libMesh
Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
libMesh::TetGenMeshInterface Class Reference

Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen library. More...

#include <mesh_tetgen_interface.h>

Inheritance diagram for libMesh::TetGenMeshInterface:
[legend]

Public Member Functions

 TetGenMeshInterface (UnstructuredMesh &mesh)
 Constructor. More...
 
virtual ~TetGenMeshInterface () override=default
 Empty destructor. More...
 
void set_switches (std::string new_switches)
 Method to set switches to tetgen, allowing for different behaviours. More...
 
virtual void triangulate () override
 Method invokes TetGen library to compute a Delaunay tetrahedralization. More...
 
void triangulate_pointset ()
 Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set. More...
 
void pointset_convexhull ()
 Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set. More...
 
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. More...
 
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. More...
 
Realdesired_volume ()
 Sets and/or gets the desired tetrahedron volume. More...
 
bool & smooth_after_generating ()
 Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid. More...
 
ElemTypeelem_type ()
 Sets and/or gets the desired element type. More...
 
void attach_hole_list (std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh >>> holes)
 Attaches a vector of Mesh pointers defining holes which will be meshed around. More...
 

Protected Member Functions

void fill_pointlist (TetGenWrapper &wrapper)
 This function copies nodes from the _mesh into TetGen's pointlist. More...
 
void assign_nodes_to_elem (unsigned *node_labels, Elem *elem)
 Assigns the node IDs contained in the 'node_labels' array to 'elem'. More...
 
unsigned check_hull_integrity ()
 This function checks the integrity of the current set of elements in the Mesh to see if they comprise a hull, that is: More...
 
void process_hull_integrity_result (unsigned result)
 This function prints an informative message and crashes based on the output of the check_hull_integrity() function. More...
 
void delete_2D_hull_elements ()
 Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization. More...
 

Static Protected Member Functions

static BoundingBox volume_to_surface_mesh (UnstructuredMesh &mesh)
 Remove volume elements from the given mesh, after converting their outer boundary faces to surface elements. More...
 

Protected Attributes

std::vector< unsigned > _sequential_to_libmesh_node_map
 We should not assume libmesh nodes are numbered sequentially... More...
 
MeshSerializer _serializer
 Tetgen only operates on serial meshes. More...
 
std::string _switches
 Parameter controlling the behaviour of tetgen. More...
 
Real _desired_volume
 The desired volume for the elements in the resulting mesh. More...
 
bool _smooth_after_generating
 Flag which tells whether we should smooth the mesh after it is generated. More...
 
ElemType _elem_type
 The exact type of tetrahedra we intend to construct. More...
 
UnstructuredMesh_mesh
 Local reference to the mesh we are working with. More...
 
std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh > > > _holes
 A pointer to a vector of meshes each defining a hole. More...
 

Detailed Description

Class TetGenMeshInterface provides an interface for tetrahedralization of meshes using the TetGen library.

For information about TetGen cf. TetGen home page.

Author
Steffen Petersen
Date
2004
Author
John W. Peterson
Date
2011

Definition at line 55 of file mesh_tetgen_interface.h.

Constructor & Destructor Documentation

◆ TetGenMeshInterface()

libMesh::TetGenMeshInterface::TetGenMeshInterface ( UnstructuredMesh mesh)
explicit

Constructor.

Takes a reference to the mesh.

Definition at line 41 of file mesh_tetgen_interface.C.

41  :
44  _switches("Q")
45 {
46 }
MeshSerializer _serializer
Tetgen only operates on serial meshes.
std::string _switches
Parameter controlling the behaviour of tetgen.
MeshBase & mesh
MeshTetInterface(UnstructuredMesh &mesh)
Constructor.
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.

◆ ~TetGenMeshInterface()

virtual libMesh::TetGenMeshInterface::~TetGenMeshInterface ( )
overridevirtualdefault

Empty destructor.

Member Function Documentation

◆ assign_nodes_to_elem()

void libMesh::TetGenMeshInterface::assign_nodes_to_elem ( unsigned *  node_labels,
Elem elem 
)
protected

Assigns the node IDs contained in the 'node_labels' array to 'elem'.

Definition at line 408 of file mesh_tetgen_interface.C.

