libMesh::Poly2TriTriangulator Class Reference

A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertion. More...

#include <poly2tri_triangulator.h>

Inheritance diagram for libMesh::Poly2TriTriangulator:

Public Types
enum	TriangulationType { GENERATE_CONVEX_HULL = 0, PSLG = 1, INVALID_TRIANGULATION_TYPE }
	The TriangulationType is used with the general triangulate function defined below. More...

Public Member Functions
	Poly2TriTriangulator (UnstructuredMesh &mesh, dof_id_type n_boundary_nodes=DofObject::invalid_id)
	The constructor. More...
virtual	~Poly2TriTriangulator ()
	Empty destructor. More...
virtual void	triangulate () override
	Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points as necessary to promote mesh quality. More...
virtual void	set_desired_area_function (FunctionBase< Real > *desired) override
	Set a function giving desired triangle area as a function of position. More...
virtual FunctionBase< Real > *	get_desired_area_function () override
	Get the function giving desired triangle area as a function of position, or nullptr if no such function has been set. More...
virtual void	set_refine_boundary_allowed (bool refine_bdy_allowed) override
	Set whether or not the triangulation is allowed to refine the mesh boundary when refining the interior. More...
virtual bool	refine_boundary_allowed () const override
	Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interior. More...
ElemType &	elem_type ()
	Sets and/or gets the desired element type. More...
Real &	desired_area ()
	Sets and/or gets the desired triangle area. More...
Real &	minimum_angle ()
	Sets and/or gets the minimum desired angle. More...
TriangulationType &	triangulation_type ()
	Sets and/or gets the desired triangulation type. More...
bool &	insert_extra_points ()
	Sets and/or gets the flag for inserting add'l points. More...
void	set_interpolate_boundary_points (int n_points)
	Complicated setter, for compatibility with insert_extra_points() More...
int	get_interpolate_boundary_points () const
	Complicated getter, for compatibility with insert_extra_points() More...
bool &	smooth_after_generating ()
	Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the grid. More...
bool &	quiet ()
	Whether not to silence internal messages to stdout. More...
void	attach_hole_list (const std::vector< Hole *> *holes)
	Attaches a vector of Hole* pointers which will be meshed around. More...
void	set_verify_hole_boundaries (bool v)
	Verifying that hole boundaries don't cross the outer boundary or each other is something like O(N_bdys^2*N_points_per_bdy^2), so we only do it if requested. More...
bool	get_verify_hole_boundaries () const
void	attach_boundary_marker (const std::vector< int > *markers)
	Attaches boundary markers. More...
void	attach_region_list (const std::vector< Region *> *regions)
	Attaches regions for using attribute to set subdomain IDs and better controlling the triangle sizes within the regions. More...
void	set_auto_area_function (const Parallel::Communicator &comm, const unsigned int num_nearest_pts, const unsigned int power, const Real background_value, const Real background_eff_dist)
	Generate an auto area function based on spacing of boundary points. More...
bool	has_auto_area_function ()
	Whether or not an auto area function has been set. More...
FunctionBase< Real > *	get_auto_area_function ()
	Get the auto area function. More...
void	calculate_auto_desired_area_samples (std::vector< Point > &function_points, std::vector< Real > &function_sizes, const Real &area_factor=1.5)
	The external boundary and all hole boundaries are collected. More...
void	set_outer_boundary_ids (std::set< std::size_t > bdy_ids)
	A set of ids to allow on the outer boundary loop: interpreted as boundary ids of 2D elements and/or subdomain ids of 1D edges. More...
const std::set< std::size_t > &	get_outer_boundary_ids () const

Public Attributes
std::vector< std::pair< unsigned int, unsigned int > >	segments
	When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly through the ordering of the points, you can use the segments vector to specify the segments explicitly, Ex: unit square numbered counter-clockwise starting from origin segments[0] = (0,1) segments[1] = (1,2) segments[2] = (2,3) segments[3] = (3,0) (For the above case you could actually use the implicit ordering!) More...
std::vector< Point >	segment_midpoints
	When constructing a second-order triangulation from a second-order boundary, we may do the triangulation using first-order elements, in which case we need to save midpoint location data in order to reconstruct curvature along boundaries. More...

Protected Member Functions
bool	is_refine_boundary_allowed (const BoundaryInfo &boundary_info, const Elem &elem, unsigned int side)
	Is refining this element's boundary side allowed? More...
void	triangulate_current_points ()
	Triangulate the current mesh and hole points. More...
bool	insert_refinement_points ()
	Add Steiner points as new mesh nodes, as necessary to refine an existing trangulation. More...
bool	should_refine_elem (Elem &elem)
	Returns true if the given element ought to be refined according to current criteria. More...
void	elems_to_segments ()
	Helper function to create PSLG segments from our other boundary-defining options (1D mesh edges, 2D mesh boundary sides), if no segments already exist. More...
void	nodes_to_segments (dof_id_type max_node_id)
	Helper function to create PSLG segments from our node ordering, up to the maximum node id, if no segments already exist. More...
void	insert_any_extra_boundary_points ()
	Helper function to add extra points (midpoints of initial segments) to a PSLG triangulation. More...
void	increase_triangle_order ()
	Helper function to upconvert Tri3 to any higher order triangle type if requested via _elem_type. More...
void	verify_holes (const Hole &outer_bdy)
	Helper function to check holes for intersections if requested. More...
unsigned int	total_hole_points ()
	Helper function to count points in and verify holes. More...

Protected Attributes
UnstructuredMesh &	_mesh
	Reference to the mesh which is to be created by triangle. More...
const std::vector< Hole * > *	_holes
	A pointer to a vector of Hole*s. More...
const std::vector< int > *	_markers
	Boundary markers. More...
const std::vector< Region * > *	_regions
	A pointer to a vector of Regions*s. More...
std::set< std::size_t >	_bdy_ids
	A set of ids to allow on the outer boundary loop. More...
ElemType	_elem_type
	The type of elements to generate. More...
Real	_desired_area
	The desired area for the elements in the resulting mesh. More...
Real	_minimum_angle
	Minimum angle in triangles. More...
TriangulationType	_triangulation_type
	The type of triangulation to perform: choices are: convex hull PSLG. More...
bool	_insert_extra_points
	Flag which tells whether or not to insert additional nodes before triangulation. More...
int	_interpolate_boundary_points
	Flag which tells how many additional nodes should be inserted between each pair of original mesh points. More...
bool	_smooth_after_generating
	Flag which tells whether we should smooth the mesh after it is generated. More...
bool	_quiet
	Flag which tells if we want to suppress stdout outputs. More...
bool	_verify_hole_boundaries
	Flag which tells if we want to check hole geometry. More...
std::unique_ptr< AutoAreaFunction >	_auto_area_function
	The auto area function based on the spacing of boundary points. More...

Private Attributes
std::map< const Hole *, std::unique_ptr< ArbitraryHole > >	replaced_holes
	We might have to replace the user-provided holes with refined versions. More...
dof_id_type	_n_boundary_nodes
	Keep track of how many mesh nodes are boundary nodes. More...
std::unique_ptr< FunctionBase< Real > >	_desired_area_func
	Location-dependent area requirements. More...
bool	_refine_bdy_allowed
	Whether to allow boundary refinement. More...

Detailed Description

A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertion.

Author

Roy H. Stogner

Date

2022

Definition at line 46 of file poly2tri_triangulator.h.

Member Enumeration Documentation

TriangulationType

enum libMesh::TriangulatorInterface::TriangulationType

inherited

The TriangulationType is used with the general triangulate function defined below.

Enumerator
GENERATE_CONVEX_HULL	First generate a convex hull from the set of points passed in, and then triangulate this set of points. This is probably the most common type of usage.
PSLG	Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of the "points" vector: a straight line is assumed to lie between each successive pair of points, with an additional line joining the final and first points. Explicitly telling the triangulator to add additional points may be important for this option.
INVALID_TRIANGULATION_TYPE	Does nothing, used as a "null" value.

Definition at line 111 of file triangulator_interface.h.

    {
      GENERATE_CONVEX_HULL = 0,

      PSLG = 1,

      INVALID_TRIANGULATION_TYPE
    };

Constructor & Destructor Documentation

Poly2TriTriangulator()

libMesh::Poly2TriTriangulator::Poly2TriTriangulator ( UnstructuredMesh &  mesh, dof_id_type  n_boundary_nodes = DofObject::invalid_id  )

explicit

The constructor.

A reference to the mesh containing the points which are to be triangulated must be provided. The first n_boundary_nodes are expected to form a closed loop around the mesh domain; any subsequent nodes are expected to be interior nodes or in the middle of (internal hole or external) boundary segments.

If n_boundary_nodes is not supplied or is invalid_id then all mesh points are expected to be boundary polyline points.

Definition at line 354 of file poly2tri_triangulator.C.

  : TriangulatorInterface(mesh),
    _n_boundary_nodes(n_boundary_nodes),
    _refine_bdy_allowed(true)
{
}

~Poly2TriTriangulator()

libMesh::Poly2TriTriangulator::~Poly2TriTriangulator ( )

virtualdefault

Empty destructor.

Defaulted in the .C so we can forward declare unique_ptr contents.

Member Function Documentation

attach_boundary_marker()

void libMesh::TriangulatorInterface::attach_boundary_marker ( const std::vector< int > *  markers )

inlineinherited

Attaches boundary markers.

If segments is set, the number of markers must be equal to the size of segments, otherwise, it is equal to the number of points.

Definition at line 294 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_markers.

{ _markers = markers; }

attach_hole_list()

void libMesh::TriangulatorInterface::attach_hole_list ( const std::vector< Hole *> *  holes )

inlineinherited

Attaches a vector of Hole* pointers which will be meshed around.

Definition at line 255 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_holes.

Referenced by libMesh::MeshTools::Generation::build_delaunay_square(), MeshTriangulationTest::testPoly2TriHolesInterpRefined(), MeshTriangulationTest::testPoly2TriRefinementBase(), MeshTriangulationTest::testTriangulatorHoles(), MeshTriangulationTest::testTriangulatorMeshedHoles(), MeshTriangulationTest::testTriangulatorRoundHole(), and triangulate_domain().

{_holes = holes;}

attach_region_list()

void libMesh::TriangulatorInterface::attach_region_list ( const std::vector< Region *> *  regions )

