#include <mesh_tetgen_interface.h>

Public Member Functions
	TetGenMeshInterface (UnstructuredMesh &mesh)
	~TetGenMeshInterface ()
void	triangulate_pointset ()
void	pointset_convexhull ()
void	triangulate_conformingDelaunayMesh (double quality_constraint=0., double volume_constraint=0.)
void	triangulate_conformingDelaunayMesh_carvehole (const std::vector< Point > &holes, double quality_constraint=0., double volume_constraint=0.)
Protected Member Functions
void	fill_pointlist (TetGenWrapper &wrapper)
void	assign_nodes_to_elem (unsigned node_labels, Elem elem)
unsigned	check_hull_integrity ()
void	process_hull_integrity_result (unsigned result)
void	delete_2D_hull_elements ()
Protected Attributes
UnstructuredMesh &	_mesh
std::vector< unsigned >	_sequential_to_libmesh_node_map
MeshSerializer	_serializer

Detailed Description

Class TetGenMeshInterface provides an interface for tetrahedrization of meshes using the TetGen library. For information about TetGen cf. TetGen home page.

Author:: , Steffen Petersen, 2004 Refactoring, John W. Peterson, 2011

Definition at line 53 of file mesh_tetgen_interface.h.

Constructor & Destructor Documentation

libMesh::TetGenMeshInterface::TetGenMeshInterface ( UnstructuredMesh & mesh ) [explicit]

Constructor. Takes a reference to the mesh.

Definition at line 38 of file mesh_tetgen_interface.C.

                                                                :
  _mesh(mesh),
  _serializer(_mesh)
{
}

libMesh::TetGenMeshInterface::~TetGenMeshInterface ( ) [inline]

Empty destructor.

Definition at line 66 of file mesh_tetgen_interface.h.

{}

Member Function Documentation

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 367 of file mesh_tetgen_interface.C.

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

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

{
  for (unsigned int j=0; j<elem->n_nodes(); ++j)
    {
      // Get the mapped node index to ask the Mesh for
      unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];

      // Parallel mesh can return NULL pointers, this is bad...
      Node* current_node = this->_mesh.node_ptr( mapped_node_id );

      if (current_node == NULL)
        libmesh_error_msg("Error! Mesh returned NULL node pointer!");

      elem->set_node(j) = current_node;
    }
}

unsigned libMesh::TetGenMeshInterface::check_hull_integrity ( ) [protected]

This function checks the integrity of the current set of elements in the Mesh to see if they comprise a convex hull, that is: .) If they are all TRI3 elements .) They all have non-NULL neighbors

Returns: .) 0 if the mesh forms a valid convex hull .) 1 if a non-TRI3 element is found .) 2 if an element with a NULL-neighbor is found

Definition at line 388 of file mesh_tetgen_interface.C.

References _mesh, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::MeshBase::n_elem(), libMesh::Elem::n_neighbors(), libMesh::Elem::neighbor(), libMesh::TRI3, and libMesh::Elem::type().

Referenced by triangulate_conformingDelaunayMesh_carvehole().

{
  // Check for easy return: if the Mesh is empty (i.e. if
  // somebody called triangulate_conformingDelaunayMesh on
  // a Mesh with no elements, then hull integrity check must
  // fail...
  if (_mesh.n_elem() == 0)
    return 3;

  MeshBase::element_iterator it        = this->_mesh.elements_begin();
  const MeshBase::element_iterator end = this->_mesh.elements_end();

  for (; it != end ; ++it)
    {
      Elem* elem = *it;

      // Check for proper element type
      if (elem->type() != TRI3)
        {
          //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!");
          return 1;
        }

      for (unsigned int i=0; i<elem->n_neighbors(); ++i)
        {
          if (elem->neighbor(i) == NULL)
            {
              // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
              return 2;
            }
        }
    }

  // If we made it here, return success!
  return 0;
}

void libMesh::TetGenMeshInterface::delete_2D_hull_elements ( ) [protected]

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

Definition at line 448 of file mesh_tetgen_interface.C.

References _mesh, libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::TRI3, and libMesh::Elem::type().

Referenced by triangulate_conformingDelaunayMesh_carvehole().

{
  MeshBase::element_iterator it        = this->_mesh.elements_begin();
  const MeshBase::element_iterator end = this->_mesh.elements_end();

  for (; it != end ; ++it)
    {
      Elem* elem = *it;

      // Check for proper element type. Yes, we legally delete elements while
      // iterating over them because no entries from the underlying container
      // are actually erased.
      if (elem->type() == TRI3)
        _mesh.delete_elem(elem);
    }
}

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. ParallelMesh) are handled.

