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