inlineinherited

Attaches regions for using attribute to set subdomain IDs and better controlling the triangle sizes within the regions.

Definition at line 300 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_regions.

{ _regions = regions; }

calculate_auto_desired_area_samples()

void libMesh::TriangulatorInterface::calculate_auto_desired_area_samples ( std::vector< Point > &  function_points, std::vector< Real > &  function_sizes, const Real &  area_factor = 1.5  )

inherited

The external boundary and all hole boundaries are collected.

The centroid of each EDGE element is used as the point position and the desired area is calculated as the area of the equilateral triangle with the edge length as the length of the EDGE element times an area_factor (default is 1.5).

Definition at line 500 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_bdy_ids, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, and libMesh::Real.

Referenced by libMesh::TriangulatorInterface::get_auto_area_function().

{
  // Get the hole mesh of the outer boundary
  // Holes should already be attached if applicable when this function is called
  const TriangulatorInterface::MeshedHole bdry_mh { _mesh, this->_bdy_ids };
  // Collect all the centroid points of the outer boundary segments
  // and the corresponding element sizes
  for (unsigned int i = 0; i < bdry_mh.n_points(); i++)
  {
    function_points.push_back((bdry_mh.point(i) + bdry_mh.point((i + 1) % bdry_mh.n_points())) /
                              Real(2.0));
    function_sizes.push_back(
        (bdry_mh.point(i) - bdry_mh.point((i + 1) % bdry_mh.n_points())).norm());
  }
  // If holes are present, do the same for the hole boundaries
  if(_holes)
    for (const Hole * hole : *_holes)
    {
      for (unsigned int i = 0; i < hole->n_points(); i++)
      {
        function_points.push_back(
            (hole->point(i) + hole->point((i + 1) % hole->n_points())) / Real(2.0));
        function_sizes.push_back(
            (hole->point(i) - hole->point((i + 1) % hole->n_points())).norm());
      }
    }

  std::for_each(
      function_sizes.begin(), function_sizes.end(), [&area_factor](Real & a) { a = a * a * area_factor * std::sqrt(3.0) / 4.0; });

}

desired_area()

Real& libMesh::TriangulatorInterface::desired_area ( )

inlineinherited

Sets and/or gets the desired triangle area.

Set to zero to disable area constraint.

If a desired_area_function is set, then desired_area() should be used to set a minimum desired area; this will reduce "false negatives" by suggesting how finely to sample desired_area_function inside large triangles, where ideally the desired_area_function will be satisfied in the triangle interior and not just at the triangle vertices.

Definition at line 177 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_desired_area.

Referenced by libMesh::MeshTools::Generation::build_delaunay_square(), MeshTriangulationTest::commonSettings(), insert_refinement_points(), main(), should_refine_elem(), MeshTriangulationTest::testPoly2TriRefinementBase(), MeshTriangulationTest::testTriangulatorInterp(), and triangulate_domain().

{return _desired_area;}

elem_type()

ElemType& libMesh::TriangulatorInterface::elem_type ( )

inlineinherited

Sets and/or gets the desired element type.

Definition at line 164 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_elem_type.

Referenced by libMesh::MeshTools::Generation::build_delaunay_square(), and MeshTriangulationTest::testTriangulatorRoundHole().

{return _elem_type;}

elems_to_segments()

void libMesh::TriangulatorInterface::elems_to_segments ( )

protectedinherited

Helper function to create PSLG segments from our other boundary-defining options (1D mesh edges, 2D mesh boundary sides), if no segments already exist.

Definition at line 149 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_bdy_ids, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_verify_hole_boundaries, libMesh::MeshBase::clear_elems(), libMesh::MeshBase::delete_node(), libMesh::make_range(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::node_ptr(), libMesh::TriangulatorInterface::segment_midpoints, libMesh::TriangulatorInterface::segments, and libMesh::TriangulatorInterface::verify_holes().

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

{
  // Don't try to override manually specified segments
  if (!this->segments.empty())
    return;

  // If we have edges, they should form the polyline with the ordering
  // we want.  Let's turn them into segments for later use, because
  // we're going to delete the original elements to replace with our
  // triangulation.
  if (_mesh.n_elem())
    {
      // Mapping from points to node ids, to back those out from
      // MeshedHole results later
      std::map<Point, dof_id_type> point_id_map;

      for (Node * node : _mesh.node_ptr_range())
        {
          // We're not going to support overlapping nodes on the boundary
          libmesh_error_msg_if
            (point_id_map.count(*node),
             "TriangulatorInterface does not support overlapping nodes found at "
             << static_cast<Point&>(*node));

          point_id_map.emplace(*node, node->id());
        }

      // We don't support directly generating Tri6, so for
      // compatibility with future stitching we need to be working
      // with first-order elements.  Let's get rid of any non-vertex
      // nodes we just added.
      for (Elem * elem : _mesh.element_ptr_range())
        for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
          point_id_map.erase(elem->point(n));

      // We'll steal the ordering calculation from
      // the MeshedHole code
      const TriangulatorInterface::MeshedHole mh { _mesh, this->_bdy_ids };

      // If we've specified only a subset of the mesh as our outer
      // boundary, then we may have nodes that don't actually fall
      // inside that boundary.  Triangulator code doesn't like Steiner
      // points that aren't inside the triangulation domain, so we
      // need to get rid of them.
      //
      // Also, if we're using Edge3 elements to define our outer
      // boundary, we're only dealing with their 2 end nodes and we'll
      // need to get rid of their central nodes.
      std::unordered_set<Node *> nodes_to_delete;

      for (Elem * elem : _mesh.element_ptr_range())
        for (auto n : make_range(elem->n_vertices(), elem->n_nodes()))
          nodes_to_delete.insert(elem->node_ptr(n));

      if (!this->_bdy_ids.empty())
        {
          for (auto & node : _mesh.node_ptr_range())
            if (!mh.contains(*node))
              nodes_to_delete.insert(node);
        }

      // And now we're done with elements.  Delete them lest they have
      // dangling pointers to nodes we'll be deleting.
      _mesh.clear_elems();

      // Make segments from boundary nodes; also make sure we don't
      // delete them.
      const std::size_t np = mh.n_points();
      for (auto i : make_range(np))
        {
          const Point pt = mh.point(i);
          const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
          nodes_to_delete.erase(_mesh.node_ptr(id0));
          const Point next_pt = mh.point((np+i+1)%np);
          const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
          this->segments.emplace_back(id0, id1);
          for (auto m : make_range(mh.n_midpoints()))
            this->segment_midpoints.emplace_back(mh.midpoint(m, i));
        }

      for (Node * node : nodes_to_delete)
        _mesh.delete_node(node);

      if (this->_verify_hole_boundaries && _holes)
        this->verify_holes(mh);
    }
}

get_auto_area_function()

FunctionBase< Real > * libMesh::TriangulatorInterface::get_auto_area_function ( )

inherited

Get the auto area function.

Definition at line 487 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_auto_area_function, and libMesh::TriangulatorInterface::calculate_auto_desired_area_samples().

Referenced by should_refine_elem().

{
  if (!_auto_area_function->initialized())
  {
    // Points and target element sizes for the interpolation
    std::vector<Point> function_points;
    std::vector<Real> function_sizes;
    calculate_auto_desired_area_samples(function_points, function_sizes);
    _auto_area_function->init_mfi(function_points, function_sizes);
  }
  return _auto_area_function.get();
}

get_desired_area_function()

FunctionBase< Real > * libMesh::Poly2TriTriangulator::get_desired_area_function ( )

overridevirtual

Get the function giving desired triangle area as a function of position, or nullptr if no such function has been set.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 456 of file poly2tri_triangulator.C.

References _desired_area_func.

Referenced by insert_refinement_points(), and should_refine_elem().

{
  return _desired_area_func.get();
}

get_interpolate_boundary_points()

int libMesh::TriangulatorInterface::get_interpolate_boundary_points ( ) const

inherited

Complicated getter, for compatibility with insert_extra_points()

Definition at line 137 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_insert_extra_points, and libMesh::TriangulatorInterface::_interpolate_boundary_points.

Referenced by libMesh::TriangulatorInterface::insert_any_extra_boundary_points().

{
  // backwards compatibility - someone might have turned this off via
  // the old API
  if (!_insert_extra_points)
    return 0;

  return _interpolate_boundary_points;
}

get_outer_boundary_ids()

const std::set<std::size_t>& libMesh::TriangulatorInterface::get_outer_boundary_ids ( ) const

inlineinherited

Definition at line 352 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_bdy_ids.

{ return _bdy_ids; }

get_verify_hole_boundaries()

bool libMesh::TriangulatorInterface::get_verify_hole_boundaries ( ) const

inlineinherited

Definition at line 264 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_verify_hole_boundaries.

{return _verify_hole_boundaries;}

has_auto_area_function()

bool libMesh::TriangulatorInterface::has_auto_area_function ( )

inlineinherited

Whether or not an auto area function has been set.

Definition at line 329 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_auto_area_function.

Referenced by insert_refinement_points(), and should_refine_elem().

{return _auto_area_function != nullptr;}

increase_triangle_order()

void libMesh::TriangulatorInterface::increase_triangle_order ( )

protectedinherited

Helper function to upconvert Tri3 to any higher order triangle type if requested via _elem_type.

Should be called at the end of triangulate()

Definition at line 332 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::MeshBase::all_complete_order(), libMesh::MeshBase::all_second_order(), libMesh::FIRST, libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::point(), libMesh::TriangulatorInterface::segment_midpoints, libMesh::TriangulatorInterface::segments, libMesh::TRI3, libMesh::TRI6, and libMesh::TRI7.

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

{
  switch (_elem_type)
    {
    case TRI3:
    // Nothing to do if we're not requested to increase order
      return;
    case TRI6:
      _mesh.all_second_order();
      break;
    case TRI7:
      _mesh.all_complete_order();
      break;
    default:
      libmesh_not_implemented();
    }

  // If we have any midpoint location data, we'll want to look it up
  // by point.  all_midpoints[{p, m}] will be the mth midpoint
  // location following after point p (when traversing a triangle
  // counter-clockwise)
  std::map<std::pair<Point, unsigned int>, Point> all_midpoints;
  unsigned int n_midpoints =
    this->segment_midpoints.size() / this->segments.size();
  libmesh_assert_equal_to(this->segments.size() * n_midpoints,
                          this->segment_midpoints.size());
  for (auto m : make_range(n_midpoints))
    for (auto i : make_range(this->segments.size()))
      {
        const Point & p = _mesh.point(this->segments[i].first);
        all_midpoints[{p,m}] =
          this->segment_midpoints[i*n_midpoints+m];
      }

  if (_holes)
    for (const Hole * hole : *_holes)
      {
        if (!hole->n_midpoints())
          continue;
        if (!n_midpoints)
          n_midpoints = hole->n_midpoints();
        else if (hole->n_midpoints() != n_midpoints)
          libmesh_not_implemented_msg
            ("Differing boundary midpoint counts " <<
             hole->n_midpoints() << " and " << n_midpoints);

        // Our inner holes are expected to have points in
        // counter-clockwise order, which is backwards from how we
        // want to traverse them when iterating in counter-clockwise
        // order over a triangle, so we'll need to reverse our maps
        // carefully here.
        const auto n_hole_points = hole->n_points();
        libmesh_assert(n_hole_points);
        for (auto m : make_range(n_midpoints))
          {
            for (auto i : make_range(n_hole_points-1))
              {
                const Point & p = hole->point(i+1);
                all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, i);
              }
            const Point & p = hole->point(0);
            all_midpoints[{p,m}] = hole->midpoint(n_midpoints-m-1, n_hole_points-1);
          }
      }

  // The n_midpoints > 1 case is for future proofing, but in the
  // present we have EDGE4 and no TRI10 yet.
  if (n_midpoints > 1)
    libmesh_not_implemented_msg
      ("Cannot construct triangles with more than 1 midpoint per edge");

  if (!n_midpoints)
    return;

  for (Elem * elem : _mesh.element_ptr_range())
    {
      // This should only be called right after we've finished
      // converting a triangulation to higher order
      libmesh_assert_equal_to(elem->n_vertices(), 3);
      libmesh_assert_not_equal_to(elem->default_order(), FIRST);

      for (auto n : make_range(3))
        {
          // Only hole/outer boundary segments need adjusted midpoints
          if (elem->neighbor_ptr(n))
            continue;

          const Point & p = elem->point(n);

          if (const auto it = all_midpoints.find({p,0});
              it != all_midpoints.end())
            elem->point(n+3) = it->second;
        }
    }
}

