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:

Public Member Functions
	TetGenMeshInterface (UnstructuredMesh &mesh)
	Constructor.
virtual	~TetGenMeshInterface () override=default
	Empty destructor.
void	set_switches (std::string new_switches)
	Method to set switches to tetgen, allowing for different behaviours.
virtual void	triangulate () override
	Method invokes TetGen library to compute a Delaunay tetrahedralization.
void	triangulate_pointset ()
	Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set.
void	pointset_convexhull ()
	Method invokes TetGen library to compute a Delaunay tetrahedralization from the nodes point set.
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.
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.
Real &	desired_volume ()
	Sets and/or gets the desired tetrahedron volume.
bool &	smooth_after_generating ()
	Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid.
ElemType &	elem_type ()
	Sets and/or gets the desired element type.
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.
void	set_verbosity (unsigned int v)
	Sets a verbosity level, defaulting to 0 (print nothing), to be set as high as 100 (print everything).

Protected Types
enum	SurfaceIntegrity {   NON_TRI3 = 1 , MISSING_NEIGHBOR = 2 , EMPTY_MESH = 3 , MISSING_BACKLINK = 4 ,   BAD_NEIGHBOR_NODES = 5 , NON_ORIENTED = 6 , BAD_NEIGHBOR_LINKS = 7 , DEGENERATE_ELEMENT = 8 ,   DEGENERATE_MESH = 9 }
	Enumeration of possible surface mesh integrity issues. More...

Protected Member Functions
void	fill_pointlist (TetGenWrapper &wrapper)
	This function copies nodes from the _mesh into TetGen's pointlist.
void	assign_nodes_to_elem (unsigned *node_labels, Elem *elem)
	Assigns the node IDs contained in the 'node_labels' array to 'elem'.
