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