Definition at line 340 of file mesh_tetgen_interface.C.

References _mesh, _sequential_to_libmesh_node_map, libMesh::TetGenWrapper::allocate_pointlist(), end, libMesh::MeshBase::n_nodes(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), and libMesh::TetGenWrapper::set_node().

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

{
  // fill input structure with point set data:
  wrapper.allocate_pointlist( this->_mesh.n_nodes() );

  // Make enough space to store a mapping between the implied sequential
  // node numbering used in tetgen and libmesh's (possibly) non-sequential
  // numbering scheme.
  _sequential_to_libmesh_node_map.clear();
  _sequential_to_libmesh_node_map.resize( this->_mesh.n_nodes() );

  {
    unsigned index = 0;
    MeshBase::node_iterator it  = this->_mesh.nodes_begin();
    const MeshBase::node_iterator end = this->_mesh.nodes_end();
    for ( ; it != end; ++it)
      {
        _sequential_to_libmesh_node_map[index] = (*it)->id();
        wrapper.set_node(index++, (**it)(0), (**it)(1), (**it)(2));
      }
  }
}

void libMesh::TetGenMeshInterface::pointset_convexhull ( )

Method invokes TetGen library to compute a Delaunay tetrahedrization from the nodes point set. Stores only 2D hull surface elements.

Definition at line 100 of file mesh_tetgen_interface.C.

