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