insert_any_extra_boundary_points()

void libMesh::TriangulatorInterface::insert_any_extra_boundary_points ( )

protectedinherited

Helper function to add extra points (midpoints of initial segments) to a PSLG triangulation.

Definition at line 278 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_triangulation_type, libMesh::MeshBase::add_point(), libMesh::TriangulatorInterface::get_interpolate_boundary_points(), libMesh::DofObject::id(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), libMesh::MeshBase::node_ptr(), libMesh::TriangulatorInterface::PSLG, and libMesh::TriangulatorInterface::segments.

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

{
  // If the initial PSLG is really simple, e.g. an L-shaped domain or
  // a square/rectangle, the resulting triangulation may be very
  // "structured" looking.  Sometimes this is a problem if your
  // intention is to work with an "unstructured" looking grid.  We can
  // attempt to work around this limitation by inserting midpoints
  // into the original PSLG.  Inserting additional points into a
  // set of points meant to be a convex hull usually makes less sense.

  const int n_interpolated = this->get_interpolate_boundary_points();
  if ((_triangulation_type==PSLG) && n_interpolated)
    {
      // If we were lucky enough to start with contiguous node ids,
      // let's keep them that way.
      dof_id_type nn = _mesh.max_node_id();

      std::vector<std::pair<unsigned int, unsigned int>> old_segments =
        std::move(this->segments);

      // We expect to have converted any elems and/or nodes into
      // segments by now.
      libmesh_assert(!old_segments.empty());

      this->segments.clear();

      // Insert a new point on each segment at evenly spaced locations
      // between existing boundary points.
      // np=index into new points vector
      // n =index into original points vector
      for (auto old_segment : old_segments)
        {
          Node * begin_node = _mesh.node_ptr(old_segment.first);
          Node * end_node = _mesh.node_ptr(old_segment.second);
          dof_id_type current_id = begin_node->id();
          for (auto i : make_range(n_interpolated))
            {
              // new points are equispaced along the original segments
              const Point new_point =
                ((n_interpolated-i) * *(Point *)(begin_node) +
                 (i+1) * *(Point *)(end_node)) /
                (n_interpolated + 1);
              Node * next_node = _mesh.add_point(new_point, nn++);
              this->segments.emplace_back(current_id,
                                          next_node->id());
              current_id = next_node->id();
            }
          this->segments.emplace_back(current_id,
                                      end_node->id());
        }
    }
}

insert_extra_points()

bool& libMesh::TriangulatorInterface::insert_extra_points ( )

inlineinherited

Sets and/or gets the flag for inserting add'l points.

Definition at line 210 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_insert_extra_points.

{return _insert_extra_points;}

insert_refinement_points()

bool libMesh::Poly2TriTriangulator::insert_refinement_points ( )

protected

Add Steiner points as new mesh nodes, as necessary to refine an existing trangulation.

Returns true iff new points were added.

