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 : // C++ includes
20 :
21 : // Local includes
22 : #include "libmesh/cell_prism.h"
23 : #include "libmesh/cell_prism6.h"
24 : #include "libmesh/enum_elem_quality.h"
25 : #include "libmesh/face_quad4.h"
26 : #include "libmesh/face_tri3.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 : // ------------------------------------------------------------
33 : // Prism class static member initializations
34 : const int Prism::num_sides;
35 : const int Prism::num_edges;
36 : const int Prism::num_children;
37 :
38 : const Real Prism::_master_points[21][3] =
39 : {
40 : {0, 0, -1},
41 : {1, 0, -1},
42 : {0, 1, -1},
43 : {0, 0, 1},
44 : {1, 0, 1},
45 : {0, 1, 1},
46 : {0.5, 0, -1},
47 : {0.5, 0.5, -1},
48 : {0, 0.5, -1},
49 : {0, 0, 0},
50 : {1, 0, 0},
51 : {0, 1, 0},
52 : {0.5, 0, 1},
53 : {0.5, 0.5, 1},
54 : {0, 0.5, 1},
55 : {0.5, 0, 0},
56 : {0.5, 0.5, 0},
57 : {0, 0.5, 0},
58 : {Real(1)/3, Real(1)/3, -1},
59 : {Real(1)/3, Real(1)/3, 1},
60 : {Real(1)/3, Real(1)/3, 0}
61 : };
62 :
63 : const unsigned int Prism::edge_sides_map[9][2] =
64 : {
65 : {0, 1}, // Edge 0
66 : {0, 2}, // Edge 1
67 : {0, 3}, // Edge 2
68 : {1, 3}, // Edge 3
69 : {1, 2}, // Edge 4
70 : {2, 3}, // Edge 5
71 : {1, 4}, // Edge 6
72 : {2, 4}, // Edge 7
73 : {3, 4} // Edge 8
74 : };
75 :
76 : const unsigned int Prism::adjacent_edges_map[/*num_vertices*/6][/*n_adjacent_edges*/3] =
77 : {
78 : {0, 2, 3}, // Edges adjacent to node 0
79 : {0, 1, 4}, // Edges adjacent to node 1
80 : {1, 2, 5}, // Edges adjacent to node 2
81 : {3, 6, 8}, // Edges adjacent to node 3
82 : {4, 6, 7}, // Edges adjacent to node 4
83 : {5, 7, 8}, // Edges adjacent to node 5
84 : };
85 :
86 :
87 : // ------------------------------------------------------------
88 : // Prism class member functions
89 0 : dof_id_type Prism::key (const unsigned int s) const
90 : {
91 0 : libmesh_assert_less (s, this->n_sides());
92 :
93 0 : switch (s)
94 : {
95 0 : case 0: // the triangular face at z=0
96 : case 4: // the triangular face at z=1
97 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
98 0 : this->node_id(Prism6::side_nodes_map[s][1]),
99 0 : this->node_id(Prism6::side_nodes_map[s][2]));
100 :
101 0 : case 1: // the quad face at y=0
102 : case 2: // the other quad face
103 : case 3: // the quad face at x=0
104 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
105 0 : this->node_id(Prism6::side_nodes_map[s][1]),
106 0 : this->node_id(Prism6::side_nodes_map[s][2]),
107 0 : this->node_id(Prism6::side_nodes_map[s][3]));
108 :
109 0 : default:
110 0 : libmesh_error_msg("Invalid side " << s);
111 : }
112 : }
113 :
114 :
115 :
116 5327330 : dof_id_type Prism::low_order_key (const unsigned int s) const
117 : {
118 237760 : libmesh_assert_less (s, this->n_sides());
119 :
120 5327330 : switch (s)
121 : {
122 2119343 : case 0: // the triangular face at z=0
123 : case 4: // the triangular face at z=1
124 2214447 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
125 2119343 : this->node_id(Prism6::side_nodes_map[s][1]),
126 2214447 : this->node_id(Prism6::side_nodes_map[s][2]));
127 :
128 3207987 : case 1: // the quad face at y=0
129 : case 2: // the other quad face
130 : case 3: // the quad face at x=0
131 3350643 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
132 3207987 : this->node_id(Prism6::side_nodes_map[s][1]),
133 3207987 : this->node_id(Prism6::side_nodes_map[s][2]),
134 3350643 : this->node_id(Prism6::side_nodes_map[s][3]));
135 :
136 0 : default:
137 0 : libmesh_error_msg("Invalid side " << s);
138 : }
139 : }
140 :
141 :
142 :
143 295832854 : unsigned int Prism::local_side_node(unsigned int side,
144 : unsigned int side_node) const
145 : {
146 26436500 : libmesh_assert_less (side, this->n_sides());
147 :
148 : // Never more than 4 nodes per side.
149 26436500 : libmesh_assert_less(side_node, Prism6::nodes_per_side);
150 :
151 : // Some sides have 3 nodes.
152 26436500 : libmesh_assert(!(side==0 || side==4) || side_node < 3);
153 :
154 295832852 : return Prism6::side_nodes_map[side][side_node];
155 : }
156 :
157 :
158 :
159 35517078 : unsigned int Prism::local_edge_node(unsigned int edge,
160 : unsigned int edge_node) const
161 : {
162 3171564 : libmesh_assert_less(edge, this->n_edges());
163 3171564 : libmesh_assert_less(edge_node, Prism6::nodes_per_edge);
164 :
165 35517078 : return Prism6::edge_nodes_map[edge][edge_node];
166 : }
167 :
168 :
169 :
170 2035233 : std::unique_ptr<Elem> Prism::side_ptr (const unsigned int i)
171 : {
172 107354 : libmesh_assert_less (i, this->n_sides());
173 :
174 2035233 : std::unique_ptr<Elem> face;
175 :
176 : // Set up the type of element
177 2035233 : switch (i)
178 : {
179 968773 : case 0: // the triangular face at z=0
180 : case 4: // the triangular face at z=1
181 : {
182 968773 : face = std::make_unique<Tri3>();
183 968773 : break;
184 : }
185 1066460 : case 1: // the quad face at y=0
186 : case 2: // the other quad face
187 : case 3: // the quad face at x=0
188 : {
189 1066460 : face = std::make_unique<Quad4>();
190 1066460 : break;
191 : }
192 0 : default:
193 0 : libmesh_error_msg("Invalid side i = " << i);
194 : }
195 :
196 : // Set the nodes
197 9207392 : for (auto n : face->node_index_range())
198 7549217 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
199 :
200 2035233 : return face;
201 0 : }
202 :
203 :
204 :
205 3666521 : void Prism::side_ptr (std::unique_ptr<Elem> & side,
206 : const unsigned int i)
207 : {
208 192334 : libmesh_assert_less (i, this->n_sides());
209 :
210 3666521 : switch (i)
211 : {
212 : // the base face
213 73008 : case 0: // the triangular face at z=0
214 : case 4: // the triangular face at z=1
215 : {
216 1319992 : if (!side.get() || side->type() != TRI3)
217 : {
218 1831444 : side = this->side_ptr(i);
219 968058 : return;
220 : }
221 20672 : break;
222 : }
223 :
224 119326 : case 1: // the quad face at y=0
225 : case 2: // the other quad face
226 : case 3: // the quad face at x=0
227 : {
228 2346529 : if (!side.get() || side->type() != QUAD4)
229 : {
230 2020858 : side = this->side_ptr(i);
231 1065395 : return;
232 : }
233 64360 : break;
234 : }
235 :
236 0 : default:
237 0 : libmesh_error_msg("Invalid side i = " << i);
238 : }
239 :
240 1718100 : side->subdomain_id() = this->subdomain_id();
241 :
242 : // Set the nodes
243 7813406 : for (auto n : side->node_index_range())
244 6499794 : side->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
245 : }
246 :
247 :
248 :
249 9077437 : bool Prism::is_child_on_side(const unsigned int c,
250 : const unsigned int s) const
251 : {
252 547932 : libmesh_assert_less (c, this->n_children());
253 547932 : libmesh_assert_less (s, this->n_sides());
254 :
255 34008071 : for (unsigned int i = 0; i != 4; ++i)
256 29480391 : if (Prism6::side_elems_map[s][i] == c)
257 242584 : return true;
258 305348 : return false;
259 : }
260 :
261 :
262 :
263 3364720 : bool Prism::is_edge_on_side(const unsigned int e,
264 : const unsigned int s) const
265 : {
266 2886286 : libmesh_assert_less (e, this->n_edges());
267 2886286 : libmesh_assert_less (s, this->n_sides());
268 :
269 3364720 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
270 : }
271 :
272 :
273 :
274 155520 : std::vector<unsigned int> Prism::sides_on_edge(const unsigned int e) const
275 : {
276 12960 : libmesh_assert_less(e, this->n_edges());
277 :
278 155520 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
279 : }
280 :
281 :
282 :
283 0 : unsigned int Prism::opposite_side(const unsigned int side_in) const
284 : {
285 0 : libmesh_assert_less (side_in, 5);
286 : static const unsigned int prism_opposites[5] =
287 : {4, invalid_uint, invalid_uint, invalid_uint, 0};
288 0 : return prism_opposites[side_in];
289 : }
290 :
291 :
292 :
293 : bool
294 9428 : Prism::is_flipped() const
295 : {
296 8868 : return (triple_product(this->point(1)-this->point(0),
297 8868 : this->point(2)-this->point(0),
298 10548 : this->point(3)-this->point(0)) < 0);
299 : }
300 :
301 :
302 : std::vector<unsigned int>
303 291840 : Prism::edges_adjacent_to_node(const unsigned int n) const
304 : {
305 24320 : libmesh_assert_less(n, this->n_nodes());
306 :
307 : // For vertices, we use the Prism::adjacent_sides_map, otherwise each
308 : // of the mid-edge nodes is adjacent only to the edge it is on, and
309 : // face/internal nodes are not adjacent to any edge.
310 291840 : if (this->is_vertex(n))
311 118560 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n])};
312 182400 : else if (this->is_edge(n))
313 131328 : return {n - this->n_vertices()};
314 :
315 4256 : libmesh_assert(this->is_face(n) || this->is_internal(n));
316 46816 : return {};
317 : }
318 :
319 :
320 :
321 7568002 : bool Prism::on_reference_element(const Point & p,
322 : const Real eps) const
323 : {
324 744121 : const Real & xi = p(0);
325 744121 : const Real & eta = p(1);
326 744121 : const Real & zeta = p(2);
327 :
328 : // inside the reference triangle with zeta in [-1,1]
329 13695153 : return ((xi >= 0.-eps) &&
330 6127151 : (eta >= 0.-eps) &&
331 4044507 : (zeta >= -1.-eps) &&
332 11945364 : (zeta <= 1.+eps) &&
333 4242538 : ((xi + eta) <= 1.+eps));
334 : }
335 :
336 :
337 :
338 4800 : std::pair<Real, Real> Prism::qual_bounds (const ElemQuality) const
339 : {
340 : // No suggested ranges have yet been established for any of the
341 : // metrics Quality::valid() lists for prisms, so we always return
342 : // the (-1, -1) sentinel.
343 : libmesh_warning("Invalid quality measure chosen.");
344 4800 : return std::make_pair(Real(-1), Real(-1));
345 : }
346 :
347 :
348 :
349 : const unsigned short int Prism::_second_order_vertex_child_number[18] =
350 : {
351 : 99,99,99,99,99,99, // Vertices
352 : 0,1,0,0,1,2,3,4,3, // Edges
353 : 0,1,0 // Faces
354 : };
355 :
356 :
357 :
358 : const unsigned short int Prism::_second_order_vertex_child_index[18] =
359 : {
360 : 99,99,99,99,99,99, // Vertices
361 : 1,2,2,3,4,5,4,5,5, // Edges
362 : 4,5,5 // Faces
363 : };
364 :
365 :
366 : const unsigned short int Prism::_second_order_adjacent_vertices[9][2] =
367 : {
368 : { 0, 1}, // vertices adjacent to node 6
369 : { 1, 2}, // vertices adjacent to node 7
370 : { 0, 2}, // vertices adjacent to node 8
371 :
372 : { 0, 3}, // vertices adjacent to node 9
373 : { 1, 4}, // vertices adjacent to node 10
374 : { 2, 5}, // vertices adjacent to node 11
375 :
376 : { 3, 4}, // vertices adjacent to node 12
377 : { 4, 5}, // vertices adjacent to node 13
378 : { 3, 5} // vertices adjacent to node 14
379 : };
380 :
381 : } // namespace libMesh
|
