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 7868 : 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 54809970 : 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 287899 : 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 : * \returns A simple Polygon face for side i, but from a clone so it's const
179 : */
180 : virtual std::unique_ptr<Elem> side_ptr (const unsigned int i) const;
181 :
182 : /**
183 : * Rebuilds a simple Polygon coincident with face i.
184 : */
185 : virtual void side_ptr (std::unique_ptr<Elem> & elem,
186 : const unsigned int i) override final;
187 :
188 : /**
189 : * Copies the Polygon side coincident with side i.
190 : */
191 : virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i) override;
192 :
193 : /**
194 : * Copies the Polygon side coincident with side i.
195 : */
196 : virtual void build_side_ptr (std::unique_ptr<Elem> & elem,
197 : const unsigned int i) override;
198 :
199 : // Avoid hiding deprecated version with different signature
200 : using Elem::build_side_ptr;
201 :
202 : /**
203 : * \returns An element coincident with edge \p i wrapped in a smart pointer.
204 : */
205 : virtual std::unique_ptr<Elem> build_edge_ptr (const unsigned int i) override final;
206 :
207 : /**
208 : * Resets the loose element \p edge, which may currently point to a
209 : * different edge than \p i or even a different element than \p
210 : * this, to point to edge \p i on \p this.
211 : */
212 : virtual void build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i) override final;
213 :
214 : /**
215 : * \returns The suggested quality bounds for
216 : * the polyhedron based on quality measure q.
217 : */
218 : virtual std::pair<Real, Real> qual_bounds (const ElemQuality q) const override;
219 :
220 : /**
221 : * \returns the (local) side numbers that touch the specified edge.
222 : */
223 : virtual std::vector<unsigned int> sides_on_edge(const unsigned int e) const override final;
224 :
225 : /**
226 : * \returns \p true if the specified edge is on the specified side.
227 : */
228 : virtual bool is_edge_on_side(const unsigned int e,
229 : const unsigned int s) const override final;
230 :
231 : /**
232 : * Maybe we have non-identity permutations, but trying to figure out
233 : * how many is an exercise in applied group theory, which is a bit
234 : * much for just expanded unit test coverage.
235 : */
236 552 : virtual unsigned int n_permutations() const override { return 1; }
237 :
238 552 : virtual void permute(unsigned int libmesh_dbg_var(perm_num)) override final
239 552 : { libmesh_assert_equal_to(perm_num, 0); }
240 :
241 : virtual bool is_flipped() const override final;
242 :
243 : /**
244 : * A flip is one of those general non-identity permutations we can't
245 : * handle.
246 : *
247 : * But we'll call it "not implemented" just in cases someone someday
248 : * wants to write an implementation to handle polyhedra with
249 : * topological symmetry planes.
250 : */
251 0 : virtual void flip(BoundaryInfo *) override final { libmesh_not_implemented(); };
252 :
253 : virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
254 :
255 : /**
256 : * Create a triangulation (tetrahedralization) from the current node
257 : * locations and face triangulations.
258 : *
259 : * If this is not called by the user, it will be called
260 : * automatically when needed.
261 : *
262 : * If the user moves the polyhedron nodes, the triangulation may
263 : * need to be regenerated manually, or the changed mapping may be
264 : * lower quality or even inverted.
265 : *
266 : * Pure virtual because this strongly depends on what mid-edge or
267 : * mid-face nodes the subclass might have.
268 : */
269 : virtual void retriangulate() = 0;
270 :
271 : /**
272 : * \returns true iff the polyhedron is convex. Some Polyhedron
273 : * methods currently assume (or in debugging modes, test for and
274 : * throw errors if not finding) convexity.
275 : */
276 : bool convex();
277 :
278 : /**
279 : * \returns the size of the triangulation of the polyhedron
280 : */
281 8125074 : unsigned int n_subelements() const { return cast_int<unsigned int>(this->_triangulation.size()); }
282 :
283 : /**
284 : * \returns the local indices of points on a subelement of the
285 : * polyhedron
286 : *
287 : * Each subelement here is a tetrahedron
288 : */
289 31795656 : virtual std::array<int, 4> subelement (unsigned int i) const
290 : {
291 2743208 : libmesh_assert_less(i, this->_triangulation.size());
292 34434224 : return this->_triangulation[i];
293 : }
294 :
295 : /**
296 : * \returns the master-space points of a subelement of the
297 : * polyhedron
298 : */
299 : virtual std::array<Point, 4> master_subelement (unsigned int i) const;
300 :
301 : /**
302 : * \returns the index of a subelement containing the master-space
303 : * point \p p, along with (the first three) barycentric coordinates
304 : * for \p; or return invalid_uint if no subelement contains p to
305 : * within tolerance \p tol.
306 : */
307 : std::tuple<unsigned int, Real, Real, Real>
308 : subelement_coordinates (const Point & p,
309 : Real tol = TOLERANCE*TOLERANCE) const;
310 :
311 : virtual bool on_reference_element(const Point & p,
312 : const Real eps = TOLERANCE) const override final;
313 :
314 : protected:
315 :
316 : /**
317 : * \returns Clones of the sides of \p this, wrapped in smart
318 : * pointers.
319 : *
320 : * This factors out much of the code required by the
321 : * disconnected_clone() method of Polyhedron subclasses.
322 : */
323 : std::vector<std::shared_ptr<Polygon>> side_clones() const;
324 :
325 : /**
326 : * Helper method for finding the non-cached side that shares an
327 : * edge, by examining the local node ids there.
328 : */
329 : bool side_has_edge_nodes(unsigned int side,
330 : unsigned int min_node,
331 : unsigned int max_node) const;
332 :
333 : /**
334 : * Data for links to parent/neighbor/interior_parent elements.
335 : *
336 : * There should be num_sides neighbors, preceded by any parent link
337 : * and followed by any interior_parent link.
338 : *
339 : * We're stuck with a heap allocation here since the number of sides
340 : * in a subclass is chosen at runtime.
341 : */
342 : std::vector<Elem *> _elemlinks_data;
343 :
344 : /**
345 : * Data for links to nodes. We're stuck with a heap allocation here
346 : * since the number of nodes in a subclass is chosen at runtime.
347 : */
348 : std::vector<Node *> _nodelinks_data;
349 :
350 : /**
351 : * Data for links to sides. We use shared_ptr so we can keep sides
352 : * consistent between neighbors. We also have a bool for each
353 : * polygon to indicate whether the zeta (xi cross eta) direction for
354 : * the polygon points into this polyhedron (true) or out of it
355 : * (false), and a map from side-local node number to element-local
356 : * node number.
357 : */
358 : std::vector<std::tuple<std::shared_ptr<Polygon>,
359 : bool,
360 : std::vector<unsigned int>>> _sidelinks_data;
361 :
362 : /**
363 : * One entry for each polyhedron edge, a pair indicating the side
364 : * number and the edge-of-side number which identifies the edge.
365 : *
366 : * Although each edge can be reached from two sides, only the
367 : * lower-numbered side is in this lookup table. For edges accessed
368 : * via the same side, the side edge numbering matches our edge
369 : * numbering.
370 : */
371 : std::vector<std::pair<unsigned int, unsigned int>> _edge_lookup;
372 :
373 : /**
374 : * Data for a triangulation (tetrahedralization) of the polyhedron.
375 : *
376 : * Positive int values should correspond to local node numbers.
377 : *
378 : * Negative int values may correspond to "special" points, e.g. a
379 : * centroid or skeleton point.
380 : */
381 : std::vector<std::array<int, 4>> _triangulation;
382 : };
383 :
384 :
385 : } // namespace libMesh
386 :
387 : #endif // LIBMESH_CELL_POLYHEDRON_H
|
