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_pyramid.h"
23 : #include "libmesh/cell_pyramid5.h"
24 : #include "libmesh/enum_elem_quality.h"
25 : #include "libmesh/face_tri3.h"
26 : #include "libmesh/face_quad4.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 : // ------------------------------------------------------------
32 : // Pyramid class static member initializations
33 : const int Pyramid::num_sides;
34 : const int Pyramid::num_edges;
35 : const int Pyramid::num_children;
36 :
37 : const Real Pyramid::_master_points[18][3] =
38 : {
39 : {-1, -1, 0},
40 : {1, -1, 0},
41 : {1, 1, 0},
42 : {-1, 1, 0},
43 : {0, 0, 1},
44 : {0, -1, 0},
45 : {1, 0, 0},
46 : {0, 1, 0},
47 : {-1, 0, 0},
48 : {-0.5, -0.5, 0.5},
49 : {0.5, -0.5, 0.5},
50 : {0.5, 0.5, 0.5},
51 : {-0.5, 0.5, 0.5},
52 : {0, 0, 0},
53 : {0, -Real(2)/3, Real(1)/3},
54 : { Real(2)/3, 0, Real(1)/3},
55 : {0, Real(2)/3, Real(1)/3},
56 : {-Real(2)/3, 0, Real(1)/3}
57 : };
58 :
59 : const unsigned int Pyramid::edge_sides_map[8][2] =
60 : {
61 : {0, 4}, // Edge 0
62 : {1, 4}, // Edge 1
63 : {2, 4}, // Edge 2
64 : {3, 4}, // Edge 3
65 : {0, 3}, // Edge 4
66 : {0, 1}, // Edge 5
67 : {1, 2}, // Edge 6
68 : {2, 3} // Edge 7
69 : };
70 :
71 : const unsigned int Pyramid::adjacent_edges_map[/*num_vertices*/5][/*max_adjacent_edges*/4] =
72 : {
73 : {0, 3, 4, 99}, // Edges adjacent to node 0
74 : {0, 1, 5, 99}, // Edges adjacent to node 1
75 : {1, 2, 6, 99}, // Edges adjacent to node 2
76 : {2, 3, 7, 99}, // Edges adjacent to node 3
77 : {4, 5, 6, 7} // Edges adjacent to node 4
78 : };
79 :
80 : // ------------------------------------------------------------
81 : // Pyramid class member functions
82 0 : dof_id_type Pyramid::key (const unsigned int s) const
83 : {
84 0 : libmesh_assert_less (s, this->n_sides());
85 :
86 0 : switch (s)
87 : {
88 0 : case 0: // triangular face 1
89 : case 1: // triangular face 2
90 : case 2: // triangular face 3
91 : case 3: // triangular face 4
92 0 : return this->compute_key (this->node_id(Pyramid5::side_nodes_map[s][0]),
93 0 : this->node_id(Pyramid5::side_nodes_map[s][1]),
94 0 : this->node_id(Pyramid5::side_nodes_map[s][2]));
95 :
96 0 : case 4: // the quad face at z=0
97 0 : return this->compute_key (this->node_id(Pyramid5::side_nodes_map[s][0]),
98 0 : this->node_id(Pyramid5::side_nodes_map[s][1]),
99 0 : this->node_id(Pyramid5::side_nodes_map[s][2]),
100 0 : this->node_id(Pyramid5::side_nodes_map[s][3]));
101 :
102 0 : default:
103 0 : libmesh_error_msg("Invalid side s = " << s);
104 : }
105 : }
106 :
107 :
108 :
109 6270484 : dof_id_type Pyramid::low_order_key (const unsigned int s) const
110 : {
111 139468 : libmesh_assert_less (s, this->n_sides());
112 :
113 6270484 : switch (s)
114 : {
115 5052488 : case 0: // triangular face 1
116 : case 1: // triangular face 2
117 : case 2: // triangular face 3
118 : case 3: // triangular face 4
119 5164088 : return this->compute_key (this->node_id(Pyramid5::side_nodes_map[s][0]),
120 5052488 : this->node_id(Pyramid5::side_nodes_map[s][1]),
121 5164088 : this->node_id(Pyramid5::side_nodes_map[s][2]));
122 :
123 1217996 : case 4: // the quad face at z=0
124 1245864 : return this->compute_key (this->node_id(Pyramid5::side_nodes_map[s][0]),
125 1217996 : this->node_id(Pyramid5::side_nodes_map[s][1]),
126 1217996 : this->node_id(Pyramid5::side_nodes_map[s][2]),
127 1245864 : this->node_id(Pyramid5::side_nodes_map[s][3]));
128 :
129 0 : default:
130 0 : libmesh_error_msg("Invalid side s = " << s);
131 : }
132 : }
133 :
134 :
135 :
136 25486 : unsigned int Pyramid::local_side_node(unsigned int side,
137 : unsigned int side_node) const
138 : {
139 2116 : libmesh_assert_less (side, this->n_sides());
140 :
141 : // Never more than 4 nodes per side.
142 2116 : libmesh_assert_less(side_node, Pyramid5::nodes_per_side);
143 :
144 : // Some sides have 3 nodes.
145 2116 : libmesh_assert(side == 4 || side_node < 3);
146 :
147 25484 : return Pyramid5::side_nodes_map[side][side_node];
148 : }
149 :
150 :
151 :
152 665379 : unsigned int Pyramid::local_edge_node(unsigned int edge,
153 : unsigned int edge_node) const
154 : {
155 55370 : libmesh_assert_less(edge, this->n_edges());
156 55370 : libmesh_assert_less(edge_node, Pyramid5::nodes_per_edge);
157 :
158 665379 : return Pyramid5::edge_nodes_map[edge][edge_node];
159 : }
160 :
161 :
162 :
163 1286756 : std::unique_ptr<Elem> Pyramid::side_ptr (const unsigned int i)
164 : {
165 29132 : libmesh_assert_less (i, this->n_sides());
166 :
167 : // Return value
168 1286756 : std::unique_ptr<Elem> face;
169 :
170 : // Set up the type of element
171 1286756 : switch (i)
172 : {
173 665599 : case 0: // triangular face 1
174 : case 1: // triangular face 2
175 : case 2: // triangular face 3
176 : case 3: // triangular face 4
177 : {
178 665599 : face = std::make_unique<Tri3>();
179 665599 : break;
180 : }
181 621157 : case 4: // the quad face at z=0
182 : {
183 621157 : face = std::make_unique<Quad4>();
184 621157 : break;
185 : }
186 0 : default:
187 0 : libmesh_error_msg("Invalid side i = " << i);
188 : }
189 :
190 : // Set the nodes
191 5768181 : for (auto n : face->node_index_range())
192 4582663 : face->set_node(n, this->node_ptr(Pyramid5::side_nodes_map[i][n]));
193 :
194 1286756 : return face;
195 0 : }
196 :
197 :
198 :
199 5317032 : void Pyramid::side_ptr (std::unique_ptr<Elem> & side,
200 : const unsigned int i)
201 : {
202 125892 : libmesh_assert_less (i, this->n_sides());
203 :
204 5317032 : switch (i)
205 : {
206 111952 : case 0: // triangular face 1
207 : case 1: // triangular face 2
208 : case 2: // triangular face 3
209 : case 3: // triangular face 4
210 : {
211 4683828 : if (!side.get() || side->type() != TRI3)
212 : {
213 1298962 : side = this->side_ptr(i);
214 664747 : return;
215 : }
216 96686 : break;
217 : }
218 :
219 13940 : case 4: // the quad face at z=0
220 : {
221 633204 : if (!side.get() || side->type() != QUAD4)
222 : {
223 1214216 : side = this->side_ptr(i);
224 620944 : return;
225 : }
226 104 : break;
227 : }
228 :
229 0 : default:
230 0 : libmesh_error_msg("Invalid side i = " << i);
231 : }
232 :
233 4128131 : side->subdomain_id() = this->subdomain_id();
234 :
235 : // Set the nodes
236 16137624 : for (auto n : side->node_index_range())
237 12396757 : side->set_node(n, this->node_ptr(Pyramid5::side_nodes_map[i][n]));
238 : }
239 :
240 :
241 :
242 0 : bool Pyramid::is_child_on_side(const unsigned int c,
243 : const unsigned int s) const
244 : {
245 0 : libmesh_assert_less (c, this->n_children());
246 0 : libmesh_assert_less (s, this->n_sides());
247 :
248 0 : for (unsigned int i = 0; i != 4; ++i)
249 0 : if (Pyramid5::side_nodes_map[s][i] == c)
250 0 : return true;
251 0 : return false;
252 : }
253 :
254 :
255 :
256 36864 : bool Pyramid::is_edge_on_side(const unsigned int e,
257 : const unsigned int s) const
258 : {
259 3072 : libmesh_assert_less (e, this->n_edges());
260 3072 : libmesh_assert_less (s, this->n_sides());
261 :
262 36864 : return (edge_sides_map[e][0] == s || edge_sides_map[e][1] == s);
263 : }
264 :
265 :
266 :
267 258048 : std::vector<unsigned int> Pyramid::sides_on_edge(const unsigned int e) const
268 : {
269 21504 : libmesh_assert_less (e, this->n_edges());
270 258048 : return {edge_sides_map[e][0], edge_sides_map[e][1]};
271 : }
272 :
273 :
274 : bool
275 22431 : Pyramid::is_flipped() const
276 : {
277 21087 : return (triple_product(this->point(1)-this->point(0),
278 21087 : this->point(3)-this->point(0),
279 25119 : this->point(4)-this->point(0)) < 0);
280 : }
281 :
282 : std::vector<unsigned int>
283 432000 : Pyramid::edges_adjacent_to_node(const unsigned int n) const
284 : {
285 36000 : libmesh_assert_less(n, this->n_nodes());
286 432000 : if (this->is_vertex(n))
287 : {
288 172800 : auto trim = (n < 4) ? 1 : 0;
289 172800 : return {std::begin(adjacent_edges_map[n]), std::end(adjacent_edges_map[n]) - trim};
290 : }
291 259200 : else if (this->is_edge(n))
292 207360 : return {n - this->n_vertices()};
293 :
294 : // Not a vertex or edge node, so must be one of the face nodes.
295 4320 : libmesh_assert(this->is_face(n));
296 47520 : return {};
297 : }
298 :
299 675138 : unsigned int Pyramid::local_singular_node(const Point & p, const Real tol) const
300 : {
301 675138 : return this->node_ref(4).absolute_fuzzy_equals(p, tol) ? 4 : invalid_uint;
302 : }
303 :
304 :
305 4766180 : bool Pyramid::on_reference_element(const Point & p,
306 : const Real eps) const
307 : {
308 131100 : const Real & xi = p(0);
309 131100 : const Real & eta = p(1);
310 131100 : const Real & zeta = p(2);
311 :
312 : // Check that the point is on the same side of all the faces
313 : // by testing whether:
314 : //
315 : // n_i.(x - x_i) <= 0
316 : //
317 : // for each i, where:
318 : // n_i is the outward normal of face i,
319 : // x_i is a point on face i.
320 7920228 : return ((-eta - 1. + zeta <= 0.+eps) &&
321 3154048 : ( xi - 1. + zeta <= 0.+eps) &&
322 1712699 : ( eta - 1. + zeta <= 0.+eps) &&
323 6178471 : ( -xi - 1. + zeta <= 0.+eps) &&
324 1310900 : ( zeta >= 0.-eps));
325 : }
326 :
327 :
328 :
329 11520 : std::pair<Real, Real> Pyramid::qual_bounds (const ElemQuality) const
330 : {
331 : // No suggested ranges have yet been established for any of the
332 : // metrics Quality::valid() lists for pyramids, so we always return
333 : // the (-1, -1) sentinel.
334 : libmesh_warning("Invalid quality measure chosen.");
335 11520 : return std::make_pair(Real(-1), Real(-1));
336 : }
337 :
338 :
339 : } // namespace libMesh
|
