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