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_tri7.h"
21 : #include "libmesh/enum_io_package.h"
22 : #include "libmesh/enum_order.h"
23 :
24 : #ifdef LIBMESH_ENABLE_AMR
25 : namespace {
26 : constexpr libMesh::Real r18 = 18;
27 : }
28 : #endif
29 :
30 : namespace libMesh
31 : {
32 :
33 :
34 :
35 :
36 : // ------------------------------------------------------------
37 : // Tri7 class static member initializations
38 : const int Tri7::num_nodes;
39 : const int Tri7::nodes_per_side;
40 :
41 :
42 : #ifdef LIBMESH_ENABLE_AMR
43 :
44 : const Real Tri7::_embedding_matrix[Tri7::num_children][Tri7::num_nodes][Tri7::num_nodes] =
45 : {
46 : // embedding matrix for child 0
47 : {
48 : // 0 1 2 3 4 5 6
49 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
50 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 1
51 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 2
52 : { .375, -.125, 0.0, .75, 0.0, 0.0, 0.0}, // 3
53 : {.09375,-.03125,-.03125, .125, -.125, .125,.84375}, // 4
54 : { .375, 0.0, -.125, 0.0, 0.0, .75, 0.0}, // 5
55 : { 5/r18,-1/r18,-1/r18,4/r18,-2/r18,4/r18, 0.5} // 6
56 : },
57 :
58 : // embedding matrix for child 1
59 : {
60 : // 0 1 2 3 4 5 6
61 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 0
62 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 1
63 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 2
64 : { -.125, .375, 0.0, .75, 0.0, 0.0, 0.0}, // 3
65 : { 0.0, .375, -.125, 0.0, .75, 0.0, 0.0}, // 4
66 : {-.03125,.09375,-.03125, .125, .125, -.125,.84375}, // 5
67 : {-1/r18, 5/r18,-1/r18,4/r18,4/r18,-2/r18, 0.5} // 6
68 : },
69 :
70 : // embedding matrix for child 2
71 : {
72 : // 0 1 2 3 4 5 6
73 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 0
74 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 1
75 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 2
76 : {-.03125,-.03125,.09375, -.125, .125, .125,.84375}, // 3
77 : { 0.0, -.125, .375, 0.0, .75, 0.0, 0.0}, // 4
78 : { -.125, 0.0, .375, 0.0, 0.0, .75, 0.0}, // 5
79 : {-1/r18,-1/r18, 5/r18,-2/r18,4/r18,4/r18, 0.5} // 6
80 : },
81 :
82 : // embedding matrix for child 3
83 : {
84 : // 0 1 2 3 4 5 6
85 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 0
86 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 1
87 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 2
88 : {-.03125,.09375,-.03125, .125, .125,-.125,.84375}, // 3
89 : {-.03125,-.03125,.09375,-.125, .125, .125,.84375}, // 4
90 : {.09375,-.03125,-.03125, .125,-.125, .125,.84375}, // 5
91 : { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 6
92 : }
93 : };
94 :
95 : const std::vector<std::pair<unsigned char, unsigned char>>
96 : Tri7::_parent_bracketing_nodes[Tri7::num_children][Tri7::num_nodes] =
97 : {
98 : // Child 0
99 : { {},{{0,1}},{{0,2}},{{0,3}},{{3,5}},{{0,5}},{{0,6}} },
100 : // Child 1
101 : { {{0,1}}, {},{{1,2}},{{1,3}},{{1,4}},{{3,4}},{{1,6}} },
102 : // Child 2
103 : { {{0,2}},{{1,2}}, {},{{4,5}},{{2,4}},{{2,5}},{{2,6}} },
104 : // Child 3
105 : { {{0,1}},{{1,2}},{{0,2}},{{3,4}},{{4,5}},{{3,5}},{{0,4},{1,5},{2,3}} }
106 : };
107 : #endif
108 :
109 :
110 :
111 : // ------------------------------------------------------------
112 : // Tri7 class member functions
113 :
114 5235677 : bool Tri7::is_vertex(const unsigned int i) const
115 : {
116 5235677 : if (i < 3)
117 2244005 : return true;
118 244100 : return false;
119 : }
120 :
121 16512 : bool Tri7::is_edge(const unsigned int i) const
122 : {
123 16512 : if (i < 3 || i == 6)
124 4128 : return false;
125 1032 : return true;
126 : }
127 :
128 336 : bool Tri7::is_face(const unsigned int i) const
129 : {
130 336 : if (i > 5)
131 336 : return true;
132 0 : return false;
133 : }
134 :
135 147156 : bool Tri7::is_node_on_side(const unsigned int n,
136 : const unsigned int s) const
137 : {
138 11066 : libmesh_assert_less (s, n_sides());
139 11066 : return std::find(std::begin(side_nodes_map[s]),
140 147156 : std::end(side_nodes_map[s]),
141 147156 : n) != std::end(side_nodes_map[s]);
142 : }
143 :
144 : std::vector<unsigned>
145 102522 : Tri7::nodes_on_side(const unsigned int s) const
146 : {
147 8520 : libmesh_assert_less(s, n_sides());
148 111042 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
149 : }
150 :
151 : std::vector<unsigned>
152 3093 : Tri7::nodes_on_edge(const unsigned int e) const
153 : {
154 3093 : return nodes_on_side(e);
155 : }
156 :
157 2404693 : bool Tri7::has_affine_map() const
158 : {
159 : // Make sure edges are straight
160 392530 : Point v = this->point(2) - this->point(1);
161 2600998 : if (!v.relative_fuzzy_equals
162 2404693 : ((this->point(4) - this->point(1))*2))
163 0 : return false;
164 2600918 : v = this->point(1) - this->point(0);
165 2600998 : if (!v.relative_fuzzy_equals
166 2404693 : ((this->point(3) - this->point(0))*2))
167 0 : return false;
168 392530 : Point v20 = this->point(2) - this->point(0);
169 2600998 : if (!v20.relative_fuzzy_equals
170 2404693 : ((this->point(5) - this->point(0))*2))
171 0 : return false;
172 :
173 : // Make sure center node is centered
174 196305 : v += v20;
175 2600998 : if (!v.relative_fuzzy_equals
176 2404693 : ((this->point(6) - this->point(0))*3))
177 0 : return false;
178 :
179 196305 : return true;
180 : }
181 :
182 :
183 :
184 89876826 : Order Tri7::default_order() const
185 : {
186 89876826 : return THIRD;
187 : }
188 :
189 :
190 :
191 13496 : Order Tri7::default_side_order() const
192 : {
193 13496 : return SECOND;
194 : }
195 :
196 :
197 :
198 0 : dof_id_type Tri7::key (const unsigned int s) const
199 : {
200 0 : libmesh_assert_less (s, this->n_sides());
201 :
202 0 : switch (s)
203 : {
204 0 : case 0:
205 :
206 : return
207 0 : this->compute_key (this->node_id(3));
208 :
209 0 : case 1:
210 :
211 : return
212 0 : this->compute_key (this->node_id(4));
213 :
214 0 : case 2:
215 :
216 : return
217 0 : this->compute_key (this->node_id(5));
218 :
219 0 : default:
220 0 : libmesh_error_msg("Invalid side s = " << s);
221 : }
222 : }
223 :
224 :
225 :
226 56081595 : unsigned int Tri7::local_side_node(unsigned int side,
227 : unsigned int side_node) const
228 : {
229 4379391 : libmesh_assert_less (side, this->n_sides());
230 4379391 : libmesh_assert_less (side_node, Tri7::nodes_per_side);
231 :
232 56081595 : return Tri7::side_nodes_map[side][side_node];
233 : }
234 :
235 :
236 :
237 1485019 : std::unique_ptr<Elem> Tri7::build_side_ptr (const unsigned int i)
238 : {
239 1485019 : return this->simple_build_side_ptr<Edge3, Tri7>(i);
240 : }
241 :
242 :
243 :
244 1661 : void Tri7::build_side_ptr (std::unique_ptr<Elem> & side,
245 : const unsigned int i)
246 : {
247 1661 : this->simple_build_side_ptr<Tri7>(side, i, EDGE3);
248 1661 : }
249 :
250 :
251 :
252 0 : void Tri7::connectivity(const unsigned int sf,
253 : const IOPackage iop,
254 : std::vector<dof_id_type> & conn) const
255 : {
256 0 : libmesh_assert_less (sf, this->n_sub_elem());
257 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
258 :
259 0 : switch (iop)
260 : {
261 0 : case TECPLOT:
262 : {
263 0 : conn.resize(4);
264 0 : switch(sf)
265 : {
266 0 : case 0:
267 : // linear sub-triangle 0
268 0 : conn[0] = this->node_id(0)+1;
269 0 : conn[1] = this->node_id(3)+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 1:
276 : // linear sub-triangle 1
277 0 : conn[0] = this->node_id(3)+1;
278 0 : conn[1] = this->node_id(1)+1;
279 0 : conn[2] = this->node_id(4)+1;
280 0 : conn[3] = this->node_id(4)+1;
281 :
282 0 : return;
283 :
284 0 : case 2:
285 : // linear sub-triangle 2
286 0 : conn[0] = this->node_id(5)+1;
287 0 : conn[1] = this->node_id(4)+1;
288 0 : conn[2] = this->node_id(2)+1;
289 0 : conn[3] = this->node_id(2)+1;
290 :
291 0 : return;
292 :
293 0 : case 3:
294 : // linear sub-triangle 3
295 0 : conn[0] = this->node_id(3)+1;
296 0 : conn[1] = this->node_id(4)+1;
297 0 : conn[2] = this->node_id(5)+1;
298 0 : conn[3] = this->node_id(5)+1;
299 :
300 0 : return;
301 :
302 0 : default:
303 0 : libmesh_error_msg("Invalid sf = " << sf);
304 : }
305 : }
306 :
307 0 : case VTK:
308 : {
309 : // VTK has a vtkBiQuadraticTriangle class whose connectivity matches libMesh's
310 0 : conn.resize(Tri7::num_nodes);
311 0 : for (auto i : index_range(conn))
312 0 : conn[i] = this->node_id(i);
313 0 : return;
314 : }
315 :
316 0 : default:
317 0 : libmesh_error_msg("Unsupported IO package " << iop);
318 : }
319 : }
320 :
321 :
322 :
323 21531913 : BoundingBox Tri7::loose_bounding_box () const
324 : {
325 : // This might have curved edges, or might be a curved surface in
326 : // 3-space, in which case the full bounding box can be larger than
327 : // the bounding box of just the nodes.
328 : //
329 : //
330 : // FIXME - I haven't yet proven the formula below to be correct for
331 : // quadratics in 2D - RHS
332 : //
333 : // FIXME - This doesn't take into account curvature caused by the
334 : // center node in 3D - RHS
335 389809 : Point pmin, pmax;
336 :
337 86127652 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
338 : {
339 65765166 : Real center = this->point(0)(d);
340 387574434 : for (unsigned int p=1; p != 6; ++p)
341 322978695 : center += this->point(p)(d);
342 64595739 : center /= 6;
343 :
344 65765166 : Real hd = std::abs(center - this->point(0)(d));
345 387574434 : for (unsigned int p=1; p != 6; ++p)
346 366792506 : hd = std::max(hd, std::abs(center - this->point(p)(d)));
347 :
348 64595739 : pmin(d) = center - hd;
349 64595739 : pmax(d) = center + hd;
350 : }
351 :
352 21921722 : return BoundingBox(pmin, pmax);
353 : }
354 :
355 :
356 :
357 12572 : unsigned int Tri7::n_second_order_adjacent_vertices (const unsigned int n) const
358 : {
359 12572 : switch (n)
360 : {
361 522 : case 3:
362 : case 4:
363 : case 5:
364 522 : return 2;
365 :
366 3143 : case 6:
367 3143 : return 3;
368 :
369 0 : default:
370 0 : libmesh_error_msg("Invalid n = " << n);
371 : }
372 : }
373 :
374 :
375 :
376 28287 : unsigned short int Tri7::second_order_adjacent_vertex (const unsigned int n,
377 : const unsigned int v) const
378 : {
379 1566 : libmesh_assert_greater_equal (n, this->n_vertices());
380 1566 : libmesh_assert_less (n, this->n_nodes());
381 :
382 28287 : switch (n)
383 : {
384 9429 : case 6:
385 : {
386 522 : libmesh_assert_less (v, 3);
387 9429 : return static_cast<unsigned short int>(v);
388 : }
389 :
390 18858 : default:
391 : {
392 1044 : libmesh_assert_less (v, 2);
393 18858 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
394 : }
395 : }
396 : }
397 :
398 :
399 :
400 : const unsigned short int Tri7::_second_order_adjacent_vertices[Tri7::num_sides][2] =
401 : {
402 : {0, 1}, // vertices adjacent to node 3
403 : {1, 2}, // vertices adjacent to node 4
404 : {0, 2} // vertices adjacent to node 5
405 : };
406 :
407 :
408 :
409 : std::pair<unsigned short int, unsigned short int>
410 0 : Tri7::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 0 : return std::pair<unsigned short int, unsigned short int>
415 0 : (_second_order_vertex_child_number[n],
416 0 : _second_order_vertex_child_index[n]);
417 : }
418 :
419 :
420 :
421 : const unsigned short int Tri7::_second_order_vertex_child_number[Tri7::num_nodes] =
422 : {
423 : 99,99,99, // Vertices
424 : 0,1,0, // Edges
425 : 3 // Interior
426 : };
427 :
428 :
429 :
430 : const unsigned short int Tri7::_second_order_vertex_child_index[Tri7::num_nodes] =
431 : {
432 : 99,99,99, // Vertices
433 : 1,2,2, // Edges
434 : 6 // Interior
435 : };
436 :
437 :
438 129834 : void Tri7::permute(unsigned int perm_num)
439 : {
440 10520 : libmesh_assert_less (perm_num, 3);
441 :
442 255665 : for (unsigned int i = 0; i != perm_num; ++i)
443 : {
444 125831 : swap3nodes(0,1,2);
445 115609 : swap3nodes(3,4,5);
446 115609 : swap3neighbors(0,1,2);
447 : }
448 129834 : }
449 :
450 :
451 3713 : void Tri7::flip(BoundaryInfo * boundary_info)
452 : {
453 174 : libmesh_assert(boundary_info);
454 :
455 3713 : swap2nodes(0,1);
456 3713 : swap2nodes(4,5);
457 174 : swap2neighbors(1,2);
458 3713 : swap2boundarysides(1,2,boundary_info);
459 3713 : swap2boundaryedges(1,2,boundary_info);
460 3713 : }
461 :
462 :
463 288 : unsigned int Tri7::center_node_on_side(const unsigned short side) const
464 : {
465 24 : libmesh_assert_less (side, Tri7::num_sides);
466 288 : return side + 3;
467 : }
468 :
469 :
470 : ElemType
471 1152 : Tri7::side_type (const unsigned int libmesh_dbg_var(s)) const
472 : {
473 96 : libmesh_assert_less (s, 3);
474 1152 : return EDGE3;
475 : }
476 :
477 :
478 : } // namespace libMesh
|
