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 :
19 : // Local includes
20 : #include "libmesh/cell_tet10.h"
21 : #include "libmesh/edge_edge3.h"
22 : #include "libmesh/face_tri6.h"
23 : #include "libmesh/enum_io_package.h"
24 : #include "libmesh/enum_order.h"
25 :
26 : namespace libMesh
27 : {
28 :
29 :
30 :
31 : // ------------------------------------------------------------
32 : // Tet10 class static member initializations
33 : const int Tet10::num_nodes;
34 : const int Tet10::nodes_per_side;
35 : const int Tet10::nodes_per_edge;
36 :
37 : // ------------------------------------------------------------
38 : // Tet10 class member functions
39 :
40 20169103 : bool Tet10::is_vertex(const unsigned int i) const
41 : {
42 20169103 : if (i < 4)
43 8267948 : return true;
44 764732 : return false;
45 : }
46 :
47 1202688 : bool Tet10::is_edge(const unsigned int i) const
48 : {
49 1202688 : if (i < 4)
50 0 : return false;
51 100224 : return true;
52 : }
53 :
54 0 : bool Tet10::is_face(const unsigned int) const
55 : {
56 0 : return false;
57 : }
58 :
59 183704 : bool Tet10::is_node_on_side(const unsigned int n,
60 : const unsigned int s) const
61 : {
62 15152 : libmesh_assert_less (s, n_sides());
63 15152 : return std::find(std::begin(side_nodes_map[s]),
64 183704 : std::end(side_nodes_map[s]),
65 183704 : n) != std::end(side_nodes_map[s]);
66 : }
67 :
68 : std::vector<unsigned>
69 16269463 : Tet10::nodes_on_side(const unsigned int s) const
70 : {
71 761036 : libmesh_assert_less(s, n_sides());
72 17030499 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
73 : }
74 :
75 : std::vector<unsigned>
76 180138 : Tet10::nodes_on_edge(const unsigned int e) const
77 : {
78 14988 : libmesh_assert_less(e, n_edges());
79 195126 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
80 : }
81 :
82 4622916 : bool Tet10::is_node_on_edge(const unsigned int n,
83 : const unsigned int e) const
84 : {
85 168408 : libmesh_assert_less (e, n_edges());
86 168408 : return std::find(std::begin(edge_nodes_map[e]),
87 4622916 : std::end(edge_nodes_map[e]),
88 4622916 : n) != std::end(edge_nodes_map[e]);
89 : }
90 :
91 :
92 : #ifdef LIBMESH_ENABLE_AMR
93 :
94 : // This function only works if LIBMESH_ENABLE_AMR...
95 11473615 : bool Tet10::is_child_on_side(const unsigned int c,
96 : const unsigned int s) const
97 : {
98 : // Table of local IDs for the midege nodes on the side opposite a given node.
99 : // See the ASCII art in the header file for this class to confirm this.
100 11473615 : const unsigned int midedge_nodes_opposite[4][3] =
101 : {
102 : {5,8,9}, // midedge nodes opposite node 0
103 : {6,7,9}, // midedge nodes opposite node 1
104 : {4,7,8}, // midedge nodes opposite node 2
105 : {4,5,6} // midedge nodes opposite node 3
106 : };
107 :
108 : // Call the base class helper function
109 11926159 : return Tet::is_child_on_side_helper(c, s, midedge_nodes_opposite);
110 : }
111 :
112 : #else
113 :
114 : bool Tet10::is_child_on_side(const unsigned int /*c*/,
115 : const unsigned int /*s*/) const
116 : {
117 : libmesh_not_implemented();
118 : return false;
119 : }
120 :
121 : #endif //LIBMESH_ENABLE_AMR
122 :
123 :
124 :
125 1387764 : bool Tet10::has_affine_map() const
126 : {
127 : // Make sure edges are straight
128 221968 : Point v = this->point(1) - this->point(0);
129 1498748 : if (!v.relative_fuzzy_equals
130 1387764 : ((this->point(4) - this->point(0))*2, affine_tol))
131 4752 : return false;
132 1436972 : v = this->point(2) - this->point(1);
133 1436972 : if (!v.relative_fuzzy_equals
134 1330740 : ((this->point(5) - this->point(1))*2, affine_tol))
135 96 : return false;
136 1435724 : v = this->point(2) - this->point(0);
137 1435724 : if (!v.relative_fuzzy_equals
138 1329588 : ((this->point(6) - this->point(0))*2, affine_tol))
139 0 : return false;
140 1435724 : v = this->point(3) - this->point(0);
141 1435724 : if (!v.relative_fuzzy_equals
142 1329588 : ((this->point(7) - this->point(0))*2, affine_tol))
143 144 : return false;
144 1433852 : v = this->point(3) - this->point(1);
145 1433852 : if (!v.relative_fuzzy_equals
146 1327860 : ((this->point(8) - this->point(1))*2, affine_tol))
147 0 : return false;
148 1433852 : v = this->point(3) - this->point(2);
149 1433852 : if (!v.relative_fuzzy_equals
150 1327860 : ((this->point(9) - this->point(2))*2, affine_tol))
151 0 : return false;
152 105992 : return true;
153 : }
154 :
155 :
156 :
157 179304460 : Order Tet10::default_order() const
158 : {
159 179304460 : return SECOND;
160 : }
161 :
162 :
163 :
164 44020487 : unsigned int Tet10::local_side_node(unsigned int side,
165 : unsigned int side_node) const
166 : {
167 1883678 : libmesh_assert_less (side, this->n_sides());
168 1883678 : libmesh_assert_less (side_node, Tet10::nodes_per_side);
169 :
170 44020487 : return Tet10::side_nodes_map[side][side_node];
171 : }
172 :
173 :
174 :
175 88194489 : unsigned int Tet10::local_edge_node(unsigned int edge,
176 : unsigned int edge_node) const
177 : {
178 6930906 : libmesh_assert_less (edge, this->n_edges());
179 6930906 : libmesh_assert_less (edge_node, Tet10::nodes_per_edge);
180 :
181 88194489 : return Tet10::edge_nodes_map[edge][edge_node];
182 : }
183 :
184 :
185 :
186 3124496 : std::unique_ptr<Elem> Tet10::build_side_ptr (const unsigned int i)
187 : {
188 3124496 : return this->simple_build_side_ptr<Tri6, Tet10>(i);
189 : }
190 :
191 :
192 :
193 5716976 : void Tet10::build_side_ptr (std::unique_ptr<Elem> & side,
194 : const unsigned int i)
195 : {
196 5716976 : this->simple_build_side_ptr<Tet10>(side, i, TRI6);
197 5716976 : }
198 :
199 :
200 :
201 5246393 : std::unique_ptr<Elem> Tet10::build_edge_ptr (const unsigned int i)
202 : {
203 5246393 : return this->simple_build_edge_ptr<Edge3,Tet10>(i);
204 : }
205 :
206 :
207 :
208 15918 : void Tet10::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
209 : {
210 15918 : this->simple_build_edge_ptr<Tet10>(edge, i, EDGE3);
211 15918 : }
212 :
213 :
214 :
215 0 : void Tet10::connectivity(const unsigned int sc,
216 : const IOPackage iop,
217 : std::vector<dof_id_type> & conn) const
218 : {
219 0 : libmesh_assert(_nodes);
220 0 : libmesh_assert_less (sc, this->n_sub_elem());
221 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
222 :
223 0 : switch (iop)
224 : {
225 0 : case TECPLOT:
226 : {
227 0 : conn.resize(8);
228 0 : switch (sc)
229 : {
230 :
231 :
232 : // Linear sub-tet 0
233 0 : case 0:
234 :
235 0 : conn[0] = this->node_id(0)+1;
236 0 : conn[1] = this->node_id(4)+1;
237 0 : conn[2] = this->node_id(6)+1;
238 0 : conn[3] = this->node_id(6)+1;
239 0 : conn[4] = this->node_id(7)+1;
240 0 : conn[5] = this->node_id(7)+1;
241 0 : conn[6] = this->node_id(7)+1;
242 0 : conn[7] = this->node_id(7)+1;
243 :
244 0 : return;
245 :
246 : // Linear sub-tet 1
247 0 : case 1:
248 :
249 0 : conn[0] = this->node_id(4)+1;
250 0 : conn[1] = this->node_id(1)+1;
251 0 : conn[2] = this->node_id(5)+1;
252 0 : conn[3] = this->node_id(5)+1;
253 0 : conn[4] = this->node_id(8)+1;
254 0 : conn[5] = this->node_id(8)+1;
255 0 : conn[6] = this->node_id(8)+1;
256 0 : conn[7] = this->node_id(8)+1;
257 :
258 0 : return;
259 :
260 : // Linear sub-tet 2
261 0 : case 2:
262 :
263 0 : conn[0] = this->node_id(5)+1;
264 0 : conn[1] = this->node_id(2)+1;
265 0 : conn[2] = this->node_id(6)+1;
266 0 : conn[3] = this->node_id(6)+1;
267 0 : conn[4] = this->node_id(9)+1;
268 0 : conn[5] = this->node_id(9)+1;
269 0 : conn[6] = this->node_id(9)+1;
270 0 : conn[7] = this->node_id(9)+1;
271 :
272 0 : return;
273 :
274 : // Linear sub-tet 3
275 0 : case 3:
276 :
277 0 : conn[0] = this->node_id(7)+1;
278 0 : conn[1] = this->node_id(8)+1;
279 0 : conn[2] = this->node_id(9)+1;
280 0 : conn[3] = this->node_id(9)+1;
281 0 : conn[4] = this->node_id(3)+1;
282 0 : conn[5] = this->node_id(3)+1;
283 0 : conn[6] = this->node_id(3)+1;
284 0 : conn[7] = this->node_id(3)+1;
285 :
286 0 : return;
287 :
288 : // Linear sub-tet 4
289 0 : case 4:
290 :
291 0 : conn[0] = this->node_id(4)+1;
292 0 : conn[1] = this->node_id(8)+1;
293 0 : conn[2] = this->node_id(6)+1;
294 0 : conn[3] = this->node_id(6)+1;
295 0 : conn[4] = this->node_id(7)+1;
296 0 : conn[5] = this->node_id(7)+1;
297 0 : conn[6] = this->node_id(7)+1;
298 0 : conn[7] = this->node_id(7)+1;
299 :
300 0 : return;
301 :
302 : // Linear sub-tet 5
303 0 : case 5:
304 :
305 0 : conn[0] = this->node_id(4)+1;
306 0 : conn[1] = this->node_id(5)+1;
307 0 : conn[2] = this->node_id(6)+1;
308 0 : conn[3] = this->node_id(6)+1;
309 0 : conn[4] = this->node_id(8)+1;
310 0 : conn[5] = this->node_id(8)+1;
311 0 : conn[6] = this->node_id(8)+1;
312 0 : conn[7] = this->node_id(8)+1;
313 :
314 0 : return;
315 :
316 : // Linear sub-tet 6
317 0 : case 6:
318 :
319 0 : conn[0] = this->node_id(5)+1;
320 0 : conn[1] = this->node_id(9)+1;
321 0 : conn[2] = this->node_id(6)+1;
322 0 : conn[3] = this->node_id(6)+1;
323 0 : conn[4] = this->node_id(8)+1;
324 0 : conn[5] = this->node_id(8)+1;
325 0 : conn[6] = this->node_id(8)+1;
326 0 : conn[7] = this->node_id(8)+1;
327 :
328 0 : return;
329 :
330 : // Linear sub-tet 7
331 0 : case 7:
332 :
333 0 : conn[0] = this->node_id(7)+1;
334 0 : conn[1] = this->node_id(6)+1;
335 0 : conn[2] = this->node_id(9)+1;
336 0 : conn[3] = this->node_id(9)+1;
337 0 : conn[4] = this->node_id(8)+1;
338 0 : conn[5] = this->node_id(8)+1;
339 0 : conn[6] = this->node_id(8)+1;
340 0 : conn[7] = this->node_id(8)+1;
341 :
342 0 : return;
343 :
344 :
345 0 : default:
346 0 : libmesh_error_msg("Invalid sc = " << sc);
347 : }
348 : }
349 :
350 0 : case VTK:
351 : {
352 : // VTK connectivity for VTK_QUADRATIC_TETRA matches libMesh's own.
353 0 : conn.resize(Tet10::num_nodes);
354 0 : for (auto i : index_range(conn))
355 0 : conn[i] = this->node_id(i);
356 0 : return;
357 : }
358 :
359 0 : default:
360 0 : libmesh_error_msg("Unsupported IO package " << iop);
361 : }
362 : }
363 :
364 :
365 :
366 : const unsigned short int Tet10::_second_order_vertex_child_number[10] =
367 : {
368 : 99,99,99,99, // Vertices
369 : 0,1,0,0,1,2 // Edges
370 : };
371 :
372 :
373 :
374 : const unsigned short int Tet10::_second_order_vertex_child_index[10] =
375 : {
376 : 99,99,99,99, // Vertices
377 : 1,2,2,3,3,3 // Edges
378 : };
379 :
380 :
381 :
382 : std::pair<unsigned short int, unsigned short int>
383 0 : Tet10::second_order_child_vertex (const unsigned int n) const
384 : {
385 0 : libmesh_assert_greater_equal (n, this->n_vertices());
386 0 : libmesh_assert_less (n, this->n_nodes());
387 0 : return std::pair<unsigned short int, unsigned short int>
388 0 : (_second_order_vertex_child_number[n],
389 0 : _second_order_vertex_child_index[n]);
390 : }
391 :
392 :
393 :
394 10744320 : unsigned short int Tet10::second_order_adjacent_vertex (const unsigned int n,
395 : const unsigned int v) const
396 : {
397 318720 : libmesh_assert_greater_equal (n, this->n_vertices());
398 318720 : libmesh_assert_less (n, this->n_nodes());
399 318720 : libmesh_assert_less (v, 2);
400 10744320 : return _second_order_adjacent_vertices[n-this->n_vertices()][v];
401 : }
402 :
403 :
404 :
405 : const unsigned short int Tet10::_second_order_adjacent_vertices[6][2] =
406 : {
407 : {0, 1}, // vertices adjacent to node 4
408 : {1, 2}, // vertices adjacent to node 5
409 : {0, 2}, // vertices adjacent to node 6
410 : {0, 3}, // vertices adjacent to node 7
411 : {1, 3}, // vertices adjacent to node 8
412 : {2, 3} // vertices adjacent to node 9
413 : };
414 :
415 :
416 :
417 :
418 :
419 : #ifdef LIBMESH_ENABLE_AMR
420 :
421 : const Real Tet10::_embedding_matrix[Tet10::num_children][Tet10::num_nodes][Tet10::num_nodes] =
422 : {
423 : // embedding matrix for child 0
424 : {
425 : // 0 1 2 3 4 5 6 7 8 9
426 : { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
427 : { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 1
428 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
429 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
430 : { 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
431 : { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 5
432 : { 0.375, 0.,-0.125, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
433 : { 0.375, 0., 0.,-0.125, 0., 0., 0., 0.75, 0., 0.}, // 7
434 : { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 8
435 : { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
436 : },
437 :
438 : // embedding matrix for child 1
439 : {
440 : // 0 1 2 3 4 5 6 7 8 9
441 : { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
442 : { 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
443 : { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 2
444 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
445 : {-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
446 : { 0., 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
447 : {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 6
448 : {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
449 : { 0., 0.375, 0.,-0.125, 0., 0., 0., 0., 0.75, 0.}, // 8
450 : { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25} // 9
451 : },
452 :
453 : // embedding matrix for child 2
454 : {
455 : // 0 1 2 3 4 5 6 7 8 9
456 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
457 : { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
458 : { 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 2
459 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 3
460 : {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
461 : { 0.,-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
462 : {-0.125, 0., 0.375, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
463 : {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 7
464 : { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 8
465 : { 0., 0., 0.375,-0.125, 0., 0., 0., 0., 0., 0.75} // 9
466 : },
467 :
468 : // embedding matrix for child 3
469 : {
470 : // 0 1 2 3 4 5 6 7 8 9
471 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
472 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
473 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
474 : { 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 3
475 : {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 4
476 : { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
477 : {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5}, // 6
478 : {-0.125, 0., 0., 0.375, 0., 0., 0., 0.75, 0., 0.}, // 7
479 : { 0.,-0.125, 0., 0.375, 0., 0., 0., 0., 0.75, 0.}, // 8
480 : { 0., 0.,-0.125, 0.375, 0., 0., 0., 0., 0., 0.75} // 9
481 : },
482 :
483 : // embedding matrix for child 4
484 : {
485 : // 0 1 2 3 4 5 6 7 8 9
486 : { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
487 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
488 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
489 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
490 : {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 4
491 : {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 5
492 : { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
493 : { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 7
494 : {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
495 : { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
496 : },
497 :
498 : // embedding matrix for child 5
499 : {
500 : // 0 1 2 3 4 5 6 7 8 9
501 : { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
502 : { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
503 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
504 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
505 : {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 4
506 : {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 5
507 : { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
508 : {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
509 : { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
510 : {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25} // 9
511 : },
512 :
513 : // embedding matrix for child 6
514 : {
515 : // 0 1 2 3 4 5 6 7 8 9
516 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
517 : { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
518 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
519 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
520 : {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
521 : { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 5
522 : {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
523 : {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 7
524 : { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
525 : { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5} // 9
526 : },
527 :
528 : // embedding matrix for child 7
529 : {
530 : // 0 1 2 3 4 5 6 7 8 9
531 : { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
532 : { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
533 : { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
534 : { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
535 : {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 4
536 : { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
537 : {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
538 : { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25}, // 7
539 : {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
540 : {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5} // 9
541 : }
542 : };
543 :
544 :
545 :
546 26690145 : Real Tet10::embedding_matrix (const unsigned int i,
547 : const unsigned int j,
548 : const unsigned int k) const
549 : {
550 : // Choose an optimal diagonal, if one has not already been selected
551 26690145 : this->choose_diagonal();
552 :
553 : // Permuted j and k indices
554 : unsigned int
555 1094408 : jp=j,
556 1094408 : kp=k;
557 :
558 26690145 : if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
559 : {
560 : // Just the enum value cast to an unsigned int...
561 177840 : const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection); // == 1 or 2
562 :
563 : // Instead of doing a lot of arithmetic, use these
564 : // straightforward arrays for the permutations. Note that 3 ->
565 : // 3, and the first array consists of "forward" permutations of
566 : // the sets {0,1,2}, {4,5,6}, and {7,8,9} while the second array
567 : // consists of "reverse" permutations of the same sets.
568 4427598 : const unsigned int perms[2][10] =
569 : {
570 : {1, 2, 0, 3, 5, 6, 4, 8, 9, 7},
571 : {2, 0, 1, 3, 6, 4, 5, 9, 7, 8}
572 : };
573 :
574 : // Permute j
575 4427598 : jp = perms[ds-1][j];
576 : // if (jp<3)
577 : // jp = (jp+ds)%3;
578 : // else if (jp>3)
579 : // jp = (jp-1+ds)%3 + 1 + 3*((jp-1)/3);
580 :
581 : // Permute k
582 4427598 : kp = perms[ds-1][k];
583 : // if (kp<3)
584 : // kp = (kp+ds)%3;
585 : // else if (kp>3)
586 : // kp = (kp-1+ds)%3 + 1 + 3*((kp-1)/3);
587 : }
588 :
589 : // Debugging:
590 : // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
591 : // libMesh::err << "j=" << j << std::endl;
592 : // libMesh::err << "k=" << k << std::endl;
593 : // libMesh::err << "jp=" << jp << std::endl;
594 : // libMesh::err << "kp=" << kp << std::endl;
595 :
596 : // Call embedding matrix with permuted indices
597 26690145 : return this->_embedding_matrix[i][jp][kp];
598 : }
599 :
600 : #endif // #ifdef LIBMESH_ENABLE_AMR
601 :
602 :
603 :
604 140448 : Real Tet10::volume () const
605 : {
606 : // This specialization is good for Lagrange mappings only in general
607 140448 : if (this->mapping_type() != LAGRANGE_MAP)
608 2653 : return this->Elem::volume();
609 :
610 : // Make copies of our points. It makes the subsequent calculations a bit
611 : // shorter and avoids dereferencing the same pointer multiple times.
612 : Point
613 192945 : x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
614 181915 : x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9);
615 :
616 : // The constant components of the dx/dxi vector, linear in xi, eta, zeta.
617 : // These were copied directly from the output of a Python script.
618 : Point dx_dxi[4] =
619 : {
620 16408 : -3*x0 - x1 + 4*x4, // constant
621 16408 : 4*x0 - 4*x4 - 4*x7 + 4*x8, // zeta
622 16408 : 4*x0 - 4*x4 + 4*x5 - 4*x6, // eta
623 16408 : 4*x0 + 4*x1 - 8*x4 // xi
624 82040 : };
625 :
626 : // The constant components of the dx/deta vector, linear in xi, eta, zeta.
627 : // These were copied directly from the output of a Python script.
628 : Point dx_deta[4] =
629 : {
630 16408 : -3*x0 - x2 + 4*x6, // constant
631 16408 : 4*x0 - 4*x6 - 4*x7 + 4*x9, // zeta
632 16408 : 4*x0 + 4*x2 - 8*x6, // eta
633 16408 : 4*x0 - 4*x4 + 4*x5 - 4*x6 // xi
634 82040 : };
635 :
636 : // The constant components of the dx/dzeta vector, linear in xi, eta, zeta.
637 : // These were copied directly from the output of a Python script.
638 : Point dx_dzeta[4] =
639 : {
640 16408 : -3*x0 - x3 + 4*x7, // constant
641 16408 : 4*x0 + 4*x3 - 8*x7, // zeta
642 16408 : 4*x0 - 4*x6 - 4*x7 + 4*x9, // eta
643 16408 : 4*x0 - 4*x4 - 4*x7 + 4*x8 // xi
644 82040 : };
645 :
646 : // 2x2x2 conical quadrature rule. Note: there is also a five point
647 : // rule for tets with a negative weight which would be cheaper, but
648 : // we'll use this one to preclude any possible issues with
649 : // cancellation error.
650 16408 : const int N = 8;
651 : static const Real w[N] =
652 : {
653 : 3.6979856358852914509238091810505e-02_R,
654 : 1.6027040598476613723156741868689e-02_R,
655 : 2.1157006454524061178256145400082e-02_R,
656 : 9.1694299214797439226823542540576e-03_R,
657 : 3.6979856358852914509238091810505e-02_R,
658 : 1.6027040598476613723156741868689e-02_R,
659 : 2.1157006454524061178256145400082e-02_R,
660 : 9.1694299214797439226823542540576e-03_R
661 : };
662 :
663 : static const Real xi[N] =
664 : {
665 : 1.2251482265544137786674043037115e-01_R,
666 : 5.4415184401122528879992623629551e-01_R,
667 : 1.2251482265544137786674043037115e-01_R,
668 : 5.4415184401122528879992623629551e-01_R,
669 : 1.2251482265544137786674043037115e-01_R,
670 : 5.4415184401122528879992623629551e-01_R,
671 : 1.2251482265544137786674043037115e-01_R,
672 : 5.4415184401122528879992623629551e-01_R
673 : };
674 :
675 : static const Real eta[N] =
676 : {
677 : 1.3605497680284601717109468420738e-01_R,
678 : 7.0679724159396903069267439165167e-02_R,
679 : 5.6593316507280088053551297149570e-01_R,
680 : 2.9399880063162286589079157179842e-01_R,
681 : 1.3605497680284601717109468420738e-01_R,
682 : 7.0679724159396903069267439165167e-02_R,
683 : 5.6593316507280088053551297149570e-01_R,
684 : 2.9399880063162286589079157179842e-01_R
685 : };
686 :
687 : static const Real zeta[N] =
688 : {
689 : 1.5668263733681830907933725249176e-01_R,
690 : 8.1395667014670255076709592007207e-02_R,
691 : 6.5838687060044409936029672711329e-02_R,
692 : 3.4202793236766414300604458388142e-02_R,
693 : 5.8474756320489429588282763292971e-01_R,
694 : 3.0377276481470755305409673253211e-01_R,
695 : 2.4571332521171333166171692542182e-01_R,
696 : 1.2764656212038543100867773351792e-01_R
697 : };
698 :
699 16408 : Real vol = 0.;
700 1240155 : for (int q=0; q<N; ++q)
701 : {
702 : // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
703 : Point
704 395984 : dx_dxi_q = dx_dxi[0] + zeta[q]*dx_dxi[1] + eta[q]*dx_dxi[2] + xi[q]*dx_dxi[3],
705 131264 : dx_deta_q = dx_deta[0] + zeta[q]*dx_deta[1] + eta[q]*dx_deta[2] + xi[q]*dx_deta[3],
706 131264 : dx_dzeta_q = dx_dzeta[0] + zeta[q]*dx_dzeta[1] + eta[q]*dx_dzeta[2] + xi[q]*dx_dzeta[3];
707 :
708 : // Compute scalar triple product, multiply by weight, and accumulate volume.
709 1190600 : vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
710 : }
711 :
712 16408 : return vol;
713 : }
714 :
715 :
716 278784 : void Tet10::permute(unsigned int perm_num)
717 : {
718 12672 : libmesh_assert_less (perm_num, 12);
719 :
720 278784 : const unsigned int side = perm_num % 4;
721 278784 : const unsigned int rotate = perm_num / 4;
722 :
723 573888 : for (unsigned int i = 0; i != rotate; ++i)
724 : {
725 295104 : swap3nodes(0,1,2);
726 281952 : swap3nodes(4,5,6);
727 281952 : swap3nodes(7,8,9);
728 281952 : swap3neighbors(1,2,3);
729 : }
730 :
731 278784 : switch (side) {
732 3648 : case 0:
733 3648 : break;
734 61536 : case 1:
735 61536 : swap3nodes(0,2,3);
736 58608 : swap3nodes(4,5,8);
737 58608 : swap3nodes(6,9,7);
738 58608 : swap3neighbors(0,2,1);
739 58608 : break;
740 53376 : case 2:
741 53376 : swap3nodes(2,0,3);
742 50688 : swap3nodes(5,4,8);
743 50688 : swap3nodes(6,7,9);
744 50688 : swap3neighbors(0,1,2);
745 50688 : break;
746 77856 : case 3:
747 77856 : swap3nodes(2,1,3);
748 74448 : swap3nodes(5,8,9);
749 74448 : swap3nodes(6,4,7);
750 74448 : swap3neighbors(0,1,3);
751 74448 : break;
752 0 : default:
753 0 : libmesh_error();
754 : }
755 278784 : }
756 :
757 :
758 6912 : void Tet10::flip(BoundaryInfo * boundary_info)
759 : {
760 576 : libmesh_assert(boundary_info);
761 :
762 6912 : swap2nodes(0,2);
763 6912 : swap2nodes(4,5);
764 6912 : swap2nodes(7,9);
765 576 : swap2neighbors(1,2);
766 6912 : swap2boundarysides(1,2,boundary_info);
767 6912 : swap2boundaryedges(0,1,boundary_info);
768 6912 : swap2boundaryedges(3,5,boundary_info);
769 6912 : }
770 :
771 :
772 30078532 : ElemType Tet10::side_type (const unsigned int libmesh_dbg_var(s)) const
773 : {
774 1316904 : libmesh_assert_less (s, 4);
775 30078532 : return TRI6;
776 : }
777 :
778 :
779 : } // namespace libMesh
|
