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_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 26480889 : bool Tet4::is_vertex(const unsigned int libmesh_dbg_var(n)) const
60 : {
61 1075059 : libmesh_assert_not_equal_to (n, invalid_uint);
62 26480889 : 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 8347944 : bool Tet4::is_node_on_edge(const unsigned int n,
76 : const unsigned int e) const
77 : {
78 262416 : libmesh_assert_less (e, n_edges());
79 262416 : return std::find(std::begin(edge_nodes_map[e]),
80 8347944 : std::end(edge_nodes_map[e]),
81 8347944 : 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 11648810 : 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 11648810 : 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 12359146 : 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 233696 : bool Tet4::is_node_on_side(const unsigned int n,
129 : const unsigned int s) const
130 : {
131 17432 : libmesh_assert_less (s, n_sides());
132 17432 : return std::find(std::begin(side_nodes_map[s]),
133 233696 : std::end(side_nodes_map[s]),
134 233696 : n) != std::end(side_nodes_map[s]);
135 : }
136 :
137 : std::vector<unsigned>
138 212987 : Tet4::nodes_on_side(const unsigned int s) const
139 : {
140 16688 : libmesh_assert_less(s, n_sides());
141 229675 : 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 15395624 : Order Tet4::default_order() const
152 : {
153 15395624 : return FIRST;
154 : }
155 :
156 4772268 : std::unique_ptr<Elem> Tet4::build_side_ptr (const unsigned int i)
157 : {
158 4772268 : return this->simple_build_side_ptr<Tri3, Tet4>(i);
159 : }
160 :
161 :
162 :
163 11044464 : void Tet4::build_side_ptr (std::unique_ptr<Elem> & side,
164 : const unsigned int i)
165 : {
166 11044464 : this->simple_build_side_ptr<Tet4>(side, i, TRI3);
167 11044464 : }
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 145894 : Real Tet4::volume () const
318 : {
319 : // This specialization is good for non-Lagrange mappings too! Even
320 : // with differing vertex weights we've still got p=1 planar faces.
321 :
322 : // The volume of a tetrahedron is 1/6 the box product formed
323 : // by its base and apex vectors
324 28344 : Point a = point(3) - point(0);
325 :
326 : // b is the vector pointing from 0 to 1
327 14173 : Point b = point(1) - point(0);
328 :
329 : // c is the vector pointing from 0 to 2
330 14173 : Point c = point(2) - point(0);
331 :
332 145894 : return triple_product(a, b, c) / 6.;
333 : }
334 :
335 :
336 :
337 :
338 0 : std::pair<Real, Real> Tet4::min_and_max_angle() const
339 : {
340 0 : Point n[4];
341 :
342 : // Compute the outward normal vectors on each face
343 0 : n[0] = (this->point(2) - this->point(0)).cross(this->point(1) - this->point(0));
344 0 : n[1] = (this->point(1) - this->point(0)).cross(this->point(3) - this->point(0));
345 0 : n[2] = (this->point(2) - this->point(1)).cross(this->point(3) - this->point(1));
346 0 : n[3] = (this->point(0) - this->point(2)).cross(this->point(3) - this->point(2));
347 :
348 : Real dihedral_angles[6]; // 01, 02, 03, 12, 13, 23
349 :
350 : // Compute dihedral angles
351 0 : for (unsigned int k=0,i=0; i<4; ++i)
352 0 : for (unsigned int j=i+1; j<4; ++j,k+=1)
353 0 : dihedral_angles[k] = std::acos(n[i]*n[j] / n[i].norm() / n[j].norm()); // return value is between 0 and PI
354 :
355 : // Return max/min dihedral angles
356 0 : return std::make_pair(*std::min_element(dihedral_angles, dihedral_angles+6),
357 0 : *std::max_element(dihedral_angles, dihedral_angles+6));
358 :
359 : }
360 :
361 :
362 :
363 0 : dof_id_type Tet4::key () const
364 : {
365 0 : return this->compute_key(this->node_id(0),
366 : this->node_id(1),
367 : this->node_id(2),
368 0 : this->node_id(3));
369 : }
370 :
371 :
372 :
373 13165743 : bool Tet4::contains_point (const Point & p, Real tol) const
374 : {
375 : // See the response by Tony Noe on this thread.
376 : // http://bbs.dartmouth.edu/~fangq/MATH/download/source/Point_in_Tetrahedron.htm
377 : Point
378 777108 : col1 = point(1) - point(0),
379 388554 : col2 = point(2) - point(0),
380 388554 : col3 = point(3) - point(0);
381 :
382 388554 : Point r;
383 :
384 : libmesh_try
385 : {
386 13165743 : RealTensorValue(col1(0), col2(0), col3(0),
387 388554 : col1(1), col2(1), col3(1),
388 13554297 : col1(2), col2(2), col3(2)).solve(p - point(0), r);
389 : }
390 0 : libmesh_catch (ConvergenceFailure &)
391 : {
392 : // The Jacobian was singular, therefore the Tet had
393 : // zero volume. In this degenerate case, it is not
394 : // possible for the Tet to contain any points.
395 0 : return false;
396 0 : }
397 :
398 : // The point p is inside the tetrahedron if r1 + r2 + r3 < 1 and
399 : // 0 <= ri <= 1 for i = 1,2,3.
400 : return
401 20135166 : r(0) > -tol &&
402 7185859 : r(1) > -tol &&
403 16465777 : r(2) > -tol &&
404 1846367 : r(0) + r(1) + r(2) < 1.0 + tol;
405 : }
406 :
407 :
408 :
409 : #ifdef LIBMESH_ENABLE_AMR
410 18001292 : Real Tet4::embedding_matrix (const unsigned int i,
411 : const unsigned int j,
412 : const unsigned int k) const
413 : {
414 : // Choose an optimal diagonal, if one has not already been selected
415 18001292 : this->choose_diagonal();
416 :
417 : // Permuted j and k indices
418 : unsigned int
419 1126480 : jp=j,
420 1126480 : kp=k;
421 :
422 18001292 : if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
423 : {
424 : // Just the enum value cast to an unsigned int...
425 402740 : const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection);
426 :
427 : // Permute j, k:
428 : // ds==1 ds==2
429 : // 0 -> 1 0 -> 2
430 : // 1 -> 2 1 -> 0
431 : // 2 -> 0 2 -> 1
432 4382382 : if (jp != 3)
433 3308759 : jp = (jp+ds)%3;
434 :
435 4382382 : if (kp != 3)
436 3220869 : kp = (kp+ds)%3;
437 : }
438 :
439 : // Debugging
440 : // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
441 : // libMesh::err << "j=" << j << std::endl;
442 : // libMesh::err << "k=" << k << std::endl;
443 : // libMesh::err << "jp=" << jp << std::endl;
444 : // libMesh::err << "kp=" << kp << std::endl;
445 :
446 : // Call embedding matrix with permuted indices
447 18001292 : return this->_embedding_matrix[i][jp][kp];
448 : }
449 :
450 :
451 :
452 : // void Tet4::reselect_diagonal (const Diagonal diag)
453 : // {
454 : // /* Make sure that the element has just been refined. */
455 : // libmesh_assert(_children);
456 : // libmesh_assert_equal_to (n_children(), 8);
457 : // libmesh_assert_equal_to (_children[0]->refinement_flag(), JUST_REFINED);
458 : // libmesh_assert_equal_to (_children[1]->refinement_flag(), JUST_REFINED);
459 : // libmesh_assert_equal_to (_children[2]->refinement_flag(), JUST_REFINED);
460 : // libmesh_assert_equal_to (_children[3]->refinement_flag(), JUST_REFINED);
461 : // libmesh_assert_equal_to (_children[4]->refinement_flag(), JUST_REFINED);
462 : // libmesh_assert_equal_to (_children[5]->refinement_flag(), JUST_REFINED);
463 : // libmesh_assert_equal_to (_children[6]->refinement_flag(), JUST_REFINED);
464 : // libmesh_assert_equal_to (_children[7]->refinement_flag(), JUST_REFINED);
465 : //
466 : // /* Check whether anything has to be changed. */
467 : // if (_diagonal_selection!=diag)
468 : // {
469 : // /* Set new diagonal selection. */
470 : // _diagonal_selection = diag;
471 : //
472 : // /* The first four children do not have to be changed. For the
473 : // others, only the nodes have to be changed. Note that we have
474 : // to keep track of the nodes ourselves since there is no \p
475 : // MeshRefinement object with a valid \p _new_nodes_map
476 : // available. */
477 : // for (unsigned int c=4; c<this->n_children(); c++)
478 : // {
479 : // Elem * child = this->child_ptr(c);
480 : // for (unsigned int nc=0; nc<child->n_nodes(); nc++)
481 : // {
482 : // /* Unassign the current node. */
483 : // child->set_node(nc, nullptr);
484 : //
485 : // /* We have to find the correct new node now. We know
486 : // that it exists somewhere. We make use of the fact
487 : // that the embedding matrix for these children consists
488 : // of entries 0.0 and 0.5 only. Also, we make use of
489 : // the properties of the embedding matrix for the first
490 : // (unchanged) four children, which allow us to use a
491 : // simple mechanism to find the required node. */
492 : //
493 : //
494 : // unsigned int first_05_in_embedding_matrix = libMesh::invalid_uint;
495 : // for (unsigned int n=0; n<this->n_nodes(); n++)
496 : // {
497 : // if (this->embedding_matrix(c,nc,n) != 0.0)
498 : // {
499 : // /* It must be 0.5 then. Check whether it's the
500 : // first or second time that we get a 0.5
501 : // value. */
502 : // if (first_05_in_embedding_matrix==libMesh::invalid_uint)
503 : // {
504 : // /* First time, so just remember this position. */
505 : // first_05_in_embedding_matrix = n;
506 : // }
507 : // else
508 : // {
509 : // /* Second time, so we know now which node to
510 : // use. */
511 : // child->set_node(nc, this->child_ptr(n)->node_ptr(first_05_in_embedding_matrix));
512 : // }
513 : //
514 : // }
515 : // }
516 : //
517 : // /* Make sure that a node has been found. */
518 : // libmesh_assert(child->node_ptr(nc));
519 : // }
520 : // }
521 : // }
522 : // }
523 :
524 :
525 :
526 : // void Tet4::reselect_optimal_diagonal (const Diagonal exclude_this)
527 : // {
528 : // Real diag_01_23 = (this->point(0)+this->point(1)-this->point(2)-this->point(3)).norm_sq();
529 : // Real diag_02_13 = (this->point(0)-this->point(1)+this->point(2)-this->point(3)).norm_sq();
530 : // Real diag_03_12 = (this->point(0)-this->point(1)-this->point(2)+this->point(3)).norm_sq();
531 : //
532 : // Diagonal use_this = INVALID_DIAG;
533 : // switch (exclude_this)
534 : // {
535 : // case DIAG_01_23:
536 : // use_this = DIAG_02_13;
537 : // if (diag_03_12 < diag_02_13)
538 : // {
539 : // use_this = DIAG_03_12;
540 : // }
541 : // break;
542 : //
543 : // case DIAG_02_13:
544 : // use_this = DIAG_03_12;
545 : // if (diag_01_23 < diag_03_12)
546 : // {
547 : // use_this = DIAG_01_23;
548 : // }
549 : // break;
550 : //
551 : // case DIAG_03_12:
552 : // use_this = DIAG_02_13;
553 : // if (diag_01_23 < diag_02_13)
554 : // {
555 : // use_this = DIAG_01_23;
556 : // }
557 : // break;
558 : //
559 : // default:
560 : // use_this = DIAG_02_13;
561 : // if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
562 : // {
563 : // if (diag_01_23 < diag_03_12)
564 : // {
565 : // use_this = DIAG_01_23;
566 : // }
567 : // else
568 : // {
569 : // use_this = DIAG_03_12;
570 : // }
571 : // }
572 : // break;
573 : // }
574 : //
575 : // reselect_diagonal (use_this);
576 : // }
577 : #endif // #ifdef LIBMESH_ENABLE_AMR
578 :
579 155952 : void Tet4::permute(unsigned int perm_num)
580 : {
581 6192 : libmesh_assert_less (perm_num, 12);
582 :
583 155952 : const unsigned int side = perm_num % 4;
584 155952 : const unsigned int rotate = perm_num / 4;
585 :
586 322596 : for (unsigned int i = 0; i != rotate; ++i)
587 : {
588 166644 : swap3nodes(0,1,2);
589 160128 : swap3neighbors(1,2,3);
590 : }
591 :
592 155952 : switch (side) {
593 1872 : case 0:
594 1872 : break;
595 33642 : case 1:
596 33642 : swap3nodes(0,2,3);
597 32256 : swap3neighbors(0,2,1);
598 32256 : break;
599 28296 : case 2:
600 28296 : swap3nodes(2,0,3);
601 27072 : swap3neighbors(0,1,2);
602 27072 : break;
603 44334 : case 3:
604 44334 : swap3nodes(2,1,3);
605 42624 : swap3neighbors(0,1,3);
606 42624 : break;
607 0 : default:
608 0 : libmesh_error();
609 : }
610 155952 : }
611 :
612 :
613 146710 : void Tet4::flip(BoundaryInfo * boundary_info)
614 : {
615 6576 : libmesh_assert(boundary_info);
616 :
617 146710 : swap2nodes(0,2);
618 6576 : swap2neighbors(1,2);
619 146710 : swap2boundarysides(1,2,boundary_info);
620 146710 : swap2boundaryedges(0,1,boundary_info);
621 146710 : swap2boundaryedges(3,5,boundary_info);
622 146710 : }
623 :
624 :
625 182620 : ElemType Tet4::side_type (const unsigned int libmesh_dbg_var(s)) const
626 : {
627 15096 : libmesh_assert_less (s, 4);
628 182620 : return TRI3;
629 : }
630 :
631 :
632 : Point
633 284 : Tet4::side_vertex_average_normal(const unsigned int s) const
634 : {
635 8 : libmesh_assert_less (s, 4);
636 300 : const Point n1 = this->point(side_nodes_map[s][1]) -
637 284 : this->point(side_nodes_map[s][0]);
638 300 : const Point n2 = this->point(side_nodes_map[s][2]) -
639 8 : this->point(side_nodes_map[s][1]);
640 276 : const Point pointing_out = n1.cross(n2);
641 284 : return pointing_out.unit();
642 : }
643 :
644 :
645 : } // namespace libMesh
|
