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