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_tet4.h"
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/face_tri3.h"
23 : #include "libmesh/tensor_value.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 :
32 : // ------------------------------------------------------------
33 : // Tet4 class static member initializations
34 : const int Tet4::num_nodes;
35 : const int Tet4::nodes_per_side;
36 : const int Tet4::nodes_per_edge;
37 :
38 : const unsigned int Tet4::side_nodes_map[Tet4::num_sides][Tet4::nodes_per_side] =
39 : {
40 : {0, 2, 1}, // Side 0
41 : {0, 1, 3}, // Side 1
42 : {1, 2, 3}, // Side 2
43 : {2, 0, 3} // Side 3
44 : };
45 :
46 : const unsigned int Tet4::edge_nodes_map[Tet4::num_edges][Tet4::nodes_per_edge] =
47 : {
48 : {0, 1}, // Edge 0
49 : {1, 2}, // Edge 1
50 : {0, 2}, // Edge 2
51 : {0, 3}, // Edge 3
52 : {1, 3}, // Edge 4
53 : {2, 3} // Edge 5
54 : };
55 :
56 : // ------------------------------------------------------------
57 : // Tet4 class member functions
58 :
59 15034297 : bool Tet4::is_vertex(const unsigned int libmesh_dbg_var(n)) const
60 : {
61 702873 : libmesh_assert_not_equal_to (n, invalid_uint);
62 15034297 : return true;
63 : }
64 :
65 0 : bool Tet4::is_edge(const unsigned int) const
66 : {
67 0 : return false;
68 : }
69 :
70 0 : bool Tet4::is_face(const unsigned int) const
71 : {
72 0 : return false;
73 : }
74 :
75 29352 : bool Tet4::is_node_on_edge(const unsigned int n,
76 : const unsigned int e) const
77 : {
78 2352 : libmesh_assert_less (e, n_edges());
79 2352 : return std::find(std::begin(edge_nodes_map[e]),
80 29352 : std::end(edge_nodes_map[e]),
81 29352 : n) != std::end(edge_nodes_map[e]);
82 : }
83 :
84 :
85 :
86 :
87 : #ifdef LIBMESH_ENABLE_AMR
88 :
89 : // This function only works if LIBMESH_ENABLE_AMR...
90 1862556 : bool Tet4::is_child_on_side(const unsigned int c,
91 : const unsigned int s) const
92 : {
93 : // OK, for the Tet4, this is pretty obvious... it is sets of nodes
94 : // not equal to the current node. But if we want this algorithm to
95 : // be generic and work for Tet10 also it helps to do it this way.
96 1862556 : const unsigned int nodes_opposite[4][3] =
97 : {
98 : {1,2,3}, // nodes opposite node 0
99 : {0,2,3}, // nodes opposite node 1
100 : {0,1,3}, // nodes opposite node 2
101 : {0,1,2} // nodes opposite node 3
102 : };
103 :
104 : // Call the base class helper function
105 2572892 : return Tet::is_child_on_side_helper(c, s, nodes_opposite);
106 : }
107 :
108 : #else
109 :
110 : bool Tet4::is_child_on_side(const unsigned int /*c*/,
111 : const unsigned int /*s*/) const
112 : {
113 : libmesh_not_implemented();
114 : return false;
115 : }
116 :
117 : #endif //LIBMESH_ENABLE_AMR
118 :
119 :
120 :
121 27648 : bool Tet4::has_invertible_map(Real tol) const
122 : {
123 27648 : return this->volume() > tol;
124 : }
125 :
126 :
127 :
128 130016 : bool Tet4::is_node_on_side(const unsigned int n,
129 : const unsigned int s) const
130 : {
131 8792 : libmesh_assert_less (s, n_sides());
132 8792 : return std::find(std::begin(side_nodes_map[s]),
133 130016 : std::end(side_nodes_map[s]),
134 130016 : n) != std::end(side_nodes_map[s]);
135 : }
136 :
137 : std::vector<unsigned>
138 176614 : Tet4::nodes_on_side(const unsigned int s) const
139 : {
140 13600 : libmesh_assert_less(s, n_sides());
141 190214 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
142 : }
143 :
144 : std::vector<unsigned>
145 69546 : Tet4::nodes_on_edge(const unsigned int e) const
146 : {
147 5772 : libmesh_assert_less(e, n_edges());
148 75318 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
149 : }
150 :
151 33258494 : Order Tet4::default_order() const
152 : {
153 33258494 : return FIRST;
154 : }
155 :
156 4344370 : std::unique_ptr<Elem> Tet4::build_side_ptr (const unsigned int i)
157 : {
158 4344370 : return this->simple_build_side_ptr<Tri3, Tet4>(i);
159 : }
160 :
161 :
162 :
163 9707604 : void Tet4::build_side_ptr (std::unique_ptr<Elem> & side,
164 : const unsigned int i)
165 : {
166 9707604 : this->simple_build_side_ptr<Tet4>(side, i, TRI3);
167 9707604 : }
168 :
169 :
170 :
171 5721818 : std::unique_ptr<Elem> Tet4::build_edge_ptr (const unsigned int i)
172 : {
173 5721818 : return this->simple_build_edge_ptr<Edge2,Tet4>(i);
174 : }
175 :
176 :
177 :
178 0 : void Tet4::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
179 : {
180 0 : this->simple_build_edge_ptr<Tet4>(edge, i, EDGE2);
181 0 : }
182 :
183 :
184 :
185 3474 : void Tet4::connectivity(const unsigned int libmesh_dbg_var(sc),
186 : const IOPackage iop,
187 : std::vector<dof_id_type> & conn) const
188 : {
189 1737 : libmesh_assert(_nodes);
190 1737 : libmesh_assert_less (sc, this->n_sub_elem());
191 1737 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
192 :
193 :
194 3474 : switch (iop)
195 : {
196 3474 : case TECPLOT:
197 : {
198 3474 : conn.resize(8);
199 5211 : conn[0] = this->node_id(0)+1;
200 3474 : conn[1] = this->node_id(1)+1;
201 3474 : conn[2] = this->node_id(2)+1;
202 3474 : conn[3] = this->node_id(2)+1;
203 3474 : conn[4] = this->node_id(3)+1;
204 3474 : conn[5] = this->node_id(3)+1;
205 3474 : conn[6] = this->node_id(3)+1;
206 3474 : conn[7] = this->node_id(3)+1;
207 3474 : return;
208 : }
209 :
210 0 : case VTK:
211 : {
212 0 : conn.resize(4);
213 0 : conn[0] = this->node_id(0);
214 0 : conn[1] = this->node_id(1);
215 0 : conn[2] = this->node_id(2);
216 0 : conn[3] = this->node_id(3);
217 0 : return;
218 : }
219 :
220 0 : default:
221 0 : libmesh_error_msg("Unsupported IO package " << iop);
222 : }
223 : }
224 :
225 :
226 :
227 : #ifdef LIBMESH_ENABLE_AMR
228 :
229 : const Real Tet4::_embedding_matrix[Tet4::num_children][Tet4::num_nodes][Tet4::num_nodes] =
230 : {
231 : // embedding matrix for child 0
232 : {
233 : // 0 1 2 3
234 : {1.0, 0.0, 0.0, 0.0}, // 0
235 : {0.5, 0.5, 0.0, 0.0}, // 1
236 : {0.5, 0.0, 0.5, 0.0}, // 2
237 : {0.5, 0.0, 0.0, 0.5} // 3
238 : },
239 :
240 : // embedding matrix for child 1
241 : {
242 : // 0 1 2 3
243 : {0.5, 0.5, 0.0, 0.0}, // 0
244 : {0.0, 1.0, 0.0, 0.0}, // 1
245 : {0.0, 0.5, 0.5, 0.0}, // 2
246 : {0.0, 0.5, 0.0, 0.5} // 3
247 : },
248 :
249 : // embedding matrix for child 2
250 : {
251 : // 0 1 2 3
252 : {0.5, 0.0, 0.5, 0.0}, // 0
253 : {0.0, 0.5, 0.5, 0.0}, // 1
254 : {0.0, 0.0, 1.0, 0.0}, // 2
255 : {0.0, 0.0, 0.5, 0.5} // 3
256 : },
257 :
258 : // embedding matrix for child 3
259 : {
260 : // 0 1 2 3
261 : {0.5, 0.0, 0.0, 0.5}, // 0
262 : {0.0, 0.5, 0.0, 0.5}, // 1
263 : {0.0, 0.0, 0.5, 0.5}, // 2
264 : {0.0, 0.0, 0.0, 1.0} // 3
265 : },
266 :
267 : // embedding matrix for child 4
268 : {
269 : // 0 1 2 3
270 : {0.5, 0.5, 0.0, 0.0}, // 0
271 : {0.0, 0.5, 0.0, 0.5}, // 1
272 : {0.5, 0.0, 0.5, 0.0}, // 2
273 : {0.5, 0.0, 0.0, 0.5} // 3
274 : },
275 :
276 : // embedding matrix for child 5
277 : {
278 : // 0 1 2 3
279 : {0.5, 0.5, 0.0, 0.0}, // 0
280 : {0.0, 0.5, 0.5, 0.0}, // 1
281 : {0.5, 0.0, 0.5, 0.0}, // 2
282 : {0.0, 0.5, 0.0, 0.5} // 3
283 : },
284 :
285 : // embedding matrix for child 6
286 : {
287 : // 0 1 2 3
288 : {0.5, 0.0, 0.5, 0.0}, // 0
289 : {0.0, 0.5, 0.5, 0.0}, // 1
290 : {0.0, 0.0, 0.5, 0.5}, // 2
291 : {0.0, 0.5, 0.0, 0.5} // 3
292 : },
293 :
294 : // embedding matrix for child 7
295 : {
296 : // 0 1 2 3
297 : {0.5, 0.0, 0.5, 0.0}, // 0
298 : {0.0, 0.5, 0.0, 0.5}, // 1
299 : {0.0, 0.0, 0.5, 0.5}, // 2
300 : {0.5, 0.0, 0.0, 0.5} // 3
301 : }
302 : };
303 :
304 : #endif // #ifdef LIBMESH_ENABLE_AMR
305 :
306 :
307 :
308 :
309 :
310 29952 : Point Tet4::true_centroid () const
311 : {
312 29952 : return Elem::vertex_average();
313 : }
314 :
315 :
316 :
317 123779 : Real Tet4::volume () const
318 : {
319 : // The volume of a tetrahedron is 1/6 the box product formed
320 : // by its base and apex vectors
321 24920 : Point a = point(3) - point(0);
322 :
323 : // b is the vector pointing from 0 to 1
324 12460 : Point b = point(1) - point(0);
325 :
326 : // c is the vector pointing from 0 to 2
327 12460 : Point c = point(2) - point(0);
328 :
329 123779 : return triple_product(a, b, c) / 6.;
330 : }
331 :
332 :
333 :
334 :
335 0 : std::pair<Real, Real> Tet4::min_and_max_angle() const
336 : {
337 0 : Point n[4];
338 :
339 : // Compute the outward normal vectors on each face
340 0 : n[0] = (this->point(2) - this->point(0)).cross(this->point(1) - this->point(0));
341 0 : n[1] = (this->point(1) - this->point(0)).cross(this->point(3) - this->point(0));
342 0 : n[2] = (this->point(2) - this->point(1)).cross(this->point(3) - this->point(1));
343 0 : n[3] = (this->point(0) - this->point(2)).cross(this->point(3) - this->point(2));
344 :
345 : Real dihedral_angles[6]; // 01, 02, 03, 12, 13, 23
346 :
347 : // Compute dihedral angles
348 0 : for (unsigned int k=0,i=0; i<4; ++i)
349 0 : for (unsigned int j=i+1; j<4; ++j,k+=1)
350 0 : dihedral_angles[k] = std::acos(n[i]*n[j] / n[i].norm() / n[j].norm()); // return value is between 0 and PI
351 :
352 : // Return max/min dihedral angles
353 0 : return std::make_pair(*std::min_element(dihedral_angles, dihedral_angles+6),
354 0 : *std::max_element(dihedral_angles, dihedral_angles+6));
355 :
356 : }
357 :
358 :
359 :
360 0 : dof_id_type Tet4::key () const
361 : {
362 0 : return this->compute_key(this->node_id(0),
363 : this->node_id(1),
364 : this->node_id(2),
365 0 : this->node_id(3));
366 : }
367 :
368 :
369 :
370 6604739 : bool Tet4::contains_point (const Point & p, Real tol) const
371 : {
372 : // See the response by Tony Noe on this thread.
373 : // http://bbs.dartmouth.edu/~fangq/MATH/download/source/Point_in_Tetrahedron.htm
374 : Point
375 394388 : col1 = point(1) - point(0),
376 197194 : col2 = point(2) - point(0),
377 197194 : col3 = point(3) - point(0);
378 :
379 197194 : Point r;
380 :
381 : libmesh_try
382 : {
383 6604739 : RealTensorValue(col1(0), col2(0), col3(0),
384 197194 : col1(1), col2(1), col3(1),
385 6801933 : col1(2), col2(2), col3(2)).solve(p - point(0), r);
386 : }
387 0 : libmesh_catch (ConvergenceFailure &)
388 : {
389 : // The Jacobian was singular, therefore the Tet had
390 : // zero volume. In this degenerate case, it is not
391 : // possible for the Tet to contain any points.
392 0 : return false;
393 0 : }
394 :
395 : // The point p is inside the tetrahedron if r1 + r2 + r3 < 1 and
396 : // 0 <= ri <= 1 for i = 1,2,3.
397 : return
398 10247006 : r(0) > -tol &&
399 3758703 : r(1) > -tol &&
400 8568559 : r(2) > -tol &&
401 1063637 : r(0) + r(1) + r(2) < 1.0 + tol;
402 : }
403 :
404 :
405 :
406 : #ifdef LIBMESH_ENABLE_AMR
407 3388382 : Real Tet4::embedding_matrix (const unsigned int i,
408 : const unsigned int j,
409 : const unsigned int k) const
410 : {
411 : // Choose an optimal diagonal, if one has not already been selected
412 3388382 : this->choose_diagonal();
413 :
414 : // Permuted j and k indices
415 : unsigned int
416 1126480 : jp=j,
417 1126480 : kp=k;
418 :
419 3388382 : if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
420 : {
421 : // Just the enum value cast to an unsigned int...
422 402740 : const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection);
423 :
424 : // Permute j, k:
425 : // ds==1 ds==2
426 : // 0 -> 1 0 -> 2
427 : // 1 -> 2 1 -> 0
428 : // 2 -> 0 2 -> 1
429 886362 : if (jp != 3)
430 664721 : jp = (jp+ds)%3;
431 :
432 886362 : if (kp != 3)
433 664923 : kp = (kp+ds)%3;
434 : }
435 :
436 : // Debugging
437 : // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
438 : // libMesh::err << "j=" << j << std::endl;
439 : // libMesh::err << "k=" << k << std::endl;
440 : // libMesh::err << "jp=" << jp << std::endl;
441 : // libMesh::err << "kp=" << kp << std::endl;
442 :
443 : // Call embedding matrix with permuted indices
444 3388382 : return this->_embedding_matrix[i][jp][kp];
445 : }
446 :
447 :
448 :
449 : // void Tet4::reselect_diagonal (const Diagonal diag)
450 : // {
451 : // /* Make sure that the element has just been refined. */
452 : // libmesh_assert(_children);
453 : // libmesh_assert_equal_to (n_children(), 8);
454 : // libmesh_assert_equal_to (_children[0]->refinement_flag(), JUST_REFINED);
455 : // libmesh_assert_equal_to (_children[1]->refinement_flag(), JUST_REFINED);
456 : // libmesh_assert_equal_to (_children[2]->refinement_flag(), JUST_REFINED);
457 : // libmesh_assert_equal_to (_children[3]->refinement_flag(), JUST_REFINED);
458 : // libmesh_assert_equal_to (_children[4]->refinement_flag(), JUST_REFINED);
459 : // libmesh_assert_equal_to (_children[5]->refinement_flag(), JUST_REFINED);
460 : // libmesh_assert_equal_to (_children[6]->refinement_flag(), JUST_REFINED);
461 : // libmesh_assert_equal_to (_children[7]->refinement_flag(), JUST_REFINED);
462 : //
463 : // /* Check whether anything has to be changed. */
464 : // if (_diagonal_selection!=diag)
465 : // {
466 : // /* Set new diagonal selection. */
467 : // _diagonal_selection = diag;
468 : //
469 : // /* The first four children do not have to be changed. For the
470 : // others, only the nodes have to be changed. Note that we have
471 : // to keep track of the nodes ourselves since there is no \p
472 : // MeshRefinement object with a valid \p _new_nodes_map
473 : // available. */
474 : // for (unsigned int c=4; c<this->n_children(); c++)
475 : // {
476 : // Elem * child = this->child_ptr(c);
477 : // for (unsigned int nc=0; nc<child->n_nodes(); nc++)
478 : // {
479 : // /* Unassign the current node. */
480 : // child->set_node(nc, nullptr);
481 : //
482 : // /* We have to find the correct new node now. We know
483 : // that it exists somewhere. We make use of the fact
484 : // that the embedding matrix for these children consists
485 : // of entries 0.0 and 0.5 only. Also, we make use of
486 : // the properties of the embedding matrix for the first
487 : // (unchanged) four children, which allow us to use a
488 : // simple mechanism to find the required node. */
489 : //
490 : //
491 : // unsigned int first_05_in_embedding_matrix = libMesh::invalid_uint;
492 : // for (unsigned int n=0; n<this->n_nodes(); n++)
493 : // {
494 : // if (this->embedding_matrix(c,nc,n) != 0.0)
495 : // {
496 : // /* It must be 0.5 then. Check whether it's the
497 : // first or second time that we get a 0.5
498 : // value. */
499 : // if (first_05_in_embedding_matrix==libMesh::invalid_uint)
500 : // {
501 : // /* First time, so just remember this position. */
502 : // first_05_in_embedding_matrix = n;
503 : // }
504 : // else
505 : // {
506 : // /* Second time, so we know now which node to
507 : // use. */
508 : // child->set_node(nc, this->child_ptr(n)->node_ptr(first_05_in_embedding_matrix));
509 : // }
510 : //
511 : // }
512 : // }
513 : //
514 : // /* Make sure that a node has been found. */
515 : // libmesh_assert(child->node_ptr(nc));
516 : // }
517 : // }
518 : // }
519 : // }
520 :
521 :
522 :
523 : // void Tet4::reselect_optimal_diagonal (const Diagonal exclude_this)
524 : // {
525 : // Real diag_01_23 = (this->point(0)+this->point(1)-this->point(2)-this->point(3)).norm_sq();
526 : // Real diag_02_13 = (this->point(0)-this->point(1)+this->point(2)-this->point(3)).norm_sq();
527 : // Real diag_03_12 = (this->point(0)-this->point(1)-this->point(2)+this->point(3)).norm_sq();
528 : //
529 : // Diagonal use_this = INVALID_DIAG;
530 : // switch (exclude_this)
531 : // {
532 : // case DIAG_01_23:
533 : // use_this = DIAG_02_13;
534 : // if (diag_03_12 < diag_02_13)
535 : // {
536 : // use_this = DIAG_03_12;
537 : // }
538 : // break;
539 : //
540 : // case DIAG_02_13:
541 : // use_this = DIAG_03_12;
542 : // if (diag_01_23 < diag_03_12)
543 : // {
544 : // use_this = DIAG_01_23;
545 : // }
546 : // break;
547 : //
548 : // case DIAG_03_12:
549 : // use_this = DIAG_02_13;
550 : // if (diag_01_23 < diag_02_13)
551 : // {
552 : // use_this = DIAG_01_23;
553 : // }
554 : // break;
555 : //
556 : // default:
557 : // use_this = DIAG_02_13;
558 : // if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
559 : // {
560 : // if (diag_01_23 < diag_03_12)
561 : // {
562 : // use_this = DIAG_01_23;
563 : // }
564 : // else
565 : // {
566 : // use_this = DIAG_03_12;
567 : // }
568 : // }
569 : // break;
570 : // }
571 : //
572 : // reselect_diagonal (use_this);
573 : // }
574 : #endif // #ifdef LIBMESH_ENABLE_AMR
575 :
576 155952 : void Tet4::permute(unsigned int perm_num)
577 : {
578 6192 : libmesh_assert_less (perm_num, 12);
579 :
580 155952 : const unsigned int side = perm_num % 4;
581 155952 : const unsigned int rotate = perm_num / 4;
582 :
583 322596 : for (unsigned int i = 0; i != rotate; ++i)
584 : {
585 166644 : swap3nodes(0,1,2);
586 160128 : swap3neighbors(1,2,3);
587 : }
588 :
589 155952 : switch (side) {
590 1872 : case 0:
591 1872 : break;
592 33642 : case 1:
593 33642 : swap3nodes(0,2,3);
594 32256 : swap3neighbors(0,2,1);
595 32256 : break;
596 28296 : case 2:
597 28296 : swap3nodes(2,0,3);
598 27072 : swap3neighbors(0,1,2);
599 27072 : break;
600 44334 : case 3:
601 44334 : swap3nodes(2,1,3);
602 42624 : swap3neighbors(0,1,3);
603 42624 : break;
604 0 : default:
605 0 : libmesh_error();
606 : }
607 155952 : }
608 :
609 :
610 146743 : void Tet4::flip(BoundaryInfo * boundary_info)
611 : {
612 6572 : libmesh_assert(boundary_info);
613 :
614 146743 : swap2nodes(0,2);
615 6572 : swap2neighbors(1,2);
616 146743 : swap2boundarysides(1,2,boundary_info);
617 146743 : swap2boundaryedges(0,1,boundary_info);
618 146743 : swap2boundaryedges(3,5,boundary_info);
619 146743 : }
620 :
621 :
622 157912 : ElemType Tet4::side_type (const unsigned int libmesh_dbg_var(s)) const
623 : {
624 13248 : libmesh_assert_less (s, 4);
625 157912 : return TRI3;
626 : }
627 :
628 :
629 : } // namespace libMesh
|