References _mesh, libMesh::MeshBase::add_elem(), assign_nodes_to_elem(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, fill_pointlist(), libMesh::TetGenWrapper::get_numberoftrifaces(), libMesh::TetGenWrapper::get_triface_node(), libMesh::Elem::n_nodes(), libMesh::TetGenWrapper::run_tetgen(), and libMesh::TetGenWrapper::set_switches().

{
  // class tetgen_wrapper allows library access on a basic level
  TetGenWrapper tetgen_wrapper;

  // Copy Mesh's node points into TetGen data structure
  this->fill_pointlist(tetgen_wrapper);

  // Run TetGen triangulation method:
  // Q = quiet, no terminal output
  // Note: if no switch is used, the input must be a list of 3D points
  // (.node file) and the Delaunay tetrahedralization of this point set
  // will be generated.  In this particular function, we are throwing
  // away the tetrahedra generated by TetGen, and keeping only the
  // convex hull...
  tetgen_wrapper.set_switches("Q");
  tetgen_wrapper.run_tetgen();
  unsigned int num_elements   = tetgen_wrapper.get_numberoftrifaces();

  // Delete *all* old elements.  Yes, we legally delete elements while
  // iterating over them because no entries from the underlying container
  // are actually erased.
  {
    MeshBase::element_iterator       it  = this->_mesh.elements_begin();
    const MeshBase::element_iterator end = this->_mesh.elements_end();
    for ( ; it != end; ++it)
      this->_mesh.delete_elem (*it);
  }


  // Add the 2D elements which comprise the convex hull back to the mesh.
  // Vector that temporarily holds the node labels defining element.
  unsigned int node_labels[3];

  for (unsigned int i=0; i<num_elements; ++i)
    {
      Elem* elem = new Tri3;

      // Get node labels associated with this element
      for (unsigned int j=0; j<elem->n_nodes(); ++j)
        node_labels[j] = tetgen_wrapper.get_triface_node(i,j);

      this->assign_nodes_to_elem(node_labels, elem);

      // Finally, add this element to the mesh.
      this->_mesh.add_elem(elem);
    }
}

void libMesh::TetGenMeshInterface::process_hull_integrity_result ( unsigned result ) [protected]

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 429 of file mesh_tetgen_interface.C.

References libMesh::err.

Referenced by triangulate_conformingDelaunayMesh_carvehole().

{
  if (result != 0)
    {
      libMesh::err << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;

      if (result==1)
        {
          libMesh::err << "Non-TRI3 elements were found in the input Mesh.  ";
          libMesh::err << "A constrained Delaunay triangulation requires a convex hull of TRI3 elements." << std::endl;
        }

      libmesh_error_msg("Consider calling TetGenMeshInterface::pointset_convexhull() followed by Mesh::find_neighbors() first.");
    }
}

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

Method invokes TetGen library to compute a Delaunay tetrahedrization from the nodes point set. Boundary constraints are taken from elements array.

Definition at line 153 of file mesh_tetgen_interface.C.

References triangulate_conformingDelaunayMesh_carvehole().

{
  // start triangulation method with empty holes list:
  std::vector<Point> noholes;
  triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
}

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 tetrahedrization from the nodes point set. Boundary constraints are taken from elements array. Include carve-out functionality.

Definition at line 163 of file mesh_tetgen_interface.C.

Referenced by triangulate_conformingDelaunayMesh().

{
  // Before calling this function, the Mesh must contain a convex hull
  // of TRI3 elements which define the boundary.
  unsigned hull_integrity_check = check_hull_integrity();

  // Possibly die if hull integrity check failed
  this->process_hull_integrity_result(hull_integrity_check);

  // class tetgen_wrapper allows library access on a basic level
  TetGenWrapper tetgen_wrapper;

  // Copy Mesh's node points into TetGen data structure
  this->fill_pointlist(tetgen_wrapper);

  // >>> fill input structure "tetgenio" with facet data:
  int facet_num = this->_mesh.n_elem();

  // allocate memory in "tetgenio" structure:
  tetgen_wrapper.allocate_facetlist
    (facet_num, cast_int<int>(holes.size()));


  // Set up tetgen data structures with existing facet information
  // from the convex hull.
  {
    int insertnum = 0;
    MeshBase::element_iterator it        = this->_mesh.elements_begin();
    const MeshBase::element_iterator end = this->_mesh.elements_end();
    for (; it != end ; ++it)
      {
        tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
        tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);

        Elem* elem = *it;

        for (unsigned int j=0; j<elem->n_nodes(); ++j)
          {
            // We need to get the sequential index of elem->get_node(j), but
            // it should already be stored in _sequential_to_libmesh_node_map...
            unsigned libmesh_node_id = elem->node(j);

            // The libmesh node IDs may not be sequential, but can we assume
            // they are at least in order???  We will do so here.
            std::vector<unsigned>::iterator node_iter =
              Utility::binary_find(_sequential_to_libmesh_node_map.begin(),
                                   _sequential_to_libmesh_node_map.end(),
                                   libmesh_node_id);

            // Check to see if not found: this could also indicate the sequential
            // node map is not sorted...
            if (node_iter == _sequential_to_libmesh_node_map.end())
              libmesh_error_msg("Global node " << libmesh_node_id << " not found in sequential node map!");

            int sequential_index = cast_int<int>
              (std::distance(_sequential_to_libmesh_node_map.begin(),
                             node_iter));

            // Debugging:
            //    libMesh::out << "libmesh_node_id=" << libmesh_node_id
            //         << ", sequential_index=" << sequential_index
            //         << std::endl;

            tetgen_wrapper.set_vertex(insertnum, // facet number
                                      0,         // polygon (always 0)
                                      j,         // local vertex index in tetgen input
                                      sequential_index);
          }

        // Go to next facet in polygonlist
        insertnum++;
      }
  }



  // fill hole list (if there are holes):
  if (holes.size() > 0)
    {
      std::vector<Point>::const_iterator ihole;
      unsigned hole_index = 0;
      for (ihole=holes.begin(); ihole!=holes.end(); ++ihole)
        tetgen_wrapper.set_hole(hole_index++, (*ihole)(0), (*ihole)(1), (*ihole)(2));
    }


  // Run TetGen triangulation method:
  // p = tetrahedralizes a piecewise linear complex (see definition in user manual)
  // Q = quiet, no terminal output
  // q = specify a minimum radius/edge ratio
  // a = tetrahedron volume constraint

  // assemble switches:
  std::ostringstream oss; // string holding switches
  oss << "pQ";

  if (quality_constraint != 0)
    oss << "q" << std::fixed << quality_constraint;

  if (volume_constraint != 0)
    oss << "a" << std::fixed << volume_constraint;

  std::string params = oss.str();

  tetgen_wrapper.set_switches(params); // TetGen switches: Piecewise linear complex, Quiet mode
  tetgen_wrapper.run_tetgen();

  // => nodes:
  unsigned int old_nodesnum = this->_mesh.n_nodes();
  REAL x=0., y=0., z=0.;
  const unsigned int num_nodes = tetgen_wrapper.get_numberofpoints();

  // Debugging:
  // libMesh::out << "Original mesh had " << old_nodesnum << " nodes." << std::endl;
  // libMesh::out << "Reserving space for " << num_nodes << " total nodes." << std::endl;

  // Reserve space for additional nodes in the node map
  _sequential_to_libmesh_node_map.reserve(num_nodes);

  // Add additional nodes to the Mesh.
  // Original code had i<=num_nodes here (Note: the indexing is:
  // foo[3*i], [3*i+1], [3*i+2]) But according to the TetGen docs, "In
  // all cases, the first item in any array is stored starting at
  // index [0]."
  for (unsigned int i=old_nodesnum; i<num_nodes; i++)
    {
      // Fill in x, y, z values
      tetgen_wrapper.get_output_node(i, x,y,z);

      // Catch the node returned by add_point()... this will tell us the ID
      // assigned by the Mesh.
      Node* new_node = this->_mesh.add_point ( Point(x,y,z) );

      // Store this new ID in our sequential-to-libmesh node mapping array
      _sequential_to_libmesh_node_map.push_back( new_node->id() );
    }

  // Debugging:
  //  std::copy(_sequential_to_libmesh_node_map.begin(),
  //    _sequential_to_libmesh_node_map.end(),
  //    std::ostream_iterator<unsigned>(std::cout, " "));
  //  std::cout << std::endl;


  // => tetrahedra:
  const unsigned int num_elements = tetgen_wrapper.get_numberoftetrahedra();

  // Vector that temporarily holds the node labels defining element connectivity.
  unsigned int node_labels[4];

  for (unsigned int i=0; i<num_elements; i++)
    {
      // TetGen only supports Tet4 elements.
      Elem* elem = new Tet4;

      // Fill up the the node_labels vector
      for (unsigned int j=0; j<elem->n_nodes(); j++)
        node_labels[j] = tetgen_wrapper.get_element_node(i,j);

      // Associate nodes with this element
      this->assign_nodes_to_elem(node_labels, elem);

      // Finally, add this element to the mesh
      this->_mesh.add_elem(elem);
    }

  // Delete original convex hull elements.  Is there ever a case where
  // we should not do this?
  this->delete_2D_hull_elements();
}

