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 : // Local includes
20 : #include "libmesh/cell_pyramid5.h"
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/face_tri3.h"
23 : #include "libmesh/face_quad4.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 : #include "libmesh/cell_hex8.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 :
32 :
33 :
34 : // ------------------------------------------------------------
35 : // Pyramid5 class static member initializations
36 : const int Pyramid5::num_nodes;
37 : const int Pyramid5::nodes_per_side;
38 : const int Pyramid5::nodes_per_edge;
39 :
40 : const unsigned int Pyramid5::side_nodes_map[Pyramid5::num_sides][Pyramid5::nodes_per_side] =
41 : {
42 : {0, 1, 4, 99}, // Side 0
43 : {1, 2, 4, 99}, // Side 1
44 : {2, 3, 4, 99}, // Side 2
45 : {3, 0, 4, 99}, // Side 3
46 : {0, 3, 2, 1} // Side 4
47 : };
48 :
49 : const unsigned int Pyramid5::edge_nodes_map[Pyramid5::num_edges][Pyramid5::nodes_per_edge] =
50 : {
51 : {0, 1}, // Edge 0
52 : {1, 2}, // Edge 1
53 : {2, 3}, // Edge 2
54 : {0, 3}, // Edge 3
55 : {0, 4}, // Edge 4
56 : {1, 4}, // Edge 5
57 : {2, 4}, // Edge 6
58 : {3, 4} // Edge 7
59 : };
60 :
61 : // ------------------------------------------------------------
62 : // Pyramid5 class member functions
63 :
64 42120 : bool Pyramid5::is_vertex(const unsigned int libmesh_dbg_var(n)) const
65 : {
66 2520 : libmesh_assert_not_equal_to (n, invalid_uint);
67 42120 : return true;
68 : }
69 :
70 0 : bool Pyramid5::is_edge(const unsigned int) const
71 : {
72 0 : return false;
73 : }
74 :
75 0 : bool Pyramid5::is_face(const unsigned int) const
76 : {
77 0 : return false;
78 : }
79 :
80 29423 : bool Pyramid5::is_node_on_side(const unsigned int n,
81 : const unsigned int s) const
82 : {
83 2354 : libmesh_assert_less (s, n_sides());
84 2354 : return std::find(std::begin(side_nodes_map[s]),
85 29423 : std::end(side_nodes_map[s]),
86 29423 : n) != std::end(side_nodes_map[s]);
87 : }
88 :
89 : std::vector<unsigned>
90 34891 : Pyramid5::nodes_on_side(const unsigned int s) const
91 : {
92 2794 : libmesh_assert_less(s, n_sides());
93 34891 : auto trim = (s == 4) ? 0 : 1;
94 34891 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
95 : }
96 :
97 : std::vector<unsigned>
98 23608 : Pyramid5::nodes_on_edge(const unsigned int e) const
99 : {
100 1936 : libmesh_assert_less(e, n_edges());
101 25544 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
102 : }
103 :
104 12056 : bool Pyramid5::is_node_on_edge(const unsigned int n,
105 : const unsigned int e) const
106 : {
107 848 : libmesh_assert_less (e, n_edges());
108 848 : return std::find(std::begin(edge_nodes_map[e]),
109 12056 : std::end(edge_nodes_map[e]),
110 12056 : n) != std::end(edge_nodes_map[e]);
111 : }
112 :
113 :
114 :
115 41602 : bool Pyramid5::has_affine_map() const
116 : {
117 : // Point v = this->point(3) - this->point(0);
118 : // return (v.relative_fuzzy_equals(this->point(2) - this->point(1)));
119 41602 : return false;
120 : }
121 :
122 :
123 :
124 2316401 : Order Pyramid5::default_order() const
125 : {
126 2316401 : return FIRST;
127 : }
128 :
129 :
130 :
131 10060 : std::unique_ptr<Elem> Pyramid5::build_side_ptr (const unsigned int i)
132 : {
133 760 : libmesh_assert_less (i, this->n_sides());
134 :
135 10060 : std::unique_ptr<Elem> face;
136 :
137 10060 : switch (i)
138 : {
139 8048 : case 0: // triangular face 1
140 : case 1: // triangular face 2
141 : case 2: // triangular face 3
142 : case 3: // triangular face 4
143 : {
144 8048 : face = std::make_unique<Tri3>();
145 8048 : break;
146 : }
147 2012 : case 4: // the quad face at z=0
148 : {
149 2012 : face = std::make_unique<Quad4>();
150 2012 : break;
151 : }
152 0 : default:
153 0 : libmesh_error_msg("Invalid side i = " << i);
154 : }
155 :
156 : // Set the nodes
157 42252 : for (auto n : face->node_index_range())
158 34624 : face->set_node(n, this->node_ptr(Pyramid5::side_nodes_map[i][n]));
159 :
160 10060 : face->set_interior_parent(this);
161 9300 : face->inherit_data_from(*this);
162 :
163 10060 : return face;
164 0 : }
165 :
166 :
167 :
168 8995 : void Pyramid5::build_side_ptr (std::unique_ptr<Elem> & side,
169 : const unsigned int i)
170 : {
171 8995 : this->side_ptr(side, i);
172 8995 : side->set_interior_parent(this);
173 8265 : side->inherit_data_from(*this);
174 8995 : }
175 :
176 :
177 :
178 1704 : std::unique_ptr<Elem> Pyramid5::build_edge_ptr (const unsigned int i)
179 : {
180 1704 : return this->simple_build_edge_ptr<Edge2,Pyramid5>(i);
181 : }
182 :
183 :
184 :
185 0 : void Pyramid5::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
186 : {
187 0 : this->simple_build_edge_ptr<Pyramid5>(edge, i, EDGE2);
188 0 : }
189 :
190 :
191 :
192 0 : void Pyramid5::connectivity(const unsigned int libmesh_dbg_var(sc),
193 : const IOPackage iop,
194 : std::vector<dof_id_type> & conn) const
195 : {
196 0 : libmesh_assert(_nodes);
197 0 : libmesh_assert_less (sc, this->n_sub_elem());
198 0 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
199 :
200 0 : switch (iop)
201 : {
202 0 : case TECPLOT:
203 : {
204 0 : conn.resize(8);
205 0 : conn[0] = this->node_id(0)+1;
206 0 : conn[1] = this->node_id(1)+1;
207 0 : conn[2] = this->node_id(2)+1;
208 0 : conn[3] = this->node_id(3)+1;
209 0 : conn[4] = this->node_id(4)+1;
210 0 : conn[5] = this->node_id(4)+1;
211 0 : conn[6] = this->node_id(4)+1;
212 0 : conn[7] = this->node_id(4)+1;
213 0 : return;
214 : }
215 :
216 0 : case VTK:
217 : {
218 0 : conn.resize(5);
219 0 : conn[0] = this->node_id(3);
220 0 : conn[1] = this->node_id(2);
221 0 : conn[2] = this->node_id(1);
222 0 : conn[3] = this->node_id(0);
223 0 : conn[4] = this->node_id(4);
224 0 : return;
225 : }
226 :
227 0 : default:
228 0 : libmesh_error_msg("Unsupported IO package " << iop);
229 : }
230 : }
231 :
232 :
233 6359 : Point Pyramid5::true_centroid () const
234 : {
235 : // Call Hex8 static helper function, passing 4 copies of the final
236 : // vertex point, effectively treating the Pyramid as a degenerate
237 : // hexahedron. In my testing, this still seems to give correct
238 : // results.
239 5683 : return Hex8::centroid_from_points(
240 : point(0), point(1), point(2), point(3),
241 6697 : point(4), point(4), point(4), point(4));
242 : }
243 :
244 :
245 :
246 3408 : Real Pyramid5::volume () const
247 : {
248 : // The pyramid with a bilinear base has volume given by the
249 : // formula in: "Calculation of the Volume of a General Hexahedron
250 : // for Flow Predictions", AIAA Journal v.23, no.6, 1984, p.954-
251 : Point
252 3696 : x0 = point(0), x1 = point(1), x2 = point(2),
253 3504 : x3 = point(3), x4 = point(4);
254 :
255 : // Construct various edge and diagonal vectors.
256 96 : Point v40 = x0 - x4;
257 96 : Point v13 = x3 - x1;
258 96 : Point v02 = x2 - x0;
259 96 : Point v03 = x3 - x0;
260 96 : Point v01 = x1 - x0;
261 :
262 : // Finally, ready to return the volume!
263 : return
264 3408 : triple_product(v40, v13, v02) / 6. +
265 3408 : triple_product(v02, v01, v03) / 12.;
266 : }
267 :
268 : BoundingBox
269 133080 : Pyramid5::loose_bounding_box () const
270 : {
271 133080 : return Elem::loose_bounding_box();
272 : }
273 :
274 4788 : void Pyramid5::permute(unsigned int perm_num)
275 : {
276 300 : libmesh_assert_less (perm_num, 4);
277 :
278 11556 : for (unsigned int i = 0; i != perm_num; ++i)
279 : {
280 6768 : swap4nodes(0,1,2,3);
281 6336 : swap4neighbors(0,1,2,3);
282 : }
283 4788 : }
284 :
285 :
286 1728 : void Pyramid5::flip(BoundaryInfo * boundary_info)
287 : {
288 144 : libmesh_assert(boundary_info);
289 :
290 1728 : swap2nodes(0,1);
291 1728 : swap2nodes(2,3);
292 144 : swap2neighbors(1,3);
293 1728 : swap2boundarysides(1,3,boundary_info);
294 1728 : swap2boundaryedges(1,3,boundary_info);
295 1728 : swap2boundaryedges(4,5,boundary_info);
296 1728 : swap2boundaryedges(6,7,boundary_info);
297 1728 : }
298 :
299 :
300 8640 : ElemType Pyramid5::side_type (const unsigned int s) const
301 : {
302 720 : libmesh_assert_less (s, 5);
303 8640 : if (s < 4)
304 6912 : return TRI3;
305 144 : return QUAD4;
306 : }
307 :
308 :
309 : } // namespace libMesh
|
