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 :
20 : #ifndef LIBMESH_CELL_POLYHEDRON_H
21 : #define LIBMESH_CELL_POLYHEDRON_H
22 :
23 :
24 : // Local includes
25 : #include "libmesh/libmesh_common.h"
26 : #include "libmesh/cell.h"
27 :
28 : namespace libMesh
29 : {
30 :
31 : // Forward declarations
32 : class Polygon;
33 :
34 : /**
35 : * The \p Polyhedron is an element in 3D with an arbitrary
36 : * number of polygonal faces.
37 : *
38 : * \author Roy H. Stogner
39 : * \date 2025
40 : * \brief A 3D element with an arbitrary number of polygonal faces
41 : */
42 : class Polyhedron : public Cell
43 : {
44 : public:
45 :
46 : /**
47 : * Arbitrary polyhedral element, takes a vector of shared pointers
48 : * to sides (which should already be constructed), and a (probably
49 : * unused) parent link. Derived classes implement 'true' elements.
50 : */
51 : Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
52 : Elem * p);
53 :
54 : Polyhedron (Polyhedron &&) = delete;
55 : Polyhedron (const Polyhedron &) = delete;
56 : Polyhedron & operator= (const Polyhedron &) = delete;
57 : Polyhedron & operator= (Polyhedron &&) = delete;
58 7588 : virtual ~Polyhedron() = default;
59 :
60 : /**
61 : * \returns true - polyhedron subclasses can have numbers of sides
62 : * and nodes which vary at runtime.
63 : */
64 12 : virtual bool runtime_topology() const override { return true; }
65 :
66 : /**
67 : * \returns The \p Point associated with local \p Node \p i,
68 : * in master element rather than physical coordinates.
69 : *
70 : * This implementation returns the master vertices; subclasses with
71 : * more points will need to further override.
72 : *
73 : * Trying to come up with a "master" shape for an arbitrary
74 : * polyhedron is *too* arbitrary, so we're just returning physical
75 : * points here.
76 : */
77 : virtual Point master_point (const unsigned int i) const override;
78 :
79 : static const int num_children = 0; // Refinement not yet supported
80 :
81 : /**
82 : * \returns the number of nodes in the polyhedron.
83 : */
84 54809770 : virtual unsigned int n_nodes() const override { return this->_nodelinks_data.size(); }
85 :
86 : /**
87 : * \returns the number of sides, which is the number of links we had
88 : * to save minus one parent link and minus one interior_parent link.
89 : */
90 287701 : virtual unsigned int n_sides() const override final { return _elemlinks_data.size()-2; }
91 :
92 : /**
93 : * \returns the number of edges.
94 : */
95 8651 : virtual unsigned int n_edges() const override final { return _edge_lookup.size(); }
96 :
97 : /**
98 : * \returns the number of faces, which is the number of links we had
99 : * to save minus one parent link and minus one interior_parent link.
100 : */
101 73 : virtual unsigned int n_faces() const override final { return _elemlinks_data.size()-2; }
102 :
103 : /**
104 : * \returns num_children.
105 : */
106 0 : virtual unsigned int n_children() const override final { return num_children; }
107 :
108 : /**
109 : * \returns \p true if the specified child is on the
110 : * specified side.
111 : *
112 : * Not implemented ... indefinitely? I don't think we'll be doing
113 : * hierarchic refinement for general polyhedra.
114 : */
115 : virtual bool is_child_on_side(const unsigned int c,
116 : const unsigned int s) const override;
117 :
118 : /**
119 : * Throws an error. I'm not sure how to define the "opposite
120 : * side" of a polyhedron.
121 : */
122 : virtual unsigned int opposite_side(const unsigned int s) const override final;
123 :
124 : /**
125 : * Throws an error - \p opposite_side(s) is too hard to define in
126 : * general on polyhedra.
127 : */
128 : virtual unsigned int opposite_node(const unsigned int n,
129 : const unsigned int s) const override final;
130 :
131 : /**
132 : * Don't hide Elem::key() defined in the base class.
133 : */
134 : using Elem::key;
135 :
136 : /**
137 : * \returns An id associated with the global node ids of this
138 : * element. The id is not necessarily unique, but should be
139 : * close.
140 : */
141 : virtual dof_id_type key () const override;
142 :
143 : /**
144 : * \returns An id associated with the \p s side of this element.
145 : * The id is not necessarily unique, but should be close.
146 : */
147 : virtual dof_id_type key (const unsigned int s) const override;
148 :
149 : /**
150 : * \returns An id associated with the \p s side of this element, as
151 : * defined solely by element vertices. The id is not necessarily
152 : * unique, but should be close. This is particularly useful in the
153 : * \p MeshBase::find_neighbors() routine.
154 : */
155 : virtual dof_id_type low_order_key (const unsigned int s) const override;
156 :
157 : /**
158 : * \returns The local node id for node \p side_node on side \p side of
159 : * this Elem.
160 : */
161 : virtual unsigned int local_side_node(unsigned int side,
162 : unsigned int side_node) const override;
163 :
164 : /**
165 : * Similar to Elem::local_side_node(), but instead of a side id, takes
166 : * an edge id and a node id on that edge and returns a local node number
167 : * for the Elem.
168 : */
169 : virtual unsigned int local_edge_node(unsigned int edge,
170 : unsigned int edge_node) const override;
171 :
172 : /**
173 : * \returns A simple Polygon face for side i.
174 : */
175 : virtual std::unique_ptr<Elem> side_ptr (const unsigned int i) override final;
176 :
177 : /**
178 : * Rebuilds a simple Polygon coincident with face i.
179 : */
180 : virtual void side_ptr (std::unique_ptr<Elem> & elem,
181 : const unsigned int i) override final;
182 :
183 : /**
184 : * Copies the Polygon side coincident with side i.
185 : */
186 : virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i) override;
187 :
188 : /**
189 : * Copies the Polygon side coincident with side i.
190 : */
191 : virtual void build_side_ptr (std::unique_ptr<Elem> & elem,
192 : const unsigned int i) override;
193 :
194 : // Avoid hiding deprecated version with different signature
195 : using Elem::build_side_ptr;
196 :
197 : /**
198 : * \returns An element coincident with edge \p i wrapped in a smart pointer.
199 : */
200 : virtual std::unique_ptr<Elem> build_edge_ptr (const unsigned int i) override final;
201 :
202 : /**
203 : * Resets the loose element \p edge, which may currently point to a
204 : * different edge than \p i or even a different element than \p
205 : * this, to point to edge \p i on \p this.
206 : */
207 : virtual void build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i) override final;
208 :
209 : /**
210 : * \returns The suggested quality bounds for
211 : * the polyhedron based on quality measure q.
212 : */
213 : virtual std::pair<Real, Real> qual_bounds (const ElemQuality q) const override;
214 :
215 : /**
216 : * \returns the (local) side numbers that touch the specified edge.
217 : */
218 : virtual std::vector<unsigned int> sides_on_edge(const unsigned int e) const override final;
219 :
220 : /**
221 : * \returns \p true if the specified edge is on the specified side.
222 : */
223 : virtual bool is_edge_on_side(const unsigned int e,
224 : const unsigned int s) const override final;
225 :
226 : /**
227 : * Maybe we have non-identity permutations, but trying to figure out
228 : * how many is an exercise in applied group theory, which is a bit
229 : * much for just expanded unit test coverage.
230 : */
231 552 : virtual unsigned int n_permutations() const override { return 1; }
232 :
233 552 : virtual void permute(unsigned int libmesh_dbg_var(perm_num)) override final
234 552 : { libmesh_assert_equal_to(perm_num, 0); }
235 :
236 : virtual bool is_flipped() const override final;
237 :
238 : /**
239 : * A flip is one of those general non-identity permutations we can't
240 : * handle.
241 : *
242 : * But we'll call it "not implemented" just in cases someone someday
243 : * wants to write an implementation to handle polyhedra with
244 : * topological symmetry planes.
245 : */
246 0 : virtual void flip(BoundaryInfo *) override final { libmesh_not_implemented(); };
247 :
248 : virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
249 :
250 : /**
251 : * Create a triangulation (tetrahedralization) from the current node
252 : * locations and face triangulations.
253 : *
254 : * If this is not called by the user, it will be called
255 : * automatically when needed.
256 : *
257 : * If the user moves the polyhedron nodes, the triangulation may
258 : * need to be regenerated manually, or the changed mapping may be
259 : * lower quality or even inverted.
260 : *
261 : * Pure virtual because this strongly depends on what mid-edge or
262 : * mid-face nodes the subclass might have.
263 : */
264 : virtual void retriangulate() = 0;
265 :
266 : /**
267 : * \returns true iff the polyhedron is convex. Some Polyhedron
268 : * methods currently assume (or in debugging modes, test for and
269 : * throw errors if not finding) convexity.
270 : */
271 : bool convex();
272 :
273 : /**
274 : * \returns the size of the triangulation of the polyhedron
275 : */
276 8125074 : unsigned int n_subelements() const { return cast_int<unsigned int>(this->_triangulation.size()); }
277 :
278 : /**
279 : * \returns the local indices of points on a subelement of the
280 : * polyhedron
281 : *
282 : * Each subelement here is a tetrahedron
283 : */
284 31795656 : virtual std::array<int, 4> subelement (unsigned int i) const
285 : {
286 2743208 : libmesh_assert_less(i, this->_triangulation.size());
287 34434224 : return this->_triangulation[i];
288 : }
289 :
290 : /**
291 : * \returns the master-space points of a subelement of the
292 : * polyhedron
293 : */
294 : virtual std::array<Point, 4> master_subelement (unsigned int i) const;
295 :
296 : /**
297 : * \returns the index of a subelement containing the master-space
298 : * point \p p, along with (the first three) barycentric coordinates
299 : * for \p; or return invalid_uint if no subelement contains p to
300 : * within tolerance \p tol.
301 : */
302 : std::tuple<unsigned int, Real, Real, Real>
303 : subelement_coordinates (const Point & p,
304 : Real tol = TOLERANCE*TOLERANCE) const;
305 :
306 : virtual bool on_reference_element(const Point & p,
307 : const Real eps = TOLERANCE) const override final;
308 :
309 : protected:
310 :
311 : /**
312 : * \returns Clones of the sides of \p this, wrapped in smart
313 : * pointers.
314 : *
315 : * This factors out much of the code required by the
316 : * disconnected_clone() method of Polyhedron subclasses.
317 : */
318 : std::vector<std::shared_ptr<Polygon>> side_clones() const;
319 :
320 : /**
321 : * Helper method for finding the non-cached side that shares an
322 : * edge, by examining the local node ids there.
323 : */
324 : bool side_has_edge_nodes(unsigned int side,
325 : unsigned int min_node,
326 : unsigned int max_node) const;
327 :
328 : /**
329 : * Data for links to parent/neighbor/interior_parent elements.
330 : *
331 : * There should be num_sides neighbors, preceded by any parent link
332 : * and followed by any interior_parent link.
333 : *
334 : * We're stuck with a heap allocation here since the number of sides
335 : * in a subclass is chosen at runtime.
336 : */
337 : std::vector<Elem *> _elemlinks_data;
338 :
339 : /**
340 : * Data for links to nodes. We're stuck with a heap allocation here
341 : * since the number of nodes in a subclass is chosen at runtime.
342 : */
343 : std::vector<Node *> _nodelinks_data;
344 :
345 : /**
346 : * Data for links to sides. We use shared_ptr so we can keep sides
347 : * consistent between neighbors. We also have a bool for each
348 : * polygon to indicate whether the zeta (xi cross eta) direction for
349 : * the polygon points into this polyhedron (true) or out of it
350 : * (false), and a map from side-local node number to element-local
351 : * node number.
352 : */
353 : std::vector<std::tuple<std::shared_ptr<Polygon>,
354 : bool,
355 : std::vector<unsigned int>>> _sidelinks_data;
356 :
357 : /**
358 : * One entry for each polyhedron edge, a pair indicating the side
359 : * number and the edge-of-side number which identifies the edge.
360 : *
361 : * Although each edge can be reached from two sides, only the
362 : * lower-numbered side is in this lookup table. For edges accessed
363 : * via the same side, the side edge numbering matches our edge
364 : * numbering.
365 : */
366 : std::vector<std::pair<unsigned int, unsigned int>> _edge_lookup;
367 :
368 : /**
369 : * Data for a triangulation (tetrahedralization) of the polyhedron.
370 : *
371 : * Positive int values should correspond to local node numbers.
372 : *
373 : * Negative int values may correspond to "special" points, e.g. a
374 : * centroid or skeleton point.
375 : */
376 : std::vector<std::array<int, 4>> _triangulation;
377 : };
378 :
379 :
380 : } // namespace libMesh
381 :
382 : #endif // LIBMESH_CELL_POLYHEDRON_H
|
