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