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