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/face_quad4.h"
25 : #include "libmesh/face_tri3.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 :
31 : // ------------------------------------------------------------
32 : // Prism class static member initializations
33 : const int Prism::num_sides;
34 : const int Prism::num_edges;
35 : const int Prism::num_children;
36 :
37 : const Real Prism::_master_points[21][3] =
38 : {
39 : {0, 0, -1},
40 : {1, 0, -1},
41 : {0, 1, -1},
42 : {0, 0, 1},
43 : {1, 0, 1},
44 : {0, 1, 1},
45 : {0.5, 0, -1},
46 : {0.5, 0.5, -1},
47 : {0, 0.5, -1},
48 : {0, 0, 0},
49 : {1, 0, 0},
50 : {0, 1, 0},
51 : {0.5, 0, 1},
52 : {0.5, 0.5, 1},
53 : {0, 0.5, 1},
54 : {0.5, 0, 0},
55 : {0.5, 0.5, 0},
56 : {0, 0.5, 0},
57 : {Real(1)/3, Real(1)/3, -1},
58 : {Real(1)/3, Real(1)/3, 1},
59 : {Real(1)/3, Real(1)/3, 0}
60 : };
61 :
62 : const unsigned int Prism::edge_sides_map[9][2] =
63 : {
64 : {0, 1}, // Edge 0
65 : {0, 2}, // Edge 1
66 : {0, 3}, // Edge 2
67 : {1, 3}, // Edge 3
68 : {1, 2}, // Edge 4
69 : {2, 3}, // Edge 5
70 : {1, 4}, // Edge 6
71 : {2, 4}, // Edge 7
72 : {3, 4} // Edge 8
73 : };
74 :
75 : const unsigned int Prism::adjacent_edges_map[/*num_vertices*/6][/*n_adjacent_edges*/3] =
76 : {
77 : {0, 2, 3}, // Edges adjacent to node 0
78 : {0, 1, 4}, // Edges adjacent to node 1
79 : {1, 2, 5}, // Edges adjacent to node 2
80 : {3, 6, 8}, // Edges adjacent to node 3
81 : {4, 6, 7}, // Edges adjacent to node 4
82 : {5, 7, 8}, // Edges adjacent to node 5
83 : };
84 :
85 :
86 : // ------------------------------------------------------------
87 : // Prism class member functions
88 0 : dof_id_type Prism::key (const unsigned int s) const
89 : {
90 0 : libmesh_assert_less (s, this->n_sides());
91 :
92 0 : switch (s)
93 : {
94 0 : case 0: // the triangular face at z=0
95 : case 4: // the triangular face at z=1
96 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
97 0 : this->node_id(Prism6::side_nodes_map[s][1]),
98 0 : this->node_id(Prism6::side_nodes_map[s][2]));
99 :
100 0 : case 1: // the quad face at y=0
101 : case 2: // the other quad face
102 : case 3: // the quad face at x=0
103 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
104 0 : this->node_id(Prism6::side_nodes_map[s][1]),
105 0 : this->node_id(Prism6::side_nodes_map[s][2]),
106 0 : this->node_id(Prism6::side_nodes_map[s][3]));
107 :
108 0 : default:
109 0 : libmesh_error_msg("Invalid side " << s);
110 : }
111 : }
112 :
113 :
114 :
115 5008295 : dof_id_type Prism::low_order_key (const unsigned int s) const
116 : {
117 230120 : libmesh_assert_less (s, this->n_sides());
118 :
119 5008295 : switch (s)
120 : {
121 2003318 : case 0: // the triangular face at z=0
122 : case 4: // the triangular face at z=1
123 2095366 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
124 2003318 : this->node_id(Prism6::side_nodes_map[s][1]),
125 2095366 : this->node_id(Prism6::side_nodes_map[s][2]));
126 :
127 3004977 : case 1: // the quad face at y=0
128 : case 2: // the other quad face
129 : case 3: // the quad face at x=0
130 3143049 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
131 3004977 : this->node_id(Prism6::side_nodes_map[s][1]),
132 3004977 : this->node_id(Prism6::side_nodes_map[s][2]),
133 3143049 : this->node_id(Prism6::side_nodes_map[s][3]));
134 :
135 0 : default:
136 0 : libmesh_error_msg("Invalid side " << s);
137 : }
138 : }
139 :
140 :
141 :
142 91870 : unsigned int Prism::local_side_node(unsigned int side,
143 : unsigned int side_node) const
144 : {
145 7648 : libmesh_assert_less (side, this->n_sides());
146 :
147 : // Never more than 4 nodes per side.
148 7648 : libmesh_assert_less(side_node, Prism6::nodes_per_side);
149 :
150 : // Some sides have 3 nodes.
151 7648 : libmesh_assert(!(side==0 || side==4) || side_node < 3);
152 :
153 91868 : return Prism6::side_nodes_map[side][side_node];
154 : }
155 :
156 :
157 :
158 31621086 : unsigned int Prism::local_edge_node(unsigned int edge,
159 : unsigned int edge_node) const
160 : {
161 2820060 : libmesh_assert_less(edge, this->n_edges());
162 2820060 : libmesh_assert_less(edge_node, Prism6::nodes_per_edge);
163 :
164 31621086 : return Prism6::edge_nodes_map[edge][edge_node];
165 : }
166 :
167 :
168 :
169 1947193 : std::unique_ptr<Elem> Prism::side_ptr (const unsigned int i)
170 : {
171 104874 : libmesh_assert_less (i, this->n_sides());
172 :
173 1947193 : std::unique_ptr<Elem> face;
174 :
175 : // Set up the type of element
176 1947193 : switch (i)
177 : {
178 934693 : case 0: // the triangular face at z=0
179 : case 4: // the triangular face at z=1
180 : {
181 934693 : face = std::make_unique<Tri3>();
182 934693 : break;
183 : }
184 1012500 : case 1: // the quad face at y=0
185 : case 2: // the other quad face
186 : case 3: // the quad face at x=0
187 : {
188 1012500 : face = std::make_unique<Quad4>();
189 1012500 : break;
190 : }
191 0 : default:
192 0 : libmesh_error_msg("Invalid side i = " << i);
193 : }
194 :
195 : // Set the nodes
196 8801272 : for (auto n : face->node_index_range())
197 7222177 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
198 :
199 1947193 : return face;
200 0 : }
201 :
202 :
203 :
204 3534461 : void Prism::side_ptr (std::unique_ptr<Elem> & side,
205 : const unsigned int i)
206 : {
207 188614 : libmesh_assert_less (i, this->n_sides());
208 :
209 3534461 : switch (i)
210 : {
211 : // the base face
212 71912 : case 0: // the triangular face at z=0
213 : case 4: // the triangular face at z=1
214 : {
215 1281084 : if (!side.get() || side->type() != TRI3)
216 : {
217 1765204 : side = this->side_ptr(i);
218 933978 : return;
219 : }
220 20536 : break;
221 : }
222 :
223 116702 : case 1: // the quad face at y=0
224 : case 2: // the other quad face
225 : case 3: // the quad face at x=0
226 : {
227 2253377 : if (!side.get() || side->type() != QUAD4)
228 : {
229 1915978 : side = this->side_ptr(i);
230 1011435 : return;
231 : }
232 63256 : break;
233 : }
234 :
235 0 : default:
236 0 : libmesh_error_msg("Invalid side i = " << i);
237 : }
238 :
239 1672840 : side->subdomain_id() = this->subdomain_id();
240 :
241 : // Set the nodes
242 7598134 : for (auto n : side->node_index_range())
243 6323718 : side->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
244 : }
245 :
246 :
247 :
248 9050449 : bool Prism::is_child_on_side(const unsigned int c,
249 : const unsigned int s) const
250 : {
251 547932 : libmesh_assert_less (c, this->n_children());
252 547932 : libmesh_assert_less (s, this->n_sides());
253 :
254 33913158 : for (unsigned int i = 0; i != 4; ++i)
255 29396451 : if (Prism6::side_elems_map[s][i] == c)
256 242584 : return true;
257 305348 : return false;
258 : }
259 :
260 :
261 :
262 3364720 : bool Prism::is_edge_on_side(const unsigned int e,
263 : const unsigned int s) const
264 : {
265 2886286 : libmesh_assert_less (e, this->n_edges());
266 2886286 : libmesh_assert_less (s, this->n_sides());
267 :
268 3364720 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
269 : }
270 :
271 :
272 :
273 34560 : std::vector<unsigned int> Prism::sides_on_edge(const unsigned int e) const
274 : {
275 2880 : libmesh_assert_less(e, this->n_edges());
276 :
277 34560 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
278 : }
279 :
280 :
281 :
282 0 : unsigned int Prism::opposite_side(const unsigned int side_in) const
283 : {
284 0 : libmesh_assert_less (side_in, 5);
285 : static const unsigned int prism_opposites[5] =
286 : {4, invalid_uint, invalid_uint, invalid_uint, 0};
287 0 : return prism_opposites[side_in];
288 : }
289 :
290 :
291 :
292 : bool
293 9428 : Prism::is_flipped() const
294 : {
295 8868 : return (triple_product(this->point(1)-this->point(0),
296 8868 : this->point(2)-this->point(0),
297 10548 : this->point(3)-this->point(0)) < 0);
298 : }
299 :
300 :
301 : std::vector<unsigned int>
302 76800 : Prism::edges_adjacent_to_node(const unsigned int n) const
303 : {
304 6400 : libmesh_assert_less(n, this->n_nodes());
305 :
306 : // For vertices, we use the Prism::adjacent_sides_map, otherwise each
307 : // of the mid-edge nodes is adjacent only to the edge it is on, and
308 : // face/internal nodes are not adjacent to any edge.
309 76800 : if (this->is_vertex(n))
310 31200 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n])};
311 48000 : else if (this->is_edge(n))
312 34560 : return {n - this->n_vertices()};
313 :
314 1120 : libmesh_assert(this->is_face(n) || this->is_internal(n));
315 12320 : return {};
316 : }
317 :
318 :
319 :
320 7558846 : bool Prism::on_reference_element(const Point & p,
321 : const Real eps) const
322 : {
323 734857 : const Real & xi = p(0);
324 734857 : const Real & eta = p(1);
325 734857 : const Real & zeta = p(2);
326 :
327 : // inside the reference triangle with zeta in [-1,1]
328 13676841 : return ((xi >= 0.-eps) &&
329 6117995 : (eta >= 0.-eps) &&
330 4035351 : (zeta >= -1.-eps) &&
331 11917788 : (zeta <= 1.+eps) &&
332 4224118 : ((xi + eta) <= 1.+eps));
333 : }
334 :
335 :
336 :
337 : const unsigned short int Prism::_second_order_vertex_child_number[18] =
338 : {
339 : 99,99,99,99,99,99, // Vertices
340 : 0,1,0,0,1,2,3,4,3, // Edges
341 : 0,1,0 // Faces
342 : };
343 :
344 :
345 :
346 : const unsigned short int Prism::_second_order_vertex_child_index[18] =
347 : {
348 : 99,99,99,99,99,99, // Vertices
349 : 1,2,2,3,4,5,4,5,5, // Edges
350 : 4,5,5 // Faces
351 : };
352 :
353 :
354 : const unsigned short int Prism::_second_order_adjacent_vertices[9][2] =
355 : {
356 : { 0, 1}, // vertices adjacent to node 6
357 : { 1, 2}, // vertices adjacent to node 7
358 : { 0, 2}, // vertices adjacent to node 8
359 :
360 : { 0, 3}, // vertices adjacent to node 9
361 : { 1, 4}, // vertices adjacent to node 10
362 : { 2, 5}, // vertices adjacent to node 11
363 :
364 : { 3, 4}, // vertices adjacent to node 12
365 : { 4, 5}, // vertices adjacent to node 13
366 : { 3, 5} // vertices adjacent to node 14
367 : };
368 :
369 : } // namespace libMesh
|
