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