Definition at line 782 of file poly2tri_triangulator.C.

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build_with_id(), libMesh::MeshBase::delete_elem(), libMesh::TriangulatorInterface::desired_area(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::MeshBase::get_boundary_info(), get_desired_area_function(), libMesh::TriangulatorInterface::ArbitraryHole::get_points(), libMesh::TriangulatorInterface::has_auto_area_function(), libMesh::DofObject::id(), libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::invalid_uint, is_refine_boundary_allowed(), libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::libmesh_merge_move(), libMesh::make_range(), libMesh::MeshBase::max_elem_id(), libMesh::MeshBase::max_node_id(), mesh, libMesh::TriangulatorInterface::minimum_angle(), libMesh::Elem::neighbor_ptr(), libMesh::MeshBase::node_ptr(), libMesh::Elem::quasicircumcenter(), libMesh::Real, libMesh::BoundaryInfo::remove(), replaced_holes, libMesh::TriangulatorInterface::segments, libMesh::Elem::set_neighbor(), libMesh::TriangulatorInterface::ArbitraryHole::set_points(), should_refine_elem(), libMesh::TOLERANCE, libMesh::TRI3, libMesh::Elem::type(), libMesh::DofObject::valid_id(), libMesh::Elem::vertex_average(), and libMesh::Elem::which_neighbor_am_i().

Referenced by triangulate().

{
  LOG_SCOPE("insert_refinement_points()", "Poly2TriTriangulator");

  if (this->minimum_angle() != 0)
    libmesh_not_implemented();

  // We need neighbor pointers for ray casting and cavity finding
  UnstructuredMesh & mesh = dynamic_cast<UnstructuredMesh &>(this->_mesh);
  mesh.find_neighbors();

  if (this->desired_area() == 0 &&
      this->get_desired_area_function() == nullptr &&
      !this->has_auto_area_function())
    return false;

  BoundaryInfo & boundary_info = _mesh.get_boundary_info();

  // We won't immediately add these, lest we invalidate iterators on a
  // ReplicatedMesh.  They'll still be in the mesh neighbor topology
  // for the purpose of doing Delaunay cavity stuff, so we need to
  // manage memory here, but there's no point in adding them to the
  // Mesh just to remove them again afterward when we hit up poly2tri.

  // We'll need to be able to remove new elems from new_elems, in
  // cases where a later refinement insertion has a not-yet-added
  // element in its cavity, so we'll use a map here to make searching
  // possible.
  //
  // For parallel consistency, we can't order a container we plan to
  // iterate through based on Elem * or a hash of it.  We'll be doing
  // Delaunay swaps so we can't iterate based on geometry.  These are
  // not-yet-added elements so we can't iterate based on proper
  // element ids ... but we can set a temporary element id to use for
  // the purpose.
  struct cmp {
    bool operator()(Elem * a, Elem * b) const {
      libmesh_assert(a == b || a->id() != b->id());
      return (a->id() < b->id());
    }
  } comp;

  std::map<Elem *, std::unique_ptr<Elem>, decltype(comp)> new_elems(comp);

  // We should already be Delaunay when we get here, otherwise we
  // won't be able to stay Delaunay later.  But we're *not* always
  // Delaunay when we get here?  What the hell, poly2tri?  Fixing this
  // is expensive!
  {
    // restore_delaunay should get to the same Delaunay triangulation up to
    // isomorphism regardless of ordering ... but we actually care
    // about the isomorphisms!  If a triangle's nodes are permuted on
    // one processor vs another that's an issue.  So sort our input
    // carefully.
    std::set<Elem *, decltype(comp)> all_elems
      { mesh.elements_begin(), mesh.elements_end(), comp };

    restore_delaunay(all_elems, boundary_info);

    libmesh_assert_delaunay(mesh, new_elems);
  }

  // Map of which points follow which in the boundary polylines.  If
  // we have to add new boundary points, we'll use this to construct
  // an updated this->segments to retriangulate with.  If we have to
  // add new hole points, we'll use this to insert points into an
  // ArbitraryHole.
  std::unordered_map<Point, Node *> next_boundary_node;

  // In cases where we've been working with contiguous node id ranges;
  // let's keep it that way.
  dof_id_type nn = _mesh.max_node_id();
  dof_id_type ne = _mesh.max_elem_id();

  // We can't handle duplicated nodes.  We shouldn't ever create one,
  // but let's make sure of that.
#ifdef DEBUG
  std::unordered_set<Point> mesh_points;
  for (const Node * node : mesh.node_ptr_range())
    {
      libmesh_assert(!mesh_points.count(*node));
      mesh_points.insert(*node);
    }
#endif

  auto add_point = [&mesh,
#ifdef DEBUG
                    &mesh_points,
#endif
                    &nn](const Point & p)
    {
#ifdef DEBUG
      libmesh_assert(!mesh_points.count(p));
      mesh_points.insert(p);
#endif
      return mesh.add_point(p, nn++);
    };

  for (auto & elem : mesh.element_ptr_range())
    {
      // element_ptr_range skips deleted elements ... right?
      libmesh_assert(elem);
      libmesh_assert(elem->valid_id());

      // We only handle triangles in our triangulation
      libmesh_assert_equal_to(elem->level(), 0u);
      libmesh_assert_equal_to(elem->type(), TRI3);

      // If this triangle is as small as we desire, move along
      if (!should_refine_elem(*elem))
        continue;

      // Otherwise add a Steiner point.  We'd like to add the
      // circumcenter ...
      Point new_pt = elem->quasicircumcenter();

      // And to give it a node;
      Node * new_node = nullptr;

      // But that might be outside our triangle, or even outside the
      // boundary.  We'll find a triangle that should contain our new
      // point
      Elem * cavity_elem = elem; // Start looking at elem anyway

      // We'll refine a boundary later if necessary.
      auto boundary_refine = [this, &next_boundary_node,
                              &cavity_elem, &new_node]
        (unsigned int side)
      {
        libmesh_ignore(this); // Only used in dbg/devel
        libmesh_assert(new_node);
        libmesh_assert(new_node->valid_id());

        Node * old_segment_start = cavity_elem->node_ptr(side),
             * old_segment_end   = cavity_elem->node_ptr((side+1)%3);
        libmesh_assert(old_segment_start);
        libmesh_assert(old_segment_start->valid_id());
        libmesh_assert(old_segment_end);
        libmesh_assert(old_segment_end->valid_id());

        if (auto it = next_boundary_node.find(*old_segment_start);
            it != next_boundary_node.end())
          {
            libmesh_assert(it->second == old_segment_end);
            it->second = new_node;
          }
        else
          {
            // This would be an O(N) sanity check if we already
            // have a segments vector or any holes.  :-P
            libmesh_assert(!this->segments.empty() ||
                           (_holes && !_holes->empty()) ||
                           (old_segment_end->id() ==
                            old_segment_start->id() + 1));
            next_boundary_node[*old_segment_start] = new_node;
          }

        next_boundary_node[*new_node] = old_segment_end;
      };

      // Let's find a triangle containing our new point, or at least
      // containing the end of a ray leading from our current triangle
      // to the new point.
      Point ray_start = elem->vertex_average();

      // What side are we coming from, and what side are we going to?
      unsigned int source_side = invalid_uint;
      unsigned int side = invalid_uint;

      while (!cavity_elem->contains_point(new_pt))
        {
          side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side);

          libmesh_assert_not_equal_to (side, invalid_uint);

          Elem * neigh = cavity_elem->neighbor_ptr(side);
          // If we're on a boundary, stop there.  Refine the boundary
          // if we're allowed, the boundary element otherwise.
          if (!neigh)
            {
              if (this->is_refine_boundary_allowed(boundary_info,
                                                   *cavity_elem,
                                                   side))
                {
                  new_pt = ray_start;
                  new_node = add_point(new_pt);
                  boundary_refine(side);
                }
              else
                {
                  // Should we just add the vertex average of the
                  // boundary element, to minimize the number of
                  // slivers created?
                  //
                  // new_pt = cavity_elem->vertex_average();
                  //
                  // That works for a while, but it
                  // seems to be able to "run away" and leave us with
                  // crazy slivers on boundaries if we push interior
                  // refinement too far while disabling boundary
                  // refinement.
                  //
                  // Let's go back to refining our original problem
                  // element.
                  cavity_elem = elem;
                  new_pt = cavity_elem->vertex_average();
                  new_node = add_point(new_pt);

                  // This was going to be a side refinement but it's
                  // now an internal refinement
                  side = invalid_uint;
                }

              break;
            }

          source_side = neigh->which_neighbor_am_i(cavity_elem);
          cavity_elem = neigh;
          side = invalid_uint;
        }

      // If we're ready to create a new node and we're not on a
      // boundary ... should we be?  We don't want to create any
      // sliver elements or confuse poly2tri or anything.
      if (side == invalid_uint && !new_node)
        {
          unsigned int worst_side = libMesh::invalid_uint;
          Real worst_cos = 1;
          for (auto s : make_range(3u))
            {
              // We never snap to a non-domain-boundary
              if (cavity_elem->neighbor_ptr(s))
                continue;

              Real ax = cavity_elem->point(s)(0) - new_pt(0),
                   ay = cavity_elem->point(s)(1) - new_pt(1),
                   bx = cavity_elem->point((s+1)%3)(0) - new_pt(0),
                   by = cavity_elem->point((s+1)%3)(1) - new_pt(1);
              const Real my_cos = (ax*bx+ay*by) /
                                  std::sqrt(ax*ax+ay*ay) /
                                  std::sqrt(bx*bx+by*by);

              if (my_cos < worst_cos)
                {
                  worst_side = s;
                  worst_cos = my_cos;
                }
            }

          // If we'd create a sliver element on the side, let's just
          // refine the side instead, if we're allowed.
          if (worst_cos < -0.6) // -0.5 is the best we could enforce?
            {
              side = worst_side;

              if (this->is_refine_boundary_allowed(boundary_info,
                                                   *cavity_elem,
                                                   side))
                {
                  // Let's just try bisecting for now
                  new_pt = (cavity_elem->point(side) +
                            cavity_elem->point((side+1)%3)) / 2;
                  new_node = add_point(new_pt);
                  boundary_refine(side);
                }
              else // Do the best we can under these restrictions
                {
                  new_pt = cavity_elem->vertex_average();
                  new_node = add_point(new_pt);

                  // This was going to be a side refinement but it's
                  // now an internal refinement
                  side = invalid_uint;
                }
            }
          else
            new_node = add_point(new_pt);
        }
      else
        libmesh_assert(new_node);

      // Find the Delaunay cavity around the new point.
      std::set<Elem *, decltype(comp)> cavity(comp);

      std::set<Elem *, decltype(comp)> unchecked_cavity ({cavity_elem}, comp);
      while (!unchecked_cavity.empty())
        {
          std::set<Elem *, decltype(comp)> checking_cavity(comp);
          checking_cavity.swap(unchecked_cavity);
          for (Elem * checking_elem : checking_cavity)
            {
              for (auto s : make_range(3u))
                {
                  Elem * neigh = checking_elem->neighbor_ptr(s);
                  if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh))
                    continue;

                  if (in_circumcircle(*neigh, new_pt, TOLERANCE*TOLERANCE))
                    unchecked_cavity.insert(neigh);
                }
            }

          libmesh_merge_move(cavity, checking_cavity);
        }

      // Retriangulate the Delaunay cavity.
      // Each of our cavity triangle edges that are exterior to the
      // cavity will be a source of one new triangle.

      // Set of elements that might need Delaunay swaps
      std::set<Elem *, decltype(comp)> check_delaunay_on(comp);

      // Keep maps for doing neighbor pointer assignment.  Not going
      // to iterate through these so hashing pointers is fine.
      std::unordered_map<Node *, std::pair<Elem *, boundary_id_type>>
        neighbors_CCW, neighbors_CW;

      for (Elem * old_elem : cavity)
        {
          for (auto s : make_range(3u))
            {
              Elem * neigh = old_elem->neighbor_ptr(s);
              if (!neigh || !cavity.count(neigh))
                {
                  Node * node_CW = old_elem->node_ptr(s),
                       * node_CCW = old_elem->node_ptr((s+1)%3);

                  auto set_neighbors =
                    [&neighbors_CW, &neighbors_CCW, &node_CW,
                      &node_CCW, &boundary_info]
                    (Elem * new_neigh, boundary_id_type bcid)
                  {
                    // Set clockwise neighbor and vice-versa if possible
                    if (const auto CW_it = neighbors_CW.find(node_CW);
                        CW_it == neighbors_CW.end())
                      {
                        libmesh_assert(!neighbors_CCW.count(node_CW));
                        neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid);
                      }
                    else
                      {
                        Elem * neigh_CW = CW_it->second.first;
                        if (new_neigh)
                          {
                            new_neigh->set_neighbor(0, neigh_CW);
                            boundary_id_type bcid_CW = CW_it->second.second;
                            if (bcid_CW != BoundaryInfo::invalid_id)
                              boundary_info.add_side(new_neigh, 0, bcid_CW);

                          }
                        if (neigh_CW)
                          {
                            neigh_CW->set_neighbor(2, new_neigh);
                            if (bcid != BoundaryInfo::invalid_id)
                              boundary_info.add_side(neigh_CW, 2, bcid);
                          }
                        neighbors_CW.erase(CW_it);
                      }

                    // Set counter-CW neighbor and vice-versa if possible
                    if (const auto CCW_it = neighbors_CCW.find(node_CCW);
                        CCW_it == neighbors_CCW.end())
                      {
                        libmesh_assert(!neighbors_CW.count(node_CCW));
                        neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid);
                      }
                    else
                      {
                        Elem * neigh_CCW = CCW_it->second.first;
                        if (new_neigh)
                          {
                            boundary_id_type bcid_CCW = CCW_it->second.second;
                            new_neigh->set_neighbor(2, neigh_CCW);
                            if (bcid_CCW != BoundaryInfo::invalid_id)
                              boundary_info.add_side(new_neigh, 2, bcid_CCW);
                          }
                        if (neigh_CCW)
                          {
                            neigh_CCW->set_neighbor(0, new_neigh);
                            if (bcid != BoundaryInfo::invalid_id)
                              boundary_info.add_side(neigh_CCW, 0, bcid);
                          }
                        neighbors_CCW.erase(CCW_it);
                      }
                  };

                  // We aren't going to try to add a sliver element if we
                  // have a new boundary node here.  We do need to
                  // keep track of other elements' neighbors, though.
                  if (old_elem == cavity_elem &&
                      s == side)
                    {
                      std::vector<boundary_id_type> bcids;
                      boundary_info.boundary_ids(old_elem, s, bcids);
                      libmesh_assert_equal_to(bcids.size(), 1);
                      set_neighbors(nullptr, bcids[0]);
                      continue;
                    }

                  auto new_elem = Elem::build_with_id(TRI3, ne++);
                  new_elem->set_node(0, new_node);
                  new_elem->set_node(1, node_CW);
                  new_elem->set_node(2, node_CCW);
                  libmesh_assert(!new_elem->is_flipped());

                  // Set in-and-out-of-cavity neighbor pointers
                  new_elem->set_neighbor(1, neigh);
                  if (neigh)
                    {
                      const unsigned int neigh_s =
                        neigh->which_neighbor_am_i(old_elem);
                      neigh->set_neighbor(neigh_s, new_elem.get());
                    }
                  else
                    {
                      std::vector<boundary_id_type> bcids;
                      boundary_info.boundary_ids(old_elem, s, bcids);
                      boundary_info.add_side(new_elem.get(), 1, bcids);
                    }

                  // Set in-cavity neighbors' neighbor pointers
                  set_neighbors(new_elem.get(), BoundaryInfo::invalid_id);

                  // C++ allows function argument evaluation in any
                  // order, but we need get() to precede move
                  Elem * new_elem_ptr = new_elem.get();
                  new_elems.emplace(new_elem_ptr, std::move(new_elem));

                  check_delaunay_on.insert(new_elem_ptr);
                }
            }

          boundary_info.remove(old_elem);
        }

      // Now that we're done using our cavity elems (including with a
      // cavity.find() that used a comparator that dereferences the
      // elements!) it's safe to delete them.
      for (Elem * old_elem : cavity)
        {
          if (const auto it = new_elems.find(old_elem);
              it == new_elems.end())
            mesh.delete_elem(old_elem);
          else
            new_elems.erase(it);
        }

      // Everybody found their match?
      libmesh_assert(neighbors_CW.empty());
      libmesh_assert(neighbors_CCW.empty());

      // Because we're preserving boundaries here, our naive cavity
      // triangulation might not be a Delaunay triangulation.  Let's
      // check and if necessary fix that; we depend on it when doing
      // future point insertions.
      restore_delaunay(check_delaunay_on, boundary_info);

      // This is too expensive to do on every cavity in devel mode
