face_quad9.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // Local includes
19 #include "libmesh/side.h"
20 #include "libmesh/edge_edge3.h"
21 #include "libmesh/face_quad9.h"
23 #include "libmesh/enum_order.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Quad9 class static member initializations
33 const int Quad9::num_nodes;
34 const int Quad9::num_sides;
35 const int Quad9::num_children;
36 const int Quad9::nodes_per_side;
37 
39  {
40  {0, 1, 4}, // Side 0
41  {1, 2, 5}, // Side 1
42  {2, 3, 6}, // Side 2
43  {3, 0, 7} // Side 3
44  };
45 
46 
47 #ifdef LIBMESH_ENABLE_AMR
48 
50  {
51  // embedding matrix for child 0
52  {
53  // 0 1 2 3 4 5 6 7 8
54  { 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
55  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
56  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 2
57  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 3
58  { 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
59  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 5
60  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 6
61  { 0.375000, 0.00000, 0.00000, -0.125000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
62  { 0.140625, -0.0468750, 0.0156250, -0.0468750, 0.281250, -0.0937500, -0.0937500, 0.281250, 0.562500 } // 8
63  },
64 
65  // embedding matrix for child 1
66  {
67  // 0 1 2 3 4 5 6 7 8
68  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
69  { 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
70  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 2
71  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 3
72  { -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
73  { 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
74  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 6
75  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 7
76  { -0.0468750, 0.140625, -0.0468750, 0.0156250, 0.281250, 0.281250, -0.0937500, -0.0937500, 0.562500 } // 8
77  },
78 
79  // embedding matrix for child 2
80  {
81  // 0 1 2 3 4 5 6 7 8
82  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 0
83  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 1
84  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 2
85  { 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 3
86  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 4
87  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 5
88  { 0.00000, 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
89  { -0.125000, 0.00000, 0.00000, 0.375000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
90  { -0.0468750, 0.0156250, -0.0468750, 0.140625, -0.0937500, -0.0937500, 0.281250, 0.281250, 0.562500 } // 8
91  },
92 
93  // embedding matrix for child 3
94  {
95  // 0 1 2 3 4 5 6 7 8
96  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 0
97  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 1
98  { 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 2
99  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 3
100  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 4
101  { 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
102  { 0.00000, 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
103  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 7
104  { 0.0156250, -0.0468750, 0.140625, -0.0468750, -0.0937500, 0.281250, 0.281250, -0.0937500, 0.562500 } // 8
105  }
106  };
107 
108 #endif
109 
110 
111 
112 // ------------------------------------------------------------
113 // Quad9 class member functions
114 
115 bool Quad9::is_vertex(const unsigned int i) const
116 {
117  if (i < 4)
118  return true;
119  return false;
120 }
121 
122 bool Quad9::is_edge(const unsigned int i) const
123 {
124  if (i < 4)
125  return false;
126  if (i > 7)
127  return false;
128  return true;
129 }
130 
131 bool Quad9::is_face(const unsigned int i) const
132 {
133  if (i > 7)
134  return true;
135  return false;
136 }
137 
138 bool Quad9::is_node_on_side(const unsigned int n,
139  const unsigned int s) const
140 {
141  libmesh_assert_less (s, n_sides());
142  return std::find(std::begin(side_nodes_map[s]),
144  n) != std::end(side_nodes_map[s]);
145 }
146 
147 std::vector<unsigned>
148 Quad9::nodes_on_side(const unsigned int s) const
149 {
150  libmesh_assert_less(s, n_sides());
151  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
152 }
153 
155 {
156  // make sure corners form a parallelogram
157  Point v = this->point(1) - this->point(0);
158  if (!v.relative_fuzzy_equals(this->point(2) - this->point(3)))
159  return false;
160  // make sure "horizontal" sides are straight
161  v /= 2;
162  if (!v.relative_fuzzy_equals(this->point(4) - this->point(0)) ||
163  !v.relative_fuzzy_equals(this->point(6) - this->point(3)))
164  return false;
165  // make sure "vertical" sides are straight
166  // and the center node is centered
167  v = (this->point(3) - this->point(0))/2;
168  if (!v.relative_fuzzy_equals(this->point(7) - this->point(0)) ||
169  !v.relative_fuzzy_equals(this->point(5) - this->point(1)) ||
170  !v.relative_fuzzy_equals(this->point(8) - this->point(4)))
171  return false;
172  return true;
173 }
174 
175 
176 
178 {
179  return SECOND;
180 }
181 
182 
183 
184 dof_id_type Quad9::key (const unsigned int s) const
185 {
186  libmesh_assert_less (s, this->n_sides());
187 
188  switch (s)
189  {
190  case 0:
191 
192  return
193  this->compute_key (this->node_id(4));
194 
195  case 1:
196 
197  return
198  this->compute_key (this->node_id(5));
199 
200  case 2:
201 
202  return
203  this->compute_key (this->node_id(6));
204 
205  case 3:
206 
207  return
208  this->compute_key (this->node_id(7));
209 
210  default:
211  libmesh_error_msg("Invalid side s = " << s);
212  }
213 }
214 
215 
216 
218 {
219  return this->compute_key(this->node_id(8));
220 }
221 
222 
223 
224 unsigned int Quad9::which_node_am_i(unsigned int side,
225  unsigned int side_node) const
226 {
227  libmesh_assert_less (side, this->n_sides());
228  libmesh_assert_less (side_node, Quad9::nodes_per_side);
229 
230  return Quad9::side_nodes_map[side][side_node];
231 }
232 
233 
234 
235 std::unique_ptr<Elem> Quad9::build_side_ptr (const unsigned int i,
236  bool proxy)
237 {
238  libmesh_assert_less (i, this->n_sides());
239 
240  if (proxy)
241  return libmesh_make_unique<Side<Edge3,Quad9>>(this,i);
242 
243  else
244  {
245  std::unique_ptr<Elem> edge = libmesh_make_unique<Edge3>();
246  edge->subdomain_id() = this->subdomain_id();
247 
248  // Set the nodes
249  for (unsigned n=0; n<edge->n_nodes(); ++n)
250  edge->set_node(n) = this->node_ptr(Quad9::side_nodes_map[i][n]);
251 
252  return edge;
253  }
254 }
255 
256 
257 
258 void Quad9::build_side_ptr (std::unique_ptr<Elem> & side,
259  const unsigned int i)
260 {
261  this->simple_build_side_ptr<Quad9>(side, i, EDGE3);
262 }
263 
264 
265 
266 void Quad9::connectivity(const unsigned int sf,
267  const IOPackage iop,
268  std::vector<dof_id_type> & conn) const
269 {
270  libmesh_assert_less (sf, this->n_sub_elem());
271  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
272 
273  conn.resize(4);
274 
275  switch (iop)
276  {
277  case TECPLOT:
278  {
279  switch(sf)
280  {
281  case 0:
282  // linear sub-quad 0
283  conn[0] = this->node_id(0)+1;
284  conn[1] = this->node_id(4)+1;
285  conn[2] = this->node_id(8)+1;
286  conn[3] = this->node_id(7)+1;
287  return;
288 
289  case 1:
290  // linear sub-quad 1
291  conn[0] = this->node_id(4)+1;
292  conn[1] = this->node_id(1)+1;
293  conn[2] = this->node_id(5)+1;
294  conn[3] = this->node_id(8)+1;
295  return;
296 
297  case 2:
298  // linear sub-quad 2
299  conn[0] = this->node_id(7)+1;
300  conn[1] = this->node_id(8)+1;
301  conn[2] = this->node_id(6)+1;
302  conn[3] = this->node_id(3)+1;
303  return;
304 
305  case 3:
306  // linear sub-quad 3
307  conn[0] = this->node_id(8)+1;
308  conn[1] = this->node_id(5)+1;
309  conn[2] = this->node_id(2)+1;
310  conn[3] = this->node_id(6)+1;
311  return;
312 
313  default:
314  libmesh_error_msg("Invalid sf = " << sf);
315  }
316  }
317 
318  case VTK:
319  {
320  conn.resize(9);
321  conn[0] = this->node_id(0);
322  conn[1] = this->node_id(1);
323  conn[2] = this->node_id(2);
324  conn[3] = this->node_id(3);
325  conn[4] = this->node_id(4);
326  conn[5] = this->node_id(5);
327  conn[6] = this->node_id(6);
328  conn[7] = this->node_id(7);
329  conn[8] = this->node_id(8);
330  return;
331 
332  /*
333  switch(sf)
334  {
335  case 0:
336  // linear sub-quad 0
337  conn[0] = this->node_id(0);
338  conn[1] = this->node_id(4);
339  conn[2] = this->node_id(8);
340  conn[3] = this->node_id(7);
341 
342  return;
343 
344  case 1:
345  // linear sub-quad 1
346  conn[0] = this->node_id(4);
347  conn[1] = this->node_id(1);
348  conn[2] = this->node_id(5);
349  conn[3] = this->node_id(8);
350 
351  return;
352 
353  case 2:
354  // linear sub-quad 2
355  conn[0] = this->node_id(7);
356  conn[1] = this->node_id(8);
357  conn[2] = this->node_id(6);
358  conn[3] = this->node_id(3);
359 
360  return;
361 
362  case 3:
363  // linear sub-quad 3
364  conn[0] = this->node_id(8);
365  conn[1] = this->node_id(5);
366  conn[2] = this->node_id(2);
367  conn[3] = this->node_id(6);
368 
369  return;
370 
371  default:
372  libmesh_error_msg("Invalid sf = " << sf);
373  }*/
374  }
375 
376  default:
377  libmesh_error_msg("Unsupported IO package " << iop);
378  }
379 }
380 
381 
382 
384 {
385  // This might have curved edges, or might be a curved surface in
386  // 3-space, in which case the full bounding box can be larger than
387  // the bounding box of just the nodes.
388  //
389  //
390  // FIXME - I haven't yet proven the formula below to be correct for
391  // biquadratics - RHS
392  Point pmin, pmax;
393 
394  for (unsigned d=0; d<LIBMESH_DIM; ++d)
395  {
396  const Real center = this->point(8)(d);
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 
409 
411 {
412  // Make copies of our points. It makes the subsequent calculations a bit
413  // shorter and avoids dereferencing the same pointer multiple times.
414  Point
415  x0 = point(0), x1 = point(1), x2 = point(2),
416  x3 = point(3), x4 = point(4), x5 = point(5),
417  x6 = point(6), x7 = point(7), x8 = point(8);
418 
419  // Construct constant data vectors.
420  // \vec{x}_{\xi} = \vec{a1}*xi*eta^2 + \vec{b1}*eta**2 + \vec{c1}*xi*eta + \vec{d1}*xi + \vec{e1}*eta + \vec{f1}
421  // \vec{x}_{\eta} = \vec{a2}*xi^2*eta + \vec{b2}*xi**2 + \vec{c2}*xi*eta + \vec{d2}*xi + \vec{e2}*eta + \vec{f2}
422  // This is copy-pasted directly from the output of a Python script.
423  Point
424  a1 = x0/2 + x1/2 + x2/2 + x3/2 - x4 - x5 - x6 - x7 + 2*x8,
425  b1 = -x0/4 + x1/4 + x2/4 - x3/4 - x5/2 + x7/2,
426  c1 = -x0/2 - x1/2 + x2/2 + x3/2 + x4 - x6,
427  d1 = x5 + x7 - 2*x8,
428  e1 = x0/4 - x1/4 + x2/4 - x3/4,
429  f1 = x5/2 - x7/2,
430  a2 = a1,
431  b2 = -x0/4 - x1/4 + x2/4 + x3/4 + x4/2 - x6/2,
432  c2 = -x0/2 + x1/2 + x2/2 - x3/2 - x5 + x7,
433  d2 = x0/4 - x1/4 + x2/4 - x3/4,
434  e2 = x4 + x6 - 2*x8,
435  f2 = -x4/2 + x6/2;
436 
437  // 3x3 quadrature, exact for bi-quintics
438  const unsigned int N = 3;
439  const Real q[N] = {-std::sqrt(15)/5., 0., std::sqrt(15)/5.};
440  const Real w[N] = {5./9, 8./9, 5./9};
441 
442  Real vol=0.;
443  for (unsigned int i=0; i<N; ++i)
444  for (unsigned int j=0; j<N; ++j)
445  vol += w[i] * w[j] *
446  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,
447  q[i]*q[i]*q[j]*a2 + q[i]*q[i]*b2 + q[j]*q[i]*c2 + q[i]*d2 + q[j]*e2 + f2);
448 
449  return vol;
450 }
451 
452 
453 
454 
455 unsigned int Quad9::n_second_order_adjacent_vertices (const unsigned int n) const
456 {
457  switch (n)
458  {
459  case 4:
460  case 5:
461  case 6:
462  case 7:
463  return 2;
464 
465  case 8:
466  return 4;
467 
468  default:
469  libmesh_error_msg("Invalid n = " << n);
470  }
471 }
472 
473 
474 
475 unsigned short int Quad9::second_order_adjacent_vertex (const unsigned int n,
476  const unsigned int v) const
477 {
478  libmesh_assert_greater_equal (n, this->n_vertices());
479  libmesh_assert_less (n, this->n_nodes());
480 
481  switch (n)
482  {
483  case 8:
484  {
485  libmesh_assert_less (v, 4);
486  return static_cast<unsigned short int>(v);
487  }
488 
489  default:
490  {
491  libmesh_assert_less (v, 2);
492  // use the matrix that we inherited from \p Quad
493  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
494  }
495  }
496 }
497 
498 
499 
500 std::pair<unsigned short int, unsigned short int>
501 Quad9::second_order_child_vertex (const unsigned int n) const
502 {
503  libmesh_assert_greater_equal (n, this->n_vertices());
504  libmesh_assert_less (n, this->n_nodes());
505  /*
506  * the _second_order_vertex_child_* vectors are
507  * stored in face_quad.C, since they are identical
508  * for Quad8 and Quad9 (for the first 4 higher-order nodes)
509  */
510  return std::pair<unsigned short int, unsigned short int>
513 }
514 
515 } // namespace libMesh
double abs(double a)
virtual unsigned int n_sub_elem() const override
Definition: face_quad9.h:81
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: face_quad9.C:266
virtual unsigned int n_vertices() const override final
Definition: face_quad.h:96
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
Definition: face_quad9.h:203
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const override
Definition: face_quad9.C:475
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const override
Definition: face_quad9.C:501
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: face_quad9.C:148
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1097
unsigned short int side
Definition: xdr_io.C:50
virtual unsigned int which_node_am_i(unsigned int side, unsigned int side_node) const override
Definition: face_quad9.C:224
IterBase * end
static const int num_children
Definition: face_quad9.h:196
long double max(long double a, double b)
virtual bool is_edge(const unsigned int i) const override
Definition: face_quad9.C:122
virtual dof_id_type key() const override
Definition: face_quad9.C:217
virtual bool has_affine_map() const override
Definition: face_quad9.C:154
static const int num_sides
Definition: face_quad9.h:195
virtual bool is_vertex(const unsigned int i) const override
Definition: face_quad9.C:115
virtual Order default_order() const override
Definition: face_quad9.C:177
static const float _embedding_matrix[num_children][num_nodes][num_nodes]
Definition: face_quad9.h:240
static const unsigned short int _second_order_vertex_child_index[9]
Definition: face_quad.h:202
static const unsigned short int _second_order_vertex_child_number[9]
Definition: face_quad.h:197
virtual unsigned int n_nodes() const override
Definition: face_quad9.h:76
static const unsigned short int _second_order_adjacent_vertices[4][2]
Definition: face_quad.h:192
virtual bool is_face(const unsigned int i) const override
Definition: face_quad9.C:131
virtual unsigned int n_sides() const override final
Definition: face_quad.h:91
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const override
Definition: face_quad9.C:455
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1957
virtual Real volume() const override
Definition: face_quad9.C:410
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: face_quad9.C:138
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:2754
static const int num_nodes
Definition: face_quad9.h:194
A geometric point in (x,y,z) space.
Definition: point.h:38
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1914
const Point & point(const unsigned int i) const
Definition: elem.h:1892
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:990
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy) override
Definition: face_quad9.C:235
static const int nodes_per_side
Definition: face_quad9.h:197
std::unique_ptr< Elem > side(const unsigned int i) const
Definition: elem.h:2202
uint8_t dof_id_type
Definition: id_types.h:64
virtual BoundingBox loose_bounding_box() const override
Definition: face_quad9.C:383