std::set< SurfaceIntegrity >	check_hull_integrity () const
	This function checks the integrity of the current set of elements in the Mesh to see if they comprise a topological manifold that (if it's also geometrically valid) would define valid boundary for a tetrahedralized volume.
std::set< SurfaceIntegrity >	improve_hull_integrity ()
	This function checks the integrity of the current set of elements in the Mesh, and corrects what it can.
void	process_hull_integrity_result (const std::set< SurfaceIntegrity > &result) const
	This function prints an informative message and throws an exception based on the output of the check_hull_integrity() function.
void	delete_2D_hull_elements ()
	Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization.
void	increase_tet_order ()
	Converts all linear tet elements to the type requested by _elem_type, if that type differs from TET4.

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.

Protected Attributes
std::vector< unsigned >	_sequential_to_libmesh_node_map
	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.
MeshSerializer	_serializer
	Tetgen only operates on serial meshes.
std::string	_switches
	Parameter controlling the behaviour of tetgen.
unsigned int	_verbosity
	verbosity setting
Real	_desired_volume
	The desired volume for the elements in the resulting mesh.
bool	_smooth_after_generating
	Flag which tells whether we should smooth the mesh after it is generated.
ElemType	_elem_type
	The exact type of tetrahedra we intend to construct.
UnstructuredMesh &	_mesh
	Local reference to the mesh we are working with.
std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh > > >	_holes
	A pointer to a vector of meshes each defining a hole.

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.

Member Enumeration Documentation

SurfaceIntegrity

enum libMesh::MeshTetInterface::SurfaceIntegrity

protectedinherited

Enumeration of possible surface mesh integrity issues.

Enumerator
NON_TRI3
MISSING_NEIGHBOR
EMPTY_MESH
MISSING_BACKLINK
BAD_NEIGHBOR_NODES
NON_ORIENTED
BAD_NEIGHBOR_LINKS
DEGENERATE_ELEMENT
DEGENERATE_MESH

Definition at line 124 of file mesh_tet_interface.h.

                        {
    NON_TRI3 = 1,           // a non-TRI3 element is found
    MISSING_NEIGHBOR = 2,   // an element with a nullptr-neighbor is found
    EMPTY_MESH = 3,         // the mesh is empty
    MISSING_BACKLINK = 4,   // an element neighbor isn't linked back to it
    BAD_NEIGHBOR_NODES = 5, // an element neighbor isn't linked to expected nodes
    NON_ORIENTED = 6,       // an element neighbor has inconsistent orientation
    BAD_NEIGHBOR_LINKS = 7, // an element neighbor has other inconsistent links
    DEGENERATE_ELEMENT = 8, // an element has zero area
    DEGENERATE_MESH = 9     // the mesh clearly bounds zero volume
  };

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.

                                                                 :
  MeshTetInterface(mesh),
  _serializer(_mesh),
  _switches("Q")
{
}

~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.

{
  for (auto j : elem->node_index_range())
    {
      // Get the mapped node index to ask the Mesh for
      unsigned mapped_node_id = _sequential_to_libmesh_node_map[ node_labels[j] ];
 
      Node * current_node = this->_mesh.node_ptr( mapped_node_id );
 
      elem->set_node(j, current_node);
    }
}

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().

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 127 of file mesh_tet_interface.C.

{
  _holes = std::move(holes);
}

References libMesh::MeshTetInterface::_holes.

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

check_hull_integrity()

std::set< MeshTetInterface::SurfaceIntegrity > libMesh::MeshTetInterface::check_hull_integrity ( ) const

protectedinherited

This function checks the integrity of the current set of elements in the Mesh to see if they comprise a topological manifold that (if it's also geometrically valid) would define valid boundary for a tetrahedralized volume.

Named check_hull_integrity() for backward compatibility, but now accepts non-convex manifolds.

Returns

a set of enums describing problems found, or an empty set if no problems are found.

Definition at line 341 of file mesh_tet_interface.C.

{
  // 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 {EMPTY_MESH};
 
  std::set<MeshTetInterface::SurfaceIntegrity> returnval;
 
  const BoundingBox bb = MeshTools::create_bounding_box(this->_mesh);
  const Point extents = bb.max() - bb.min();
  if (extents(0) == 0 ||
      extents(1) == 0 ||
      extents(2) == 0)
    returnval.insert(DEGENERATE_MESH);
 
  // Figure a area to use for relative tolerances when detecting
  // degenerate elements
  const Real ref_area = std::abs(extents(0) * extents(1)) +
                        std::abs(extents(0) * extents(2)) +
                        std::abs(extents(1) * extents(2));
 
  struct TriChecker {
    std::set<MeshTetInterface::SurfaceIntegrity> my_returnval;
    const Real my_ref_area;
    const unsigned int my_verbosity;
 
    TriChecker (Real ref_area, unsigned int verbosity) :
      my_returnval(), my_ref_area(ref_area),
      my_verbosity(verbosity) {}
    TriChecker (TriChecker & other, Threads::split) :
      my_returnval(), my_ref_area(other.my_ref_area),
      my_verbosity(other.my_verbosity) {}
 
    void operator()(const ConstElemRange & range) {
 
      for (const Elem * elem : range)
        {
          // Check for proper element type
          if (elem->type() != TRI3)
            {
              if (my_verbosity >= 50)
                std::cerr << "Non-Tri3: " << elem->get_info() << std::endl;
              my_returnval.insert(NON_TRI3);
            }
 
          // Make sure it's a decent element.
          if (elem->volume() < my_ref_area * TOLERANCE * TOLERANCE)
            {
              if (my_verbosity >= 50)
                std::cerr << "Degenerate element: " << elem->get_info() << std::endl;
              my_returnval.insert(DEGENERATE_ELEMENT);
            }
 
          for (auto s : elem->side_index_range())
            {
              const Elem * const neigh = elem->neighbor_ptr(s);
 
              if (neigh == nullptr)
                {
                  if (my_verbosity >= 50)
                    std::cerr << "Element missing neighbor " << s << ": " << elem->get_info() << std::endl;
                  my_returnval.insert(MISSING_NEIGHBOR);
                  continue;
                }
 
              // Make sure our neighbor points back to us
              const unsigned int nn = neigh->which_neighbor_am_i(elem);
 
              if (nn >= 3)
                {
                  if (my_verbosity >= 50)
                    std::cerr << "Element missing backlink " << s << ": " << elem->get_info() << std::endl;
                  my_returnval.insert(MISSING_BACKLINK);
                  continue;
                }
 
              // Our neighbor should have the same the edge nodes we do on
              // the neighboring edgei
              const Node * const n1 = elem->node_ptr(s);
              const Node * const n2 = elem->node_ptr((s+1)%3);
 
              const unsigned int i1 = neigh->local_node(n1->id());
              const unsigned int i2 = neigh->local_node(n2->id());
              if (i1 >= 3 || i2 >= 3)
                {
                  if (my_verbosity >= 50)
                    std::cerr << "Element with bad neighbor " << s << " nodes: " << elem->get_info() << std::endl;
                  my_returnval.insert(BAD_NEIGHBOR_NODES);
                  continue;
                }
 
              // It should have those edge nodes in the opposite order
              // (because they have the same orientation we do)
              if ((i2 + 1)%3 != i1)
                {
                  if (my_verbosity >= 50)
                    std::cerr << "Element orientation mismatch with neighbor " << s << ": " << elem->get_info() << std::endl;
                  my_returnval.insert(NON_ORIENTED);
                  continue;
                }
 
              // And it should have those edge nodes in the expected
              // places relative to its neighbor link
              if (i2 != nn)
                {
                  if (my_verbosity >= 50)
                    std::cerr << "Element with bad links on neighbor " << s << ": " << elem->get_info() << std::endl;
                  my_returnval.insert(BAD_NEIGHBOR_LINKS);
                  continue;
                }
            }
        }
    }
 
    void join(TriChecker & other) {
      my_returnval.merge(other.my_returnval);
    }
  };
 
  TriChecker checker (ref_area, this->_verbosity);
 
  Threads::parallel_reduce
    (this->_mesh.active_local_element_stored_range(), checker);
 
  // Join problems found in threaded loop
  returnval.merge(checker.my_returnval);
 
  // Join problems found on other ranks
  std::set<char> int_set;
  std::transform
    (returnval.begin(), returnval.end(),
     std::inserter<std::set<char>>(int_set, int_set.end()),
     [](SurfaceIntegrity i){return int(i);});
  _mesh.comm().set_union(int_set);
  std::transform
    (int_set.begin(), int_set.end(),
     std::inserter<std::set<SurfaceIntegrity>>(returnval, returnval.end()),
     [](int i){return SurfaceIntegrity(i);});
 
  return returnval;
}

References libMesh::MeshTetInterface::_mesh, libMesh::MeshTetInterface::_verbosity, libMesh::MeshBase::active_local_element_stored_range(), libMesh::ParallelObject::comm(), libMesh::MeshTools::create_bounding_box(), libMesh::MeshTetInterface::DEGENERATE_MESH, libMesh::MeshTetInterface::EMPTY_MESH, libMesh::Elem::get_info(), libMesh::DofObject::id(), libMesh::Elem::local_node(), libMesh::BoundingBox::max(), libMesh::BoundingBox::min(), libMesh::MeshBase::n_elem(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_ptr(), libMesh::Threads::parallel_reduce(), libMesh::Real, libMesh::Parallel::Communicator::set_union(), libMesh::Elem::side_index_range(), libMesh::TOLERANCE, libMesh::TRI3, libMesh::Elem::type(), libMesh::Elem::volume(), and libMesh::Elem::which_neighbor_am_i().

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

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 685 of file mesh_tet_interface.C.

{
  for (auto & elem : this->_mesh.element_ptr_range())
    {
      // 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);
    }
 
  // We just removed any boundary info associated with hull element
  // edges, so let's update the boundary id caches.
  this->_mesh.get_boundary_info().regenerate_id_sets();
}

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().

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.

{return _desired_volume;}

References libMesh::MeshTetInterface::_desired_volume.

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

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.

{return _elem_type;}

References libMesh::MeshTetInterface::_elem_type.

Referenced by MeshTetTest::testNetGenQuadratic(), MeshTetTest::testNetGenQuadraticCurved(), MeshTetTest::testNetGenTet14(), and MeshTetTest::testNetGenTet14Curved().

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.

{
  // 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;
    for (auto & node : this->_mesh.node_ptr_range())
      {
        _sequential_to_libmesh_node_map[index] = node->id();
        wrapper.set_node(index++,
                         REAL((*node)(0)),
                         REAL((*node)(1)),
                         REAL((*node)(2)));
      }
  }
}

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().

improve_hull_integrity()

std::set< MeshTetInterface::SurfaceIntegrity > libMesh::MeshTetInterface::improve_hull_integrity ( )

protectedinherited

This function checks the integrity of the current set of elements in the Mesh, and corrects what it can.

Returns

A set of SurfaceIntegrity codes from check_hull_integrity() if there are problems it can't fix, or an empty set otherwise.

Definition at line 488 of file mesh_tet_interface.C.

{
  // We don't really do anything parallel here, but we aspire to.
  libmesh_parallel_only(this->_mesh.comm());
 
  std::set<MeshTetInterface::SurfaceIntegrity> integrityproblems =
    this->check_hull_integrity();
 
  // If we have no problem, or a problem we can't fix, we're done.
  if (integrityproblems.empty() ||
      integrityproblems.count(NON_TRI3) ||
      integrityproblems.count(EMPTY_MESH))
    return integrityproblems;
 
  // Possibly the user gave us an unprepared mesh with missing or bad
  // neighbor links?
  if (integrityproblems.count(MISSING_NEIGHBOR) ||
      integrityproblems.count(MISSING_BACKLINK) ||
      integrityproblems.count(BAD_NEIGHBOR_LINKS))
  {
    this->_mesh.find_neighbors();
    integrityproblems = this->check_hull_integrity();
  }
 
  // If find_neighbors() doesn't fix these, I give up.
  if (integrityproblems.count(MISSING_NEIGHBOR) ||
      integrityproblems.count(MISSING_BACKLINK) ||
      integrityproblems.count(BAD_NEIGHBOR_LINKS))
    return integrityproblems;
 
  // find_neighbors() might have fixed everything
  if (integrityproblems.empty())
    return integrityproblems;
 
  // A non-oriented (but orientable!) surface is the only thing we
  // shouldn't have fixed or given up on by now.
  libmesh_assert_equal_to(integrityproblems.size(), 1);
  libmesh_assert_equal_to(integrityproblems.count(NON_ORIENTED), 1);
 
  // We need one known-good triangle to start from.  We'll pick the
  // most-negative-x normal among the triangles on the most-negative-x
  // point.
 
  // We'll just implement this in serial for now.
  MeshSerializer mesh_serializer(this->_mesh);
 
  // I don't see why we'd need boundary info here, but maybe we'll
  // want to preserve edge/node conditions eventually?
  BoundaryInfo & bi = this->_mesh.get_boundary_info();
 
  const Node * lowest_point = (*this->_mesh.elements_begin())->node_ptr(0);
 
  // Index by ids, not pointers, for consistency in parallel
  std::unordered_set<dof_id_type> attached_elements;
 
  for (Elem * elem : this->_mesh.element_ptr_range())
    {
      for (const Node & node : elem->node_ref_range())
        {
          if (node(0) < (*lowest_point)(0))
            {
              lowest_point = &node;
              attached_elements.clear();
            }
          if (&node == lowest_point)
            attached_elements.insert(elem->id());
        }
    }
 
  Elem * best_elem = nullptr;
  Real best_abs_normal_0 = 0;
 
  for (dof_id_type id : attached_elements)
    {
      Elem * elem = this->_mesh.elem_ptr(id);
      const Point e01 = elem->point(1) - elem->point(0);
      const Point e02 = elem->point(2) - elem->point(0);
      const Point normal = e01.cross(e02).unit();
      const Real abs_normal_0 = std::abs(normal(0));
 
      if (!best_elem || abs_normal_0 > best_abs_normal_0)
        {
          best_elem = elem;
          best_abs_normal_0 = abs_normal_0;
 
          // Make sure that element is actually a good one, by
          // flipping it if it's not.
          if (abs_normal_0 == normal(0))
            elem->flip(&bi);
        }
    }
 
  // Now flood-fill from that element to get a consistent orientation
  // for the others.
  std::unordered_set<dof_id_type> frontier_elements{best_elem->id()},
                                  finished_elements{};
 
  while (!frontier_elements.empty())
    {
      const dof_id_type elem_id = *frontier_elements.begin();
      Elem & elem = this->_mesh.elem_ref(elem_id);
      for (auto s : elem.side_index_range())
        {
          Elem * neigh = elem.neighbor_ptr(s);
          libmesh_assert(neigh);
          libmesh_assert_less(neigh->which_neighbor_am_i(&elem), 3);
 
          const Node * const n1 = elem.node_ptr(s);
          const Node * const n2 = elem.node_ptr((s+1)%3);
          const unsigned int i1 = neigh->local_node(n1->id());
          const unsigned int i2 = neigh->local_node(n2->id());
          libmesh_assert_less(i1, 3);
          libmesh_assert_less(i2, 3);
 
          const dof_id_type neigh_id = neigh->id();
 
          const bool frontier_neigh = frontier_elements.count(neigh_id);
          const bool finished_neigh = finished_elements.count(neigh_id);
 
          // Are we flipped?
          if ((i2 + 1)%3 != i1)
            {
              // Are we a Moebius strip???  We give up.
              if (frontier_neigh || finished_neigh)
                return integrityproblems;
 
              neigh->flip(&bi);
            }
 
          if (!frontier_neigh && !finished_neigh)
            frontier_elements.insert(neigh_id);
        }
 
      finished_elements.insert(elem_id);
      frontier_elements.erase(elem_id);
    }
 
  this->_mesh.find_neighbors();
 
  libmesh_assert(this->check_hull_integrity().empty());
 
  return {};
}

References libMesh::MeshTetInterface::_mesh, libMesh::MeshTetInterface::BAD_NEIGHBOR_LINKS, libMesh::MeshTetInterface::check_hull_integrity(), libMesh::ParallelObject::comm(), libMesh::TypeVector< T >::cross(), libMesh::MeshBase::elem_ptr(), libMesh::MeshBase::elem_ref(), libMesh::MeshTetInterface::EMPTY_MESH, libMesh::UnstructuredMesh::find_neighbors(), libMesh::Elem::flip(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::libmesh_assert(), libMesh::Elem::local_node(), libMesh::MeshTetInterface::MISSING_BACKLINK, libMesh::MeshTetInterface::MISSING_NEIGHBOR, libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_ptr(), libMesh::Elem::node_ref_range(), libMesh::MeshTetInterface::NON_ORIENTED, libMesh::MeshTetInterface::NON_TRI3, libMesh::Elem::point(), libMesh::Real, libMesh::Elem::side_index_range(), and libMesh::Elem::which_neighbor_am_i().

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

increase_tet_order()

void libMesh::MeshTetInterface::increase_tet_order ( )

protectedinherited

Converts all linear tet elements to the type requested by _elem_type, if that type differs from TET4.

Definition at line 703 of file mesh_tet_interface.C.

{
  switch (_elem_type)
    {
    case TET4:
      return;
    case TET10:
      _mesh.all_second_order();
      break;
    case TET14:
      _mesh.all_complete_order();
      break;
    default:
      libmesh_not_implemented();
    }
}

References libMesh::MeshTetInterface::_elem_type, libMesh::MeshTetInterface::_mesh, libMesh::MeshBase::all_complete_order(), libMesh::MeshBase::all_second_order(), libMesh::TET10, libMesh::TET14, and libMesh::TET4.

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

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.

{
  // 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(_switches);
  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.
  for (auto & elem : this->_mesh.element_ptr_range())
    this->_mesh.delete_elem (elem);
 
  // We just removed any boundary info associated with element faces
  // or edges, so let's update the boundary id caches.
  this->_mesh.get_boundary_info().regenerate_id_sets();
 
  // 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)
    {
      auto elem = Elem::build(TRI3);
 
      // Get node labels associated with this element
      for (auto j : elem->node_index_range())
        node_labels[j] = tetgen_wrapper.get_triface_node(i,j);
 
      this->assign_nodes_to_elem(node_labels, elem.get());
 
      // Finally, add this element to the mesh.
      this->_mesh.add_elem(std::move(elem));
    }
 
  // To the naked eye, a few smoothing iterations usually looks better.
  // We don't do this by default.
  if (this->_smooth_after_generating)
    LaplaceMeshSmoother(this->_mesh, 2).smooth();
}

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().

process_hull_integrity_result()

void libMesh::MeshTetInterface::process_hull_integrity_result ( const std::set< SurfaceIntegrity > &  result ) const

protectedinherited

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

It is a separate function so that you can check hull integrity without exiting or catching an exception if desired.

Definition at line 633 of file mesh_tet_interface.C.

{
  std::ostringstream err_msg;
 
  if (result.empty()) // success
    return;
 
  err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
 
  if (result.count(NON_TRI3))
    {
      err_msg << "At least one non-Tri3 element was found in the input boundary mesh.  ";
      err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
    }
  if (result.count(MISSING_NEIGHBOR))
    {
      err_msg << "At least one triangle without three neighbors was found in the input boundary mesh.  ";
      err_msg << "A constrained Delaunay tetrahedralization boundary must be a triangular manifold without boundary." << std::endl;
    }
  if (result.count(EMPTY_MESH))
    {
      err_msg << "The input boundary mesh was empty!" << std::endl;
      err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
    }
  if (result.count(MISSING_BACKLINK))
    {
      err_msg << "At least one triangle neighbor without a return neighbor link was found in the input boundary mesh.  ";
      err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
    }
  if (result.count(BAD_NEIGHBOR_NODES))
    {
      err_msg << "At least one triangle neighbor without expected node links was found in the input boundary mesh.  ";
      err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
    }
  if (result.count(NON_ORIENTED))
    {
      err_msg << "At least one triangle neighbor with an inconsistent orientation was found in the input boundary mesh.  ";
      err_msg << "A constrained Delaunay tetrahedralization boundary must be an oriented Tri3 mesh." << std::endl;
    }
  if (result.count(BAD_NEIGHBOR_LINKS))
    err_msg << "At least one triangle neighbor with inconsistent node and neighbor links was found in the input boundary mesh." << std::endl;
  if (result.count(DEGENERATE_ELEMENT))
    err_msg << "At least one input triangle is degenerate, with near-zero area relative to the manifold." << std::endl;
  if (result.count(DEGENERATE_MESH))
    err_msg << "The input mesh is degenerate, with zero thickness in at least one direction." << std::endl;
 
  libmesh_error_msg(err_msg.str());
}

References libMesh::MeshTetInterface::BAD_NEIGHBOR_LINKS, libMesh::MeshTetInterface::BAD_NEIGHBOR_NODES, libMesh::MeshTetInterface::DEGENERATE_ELEMENT, libMesh::MeshTetInterface::DEGENERATE_MESH, libMesh::MeshTetInterface::EMPTY_MESH, libMesh::MeshTetInterface::MISSING_BACKLINK, libMesh::MeshTetInterface::MISSING_NEIGHBOR, libMesh::MeshTetInterface::NON_ORIENTED, and libMesh::MeshTetInterface::NON_TRI3.

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

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.

{
  // set the tetgen switch manually:
  // 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
  // V = verbose output
  // for full list of options and their meaning: see the tetgen manual
  // (http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html)
  _switches = std::move(switches);
}

References _switches.

Referenced by tetrahedralize_domain().

set_verbosity()

void libMesh::MeshTetInterface::set_verbosity ( unsigned int  v )

inlineinherited

Sets a verbosity level, defaulting to 0 (print nothing), to be set as high as 100 (print everything).

For verbosity >= 50, print all detected surface mesh integrity issues as they're found. Subclasses may add other output at other verbosity levels.

Definition at line 214 of file mesh_tet_interface.h.

{
  this->_verbosity = v;
}

References libMesh::MeshTetInterface::_verbosity.

Referenced by MeshTetTest::testNetGen().

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.

{return _smooth_after_generating;}

References libMesh::MeshTetInterface::_smooth_after_generating.

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.

{
  this->triangulate_pointset();
}

References triangulate_pointset().

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.

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

References triangulate_conformingDelaunayMesh_carvehole().

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.

{
  // Before calling this function, the Mesh must contain a convex hull
  // of TRI3 elements which define the boundary.
  auto hull_integrity_check = this->improve_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;
    for (auto & elem : this->_mesh.element_ptr_range())
      {
        tetgen_wrapper.allocate_facet_polygonlist(insertnum, 1);
        tetgen_wrapper.allocate_polygon_vertexlist(insertnum, 0, 3);
 
        for (auto j : elem->node_index_range())
          {
            // We need to get the sequential index of elem->node_ptr(j), but
            // it should already be stored in _sequential_to_libmesh_node_map...
            unsigned libmesh_node_id = elem->node_id(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...
            libmesh_error_msg_if(node_iter == _sequential_to_libmesh_node_map.end(),
                                 "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)
    {
      unsigned hole_index = 0;
      for (Point ihole : holes)
        tetgen_wrapper.set_hole(hole_index++,
                                REAL(ihole(0)),
                                REAL(ihole(1)),
                                REAL(ihole(2)));
    }
 
 
  // Run TetGen triangulation method
 
  // Assemble switches: we append the user's switches (if any) to
  //  - 'p'  tetrahedralize a piecewise linear complex
  //  - 'C'  check consistency of mesh (avoid inverted elements)
  //  (see definition and further options in user manual
  //  http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html )
  std::ostringstream oss;
  oss << "pC";
  oss << _switches;
 
  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.
      auto elem = Elem::build(TET4);
 
      // Fill up the the node_labels vector
      for (auto j : elem->node_index_range())
        node_labels[j] = tetgen_wrapper.get_element_node(i,j);
 
      // Associate nodes with this element
      this->assign_nodes_to_elem(node_labels, elem.get());
 
      // Finally, add this element to the mesh
      this->_mesh.add_elem(std::move(elem));
    }
 
  // Delete original convex hull elements.  Is there ever a case where
  // we should not do this?
  this->delete_2D_hull_elements();
 
  // To the naked eye, a few smoothing iterations usually looks better.
  // We don't do this by default.
  if (this->_smooth_after_generating)
    LaplaceMeshSmoother(_mesh, 2).smooth();
}

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::delete_2D_hull_elements(), fill_pointlist(), libMesh::TetGenWrapper::get_element_node(), libMesh::TetGenWrapper::get_numberofpoints(), libMesh::TetGenWrapper::get_numberoftetrahedra(), libMesh::TetGenWrapper::get_output_node(), libMesh::DofObject::id(), libMesh::MeshTetInterface::improve_hull_integrity(), 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().

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.

{
  // 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 nullptr is passed
  // to the routine.
  std::ostringstream oss;
  oss << _switches;
  // oss << "V"; // verbose operation
  //oss  << "q" << std::fixed << 2.0;  // quality constraint
  //oss  << "a" << std::fixed << 100.; // volume constraint
 
  // But if the user wants refinement, let's do our best.
  if (_desired_volume)
    oss  << "a" << std::fixed << _desired_volume; // 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)
    {
      auto elem = Elem::build(TET4);
 
      // Get the nodes associated with this element
      for (auto j : elem->node_index_range())
        node_labels[j] = tetgen_wrapper.get_element_node(i,j);
 
      // Associate the nodes with this element
      this->assign_nodes_to_elem(node_labels, elem.get());
 
      // Finally, add this element to the mesh.
      this->_mesh.add_elem(std::move(elem));
    }
 
  // To the naked eye, a few smoothing iterations usually looks better.
  // We don't do this by default.
  if (this->_smooth_after_generating)
    LaplaceMeshSmoother(this->_mesh, 2).smooth();
}

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().

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 134 of file mesh_tet_interface.C.

{
  // If we've been handed an unprepared mesh then we need to be made
  // aware of that and fix that; we're relying on neighbor pointers.
  libmesh_assert(MeshTools::valid_is_prepared(mesh));
 
  if (!mesh.is_prepared())
    mesh.prepare_for_use();
 
  // We'll return a bounding box for use by subclasses in basic sanity checks.
  BoundingBox surface_bb;
 
  // First convert all volume boundaries to surface elements; this
  // gives us a manifold bounding the mesh, though it may not be a
  // connected manifold even if the volume mesh was connected.
  {
    // Make sure ids are in sync and valid on a DistributedMesh
    const dof_id_type max_orig_id = mesh.max_elem_id();
#ifdef LIBMESH_ENABLE_UNIQUE_ID
    const unique_id_type max_unique_id = mesh.parallel_max_unique_id();
#endif
 
    // Change this if we add arbitrary polyhedra...
    const dof_id_type max_sides = 6;
 
    std::unordered_set<Elem *> elems_to_delete;
 
    std::vector<std::unique_ptr<Elem>> elems_to_add;
 
    // Convert all faces to surface elements
    for (auto * elem : mesh.active_element_ptr_range())
      {
        libmesh_error_msg_if (elem->dim() < 2,
          "Cannot use meshes with 0D or 1D elements to define a volume");
 
        // If we've already got 2D elements then those are (part of)
        // our surface.
        if (elem->dim() == 2)
          continue;
 
        // 3D elements will be removed after we've extracted their
        // surface faces.
        elems_to_delete.insert(elem);
 
        for (auto s : make_range(elem->n_sides()))
          {
            // If there's a neighbor on this side then there's not a
            // boundary
            if (elem->neighbor_ptr(s))
              {
                // We're not supporting AMR meshes here yet
                if (elem->level() != elem->neighbor_ptr(s)->level())
                  libmesh_not_implemented_msg
                    ("Tetrahedralizaton of adapted meshes is not currently supported");
                continue;
              }
 
            elems_to_add.push_back(elem->build_side_ptr(s));
            Elem * side_elem = elems_to_add.back().get();
 
            // Wipe the interior_parent before it can become a
            // dangling pointer later
            side_elem->set_interior_parent(nullptr);
 
            // If the mesh is replicated then its automatic id
            // setting is fine.  If not, then we need unambiguous ids
            // independent of element traversal.
            if (!mesh.is_replicated())
              {
                side_elem->set_id(max_orig_id + max_sides*elem->id() + s);
#ifdef LIBMESH_ENABLE_UNIQUE_ID
                side_elem->set_unique_id(max_unique_id + max_sides*elem->id() + s);
#endif
              }
          }
      }
 
    // If the mesh is replicated then its automatic neighbor finding
    // is fine.  If not, then we need to insert them ourselves, but
    // it's easy because we can use the fact (from our implementation
    // above) that our new elements have no parents or children, plus
    // the fact (from the tiny fraction of homology I understand) that
    // a manifold boundary is a manifold with no boundary.
    //
    // See UnstructuredMesh::find_neighbors() for more explanation of
    // (a more complicated version of) the algorithm here.
    if (!mesh.is_replicated())
      {
        typedef dof_id_type                     key_type;
        typedef std::pair<Elem *, unsigned char> val_type;
        typedef std::unordered_multimap<key_type, val_type> map_type;
        map_type side_to_elem_map;
 
        std::unique_ptr<Elem> my_side, their_side;
 
        for (auto & elem : elems_to_add)
          {
            for (auto s : elem->side_index_range())
              {
                if (elem->neighbor_ptr(s))
                  continue;
                const dof_id_type key = elem->low_order_key(s);
                auto bounds = side_to_elem_map.equal_range(key);
                if (bounds.first != bounds.second)
                  {
                    elem->side_ptr(my_side, s);
                    while (bounds.first != bounds.second)
                      {
                        Elem * potential_neighbor = bounds.first->second.first;
                        const unsigned int ns = bounds.first->second.second;
                        potential_neighbor->side_ptr(their_side, ns);
                        if (*my_side == *their_side)
                          {
                            elem->set_neighbor(s, potential_neighbor);
                            potential_neighbor->set_neighbor(ns, elem.get());
                            side_to_elem_map.erase (bounds.first);
                            break;
                          }
                        ++bounds.first;
                      }
 
                    if (!elem->neighbor_ptr(s))
                      side_to_elem_map.emplace
                        (key, std::make_pair(elem.get(), cast_int<unsigned char>(s)));
                  }
              }
          }
 
        // At this point we *should* have a match for everything, so
        // anything we don't have a match for is remote.
        for (auto & elem : elems_to_add)
          for (auto s : elem->side_index_range())
            if (!elem->neighbor_ptr(s))
              elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
      }
 
    // Remove volume and edge elements
    for (Elem * elem : elems_to_delete)
      mesh.delete_elem(elem);
 
    // Add the new elements outside the loop so we don't risk
    // invalidating iterators.
    for (auto & elem : elems_to_add)
      mesh.add_elem(std::move(elem));
  }
 
  // Fix up neighbor pointers, element counts, etc.
  mesh.prepare_for_use();
 
  // We're making tets; we need to start with tris
  MeshTools::Modification::all_tri(mesh);
 
  // Partition surface into connected components.  At this point I'm
  // finally going to give up and serialize, because at least we got
  // from 3D down to 2D first, and because I don't want to have to
  // turn flood_component into a while loop with a parallel sync in
  // the middle, and because we do have to serialize *eventually*
  // anyways unless we get a parallel tetrahedralizer backend someday.
  MeshSerializer mesh_serializer(mesh);
 
  std::vector<std::unordered_set<Elem *>> components;
  std::unordered_set<Elem *> in_component;
 
  for (auto * elem : mesh.element_ptr_range())
    if (!in_component.count(elem))
      components.emplace_back(flood_component(in_component, elem));
 
  const std::unordered_set<Elem *> * biggest_component = nullptr;
  Real biggest_six_vol = 0;
  for (const auto & component : components)
    {
      Real six_vol = six_times_signed_volume(component);
      if (std::abs(six_vol) > std::abs(biggest_six_vol))
        {
          biggest_six_vol = six_vol;
          biggest_component = &component;
        }
    }
 
  if (!biggest_component)
    libmesh_error_msg("No non-zero-volume component found among " <<
                      components.size() << " boundary components");
 
  for (const auto & component : components)
    if (&component != biggest_component)
      {
        for (Elem * elem: component)
          mesh.delete_elem(elem);
      }
    else
      {
        for (Elem * elem: component)
          {
            if (biggest_six_vol < 0)
              elem->flip(&mesh.get_boundary_info());
 
            for (auto & node : elem->node_ref_range())
              surface_bb.union_with(node);
          }
      }
 
  mesh.prepare_for_use();
 
  return surface_bb;
}

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().

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 184 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 195 of file mesh_tet_interface.h.

Referenced by libMesh::MeshTetInterface::elem_type(), libMesh::MeshTetInterface::increase_tet_order(), and libMesh::NetGenMeshInterface::triangulate().

_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 206 of file mesh_tet_interface.h.

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

_mesh

UnstructuredMesh& libMesh::MeshTetInterface::_mesh

protectedinherited

Local reference to the mesh we are working with.

Definition at line 200 of file mesh_tet_interface.h.

Referenced by assign_nodes_to_elem(), libMesh::MeshTetInterface::check_hull_integrity(), libMesh::MeshTetInterface::delete_2D_hull_elements(), fill_pointlist(), libMesh::MeshTetInterface::improve_hull_integrity(), libMesh::MeshTetInterface::increase_tet_order(), pointset_convexhull(), libMesh::NetGenMeshInterface::triangulate(), triangulate_conformingDelaunayMesh_carvehole(), and triangulate_pointset().

_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 190 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().

_verbosity

unsigned int libMesh::MeshTetInterface::_verbosity

protectedinherited

verbosity setting

Definition at line 110 of file mesh_tet_interface.h.

Referenced by libMesh::MeshTetInterface::check_hull_integrity(), and libMesh::MeshTetInterface::set_verbosity().

---

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

• include/mesh/mesh_tetgen_interface.h
• src/mesh/mesh_tetgen_interface.C