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 : // 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[18][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 : };
58 :
59 : const unsigned int Prism::edge_sides_map[9][2] =
60 : {
61 : {0, 1}, // Edge 0
62 : {0, 2}, // Edge 1
63 : {0, 3}, // Edge 2
64 : {1, 3}, // Edge 3
65 : {1, 2}, // Edge 4
66 : {2, 3}, // Edge 5
67 : {1, 4}, // Edge 6
68 : {2, 4}, // Edge 7
69 : {3, 4} // Edge 8
70 : };
71 :
72 : const unsigned int Prism::adjacent_edges_map[/*num_vertices*/6][/*n_adjacent_edges*/3] =
73 : {
74 : {0, 2, 3}, // Edges adjacent to node 0
75 : {0, 1, 4}, // Edges adjacent to node 1
76 : {1, 2, 5}, // Edges adjacent to node 2
77 : {3, 6, 8}, // Edges adjacent to node 3
78 : {4, 6, 7}, // Edges adjacent to node 4
79 : {5, 7, 8}, // Edges adjacent to node 5
80 : };
81 :
82 :
83 : // ------------------------------------------------------------
84 : // Prism class member functions
85 0 : dof_id_type Prism::key (const unsigned int s) const
86 : {
87 0 : libmesh_assert_less (s, this->n_sides());
88 :
89 0 : switch (s)
90 : {
91 0 : case 0: // the triangular face at z=0
92 : case 4: // the triangular face at z=1
93 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
94 0 : this->node_id(Prism6::side_nodes_map[s][1]),
95 0 : this->node_id(Prism6::side_nodes_map[s][2]));
96 :
97 0 : case 1: // the quad face at y=0
98 : case 2: // the other quad face
99 : case 3: // the quad face at x=0
100 0 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
101 0 : this->node_id(Prism6::side_nodes_map[s][1]),
102 0 : this->node_id(Prism6::side_nodes_map[s][2]),
103 0 : this->node_id(Prism6::side_nodes_map[s][3]));
104 :
105 0 : default:
106 0 : libmesh_error_msg("Invalid side " << s);
107 : }
108 : }
109 :
110 :
111 :
112 4348745 : dof_id_type Prism::low_order_key (const unsigned int s) const
113 : {
114 217220 : libmesh_assert_less (s, this->n_sides());
115 :
116 4348745 : switch (s)
117 : {
118 1739498 : case 0: // the triangular face at z=0
119 : case 4: // the triangular face at z=1
120 1826386 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
121 1739498 : this->node_id(Prism6::side_nodes_map[s][1]),
122 1826386 : this->node_id(Prism6::side_nodes_map[s][2]));
123 :
124 2609247 : case 1: // the quad face at y=0
125 : case 2: // the other quad face
126 : case 3: // the quad face at x=0
127 2739579 : return this->compute_key (this->node_id(Prism6::side_nodes_map[s][0]),
128 2609247 : this->node_id(Prism6::side_nodes_map[s][1]),
129 2609247 : this->node_id(Prism6::side_nodes_map[s][2]),
130 2739579 : this->node_id(Prism6::side_nodes_map[s][3]));
131 :
132 0 : default:
133 0 : libmesh_error_msg("Invalid side " << s);
134 : }
135 : }
136 :
137 :
138 :
139 91870 : unsigned int Prism::local_side_node(unsigned int side,
140 : unsigned int side_node) const
141 : {
142 7648 : libmesh_assert_less (side, this->n_sides());
143 :
144 : // Never more than 4 nodes per side.
145 7648 : libmesh_assert_less(side_node, Prism6::nodes_per_side);
146 :
147 : // Some sides have 3 nodes.
148 7648 : libmesh_assert(!(side==0 || side==4) || side_node < 3);
149 :
150 91868 : return Prism6::side_nodes_map[side][side_node];
151 : }
152 :
153 :
154 :
155 31575726 : unsigned int Prism::local_edge_node(unsigned int edge,
156 : unsigned int edge_node) const
157 : {
158 2816280 : libmesh_assert_less(edge, this->n_edges());
159 2816280 : libmesh_assert_less(edge_node, Prism6::nodes_per_edge);
160 :
161 31575726 : return Prism6::edge_nodes_map[edge][edge_node];
162 : }
163 :
164 :
165 :
166 264791983 : std::unique_ptr<Elem> Prism::side_ptr (const unsigned int i)
167 : {
168 23560998 : libmesh_assert_less (i, this->n_sides());
169 :
170 264791983 : std::unique_ptr<Elem> face;
171 :
172 : // Set up the type of element
173 264791983 : switch (i)
174 : {
175 819381 : case 0: // the triangular face at z=0
176 : case 4: // the triangular face at z=1
177 : {
178 819381 : face = std::make_unique<Tri3>();
179 819381 : break;
180 : }
181 263972602 : case 1: // the quad face at y=0
182 : case 2: // the other quad face
183 : case 3: // the quad face at x=0
184 : {
185 263972602 : face = std::make_unique<Quad4>();
186 263972602 : break;
187 : }
188 0 : default:
189 0 : libmesh_error_msg("Invalid side i = " << i);
190 : }
191 :
192 : // Set the nodes
193 1323140534 : for (auto n : face->node_index_range())
194 1152542949 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
195 :
196 264791983 : return face;
197 0 : }
198 :
199 :
200 :
201 3123892 : void Prism::side_ptr (std::unique_ptr<Elem> & side,
202 : const unsigned int i)
203 : {
204 180456 : libmesh_assert_less (i, this->n_sides());
205 :
206 3123892 : switch (i)
207 : {
208 : // the base face
209 69128 : case 0: // the triangular face at z=0
210 : case 4: // the triangular face at z=1
211 : {
212 1141932 : if (!side.get() || side->type() != TRI3)
213 : {
214 1523572 : side = this->side_ptr(i);
215 810674 : return;
216 : }
217 20240 : break;
218 : }
219 :
220 111328 : case 1: // the quad face at y=0
221 : case 2: // the other quad face
222 : case 3: // the quad face at x=0
223 : {
224 1981960 : if (!side.get() || side->type() != QUAD4)
225 : {
226 1646890 : side = this->side_ptr(i);
227 873979 : return;
228 : }
229 60794 : break;
230 : }
231 :
232 0 : default:
233 0 : libmesh_error_msg("Invalid side i = " << i);
234 : }
235 :
236 1520273 : side->subdomain_id() = this->subdomain_id();
237 :
238 : // Set the nodes
239 6864937 : for (auto n : side->node_index_range())
240 5729594 : side->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
241 : }
242 :
243 :
244 :
245 1248672 : bool Prism::is_child_on_side(const unsigned int c,
246 : const unsigned int s) const
247 : {
248 550492 : libmesh_assert_less (c, this->n_children());
249 550492 : libmesh_assert_less (s, this->n_sides());
250 :
251 4770226 : for (unsigned int i = 0; i != 4; ++i)
252 4110115 : if (Prism6::side_elems_map[s][i] == c)
253 243224 : return true;
254 307268 : return false;
255 : }
256 :
257 :
258 :
259 3390736 : bool Prism::is_edge_on_side(const unsigned int e,
260 : const unsigned int s) const
261 : {
262 2912302 : libmesh_assert_less (e, this->n_edges());
263 2912302 : libmesh_assert_less (s, this->n_sides());
264 :
265 3390736 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
266 : }
267 :
268 :
269 :
270 34560 : std::vector<unsigned int> Prism::sides_on_edge(const unsigned int e) const
271 : {
272 2880 : libmesh_assert_less(e, this->n_edges());
273 :
274 34560 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
275 : }
276 :
277 :
278 :
279 0 : unsigned int Prism::opposite_side(const unsigned int side_in) const
280 : {
281 0 : libmesh_assert_less (side_in, 5);
282 : static const unsigned int prism_opposites[5] =
283 : {4, invalid_uint, invalid_uint, invalid_uint, 0};
284 0 : return prism_opposites[side_in];
285 : }
286 :
287 :
288 :
289 : bool
290 9428 : Prism::is_flipped() const
291 : {
292 8868 : return (triple_product(this->point(1)-this->point(0),
293 8868 : this->point(2)-this->point(0),
294 10548 : this->point(3)-this->point(0)) < 0);
295 : }
296 :
297 :
298 : std::vector<unsigned int>
299 76800 : Prism::edges_adjacent_to_node(const unsigned int n) const
300 : {
301 6400 : libmesh_assert_less(n, this->n_nodes());
302 :
303 : // For vertices, we use the Prism::adjacent_sides_map, otherwise each
304 : // of the mid-edge nodes is adjacent only to the edge it is on, and
305 : // face/internal nodes are not adjacent to any edge.
306 76800 : if (this->is_vertex(n))
307 31200 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n])};
308 48000 : else if (this->is_edge(n))
309 34560 : return {n - this->n_vertices()};
310 :
311 1120 : libmesh_assert(this->is_face(n) || this->is_internal(n));
312 12320 : return {};
313 : }
314 :
315 :
316 :
317 7028703 : bool Prism::on_reference_element(const Point & p,
318 : const Real eps) const
319 : {
320 720231 : const Real & xi = p(0);
321 720231 : const Real & eta = p(1);
322 720231 : const Real & zeta = p(2);
323 :
324 : // inside the reference triangle with zeta in [-1,1]
325 12818967 : return ((xi >= 0.-eps) &&
326 5790264 : (eta >= 0.-eps) &&
327 3752311 : (zeta >= -1.-eps) &&
328 11111463 : (zeta <= 1.+eps) &&
329 3981476 : ((xi + eta) <= 1.+eps));
330 : }
331 :
332 :
333 :
334 : const unsigned short int Prism::_second_order_vertex_child_number[18] =
335 : {
336 : 99,99,99,99,99,99, // Vertices
337 : 0,1,0,0,1,2,3,4,3, // Edges
338 : 0,1,0 // Faces
339 : };
340 :
341 :
342 :
343 : const unsigned short int Prism::_second_order_vertex_child_index[18] =
344 : {
345 : 99,99,99,99,99,99, // Vertices
346 : 1,2,2,3,4,5,4,5,5, // Edges
347 : 4,5,5 // Faces
348 : };
349 :
350 :
351 : const unsigned short int Prism::_second_order_adjacent_vertices[9][2] =
352 : {
353 : { 0, 1}, // vertices adjacent to node 6
354 : { 1, 2}, // vertices adjacent to node 7
355 : { 0, 2}, // vertices adjacent to node 8
356 :
357 : { 0, 3}, // vertices adjacent to node 9
358 : { 1, 4}, // vertices adjacent to node 10
359 : { 2, 5}, // vertices adjacent to node 11
360 :
361 : { 3, 4}, // vertices adjacent to node 12
362 : { 4, 5}, // vertices adjacent to node 13
363 : { 3, 5} // vertices adjacent to node 14
364 : };
365 :
366 : } // namespace libMesh
|