#ifdef DEBUG
      libmesh_assert_delaunay(mesh, new_elems);
#endif
    }

  // If we added any new boundary nodes, we're going to need to keep
  // track of the changes they made to the outer polyline and/or to
  // any holes.
  if (!next_boundary_node.empty())
    {
      auto checked_emplace = [this](dof_id_type new_first,
                                    dof_id_type new_second)
      {
#ifdef DEBUG
        for (auto [first, second] : this->segments)
          {
            libmesh_assert_not_equal_to(first, new_first);
            libmesh_assert_not_equal_to(second, new_second);
          }
        if (!this->segments.empty())
          libmesh_assert_equal_to(this->segments.back().second, new_first);
#endif
        libmesh_assert_not_equal_to(new_first, new_second);

        this->segments.emplace_back (new_first, new_second);
      };

      // Keep track of the outer polyline
      if (this->segments.empty())
        {
          dof_id_type last_id = DofObject::invalid_id;

          // Custom loop because we increment node_it 1+ times inside
          for (auto node_it = _mesh.nodes_begin(),
               node_end = _mesh.nodes_end();
               node_it != node_end;)
            {
              Node & node = **node_it;
              ++node_it;

              const dof_id_type node_id = node.id();

              // Don't add Steiner points
              if (node_id >= _n_boundary_nodes)
                break;

              // Connect up the previous node, if we didn't already
              // connect it after some newly inserted nodes
              if (!this->segments.empty())
                last_id = this->segments.back().second;

              if (last_id != DofObject::invalid_id &&
                  last_id != node_id)
                checked_emplace(last_id, node_id);

              last_id = node_id;

              // Connect to any newly inserted nodes
              Node * this_node = &node;
              auto it = next_boundary_node.find(*this_node);
              while (it != next_boundary_node.end())
                {
                  libmesh_assert(this_node->valid_id());
                  Node * next_node = it->second;
                  libmesh_assert(next_node->valid_id());

                  if (node_it != node_end &&
                      next_node == *node_it)
                    ++node_it;

                  checked_emplace(this_node->id(), next_node->id());

                  this_node = next_node;
                  if (this_node->id() == this->segments.front().first)
                    break;

                  it = next_boundary_node.find(*this_node);
                }
            }

          // We expect a closed loop here
          if (this->segments.back().second != this->segments.front().first)
            checked_emplace(this->segments.back().second,
                            this->segments.front().first);
        }
      else
        {
          std::vector<std::pair<unsigned int, unsigned int>> old_segments;
          old_segments.swap(this->segments);

          auto old_it  = old_segments.begin();

          const Node * node = _mesh.node_ptr(old_it->first);
          const Node * const first_node = node;

          do
            {
              const dof_id_type node_id = node->id();
              if (const auto it = next_boundary_node.find(*node);
                  it == next_boundary_node.end())
                {
                  while (node_id != old_it->first)
                    {
                      ++old_it;
                      libmesh_assert(old_it != old_segments.end());
                    }
                  node = mesh.node_ptr(old_it->second);
                }
              else
                {
                  node = it->second;
                }

              checked_emplace(node_id, node->id());
            }
          while (node != first_node);
        }

      // Keep track of any holes
      if (this->_holes)
        {
          // Do we have any holes that need to be newly replaced?
          for (const Hole * hole : *this->_holes)
            {
              if (this->replaced_holes.count(hole))
                continue;

              bool hole_point_insertion = false;
                for (auto p : make_range(hole->n_points()))
                  if (next_boundary_node.count(hole->point(p)))
                    {
                      hole_point_insertion = true;
                      break;
                    }
              if (hole_point_insertion)
                this->replaced_holes.emplace
                  (hole, std::make_unique<ArbitraryHole>(*hole));
            }

          // If we have any holes that are being replaced, make sure
          // their replacements are up to date.
          for (const Hole * hole : *this->_holes)
            {
              auto hole_it = replaced_holes.find(hole);
              if (hole_it == replaced_holes.end())
                continue;

              ArbitraryHole & arb = *hole_it->second;

              // We only need to update a replacement that's just had
              // new points inserted
              bool point_inserted = false;
              for (const Point & point : arb.get_points())
                if (next_boundary_node.count(point))
                  {
                    point_inserted = true;
                    break;
                  }

              if (!point_inserted)
                continue;

              // Find all points in the replacement hole
              std::vector<Point> new_points;

              // Our outer polyline is expected to have points in
              // counter-clockwise order, so it proceeds "to the left"
              // from the point of view of rays inside the domain
              // pointing outward, and our next_boundary_node ordering
              // was filled accordingly.
              //
              // Our inner holes are expected to have points in
              // counter-clockwise order, but for holes "to the left
              // as viewed from the hole interior is the *opposite* of
              // "to the left as viewed from the domain interior".  We
              // need to build the updated hole ordering "backwards".

              // We should never see duplicate points when we add one
              // to a hole; if we do then we did something wrong.
              auto push_back_new_point = [&new_points](const Point & p) {
                // O(1) assert in devel
                libmesh_assert(new_points.empty() ||
                               new_points.back() != p);
#ifdef DEBUG
                // O(N) asserts in dbg
                for (auto old_p : new_points)
                  libmesh_assert_not_equal_to(old_p, p);
#endif
                new_points.push_back(p);
              };

              for (auto point_it = arb.get_points().rbegin(),
                   point_end = arb.get_points().rend();
                   point_it != point_end;)
                {
                  Point point = *point_it;
                  ++point_it;

                  if (new_points.empty() ||
                      (point != new_points.back() &&
                       point != new_points.front()))
                    push_back_new_point(point);

                  auto it = next_boundary_node.find(point);
                  while (it != next_boundary_node.end())
                    {
                      point = *it->second;
                      if (point == new_points.front())
                        break;
                      if (point_it != point_end &&
                          point == *point_it)
                        ++point_it;
                      push_back_new_point(point);
                      it = next_boundary_node.find(point);
                    }
                }

              std::reverse(new_points.begin(), new_points.end());

              arb.set_points(std::move(new_points));
            }
        }
    }

  // Okay, *now* we can add the new elements.
  for (auto & [raw_elem, unique_elem] : new_elems)
    {
      libmesh_assert_equal_to(raw_elem, unique_elem.get());
      libmesh_assert(!raw_elem->is_flipped());
      libmesh_ignore(raw_elem); // Old gcc warns "unused variable"
      mesh.add_elem(std::move(unique_elem));
    }

  // Did we add anything?
  return !new_elems.empty();
}

is_refine_boundary_allowed()

bool libMesh::Poly2TriTriangulator::is_refine_boundary_allowed ( const BoundaryInfo &  boundary_info, const Elem &  elem, unsigned int  side  )

protected

Is refining this element's boundary side allowed?

Definition at line 463 of file poly2tri_triangulator.C.

References libMesh::BoundaryInfo::boundary_ids(), libMesh::libmesh_assert(), libMesh::Elem::neighbor_ptr(), and libMesh::TriangulatorInterface::Hole::refine_boundary_allowed().

Referenced by insert_refinement_points().

{
  // We should only be calling this on a boundary side
  libmesh_assert(!elem.neighbor_ptr(side));

  std::vector<boundary_id_type> bcids;
  boundary_info.boundary_ids(&elem, side, bcids);

  // We should have one bcid on every boundary side.
  libmesh_assert_equal_to(bcids.size(), 1);

  if (bcids[0] == 0)
    return this->refine_boundary_allowed();

  // If we're not on an outer boundary side we'd better be on a hole
  // side
  libmesh_assert(this->_holes);

  const boundary_id_type hole_num = bcids[0]-1;
  libmesh_assert_less(hole_num, this->_holes->size());
  const Hole * hole = (*this->_holes)[hole_num];
  return hole->refine_boundary_allowed();
}

minimum_angle()

Real& libMesh::TriangulatorInterface::minimum_angle ( )

inlineinherited

Sets and/or gets the minimum desired angle.

Set to zero to disable angle constraint.

Definition at line 200 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_minimum_angle.

Referenced by MeshTriangulationTest::commonSettings(), insert_refinement_points(), and main().

{return _minimum_angle;}

nodes_to_segments()

void libMesh::TriangulatorInterface::nodes_to_segments ( dof_id_type  max_node_id )

protectedinherited

Helper function to create PSLG segments from our node ordering, up to the maximum node id, if no segments already exist.

Definition at line 239 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, libMesh::TriangulatorInterface::_triangulation_type, libMesh::TriangulatorInterface::_verify_hole_boundaries, libMesh::DofObject::id(), libMesh::MeshBase::point(), libMesh::TriangulatorInterface::PSLG, libMesh::TriangulatorInterface::segments, and libMesh::TriangulatorInterface::verify_holes().

Referenced by libMesh::TriangleInterface::triangulate(), and triangulate().

