mesh_tetgen_interface.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include "libmesh/libmesh_config.h"
#ifdef LIBMESH_HAVE_TETGEN


// C++ includes
#include <sstream>

// Local includes
#include "libmesh/cell_tet4.h"
#include "libmesh/face_tri3.h"
#include "libmesh/unstructured_mesh.h"
#include "libmesh/mesh_tetgen_interface.h"
#include "libmesh/utility.h" // binary_find
#include "libmesh/mesh_tetgen_wrapper.h"

namespace libMesh
{

//----------------------------------------------------------------------
// TetGenMeshInterface class members
TetGenMeshInterface::TetGenMeshInterface (UnstructuredMesh& mesh) :
  _mesh(mesh),
  _serializer(_mesh)
{
}



void TetGenMeshInterface::triangulate_pointset ()
{
  // 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);
    }
}



void TetGenMeshInterface::pointset_convexhull ()
{
  // 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 TetGenMeshInterface::triangulate_conformingDelaunayMesh (double quality_constraint,
                                                              double volume_constraint)
{
  // start triangulation method with empty holes list:
  std::vector<Point> noholes;
  triangulate_conformingDelaunayMesh_carvehole(noholes, quality_constraint, volume_constraint);
}



void TetGenMeshInterface::triangulate_conformingDelaunayMesh_carvehole  (const std::vector<Point>& holes,
                                                                         double quality_constraint,
                                                                         double volume_constraint)
{
  // 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 TetGenMeshInterface::fill_pointlist(TetGenWrapper& wrapper)
{
  // 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 TetGenMeshInterface::assign_nodes_to_elem(unsigned* node_labels, Elem* elem)
{
  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 TetGenMeshInterface::check_hull_integrity()
{
  // 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 TetGenMeshInterface::process_hull_integrity_result(unsigned result)
{
  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 TetGenMeshInterface::delete_2D_hull_elements()
{
  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);
    }
}



} // namespace libMesh


#endif // #ifdef LIBMESH_HAVE_TETGEN