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