{
  // Don't try to override manually specified segments, or try to add
  // segments if we're doing a convex hull
  if (!this->segments.empty() || _triangulation_type != PSLG)
    return;

  for (auto node_it = _mesh.nodes_begin(),
       node_end = _mesh.nodes_end();
       node_it != node_end;)
    {
      Node * node = *node_it;

      // If we're out of boundary nodes, the rest are going to be
      // Steiner points or hole points
      if (node->id() >= max_node_id)
        break;

      ++node_it;

      Node * next_node = (node_it == node_end) ?
        *_mesh.nodes_begin() : *node_it;

      this->segments.emplace_back(node->id(), next_node->id());
    }

  if (this->_verify_hole_boundaries && _holes)
    {
      std::vector<Point> outer_pts;
      for (auto segment : this->segments)
        outer_pts.push_back(_mesh.point(segment.first));

      ArbitraryHole ah(outer_pts);
      this->verify_holes(ah);
    }
}

quiet()

bool& libMesh::TriangulatorInterface::quiet ( )

inlineinherited

Whether not to silence internal messages to stdout.

Definition at line 249 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_quiet.

{return _quiet;}

refine_boundary_allowed()

virtual bool libMesh::Poly2TriTriangulator::refine_boundary_allowed ( ) const

inlineoverridevirtual

Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.

True by default.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 104 of file poly2tri_triangulator.h.

References _refine_bdy_allowed.

  { return _refine_bdy_allowed; }

set_auto_area_function()

void libMesh::TriangulatorInterface::set_auto_area_function ( const Parallel::Communicator &  comm, const unsigned int  num_nearest_pts, const unsigned int  power, const Real  background_value, const Real  background_eff_dist  )

inherited

Generate an auto area function based on spacing of boundary points.

The external boundary as well as the hole boundaries are taken into consideration to generate the auto area function based on inverse distance interpolation. For each EDGE element on these boundaries, its centroid (midpoint) is used as the point position and the desired area is calculated as 1.5 times of the area of the equilateral triangle with the edge length as the length of the EDGE element. For a given position, the inverse distance interpolation only considers a number of nearest points (set by num_nearest_pts) to calculate the desired area. The weight of the value at each point is calculated as 1/distance^power.

In addition to these conventional inverse distance interpolation features, a concept of "background value" and "background effective distance" is introduced. The background value belongs to a virtual point located at a constant distance (background effective distance) from the given position. The weight of the value at this virtual point is calculated as 1/background_effective_distance^power. Effectively, the background value is the value when the given position is far away from the boundary points.

Definition at line 478 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_auto_area_function.

{
   _auto_area_function = std::make_unique<AutoAreaFunction>(comm, num_nearest_pts, power, background_value, background_eff_dist);
}

set_desired_area_function()

void libMesh::Poly2TriTriangulator::set_desired_area_function ( FunctionBase< Real > *  desired )

overridevirtual

Set a function giving desired triangle area as a function of position.

Set this to nullptr to disable position-dependent area constraint (falling back on desired_area()).

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 447 of file poly2tri_triangulator.C.

References libMesh::FunctionBase< Output >::clone().

Referenced by MeshTriangulationTest::testPoly2TriRefinementBase().

{
  if (desired)
    _desired_area_func = desired->clone();
  else
    _desired_area_func.reset();
}

set_interpolate_boundary_points()

void libMesh::TriangulatorInterface::set_interpolate_boundary_points ( int  n_points )

inherited

Complicated setter, for compatibility with insert_extra_points()

Definition at line 123 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_insert_extra_points, libMesh::TriangulatorInterface::_interpolate_boundary_points, and libMesh::libmesh_assert().

Referenced by MeshTriangulationTest::testTriangulatorInterp().

{
  // Maybe we'll reserve a meaning for negatives later?
  libmesh_assert(n_points >= 0);

  _interpolate_boundary_points = n_points;

  // backwards compatibility - someone (including us) might want to
  // query this via the old API.
  _insert_extra_points = n_points;
}

set_outer_boundary_ids()

void libMesh::TriangulatorInterface::set_outer_boundary_ids ( std::set< std::size_t >  bdy_ids )

inlineinherited

A set of ids to allow on the outer boundary loop: interpreted as boundary ids of 2D elements and/or subdomain ids of 1D edges.

If this is empty, then the outer boundary may be constructed from boundary edges of any id!

Definition at line 351 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_bdy_ids.

Referenced by MeshTriangulationTest::testHalfDomain().

{ _bdy_ids = std::move(bdy_ids); }

set_refine_boundary_allowed()

virtual void libMesh::Poly2TriTriangulator::set_refine_boundary_allowed ( bool  refine_bdy_allowed )

inlineoverridevirtual

Set whether or not the triangulation is allowed to refine the mesh boundary when refining the interior.

This is true by default, but may be set to false to make the mesh boundary more predictable (and so easier to stitch to other meshes) later.

Reimplemented from libMesh::TriangulatorInterface.

Definition at line 97 of file poly2tri_triangulator.h.

References _refine_bdy_allowed.

Referenced by MeshTriangulationTest::testPoly2TriHolesInterpRefined().

  { _refine_bdy_allowed = refine_bdy_allowed; }

set_verify_hole_boundaries()

void libMesh::TriangulatorInterface::set_verify_hole_boundaries ( bool  v )

inlineinherited

Verifying that hole boundaries don't cross the outer boundary or each other is something like O(N_bdys^2*N_points_per_bdy^2), so we only do it if requested.

Definition at line 262 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_verify_hole_boundaries.

Referenced by MeshTriangulationTest::commonSettings().

{_verify_hole_boundaries = v;}

should_refine_elem()

bool libMesh::Poly2TriTriangulator::should_refine_elem ( Elem &  elem )

protected

Returns true if the given element ought to be refined according to current criteria.

Definition at line 1478 of file poly2tri_triangulator.C.

