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 : #ifndef LIBMESH_MESH_TRIANGULATOR_INTERFACE_H
20 : #define LIBMESH_MESH_TRIANGULATOR_INTERFACE_H
21 :
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : // Local Includes
26 : #include "libmesh/libmesh.h"
27 : #include "libmesh/point.h"
28 :
29 : #include "libmesh/function_base.h"
30 :
31 : // C++ includes
32 : #include <set>
33 : #include <vector>
34 :
35 : namespace libMesh
36 : {
37 :
38 : // Forward Declarations
39 : class UnstructuredMesh;
40 : template <unsigned int KDDim>
41 : class InverseDistanceInterpolation;
42 : enum ElemType : int;
43 :
44 : // Helper function for auto area calculation
45 0 : class AutoAreaFunction : public FunctionBase<Real>
46 : {
47 : public:
48 : AutoAreaFunction (const Parallel::Communicator &comm,
49 : const unsigned int num_nearest_pts,
50 : const unsigned int power,
51 : const Real background_value,
52 : const Real background_eff_dist);
53 :
54 : AutoAreaFunction (const AutoAreaFunction &);
55 : AutoAreaFunction & operator= (const AutoAreaFunction &);
56 :
57 : AutoAreaFunction (AutoAreaFunction &&) = delete;
58 : AutoAreaFunction & operator= (AutoAreaFunction &&) = delete;
59 : virtual ~AutoAreaFunction ();
60 :
61 : void init_mfi (const std::vector<Point> & input_pts,
62 : const std::vector<Real> & input_vals);
63 :
64 : virtual Real operator() (const Point & p,
65 : const Real /*time*/) override;
66 :
67 0 : virtual void operator() (const Point & p,
68 : const Real time,
69 : DenseVector<Real> & output) override
70 : {
71 0 : output.resize(1);
72 0 : output(0) = (*this)(p,time);
73 0 : return;
74 : }
75 :
76 0 : virtual std::unique_ptr<FunctionBase<Real>> clone() const override
77 : {
78 0 : return std::make_unique<AutoAreaFunction>(_comm, _num_nearest_pts, _power, _background_value, _background_eff_dist);
79 : }
80 :
81 : private:
82 : const Parallel::Communicator & _comm;
83 : const unsigned int _num_nearest_pts;
84 : const unsigned int _power;
85 : const Real _background_value;
86 : const Real _background_eff_dist;
87 : std::unique_ptr<InverseDistanceInterpolation<3>> _auto_area_mfi;
88 : };
89 :
90 : class TriangulatorInterface
91 : {
92 : public:
93 : /**
94 : * The constructor. A reference to the mesh containing the points
95 : * which are to be triangulated must be provided. Unless otherwise
96 : * specified, a convex hull will be computed for the set of input points
97 : * and the convex hull will be meshed.
98 : */
99 : explicit
100 : TriangulatorInterface(UnstructuredMesh & mesh);
101 :
102 : /**
103 : * Empty destructor.
104 : */
105 1810 : virtual ~TriangulatorInterface() = default;
106 :
107 : /**
108 : * The TriangulationType is used with the general triangulate function
109 : * defined below.
110 : */
111 : enum TriangulationType
112 : {
113 : /**
114 : * First generate a convex hull from the set of points passed
115 : * in, and then triangulate this set of points. This is
116 : * probably the most common type of usage.
117 : */
118 : GENERATE_CONVEX_HULL = 0,
119 :
120 : /**
121 : * Triangulate the interior of a Planar Straight Line Graph,
122 : * which is defined implicitly by the order of the "points"
123 : * vector: a straight line is assumed to lie between each
124 : * successive pair of points, with an additional line joining
125 : * the final and first points.
126 : *
127 : * Explicitly telling the triangulator to add additional points
128 : * may be important for this option.
129 : */
130 : PSLG = 1,
131 :
132 : /**
133 : * Does nothing, used as a "null" value.
134 : */
135 : INVALID_TRIANGULATION_TYPE
136 : };
137 :
138 : /**
139 : * The hole class and its several subclasses define the interface
140 : * and functionality of a "hole" which appears in a 2D mesh.
141 : * See mesh_triangle_holes.C/h for definitions.
142 : */
143 : class Hole;
144 : class AffineHole;
145 : class PolygonHole;
146 : class ArbitraryHole;
147 : class MeshedHole;
148 :
149 : /**
150 : * The region class defines the interface
151 : * and functionality of a "region" which appears in a 2D mesh.
152 : * See mesh_triangle_holes.C/h for definitions.
153 : */
154 : class Region;
155 :
156 : /**
157 : * This is the main public interface for this function.
158 : */
159 : virtual void triangulate() = 0;
160 :
161 : /**
162 : * Sets and/or gets the desired element type.
163 : */
164 0 : ElemType & elem_type() {return _elem_type;}
165 :
166 : /**
167 : * Sets and/or gets the desired triangle area. Set to zero to disable
168 : * area constraint.
169 : *
170 : * If a \p desired_area_function is set, then \p desired_area()
171 : * should be used to set a *minimum* desired area; this will reduce
172 : * "false negatives" by suggesting how finely to sample \p
173 : * desired_area_function inside large triangles, where ideally the
174 : * \p desired_area_function will be satisfied in the triangle
175 : * interior and not just at the triangle vertices.
176 : */
177 230 : Real & desired_area() {return _desired_area;}
178 :
179 : /**
180 : * Set a function giving desired triangle area as a function of
181 : * position. Set this to nullptr to disable position-dependent area
182 : * constraint (falling back on desired_area()).
183 : *
184 : * This may not be implemented in all subclasses.
185 : */
186 0 : virtual void set_desired_area_function (FunctionBase<Real> *)
187 0 : { libmesh_not_implemented(); }
188 :
189 : /**
190 : * Get the function giving desired triangle area as a function of
191 : * position, or \p nullptr if no such function has been set.
192 : */
193 0 : virtual FunctionBase<Real> * get_desired_area_function ()
194 0 : { return nullptr; }
195 :
196 : /**
197 : * Sets and/or gets the minimum desired angle. Set to zero to
198 : * disable angle constraint.
199 : */
200 230 : Real & minimum_angle() {return _minimum_angle;}
201 :
202 : /**
203 : * Sets and/or gets the desired triangulation type.
204 : */
205 0 : TriangulationType & triangulation_type() {return _triangulation_type;}
206 :
207 : /**
208 : * Sets and/or gets the flag for inserting add'l points.
209 : */
210 : bool & insert_extra_points() {return _insert_extra_points;}
211 :
212 : /**
213 : * Complicated setter, for compatibility with insert_extra_points()
214 : */
215 : void set_interpolate_boundary_points (int n_points);
216 :
217 : /**
218 : * Complicated getter, for compatibility with insert_extra_points()
219 : */
220 : int get_interpolate_boundary_points () const;
221 :
222 : /**
223 : * Set whether or not the triangulation is allowed to refine the
224 : * mesh boundary when refining the interior. This is true by
225 : * default, but may be set to false to make the mesh boundary more
226 : * predictable (and so easier to stitch to other meshes) later.
227 : *
228 : * This may not be implemented in all subclasses.
229 : */
230 0 : virtual void set_refine_boundary_allowed (bool)
231 0 : { libmesh_not_implemented(); }
232 :
233 : /**
234 : * Get whether or not the triangulation is allowed to refine the
235 : * mesh boundary when refining the interior. True by default.
236 : */
237 0 : virtual bool refine_boundary_allowed () const
238 0 : { return true; }
239 :
240 : /**
241 : * Sets/gets flag which tells whether to do two steps of Laplace
242 : * mesh smoothing after generating the grid.
243 : */
244 230 : bool & smooth_after_generating() {return _smooth_after_generating;}
245 :
246 : /**
247 : * Whether not to silence internal messages to stdout
248 : */
249 : bool & quiet() {return _quiet;}
250 :
251 : /**
252 : * Attaches a vector of Hole* pointers which will be
253 : * meshed around.
254 : */
255 0 : void attach_hole_list(const std::vector<Hole*> * holes) {_holes = holes;}
256 :
257 : /**
258 : * Verifying that hole boundaries don't cross the outer boundary or
259 : * each other is something like O(N_bdys^2*N_points_per_bdy^2), so
260 : * we only do it if requested
261 : */
262 : void set_verify_hole_boundaries(bool v) {_verify_hole_boundaries = v;}
263 :
264 : bool get_verify_hole_boundaries() const {return _verify_hole_boundaries;}
265 :
266 : /**
267 : * When constructing a PSLG, if the node numbers do not define the
268 : * desired boundary segments implicitly through the ordering of the
269 : * points, you can use the segments vector to specify the segments
270 : * explicitly, Ex: unit square numbered counter-clockwise starting
271 : * from origin
272 : * segments[0] = (0,1)
273 : * segments[1] = (1,2)
274 : * segments[2] = (2,3)
275 : * segments[3] = (3,0)
276 : * (For the above case you could actually use the implicit
277 : * ordering!)
278 : */
279 : std::vector<std::pair<unsigned int, unsigned int>> segments;
280 :
281 : /**
282 : * When constructing a second-order triangulation from a
283 : * second-order boundary, we may do the triangulation using
284 : * first-order elements, in which case we need to save midpoint
285 : * location data in order to reconstruct curvature along boundaries.
286 : */
287 : std::vector<Point> segment_midpoints;
288 :
289 : /**
290 : * Attaches boundary markers.
291 : * If segments is set, the number of markers must be equal to the size of segments,
292 : * otherwise, it is equal to the number of points.
293 : */
294 : void attach_boundary_marker(const std::vector<int> * markers) { _markers = markers; }
295 :
296 : /**
297 : * Attaches regions for using attribute to set subdomain IDs and better
298 : * controlling the triangle sizes within the regions.
299 : */
300 : void attach_region_list(const std::vector<Region*> * regions) { _regions = regions; }
301 :
302 : /**
303 : * Generate an auto area function based on spacing of boundary points.
304 : * The external boundary as well as the hole boundaries are taken into consideration
305 : * to generate the auto area function based on inverse distance interpolation.
306 : * For each EDGE element on these boundaries, its centroid (midpoint) is used as the point
307 : * position and the desired area is calculated as 1.5 times of the area of the equilateral
308 : * triangle with the edge length as the length of the EDGE element.
309 : * For a given position, the inverse distance interpolation only considers a number of nearest
310 : * points (set by num_nearest_pts) to calculate the desired area. The weight of the value at
311 : * each point is calculated as 1/distance^power.
312 : *
313 : * In addition to these conventional inverse distance interpolation features, a concept of
314 : * "background value" and "background effective distance" is introduced. The background value
315 : * belongs to a virtual point located at a constant distance (background effective distance)
316 : * from the given position. The weight of the value at this virtual point is calculated as
317 : * 1/background_effective_distance^power. Effectively, the background value is the value when
318 : * the given position is far away from the boundary points.
319 : */
320 : void set_auto_area_function(const Parallel::Communicator &comm,
321 : const unsigned int num_nearest_pts,
322 : const unsigned int power,
323 : const Real background_value,
324 : const Real background_eff_dist);
325 :
326 : /**
327 : * Whether or not an auto area function has been set
328 : */
329 42918 : bool has_auto_area_function() {return _auto_area_function != nullptr;}
330 :
331 : /**
332 : * Get the auto area function
333 : */
334 : FunctionBase<Real> * get_auto_area_function ();
335 :
336 : /**
337 : * The external boundary and all hole boundaries are collected. The centroid of each EDGE element is
338 : * used as the point position and the desired area is calculated as the area of the equilateral triangle
339 : * with the edge length as the length of the EDGE element times an area_factor (default is 1.5).
340 : */
341 : void calculate_auto_desired_area_samples(std::vector<Point> & function_points,
342 : std::vector<Real> & function_sizes,
343 0 : const Real & area_factor = 1.5);
344 :
345 : /**
346 : * A set of ids to allow on the outer boundary loop: interpreted as
347 : * boundary ids of 2D elements and/or subdomain ids of 1D edges. If
348 : * this is empty, then the outer boundary may be constructed from
349 : * boundary edges of any id!
350 : */
351 : void set_outer_boundary_ids(std::set<std::size_t> bdy_ids) { _bdy_ids = std::move(bdy_ids); }
352 : const std::set<std::size_t> & get_outer_boundary_ids() const { return _bdy_ids; }
353 :
354 : protected:
355 : /**
356 : * Helper function to create PSLG segments from our other
357 : * boundary-defining options (1D mesh edges, 2D mesh
358 : * boundary sides), if no segments already exist.
359 : */
360 : void elems_to_segments();
361 :
362 : /**
363 : * Helper function to create PSLG segments from our node ordering,
364 : * up to the maximum node id, if no segments already exist.
365 : */
366 : void nodes_to_segments(dof_id_type max_node_id);
367 :
368 : /**
369 : * Helper function to add extra points (midpoints of initial
370 : * segments) to a PSLG triangulation
371 : */
372 : void insert_any_extra_boundary_points();
373 :
374 : /**
375 : * Helper function to upconvert Tri3 to any higher order triangle
376 : * type if requested via \p _elem_type. Should be called at the end
377 : * of triangulate()
378 : */
379 : void increase_triangle_order();
380 :
381 : /**
382 : * Helper function to check holes for intersections if requested.
383 : */
384 : void verify_holes(const Hole & outer_bdy);
385 :
386 : /**
387 : * Helper function to count points in and verify holes
388 : */
389 : unsigned int total_hole_points();
390 :
391 : /**
392 : * Reference to the mesh which is to be created by triangle.
393 : */
394 : UnstructuredMesh & _mesh;
395 :
396 : /**
397 : * A pointer to a vector of Hole*s. If this is nullptr, there
398 : * are no holes!
399 : */
400 : const std::vector<Hole*> * _holes;
401 :
402 : /**
403 : * Boundary markers
404 : */
405 : const std::vector<int> * _markers;
406 :
407 : /**
408 : * A pointer to a vector of Regions*s. If this is nullptr, there
409 : * are no regions!
410 : */
411 : const std::vector<Region*> * _regions;
412 :
413 : /**
414 : * A set of ids to allow on the outer boundary loop.
415 : */
416 : std::set<std::size_t> _bdy_ids;
417 :
418 : /**
419 : * The type of elements to generate. (Defaults to
420 : * TRI3).
421 : */
422 : ElemType _elem_type;
423 :
424 : /**
425 : * The desired area for the elements in the resulting mesh.
426 : */
427 : Real _desired_area;
428 :
429 : /**
430 : * Minimum angle in triangles
431 : */
432 : Real _minimum_angle;
433 :
434 : /**
435 : * The type of triangulation to perform: choices are:
436 : * convex hull
437 : * PSLG
438 : */
439 : TriangulationType _triangulation_type;
440 :
441 : /**
442 : * Flag which tells whether or not to insert additional nodes
443 : * before triangulation. This can sometimes be used to "de-regularize"
444 : * the resulting triangulation.
445 : *
446 : * This flag is supported for backwards compatibility; setting
447 : * _interpolate_boundary_points = 1 is equivalent.
448 : */
449 : bool _insert_extra_points;
450 :
451 : /**
452 : * Flag which tells how many additional nodes should be inserted
453 : * between each pair of original mesh points.
454 : */
455 : int _interpolate_boundary_points;
456 :
457 : /**
458 : * Flag which tells whether we should smooth the mesh after
459 : * it is generated. True by default.
460 : */
461 : bool _smooth_after_generating;
462 :
463 : /**
464 : * Flag which tells if we want to suppress stdout outputs
465 : */
466 : bool _quiet;
467 :
468 : /**
469 : * Flag which tells if we want to check hole geometry
470 : */
471 : bool _verify_hole_boundaries;
472 :
473 : /**
474 : * The auto area function based on the spacing of boundary points
475 : */
476 : std::unique_ptr<AutoAreaFunction> _auto_area_function;
477 : };
478 :
479 : } // namespace libMesh
480 :
481 : #endif // ifndef LIBMESH_MESH_TRIANGULATOR_INTERFACE_H
|
