libMesh
face_quad9.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 
20 // Local includes
21 #include "libmesh/side.h"
22 #include "libmesh/edge_edge3.h"
23 #include "libmesh/face_quad9.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Quad9 class static member initializations
33 const unsigned int Quad9::side_nodes_map[4][3] =
34  {
35  {0, 1, 4}, // Side 0
36  {1, 2, 5}, // Side 1
37  {2, 3, 6}, // Side 2
38  {3, 0, 7} // Side 3
39  };
40 
41 
42 #ifdef LIBMESH_ENABLE_AMR
43 
44 const float Quad9::_embedding_matrix[4][9][9] =
45  {
46  // embedding matrix for child 0
47  {
48  // 0 1 2 3 4 5 6 7 8
49  { 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
50  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
51  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 2
52  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 3
53  { 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
54  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 5
55  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 6
56  { 0.375000, 0.00000, 0.00000, -0.125000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
57  { 0.140625, -0.0468750, 0.0156250, -0.0468750, 0.281250, -0.0937500, -0.0937500, 0.281250, 0.562500 } // 8
58  },
59 
60  // embedding matrix for child 1
61  {
62  // 0 1 2 3 4 5 6 7 8
63  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
64  { 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
65  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 2
66  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 3
67  { -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
68  { 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
69  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 6
70  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 7
71  { -0.0468750, 0.140625, -0.0468750, 0.0156250, 0.281250, 0.281250, -0.0937500, -0.0937500, 0.562500 } // 8
72  },
73 
74  // embedding matrix for child 2
75  {
76  // 0 1 2 3 4 5 6 7 8
77  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 0
78  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 1
79  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 2
80  { 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 3
81  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 4
82  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 5
83  { 0.00000, 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
84  { -0.125000, 0.00000, 0.00000, 0.375000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
85  { -0.0468750, 0.0156250, -0.0468750, 0.140625, -0.0937500, -0.0937500, 0.281250, 0.281250, 0.562500 } // 8
86  },
87 
88  // embedding matrix for child 3
89  {
90  // 0 1 2 3 4 5 6 7 8
91  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 0
92  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 1
93  { 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 2
94  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 3
95  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 4
96  { 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
97  { 0.00000, 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
98  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 7
99  { 0.0156250, -0.0468750, 0.140625, -0.0468750, -0.0937500, 0.281250, 0.281250, -0.0937500, 0.562500 } // 8
100  }
101  };
102 
103 #endif
104 
105 
106 
107 // ------------------------------------------------------------
108 // Quad9 class member functions
109 
110 bool Quad9::is_vertex(const unsigned int i) const
111 {
112  if (i < 4)
113  return true;
114  return false;
115 }
116 
117 bool Quad9::is_edge(const unsigned int i) const
118 {
119  if (i < 4)
120  return false;
121  if (i > 7)
122  return false;
123  return true;
124 }
125 
126 bool Quad9::is_face(const unsigned int i) const
127 {
128  if (i > 7)
129  return true;
130  return false;
131 }
132 
133 bool Quad9::is_node_on_side(const unsigned int n,
134  const unsigned int s) const
135 {
136  libmesh_assert_less (s, n_sides());
137  for (unsigned int i = 0; i != 3; ++i)
138  if (side_nodes_map[s][i] == n)
139  return true;
140  return false;
141 }
142 
143 
144 
146 {
147  // make sure corners form a parallelogram
148  Point v = this->point(1) - this->point(0);
149  if (!v.relative_fuzzy_equals(this->point(2) - this->point(3)))
150  return false;
151  // make sure "horizontal" sides are straight
152  v /= 2;
153  if (!v.relative_fuzzy_equals(this->point(4) - this->point(0)) ||
154  !v.relative_fuzzy_equals(this->point(6) - this->point(3)))
155  return false;
156  // make sure "vertical" sides are straight
157  // and the center node is centered
158  v = (this->point(3) - this->point(0))/2;
159  if (!v.relative_fuzzy_equals(this->point(7) - this->point(0)) ||
160  !v.relative_fuzzy_equals(this->point(5) - this->point(1)) ||
161  !v.relative_fuzzy_equals(this->point(8) - this->point(4)))
162  return false;
163  return true;
164 }
165 
166 
167 
168 dof_id_type Quad9::key (const unsigned int s) const
169 {
170  libmesh_assert_less (s, this->n_sides());
171 
172  switch (s)
173  {
174  case 0:
175 
176  return
177  this->compute_key (this->node_id(4));
178 
179  case 1:
180 
181  return
182  this->compute_key (this->node_id(5));
183 
184  case 2:
185 
186  return
187  this->compute_key (this->node_id(6));
188 
189  case 3:
190 
191  return
192  this->compute_key (this->node_id(7));
193 
194  default:
195  libmesh_error_msg("Invalid side s = " << s);
196  }
197 
198  libmesh_error_msg("We'll never get here!");
199  return 0;
200 }
201 
202 
203 
205 {
206  return this->compute_key(this->node_id(8));
207 }
208 
209 
210 
211 unsigned int Quad9::which_node_am_i(unsigned int side,
212  unsigned int side_node) const
213 {
214  libmesh_assert_less (side, this->n_sides());
215  libmesh_assert_less (side_node, 3);
216 
217  return Quad9::side_nodes_map[side][side_node];
218 }
219 
220 
221 
223  bool proxy)
224 {
225  libmesh_assert_less (i, this->n_sides());
226 
227  if (proxy)
228  return UniquePtr<Elem>(new Side<Edge3,Quad9>(this,i));
229 
230  else
231  {
232  Elem * edge = new Edge3;
233  edge->subdomain_id() = this->subdomain_id();
234 
235  // Set the nodes
236  for (unsigned n=0; n<edge->n_nodes(); ++n)
237  edge->set_node(n) = this->node_ptr(Quad9::side_nodes_map[i][n]);
238 
239  return UniquePtr<Elem>(edge);
240  }
241 
242  libmesh_error_msg("We'll never get here!");
243  return UniquePtr<Elem>();
244 }
245 
246 
247 
248 
249 
250 
251 
252 void Quad9::connectivity(const unsigned int sf,
253  const IOPackage iop,
254  std::vector<dof_id_type> & conn) const
255 {
256  libmesh_assert_less (sf, this->n_sub_elem());
257  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
258 
259  conn.resize(4);
260 
261  switch (iop)
262  {
263  case TECPLOT:
264  {
265  switch(sf)
266  {
267  case 0:
268  // linear sub-quad 0
269  conn[0] = this->node_id(0)+1;
270  conn[1] = this->node_id(4)+1;
271  conn[2] = this->node_id(8)+1;
272  conn[3] = this->node_id(7)+1;
273  return;
274 
275  case 1:
276  // linear sub-quad 1
277  conn[0] = this->node_id(4)+1;
278  conn[1] = this->node_id(1)+1;
279  conn[2] = this->node_id(5)+1;
280  conn[3] = this->node_id(8)+1;
281  return;
282 
283  case 2:
284  // linear sub-quad 2
285  conn[0] = this->node_id(7)+1;
286  conn[1] = this->node_id(8)+1;
287  conn[2] = this->node_id(6)+1;
288  conn[3] = this->node_id(3)+1;
289  return;
290 
291  case 3:
292  // linear sub-quad 3
293  conn[0] = this->node_id(8)+1;
294  conn[1] = this->node_id(5)+1;
295  conn[2] = this->node_id(2)+1;
296  conn[3] = this->node_id(6)+1;
297  return;
298 
299  default:
300  libmesh_error_msg("Invalid sf = " << sf);
301  }
302  }
303 
304  case VTK:
305  {
306  conn.resize(9);
307  conn[0] = this->node_id(0);
308  conn[1] = this->node_id(1);
309  conn[2] = this->node_id(2);
310  conn[3] = this->node_id(3);
311  conn[4] = this->node_id(4);
312  conn[5] = this->node_id(5);
313  conn[6] = this->node_id(6);
314  conn[7] = this->node_id(7);
315  conn[8] = this->node_id(8);
316  return;
317 
318  /*
319  switch(sf)
320  {
321  case 0:
322  // linear sub-quad 0
323  conn[0] = this->node_id(0);
324  conn[1] = this->node_id(4);
325  conn[2] = this->node_id(8);
326  conn[3] = this->node_id(7);
327 
328  return;
329 
330  case 1:
331  // linear sub-quad 1
332  conn[0] = this->node_id(4);
333  conn[1] = this->node_id(1);
334  conn[2] = this->node_id(5);
335  conn[3] = this->node_id(8);
336 
337  return;
338 
339  case 2:
340  // linear sub-quad 2
341  conn[0] = this->node_id(7);
342  conn[1] = this->node_id(8);
343  conn[2] = this->node_id(6);
344  conn[3] = this->node_id(3);
345 
346  return;
347 
348  case 3:
349  // linear sub-quad 3
350  conn[0] = this->node_id(8);
351  conn[1] = this->node_id(5);
352  conn[2] = this->node_id(2);
353  conn[3] = this->node_id(6);
354 
355  return;
356 
357  default:
358  libmesh_error_msg("Invalid sf = " << sf);
359  }*/
360  }
361 
362  default:
363  libmesh_error_msg("Unsupported IO package " << iop);
364  }
365 }
366 
367 
368 
370 {
371  // This might have curved edges, or might be a curved surface in
372  // 3-space, in which case the full bounding box can be larger than
373  // the bounding box of just the nodes.
374  //
375  //
376  // FIXME - I haven't yet proven the formula below to be correct for
377  // biquadratics - RHS
378  Point pmin, pmax;
379 
380  for (unsigned d=0; d<LIBMESH_DIM; ++d)
381  {
382  const Real center = this->point(8)(d);
383  Real hd = std::abs(center - this->point(0)(d));
384  for (unsigned int p=0; p != 8; ++p)
385  hd = std::max(hd, std::abs(center - this->point(p)(d)));
386 
387  pmin(d) = center - hd;
388  pmax(d) = center + hd;
389  }
390 
391  return BoundingBox(pmin, pmax);
392 }
393 
394 
395 
397 {
398  // Make copies of our points. It makes the subsequent calculations a bit
399  // shorter and avoids dereferencing the same pointer multiple times.
400  Point
401  x0 = point(0), x1 = point(1), x2 = point(2),
402  x3 = point(3), x4 = point(4), x5 = point(5),
403  x6 = point(6), x7 = point(7), x8 = point(8);
404 
405  // Construct constant data vectors.
406  // \vec{x}_{\xi} = \vec{a1}*xi*eta^2 + \vec{b1}*eta**2 + \vec{c1}*xi*eta + \vec{d1}*xi + \vec{e1}*eta + \vec{f1}
407  // \vec{x}_{\eta} = \vec{a2}*xi^2*eta + \vec{b2}*xi**2 + \vec{c2}*xi*eta + \vec{d2}*xi + \vec{e2}*eta + \vec{f2}
408  // This is copy-pasted directly from the output of a Python script.
409  Point
410  a1 = x0/2 + x1/2 + x2/2 + x3/2 - x4 - x5 - x6 - x7 + 2*x8,
411  b1 = -x0/4 + x1/4 + x2/4 - x3/4 - x5/2 + x7/2,
412  c1 = -x0/2 - x1/2 + x2/2 + x3/2 + x4 - x6,
413  d1 = x5 + x7 - 2*x8,
414  e1 = x0/4 - x1/4 + x2/4 - x3/4,
415  f1 = x5/2 - x7/2,
416  a2 = a1,
417  b2 = -x0/4 - x1/4 + x2/4 + x3/4 + x4/2 - x6/2,
418  c2 = -x0/2 + x1/2 + x2/2 - x3/2 - x5 + x7,
419  d2 = x0/4 - x1/4 + x2/4 - x3/4,
420  e2 = x4 + x6 - 2*x8,
421  f2 = -x4/2 + x6/2;
422 
423  // 3x3 quadrature, exact for bi-quintics
424  const unsigned int N = 3;
425  const Real q[N] = {-std::sqrt(15)/5., 0., std::sqrt(15)/5.};
426  const Real w[N] = {5./9, 8./9, 5./9};
427 
428  Real vol=0.;
429  for (unsigned int i=0; i<N; ++i)
430  for (unsigned int j=0; j<N; ++j)
431  vol += w[i] * w[j] *
432  cross_norm(q[i]*q[j]*q[j]*a1 + q[j]*q[j]*b1 + q[j]*q[i]*c1 + q[i]*d1 + q[j]*e1 + f1,
433  q[i]*q[i]*q[j]*a2 + q[i]*q[i]*b2 + q[j]*q[i]*c2 + q[i]*d2 + q[j]*e2 + f2);
434 
435  return vol;
436 }
437 
438 
439 
440 
441 unsigned int Quad9::n_second_order_adjacent_vertices (const unsigned int n) const
442 {
443  switch (n)
444  {
445  case 4:
446  case 5:
447  case 6:
448  case 7:
449  return 2;
450 
451  case 8:
452  return 4;
453 
454  default:
455  libmesh_error_msg("Invalid n = " << n);
456  }
457 
458  libmesh_error_msg("We'll never get here!");
459  return libMesh::invalid_uint;
460 }
461 
462 
463 
464 unsigned short int Quad9::second_order_adjacent_vertex (const unsigned int n,
465  const unsigned int v) const
466 {
467  libmesh_assert_greater_equal (n, this->n_vertices());
468  libmesh_assert_less (n, this->n_nodes());
469 
470  switch (n)
471  {
472  case 8:
473  {
474  libmesh_assert_less (v, 4);
475  return static_cast<unsigned short int>(v);
476  }
477 
478  default:
479  {
480  libmesh_assert_less (v, 2);
481  // use the matrix that we inherited from \p Quad
482  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
483  }
484  }
485 }
486 
487 
488 
489 std::pair<unsigned short int, unsigned short int>
490 Quad9::second_order_child_vertex (const unsigned int n) const
491 {
492  libmesh_assert_greater_equal (n, this->n_vertices());
493  libmesh_assert_less (n, this->n_nodes());
494  /*
495  * the _second_order_vertex_child_* vectors are
496  * stored in face_quad.C, since they are identical
497  * for Quad8 and Quad9 (for the first 4 higher-order nodes)
498  */
499  return std::pair<unsigned short int, unsigned short int>
502 }
503 
504 } // namespace libMesh
virtual unsigned int n_nodes() const libmesh_override
Definition: face_quad9.h:68
double abs(double a)
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
virtual dof_id_type key() const libmesh_override
Compute a unique key for this element which is suitable for hashing (not necessarily unique...
Definition: face_quad9.C:204
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
virtual unsigned int n_sub_elem() const libmesh_override
Definition: face_quad9.h:73
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_quad9.C:133
UniquePtr< Elem > side(const unsigned int i) const
Definition: elem.h:2105
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Calls cross_norm_sq() and takes the square root of the result.
Definition: type_vector.h:1092
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_quad9.C:222
virtual Real volume() const libmesh_override
An optimized method for approximating the area of a QUAD9 using quadrature.
Definition: face_quad9.C:396
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static const unsigned int side_nodes_map[4][3]
This maps the node of the side to element node numbers.
Definition: face_quad9.h:179
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
static const float _embedding_matrix[4][9][9]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition: face_quad9.h:216
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
virtual unsigned int n_vertices() const libmesh_override
Definition: face_quad.h:90
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
virtual BoundingBox loose_bounding_box() const libmesh_override
Definition: face_quad9.C:369
static const unsigned short int _second_order_vertex_child_index[9]
Vector that names the child vertex index for each second order node.
Definition: face_quad.h:190
static const unsigned short int _second_order_vertex_child_number[9]
Vector that names a child sharing each second order node.
Definition: face_quad.h:185
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
The Edge3 is an element in 1D composed of 3 nodes.
Definition: edge_edge3.h:43
static const unsigned short int _second_order_adjacent_vertices[4][2]
Matrix that tells which vertices define the location of mid-side (or second-order) nodes...
Definition: face_quad.h:180
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const libmesh_override
Definition: face_quad9.C:441
This defines the Side class.
Definition: side.h:48
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const libmesh_override
Definition: face_quad9.C:490
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const libmesh_override
Definition: face_quad9.C:464
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const libmesh_override
Definition: face_quad9.C:211
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual unsigned int n_sides() const libmesh_override
Definition: face_quad.h:85
virtual bool is_edge(const unsigned int i) const libmesh_override
Definition: face_quad9.C:117
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_quad9.C:126
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_quad9.C:110
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:2621
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_quad9.C:252
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:987
virtual bool has_affine_map() const libmesh_override
Definition: face_quad9.C:145
uint8_t dof_id_type
Definition: id_types.h:64