libMesh::NetGenMeshInterface Class Reference

Class NetGenMeshInterface provides an interface for tetrahedralization of meshes using the NetGen library. More...

#include <mesh_netgen_interface.h>

Inheritance diagram for libMesh::NetGenMeshInterface:

Public Member Functions
	NetGenMeshInterface (UnstructuredMesh &mesh)
	Constructor. More...
virtual	~NetGenMeshInterface () override=default
	Empty destructor. More...
virtual void	triangulate () override
	Method invokes NetGen library to compute a tetrahedralization. More...
Real &	desired_volume ()
	Sets and/or gets the desired tetrahedron volume. More...
bool &	smooth_after_generating ()
	Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid. More...
ElemType &	elem_type ()
	Sets and/or gets the desired element type. More...
void	attach_hole_list (std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh >>> holes)
	Attaches a vector of Mesh pointers defining holes which will be meshed around. More...

Protected Member Functions
unsigned	check_hull_integrity ()
	This function checks the integrity of the current set of elements in the Mesh to see if they comprise a hull, that is: More...
void	process_hull_integrity_result (unsigned result)
	This function prints an informative message and crashes based on the output of the check_hull_integrity() function. More...
void	delete_2D_hull_elements ()
	Delete original convex hull elements from the Mesh after performing a Delaunay tetrahedralization. More...

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

Protected Attributes
MeshSerializer	_serializer
	Tetgen only operates on serial meshes. More...
Real	_desired_volume
	The desired volume for the elements in the resulting mesh. More...
bool	_smooth_after_generating
	Flag which tells whether we should smooth the mesh after it is generated. More...
ElemType	_elem_type
	The exact type of tetrahedra we intend to construct. More...
UnstructuredMesh &	_mesh
	Local reference to the mesh we are working with. More...
std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh > > >	_holes
	A pointer to a vector of meshes each defining a hole. More...

Detailed Description

Class NetGenMeshInterface provides an interface for tetrahedralization of meshes using the NetGen library.

For information about TetGen cf. NetGen github repository.

Author

Roy H. Stogner

Date

2024

Definition at line 52 of file mesh_netgen_interface.h.

Constructor & Destructor Documentation

NetGenMeshInterface()

libMesh::NetGenMeshInterface::NetGenMeshInterface ( UnstructuredMesh &  mesh )

explicit

Constructor.

Takes a reference to the mesh containing the triangulated surface which is to be tetrahedralized.

Definition at line 74 of file mesh_netgen_interface.C.

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

~NetGenMeshInterface()

virtual libMesh::NetGenMeshInterface::~NetGenMeshInterface ( )

overridevirtualdefault

Empty destructor.

Member Function Documentation

attach_hole_list()

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

inherited

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

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

Definition at line 125 of file mesh_tet_interface.C.

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

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

check_hull_integrity()

unsigned libMesh::MeshTetInterface::check_hull_integrity ( )

protectedinherited

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

• If they are all TRI3 elements
• They all have non-nullptr neighbors

Returns

• 0 if the mesh forms a valid hull
• 1 if a non-TRI3 element is found
• 2 if an element with a nullptr-neighbor is found

Definition at line 338 of file mesh_tet_interface.C.

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

Referenced by triangulate(), and libMesh::TetGenMeshInterface::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;

  for (auto & elem : this->_mesh.element_ptr_range())
    {
      // 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 (auto neigh : elem->neighbor_ptr_range())
        {
          if (neigh == nullptr)
            {
              // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
              return 2;
            }
        }
    }

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

delete_2D_hull_elements()

void libMesh::MeshTetInterface::delete_2D_hull_elements ( )

protectedinherited

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

Definition at line 390 of file mesh_tet_interface.C.

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

Referenced by libMesh::TetGenMeshInterface::triangulate_conformingDelaunayMesh_carvehole().

{
  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();
}

desired_volume()

Real& libMesh::MeshTetInterface::desired_volume ( )

inlineinherited

Sets and/or gets the desired tetrahedron volume.

Set to zero to disable volume constraint.

Definition at line 68 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_desired_volume.

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

{return _desired_volume;}

elem_type()

ElemType& libMesh::MeshTetInterface::elem_type ( )

inlineinherited

Sets and/or gets the desired element type.

This should be a Tet type.

Definition at line 81 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_elem_type.

{return _elem_type;}

process_hull_integrity_result()

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

protectedinherited

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

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

Definition at line 372 of file mesh_tet_interface.C.

References libMesh::err.

Referenced by libMesh::TetGenMeshInterface::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.");
    }
}

