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