void libMesh::TetGenMeshInterface::triangulate_pointset ( )

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

Definition at line 46 of file mesh_tetgen_interface.C.

References _mesh, libMesh::MeshBase::add_elem(), assign_nodes_to_elem(), fill_pointlist(), libMesh::TetGenWrapper::get_element_node(), libMesh::TetGenWrapper::get_numberoftetrahedra(), libMesh::Elem::n_nodes(), libMesh::TetGenWrapper::run_tetgen(), and libMesh::TetGenWrapper::set_switches().

{
  // class tetgen_wrapper allows library access on a basic level:
  TetGenWrapper tetgen_wrapper;

  // fill input structure with point set data:
  this->fill_pointlist(tetgen_wrapper);

  // Run TetGen triangulation method:
  // Q = quiet, no terminal output
  // V = verbose, more terminal output
  // Note: if no switch is used, the input must be a list of 3D points
  // (.node file) and the Delaunay tetrahedralization of this point set
  // will be generated.

  // Can we apply quality and volume constraints in
  // triangulate_pointset()?.  On at least one test problem,
  // specifying any quality or volume constraints here causes tetgen
  // to segfault down in the insphere method: a NULL pointer is passed
  // to the routine.
  std::ostringstream oss;
  oss << "Q"; // quiet operation
  // oss << "V"; // verbose operation
  //oss  << "q" << std::fixed << 2.0;  // quality constraint
  //oss  << "a" << std::fixed << 100.; // volume constraint
  tetgen_wrapper.set_switches(oss.str());

  // Run tetgen
  tetgen_wrapper.run_tetgen();

  // save elements to mesh structure, nodes will not be changed:
  const unsigned int num_elements   = tetgen_wrapper.get_numberoftetrahedra();

  // Vector that temporarily holds the node labels defining element.
  unsigned int node_labels[4];

  for (unsigned int i=0; i<num_elements; ++i)
    {
      Elem* elem = new Tet4;

      // Get the nodes associated with this element
      for (unsigned int j=0; j<elem->n_nodes(); ++j)
        node_labels[j] = tetgen_wrapper.get_element_node(i,j);

      // Associate the nodes with this element
      this->assign_nodes_to_elem(node_labels, elem);

      // Finally, add this element to the mesh.
      this->_mesh.add_elem(elem);
    }
}

Member Data Documentation

UnstructuredMesh& libMesh::TetGenMeshInterface::_mesh [protected]

Local reference to the mesh we are working with.

Definition at line 145 of file mesh_tetgen_interface.h.

Referenced by assign_nodes_to_elem(), check_hull_integrity(), delete_2D_hull_elements(), fill_pointlist(), pointset_convexhull(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().

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 ParallelMesh, 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 154 of file mesh_tetgen_interface.h.

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

MeshSerializer libMesh::TetGenMeshInterface::_serializer [protected]

Tetgen only operates on serial meshes.

Definition at line 159 of file mesh_tetgen_interface.h.

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

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation