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_tet.h"
21 : #include "libmesh/cell_tet4.h"
22 : #include "libmesh/face_tri3.h"
23 : #include "libmesh/enum_elem_quality.h"
24 :
25 : namespace libMesh
26 : {
27 :
28 :
29 :
30 : // ------------------------------------------------------------
31 : // Tet class static member initializations
32 : const int Tet::num_sides;
33 : const int Tet::num_edges;
34 : const int Tet::num_children;
35 :
36 : const Real Tet::_master_points[14][3] =
37 : {
38 : {0, 0, 0},
39 : {1, 0, 0},
40 : {0, 1, 0},
41 : {0, 0, 1},
42 : {0.5, 0, 0},
43 : {0.5, 0.5, 0},
44 : {0, 0.5, 0},
45 : {0, 0, 0.5},
46 : {0.5, 0, 0.5},
47 : {0, 0.5, 0.5},
48 : {1/Real(3), 1/Real(3), 0},
49 : {1/Real(3), 0, 1/Real(3)},
50 : {1/Real(3), 1/Real(3), 1/Real(3)},
51 : {0, 1/Real(3), 1/Real(3)}
52 : };
53 :
54 : const unsigned int Tet::edge_sides_map[6][2] =
55 : {
56 : {0, 1}, // Edge 0
57 : {0, 2}, // Edge 1
58 : {0, 3}, // Edge 2
59 : {1, 3}, // Edge 3
60 : {1, 2}, // Edge 4
61 : {2, 3} // Edge 5
62 : };
63 :
64 : const unsigned int Tet::adjacent_edges_map[/*num_vertices*/4][/*n_adjacent_edges*/3] =
65 : {
66 : {0, 2, 3}, // Edges adjacent to node 0
67 : {0, 1, 4}, // Edges adjacent to node 1
68 : {1, 2, 5}, // Edges adjacent to node 2
69 : {3, 4, 5}, // Edges adjacent to node 3
70 : };
71 :
72 : // ------------------------------------------------------------
73 : // Tet class member functions
74 288 : dof_id_type Tet::key (const unsigned int s) const
75 : {
76 24 : libmesh_assert_less (s, this->n_sides());
77 :
78 312 : return this->compute_key(this->node_id(Tet4::side_nodes_map[s][0]),
79 288 : this->node_id(Tet4::side_nodes_map[s][1]),
80 312 : this->node_id(Tet4::side_nodes_map[s][2]));
81 : }
82 :
83 :
84 :
85 44529196 : dof_id_type Tet::low_order_key (const unsigned int s) const
86 : {
87 1611664 : libmesh_assert_less (s, this->n_sides());
88 :
89 46097820 : return this->compute_key(this->node_id(Tet4::side_nodes_map[s][0]),
90 44529196 : this->node_id(Tet4::side_nodes_map[s][1]),
91 46140860 : this->node_id(Tet4::side_nodes_map[s][2]));
92 : }
93 :
94 :
95 :
96 4823 : unsigned int Tet::local_side_node(unsigned int side,
97 : unsigned int side_node) const
98 : {
99 398 : libmesh_assert_less (side, this->n_sides());
100 398 : libmesh_assert_less (side_node, Tet4::nodes_per_side);
101 :
102 4823 : return Tet4::side_nodes_map[side][side_node];
103 : }
104 :
105 :
106 :
107 9067857 : unsigned int Tet::local_edge_node(unsigned int edge,
108 : unsigned int edge_node) const
109 : {
110 410286 : libmesh_assert_less (edge, this->n_edges());
111 410286 : libmesh_assert_less (edge_node, Tet4::nodes_per_edge);
112 :
113 9067857 : return Tet4::edge_nodes_map[edge][edge_node];
114 : }
115 :
116 :
117 93197730 : std::unique_ptr<Elem> Tet::side_ptr (const unsigned int i)
118 : {
119 7758188 : libmesh_assert_less (i, this->n_sides());
120 :
121 93197730 : std::unique_ptr<Elem> face = std::make_unique<Tri3>();
122 :
123 372790920 : for (auto n : face->node_index_range())
124 302867718 : face->set_node(n, this->node_ptr(Tet4::side_nodes_map[i][n]));
125 :
126 93197730 : return face;
127 0 : }
128 :
129 :
130 :
131 40215766 : void Tet::side_ptr (std::unique_ptr<Elem> & side,
132 : const unsigned int i)
133 : {
134 40215766 : this->simple_side_ptr<Tet,Tet4>(side, i, TRI3);
135 40215766 : }
136 :
137 :
138 :
139 0 : void Tet::select_diagonal (const Diagonal diag) const
140 : {
141 0 : libmesh_assert_equal_to (_diagonal_selection, INVALID_DIAG);
142 0 : _diagonal_selection = diag;
143 0 : }
144 :
145 :
146 :
147 :
148 :
149 : #ifdef LIBMESH_ENABLE_AMR
150 :
151 :
152 5086380 : bool Tet::is_child_on_side_helper(const unsigned int c,
153 : const unsigned int s,
154 : const unsigned int checked_nodes[][3]) const
155 : {
156 1615424 : libmesh_assert_less (c, this->n_children());
157 1615424 : libmesh_assert_less (s, this->n_sides());
158 :
159 : // For the 4 vertices, child c touches vertex c, so we can return
160 : // true if that vertex is on side s
161 16408620 : for (unsigned int i = 0; i != 3; ++i)
162 13252551 : if (Tet4::side_nodes_map[s][i] == c)
163 583776 : return true;
164 :
165 : // If we are a "vertex child" and we didn't already return true,
166 : // we must not be on the side in question
167 3156069 : if (c < 4)
168 209264 : return false;
169 :
170 : // For the 4 non-vertex children, the child ordering depends on the
171 : // diagonal selection. We'll let the embedding matrix figure that
172 : // out: if this child has three nodes that don't depend on the
173 : // position of the node_facing_side[s], then we're on side s. Which
174 : // three nodes those are depends on the subclass, so their responsibility
175 : // is to call this function with the proper check_nodes array
176 2523304 : const unsigned int node_facing_side[4] = {3, 2, 0, 1};
177 2523304 : const unsigned int n = node_facing_side[s];
178 :
179 : // Add up the absolute values of the entries of the embedding matrix for the
180 : // nodes opposite node n. If it is equal to zero, then the child in question is
181 : // on side s, so return true.
182 822384 : Real embedding_sum = 0.;
183 10093216 : for (unsigned i=0; i<3; ++i)
184 7569912 : embedding_sum += std::abs(this->embedding_matrix(c, checked_nodes[n][i], n));
185 :
186 2523304 : return ( std::abs(embedding_sum) < 1.e-3 );
187 : }
188 :
189 : #else
190 :
191 : bool Tet::is_child_on_side_helper(const unsigned int /*c*/,
192 : const unsigned int /*s*/,
193 : const unsigned int /*checked_nodes*/[][3]) const
194 : {
195 : libmesh_not_implemented();
196 : return false;
197 : }
198 :
199 : #endif //LIBMESH_ENABLE_AMR
200 :
201 :
202 :
203 :
204 92862639 : void Tet::choose_diagonal() const
205 : {
206 : // If uninitialized diagonal selection, select the shortest octahedron diagonal
207 92862639 : if (this->_diagonal_selection==INVALID_DIAG)
208 : {
209 28012 : Real diag_01_23 = (this->point(0) + this->point(1) - this->point(2) - this->point(3)).norm_sq();
210 14006 : Real diag_02_13 = (this->point(0) - this->point(1) + this->point(2) - this->point(3)).norm_sq();
211 14006 : Real diag_03_12 = (this->point(0) - this->point(1) - this->point(2) + this->point(3)).norm_sq();
212 :
213 269245 : std::array<Real, 3> D = {diag_02_13, diag_03_12, diag_01_23};
214 :
215 269245 : this->_diagonal_selection =
216 269245 : Diagonal(std::distance(D.begin(), std::min_element(D.begin(), D.end())));
217 : }
218 92862639 : }
219 :
220 :
221 :
222 5221520 : bool Tet::is_edge_on_side(const unsigned int e,
223 : const unsigned int s) const
224 : {
225 5133212 : libmesh_assert_less (e, this->n_edges());
226 5133212 : libmesh_assert_less (s, this->n_sides());
227 :
228 5221520 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
229 : }
230 :
231 :
232 :
233 166740 : std::vector<unsigned int> Tet::sides_on_edge(const unsigned int e) const
234 : {
235 13848 : libmesh_assert_less (e, this->n_edges());
236 166740 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
237 : }
238 :
239 :
240 :
241 : bool
242 377146 : Tet::is_flipped() const
243 : {
244 360502 : return (triple_product(this->point(1)-this->point(0),
245 360502 : this->point(2)-this->point(0),
246 410434 : this->point(3)-this->point(0)) < 0);
247 : }
248 :
249 :
250 : std::vector<unsigned int>
251 324264 : Tet::edges_adjacent_to_node(const unsigned int n) const
252 : {
253 26928 : libmesh_assert_less(n, this->n_nodes());
254 :
255 : // For vertices, we use the Tet::adjacent_sides_map, otherwise each
256 : // of the mid-edge nodes is adjacent only to the edge it is on, and the
257 : // mid-face nodes are not adjacent to any edges.
258 324264 : if (this->is_vertex(n))
259 151512 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n])};
260 184320 : else if (this->is_edge(n))
261 138240 : return {n - this->n_vertices()};
262 :
263 : // Current Tets have only vertex, edge, and face nodes.
264 3840 : libmesh_assert(this->is_face(n));
265 42240 : return {};
266 : }
267 :
268 :
269 48952 : Real Tet::quality(const ElemQuality q) const
270 : {
271 48952 : return Elem::quality(q); // Not implemented
272 : }
273 :
274 :
275 :
276 :
277 0 : std::pair<Real, Real> Tet::qual_bounds (const ElemQuality q) const
278 : {
279 0 : std::pair<Real, Real> bounds;
280 :
281 0 : switch (q)
282 : {
283 :
284 0 : case ASPECT_RATIO_BETA:
285 : case ASPECT_RATIO_GAMMA:
286 0 : bounds.first = 1.;
287 0 : bounds.second = 3.;
288 0 : break;
289 :
290 0 : case SIZE:
291 : case SHAPE:
292 0 : bounds.first = 0.2;
293 0 : bounds.second = 1.;
294 0 : break;
295 :
296 0 : case CONDITION:
297 0 : bounds.first = 1.;
298 0 : bounds.second = 3.;
299 0 : break;
300 :
301 0 : case DISTORTION:
302 0 : bounds.first = 0.6;
303 0 : bounds.second = 1.;
304 0 : break;
305 :
306 0 : case JACOBIAN:
307 : case SCALED_JACOBIAN:
308 0 : bounds.first = 0.5;
309 0 : bounds.second = 1.414;
310 0 : break;
311 :
312 0 : default:
313 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
314 0 : bounds.first = -1;
315 0 : bounds.second = -1;
316 : }
317 :
318 0 : return bounds;
319 : }
320 :
321 :
322 :
323 19036304 : bool Tet::on_reference_element(const Point & p,
324 : const Real eps) const
325 : {
326 521265 : const Real & xi = p(0);
327 521265 : const Real & eta = p(1);
328 521265 : const Real & zeta = p(2);
329 : // The reference tetrahedral is isosceles
330 : // and is bound by xi=0, eta=0, zeta=0,
331 : // and xi+eta+zeta=1.
332 30099408 : return ((xi >= 0.-eps) &&
333 11063104 : (eta >= 0.-eps) &&
334 24655543 : (zeta >= 0.-eps) &&
335 3172270 : ((xi + eta + zeta) <= 1.+eps));
336 : }
337 :
338 :
339 : } // namespace libMesh
|