References libMesh::TriangulatorInterface::desired_area(), libMesh::TriangulatorInterface::get_auto_area_function(), get_desired_area_function(), libMesh::TriangulatorInterface::has_auto_area_function(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::Elem::n_vertices(), libMesh::Elem::point(), libMesh::Real, and libMesh::Elem::volume().

Referenced by insert_refinement_points().

{
  const Real min_area_target = this->desired_area();
  FunctionBase<Real> *area_func = this->has_auto_area_function() ? this->get_auto_area_function() : this->get_desired_area_function();

  // If this isn't a question, why are we here?
  libmesh_assert(min_area_target > 0 ||
                 area_func != nullptr ||
                 this->has_auto_area_function());

  const Real area = elem.volume();

  // If we don't have position-dependent area targets we can make a
  // decision quickly
  if (!area_func && !this->has_auto_area_function())
    return (area > min_area_target);
  else if(area_func && this->has_auto_area_function())
    libmesh_warning("WARNING:  both desired are function and automatic area function are set.  Using automatic area function.");

  // If we do?
  //
  // See if we're meeting the local area target at all the elem
  // vertices first
  for (auto v : make_range(elem.n_vertices()))
    {
      // If we have an auto area function, we'll use it and override other area options
      const Real local_area_target = (*area_func)(elem.point(v));
      libmesh_error_msg_if
        (local_area_target <= 0,
         "Non-positive desired element areas are unachievable");
      if (area > local_area_target)
        return true;
    }

  // If our vertices are happy, it's still possible that our interior
  // isn't.  Are we allowed not to bother checking it?
  if (!min_area_target)
    return false;

  libmesh_not_implemented_msg
    ("Combining a minimum desired_area with an area function isn't yet supported.");
}

smooth_after_generating()

bool& libMesh::TriangulatorInterface::smooth_after_generating ( )

inlineinherited

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

Definition at line 244 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_smooth_after_generating.

Referenced by MeshTriangulationTest::commonSettings(), and triangulate_domain().

{return _smooth_after_generating;}

total_hole_points()

unsigned int libMesh::TriangulatorInterface::total_hole_points ( )

protectedinherited

Helper function to count points in and verify holes.

Definition at line 454 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_holes.

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

{
  // If the holes vector is non-nullptr (and non-empty) we need to determine
  // the number of additional points which the holes will add to the
  // triangulation.
  // Note that the number of points is always equal to the number of segments
  // that form the holes.
  unsigned int n_hole_points = 0;

  if (_holes)
    for (const auto & hole : *_holes)
    {
      n_hole_points += hole->n_points();
      // A hole at least has one enclosure.
      // Points on enclosures are ordered so that we can add segments implicitly.
      // Elements in segment_indices() indicates the starting points of all enclosures.
      // The last element in segment_indices() is the number of total points.
      libmesh_assert_greater(hole->segment_indices().size(), 1);
      libmesh_assert_equal_to(hole->segment_indices().back(), hole->n_points());
    }

  return n_hole_points;
}

triangulate()

void libMesh::Poly2TriTriangulator::triangulate ( )

overridevirtual

Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points as necessary to promote mesh quality.

Implements libMesh::TriangulatorInterface.

Definition at line 367 of file poly2tri_triangulator.C.

References libMesh::TriangulatorInterface::_elem_type, libMesh::TriangulatorInterface::_markers, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::TriangulatorInterface::_regions, libMesh::TriangulatorInterface::_smooth_after_generating, libMesh::TriangulatorInterface::_triangulation_type, libMesh::ParallelObject::comm(), libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), libMesh::MeshBase::prepare_for_use(), libMesh::TriangulatorInterface::PSLG, libMesh::MeshBase::set_mesh_dimension(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and triangulate_current_points().

Referenced by main(), and MeshTriangulationTest::testPoly2TriRefinementBase().

{
  LOG_SCOPE("triangulate()", "Poly2TriTriangulator");

  // We only operate on serialized meshes.  And it's not safe to
  // serialize earlier, because it would then be possible for the user
  // to re-parallelize the mesh in between there and here.
  MeshSerializer serializer(_mesh);

  // We don't yet support every set of Triangulator options in the
  // poly2tri implementation

  // We don't support convex hull triangulation, only triangulation of
  // (implicitly defined, by node ordering) polygons (with holes if
  // requested)
  if (_triangulation_type != PSLG)
    libmesh_not_implemented();

  // We currently don't handle region specifications
  if (_regions)
    libmesh_not_implemented();

  // We won't support quads any time soon, or 1D/3D in this interface
  // ever.
  if (_elem_type != TRI3 &&
      _elem_type != TRI6 &&
      _elem_type != TRI7)
    libmesh_not_implemented();

  // If we have no explicit segments defined, we may get them from
  // mesh elements
  this->elems_to_segments();

  // If we *still* have no explicit segments defined, we get them from
  // the order of nodes.
  this->nodes_to_segments(_n_boundary_nodes);

  // Insert additional new points in between existing boundary points,
  // if that is requested and reasonable
  this->insert_any_extra_boundary_points();

  // Triangulate the points we have, then see if we need to add more;
  // repeat until we don't need to add more.
  //
  // This is currently done redundantly in parallel; make sure no
  // processor quits before the others.
  do
    {
      libmesh_parallel_only(_mesh.comm());
      this->triangulate_current_points();
    }
  while (this->insert_refinement_points());

  libmesh_parallel_only(_mesh.comm());

  // Okay, we really do need to support boundary ids soon, but we
  // don't yet
  if (_markers)
    libmesh_not_implemented();

  _mesh.set_mesh_dimension(2);

  // To the naked eye, a few smoothing iterations usually looks better,
  // so we do this by default unless the user says not to.
  if (this->_smooth_after_generating)
    LaplaceMeshSmoother(_mesh).smooth(2);

  // The user might have requested TRI6 or higher instead of TRI3.  We
  // can do this before prepare_for_use() because all we need for it
  // is find_neighbors(), which we did in insert_refinement_points()
  this->increase_triangle_order();

  // Prepare the mesh for use before returning.  This ensures (among
  // other things) that it is partitioned and therefore users can
  // iterate over local elements, etc.
  _mesh.prepare_for_use();
}

triangulate_current_points()

void libMesh::Poly2TriTriangulator::triangulate_current_points ( )

protected

Triangulate the current mesh and hole points.

Definition at line 490 of file poly2tri_triangulator.C.

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::_mesh, _n_boundary_nodes, libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::Elem::build_with_id(), libMesh::BoundaryInfo::clear(), libMesh::MeshBase::clear_elems(), distance(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::DofObject::invalid_id, libMesh::Elem::is_flipped(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), libMesh::MeshBase::n_nodes(), libMesh::TriangulatorInterface::Hole::n_points(), libMesh::TriangulatorInterface::Hole::point(), libMesh::MeshBase::query_node_ptr(), replaced_holes, libMesh::TriangulatorInterface::segments, libMesh::Elem::set_node(), and libMesh::TRI3.

Referenced by triangulate().

{
  LOG_SCOPE("triangulate_current_points()", "Poly2TriTriangulator");

  // Will the triangulation have holes?
  const std::size_t n_holes = _holes != nullptr ? _holes->size() : 0;

  // Mapping from Poly2Tri points to libMesh nodes, so we can get the
  // connectivity translated back later.
  std::map<const p2t::Point, Node *, P2TPointCompare> point_node_map;

  // Poly2Tri data structures
  // Poly2Tri takes vectors of pointers-to-Point for some reason, but
  // we'll just make those shims to vectors of Point rather than
  // individually/manually heap allocating everything.
  std::vector<p2t::Point> outer_boundary_points;
  std::vector<std::vector<p2t::Point>> inner_hole_points(n_holes);

  dof_id_type nn = _mesh.max_node_id();
  libmesh_error_msg_if
    (!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!");

  // Unless we're using an explicit segments list, we assume node ids
  // are contiguous here.
  if (this->segments.empty())
    libmesh_error_msg_if
      (_mesh.n_nodes() != nn,
       "Poly2TriTriangulator needs contiguous node ids or explicit segments!");

  // And if we have more nodes than outer boundary points, the rest
  // may be interior "Steiner points".  We use a set here so we can
  // cheaply search and erase from it later, when we're identifying
  // hole points.
  std::set<p2t::Point, P2TPointCompare> steiner_points;

  // If we were asked to use all mesh nodes as boundary nodes, now's
  // the time to see how many that is.
  if (_n_boundary_nodes == DofObject::invalid_id)
    {
      _n_boundary_nodes = _mesh.n_nodes();
      libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes),
                              std::distance(_mesh.nodes_begin(),
                                            _mesh.nodes_end()));

    }
  else
    libmesh_assert_less_equal(_n_boundary_nodes,
                              _mesh.n_nodes());

  // Prepare poly2tri points for our nodes, sorted into outer boundary
  // points and interior Steiner points.

  if (this->segments.empty())
    {
      // If we have no segments even after taking elems into account,
      // the nodal id ordering defines our outer polyline ordering
      for (auto & node : _mesh.node_ptr_range())
        {
          const p2t::Point pt = to_p2t(*node);

          // If we're out of boundary nodes, the rest are going to be
          // Steiner points or hole points
          if (node->id() < _n_boundary_nodes)
            outer_boundary_points.push_back(pt);
          else
            steiner_points.insert(pt);

          // We're not going to support overlapping nodes on the boundary
          if (point_node_map.count(pt))
            libmesh_not_implemented();

          point_node_map.emplace(pt, node);
        }
    }
  // If we have explicit segments defined, that's our outer polyline
  // ordering:
  else
    {
      // Let's make sure our segments are in order
      dof_id_type last_id = DofObject::invalid_id;

      // Add nodes from every segment, in order, to the outer polyline
      for (auto [segment_start, segment_end] : this->segments)
        {
          if (last_id != DofObject::invalid_id)
            libmesh_error_msg_if(segment_start != last_id,
                                 "Disconnected triangulator segments");
          last_id = segment_end;

          Node * node = _mesh.query_node_ptr(segment_start);

          libmesh_error_msg_if(!node,
                               "Triangulator segments reference nonexistent node id " <<
                               segment_start);

          outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1)));
          p2t::Point * pt = &outer_boundary_points.back();

          // We're not going to support overlapping nodes on the boundary
          if (point_node_map.count(*pt))
            libmesh_not_implemented_msg
              ("Triangulating overlapping boundary nodes is unsupported");

          point_node_map.emplace(*pt, node);
        }

        libmesh_error_msg_if(last_id != this->segments[0].first,
                             "Non-closed-loop triangulator segments");

      // If we have points that aren't in any segments, those will be
      // Steiner points
      for (auto & node : _mesh.node_ptr_range())
        {
          const p2t::Point pt = to_p2t(*node);
          if (const auto it = point_node_map.find(pt);
              it == point_node_map.end())
            {
              steiner_points.insert(pt);
              point_node_map.emplace(pt, node);
            }
          else
            libmesh_assert_equal_to(it->second, node);
        }
    }

  // If we have any elements from a previous triangulation, we're
  // going to replace them with our own.  If we have any elements that
  // were used to create our segments, we're done creating and we no
  // longer need them.
  _mesh.clear_elems();

  // 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've already got the point_node_map to find
  // nodes, so we can just key on pairs of node ids to identify a side.
  std::unordered_map<std::pair<dof_id_type,dof_id_type>,
                     boundary_id_type, libMesh::hash> side_boundary_id;

  const boundary_id_type outer_bcid = 0;
  const std::size_t n_outer = outer_boundary_points.size();

  for (auto i : make_range(n_outer))
    {
      const Node * node1 =
        libmesh_map_find(point_node_map, outer_boundary_points[i]),
                 * node2 =
        libmesh_map_find(point_node_map, outer_boundary_points[(i+1)%n_outer]);

      side_boundary_id.emplace(std::make_pair(node1->id(),
                                              node2->id()),
                               outer_bcid);
    }

  // Create poly2tri triangulator with our mesh points
  std::vector<p2t::Point *> outer_boundary_pointers(n_outer);
  std::transform(outer_boundary_points.begin(),
                 outer_boundary_points.end(),
                 outer_boundary_pointers.begin(),
                 [](p2t::Point & p) { return &p; });


  // Make sure shims for holes last as long as the CDT does; the
  // poly2tri headers don't make clear whether or not they're hanging
  // on to references to these vectors, and it would be reasonable for
  // them to do so.
  std::vector<std::vector<p2t::Point *>> inner_hole_pointers(n_holes);

  p2t::CDT cdt{outer_boundary_pointers};

  // Add any holes
  for (auto h : make_range(n_holes))
    {
      const Hole * initial_hole = (*_holes)[h];
      auto it = replaced_holes.find(initial_hole);
      const Hole & our_hole =
        (it == replaced_holes.end()) ?
        *initial_hole : *it->second;
      auto & poly2tri_hole = inner_hole_points[h];

      for (auto i : make_range(our_hole.n_points()))
        {
          Point p = our_hole.point(i);
          poly2tri_hole.emplace_back(to_p2t(p));

          const auto & pt = poly2tri_hole.back();

          // This won't be a steiner point.
          steiner_points.erase(pt);

          // If we see a hole point already in the mesh, we'll share
          // that node.  This might be a problem if it's a boundary
          // node, but it might just be the same hole point already
          // added during a previous triangulation refinement step.
          if (point_node_map.count(pt))
            {
              libmesh_assert_equal_to
                (point_node_map[pt],
                 _mesh.query_node_ptr(point_node_map[pt]->id()));
            }
          else
            {
              Node * node = _mesh.add_point(p, nn++);
              point_node_map[pt] = node;
            }
        }

      const boundary_id_type inner_bcid = h+1;
      const std::size_t n_inner = poly2tri_hole.size();

      for (auto i : make_range(n_inner))
        {
          const Node * node1 =
            libmesh_map_find(point_node_map, poly2tri_hole[i]),
                     * node2 =
            libmesh_map_find(point_node_map, poly2tri_hole[(i+1)%n_inner]);

          side_boundary_id.emplace(std::make_pair(node1->id(),
                                                  node2->id()),
                                   inner_bcid);
        }

      auto & poly2tri_ptrs = inner_hole_pointers[h];
      poly2tri_ptrs.resize(n_inner);

      std::transform(poly2tri_hole.begin(),
                     poly2tri_hole.end(),
                     poly2tri_ptrs.begin(),
                     [](p2t::Point & p) { return &p; });

      cdt.AddHole(poly2tri_ptrs);
    }

  // Add any steiner points.  We had them in a set, but post-C++11
  // that won't give us non-const element access (even if we
  // pinky-promise not to change the elements in any way that affects
  // our Comparator), and Poly2Tri wants non-const elements (to store
  // edge data?), so we need to move them here.
  std::vector<p2t::Point> steiner_vector(steiner_points.begin(), steiner_points.end());
  steiner_points.clear();
  for (auto & p : steiner_vector)
    cdt.AddPoint(&p);

  // Triangulate!
  cdt.Triangulate();

  // Get poly2tri triangles, turn them into libMesh triangles
  std::vector<p2t::Triangle *> triangles = cdt.GetTriangles();

  // Do our own numbering, even on DistributedMesh
  dof_id_type next_id = 0;

  BoundaryInfo & boundary_info = _mesh.get_boundary_info();
  boundary_info.clear();

  // Add the triangles to our Mesh data structure.
  for (auto ptri_ptr : triangles)
    {
      p2t::Triangle & ptri = *ptri_ptr;

      // We always use TRI3 here, since that's what we have nodes for;
      // if we need a higher order we can convert at the end.
      auto elem = Elem::build_with_id(TRI3, next_id++);
      for (auto v : make_range(3))
        {
          const p2t::Point & vertex = *ptri.GetPoint(v);

          Node * node = libmesh_map_find(point_node_map, vertex);
          libmesh_assert(node);
          elem->set_node(v, node);
        }

      // We expect a consistent triangle orientation
      libmesh_assert(!elem->is_flipped());

      Elem * added_elem = _mesh.add_elem(std::move(elem));

      for (auto v : make_range(3))
        {
          const Node & node1 = added_elem->node_ref(v),
                     & node2 = added_elem->node_ref((v+1)%3);

          auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id()));
          if (it == side_boundary_id.end())
            it = side_boundary_id.find(std::make_pair(node2.id(), node1.id()));
          if (it != side_boundary_id.end())
            boundary_info.add_side(added_elem, v, it->second);
        }
    }
}