References libMesh::MeshTetInterface::_mesh, _sequential_to_libmesh_node_map, libMesh::Elem::node_index_range(), libMesh::MeshBase::node_ptr(), and libMesh::Elem::set_node().

Referenced by pointset_convexhull(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().

409 {
410  for (auto j : elem->node_index_range())
411  {
412  // Get the mapped node index to ask the Mesh for
413  unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];
414 
415  Node * current_node = this->_mesh.node_ptr( mapped_node_id );
416 
417  elem->set_node(j, current_node);
418  }
419 }
std::vector< unsigned > _sequential_to_libmesh_node_map
We should not assume libmesh nodes are numbered sequentially...
virtual const Node * node_ptr(const dof_id_type i) const =0
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.

◆ attach_hole_list()

void libMesh::MeshTetInterface::attach_hole_list ( std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh >>>  holes)
inherited

Attaches a vector of Mesh pointers defining holes which will be meshed around.

We use unique_ptr here because we expect that we may need to modify these meshes internally.

Definition at line 125 of file mesh_tet_interface.C.

Referenced by MeshTetTest::testHole(), and MeshTetTest::testSphereShell().

126 {
127  _holes = std::move(holes);
128 }
std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh > > > _holes
A pointer to a vector of meshes each defining a hole.

◆ check_hull_integrity()

unsigned libMesh::MeshTetInterface::check_hull_integrity ( )
protectedinherited

This function checks the integrity of the current set of elements in the Mesh to see if they comprise a hull, that is:

  • If they are all TRI3 elements
  • They all have non-nullptr neighbors
Returns
  • 0 if the mesh forms a valid hull
  • 1 if a non-TRI3 element is found
  • 2 if an element with a nullptr-neighbor is found

Definition at line 338 of file mesh_tet_interface.C.

References libMesh::MeshTetInterface::_mesh, libMesh::MeshBase::n_elem(), libMesh::Elem::neighbor_ptr_range(), libMesh::TRI3, and libMesh::Elem::type().

Referenced by libMesh::NetGenMeshInterface::triangulate(), and triangulate_conformingDelaunayMesh_carvehole().

339 {
340  // Check for easy return: if the Mesh is empty (i.e. if
341  // somebody called triangulate_conformingDelaunayMesh on
342  // a Mesh with no elements, then hull integrity check must
343  // fail...
344  if (_mesh.n_elem() == 0)
345  return 3;
346 
347  for (auto & elem : this->_mesh.element_ptr_range())
348  {
349  // Check for proper element type
350  if (elem->type() != TRI3)
351  {
352  //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!");
353  return 1;
354  }
355 
356  for (auto neigh : elem->neighbor_ptr_range())
357  {
358  if (neigh == nullptr)
359  {
360  // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
361  return 2;
362  }
363  }
364  }
365 
366  // If we made it here, return success!
367  return 0;
368 }
virtual dof_id_type n_elem() const =0
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.

◆ delete_2D_hull_elements()

void libMesh::MeshTetInterface::delete_2D_hull_elements ( )
protectedinherited

Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization.

Definition at line 390 of file mesh_tet_interface.C.

References libMesh::MeshTetInterface::_mesh, libMesh::MeshBase::delete_elem(), libMesh::MeshBase::get_boundary_info(), libMesh::BoundaryInfo::regenerate_id_sets(), libMesh::TRI3, and libMesh::Elem::type().

Referenced by triangulate_conformingDelaunayMesh_carvehole().

391 {
392  for (auto & elem : this->_mesh.element_ptr_range())
393  {
394  // Check for proper element type. Yes, we legally delete elements while
395  // iterating over them because no entries from the underlying container
396  // are actually erased.
397  if (elem->type() == TRI3)
398  _mesh.delete_elem(elem);
399  }
400 
401  // We just removed any boundary info associated with hull element
402  // edges, so let's update the boundary id caches.
404 }
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
void regenerate_id_sets()
Clears and regenerates the cached sets of ids.
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.

◆ desired_volume()

Real& libMesh::MeshTetInterface::desired_volume ( )
inlineinherited

Sets and/or gets the desired tetrahedron volume.

Set to zero to disable volume constraint.

Definition at line 68 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_desired_volume.

Referenced by main(), and MeshTetTest::testSphereShell().

68 {return _desired_volume;}
Real _desired_volume
The desired volume for the elements in the resulting mesh.