smooth_after_generating()

bool& libMesh::MeshTetInterface::smooth_after_generating ( )

inlineinherited

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

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

Definition at line 75 of file mesh_tet_interface.h.

References libMesh::MeshTetInterface::_smooth_after_generating.

{return _smooth_after_generating;}

triangulate()

void libMesh::NetGenMeshInterface::triangulate ( )

overridevirtual

Method invokes NetGen library to compute a tetrahedralization.

Implements libMesh::MeshTetInterface.

Definition at line 82 of file mesh_netgen_interface.C.

References libMesh::MeshTetInterface::_desired_volume, libMesh::MeshTetInterface::_holes, libMesh::MeshTetInterface::_mesh, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshCommunication::broadcast(), libMesh::Elem::build_with_id(), libMesh::MeshTetInterface::check_hull_integrity(), libMesh::BoundingBox::contains(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::DofObject::invalid_id, libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::make_range(), libMesh::MeshTools::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::Utility::pow(), libMesh::MeshBase::prepare_for_use(), libMesh::ParallelObject::processor_id(), libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and libMesh::MeshTetInterface::volume_to_surface_mesh().

Referenced by main().

{
  using namespace nglib;

  LOG_SCOPE("triangulate()", "NetGenMeshInterface");

  // We're hoping to do volume_to_surface_mesh in parallel at least,
  // but then we'll need to serialize any hole meshes to rank 0 so it
  // can use them in serial.

  const BoundingBox mesh_bb =
    MeshTetInterface::volume_to_surface_mesh(this->_mesh);

  std::vector<MeshSerializer> hole_serializers;
  if (_holes)
    for (std::unique_ptr<UnstructuredMesh> & hole : *_holes)
      {
        const BoundingBox hole_bb =
          MeshTetInterface::volume_to_surface_mesh(*hole);

        libmesh_error_msg_if
          (!mesh_bb.contains(hole_bb),
           "Found hole with bounding box " << hole_bb <<
           "\nextending outside of mesh bounding box " << mesh_bb);

        hole_serializers.emplace_back
          (*hole, /* need_serial */ true,
           /* serial_only_needed_on_proc_0 */ true);
      }

  // If we're not rank 0, we're just going to wait for rank 0 to call
  // Netgen, then receive its data afterward, we're not going to hope
  // that Netgen does the exact same thing on every processor.
  if (this->_mesh.processor_id() != 0)
    {
      // We don't need our holes anymore.  Delete their serializers
      // first to avoid dereferencing dangling pointers.
      hole_serializers.clear();
      if (_holes)
        _holes->clear();

      // Receive the mesh data rank 0 will send later, then fix it up
      // together
      MeshCommunication().broadcast(this->_mesh);
      this->_mesh.prepare_for_use();
      return;
    }

  this->check_hull_integrity();

  Ng_Meshing_Parameters params;

  // Override any default parameters we might need to, to avoid
  // inserting nodes we don't want.
  params.uselocalh = false;
  params.minh = 0;
  params.elementsperedge = 1;
  params.elementspercurve = 1;
  params.closeedgeenable = false;
  params.closeedgefact = 0;
  params.minedgelenenable = false;
  params.minedgelen = 0;

  // Try to get a no-extra-nodes mesh if we're asked to, or try to
  // translate our desired volume into NetGen terms otherwise.
  //
  // Spoiler alert: all we can do is try; NetGen uses a marching front
  // algorithm that can insert extra nodes despite all my best
  // efforts.
  if (_desired_volume == 0) // shorthand for "no refinement"
    {
      params.maxh = std::numeric_limits<double>::max();
      params.fineness = 0; // "coarse" in the docs
      params.grading = 1;  // "aggressive local grading" to avoid smoothing??

      // Turning off optimization steps avoids another opportunity for
      // Netgen to try to add more nodes.
      params.optsteps_3d = 0;
    }
  else
    params.maxh = double(std::pow(_desired_volume, 1./3.));

  // Keep track of how NetGen copies of nodes map back to our original
  // nodes, so we can connect new elements to nodes correctly.
  std::unordered_map<int, dof_id_type> ng_to_libmesh_id;

  auto handle_ng_result = [](Ng_Result result) {
    static const std::vector<std::string> result_types =
      {"Netgen error", "Netgen success", "Netgen surface input error",
       "Netgen volume failure", "Netgen STL input error",
       "Netgen surface failure", "Netgen file not found"};

    if (result+1 >= 0 &&
        std::size_t(result+1) < result_types.size())
      libmesh_error_msg_if
        (result, "Ng_GenerateVolumeMesh failed: " <<
         result_types[result+1]);
    else
      libmesh_error_msg
        ("Ng_GenerateVolumeMesh failed with an unknown error code");
  };

  // Keep track of what boundary ids we want to assign to each new
  // triangle.  We'll give the outer boundary BC 0, and give holes ids
  // starting from 1.
  // We key on sorted tuples of node ids to identify a side.
  std::unordered_map<std::array<dof_id_type,3>,
                     boundary_id_type, libMesh::hash> side_boundary_id;

  auto insert_id = []
    (std::array<dof_id_type,3> & array,
     dof_id_type n_id)
  {
    libmesh_assert_less(n_id, DofObject::invalid_id);
    unsigned int i=0;
    while (array[i] < n_id)
      ++i;
    while (i < 3)
      std::swap(array[i++], n_id);
  };

  WrappedNgMesh ngmesh;

  // Create surface mesh in the WrappedNgMesh
  {
    // NetGen appears to use ONE-BASED numbering for its nodes, and
    // since it doesn't return an id when adding nodes we'll have to
    // track the numbering ourselves.
    int ng_id = 1;

    auto create_surface_component =
      [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
      (UnstructuredMesh & srcmesh, bool hole_mesh,
       boundary_id_type bcid)
    {
      LOG_SCOPE("create_surface_component()", "NetGenMeshInterface");

      // Keep track of what nodes we've already added to the Netgen
      // mesh vs what nodes we need to add.  We'll keep track by id,
      // not by point location.  I don't know if Netgen can handle
      // multiple nodes with the same point location, but if they can
      // it's not going to be *us* who breaks that feature.
      std::unordered_map<dof_id_type, int> libmesh_to_ng_id;

      // Keep track of what nodes we've already added to the main
      // mesh from a hole mesh.
      std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;

      // Use a separate array for passing points to NetGen, just in case
      // we're not using double-precision ourselves.
      std::array<double, 3> point_val;

      // And an array for element vertices
      std::array<int, 3> elem_nodes;

      for (const auto * elem : srcmesh.element_ptr_range())
        {
          // If someone has triangles we can't triangulate, we have a
          // problem
          if (elem->type() == TRI6 ||
              elem->type() == TRI7)
            libmesh_not_implemented_msg
              ("Netgen tetrahedralization currently only supports TRI3 boundaries");

          // If someone has non-triangles, let's just ignore them.
          if (elem->type() != TRI3)
            continue;

          std::array<dof_id_type,3> sorted_ids =
            {DofObject::invalid_id, DofObject::invalid_id,
             DofObject::invalid_id};

          for (int ni : make_range(3))
            {
              // Just using the "invert_trigs" option in NetGen params
              // doesn't work for me, so we'll have to have properly
              // oriented the tris earlier.
              auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];

              const Node & n = elem->node_ref(ni);
              auto n_id = n.id();
              if (hole_mesh)
                {
                  if (auto it = hole_to_main_mesh_id.find(n_id);
                      it != hole_to_main_mesh_id.end())
                    {
                      n_id = it->second;
                    }
                  else
                    {
                      Node * n_new = this->_mesh.add_point(n);
                      const dof_id_type n_new_id = n_new->id();
                      hole_to_main_mesh_id.emplace(n_id, n_new_id);
                      n_id = n_new_id;
                    }
                }

              if (auto it = libmesh_to_ng_id.find(n_id);
                  it != libmesh_to_ng_id.end())
                {
                  const int existing_ng_id = it->second;
                  elem_node = existing_ng_id;
                }
              else
                {
                  for (auto i : make_range(3))
                    point_val[i] = double(n(i));

                  Ng_AddPoint(ngmesh, point_val.data());

                  ng_to_libmesh_id[ng_id] = n_id;
                  libmesh_to_ng_id[n_id] = ng_id;
                  elem_node = ng_id;
                  ++ng_id;
                }

              insert_id(sorted_ids, n_id);
            }

          side_boundary_id[sorted_ids] = bcid;

          Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
        }
    };

    // Number the outer boundary 0, and the holes starting from 1
    boundary_id_type bcid = 0;

    create_surface_component(this->_mesh, false, bcid);

    if (_holes)
      for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
        create_surface_component(*h, true, ++bcid);
  }

  {
    LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");

    auto result = Ng_GenerateVolumeMesh(ngmesh, &params);
    handle_ng_result(result);
  }

  const int n_elem = Ng_GetNE(ngmesh);

  libmesh_error_msg_if (n_elem <= 0,
                        "NetGen failed to generate any tetrahedra");

  const dof_id_type n_points = Ng_GetNP(ngmesh);
  const dof_id_type old_nodes = this->_mesh.n_nodes();

  // Netgen may have generated new interior nodes
  if (n_points != old_nodes)
    {
      std::array<double, 3> point_val;

      // We should only be getting new nodes if we asked for them
      if (!_desired_volume)
        {
          std::cout <<
            "NetGen output " << n_points <<
            " points when we gave it " <<
            old_nodes << " and disabled refinement\n" <<
            "If new interior points are acceptable in your mesh, please set\n" <<
            "a non-zero desired_volume to indicate that.  If new interior\n" <<
            "points are not acceptable in your mesh, you may need a different\n" <<
            "(non-advancing-front?) mesh generator." << std::endl;
          libmesh_error();
        }
      else
        for (auto i : make_range(old_nodes, n_points))
          {
            // i+1 since ng uses ONE-BASED numbering
            Ng_GetPoint (ngmesh, i+1, point_val.data());
            const Point p(point_val[0], point_val[1], point_val[2]);
            Node * n_new = this->_mesh.add_point(p);
            const dof_id_type n_new_id = n_new->id();
            ng_to_libmesh_id[i+1] = n_new_id;
          }
    }

  for (auto * elem : this->_mesh.element_ptr_range())
    this->_mesh.delete_elem(elem);

  BoundaryInfo * bi = & this->_mesh.get_boundary_info();

  for (auto i : make_range(n_elem))
  {
    // Enough data to return even a Tet10 without a segfault if nglib
    // went nuts
    int ngnodes[11];

    // i+1 since we must be 1-based with these ids too...
    Ng_Volume_Element_Type ngtype =
      Ng_GetVolumeElement(ngmesh, i+1, ngnodes);

    // But really nglib shouldn't go nuts
    libmesh_assert(ngtype == NG_TET);
    libmesh_ignore(ngtype);

    auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
    for (auto n : make_range(4))
      {
        const dof_id_type node_id =
          libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
        elem->set_node(n, this->_mesh.node_ptr(node_id));
      }

    // NetGen and we disagree about node numbering orientation
    elem->orient(bi);

    for (auto s : make_range(4))
      {
        std::array<dof_id_type,3> sorted_ids =
          {DofObject::invalid_id, DofObject::invalid_id,
           DofObject::invalid_id};

        std::vector<unsigned int> nos = elem->nodes_on_side(s);
        for (auto n : nos)
          insert_id(sorted_ids, elem->node_id(n));

        if (auto it = side_boundary_id.find(sorted_ids);
            it != side_boundary_id.end())
          bi->add_side(elem, s, it->second);
      }
  }

  // We don't need our holes anymore.  Delete their serializers
  // first to avoid dereferencing dangling pointers.
  hole_serializers.clear();
  if (_holes)
    _holes->clear();

  // Send our data to other ranks, then fix it up together
  MeshCommunication().broadcast(this->_mesh);
  this->_mesh.prepare_for_use();
}

volume_to_surface_mesh()

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

staticprotectedinherited

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

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

Definition at line 131 of file mesh_tet_interface.C.

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

Referenced by triangulate().

{
  // 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;
}

Member Data Documentation

_desired_volume

Real libMesh::MeshTetInterface::_desired_volume

protectedinherited

The desired volume for the elements in the resulting mesh.

Unlimited (indicated by 0) by default

Definition at line 138 of file mesh_tet_interface.h.

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

_elem_type

ElemType libMesh::MeshTetInterface::_elem_type

protectedinherited

The exact type of tetrahedra we intend to construct.

Definition at line 149 of file mesh_tet_interface.h.

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

_holes

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

protectedinherited

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

If this is nullptr, there are no holes!

Definition at line 160 of file mesh_tet_interface.h.

Referenced by triangulate().

_mesh

UnstructuredMesh& libMesh::MeshTetInterface::_mesh

protectedinherited

Local reference to the mesh we are working with.

Definition at line 154 of file mesh_tet_interface.h.

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

_serializer

MeshSerializer libMesh::NetGenMeshInterface::_serializer

protected

Tetgen only operates on serial meshes.

Definition at line 78 of file mesh_netgen_interface.h.

_smooth_after_generating

bool libMesh::MeshTetInterface::_smooth_after_generating

protectedinherited

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

False by default.

Definition at line 144 of file mesh_tet_interface.h.

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

---

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

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