triangulation_type()

TriangulationType& libMesh::TriangulatorInterface::triangulation_type ( )

inlineinherited

Sets and/or gets the desired triangulation type.

Definition at line 205 of file triangulator_interface.h.

References libMesh::TriangulatorInterface::_triangulation_type.

Referenced by libMesh::MeshTools::Generation::build_delaunay_square(), MeshTriangulationTest::commonSettings(), main(), and triangulate_domain().

{return _triangulation_type;}

verify_holes()

void libMesh::TriangulatorInterface::verify_holes ( const Hole &  outer_bdy )

protectedinherited

Helper function to check holes for intersections if requested.

Definition at line 429 of file triangulator_interface.C.

References libMesh::TriangulatorInterface::_holes, libMesh::TriangulatorInterface::Hole::contains(), and libMesh::make_range().

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), and libMesh::TriangulatorInterface::nodes_to_segments().

{
  for (const Hole * hole : *_holes)
    {
      for (const Hole * hole2 : *_holes)
        {
          if (hole == hole2)
            continue;

          for (auto i : make_range(hole2->n_points()))
            if (hole->contains(hole2->point(i)))
              libmesh_error_msg
                ("Found point " << hole2->point(i) <<
                 " on one hole boundary and another's interior");
        }

      for (auto i : make_range(hole->n_points()))
        if (!outer_bdy.contains(hole->point(i)))
          libmesh_error_msg
            ("Found point " << hole->point(i) <<
             " on hole boundary but outside outer boundary");
    }
}

Member Data Documentation

_auto_area_function

std::unique_ptr<AutoAreaFunction> libMesh::TriangulatorInterface::_auto_area_function

protectedinherited

The auto area function based on the spacing of boundary points.

Definition at line 476 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_auto_area_function(), libMesh::TriangulatorInterface::has_auto_area_function(), and libMesh::TriangulatorInterface::set_auto_area_function().

_bdy_ids

std::set<std::size_t> libMesh::TriangulatorInterface::_bdy_ids

protectedinherited

A set of ids to allow on the outer boundary loop.

Definition at line 416 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::calculate_auto_desired_area_samples(), libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::get_outer_boundary_ids(), and libMesh::TriangulatorInterface::set_outer_boundary_ids().

_desired_area

Real libMesh::TriangulatorInterface::_desired_area

protectedinherited

The desired area for the elements in the resulting mesh.

Definition at line 427 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::desired_area(), and libMesh::TriangleInterface::triangulate().

_desired_area_func

std::unique_ptr<FunctionBase<Real> > libMesh::Poly2TriTriangulator::_desired_area_func

private

Location-dependent area requirements.

Definition at line 149 of file poly2tri_triangulator.h.

Referenced by get_desired_area_function().

_elem_type

ElemType libMesh::TriangulatorInterface::_elem_type

protectedinherited

The type of elements to generate.

(Defaults to TRI3).

Definition at line 422 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elem_type(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangleInterface::triangulate(), and triangulate().

_holes

const std::vector<Hole*>* libMesh::TriangulatorInterface::_holes

protectedinherited

A pointer to a vector of Hole*s.

If this is nullptr, there are no holes!

Definition at line 400 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::attach_hole_list(), libMesh::TriangulatorInterface::calculate_auto_desired_area_samples(), libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), libMesh::TriangulatorInterface::total_hole_points(), libMesh::TriangleInterface::triangulate(), triangulate_current_points(), and libMesh::TriangulatorInterface::verify_holes().

_insert_extra_points

bool libMesh::TriangulatorInterface::_insert_extra_points

protectedinherited

Flag which tells whether or not to insert additional nodes before triangulation.

This can sometimes be used to "de-regularize" the resulting triangulation.

This flag is supported for backwards compatibility; setting _interpolate_boundary_points = 1 is equivalent.

Definition at line 449 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_interpolate_boundary_points(), libMesh::TriangulatorInterface::insert_extra_points(), and libMesh::TriangulatorInterface::set_interpolate_boundary_points().

_interpolate_boundary_points

int libMesh::TriangulatorInterface::_interpolate_boundary_points

protectedinherited

Flag which tells how many additional nodes should be inserted between each pair of original mesh points.

Definition at line 455 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::get_interpolate_boundary_points(), and libMesh::TriangulatorInterface::set_interpolate_boundary_points().

_markers

const std::vector<int>* libMesh::TriangulatorInterface::_markers

protectedinherited

Boundary markers.

Definition at line 405 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::attach_boundary_marker(), libMesh::TriangleInterface::triangulate(), and triangulate().

_mesh

UnstructuredMesh& libMesh::TriangulatorInterface::_mesh

protectedinherited

Reference to the mesh which is to be created by triangle.

Definition at line 394 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::calculate_auto_desired_area_samples(), libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), libMesh::TriangleInterface::triangulate(), triangulate(), and triangulate_current_points().

_minimum_angle

Real libMesh::TriangulatorInterface::_minimum_angle

protectedinherited

Minimum angle in triangles.

Definition at line 432 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::minimum_angle(), and libMesh::TriangleInterface::triangulate().

_n_boundary_nodes

dof_id_type libMesh::Poly2TriTriangulator::_n_boundary_nodes

private

Keep track of how many mesh nodes are boundary nodes.

Definition at line 144 of file poly2tri_triangulator.h.

Referenced by insert_refinement_points(), triangulate(), and triangulate_current_points().

_quiet

bool libMesh::TriangulatorInterface::_quiet

protectedinherited

Flag which tells if we want to suppress stdout outputs.

Definition at line 466 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::quiet(), and libMesh::TriangleInterface::triangulate().

_refine_bdy_allowed

bool libMesh::Poly2TriTriangulator::_refine_bdy_allowed

private

Whether to allow boundary refinement.

Definition at line 154 of file poly2tri_triangulator.h.

Referenced by refine_boundary_allowed(), and set_refine_boundary_allowed().

_regions

const std::vector<Region*>* libMesh::TriangulatorInterface::_regions

protectedinherited

A pointer to a vector of Regions*s.

If this is nullptr, there are no regions!

Definition at line 411 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::attach_region_list(), libMesh::TriangleInterface::triangulate(), and triangulate().

_smooth_after_generating

bool libMesh::TriangulatorInterface::_smooth_after_generating

protectedinherited

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

True by default.

Definition at line 461 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::smooth_after_generating(), libMesh::TriangleInterface::triangulate(), and triangulate().

_triangulation_type

TriangulationType libMesh::TriangulatorInterface::_triangulation_type

protectedinherited

The type of triangulation to perform: choices are: convex hull PSLG.

Definition at line 439 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), libMesh::TriangulatorInterface::nodes_to_segments(), libMesh::TriangleInterface::triangulate(), triangulate(), and libMesh::TriangulatorInterface::triangulation_type().

_verify_hole_boundaries

bool libMesh::TriangulatorInterface::_verify_hole_boundaries

protectedinherited

Flag which tells if we want to check hole geometry.

Definition at line 471 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::get_verify_hole_boundaries(), libMesh::TriangulatorInterface::nodes_to_segments(), and libMesh::TriangulatorInterface::set_verify_hole_boundaries().

replaced_holes

std::map<const Hole *, std::unique_ptr<ArbitraryHole> > libMesh::Poly2TriTriangulator::replaced_holes

private

We might have to replace the user-provided holes with refined versions.

Definition at line 139 of file poly2tri_triangulator.h.

Referenced by insert_refinement_points(), and triangulate_current_points().

segment_midpoints

std::vector<Point> libMesh::TriangulatorInterface::segment_midpoints

inherited

When constructing a second-order triangulation from a second-order boundary, we may do the triangulation using first-order elements, in which case we need to save midpoint location data in order to reconstruct curvature along boundaries.

Definition at line 287 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), and libMesh::TriangulatorInterface::increase_triangle_order().

segments

std::vector<std::pair<unsigned int, unsigned int> > libMesh::TriangulatorInterface::segments

inherited

When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly through the ordering of the points, you can use the segments vector to specify the segments explicitly, Ex: unit square numbered counter-clockwise starting from origin segments[0] = (0,1) segments[1] = (1,2) segments[2] = (2,3) segments[3] = (3,0) (For the above case you could actually use the implicit ordering!)

Definition at line 279 of file triangulator_interface.h.

Referenced by libMesh::TriangulatorInterface::elems_to_segments(), libMesh::TriangulatorInterface::increase_triangle_order(), libMesh::TriangulatorInterface::insert_any_extra_boundary_points(), insert_refinement_points(), libMesh::TriangulatorInterface::nodes_to_segments(), MeshTriangulationTest::testTriangulatorSegments(), libMesh::TriangleInterface::triangulate(), and triangulate_current_points().

---

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

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