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