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 1879 : 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 : * When saving the midpoint location data, we need to save the
291 : * corresponding segment information too. Here the first point
292 : * of the segment is saved so that it can be used as a key to
293 : * find the corresponding segment midpoint.
294 : */
295 : std::vector<Point> segment_midpoints_keys;
296 :
297 : /**
298 : * Attaches boundary markers.
299 : * If segments is set, the number of markers must be equal to the size of segments,
300 : * otherwise, it is equal to the number of points.
301 : */
302 : void attach_boundary_marker(const std::vector<int> * markers) { _markers = markers; }
303 :
304 : /**
305 : * Attaches regions for using attribute to set subdomain IDs and better
306 : * controlling the triangle sizes within the regions.
307 : */
308 : void attach_region_list(const std::vector<Region*> * regions) { _regions = regions; }
309 :
310 : /**
311 : * Generate an auto area function based on spacing of boundary points.
312 : * The external boundary as well as the hole boundaries are taken into consideration
313 : * to generate the auto area function based on inverse distance interpolation.
314 : * For each EDGE element on these boundaries, its centroid (midpoint) is used as the point
315 : * position and the desired area is calculated as 1.5 times of the area of the equilateral
316 : * triangle with the edge length as the length of the EDGE element.
317 : * For a given position, the inverse distance interpolation only considers a number of nearest
318 : * points (set by num_nearest_pts) to calculate the desired area. The weight of the value at
319 : * each point is calculated as 1/distance^power.
320 : *
321 : * In addition to these conventional inverse distance interpolation features, a concept of
322 : * "background value" and "background effective distance" is introduced. The background value
323 : * belongs to a virtual point located at a constant distance (background effective distance)
324 : * from the given position. The weight of the value at this virtual point is calculated as
325 : * 1/background_effective_distance^power. Effectively, the background value is the value when
326 : * the given position is far away from the boundary points.
327 : */
328 : void set_auto_area_function(const Parallel::Communicator &comm,
329 : const unsigned int num_nearest_pts,
330 : const unsigned int power,
331 : const Real background_value,
332 : const Real background_eff_dist);
333 :
334 : /**
335 : * Whether or not an auto area function has been set
336 : */
337 42926 : bool has_auto_area_function() {return _auto_area_function != nullptr;}
338 :
339 : /**
340 : * Get the auto area function
341 : */
342 : FunctionBase<Real> * get_auto_area_function ();
343 :
344 : /**
345 : * The external boundary and all hole boundaries are collected. The centroid of each EDGE element is
346 : * used as the point position and the desired area is calculated as the area of the equilateral triangle
347 : * with the edge length as the length of the EDGE element times an area_factor (default is 1.5).
348 : */
349 : void calculate_auto_desired_area_samples(std::vector<Point> & function_points,
350 : std::vector<Real> & function_sizes,
351 0 : const Real & area_factor = 1.5);
352 :
353 : /**
354 : * A set of ids to allow on the outer boundary loop: interpreted as
355 : * boundary ids of 2D elements and/or subdomain ids of 1D edges. If
356 : * this is empty, then the outer boundary may be constructed from
357 : * boundary edges of any id!
358 : */
359 : void set_outer_boundary_ids(std::set<std::size_t> bdy_ids) { _bdy_ids = std::move(bdy_ids); }
360 : const std::set<std::size_t> & get_outer_boundary_ids() const { return _bdy_ids; }
361 :
362 : protected:
363 : /**
364 : * Helper function to create PSLG segments from our other
365 : * boundary-defining options (1D mesh edges, 2D mesh
366 : * boundary sides), if no segments already exist.
367 : */
368 : void elems_to_segments();
369 :
370 : /**
371 : * Helper function to create PSLG segments from our node ordering,
372 : * up to the maximum node id, if no segments already exist.
373 : */
374 : void nodes_to_segments(dof_id_type max_node_id);
375 :
376 : /**
377 : * Helper function to add extra points (midpoints of initial
378 : * segments) to a PSLG triangulation
379 : */
380 : void insert_any_extra_boundary_points();
381 :
382 : /**
383 : * Helper function to upconvert Tri3 to any higher order triangle
384 : * type if requested via \p _elem_type. Should be called at the end
385 : * of triangulate()
386 : */
387 : void increase_triangle_order();
388 :
389 : /**
390 : * Helper function to check holes for intersections if requested.
391 : */
392 : void verify_holes(const Hole & outer_bdy);
393 :
394 : /**
395 : * Helper function to count points in and verify holes
396 : */
397 : unsigned int total_hole_points();
398 :
399 : /**
400 : * Reference to the mesh which is to be created by triangle.
401 : */
402 : UnstructuredMesh & _mesh;
403 :
404 : /**
405 : * A pointer to a vector of Hole*s. If this is nullptr, there
406 : * are no holes!
407 : */
408 : const std::vector<Hole*> * _holes;
409 :
410 : /**
411 : * Boundary markers
412 : */
413 : const std::vector<int> * _markers;
414 :
415 : /**
416 : * A pointer to a vector of Regions*s. If this is nullptr, there
417 : * are no regions!
418 : */
419 : const std::vector<Region*> * _regions;
420 :
421 : /**
422 : * A set of ids to allow on the outer boundary loop.
423 : */
424 : std::set<std::size_t> _bdy_ids;
425 :
426 : /**
427 : * The type of elements to generate. (Defaults to
428 : * TRI3).
429 : */
430 : ElemType _elem_type;
431 :
432 : /**
433 : * The desired area for the elements in the resulting mesh.
434 : */
435 : Real _desired_area;
436 :
437 : /**
438 : * Minimum angle in triangles
439 : */
440 : Real _minimum_angle;
441 :
442 : /**
443 : * The type of triangulation to perform: choices are:
444 : * convex hull
445 : * PSLG
446 : */
447 : TriangulationType _triangulation_type;
448 :
449 : /**
450 : * Flag which tells whether or not to insert additional nodes
451 : * before triangulation. This can sometimes be used to "de-regularize"
452 : * the resulting triangulation.
453 : *
454 : * This flag is supported for backwards compatibility; setting
455 : * _interpolate_boundary_points = 1 is equivalent.
456 : */
457 : bool _insert_extra_points;
458 :
459 : /**
460 : * Flag which tells how many additional nodes should be inserted
461 : * between each pair of original mesh points.
462 : */
463 : int _interpolate_boundary_points;
464 :
465 : /**
466 : * Flag which tells whether we should smooth the mesh after
467 : * it is generated. True by default.
468 : */
469 : bool _smooth_after_generating;
470 :
471 : /**
472 : * Flag which tells if we want to suppress stdout outputs
473 : */
474 : bool _quiet;
475 :
476 : /**
477 : * Flag which tells if we want to check hole geometry
478 : */
479 : bool _verify_hole_boundaries;
480 :
481 : /**
482 : * The auto area function based on the spacing of boundary points
483 : */
484 : std::unique_ptr<AutoAreaFunction> _auto_area_function;
485 : };
486 :
487 : } // namespace libMesh
488 :
489 : #endif // ifndef LIBMESH_MESH_TRIANGULATOR_INTERFACE_H
|
