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