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