libMesh
mesh_tetgen_interface.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/cell_tet4.h"
27 #include "libmesh/face_tri3.h"
28 #include "libmesh/unstructured_mesh.h"
29 #include "libmesh/mesh_tetgen_interface.h"
30 #include "libmesh/utility.h" // binary_find
31 #include "libmesh/mesh_tetgen_wrapper.h"
32 
33 namespace libMesh
34 {
35 
36 //----------------------------------------------------------------------
37 // TetGenMeshInterface class members
39  _mesh(mesh),
40  _serializer(_mesh),
41  _switches("Q")
42 {
43 }
44 
45 void TetGenMeshInterface::set_switches(const std::string & switches)
46 {
47  // set the tetgen switch manually:
48  // p = tetrahedralizes a piecewise linear complex (see definition in user manual)
49  // Q = quiet, no terminal output
50  // q = specify a minimum radius/edge ratio
51  // a = tetrahedron volume constraint
52  // V = verbose output
53  // for full list of options and their meaning: see the tetgen manual
54  // (http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html)
55  _switches = switches;
56 }
57 
58 
60 {
61  // class tetgen_wrapper allows library access on a basic level:
62  TetGenWrapper tetgen_wrapper;
63 
64  // fill input structure with point set data:
65  this->fill_pointlist(tetgen_wrapper);
66 
67  // Run TetGen triangulation method:
68  // Q = quiet, no terminal output
69  // V = verbose, more terminal output
70  // Note: if no switch is used, the input must be a list of 3D points
71  // (.node file) and the Delaunay tetrahedralization of this point set
72  // will be generated.
73 
74  // Can we apply quality and volume constraints in
75  // triangulate_pointset()?. On at least one test problem,
76  // specifying any quality or volume constraints here causes tetgen
77  // to segfault down in the insphere method: a NULL pointer is passed
78  // to the routine.
79  std::ostringstream oss;
80  oss << _switches;
81  // oss << "V"; // verbose operation
82  //oss << "q" << std::fixed << 2.0; // quality constraint
83  //oss << "a" << std::fixed << 100.; // volume constraint
84  tetgen_wrapper.set_switches(oss.str());
85 
86  // Run tetgen
87  tetgen_wrapper.run_tetgen();
88 
89  // save elements to mesh structure, nodes will not be changed:
90  const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
91 
92  // Vector that temporarily holds the node labels defining element.
93  unsigned int node_labels[4];
94 
95  for (unsigned int i=0; i<num_elements; ++i)
96  {
97  Elem * elem = new Tet4;
98 
99  // Get the nodes associated with this element
100  for (unsigned int j=0; j<elem->n_nodes(); ++j)
101  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
102 
103  // Associate the nodes with this element
104  this->assign_nodes_to_elem(node_labels, elem);
105 
106  // Finally, add this element to the mesh.
107  this->_mesh.add_elem(elem);
108  }
109 }
110 
111 
112 
114 {
115  // class tetgen_wrapper allows library access on a basic level
116  TetGenWrapper tetgen_wrapper;
117 
118  // Copy Mesh's node points into TetGen data structure
119  this->fill_pointlist(tetgen_wrapper);
120 
121  // Run TetGen triangulation method:
122  // Q = quiet, no terminal output
123  // Note: if no switch is used, the input must be a list of 3D points
124  // (.node file) and the Delaunay tetrahedralization of this point set
125  // will be generated. In this particular function, we are throwing
126  // away the tetrahedra generated by TetGen, and keeping only the
127  // convex hull...
128  tetgen_wrapper.set_switches(_switches);
129  tetgen_wrapper.run_tetgen();
130  unsigned int num_elements = tetgen_wrapper.get_numberoftrifaces();
131 
132  // Delete *all* old elements. Yes, we legally delete elements while
133  // iterating over them because no entries from the underlying container
134  // are actually erased.
135  for (auto & elem : this->_mesh.element_ptr_range())
136  this->_mesh.delete_elem (elem);
137 
138  // We just removed any boundary info associated with element faces
139  // or edges, so let's update the boundary id caches.
141 
142  // Add the 2D elements which comprise the convex hull back to the mesh.
143  // Vector that temporarily holds the node labels defining element.
144  unsigned int node_labels[3];
145 
146  for (unsigned int i=0; i<num_elements; ++i)
147  {
148  Elem * elem = new Tri3;
149 
150  // Get node labels associated with this element
151  for (unsigned int j=0; j<elem->n_nodes(); ++j)
152  node_labels[j] = tetgen_wrapper.get_triface_node(i,j);
153 
154  this->assign_nodes_to_elem(node_labels, elem);
155 
156  // Finally, add this element to the mesh.
157  this->_mesh.add_elem(elem);
158  }
159 }
160 
161 
162 
163 
164 
166  double volume_constraint)
167 {
168  // start triangulation method with empty holes list:
169  std::vector<Point> noholes;
170  triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
171 }
172 
173 
174 
176  double quality_constraint,
177  double volume_constraint)
178 {
179  // Before calling this function, the Mesh must contain a convex hull
180  // of TRI3 elements which define the boundary.
181  unsigned hull_integrity_check = check_hull_integrity();
182 
183  // Possibly die if hull integrity check failed
184  this->process_hull_integrity_result(hull_integrity_check);
185 
186  // class tetgen_wrapper allows library access on a basic level
187  TetGenWrapper tetgen_wrapper;
188 
189  // Copy Mesh's node points into TetGen data structure
190  this->fill_pointlist(tetgen_wrapper);
191 
192  // >>> fill input structure "tetgenio" with facet data:
193  int facet_num = this->_mesh.n_elem();
194 
195  // allocate memory in "tetgenio" structure:
196  tetgen_wrapper.allocate_facetlist
197  (facet_num, cast_int<int>(holes.size()));
198 
199 
200  // Set up tetgen data structures with existing facet information
201  // from the convex hull.
202  {
203  int insertnum = 0;
204  for (auto & elem : this->_mesh.element_ptr_range())
205  {
206  tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
207  tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);
208 
209  for (unsigned int j=0; j<elem->n_nodes(); ++j)
210  {
211  // We need to get the sequential index of elem->node_ptr(j), but
212  // it should already be stored in _sequential_to_libmesh_node_map...
213  unsigned libmesh_node_id = elem->node_id(j);
214 
215  // The libmesh node IDs may not be sequential, but can we assume
216  // they are at least in order??? We will do so here.
217  std::vector<unsigned>::iterator node_iter =
220  libmesh_node_id);
221 
222  // Check to see if not found: this could also indicate the sequential
223  // node map is not sorted...
224  if (node_iter == _sequential_to_libmesh_node_map.end())
225  libmesh_error_msg("Global node " << libmesh_node_id << " not found in sequential node map!");
226 
227  int sequential_index = cast_int<int>
229  node_iter));
230 
231  // Debugging:
232  // libMesh::out << "libmesh_node_id=" << libmesh_node_id
233  // << ", sequential_index=" << sequential_index
234  // << std::endl;
235 
236  tetgen_wrapper.set_vertex(insertnum, // facet number
237  0, // polygon (always 0)
238  j, // local vertex index in tetgen input
239  sequential_index);
240  }
241 
242  // Go to next facet in polygonlist
243  insertnum++;
244  }
245  }
246 
247 
248 
249  // fill hole list (if there are holes):
250  if (holes.size() > 0)
251  {
252  std::vector<Point>::const_iterator ihole;
253  unsigned hole_index = 0;
254  for (ihole=holes.begin(); ihole!=holes.end(); ++ihole)
255  tetgen_wrapper.set_hole(hole_index++, (*ihole)(0), (*ihole)(1), (*ihole)(2));
256  }
257 
258 
259  // Run TetGen triangulation method
260 
261  // Assemble switches: we append the user's switches (if any) to 'p',
262  // which tetrahedralizes a piecewise linear complex (see definition
263  // in user manual)
264  std::ostringstream oss;
265  oss << "p";
266  oss << _switches;
267 
268  if (quality_constraint != 0)
269  oss << "q" << std::fixed << quality_constraint;
270 
271  if (volume_constraint != 0)
272  oss << "a" << std::fixed << volume_constraint;
273 
274  std::string params = oss.str();
275 
276  tetgen_wrapper.set_switches(params); // TetGen switches: Piecewise linear complex, Quiet mode
277  tetgen_wrapper.run_tetgen();
278 
279  // => nodes:
280  unsigned int old_nodesnum = this->_mesh.n_nodes();
281  REAL x=0., y=0., z=0.;
282  const unsigned int num_nodes = tetgen_wrapper.get_numberofpoints();
283 
284  // Debugging:
285  // libMesh::out << "Original mesh had " << old_nodesnum << " nodes." << std::endl;
286  // libMesh::out << "Reserving space for " << num_nodes << " total nodes." << std::endl;
287 
288  // Reserve space for additional nodes in the node map
289  _sequential_to_libmesh_node_map.reserve(num_nodes);
290 
291  // Add additional nodes to the Mesh.
292  // Original code had i<=num_nodes here (Note: the indexing is:
293  // foo[3*i], [3*i+1], [3*i+2]) But according to the TetGen docs, "In
294  // all cases, the first item in any array is stored starting at
295  // index [0]."
296  for (unsigned int i=old_nodesnum; i<num_nodes; i++)
297  {
298  // Fill in x, y, z values
299  tetgen_wrapper.get_output_node(i, x,y,z);
300 
301  // Catch the node returned by add_point()... this will tell us the ID
302  // assigned by the Mesh.
303  Node * new_node = this->_mesh.add_point ( Point(x,y,z) );
304 
305  // Store this new ID in our sequential-to-libmesh node mapping array
306  _sequential_to_libmesh_node_map.push_back( new_node->id() );
307  }
308 
309  // Debugging:
310  // std::copy(_sequential_to_libmesh_node_map.begin(),
311  // _sequential_to_libmesh_node_map.end(),
312  // std::ostream_iterator<unsigned>(std::cout, " "));
313  // std::cout << std::endl;
314 
315 
316  // => tetrahedra:
317  const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();
318 
319  // Vector that temporarily holds the node labels defining element connectivity.
320  unsigned int node_labels[4];
321 
322  for (unsigned int i=0; i<num_elements; i++)
323  {
324  // TetGen only supports Tet4 elements.
325  Elem * elem = new Tet4;
326 
327  // Fill up the the node_labels vector
328  for (unsigned int j=0; j<elem->n_nodes(); j++)
329  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
330 
331  // Associate nodes with this element
332  this->assign_nodes_to_elem(node_labels, elem);
333 
334  // Finally, add this element to the mesh
335  this->_mesh.add_elem(elem);
336  }
337 
338  // Delete original convex hull elements. Is there ever a case where
339  // we should not do this?
340  this->delete_2D_hull_elements();
341 }
342 
343 
344 
345 
346 
348 {
349  // fill input structure with point set data:
350  wrapper.allocate_pointlist( this->_mesh.n_nodes() );
351 
352  // Make enough space to store a mapping between the implied sequential
353  // node numbering used in tetgen and libmesh's (possibly) non-sequential
354  // numbering scheme.
357 
358  {
359  unsigned index = 0;
360  for (auto & node : this->_mesh.node_ptr_range())
361  {
362  _sequential_to_libmesh_node_map[index] = node->id();
363  wrapper.set_node(index++, (*node)(0), (*node)(1), (*node)(2));
364  }
365  }
366 }
367 
368 
369 
370 
371 
372 void TetGenMeshInterface::assign_nodes_to_elem(unsigned * node_labels, Elem * elem)
373 {
374  for (unsigned int j=0; j<elem->n_nodes(); ++j)
375  {
376  // Get the mapped node index to ask the Mesh for
377  unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];
378 
379  Node * current_node = this->_mesh.node_ptr( mapped_node_id );
380 
381  elem->set_node(j) = current_node;
382  }
383 }
384 
385 
386 
387 
388 
390 {
391  // Check for easy return: if the Mesh is empty (i.e. if
392  // somebody called triangulate_conformingDelaunayMesh on
393  // a Mesh with no elements, then hull integrity check must
394  // fail...
395  if (_mesh.n_elem() == 0)
396  return 3;
397 
398  for (auto & elem : this->_mesh.element_ptr_range())
399  {
400  // Check for proper element type
401  if (elem->type() != TRI3)
402  {
403  //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!");
404  return 1;
405  }
406 
407  for (auto neigh : elem->neighbor_ptr_range())
408  {
409  if (neigh == libmesh_nullptr)
410  {
411  // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
412  return 2;
413  }
414  }
415  }
416 
417  // If we made it here, return success!
418  return 0;
419 }
420 
421 
422 
423 
424 
426 {
427  if (result != 0)
428  {
429  libMesh::err << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
430 
431  if (result==1)
432  {
433  libMesh::err << "Non-TRI3 elements were found in the input Mesh. ";
434  libMesh::err << "A constrained Delaunay triangulation requires a convex hull of TRI3 elements." << std::endl;
435  }
436 
437  libmesh_error_msg("Consider calling TetGenMeshInterface::pointset_convexhull() followed by Mesh::find_neighbors() first.");
438  }
439 }
440 
441 
442 
443 
445 {
446  for (auto & elem : this->_mesh.element_ptr_range())
447  {
448  // Check for proper element type. Yes, we legally delete elements while
449  // iterating over them because no entries from the underlying container
450  // are actually erased.
451  if (elem->type() == TRI3)
452  _mesh.delete_elem(elem);
453  }
454 
455  // We just removed any boundary info associated with hull element
456  // edges, so let's update the boundary id caches.
458 }
459 
460 
461 
462 } // namespace libMesh
463 
464 
465 #endif // #ifdef LIBMESH_HAVE_TETGEN
OStreamProxy err
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:54
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
void set_switches(const std::string &)
Method to set switches to tetgen, allowing for different behaviours.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
int get_element_node(unsigned i, unsigned j)
void get_output_node(unsigned i, REAL &x, REAL &y, REAL &z)
int get_triface_node(unsigned i, unsigned j)
std::string _switches
Parameter controlling the behaviour of tetgen.
void triangulate_conformingDelaunayMesh(double quality_constraint=0., double volume_constraint=0.)
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
std::vector< unsigned > _sequential_to_libmesh_node_map
We should not assume libmesh nodes are numbered sequentially...
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual const Node * node_ptr(const dof_id_type i) const =0
void set_vertex(unsigned i, unsigned j, unsigned k, int nodeindex)
Sets index of ith facet, jth polygon, kth vertex in the TetGen input.
TetGenMeshInterface(UnstructuredMesh &mesh)
Constructor.
The libMesh namespace provides an interface to certain functionality in the library.
void allocate_polygon_vertexlist(unsigned i, unsigned j, int numofvertices)
Allocates memory, sets number of vertices for polygon j, facet i in the TetGen input.
void process_hull_integrity_result(unsigned result)
This function prints an informative message and crashes based on the output of the check_hull_integri...
void allocate_pointlist(int numofpoints)
Allocates memory, sets number of nodes in the TetGen input.
Real distance(const Point &p)
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
virtual unsigned int n_nodes() const =0
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
PetscErrorCode Vec x
void triangulate_pointset()
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
The TetGenWrapper provides an interface for basic access to TetGen data structures and methods...
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
virtual SimpleRange< element_iterator > element_ptr_range()=0
void set_node(unsigned i, REAL x, REAL y, REAL z)
Sets coordinates of point i in the TetGen input.
The UnstructuredMesh class is derived from the MeshBase class.
void triangulate_conformingDelaunayMesh_carvehole(const std::vector< Point > &holes, double quality_constraint=0., double volume_constraint=0.)
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual SimpleRange< node_iterator > node_ptr_range()=0
void assign_nodes_to_elem(unsigned *node_labels, Elem *elem)
Assigns the node IDs contained in the &#39;node_labels&#39; array to &#39;elem&#39;.
void run_tetgen()
Starts the triangulation.
void set_switches(const std::string &s)
Method to set TetGen commandline switches -p Tetrahedralizes a piecewise linear complex (...
void regenerate_id_sets()
Clears and regenerates the cached sets of ids.
void allocate_facet_polygonlist(unsigned i, int numofpolygons)
Allocates memory, sets number of polygons for facet i in the TetGen input.
void set_hole(unsigned i, REAL x, REAL y, REAL z)
Sets coordinates of hole i in the TetGen input.
ForwardIterator binary_find(ForwardIterator first, ForwardIterator last, const T &value)
The STL provides std::binary_search() which returns true or false depending on whether the searched-f...
Definition: utility.h:132
unsigned check_hull_integrity()
This function checks the integrity of the current set of elements in the Mesh to see if they comprise...
void allocate_facetlist(int numoffacets, int numofholes)
Allocates memory, sets number of facets, holes in the TetGen input.
void fill_pointlist(TetGenWrapper &wrapper)
This function copies nodes from the _mesh into TetGen&#39;s pointlist.
void pointset_convexhull()
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...
dof_id_type id() const
Definition: dof_object.h:632
virtual dof_id_type n_nodes() const =0
virtual dof_id_type n_elem() const =0
void delete_2D_hull_elements()
Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38