1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
33
34 // We need to require C++11...
35 const Real Tet::_master_points[10][3] =
36 {
37 {0, 0, 0},
38 {1, 0, 0},
39 {0, 1, 0},
40 {0, 0, 1},
41 {0.5, 0, 0},
42 {0.5, 0.5, 0},
43 {0, 0.5, 0},
44 {0, 0, 0.5},
45 {0.5, 0, 0.5},
46 {0, 0.5, 0.5}
47 };
48
49
50
51
52 // ------------------------------------------------------------
53 // Tet class member functions
key(const unsigned int s)54 dof_id_type Tet::key (const unsigned int s) const
55 {
56 libmesh_assert_less (s, this->n_sides());
57
58 return this->compute_key(this->node_id(Tet4::side_nodes_map[s][0]),
59 this->node_id(Tet4::side_nodes_map[s][1]),
60 this->node_id(Tet4::side_nodes_map[s][2]));
61 }
62
63
64
local_side_node(unsigned int side,unsigned int side_node)65 unsigned int Tet::local_side_node(unsigned int side,
66 unsigned int side_node) const
67 {
68 libmesh_assert_less (side, this->n_sides());
69 libmesh_assert_less (side_node, Tet4::nodes_per_side);
70
71 return Tet4::side_nodes_map[side][side_node];
72 }
73
74
75
local_edge_node(unsigned int edge,unsigned int edge_node)76 unsigned int Tet::local_edge_node(unsigned int edge,
77 unsigned int edge_node) const
78 {
79 libmesh_assert_less (edge, this->n_edges());
80 libmesh_assert_less (edge_node, Tet4::nodes_per_edge);
81
82 return Tet4::edge_nodes_map[edge][edge_node];
83 }
84
85
side_ptr(const unsigned int i)86 std::unique_ptr<Elem> Tet::side_ptr (const unsigned int i)
87 {
88 libmesh_assert_less (i, this->n_sides());
89
90 std::unique_ptr<Elem> face = libmesh_make_unique<Tri3>();
91
92 for (auto n : face->node_index_range())
93 face->set_node(n) = this->node_ptr(Tet4::side_nodes_map[i][n]);
94
95 return face;
96 }
97
98
99
side_ptr(std::unique_ptr<Elem> & side,const unsigned int i)100 void Tet::side_ptr (std::unique_ptr<Elem> & side,
101 const unsigned int i)
102 {
103 this->simple_side_ptr<Tet,Tet4>(side, i, TRI3);
104 }
105
106
107
select_diagonal(const Diagonal diag)108 void Tet::select_diagonal (const Diagonal diag) const
109 {
110 libmesh_assert_equal_to (_diagonal_selection, INVALID_DIAG);
111 _diagonal_selection = diag;
112 }
113
114
115
116
117
118 #ifdef LIBMESH_ENABLE_AMR
119
120
is_child_on_side_helper(const unsigned int c,const unsigned int s,const unsigned int checked_nodes[][3])121 bool Tet::is_child_on_side_helper(const unsigned int c,
122 const unsigned int s,
123 const unsigned int checked_nodes[][3]) const
124 {
125 libmesh_assert_less (c, this->n_children());
126 libmesh_assert_less (s, this->n_sides());
127
128 // For the 4 vertices, child c touches vertex c, so we can return
129 // true if that vertex is on side s
130 for (unsigned int i = 0; i != 3; ++i)
131 if (Tet4::side_nodes_map[s][i] == c)
132 return true;
133
134 // If we are a "vertex child" and we didn't already return true,
135 // we must not be on the side in question
136 if (c < 4)
137 return false;
138
139 // For the 4 non-vertex children, the child ordering depends on the
140 // diagonal selection. We'll let the embedding matrix figure that
141 // out: if this child has three nodes that don't depend on the
142 // position of the node_facing_side[s], then we're on side s. Which
143 // three nodes those are depends on the subclass, so their responsibility
144 // is to call this function with the proper check_nodes array
145 const unsigned int node_facing_side[4] = {3, 2, 0, 1};
146 const unsigned int n = node_facing_side[s];
147
148 // Add up the absolute values of the entries of the embedding matrix for the
149 // nodes opposite node n. If it is equal to zero, then the child in question is
150 // on side s, so return true.
151 Real embedding_sum = 0.;
152 for (unsigned i=0; i<3; ++i)
153 embedding_sum += std::abs(this->embedding_matrix(c, checked_nodes[n][i], n));
154
155 return ( std::abs(embedding_sum) < 1.e-3 );
156 }
157
158 #else
159
is_child_on_side_helper(const unsigned int,const unsigned int,const unsigned int[][3])160 bool Tet::is_child_on_side_helper(const unsigned int /*c*/,
161 const unsigned int /*s*/,
162 const unsigned int /*checked_nodes*/[][3]) const
163 {
164 libmesh_not_implemented();
165 return false;
166 }
167
168 #endif //LIBMESH_ENABLE_AMR
169
170
171
172
choose_diagonal()173 void Tet::choose_diagonal() const
174 {
175 // Check for uninitialized diagonal selection
176 if (this->_diagonal_selection==INVALID_DIAG)
177 {
178 Real diag_01_23 = (this->point(0) + this->point(1) - this->point(2) - this->point(3)).norm_sq();
179 Real diag_02_13 = (this->point(0) - this->point(1) + this->point(2) - this->point(3)).norm_sq();
180 Real diag_03_12 = (this->point(0) - this->point(1) - this->point(2) + this->point(3)).norm_sq();
181
182 this->_diagonal_selection=DIAG_02_13;
183
184 if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
185 {
186 if (diag_01_23 < diag_03_12)
187 this->_diagonal_selection=DIAG_01_23;
188
189 else
190 this->_diagonal_selection=DIAG_03_12;
191 }
192 }
193 }
194
195
196
is_edge_on_side(const unsigned int e,const unsigned int s)197 bool Tet::is_edge_on_side(const unsigned int e,
198 const unsigned int s) const
199 {
200 libmesh_assert_less (e, this->n_edges());
201 libmesh_assert_less (s, this->n_sides());
202
203 return (Tet4::edge_sides_map[e][0] == s ||
204 Tet4::edge_sides_map[e][1] == s);
205 }
206
207
208
sides_on_edge(const unsigned int e)209 std::vector<unsigned int> Tet::sides_on_edge(const unsigned int e) const
210 {
211 return {Tet4::edge_sides_map[e][0], Tet4::edge_sides_map[e][1]};
212 }
213
214
215
quality(const ElemQuality q)216 Real Tet::quality(const ElemQuality q) const
217 {
218 return Elem::quality(q); // Not implemented
219 }
220
221
222
223
qual_bounds(const ElemQuality q)224 std::pair<Real, Real> Tet::qual_bounds (const ElemQuality q) const
225 {
226 std::pair<Real, Real> bounds;
227
228 switch (q)
229 {
230
231 case ASPECT_RATIO_BETA:
232 case ASPECT_RATIO_GAMMA:
233 bounds.first = 1.;
234 bounds.second = 3.;
235 break;
236
237 case SIZE:
238 case SHAPE:
239 bounds.first = 0.2;
240 bounds.second = 1.;
241 break;
242
243 case CONDITION:
244 bounds.first = 1.;
245 bounds.second = 3.;
246 break;
247
248 case DISTORTION:
249 bounds.first = 0.6;
250 bounds.second = 1.;
251 break;
252
253 case JACOBIAN:
254 bounds.first = 0.5;
255 bounds.second = 1.414;
256 break;
257
258 default:
259 libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
260 bounds.first = -1;
261 bounds.second = -1;
262 }
263
264 return bounds;
265 }
266
267 } // namespace libMesh
268