libMesh
face_quad8.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_quad8.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Quad8 class static member initializations
33 const unsigned int Quad8::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 Quad8::_embedding_matrix[4][8][8] =
45  {
46  // embedding matrix for child 0
47  {
48  // 0 1 2 3 4 5 6 7
49  { 1.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 }, // 1
51  { -0.250000, -0.250000, -0.250000, -0.250000, 0.500000, 0.500000, 0.500000, 0.500000 }, // 2
52  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 3
53  { 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 4
54  { -0.187500, -0.187500, -0.187500, -0.187500, 0.750000, 0.375000, 0.250000, 0.375000 }, // 5
55  { -0.187500, -0.187500, -0.187500, -0.187500, 0.375000, 0.250000, 0.375000, 0.750000 }, // 6
56  { 0.375000, 0.00000, 0.00000, -0.125000, 0.00000, 0.00000, 0.00000, 0.750000 } // 7
57  },
58 
59  // embedding matrix for child 1
60  {
61  // 0 1 2 3 4 5 6 7
62  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 0
63  { 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
64  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 2
65  { -0.250000, -0.250000, -0.250000, -0.250000, 0.500000, 0.500000, 0.500000, 0.500000 }, // 3
66  { -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 4
67  { 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 5
68  { -0.187500, -0.187500, -0.187500, -0.187500, 0.375000, 0.750000, 0.375000, 0.250000 }, // 6
69  { -0.187500, -0.187500, -0.187500, -0.187500, 0.750000, 0.375000, 0.250000, 0.375000 } // 7
70  },
71 
72  // embedding matrix for child 2
73  {
74  // 0 1 2 3 4 5 6 7
75  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 0
76  { -0.250000, -0.250000, -0.250000, -0.250000, 0.500000, 0.500000, 0.500000, 0.500000 }, // 1
77  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 2
78  { 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 3
79  { -0.187500, -0.187500, -0.187500, -0.187500, 0.375000, 0.250000, 0.375000, 0.750000 }, // 4
80  { -0.187500, -0.187500, -0.187500, -0.187500, 0.250000, 0.375000, 0.750000, 0.375000 }, // 5
81  { 0.00000, 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 6
82  { -0.125000, 0.00000, 0.00000, 0.375000, 0.00000, 0.00000, 0.00000, 0.750000 } // 7
83  },
84 
85  // embedding matrix for child 3
86  {
87  // 0 1 2 3 4 5 6 7
88  { -0.250000, -0.250000, -0.250000, -0.250000, 0.500000, 0.500000, 0.500000, 0.500000 }, // 0
89  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 1
90  { 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 2
91  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 3
92  { -0.187500, -0.187500, -0.187500, -0.187500, 0.375000, 0.750000, 0.375000, 0.250000 }, // 4
93  { 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 5
94  { 0.00000, 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 6
95  { -0.187500, -0.187500, -0.187500, -0.187500, 0.250000, 0.375000, 0.750000, 0.375000 } // 7
96  }
97  };
98 
99 
100 #endif
101 
102 
103 // ------------------------------------------------------------
104 // Quad8 class member functions
105 
106 bool Quad8::is_vertex(const unsigned int i) const
107 {
108  if (i < 4)
109  return true;
110  return false;
111 }
112 
113 bool Quad8::is_edge(const unsigned int i) const
114 {
115  if (i < 4)
116  return false;
117  return true;
118 }
119 
120 bool Quad8::is_face(const unsigned int) const
121 {
122  return false;
123 }
124 
125 bool Quad8::is_node_on_side(const unsigned int n,
126  const unsigned int s) const
127 {
128  libmesh_assert_less (s, n_sides());
129  for (unsigned int i = 0; i != 3; ++i)
130  if (side_nodes_map[s][i] == n)
131  return true;
132  return false;
133 }
134 
135 
136 
138 {
139  // make sure corners form a parallelogram
140  Point v = this->point(1) - this->point(0);
141  if (!v.relative_fuzzy_equals(this->point(2) - this->point(3)))
142  return false;
143  // make sure sides are straight
144  v /= 2;
145  if (!v.relative_fuzzy_equals(this->point(4) - this->point(0)) ||
146  !v.relative_fuzzy_equals(this->point(6) - this->point(3)))
147  return false;
148  v = (this->point(3) - this->point(0))/2;
149  if (!v.relative_fuzzy_equals(this->point(7) - this->point(0)) ||
150  !v.relative_fuzzy_equals(this->point(5) - this->point(1)))
151  return false;
152  return true;
153 }
154 
155 
156 
157 dof_id_type Quad8::key (const unsigned int s) const
158 {
159  libmesh_assert_less (s, this->n_sides());
160 
161  switch (s)
162  {
163  case 0:
164 
165  return
166  this->compute_key (this->node_id(4));
167 
168  case 1:
169 
170  return
171  this->compute_key (this->node_id(5));
172 
173  case 2:
174 
175  return
176  this->compute_key (this->node_id(6));
177 
178  case 3:
179 
180  return
181  this->compute_key (this->node_id(7));
182 
183  default:
184  libmesh_error_msg("Invalid side s = " << s);
185  }
186 
187  libmesh_error_msg("We'll never get here!");
188  return 0;
189 }
190 
191 
192 
193 unsigned int Quad8::which_node_am_i(unsigned int side,
194  unsigned int side_node) const
195 {
196  libmesh_assert_less (side, this->n_sides());
197  libmesh_assert_less (side_node, 3);
198 
199  return Quad8::side_nodes_map[side][side_node];
200 }
201 
202 
203 
205  bool proxy)
206 {
207  libmesh_assert_less (i, this->n_sides());
208 
209  if (proxy)
210  return UniquePtr<Elem>(new Side<Edge3,Quad8>(this,i));
211 
212  else
213  {
214  Elem * edge = new Edge3;
215  edge->subdomain_id() = this->subdomain_id();
216 
217  // Set the nodes
218  for (unsigned n=0; n<edge->n_nodes(); ++n)
219  edge->set_node(n) = this->node_ptr(Quad8::side_nodes_map[i][n]);
220 
221  return UniquePtr<Elem>(edge);
222  }
223 
224  libmesh_error_msg("We'll never get here!");
225  return UniquePtr<Elem>();
226 }
227 
228 
229 
230 
231 
232 
233 void Quad8::connectivity(const unsigned int sf,
234  const IOPackage iop,
235  std::vector<dof_id_type> & conn) const
236 {
237  libmesh_assert_less (sf, this->n_sub_elem());
238  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
239 
240  switch (iop)
241  {
242  // Note: TECPLOT connectivity is output as four triangles with
243  // a central quadrilateral. Therefore, the first four connectivity
244  // arrays are degenerate quads (triangles in Tecplot).
245  case TECPLOT:
246  {
247  // Create storage
248  conn.resize(4);
249 
250  switch(sf)
251  {
252  case 0:
253  // linear sub-tri 0
254  conn[0] = this->node_id(0)+1;
255  conn[1] = this->node_id(4)+1;
256  conn[2] = this->node_id(7)+1;
257  conn[3] = this->node_id(7)+1;
258 
259  return;
260 
261  case 1:
262  // linear sub-tri 1
263  conn[0] = this->node_id(4)+1;
264  conn[1] = this->node_id(1)+1;
265  conn[2] = this->node_id(5)+1;
266  conn[3] = this->node_id(5)+1;
267 
268  return;
269 
270  case 2:
271  // linear sub-tri 2
272  conn[0] = this->node_id(5)+1;
273  conn[1] = this->node_id(2)+1;
274  conn[2] = this->node_id(6)+1;
275  conn[3] = this->node_id(6)+1;
276 
277  return;
278 
279  case 3:
280  // linear sub-tri 3
281  conn[0] = this->node_id(7)+1;
282  conn[1] = this->node_id(6)+1;
283  conn[2] = this->node_id(3)+1;
284  conn[3] = this->node_id(3)+1;
285 
286  return;
287 
288  case 4:
289  // linear sub-quad
290  conn[0] = this->node_id(4)+1;
291  conn[1] = this->node_id(5)+1;
292  conn[2] = this->node_id(6)+1;
293  conn[3] = this->node_id(7)+1;
294 
295  return;
296 
297  default:
298  libmesh_error_msg("Invalid sf = " << sf);
299  }
300  }
301 
302 
303  // Note: VTK connectivity is output as four triangles with
304  // a central quadrilateral. Therefore most of the connectivity
305  // arrays have length three.
306  case VTK:
307  {
308  // Create storage
309  conn.resize(8);
310  conn[0] = this->node_id(0);
311  conn[1] = this->node_id(1);
312  conn[2] = this->node_id(2);
313  conn[3] = this->node_id(3);
314  conn[4] = this->node_id(4);
315  conn[5] = this->node_id(5);
316  conn[6] = this->node_id(6);
317  conn[7] = this->node_id(7);
318  return;
319  /*
320  conn.resize(3);
321 
322  switch (sf)
323  {
324  case 0:
325  // linear sub-tri 0
326  conn[0] = this->node_id(0);
327  conn[1] = this->node_id(4);
328  conn[2] = this->node_id(7);
329 
330  return;
331 
332  case 1:
333  // linear sub-tri 1
334  conn[0] = this->node_id(4);
335  conn[1] = this->node_id(1);
336  conn[2] = this->node_id(5);
337 
338  return;
339 
340  case 2:
341  // linear sub-tri 2
342  conn[0] = this->node_id(5);
343  conn[1] = this->node_id(2);
344  conn[2] = this->node_id(6);
345 
346  return;
347 
348  case 3:
349  // linear sub-tri 3
350  conn[0] = this->node_id(7);
351  conn[1] = this->node_id(6);
352  conn[2] = this->node_id(3);
353 
354  return;
355 
356  case 4:
357  conn.resize(4);
358 
359  // linear sub-quad
360  conn[0] = this->node_id(4);
361  conn[1] = this->node_id(5);
362  conn[2] = this->node_id(6);
363  conn[3] = this->node_id(7);
364  */
365  // return;
366 
367  // default:
368  // libmesh_error_msg("Invalid sf = " << sf);
369  // }
370  }
371 
372  default:
373  libmesh_error_msg("Unsupported IO package " << iop);
374  }
375 }
376 
377 
378 
380 {
381  // This might have curved edges, or might be a curved surface in
382  // 3-space, in which case the full bounding box can be larger than
383  // the bounding box of just the nodes.
384  //
385  //
386  // FIXME - I haven't yet proven the formula below to be correct for
387  // biquadratics - RHS
388  Point pmin, pmax;
389 
390  for (unsigned d=0; d<LIBMESH_DIM; ++d)
391  {
392  Real center = this->point(0)(d);
393  for (unsigned int p=1; p != 8; ++p)
394  center += this->point(p)(d);
395  center /= 8;
396 
397  Real hd = std::abs(center - this->point(0)(d));
398  for (unsigned int p=0; p != 8; ++p)
399  hd = std::max(hd, std::abs(center - this->point(p)(d)));
400 
401  pmin(d) = center - hd;
402  pmax(d) = center + hd;
403  }
404 
405  return BoundingBox(pmin, pmax);
406 }
407 
408 
410 {
411  // Make copies of our points. It makes the subsequent calculations a bit
412  // shorter and avoids dereferencing the same pointer multiple times.
413  Point
414  x0 = point(0),
415  x1 = point(1),
416  x2 = point(2),
417  x3 = point(3),
418  x4 = point(4),
419  x5 = point(5),
420  x6 = point(6),
421  x7 = point(7);
422 
423  // Construct constant data vectors.
424  // \vec{x}_{\xi} = \vec{a1}*eta**2 + \vec{b1}*xi*eta + \vec{c1}*xi + \vec{d1}*eta + \vec{e1}
425  // \vec{x}_{\eta} = \vec{a2}*xi**2 + \vec{b2}*xi*eta + \vec{c2}*xi + \vec{d2}*eta + \vec{e2}
426  // This is copy-pasted directly from the output of a Python script.
427  Point
428  a1 = -x0/4 + x1/4 + x2/4 - x3/4 - x5/2 + x7/2,
429  b1 = -x0/2 - x1/2 + x2/2 + x3/2 + x4 - x6,
430  c1 = x0/2 + x1/2 + x2/2 + x3/2 - x4 - x6,
431  d1 = x0/4 - x1/4 + x2/4 - x3/4,
432  e1 = x5/2 - x7/2,
433  a2 = -x0/4 - x1/4 + x2/4 + x3/4 + x4/2 - x6/2,
434  b2 = -x0/2 + x1/2 + x2/2 - x3/2 - x5 + x7,
435  c2 = x0/4 - x1/4 + x2/4 - x3/4,
436  d2 = x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
437  e2 = -x4/2 + x6/2;
438 
439  // 3x3 quadrature, exact for bi-quintics
440  const unsigned int N = 3;
441  const Real q[N] = {-std::sqrt(15)/5., 0., std::sqrt(15)/5.};
442  const Real w[N] = {5./9, 8./9, 5./9};
443 
444  Real vol=0.;
445  for (unsigned int i=0; i<N; ++i)
446  for (unsigned int j=0; j<N; ++j)
447  vol += w[i] * w[j] * cross_norm(q[j]*q[j]*a1 + q[i]*q[j]*b1 + q[i]*c1 + q[j]*d1 + e1,
448  q[i]*q[i]*a2 + q[i]*q[j]*b2 + q[i]*c2 + q[j]*d2 + e2);
449 
450  return vol;
451 }
452 
453 
454 
455 unsigned short int Quad8::second_order_adjacent_vertex (const unsigned int n,
456  const unsigned int v) const
457 {
458  libmesh_assert_greater_equal (n, this->n_vertices());
459  libmesh_assert_less (n, this->n_nodes());
460  libmesh_assert_less (v, 2);
461  // use the matrix from \p face_quad.C
462  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
463 }
464 
465 
466 
467 std::pair<unsigned short int, unsigned short int>
468 Quad8::second_order_child_vertex (const unsigned int n) const
469 {
470  libmesh_assert_greater_equal (n, this->n_vertices());
471  libmesh_assert_less (n, this->n_nodes());
472  /*
473  * the _second_order_vertex_child_* vectors are
474  * stored in face_quad.C, since they are identical
475  * for Quad8 and Quad9 (for the first 4 higher-order nodes)
476  */
477  return std::pair<unsigned short int, unsigned short int>
480 }
481 
482 } // namespace libMesh
double abs(double a)
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
virtual unsigned int n_nodes() const libmesh_override
Definition: face_quad8.h:68
static const float _embedding_matrix[4][8][8]
Matrix that computes new nodal locations/solution values from current nodes/solution.
Definition: face_quad8.h:209
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_quad8.C:125
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_quad8.C:233
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const libmesh_override
Definition: face_quad8.C:468
virtual dof_id_type key() const libmesh_override
Definition: face_quad.C:74
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 unsigned int n_sub_elem() const libmesh_override
Definition: face_quad8.h:73
static const unsigned int side_nodes_map[4][3]
This maps the node of the side to element node numbers.
Definition: face_quad8.h:172
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_quad8.C:106
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_quad8.C:204
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 bool is_edge(const unsigned int i) const libmesh_override
Definition: face_quad8.C:113
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 Real volume() const libmesh_override
An optimized method for approximating the area of a QUAD8 using quadrature.
Definition: face_quad8.C:409
This defines the Side class.
Definition: side.h:48
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const libmesh_override
Definition: face_quad8.C:455
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
virtual BoundingBox loose_bounding_box() const libmesh_override
Definition: face_quad8.C:379
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const libmesh_override
Definition: face_quad8.C:193
virtual unsigned int n_sides() const libmesh_override
Definition: face_quad.h:85
virtual bool has_affine_map() const libmesh_override
Definition: face_quad8.C:137
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.
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
uint8_t dof_id_type
Definition: id_types.h:64
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_quad8.C:120