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_tri6.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 : // Tri6 class static member initializations
32 : const int Tri6::num_nodes;
33 : const int Tri6::nodes_per_side;
34 :
35 : const unsigned int Tri6::side_nodes_map[Tri6::num_sides][Tri6::nodes_per_side] =
36 : {
37 : {0, 1, 3}, // Side 0
38 : {1, 2, 4}, // Side 1
39 : {2, 0, 5} // Side 2
40 : };
41 :
42 :
43 : #ifdef LIBMESH_ENABLE_AMR
44 :
45 : const Real Tri6::_embedding_matrix[Tri6::num_children][Tri6::num_nodes][Tri6::num_nodes] =
46 : {
47 : // embedding matrix for child 0
48 : {
49 : // 0 1 2 3 4 5
50 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
51 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 1
52 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 2
53 : {.375, -.125, 0.0, .75, 0.0, 0.0}, // 3
54 : { 0.0, -.125, -.125, 0.5, .25, 0.5}, // 4
55 : {.375, 0.0, -.125, 0.0, 0.0, .75} // 5
56 : },
57 :
58 : // embedding matrix for child 1
59 : {
60 : // 0 1 2 3 4 5
61 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 0
62 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
63 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 2
64 : {-.125, .375, 0.0, .75, 0.0, 0.0}, // 3
65 : { 0.0, .375, -.125, 0.0, .75, 0.0}, // 4
66 : {-.125, 0.0, -.125, 0.5, 0.5, .25} // 5
67 : },
68 :
69 : // embedding matrix for child 2
70 : {
71 : // 0 1 2 3 4 5
72 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 0
73 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 1
74 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
75 : {-.125, -.125, 0.0, .25, 0.5, 0.5}, // 3
76 : { 0.0, -.125, .375, 0.0, .75, 0.0}, // 4
77 : {-.125, 0.0, .375, 0.0, 0.0, .75} // 5
78 : },
79 :
80 : // embedding matrix for child 3
81 : {
82 : // 0 1 2 3 4 5
83 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 0
84 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 1
85 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 2
86 : {-.125, 0.0, -.125, 0.5, 0.5, .25}, // 3
87 : {-.125, -.125, 0.0, .25, 0.5, 0.5}, // 4
88 : { 0.0, -.125, -.125, 0.5, .25, 0.5} // 5
89 : }
90 : };
91 :
92 : #endif
93 :
94 :
95 :
96 : // ------------------------------------------------------------
97 : // Tri6 class member functions
98 :
99 12586213 : bool Tri6::is_vertex(const unsigned int i) const
100 : {
101 12586213 : if (i < 3)
102 6271977 : return true;
103 545793 : return false;
104 : }
105 :
106 18864 : bool Tri6::is_edge(const unsigned int i) const
107 : {
108 18864 : if (i < 3)
109 0 : return false;
110 1700 : return true;
111 : }
112 :
113 0 : bool Tri6::is_face(const unsigned int) const
114 : {
115 0 : return false;
116 : }
117 :
118 531072 : bool Tri6::is_node_on_side(const unsigned int n,
119 : const unsigned int s) const
120 : {
121 41423 : libmesh_assert_less (s, n_sides());
122 41423 : return std::find(std::begin(side_nodes_map[s]),
123 531072 : std::end(side_nodes_map[s]),
124 531072 : n) != std::end(side_nodes_map[s]);
125 : }
126 :
127 : std::vector<unsigned>
128 24238 : Tri6::nodes_on_side(const unsigned int s) const
129 : {
130 1652 : libmesh_assert_less(s, n_sides());
131 25794 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
132 : }
133 :
134 : std::vector<unsigned>
135 1653 : Tri6::nodes_on_edge(const unsigned int e) const
136 : {
137 1653 : return nodes_on_side(e);
138 : }
139 :
140 5887019 : bool Tri6::has_affine_map() const
141 : {
142 : // Make sure edges are straight
143 1030828 : Point v = this->point(1) - this->point(0);
144 6402433 : if (!v.relative_fuzzy_equals
145 5887019 : ((this->point(3) - this->point(0))*2, affine_tol))
146 644 : return false;
147 6394061 : v = this->point(2) - this->point(1);
148 6394061 : if (!v.relative_fuzzy_equals
149 5879291 : ((this->point(4) - this->point(1))*2, affine_tol))
150 4 : return false;
151 6394009 : v = this->point(2) - this->point(0);
152 6394009 : if (!v.relative_fuzzy_equals
153 5879243 : ((this->point(5) - this->point(0))*2, affine_tol))
154 32 : return false;
155 :
156 514750 : return true;
157 : }
158 :
159 :
160 :
161 306860723 : Order Tri6::default_order() const
162 : {
163 306860723 : return SECOND;
164 : }
165 :
166 :
167 :
168 504 : dof_id_type Tri6::key (const unsigned int s) const
169 : {
170 42 : libmesh_assert_less (s, this->n_sides());
171 :
172 504 : switch (s)
173 : {
174 168 : case 0:
175 :
176 : return
177 182 : this->compute_key (this->node_id(3));
178 :
179 168 : case 1:
180 :
181 : return
182 182 : this->compute_key (this->node_id(4));
183 :
184 168 : case 2:
185 :
186 : return
187 182 : this->compute_key (this->node_id(5));
188 :
189 0 : default:
190 0 : libmesh_error_msg("Invalid side s = " << s);
191 : }
192 : }
193 :
194 :
195 :
196 72232147 : unsigned int Tri6::local_side_node(unsigned int side,
197 : unsigned int side_node) const
198 : {
199 6046926 : libmesh_assert_less (side, this->n_sides());
200 6046926 : libmesh_assert_less (side_node, Tri6::nodes_per_side);
201 :
202 72232147 : return Tri6::side_nodes_map[side][side_node];
203 : }
204 :
205 :
206 :
207 4341028 : std::unique_ptr<Elem> Tri6::build_side_ptr (const unsigned int i)
208 : {
209 4341028 : return this->simple_build_side_ptr<Edge3, Tri6>(i);
210 : }
211 :
212 :
213 :
214 16136 : void Tri6::build_side_ptr (std::unique_ptr<Elem> & side,
215 : const unsigned int i)
216 : {
217 16136 : this->simple_build_side_ptr<Tri6>(side, i, EDGE3);
218 16136 : }
219 :
220 :
221 :
222 0 : void Tri6::connectivity(const unsigned int sf,
223 : const IOPackage iop,
224 : std::vector<dof_id_type> & conn) const
225 : {
226 0 : libmesh_assert_less (sf, this->n_sub_elem());
227 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
228 :
229 0 : switch (iop)
230 : {
231 0 : case TECPLOT:
232 : {
233 0 : conn.resize(4);
234 0 : switch(sf)
235 : {
236 0 : case 0:
237 : // linear sub-triangle 0
238 0 : conn[0] = this->node_id(0)+1;
239 0 : conn[1] = this->node_id(3)+1;
240 0 : conn[2] = this->node_id(5)+1;
241 0 : conn[3] = this->node_id(5)+1;
242 :
243 0 : return;
244 :
245 0 : case 1:
246 : // linear sub-triangle 1
247 0 : conn[0] = this->node_id(3)+1;
248 0 : conn[1] = this->node_id(1)+1;
249 0 : conn[2] = this->node_id(4)+1;
250 0 : conn[3] = this->node_id(4)+1;
251 :
252 0 : return;
253 :
254 0 : case 2:
255 : // linear sub-triangle 2
256 0 : conn[0] = this->node_id(5)+1;
257 0 : conn[1] = this->node_id(4)+1;
258 0 : conn[2] = this->node_id(2)+1;
259 0 : conn[3] = this->node_id(2)+1;
260 :
261 0 : return;
262 :
263 0 : case 3:
264 : // linear sub-triangle 3
265 0 : conn[0] = this->node_id(3)+1;
266 0 : conn[1] = this->node_id(4)+1;
267 0 : conn[2] = this->node_id(5)+1;
268 0 : conn[3] = this->node_id(5)+1;
269 :
270 0 : return;
271 :
272 0 : default:
273 0 : libmesh_error_msg("Invalid sf = " << sf);
274 : }
275 : }
276 :
277 0 : case VTK:
278 : {
279 : // VTK_QUADRATIC_TRIANGLE has same numbering as libmesh TRI6
280 0 : conn.resize(Tri6::num_nodes);
281 0 : for (auto i : index_range(conn))
282 0 : conn[i] = this->node_id(i);
283 0 : return;
284 : }
285 :
286 0 : default:
287 0 : libmesh_error_msg("Unsupported IO package " << iop);
288 : }
289 : }
290 :
291 :
292 :
293 8959274 : BoundingBox Tri6::loose_bounding_box () const
294 : {
295 : // This might have curved edges, or might be a curved surface in
296 : // 3-space, in which case the full bounding box can be larger than
297 : // the bounding box of just the nodes.
298 : //
299 : //
300 : // FIXME - I haven't yet proven the formula below to be correct for
301 : // quadratics in 2D - RHS
302 357181 : Point pmin, pmax;
303 :
304 35837096 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
305 : {
306 27949365 : Real center = this->point(0)(d);
307 161266932 : for (unsigned int p=1; p != 6; ++p)
308 134389110 : center += this->point(p)(d);
309 26877822 : center /= 6;
310 :
311 27949365 : Real hd = std::abs(center - this->point(0)(d));
312 161266932 : for (unsigned int p=1; p != 6; ++p)
313 154245071 : hd = std::max(hd, std::abs(center - this->point(p)(d)));
314 :
315 26877822 : pmin(d) = center - hd;
316 26877822 : pmax(d) = center + hd;
317 : }
318 :
319 9316455 : return BoundingBox(pmin, pmax);
320 : }
321 :
322 :
323 :
324 153 : Real Tri6::volume () const
325 : {
326 : // This specialization is good for Lagrange mappings only
327 153 : if (this->mapping_type() != LAGRANGE_MAP)
328 0 : return this->Elem::volume();
329 :
330 12 : Real vol=0.;
331 :
332 : #if LIBMESH_DIM > 1
333 : // Make copies of our points. It makes the subsequent calculations a bit
334 : // shorter and avoids dereferencing the same pointer multiple times.
335 : Point
336 183 : x0 = point(0), x1 = point(1), x2 = point(2),
337 173 : x3 = point(3), x4 = point(4), x5 = point(5);
338 :
339 : // Construct constant data vectors.
340 : // \vec{x}_{\xi} = \vec{a1}*xi + \vec{b1}*eta + \vec{c1}
341 : // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}*eta + \vec{c2}
342 : Point
343 12 : a1 = 4*x0 + 4*x1 - 8*x3,
344 12 : b1 = 4*x0 - 4*x3 + 4*x4 - 4*x5, /*=a2*/
345 12 : c1 = -3*x0 - 1*x1 + 4*x3,
346 12 : b2 = 4*x0 + 4*x2 - 8*x5,
347 12 : c2 = -3*x0 - 1*x2 + 4*x5;
348 :
349 : // If a1 == b1 == a2 == b2 == 0, this is a TRI6 with straight sides,
350 : // and we can use the TRI3 formula to compute the volume.
351 153 : if (a1.relative_fuzzy_equals(Point(0,0,0)) &&
352 304 : b1.relative_fuzzy_equals(Point(0,0,0)) &&
353 22 : b2.relative_fuzzy_equals(Point(0,0,0)))
354 151 : return 0.5 * cross_norm(c1, c2);
355 :
356 : // 7-point rule, exact for quintics.
357 2 : const unsigned int N = 7;
358 :
359 : // Parameters of the quadrature rule
360 : static const Real
361 : w1 = Real(31)/480 + std::sqrt(Real(15))/2400,
362 : w2 = Real(31)/480 - std::sqrt(Real(15))/2400,
363 : q1 = Real(2)/7 + std::sqrt(Real(15))/21,
364 : q2 = Real(2)/7 - std::sqrt(Real(15))/21;
365 :
366 : static const Real xi[N] = {Real(1)/3, q1, q1, 1-2*q1, q2, q2, 1-2*q2};
367 : static const Real eta[N] = {Real(1)/3, q1, 1-2*q1, q1, q2, 1-2*q2, q2};
368 : static const Real wts[N] = {Real(9)/80, w1, w1, w1, w2, w2, w2};
369 :
370 : // Approximate the area with quadrature
371 16 : for (unsigned int q=0; q<N; ++q)
372 14 : vol += wts[q] * cross_norm(xi[q]*a1 + eta[q]*b1 + c1,
373 28 : xi[q]*b1 + eta[q]*b2 + c2);
374 : #endif // LIBMESH_DIM > 1
375 :
376 2 : return vol;
377 : }
378 :
379 :
380 :
381 12522 : unsigned short int Tri6::second_order_adjacent_vertex (const unsigned int n,
382 : const unsigned int v) const
383 : {
384 540 : libmesh_assert_greater_equal (n, this->n_vertices());
385 540 : libmesh_assert_less (n, this->n_nodes());
386 540 : libmesh_assert_less (v, 2);
387 12522 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
388 : }
389 :
390 :
391 :
392 : const unsigned short int Tri6::_second_order_adjacent_vertices[Tri6::num_sides][2] =
393 : {
394 : {0, 1}, // vertices adjacent to node 3
395 : {1, 2}, // vertices adjacent to node 4
396 : {0, 2} // vertices adjacent to node 5
397 : };
398 :
399 :
400 :
401 : std::pair<unsigned short int, unsigned short int>
402 0 : Tri6::second_order_child_vertex (const unsigned int n) const
403 : {
404 0 : libmesh_assert_greater_equal (n, this->n_vertices());
405 0 : libmesh_assert_less (n, this->n_nodes());
406 0 : return std::pair<unsigned short int, unsigned short int>
407 0 : (_second_order_vertex_child_number[n],
408 0 : _second_order_vertex_child_index[n]);
409 : }
410 :
411 :
412 :
413 : const unsigned short int Tri6::_second_order_vertex_child_number[Tri6::num_nodes] =
414 : {
415 : 99,99,99, // Vertices
416 : 0,1,0 // Edges
417 : };
418 :
419 :
420 :
421 : const unsigned short int Tri6::_second_order_vertex_child_index[Tri6::num_nodes] =
422 : {
423 : 99,99,99, // Vertices
424 : 1,2,2 // Edges
425 : };
426 :
427 :
428 158701 : void Tri6::permute(unsigned int perm_num)
429 : {
430 15288 : libmesh_assert_less (perm_num, 3);
431 :
432 312481 : for (unsigned int i = 0; i != perm_num; ++i)
433 : {
434 153780 : swap3nodes(0,1,2);
435 138918 : swap3nodes(3,4,5);
436 138918 : swap3neighbors(0,1,2);
437 : }
438 158701 : }
439 :
440 :
441 3713 : void Tri6::flip(BoundaryInfo * boundary_info)
442 : {
443 174 : libmesh_assert(boundary_info);
444 :
445 3713 : swap2nodes(0,1);
446 3713 : swap2nodes(4,5);
447 174 : swap2neighbors(1,2);
448 3713 : swap2boundarysides(1,2,boundary_info);
449 3713 : swap2boundaryedges(1,2,boundary_info);
450 3713 : }
451 :
452 :
453 288 : unsigned int Tri6::center_node_on_side(const unsigned short side) const
454 : {
455 24 : libmesh_assert_less (side, Tri6::num_sides);
456 288 : return side + 3;
457 : }
458 :
459 :
460 : ElemType
461 14844 : Tri6::side_type (const unsigned int libmesh_dbg_var(s)) const
462 : {
463 718 : libmesh_assert_less (s, 3);
464 14844 : return EDGE3;
465 : }
466 :
467 :
468 : } // namespace libMesh
|