◆ elem_type()

ElemType& libMesh::MeshTetInterface::elem_type ( )
inlineinherited

Sets and/or gets the desired element type.

This should be a Tet type.

Definition at line 81 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_elem_type.

81 {return _elem_type;}
ElemType _elem_type
The exact type of tetrahedra we intend to construct.

◆ fill_pointlist()

void libMesh::TetGenMeshInterface::fill_pointlist ( TetGenWrapper wrapper)
protected

This function copies nodes from the _mesh into TetGen's pointlist.

Takes some pains to ensure that non-sequential node numberings (which can happen with e.g. DistributedMesh) are handled.

Definition at line 380 of file mesh_tetgen_interface.C.

References libMesh::MeshTetInterface::_mesh, _sequential_to_libmesh_node_map, libMesh::TetGenWrapper::allocate_pointlist(), libMesh::MeshBase::n_nodes(), and libMesh::TetGenWrapper::set_node().

Referenced by pointset_convexhull(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().

381 {
382  // fill input structure with point set data:
383  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.
390 
391  {
392  unsigned index = 0;
393  for (auto & node : this->_mesh.node_ptr_range())
394  {
395  _sequential_to_libmesh_node_map[index] = node->id();
396  wrapper.set_node(index++,
397  REAL((*node)(0)),
398  REAL((*node)(1)),
399  REAL((*node)(2)));
400  }
401  }
402 }
std::vector< unsigned > _sequential_to_libmesh_node_map
We should not assume libmesh nodes are numbered sequentially...
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
virtual dof_id_type n_nodes() const =0

◆ pointset_convexhull()

void libMesh::TetGenMeshInterface::pointset_convexhull ( )

Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set.

Stores only 2D hull surface elements.

Definition at line 132 of file mesh_tetgen_interface.C.

References libMesh::MeshTetInterface::_mesh, libMesh::MeshTetInterface::_smooth_after_generating, _switches, libMesh::MeshBase::add_elem(), assign_nodes_to_elem(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), fill_pointlist(), libMesh::MeshBase::get_boundary_info(), libMesh::TetGenWrapper::get_numberoftrifaces(), libMesh::TetGenWrapper::get_triface_node(), libMesh::BoundaryInfo::regenerate_id_sets(), libMesh::TetGenWrapper::run_tetgen(), libMesh::TetGenWrapper::set_switches(), libMesh::LaplaceMeshSmoother::smooth(), and libMesh::TRI3.

Referenced by add_cube_convex_hull_to_mesh().

133 {
134  // class tetgen_wrapper allows library access on a basic level
135  TetGenWrapper tetgen_wrapper;
136 
137  // Copy Mesh's node points into TetGen data structure
138  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  tetgen_wrapper.set_switches(_switches);
148  tetgen_wrapper.run_tetgen();
149  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  for (auto & elem : this->_mesh.element_ptr_range())
155  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.
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  for (unsigned int i=0; i<num_elements; ++i)
166  {
167  auto elem = Elem::build(TRI3);
168 
169  // Get node labels associated with this element
170  for (auto j : elem->node_index_range())
171  node_labels[j] = tetgen_wrapper.get_triface_node(i,j);
172 
173  this->assign_nodes_to_elem(node_labels, elem.get());
174 
175  // Finally, add this element to the mesh.
176  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  if (this->_smooth_after_generating)
182  LaplaceMeshSmoother(this->_mesh).smooth(2);
183 }
std::string _switches
Parameter controlling the behaviour of tetgen.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
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 regenerate_id_sets()
Clears and regenerates the cached sets of ids.
void fill_pointlist(TetGenWrapper &wrapper)
This function copies nodes from the _mesh into TetGen&#39;s pointlist.
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.

◆ process_hull_integrity_result()

void libMesh::MeshTetInterface::process_hull_integrity_result ( unsigned  result)
protectedinherited

This function prints an informative message and crashes based on the output of the check_hull_integrity() function.

It is a separate function so that you can check hull integrity without crashing if you desire.

Definition at line 372 of file mesh_tet_interface.C.

References libMesh::err.

Referenced by triangulate_conformingDelaunayMesh_carvehole().

373 {
374  if (result != 0)
375  {
376  libMesh::err << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
377 
378  if (result==1)
379  {
380  libMesh::err << "Non-TRI3 elements were found in the input Mesh. ";
381  libMesh::err << "A constrained Delaunay triangulation requires a convex hull of TRI3 elements." << std::endl;
382  }
383 
384  libmesh_error_msg("Consider calling TetGenMeshInterface::pointset_convexhull() followed by Mesh::find_neighbors() first.");
385  }
386 }
OStreamProxy err

◆ set_switches()

void libMesh::TetGenMeshInterface::set_switches ( std::string  new_switches)

Method to set switches to tetgen, allowing for different behaviours.

Definition at line 48 of file mesh_tetgen_interface.C.

References _switches.

Referenced by tetrahedralize_domain().

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  _switches = std::move(switches);
59 }
std::string _switches
Parameter controlling the behaviour of tetgen.

◆ smooth_after_generating()

bool& libMesh::MeshTetInterface::smooth_after_generating ( )
inlineinherited

Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid.

False by default (for compatibility with old TetGenMeshInterface behavior).

Definition at line 75 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_smooth_after_generating.

bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.

◆ triangulate()

void libMesh::TetGenMeshInterface::triangulate ( )
overridevirtual

Method invokes TetGen library to compute a Delaunay tetrahedralization.

Implements libMesh::MeshTetInterface.

Definition at line 62 of file mesh_tetgen_interface.C.

References triangulate_pointset().

63 {
64  this->triangulate_pointset();
65 }
void triangulate_pointset()
Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set...

◆ triangulate_conformingDelaunayMesh()

void libMesh::TetGenMeshInterface::triangulate_conformingDelaunayMesh ( double  quality_constraint = 0.,
double  volume_constraint = 0. 
)

Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set.

Boundary constraints are taken from elements array.

Definition at line 189 of file mesh_tetgen_interface.C.

References triangulate_conformingDelaunayMesh_carvehole().

191 {
192  // start triangulation method with empty holes list:
193  std::vector<Point> noholes;
194  triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
195 }
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...

◆ triangulate_conformingDelaunayMesh_carvehole()

void libMesh::TetGenMeshInterface::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.

Boundary constraints are taken from elements array. Include carve-out functionality.

Definition at line 199 of file mesh_tetgen_interface.C.

References libMesh::MeshTetInterface::_mesh, _sequential_to_libmesh_node_map, libMesh::MeshTetInterface::_smooth_after_generating, _switches, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::TetGenWrapper::allocate_facet_polygonlist(), libMesh::TetGenWrapper::allocate_facetlist(), libMesh::TetGenWrapper::allocate_polygon_vertexlist(), assign_nodes_to_elem(), libMesh::Utility::binary_find(), libMesh::Elem::build(), libMesh::MeshTetInterface::check_hull_integrity(), libMesh::MeshTetInterface::delete_2D_hull_elements(), distance(), fill_pointlist(), libMesh::TetGenWrapper::get_element_node(), libMesh::TetGenWrapper::get_numberofpoints(), libMesh::TetGenWrapper::get_numberoftetrahedra(), libMesh::TetGenWrapper::get_output_node(), libMesh::DofObject::id(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshTetInterface::process_hull_integrity_result(), libMesh::TetGenWrapper::run_tetgen(), libMesh::TetGenWrapper::set_hole(), libMesh::TetGenWrapper::set_switches(), libMesh::TetGenWrapper::set_vertex(), libMesh::LaplaceMeshSmoother::smooth(), and libMesh::TET4.

Referenced by tetrahedralize_domain(), and triangulate_conformingDelaunayMesh().

202 {
203  // Before calling this function, the Mesh must contain a convex hull
204  // of TRI3 elements which define the boundary.
205  unsigned hull_integrity_check = check_hull_integrity();
206 
207  // Possibly die if hull integrity check failed
208  this->process_hull_integrity_result(hull_integrity_check);
209 
210  // class tetgen_wrapper allows library access on a basic level
211  TetGenWrapper tetgen_wrapper;
212 
213  // Copy Mesh's node points into TetGen data structure
214  this->fill_pointlist(tetgen_wrapper);
215 
216  // >>> fill input structure "tetgenio" with facet data:
217  int facet_num = this->_mesh.n_elem();
218 
219  // allocate memory in "tetgenio" structure:
220  tetgen_wrapper.allocate_facetlist
221  (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  int insertnum = 0;
228  for (auto & elem : this->_mesh.element_ptr_range())
229  {
230  tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
231  tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);
232 
233  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  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 =
244  libmesh_node_id);
245 
246  // Check to see if not found: this could also indicate the sequential
247  // node map is not sorted...
248  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>
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  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  insertnum++;
268  }
269  }
270 
271 
272 
273  // fill hole list (if there are holes):
274  if (holes.size() > 0)
275  {
276  unsigned hole_index = 0;
277  for (Point ihole : holes)
278  tetgen_wrapper.set_hole(hole_index++,
279  REAL(ihole(0)),
280  REAL(ihole(1)),
281  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  std::ostringstream oss;
293  oss << "pC";
294  oss << _switches;
295 
296  if (quality_constraint != 0)
297  oss << "q" << std::fixed << quality_constraint;
298 
299  if (volume_constraint != 0)
300  oss << "a" << std::fixed << volume_constraint;
301 
302  std::string params = oss.str();
303 
304  tetgen_wrapper.set_switches(params); // TetGen switches: Piecewise linear complex, Quiet mode
305  tetgen_wrapper.run_tetgen();
306 
307  // => nodes:
308  unsigned int old_nodesnum = this->_mesh.n_nodes();
309  REAL x=0., y=0., z=0.;
310  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  _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  for (unsigned int i=old_nodesnum; i<num_nodes; i++)
325  {
326  // Fill in x, y, z values
327  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  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  _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  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  for (unsigned int i=0; i<num_elements; i++)
351  {
352  // TetGen only supports Tet4 elements.
353  auto elem = Elem::build(TET4);
354 
355  // Fill up the the node_labels vector
356  for (auto j : elem->node_index_range())
357  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
358 
359  // Associate nodes with this element
360  this->assign_nodes_to_elem(node_labels, elem.get());
361 
362  // Finally, add this element to the mesh
363  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  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  if (this->_smooth_after_generating)
373  LaplaceMeshSmoother(_mesh).smooth(2);
374 }
unsigned check_hull_integrity()
This function checks the integrity of the current set of elements in the Mesh to see if they comprise...
std::string _switches
Parameter controlling the behaviour of tetgen.
std::vector< unsigned > _sequential_to_libmesh_node_map
We should not assume libmesh nodes are numbered sequentially...
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 Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
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 process_hull_integrity_result(unsigned result)
This function prints an informative message and crashes based on the output of the check_hull_integri...
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:265
void fill_pointlist(TetGenWrapper &wrapper)
This function copies nodes from the _mesh into TetGen&#39;s pointlist.
virtual dof_id_type n_elem() const =0
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
void delete_2D_hull_elements()
Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization.
virtual dof_id_type n_nodes() const =0
bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.

◆ triangulate_pointset()

void libMesh::TetGenMeshInterface::triangulate_pointset ( )

Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set.

Definition at line 68 of file mesh_tetgen_interface.C.

References libMesh::MeshTetInterface::_desired_volume, libMesh::MeshTetInterface::_mesh, libMesh::MeshTetInterface::_smooth_after_generating, _switches, libMesh::MeshBase::add_elem(), assign_nodes_to_elem(), libMesh::Elem::build(), fill_pointlist(), libMesh::TetGenWrapper::get_element_node(), libMesh::TetGenWrapper::get_numberoftetrahedra(), libMesh::TetGenWrapper::run_tetgen(), libMesh::TetGenWrapper::set_switches(), libMesh::LaplaceMeshSmoother::smooth(), and libMesh::TET4.

Referenced by triangulate().

69 {
70  // class tetgen_wrapper allows library access on a basic level:
71  TetGenWrapper tetgen_wrapper;
72 
73  // fill input structure with point set data:
74  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  std::ostringstream oss;
89  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  if (_desired_volume)
96  oss << "a" << std::fixed << _desired_volume; // volume constraint
97 
98  tetgen_wrapper.set_switches(oss.str());
99 
100  // Run tetgen
101  tetgen_wrapper.run_tetgen();
102 
103  // save elements to mesh structure, nodes will not be changed:
104  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  for (unsigned int i=0; i<num_elements; ++i)
110  {
111  auto elem = Elem::build(TET4);
112 
113  // Get the nodes associated with this element
114  for (auto j : elem->node_index_range())
115  node_labels[j] = tetgen_wrapper.get_element_node(i,j);
116 
117  // Associate the nodes with this element
118  this->assign_nodes_to_elem(node_labels, elem.get());
119 
120  // Finally, add this element to the mesh.
121  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  if (this->_smooth_after_generating)
127  LaplaceMeshSmoother(this->_mesh).smooth(2);
128 }
std::string _switches
Parameter controlling the behaviour of tetgen.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
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;.
Real _desired_volume
The desired volume for the elements in the resulting mesh.
void fill_pointlist(TetGenWrapper &wrapper)
This function copies nodes from the _mesh into TetGen&#39;s pointlist.
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
bool _smooth_after_generating
Flag which tells whether we should smooth the mesh after it is generated.

◆ volume_to_surface_mesh()

BoundingBox libMesh::MeshTetInterface::volume_to_surface_mesh ( UnstructuredMesh mesh)
staticprotectedinherited

Remove volume elements from the given mesh, after converting their outer boundary faces to surface elements.

Returns the bounding box of the mesh; this is useful for detecting misplaced holes later.

Definition at line 131 of file mesh_tet_interface.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshTools::Modification::all_tri(), libMesh::Elem::build_side_ptr(), libMesh::MeshBase::delete_elem(), libMesh::Elem::dim(), libMesh::Elem::flip(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::Elem::low_order_key(), libMesh::make_range(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::Elem::n_sides(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_ref_range(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::MeshBase::prepare_for_use(), libMesh::Real, libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_interior_parent(), libMesh::Elem::set_neighbor(), libMesh::DofObject::set_unique_id(), libMesh::Elem::side_index_range(), libMesh::Elem::side_ptr(), libMesh::BoundingBox::union_with(), and libMesh::MeshTools::valid_is_prepared().

Referenced by libMesh::NetGenMeshInterface::triangulate().

132 {
133  // If we've been handed an unprepared mesh then we need to be made
134  // aware of that and fix that; we're relying on neighbor pointers.
136 
137  if (!mesh.is_prepared())
139 
140  // We'll return a bounding box for use by subclasses in basic sanity checks.
141  BoundingBox surface_bb;
142 
143  // First convert all volume boundaries to surface elements; this
144  // gives us a manifold bounding the mesh, though it may not be a
145  // connected manifold even if the volume mesh was connected.
146  {
147  // Make sure ids are in sync and valid on a DistributedMesh
148  const dof_id_type max_orig_id = mesh.max_elem_id();
149 #ifdef LIBMESH_ENABLE_UNIQUE_ID
150  const unique_id_type max_unique_id = mesh.parallel_max_unique_id();
151 #endif
152 
153  // Change this if we add arbitrary polyhedra...
154  const dof_id_type max_sides = 6;
155 
156  std::unordered_set<Elem *> elems_to_delete;
157 
158  std::vector<std::unique_ptr<Elem>> elems_to_add;
159 
160  // Convert all faces to surface elements
161  for (auto * elem : mesh.active_element_ptr_range())
162  {
163  libmesh_error_msg_if (elem->dim() < 2,
164  "Cannot use meshes with 0D or 1D elements to define a volume");
165 
166  // If we've already got 2D elements then those are (part of)
167  // our surface.
168  if (elem->dim() == 2)
169  continue;
170 
171  // 3D elements will be removed after we've extracted their
172  // surface faces.
173  elems_to_delete.insert(elem);
174 
175  for (auto s : make_range(elem->n_sides()))
176  {
177  // If there's a neighbor on this side then there's not a
178  // boundary
179  if (elem->neighbor_ptr(s))
180  {
181  // We're not supporting AMR meshes here yet
182  if (elem->level() != elem->neighbor_ptr(s)->level())
183  libmesh_not_implemented_msg
184  ("Tetrahedralizaton of adapted meshes is not currently supported");
185  continue;
186  }
187 
188  elems_to_add.push_back(elem->build_side_ptr(s));
189  Elem * side_elem = elems_to_add.back().get();
190 
191  // Wipe the interior_parent before it can become a
192  // dangling pointer later
193  side_elem->set_interior_parent(nullptr);
194 
195  // If the mesh is replicated then its automatic id
196  // setting is fine. If not, then we need unambiguous ids
197  // independent of element traversal.
198  if (!mesh.is_replicated())
199  {
200  side_elem->set_id(max_orig_id + max_sides*elem->id() + s);
201 #ifdef LIBMESH_ENABLE_UNIQUE_ID
202  side_elem->set_unique_id(max_unique_id + max_sides*elem->id() + s);
203 #endif
204  }
205  }
206  }
207 
208  // If the mesh is replicated then its automatic neighbor finding
209  // is fine. If not, then we need to insert them ourselves, but
210  // it's easy because we can use the fact (from our implementation
211  // above) that our new elements have no parents or children, plus
212  // the fact (from the tiny fraction of homology I understand) that
213  // a manifold boundary is a manifold with no boundary.
214  //
215  // See UnstructuredMesh::find_neighbors() for more explanation of
216  // (a more complicated version of) the algorithm here.
217  if (!mesh.is_replicated())
218  {
219  typedef dof_id_type key_type;
220  typedef std::pair<Elem *, unsigned char> val_type;
221  typedef std::unordered_multimap<key_type, val_type> map_type;
222  map_type side_to_elem_map;
223 
224  std::unique_ptr<Elem> my_side, their_side;
225 
226  for (auto & elem : elems_to_add)
227  {
228  for (auto s : elem->side_index_range())
229  {
230  if (elem->neighbor_ptr(s))
231  continue;
232  const dof_id_type key = elem->low_order_key(s);
233  auto bounds = side_to_elem_map.equal_range(key);
234  if (bounds.first != bounds.second)
235  {
236  elem->side_ptr(my_side, s);
237  while (bounds.first != bounds.second)
238  {
239  Elem * potential_neighbor = bounds.first->second.first;
240  const unsigned int ns = bounds.first->second.second;
241  potential_neighbor->side_ptr(their_side, ns);
242  if (*my_side == *their_side)
243  {
244  elem->set_neighbor(s, potential_neighbor);
245  potential_neighbor->set_neighbor(ns, elem.get());
246  side_to_elem_map.erase (bounds.first);
247  break;
248  }
249  ++bounds.first;
250  }
251 
252  if (!elem->neighbor_ptr(s))
253  side_to_elem_map.emplace
254  (key, std::make_pair(elem.get(), cast_int<unsigned char>(s)));
255  }
256  }
257  }
258 
259  // At this point we *should* have a match for everything, so
260  // anything we don't have a match for is remote.
261  for (auto & elem : elems_to_add)
262  for (auto s : elem->side_index_range())
263  if (!elem->neighbor_ptr(s))
264  elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
265  }
266 
267  // Remove volume and edge elements
268  for (Elem * elem : elems_to_delete)
269  mesh.delete_elem(elem);
270 
271  // Add the new elements outside the loop so we don't risk
272  // invalidating iterators.
273  for (auto & elem : elems_to_add)
274  mesh.add_elem(std::move(elem));
275  }
276 
277  // Fix up neighbor pointers, element counts, etc.
279 
280  // We're making tets; we need to start with tris
282 
283  // Partition surface into connected components. At this point I'm
284  // finally going to give up and serialize, because at least we got
285  // from 3D down to 2D first, and because I don't want to have to
286  // turn flood_component into a while loop with a parallel sync in
287  // the middle, and because we do have to serialize *eventually*
288  // anyways unless we get a parallel tetrahedralizer backend someday.
289  MeshSerializer mesh_serializer(mesh);
290 
291  std::vector<std::unordered_set<Elem *>> components;
292  std::unordered_set<Elem *> in_component;
293 
294  for (auto * elem : mesh.element_ptr_range())
295  if (!in_component.count(elem))
296  components.emplace_back(flood_component(in_component, elem));
297 
298  const std::unordered_set<Elem *> * biggest_component = nullptr;
299  Real biggest_six_vol = 0;
300  for (const auto & component : components)
301  {
302  Real six_vol = six_times_signed_volume(component);
303  if (std::abs(six_vol) > std::abs(biggest_six_vol))
304  {
305  biggest_six_vol = six_vol;
306  biggest_component = &component;
307  }
308  }
309 
310  if (!biggest_component)
311  libmesh_error_msg("No non-zero-volume component found among " <<
312  components.size() << " boundary components");
313 
314  for (const auto & component : components)
315  if (&component != biggest_component)
316  {
317  for (Elem * elem: component)
318  mesh.delete_elem(elem);
319  }
320  else
321  {
322  for (Elem * elem: component)
323  {
324  if (biggest_six_vol < 0)
325  elem->flip(&mesh.get_boundary_info());
326 
327  for (auto & node : elem->node_ref_range())
328  surface_bb.union_with(node);
329  }
330  }
331 
333 
334  return surface_bb;
335 }
unique_id_type & set_unique_id()
Definition: dof_object.h:858
bool is_prepared() const
Definition: mesh_base.h:198
virtual unique_id_type parallel_max_unique_id() const =0
bool valid_is_prepared(const MeshBase &mesh)
A function for testing whether a mesh&#39;s cached is_prepared() setting is not a false positive...
Definition: mesh_tools.C:1182
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void set_interior_parent(Elem *p)
Sets the pointer to the element&#39;s interior_parent.
Definition: elem.C:1248
dof_id_type & set_id()
Definition: dof_object.h:836
void all_tri(MeshBase &mesh)
Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual dof_id_type max_elem_id() const =0
libmesh_assert(ctx)
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2618
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual std::unique_ptr< Elem > side_ptr(unsigned int i)=0
Temporarily serialize a DistributedMesh for non-distributed-mesh capable code paths.
virtual bool is_replicated() const
Definition: mesh_base.h:233
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
void union_with(const Point &p)
Enlarges this bounding box to include the given point.
Definition: bounding_box.h:286
uint8_t unique_id_type
Definition: id_types.h:86
uint8_t dof_id_type
Definition: id_types.h:67
const RemoteElem * remote_elem
Definition: remote_elem.C:57

Member Data Documentation

◆ _desired_volume

Real libMesh::MeshTetInterface::_desired_volume
protectedinherited

The desired volume for the elements in the resulting mesh.

Unlimited (indicated by 0) by default

Definition at line 138 of file mesh_tet_interface.h.

Referenced by libMesh::MeshTetInterface::desired_volume(), libMesh::NetGenMeshInterface::triangulate(), and triangulate_pointset().

◆ _elem_type

ElemType libMesh::MeshTetInterface::_elem_type
protectedinherited

The exact type of tetrahedra we intend to construct.

Definition at line 149 of file mesh_tet_interface.h.

Referenced by libMesh::MeshTetInterface::elem_type().

◆ _holes

std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh> > > libMesh::MeshTetInterface::_holes
protectedinherited

A pointer to a vector of meshes each defining a hole.

If this is nullptr, there are no holes!

Definition at line 160 of file mesh_tet_interface.h.

Referenced by libMesh::NetGenMeshInterface::triangulate().

◆ _mesh

UnstructuredMesh& libMesh::MeshTetInterface::_mesh
protectedinherited

◆ _sequential_to_libmesh_node_map

std::vector<unsigned> libMesh::TetGenMeshInterface::_sequential_to_libmesh_node_map
protected

We should not assume libmesh nodes are numbered sequentially...

This is not the default behavior of DistributedMesh, for example, unless you specify node IDs explicitly. So this array allows us to keep a mapping between the sequential numbering in tetgen_data.pointlist.

Definition at line 133 of file mesh_tetgen_interface.h.

Referenced by assign_nodes_to_elem(), fill_pointlist(), and triangulate_conformingDelaunayMesh_carvehole().

◆ _serializer

MeshSerializer libMesh::TetGenMeshInterface::_serializer
protected

Tetgen only operates on serial meshes.

Definition at line 138 of file mesh_tetgen_interface.h.

◆ _smooth_after_generating

bool libMesh::MeshTetInterface::_smooth_after_generating
protectedinherited

Flag which tells whether we should smooth the mesh after it is generated.

False by default.

Definition at line 144 of file mesh_tet_interface.h.

Referenced by pointset_convexhull(), libMesh::MeshTetInterface::smooth_after_generating(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().

◆ _switches

std::string libMesh::TetGenMeshInterface::_switches
protected

Parameter controlling the behaviour of tetgen.

By default quiet.

Definition at line 144 of file mesh_tetgen_interface.h.

Referenced by pointset_convexhull(), set_switches(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().


The documentation for this class was generated from the following files: