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