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