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 : #include "libmesh/libmesh_config.h"
20 :
21 : #ifdef LIBMESH_HAVE_TRIANGLE
22 :
23 : // libmesh includes
24 : #include "libmesh/mesh_triangle_interface.h"
25 : #include "libmesh/unstructured_mesh.h"
26 : #include "libmesh/face_tri3.h"
27 : #include "libmesh/face_tri6.h"
28 : #include "libmesh/mesh_generation.h"
29 : #include "libmesh/mesh_smoother_laplace.h"
30 : #include "libmesh/boundary_info.h"
31 : #include "libmesh/mesh_triangle_holes.h"
32 : #include "libmesh/mesh_triangle_wrapper.h"
33 : #include "libmesh/enum_elem_type.h"
34 : #include "libmesh/enum_to_string.h"
35 :
36 : // C/C++ includes
37 : #include <sstream>
38 :
39 :
40 : namespace libMesh
41 : {
42 :
43 : //
44 : // Function definitions for the TrianguleInterface class
45 : //
46 :
47 : // Constructor
48 56 : TriangleInterface::TriangleInterface(UnstructuredMesh & mesh)
49 : : TriangulatorInterface(mesh),
50 : _extra_flags(""),
51 56 : _serializer(_mesh)
52 56 : {}
53 :
54 :
55 :
56 :
57 : // Primary function responsible for performing the triangulation
58 500 : void TriangleInterface::triangulate()
59 : {
60 : // Will the triangulation have holes?
61 500 : const bool have_holes = ((_holes != nullptr) && (!_holes->empty()));
62 :
63 500 : unsigned int n_hole_points = this->total_hole_points();
64 :
65 : // If we have no explicit segments defined, we may get them from
66 : // mesh elements
67 500 : this->elems_to_segments();
68 :
69 : // If we're doing PSLG without segments, construct them from all our
70 : // mesh nodes
71 500 : this->nodes_to_segments(_mesh.max_node_id());
72 :
73 : // Insert additional new points in between existing boundary points,
74 : // if that is requested and reasonable
75 500 : this->insert_any_extra_boundary_points();
76 :
77 : // Regardless of whether we added additional points, the set of points to
78 : // triangulate is now sitting in the mesh.
79 :
80 : // Triangle data structure for the mesh
81 : TriangleWrapper::triangulateio initial;
82 : TriangleWrapper::triangulateio final;
83 : TriangleWrapper::triangulateio voronoi;
84 :
85 : // Pseudo-Constructor for the triangle io structs
86 500 : TriangleWrapper::init(initial);
87 500 : TriangleWrapper::init(final);
88 500 : TriangleWrapper::init(voronoi);
89 :
90 500 : initial.numberofpoints = _mesh.n_nodes() + n_hole_points;
91 500 : initial.pointlist = static_cast<REAL*>(std::malloc(initial.numberofpoints * 2 * sizeof(REAL)));
92 :
93 500 : if (_triangulation_type==PSLG)
94 : {
95 : // Implicit segment ordering: One segment per point, including hole points
96 40 : if (this->segments.empty())
97 0 : initial.numberofsegments = initial.numberofpoints;
98 :
99 : // User-defined segment ordering: One segment per entry in the segments vector
100 : else
101 40 : initial.numberofsegments = this->segments.size() + n_hole_points;
102 : }
103 :
104 460 : else if (_triangulation_type==GENERATE_CONVEX_HULL)
105 460 : initial.numberofsegments = n_hole_points; // One segment for each hole point
106 :
107 : // Allocate space for the segments (2 int per segment)
108 500 : if (initial.numberofsegments > 0)
109 : {
110 40 : initial.segmentlist = static_cast<int *> (std::malloc(initial.numberofsegments * 2 * sizeof(int)));
111 40 : if (_markers)
112 0 : initial.segmentmarkerlist = static_cast<int *> (std::malloc(initial.numberofsegments * sizeof(int)));
113 : }
114 :
115 :
116 : // Copy all the holes' points and segments into the triangle struct.
117 :
118 : // The hole_offset is a constant offset into the points vector which points
119 : // past the end of the last hole point added.
120 250 : unsigned int hole_offset=0;
121 :
122 500 : if (have_holes)
123 40 : for (const auto & hole : *_holes)
124 : {
125 60 : for (unsigned int ctr=0, h=0, i=0, hsism=hole->segment_indices().size()-1; i<hsism; ++i)
126 : {
127 24 : unsigned int begp = hole_offset + hole->segment_indices()[i];
128 36 : unsigned int endp = hole->segment_indices()[i+1];
129 :
130 716 : for (; h<endp; ctr+=2, ++h)
131 : {
132 692 : Point p = hole->point(h);
133 :
134 692 : const unsigned int index0 = 2*hole_offset+ctr;
135 692 : const unsigned int index1 = 2*hole_offset+ctr+1;
136 :
137 : // Save the x,y locations in the triangle struct.
138 692 : initial.pointlist[index0] = p(0);
139 692 : initial.pointlist[index1] = p(1);
140 :
141 : // Set the points which define the segments
142 692 : initial.segmentlist[index0] = hole_offset+h;
143 692 : initial.segmentlist[index1] = (h == endp - 1) ? begp : hole_offset + h + 1; // wrap around
144 692 : if (_markers)
145 : // 1 is reserved for boundaries of holes
146 0 : initial.segmentmarkerlist[hole_offset+h] = 1;
147 : }
148 : }
149 :
150 : // Update the hole_offset for the next hole
151 24 : hole_offset += hole->n_points();
152 : }
153 :
154 :
155 : // Copy all the non-hole points and segments into the triangle struct.
156 750 : std::vector<unsigned int> libmesh_id_to_pointlist_index(_mesh.max_node_id());
157 : {
158 250 : dof_id_type ctr=0;
159 2634 : for (auto & node : _mesh.node_ptr_range())
160 : {
161 1884 : dof_id_type index = 2*hole_offset + ctr;
162 :
163 : // Set x,y values in pointlist
164 1884 : initial.pointlist[index] = (*node)(0);
165 1884 : initial.pointlist[index+1] = (*node)(1);
166 1884 : libmesh_id_to_pointlist_index[node->id()] = ctr/2;
167 :
168 : // If the user requested a PSLG, the non-hole points are also segments
169 1884 : if (_triangulation_type==PSLG)
170 : {
171 : // Use implicit ordering to define segments
172 216 : if (this->segments.empty())
173 : {
174 0 : dof_id_type n = ctr/2; // ctr is always even
175 0 : initial.segmentlist[index] = hole_offset+n;
176 0 : initial.segmentlist[index+1] = (n==_mesh.n_nodes()-1) ? hole_offset : hole_offset+n+1; // wrap around
177 0 : if (_markers)
178 0 : initial.segmentmarkerlist[hole_offset + n] = (*_markers)[n];
179 : }
180 : }
181 :
182 1884 : ctr +=2;
183 : }
184 : }
185 :
186 :
187 : // If the user provided it, use his ordering to define the segments
188 716 : for (std::size_t ctr=0, s=0, ss=this->segments.size(); s<ss; ctr+=2, ++s)
189 : {
190 216 : const unsigned int index0 = 2*hole_offset+ctr;
191 216 : const unsigned int index1 = 2*hole_offset+ctr+1;
192 :
193 216 : initial.segmentlist[index0] = hole_offset +
194 216 : libmesh_id_to_pointlist_index[this->segments[s].first];
195 216 : initial.segmentlist[index1] = hole_offset +
196 216 : libmesh_id_to_pointlist_index[this->segments[s].second];
197 216 : if (_markers)
198 0 : initial.segmentmarkerlist[hole_offset + s] = (*_markers)[s];
199 : }
200 :
201 :
202 :
203 : // Tell the input struct about the holes
204 500 : if (have_holes)
205 : {
206 16 : initial.numberofholes = _holes->size();
207 16 : initial.holelist = static_cast<REAL*>(std::malloc(initial.numberofholes * 2 * sizeof(REAL)));
208 40 : for (std::size_t i=0, ctr=0, hs=_holes->size(); i<hs; ++i, ctr+=2)
209 : {
210 36 : Point inside_point = (*_holes)[i]->inside();
211 24 : initial.holelist[ctr] = inside_point(0);
212 24 : initial.holelist[ctr+1] = inside_point(1);
213 : }
214 : }
215 :
216 500 : if (_regions)
217 : {
218 0 : initial.numberofregions = _regions->size();
219 0 : initial.regionlist = static_cast<REAL*>(std::malloc(initial.numberofregions * 4 * sizeof(REAL)));
220 0 : for (std::size_t i=0, ctr=0, rs=_regions->size(); i<rs; ++i, ctr+=4)
221 : {
222 0 : Point inside_point = (*_regions)[i]->inside();
223 0 : initial.regionlist[ctr] = inside_point(0);
224 0 : initial.regionlist[ctr+1] = inside_point(1);
225 0 : initial.regionlist[ctr+2] = (*_regions)[i]->attribute();
226 0 : initial.regionlist[ctr+3] = (*_regions)[i]->max_area();
227 : }
228 : }
229 :
230 : // Set the triangulation flags.
231 : // c ~ enclose convex hull with segments
232 : // z ~ use zero indexing
233 : // B ~ Suppresses boundary markers in the output
234 : // Q ~ run in "quiet" mode
235 : // p ~ Triangulates a Planar Straight Line Graph
236 : // If the `p' switch is used, `segmentlist' must point to a list of
237 : // segments, `numberofsegments' must be properly set, and
238 : // `segmentmarkerlist' must either be set to nullptr (in which case all
239 : // markers default to zero), or must point to a list of markers.
240 : // D ~ Conforming Delaunay: use this switch if you want all triangles
241 : // in the mesh to be Delaunay, and not just constrained Delaunay
242 : // q ~ Quality mesh generation with no angles smaller than 20 degrees.
243 : // An alternate minimum angle may be specified after the q
244 : // a ~ Imposes a maximum triangle area constraint.
245 : // -P Suppresses the output .poly file. Saves disk space, but you lose the ability to maintain
246 : // constraining segments on later refinements of the mesh.
247 : // -e Outputs (to an .edge file) a list of edges of the triangulation.
248 : // -v Outputs the Voronoi diagram associated with the triangulation.
249 : // Create the flag strings, depends on element type
250 1000 : std::ostringstream flags;
251 :
252 : // Default flags always used
253 500 : flags << "z";
254 :
255 500 : if (_quiet)
256 500 : flags << "QP";
257 : else
258 0 : flags << "V";
259 :
260 500 : if (_markers)
261 0 : flags << "ev";
262 :
263 : // Flags which are specific to the type of triangulation
264 500 : switch (_triangulation_type)
265 : {
266 460 : case GENERATE_CONVEX_HULL:
267 : {
268 460 : flags << "c";
269 230 : break;
270 : }
271 :
272 40 : case PSLG:
273 : {
274 40 : flags << "p";
275 20 : break;
276 : }
277 :
278 0 : case INVALID_TRIANGULATION_TYPE:
279 0 : libmesh_error_msg("ERROR: INVALID_TRIANGULATION_TYPE selected!");
280 :
281 0 : default:
282 0 : libmesh_error_msg("Unrecognized _triangulation_type");
283 : }
284 :
285 :
286 : // Flags specific to the type of element
287 500 : switch (_elem_type)
288 : {
289 248 : case TRI3:
290 : {
291 : // do nothing.
292 248 : break;
293 : }
294 :
295 4 : case TRI6:
296 : {
297 4 : flags << "o2";
298 2 : break;
299 : }
300 :
301 0 : default:
302 0 : libmesh_error_msg("ERROR: Unrecognized triangular element type == " << Utility::enum_to_string(_elem_type));
303 : }
304 :
305 :
306 : // If we do have holes and the user asked to GENERATE_CONVEX_HULL,
307 : // need to add the p flag so the triangulation respects those segments.
308 500 : if ((_triangulation_type==GENERATE_CONVEX_HULL) && (have_holes))
309 0 : flags << "p";
310 :
311 : // Finally, add the area constraint
312 500 : if (_desired_area > TOLERANCE)
313 500 : flags << "a" << std::fixed << _desired_area;
314 :
315 : // add minimum angle constraint
316 500 : if (_minimum_angle > TOLERANCE)
317 464 : flags << "q" << std::fixed << _minimum_angle;
318 :
319 500 : if (_regions)
320 0 : flags << "Aa";
321 :
322 : // add user provided extra flags
323 500 : if (_extra_flags.size() > 0)
324 0 : flags << _extra_flags;
325 :
326 : // Refine the initial output to conform to the area constraint
327 500 : if (_markers)
328 : {
329 : // need Voronoi to generate boundary information
330 0 : TriangleWrapper::triangulate(const_cast<char *>(flags.str().c_str()),
331 : &initial,
332 : &final,
333 : &voronoi);
334 :
335 : // Send the information computed by Triangle to the Mesh.
336 0 : TriangleWrapper::copy_tri_to_mesh(final,
337 : _mesh,
338 : _elem_type,
339 : &voronoi);
340 : }
341 : else
342 : {
343 500 : TriangleWrapper::triangulate(const_cast<char *>(flags.str().c_str()),
344 : &initial,
345 : &final,
346 : nullptr);
347 : // Send the information computed by Triangle to the Mesh.
348 500 : TriangleWrapper::copy_tri_to_mesh(final,
349 : _mesh,
350 : _elem_type);
351 : }
352 :
353 :
354 500 : _mesh.set_mesh_dimension(2);
355 :
356 : // The user might have requested TRI6 or higher instead of TRI3. If
357 : // so then we'll need to get our neighbor pointers in order ahead of
358 : // time for increase_triangle_order() to use.
359 500 : if (_elem_type != TRI3)
360 : {
361 4 : _mesh.find_neighbors();
362 4 : this->increase_triangle_order();
363 : }
364 :
365 : // To the naked eye, a few smoothing iterations usually looks better,
366 : // so we do this by default unless the user says not to.
367 500 : if (this->_smooth_after_generating)
368 6 : LaplaceMeshSmoother(_mesh).smooth(2);
369 :
370 :
371 : // Clean up.
372 500 : TriangleWrapper::destroy(initial, TriangleWrapper::INPUT);
373 500 : TriangleWrapper::destroy(final, TriangleWrapper::OUTPUT);
374 500 : TriangleWrapper::destroy(voronoi, TriangleWrapper::OUTPUT);
375 :
376 : // Prepare the mesh for use before returning. This ensures (among
377 : // other things) that it is partitioned and therefore users can
378 : // iterate over local elements, etc.
379 500 : _mesh.prepare_for_use();
380 500 : }
381 :
382 : } // namespace libMesh
383 :
384 :
385 :
386 :
387 :
388 :
389 :
390 : #endif // LIBMESH_HAVE_TRIANGLE
|
