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