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 : #include "libmesh/libmesh_config.h"
19 : #ifdef LIBMESH_HAVE_TETGEN
20 :
21 :
22 : // C++ includes
23 : #include <sstream>
24 :
25 : // Local includes
26 : #include "libmesh/mesh_tetgen_interface.h"
27 :
28 : #include "libmesh/boundary_info.h"
29 : #include "libmesh/cell_tet4.h"
30 : #include "libmesh/face_tri3.h"
31 : #include "libmesh/mesh_smoother_laplace.h"
32 : #include "libmesh/unstructured_mesh.h"
33 : #include "libmesh/utility.h" // binary_find
34 : #include "libmesh/mesh_tetgen_wrapper.h"
35 :
36 : namespace libMesh
37 : {
38 :
39 : //----------------------------------------------------------------------
40 : // TetGenMeshInterface class members
41 28 : TetGenMeshInterface::TetGenMeshInterface (UnstructuredMesh & mesh) :
42 : MeshTetInterface(mesh),
43 28 : _serializer(_mesh),
44 28 : _switches("Q")
45 : {
46 28 : }
47 :
48 2 : void TetGenMeshInterface::set_switches(std::string switches)
49 : {
50 : // set the tetgen switch manually:
51 : // p = tetrahedralizes a piecewise linear complex (see definition in user manual)
52 : // Q = quiet, no terminal output
53 : // q = specify a minimum radius/edge ratio
54 : // a = tetrahedron volume constraint
55 : // V = verbose output
56 : // for full list of options and their meaning: see the tetgen manual
57 : // (http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html)
58 2 : _switches = std::move(switches);
59 2 : }
60 :
61 :
62 0 : void TetGenMeshInterface::triangulate ()
63 : {
64 0 : this->triangulate_pointset();
65 0 : }
66 :
67 :
68 1080 : void TetGenMeshInterface::triangulate_pointset ()
69 : {
70 : // class tetgen_wrapper allows library access on a basic level:
71 2160 : TetGenWrapper tetgen_wrapper;
72 :
73 : // fill input structure with point set data:
74 1080 : this->fill_pointlist(tetgen_wrapper);
75 :
76 : // Run TetGen triangulation method:
77 : // Q = quiet, no terminal output
78 : // V = verbose, more terminal output
79 : // Note: if no switch is used, the input must be a list of 3D points
80 : // (.node file) and the Delaunay tetrahedralization of this point set
81 : // will be generated.
82 :
83 : // Can we apply quality and volume constraints in
84 : // triangulate_pointset()?. On at least one test problem,
85 : // specifying any quality or volume constraints here causes tetgen
86 : // to segfault down in the insphere method: a nullptr is passed
87 : // to the routine.
88 2160 : std::ostringstream oss;
89 540 : oss << _switches;
90 : // oss << "V"; // verbose operation
91 : //oss << "q" << std::fixed << 2.0; // quality constraint
92 : //oss << "a" << std::fixed << 100.; // volume constraint
93 :
94 : // But if the user wants refinement, let's do our best.
95 1080 : if (_desired_volume)
96 0 : oss << "a" << std::fixed << _desired_volume; // volume constraint
97 :
98 1080 : tetgen_wrapper.set_switches(oss.str());
99 :
100 : // Run tetgen
101 1080 : tetgen_wrapper.run_tetgen();
102 :
103 : // save elements to mesh structure, nodes will not be changed:
104 1080 : const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
105 :
106 : // Vector that temporarily holds the node labels defining element.
107 : unsigned int node_labels[4];
108 :
109 4554 : for (unsigned int i=0; i<num_elements; ++i)
110 : {
111 3474 : auto elem = Elem::build(TET4);
112 :
113 : // Get the nodes associated with this element
114 17370 : for (auto j : elem->node_index_range())
115 13896 : node_labels[j] = tetgen_wrapper.get_element_node(i,j);
116 :
117 : // Associate the nodes with this element
118 3474 : this->assign_nodes_to_elem(node_labels, elem.get());
119 :
120 : // Finally, add this element to the mesh.
121 6948 : this->_mesh.add_elem(std::move(elem));
122 : }
123 :
124 : // To the naked eye, a few smoothing iterations usually looks better.
125 : // We don't do this by default.
126 1080 : if (this->_smooth_after_generating)
127 0 : LaplaceMeshSmoother(this->_mesh).smooth(2);
128 1080 : }
129 :
130 :
131 :
132 8 : void TetGenMeshInterface::pointset_convexhull ()
133 : {
134 : // class tetgen_wrapper allows library access on a basic level
135 16 : TetGenWrapper tetgen_wrapper;
136 :
137 : // Copy Mesh's node points into TetGen data structure
138 8 : this->fill_pointlist(tetgen_wrapper);
139 :
140 : // Run TetGen triangulation method:
141 : // Q = quiet, no terminal output
142 : // Note: if no switch is used, the input must be a list of 3D points
143 : // (.node file) and the Delaunay tetrahedralization of this point set
144 : // will be generated. In this particular function, we are throwing
145 : // away the tetrahedra generated by TetGen, and keeping only the
146 : // convex hull...
147 8 : tetgen_wrapper.set_switches(_switches);
148 8 : tetgen_wrapper.run_tetgen();
149 8 : unsigned int num_elements = tetgen_wrapper.get_numberoftrifaces();
150 :
151 : // Delete *all* old elements. Yes, we legally delete elements while
152 : // iterating over them because no entries from the underlying container
153 : // are actually erased.
154 20 : for (auto & elem : this->_mesh.element_ptr_range())
155 8 : this->_mesh.delete_elem (elem);
156 :
157 : // We just removed any boundary info associated with element faces
158 : // or edges, so let's update the boundary id caches.
159 12 : this->_mesh.get_boundary_info().regenerate_id_sets();
160 :
161 : // Add the 2D elements which comprise the convex hull back to the mesh.
162 : // Vector that temporarily holds the node labels defining element.
163 : unsigned int node_labels[3];
164 :
165 104 : for (unsigned int i=0; i<num_elements; ++i)
166 : {
167 96 : auto elem = Elem::build(TRI3);
168 :
169 : // Get node labels associated with this element
170 384 : for (auto j : elem->node_index_range())
171 288 : node_labels[j] = tetgen_wrapper.get_triface_node(i,j);
172 :
173 96 : this->assign_nodes_to_elem(node_labels, elem.get());
174 :
175 : // Finally, add this element to the mesh.
176 192 : this->_mesh.add_elem(std::move(elem));
177 : }
178 :
179 : // To the naked eye, a few smoothing iterations usually looks better.
180 : // We don't do this by default.
181 8 : if (this->_smooth_after_generating)
182 0 : LaplaceMeshSmoother(this->_mesh).smooth(2);
183 8 : }
184 :
185 :
186 :
187 :
188 :
189 0 : void TetGenMeshInterface::triangulate_conformingDelaunayMesh (double quality_constraint,
190 : double volume_constraint)
191 : {
192 : // start triangulation method with empty holes list:
193 0 : std::vector<Point> noholes;
194 0 : triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
195 0 : }
196 :
197 :
198 :
199 4 : void TetGenMeshInterface::triangulate_conformingDelaunayMesh_carvehole (const std::vector<Point> & holes,
200 : double quality_constraint,
201 : double volume_constraint)
202 : {
203 : // Before calling this function, the Mesh must contain a convex hull
204 : // of TRI3 elements which define the boundary.
205 4 : unsigned hull_integrity_check = check_hull_integrity();
206 :
207 : // Possibly die if hull integrity check failed
208 4 : this->process_hull_integrity_result(hull_integrity_check);
209 :
210 : // class tetgen_wrapper allows library access on a basic level
211 8 : TetGenWrapper tetgen_wrapper;
212 :
213 : // Copy Mesh's node points into TetGen data structure
214 4 : this->fill_pointlist(tetgen_wrapper);
215 :
216 : // >>> fill input structure "tetgenio" with facet data:
217 4 : int facet_num = this->_mesh.n_elem();
218 :
219 : // allocate memory in "tetgenio" structure:
220 : tetgen_wrapper.allocate_facetlist
221 6 : (facet_num, cast_int<int>(holes.size()));
222 :
223 :
224 : // Set up tetgen data structures with existing facet information
225 : // from the convex hull.
226 : {
227 2 : int insertnum = 0;
228 102 : for (auto & elem : this->_mesh.element_ptr_range())
229 : {
230 96 : tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
231 96 : tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);
232 :
233 432 : for (auto j : elem->node_index_range())
234 : {
235 : // We need to get the sequential index of elem->node_ptr(j), but
236 : // it should already be stored in _sequential_to_libmesh_node_map...
237 432 : unsigned libmesh_node_id = elem->node_id(j);
238 :
239 : // The libmesh node IDs may not be sequential, but can we assume
240 : // they are at least in order??? We will do so here.
241 : std::vector<unsigned>::iterator node_iter =
242 : Utility::binary_find(_sequential_to_libmesh_node_map.begin(),
243 : _sequential_to_libmesh_node_map.end(),
244 288 : libmesh_node_id);
245 :
246 : // Check to see if not found: this could also indicate the sequential
247 : // node map is not sorted...
248 288 : libmesh_error_msg_if(node_iter == _sequential_to_libmesh_node_map.end(),
249 : "Global node " << libmesh_node_id << " not found in sequential node map!");
250 :
251 : int sequential_index = cast_int<int>
252 144 : (std::distance(_sequential_to_libmesh_node_map.begin(),
253 : node_iter));
254 :
255 : // Debugging:
256 : // libMesh::out << "libmesh_node_id=" << libmesh_node_id
257 : // << ", sequential_index=" << sequential_index
258 : // << std::endl;
259 :
260 288 : tetgen_wrapper.set_vertex(insertnum, // facet number
261 : 0, // polygon (always 0)
262 : j, // local vertex index in tetgen input
263 : sequential_index);
264 : }
265 :
266 : // Go to next facet in polygonlist
267 96 : insertnum++;
268 : }
269 : }
270 :
271 :
272 :
273 : // fill hole list (if there are holes):
274 4 : if (holes.size() > 0)
275 : {
276 2 : unsigned hole_index = 0;
277 8 : for (Point ihole : holes)
278 8 : tetgen_wrapper.set_hole(hole_index++,
279 2 : REAL(ihole(0)),
280 2 : REAL(ihole(1)),
281 2 : REAL(ihole(2)));
282 : }
283 :
284 :
285 : // Run TetGen triangulation method
286 :
287 : // Assemble switches: we append the user's switches (if any) to
288 : // - 'p' tetrahedralize a piecewise linear complex
289 : // - 'C' check consistency of mesh (avoid inverted elements)
290 : // (see definition and further options in user manual
291 : // http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html )
292 8 : std::ostringstream oss;
293 4 : oss << "pC";
294 2 : oss << _switches;
295 :
296 4 : if (quality_constraint != 0)
297 2 : oss << "q" << std::fixed << quality_constraint;
298 :
299 4 : if (volume_constraint != 0)
300 2 : oss << "a" << std::fixed << volume_constraint;
301 :
302 4 : std::string params = oss.str();
303 :
304 4 : tetgen_wrapper.set_switches(params); // TetGen switches: Piecewise linear complex, Quiet mode
305 4 : tetgen_wrapper.run_tetgen();
306 :
307 : // => nodes:
308 4 : unsigned int old_nodesnum = this->_mesh.n_nodes();
309 4 : REAL x=0., y=0., z=0.;
310 4 : const unsigned int num_nodes = tetgen_wrapper.get_numberofpoints();
311 :
312 : // Debugging:
313 : // libMesh::out << "Original mesh had " << old_nodesnum << " nodes." << std::endl;
314 : // libMesh::out << "Reserving space for " << num_nodes << " total nodes." << std::endl;
315 :
316 : // Reserve space for additional nodes in the node map
317 4 : _sequential_to_libmesh_node_map.reserve(num_nodes);
318 :
319 : // Add additional nodes to the Mesh.
320 : // Original code had i<=num_nodes here (Note: the indexing is:
321 : // foo[3*i], [3*i+1], [3*i+2]) But according to the TetGen docs, "In
322 : // all cases, the first item in any array is stored starting at
323 : // index [0]."
324 2932 : for (unsigned int i=old_nodesnum; i<num_nodes; i++)
325 : {
326 : // Fill in x, y, z values
327 2928 : tetgen_wrapper.get_output_node(i, x,y,z);
328 :
329 : // Catch the node returned by add_point()... this will tell us the ID
330 : // assigned by the Mesh.
331 4392 : Node * new_node = this->_mesh.add_point ( Point(x,y,z) );
332 :
333 : // Store this new ID in our sequential-to-libmesh node mapping array
334 2928 : _sequential_to_libmesh_node_map.push_back( new_node->id() );
335 : }
336 :
337 : // Debugging:
338 : // std::copy(_sequential_to_libmesh_node_map.begin(),
339 : // _sequential_to_libmesh_node_map.end(),
340 : // std::ostream_iterator<unsigned>(std::cout, " "));
341 : // std::cout << std::endl;
342 :
343 :
344 : // => tetrahedra:
345 4 : const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
346 :
347 : // Vector that temporarily holds the node labels defining element connectivity.
348 : unsigned int node_labels[4];
349 :
350 9980 : for (unsigned int i=0; i<num_elements; i++)
351 : {
352 : // TetGen only supports Tet4 elements.
353 9976 : auto elem = Elem::build(TET4);
354 :
355 : // Fill up the the node_labels vector
356 49880 : for (auto j : elem->node_index_range())
357 39904 : node_labels[j] = tetgen_wrapper.get_element_node(i,j);
358 :
359 : // Associate nodes with this element
360 9976 : this->assign_nodes_to_elem(node_labels, elem.get());
361 :
362 : // Finally, add this element to the mesh
363 19952 : this->_mesh.add_elem(std::move(elem));
364 : }
365 :
366 : // Delete original convex hull elements. Is there ever a case where
367 : // we should not do this?
368 4 : this->delete_2D_hull_elements();
369 :
370 : // To the naked eye, a few smoothing iterations usually looks better.
371 : // We don't do this by default.
372 4 : if (this->_smooth_after_generating)
373 0 : LaplaceMeshSmoother(_mesh).smooth(2);
374 4 : }
375 :
376 :
377 :
378 :
379 :
380 1092 : void TetGenMeshInterface::fill_pointlist(TetGenWrapper & wrapper)
381 : {
382 : // fill input structure with point set data:
383 1092 : wrapper.allocate_pointlist( this->_mesh.n_nodes() );
384 :
385 : // Make enough space to store a mapping between the implied sequential
386 : // node numbering used in tetgen and libmesh's (possibly) non-sequential
387 : // numbering scheme.
388 1092 : _sequential_to_libmesh_node_map.clear();
389 1092 : _sequential_to_libmesh_node_map.resize( this->_mesh.n_nodes() );
390 :
391 : {
392 546 : unsigned index = 0;
393 8014 : for (auto & node : this->_mesh.node_ptr_range())
394 : {
395 6376 : _sequential_to_libmesh_node_map[index] = node->id();
396 12752 : wrapper.set_node(index++,
397 3188 : REAL((*node)(0)),
398 3188 : REAL((*node)(1)),
399 3188 : REAL((*node)(2)));
400 : }
401 : }
402 1092 : }
403 :
404 :
405 :
406 :
407 :
408 13546 : void TetGenMeshInterface::assign_nodes_to_elem(unsigned * node_labels, Elem * elem)
409 : {
410 67634 : for (auto j : elem->node_index_range())
411 : {
412 : // Get the mapped node index to ask the Mesh for
413 54088 : unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];
414 :
415 54088 : Node * current_node = this->_mesh.node_ptr( mapped_node_id );
416 :
417 54088 : elem->set_node(j, current_node);
418 : }
419 13546 : }
420 :
421 : } // namespace libMesh
422 :
423 :
424 : #endif // #ifdef LIBMESH_HAVE_TETGEN
|
