mesh_triangle_holes.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
#include "libmesh/libmesh_config.h"
 
// Local includes
#include "libmesh/mesh_triangle_holes.h"
 
#include "libmesh/boundary_info.h"
#include "libmesh/elem.h"
#include "libmesh/int_range.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_serializer.h"
#include "libmesh/node.h"
#include "libmesh/parallel_algebra.h" // Packing<Point>
#include "libmesh/simple_range.h"
 
// TIMPI includes
#include "timpi/parallel_implementation.h" // broadcast
 
// C++ includes
#include <algorithm>
 
namespace
{
  using namespace libMesh;
 
  int signof(Real val) {
    return (0 < val) - (val < 0);
  }
 
  // Return 1 iff counter-clockwise turn
  // Return -1 iff clockwise turn
  // Return 0 iff collinear
  int orientation(const Point & p0,
                  const Point & p1,
                  const Point & p2)
  {
    const Real detleft  = (p0(0)-p2(0))*(p1(1)-p2(1));
    const Real detright = (p0(1)-p2(1))*(p1(0)-p2(0));
 
    return signof(detleft - detright);
  }
 
  // Same, but for the ray target as it goes to infinity
  int ray_orientation(const Point & p0,
                      const Point & p1,
                      const Point & source,
                      const Point & ray_target)
  {
    const Point rayvec = ray_target - source;
    const Point edgevec = p1 - p0;
    const Real det = edgevec(0)*rayvec(1)-edgevec(1)*rayvec(0);
 
    return signof(det);
  }
 
  // Return true iff the ray from source toward ray_target crosses the
  // edge from edge_pt0 (non-inclusive) to edge_pt1 (inclusive) in the
  // forward ray direction.  edge_pt2 is the vertex following edge_pt1
  // along the polygon; it is used to classify a crossing that lands
  // exactly on the shared vertex edge_pt1.
  //
  // The decision is made purely from orientation predicates on the
  // vertices, never from a parametric edge coordinate compared against
  // a tolerance.  This is what guarantees consistency: a shared vertex
  // yields the exact same orientation when viewed from either of its
  // two edges, so a ray passing through a vertex is counted exactly
  // once for a "through" crossing and not at all for a "tangent"
  // (glancing) touch, regardless of floating-point round-off.
  bool is_intersection(const Point & source,
                       const Point & ray_target,
                       const Point & edge_pt0,
                       const Point & edge_pt1,
                       const Point & edge_pt2)
  {
    const int orient_st0 = orientation(source, ray_target, edge_pt0);
    const int orient_st1 = orientation(source, ray_target, edge_pt1);
 
    // Which side of the edge line the source lies on, and which way the
    // ray turns off the edge, together tell us whether the crossing is
    // ahead of the source rather than behind it.
    const int orient_edge_s = orientation(edge_pt0, edge_pt1, source);
    const int orient_edge_t = ray_orientation(edge_pt0, edge_pt1, source, ray_target);
    const bool forward = (orient_edge_s != orient_edge_t);
 
    // The ray line passes exactly through edge_pt1, the inclusive
    // endpoint of this edge.  Count it only as a genuine "through"
    // crossing: the previous vertex (edge_pt0) and the next vertex
    // (edge_pt2) must lie on opposite sides of the ray line.  A tangent
    // touch (both neighbors on the same side) or a collinear edge
    // (edge_pt0 also on the line) does not count.
    if (orient_st1 == 0)
      {
        const int orient_st2 = orientation(source, ray_target, edge_pt2);
        return forward && (orient_st0 != 0) && (orient_st0 == -orient_st2);
      }
 
    // The ray line passes through edge_pt0, the non-inclusive endpoint;
    // that vertex is counted (if at all) by the previous edge.
    //
    // A run of two or more vertices all lying exactly on the ray line
    // is therefore never counted, even when the boundary arrives at the
    // run from one side and leaves from the other, which ought to count
    // once.  Deciding that needs a look further along the polygon than
    // edge_pt2.  This is long-standing behavior, not something the
    // orientation predicates changed, and it takes exact collinearity
    // to reach.
    if (orient_st0 == 0)
      return false;
 
    // Ordinary case: count it iff the edge endpoints straddle the ray
    // line.
    return forward && (orient_st0 == -orient_st1);
  }
 
  // Returns a non-negative distance iff the ray from source in the
  // direction of ray_target intersects the edge from pt0
  // (non-inclusive) to pt1 (inclusive), -1 otherwise.
  //
  // If the intersection is a "glancing" one at a corner, return -1.
  //
  // edge_pt2 is the vertex following edge_pt1, needed to tell a
  // "glancing" corner from one the ray passes through.
  Real find_intersection(const Point & source,
                         const Point & ray_target,
                         const Point & edge_pt0,
                         const Point & edge_pt1,
                         const Point & edge_pt2)
  {
    // Decide whether the ray crosses this edge using only orientation
    // predicates on the edge vertices, so that shared vertices are
    // classified consistently between neighboring edges.
    if (!is_intersection(source, ray_target, edge_pt0, edge_pt1, edge_pt2))
      return -1;
 
    // Now that the crossing decision is settled, use parametric math
    // only to recover the distance to the intersection.
    const Real raydx = ray_target(0)-source(0),
               raydy = ray_target(1)-source(1),
               edgedx = edge_pt1(0)-edge_pt0(0),
               edgedy = edge_pt1(1)-edge_pt0(1);
    const Real denom = edgedx * raydy - edgedy * raydx;
 
    // divide-by-zero means the segments are parallel; is_intersection
    // rejects that case, but guard against it anyway.
    if (denom == 0)
      return -1;
 
    const Real one_over_denom = 1 / denom;
 
    const Real targetsdx = edge_pt1(0)-ray_target(0),
               targetsdy = edge_pt1(1)-ray_target(1);
 
    const Real u_num = targetsdx * edgedy - targetsdy * edgedx;
    const Real u = u_num * one_over_denom;
    const Real ray_fraction = (1-u);
 
    // is_intersection guarantees a forward crossing, so clamp away any
    // floating-point noise that would otherwise make a vertex hit read
    // as a tiny negative distance.
    const Real distance =
      ray_fraction * std::sqrt(raydx*raydx + raydy*raydy);
    return std::max(Real(0), distance);
  }
}
 
 
namespace libMesh
{
 
//
// Hole member functions
//
Real TriangulatorInterface::Hole::area() const
{
  return this->areavec().norm() / 2;
}
 
 
RealGradient TriangulatorInterface::Hole::areavec() const
{
  const unsigned int np = this->n_points();
 
  if (np < 3)
    return 0;
 
  const Point p0 = this->point(0);
 
  // Every segment (p_{i-1},p_i) from i=2 on defines a triangle w.r.t.
  // p_0.  Add up the cross products of those triangles.  We'll save
  // the division by 2 and the norm for the end.
  //
  // Your hole points had best be coplanar, but this should work
  // regardless of which plane they're in.  If you're in the XY plane,
  // then the standard counter-clockwise hole point ordering gives you
  // a positive areavec(2);
 
  RealGradient areavec = 0;
 
  for (unsigned int i=2; i != np; ++i)
    {
      const Point e_0im = this->point(i-1) - p0,
                  e_0i  = this->point(i) - p0;
 
      areavec += e_0i.cross(e_0im);
    }
 
  return areavec;
}
 
 
 
std::vector<Real>
TriangulatorInterface::Hole::find_ray_intersections(Point ray_start,
                                                    Point ray_target) const
{
  const auto np = this->n_points();
 
  std::vector<Real> intersection_distances;
 
  for (auto i : make_range(np))
    {
      const Point & p0 = this->point(i),
                  & p1 = this->point((i+1)%np),
                  & p2 = this->point((i+2)%np);
      const Real intersection_distance =
        find_intersection(ray_start, ray_target, p0, p1, p2);
      if (intersection_distance >= 0)
        intersection_distances.push_back
          (intersection_distance);
    }
 
  return intersection_distances;
}
 
 
 
Point TriangulatorInterface::Hole::calculate_inside_point() const
{
  // Start with the vertex average
 
  // Turns out "I'm a fully compliant C++17 compiler!" doesn't
  // mean "I have a full C++17 standard library!"
  // inside = std::reduce(points.begin(), points.end());
  Point inside = 0;
  for (auto i : make_range(this->n_points()))
    inside += this->point(i);
 
  inside /= this->n_points();
 
  // Count the number of intersections with a ray to the right,
  // keep track of how far they are
  Point ray_target = inside + Point(1);
  std::vector<Real> intersection_distances =
    this->find_ray_intersections(inside, ray_target);
 
  // The vertex average isn't on the interior, and we found no
  // intersections to the right?  Try looking to the left.
  if (!intersection_distances.size())
    {
      ray_target = inside - Point(1);
      intersection_distances =
        this->find_ray_intersections(inside, ray_target);
    }
 
  // I'd make this an assert, but I'm not 100% confident we can't
  // get here via some kind of FP error on a weird hole shape.
  libmesh_error_msg_if
    (!intersection_distances.size(),
     "Can't find a center for a MeshedHole!");
 
  if (intersection_distances.size() % 2)
    return inside;
 
  // The vertex average is outside.  So go from the vertex average to
  // the closest edge intersection, then halfway to the next-closest.
 
  // Find the nearest first.
  Real min_distance    = std::numeric_limits<Real>::max(),
       second_distance = std::numeric_limits<Real>::max();
  for (Real d : intersection_distances)
    if (d < min_distance)
      {
        second_distance = min_distance;
        min_distance = d;
      }
 
  const Point ray = ray_target - inside;
  inside += ray * (min_distance + second_distance)/2;
 
  return inside;
}
 
 
bool TriangulatorInterface::Hole::contains(Point p) const
{
  // Count the number of intersections with a ray to the right,
  // keep track of how far they are
  Point ray_target = p + Point(1);
  std::vector<Real> intersection_distances =
    this->find_ray_intersections(p, ray_target);
 
  // Odd number of intersections == we're inside
  // Even number == we're outside
  return intersection_distances.size() % 2;
}
 
 
 
//
// PolygonHole member functions
//
TriangulatorInterface::PolygonHole::PolygonHole(const Point & center,
                                                Real radius,
                                                unsigned int n_points_in) :
  _center(center),
  _radius(radius),
  _n_points(n_points_in)
{}
 
 
unsigned int TriangulatorInterface::PolygonHole::n_points() const
{
  return _n_points;
}
 
 
Point TriangulatorInterface::PolygonHole::point(const unsigned int n) const
{
  // The nth point lies at the angle theta = 2 * pi * n / _n_points
  const Real theta = static_cast<Real>(n) * 2.0 * libMesh::pi / static_cast<Real>(_n_points);
 
  return Point(_center(0) + _radius*std::cos(theta), // x=r*cos(theta)
               _center(1) + _radius*std::sin(theta), // y=r*sin(theta)
               0.);
}
 
 
Point TriangulatorInterface::PolygonHole::inside() const
{
  // The center of the hole is definitely inside.
  return _center;
}
 
 
//
// AffineHole member functions
//
Point TriangulatorInterface::AffineHole::point(const unsigned int n) const
{
  return this->transform(_underlying.point(n));
}
 
 
Point TriangulatorInterface::AffineHole::inside() const
{
  return this->transform(_underlying.inside());
}
 
 
Point TriangulatorInterface::AffineHole::transform(const Point & p) const
{
  const Real cos_a = std::cos(_angle);
  const Real sin_a = std::sin(_angle);
  return Point(p(0)*cos_a-p(1)*sin_a + _shift(0),
               p(1)*cos_a+p(1)*sin_a + _shift(1));
}
 
 
//
// ArbitraryHole member functions
//
TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Point & center,
                                                    std::vector<Point> points)
  : _center(center),
    _points(std::move(points))
{
  _segment_indices.push_back(0);
  _segment_indices.push_back(_points.size());
}
 
 
TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Point & center,
                                                    std::vector<Point> points,
                                                    std::vector<unsigned int> segment_indices)
  : _center(center),
    _points(std::move(points)),
    _segment_indices(std::move(segment_indices))
{}
 
 
TriangulatorInterface::ArbitraryHole::ArbitraryHole(std::vector<Point> points)
  : _points(std::move(points))
{
  _segment_indices.push_back(0);
  _segment_indices.push_back(_points.size());
  _center = this->calculate_inside_point();
}
 
 
TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Hole & orig)
  : _center(orig.inside())
{
  const unsigned int np = orig.n_points();
  _points.reserve(np);
  for (auto i : make_range(np))
    _points.push_back(orig.point(i));
}
 
 
unsigned int TriangulatorInterface::ArbitraryHole::n_points() const
{
  return _points.size();
}
 
 
Point TriangulatorInterface::ArbitraryHole::point(const unsigned int n) const
{
  libmesh_assert_less (n, _points.size());
  return _points[n];
}
 
 
Point TriangulatorInterface::ArbitraryHole::inside() const
{
  return _center;
}
 
 
std::vector<unsigned int> TriangulatorInterface::ArbitraryHole::segment_indices() const
{
  return _segment_indices;
}
 
 
//
// MeshedHole member functions
//
 
 
TriangulatorInterface::MeshedHole::MeshedHole(const MeshBase & mesh,
                                              std::set<std::size_t> ids)
  : _center(std::numeric_limits<Real>::max())
{
  // We'll want to do this on one processor and broadcast to the rest;
  // otherwise we can get out of sync by doing things like using
  // pointers as keys.
  libmesh_parallel_only(mesh.comm());
 
  MeshSerializer serial(const_cast<MeshBase &>(mesh),
                        /* serial */ true, /* only proc 0 */ true);
 
  // Try to keep in sync even if we throw an error on proc 0, so we
  // can examine errors in our unit tests in parallel too.
  std::string error_reported;
 
  auto report_error = [&mesh, &error_reported](std::string er) {
    error_reported = std::move(er);
    mesh.comm().broadcast(error_reported);
    libmesh_error_msg(error_reported);
  };
 
  if (mesh.processor_id() != 0)
    {
      // Make sure proc 0 didn't just fail
      mesh.comm().broadcast(error_reported);
      libmesh_error_msg_if(!error_reported.empty(), error_reported);
 
      // Receive the points proc 0 will send later
      mesh.comm().broadcast(_points);
      mesh.comm().broadcast(_midpoints);
      return;
    }
 
  // We'll find all the line segments first, then stitch them together
  // afterward.  If the line segments come from 2D element sides then
  // we'll label their edge_type as "1" for clockwise orientation
  // around the element or "2" for CCW, to make it easier to detect
  // and scream about cases where we have a disconnected outer
  // boundary.
  std::multimap<const Node *,
                std::pair<const Node *, int>> hole_edge_map;
 
  // If we're looking at higher-order elements, we have mid-edge edge
  // nodes to worry about.  hole_midpoint_map[{m,n}][i] should give us
  // the ith mid-edge node traveling from vertex m to vertex n
  std::map<std::pair<const Node *, const Node *>,
           std::vector<const Node *>> hole_midpoint_map;
 
  std::vector<boundary_id_type> bcids;
 
  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
  for (const auto & elem : mesh.active_element_ptr_range())
    {
      if (elem->dim() == 1)
        {
          if (ids.empty() || ids.count(elem->subdomain_id()))
            {
              hole_edge_map.emplace(elem->node_ptr(0),
                                    std::make_pair(elem->node_ptr(1),
                                                   /*edge*/ 0));
              hole_edge_map.emplace(elem->node_ptr(1),
                                    std::make_pair(elem->node_ptr(0),
                                                   /*edge*/ 0));
              if (elem->type() == EDGE3)
                {
                  hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(0),
                                                           elem->node_ptr(1)),
                                            std::vector<const Node *>{elem->node_ptr(2)});
                  hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(1),
                                                           elem->node_ptr(0)),
                                            std::vector<const Node *>{elem->node_ptr(2)});
                }
              else if (elem->type() == EDGE4)
                {
                  hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(0),
                                                           elem->node_ptr(1)),
                                            std::vector<const Node *>{elem->node_ptr(2),
                                                                      elem->node_ptr(3)});
                  hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(1),
                                                           elem->node_ptr(0)),
                                            std::vector<const Node *>{elem->node_ptr(3),
                                                                      elem->node_ptr(2)});
                }
              else
                libmesh_assert_equal_to(elem->default_side_order(), 1);
            }
          continue;
        }
 
      if (elem->dim() == 2)
        {
          const auto ns = elem->n_sides();
          for (auto s : make_range(ns))
            {
              boundary_info.boundary_ids(elem, s, bcids);
 
              bool add_edge = false;
              if (!elem->neighbor_ptr(s) && ids.empty())
                add_edge = true;
 
              if (!add_edge)
                for (auto b : bcids)
                  if (ids.count(b))
                    add_edge = true;
 
              if (add_edge)
                {
                  hole_edge_map.emplace(elem->node_ptr(s),
                                        std::make_pair(elem->node_ptr((s+1)%ns),
                                                       /*counter-CW*/ 2));
                  // Do we really need to support flipped 2D elements?
                  hole_edge_map.emplace(elem->node_ptr((s+1)%ns),
                                        std::make_pair(elem->node_ptr(s),
                                                       /*clockwise*/ 1));
 
                  if (elem->default_side_order() == 2)
                    {
                      hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(s),
                                                               elem->node_ptr((s+1)%ns)),
                                                std::vector<const Node *>{elem->node_ptr(s+ns)});
                      hole_midpoint_map.emplace(std::make_pair(elem->node_ptr((s+1)%ns),
                                                               elem->node_ptr(s)),
                                                std::vector<const Node *>{elem->node_ptr(s+ns)});
                    }
                  else
                    libmesh_assert_equal_to(elem->default_side_order(), 1);
 
                  continue;
                }
            }
        }
    }
 
  if (hole_edge_map.empty())
    report_error("No valid hole edges found in mesh!");
 
  // Function to pull a vector of points out of the map; a loop of
  // edges connecting these points defines a hole boundary.  If the
  // mesh has multiple boundaries (e.g. because it had holes itself),
  // then a random vector will be extracted; this function will be
  // called multiple times so that the various options can be
  // compared.  We choose the largest option.
  auto extract_edge_vector =
    [&report_error, &hole_edge_map, &hole_midpoint_map]() {
    std::tuple<std::vector<const Node *>, std::vector<const Node *>, int>
    hole_points_and_edge_type
      {{hole_edge_map.begin()->first, hole_edge_map.begin()->second.first},
       {}, hole_edge_map.begin()->second.second};
 
    auto & hole_points = std::get<0>(hole_points_and_edge_type);
    auto & midpoint_points = std::get<1>(hole_points_and_edge_type);
    int & edge_type = std::get<2>(hole_points_and_edge_type);
 
    // We won't be needing to search for this edge
    hole_edge_map.erase(hole_points.front());
 
    // Sort the remaining edges into a connected order
    for (const Node * last = hole_points.front(),
                    *    n = hole_points.back();
                         n != hole_points.front();
                      last = n,
                         n = hole_points.back())
      {
        auto [next_it_begin, next_it_end] = hole_edge_map.equal_range(n);
 
        if (std::distance(next_it_begin, next_it_end) != 2)
          report_error("Bad edge topology found by MeshedHole");
 
        const Node * next = nullptr;
        for (const auto & [key, val] : as_range(next_it_begin, next_it_end))
          {
            libmesh_assert_equal_to(key, n);
            libmesh_ignore(key);
            libmesh_assert_not_equal_to(val.first, n);
 
            // Don't go backwards on the edge we just traversed
            if (val.first == last)
              continue;
 
            // We can support mixes of Edge and Tri-side edges, but we
            // can't do proper error detection on flipped triangles.
            if (val.second != edge_type &&
                val.second != 0)
              {
                if (!edge_type)
                  edge_type = val.second;
                else
                  report_error("MeshedHole sees inconsistent triangle orientations on boundary");
              }
            next = val.first;
          }
 
        // We should never hit the same n twice!
        hole_edge_map.erase(next_it_begin, next_it_end);
 
        hole_points.push_back(next);
      }
 
    for (auto i : make_range(hole_points.size()-1))
      {
        const auto & midpoints = hole_midpoint_map[{hole_points[i],hole_points[i+1]}];
        midpoint_points.insert(midpoint_points.end(),
                               midpoints.begin(), midpoints.end());
      }
 
    hole_points.pop_back();
 
    return hole_points_and_edge_type;
  };
 
  /*
   * If it's not obvious which loop we find is really the loop we
   * want, then we should die with a nice error message.
   */
  int n_negative_areas = 0,
      n_positive_areas = 0,
      n_edgeelem_loops = 0;
 
  std::vector<const Node *> outer_hole_points, outer_mid_points;
  int outer_edge_type = -1;
  Real twice_outer_area = 0,
       abs_twice_outer_area = 0;
 
#ifdef DEBUG
  // Area and edge type, for error reporting
  std::vector<std::pair<Real, int>> areas;
#endif
 
  while (!hole_edge_map.empty()) {
    auto [hole_points, mid_points, edge_type] = extract_edge_vector();
 
    if (edge_type == 0)
    {
      ++n_edgeelem_loops;
      if (n_edgeelem_loops > 1)
        report_error("MeshedHole is confused by multiple loops of Edge elements");
      if (n_positive_areas || n_negative_areas)
        report_error("MeshedHole is confused by meshes with both Edge and 2D-side boundaries");
    }
 
    const std::size_t n_hole_points = hole_points.size();
    if (n_hole_points < 3)
      report_error("Loop with only " + std::to_string(n_hole_points) +
                   " hole edges found in mesh!");
 
    Real twice_this_area = 0;
    const Point p0 = *hole_points[0];
    for (unsigned int i=2; i != n_hole_points; ++i)
      {
        const Point e_0im = *hole_points[i-1] - p0,
                    e_0i  = *hole_points[i] - p0;
 
        twice_this_area += e_0i.cross(e_0im)(2);
      }
 
    auto abs_twice_this_area = std::abs(twice_this_area);
 
    if (((twice_this_area > 0) && edge_type == 2) ||
        ((twice_this_area < 0) && edge_type == 1))
      ++n_positive_areas;
    else if (edge_type != 0)
      ++n_negative_areas;
 
#ifdef DEBUG
    areas.push_back({twice_this_area/2,edge_type});
#endif
 
    if (abs_twice_this_area > abs_twice_outer_area)
      {
        twice_outer_area = twice_this_area;
        abs_twice_outer_area = abs_twice_this_area;
        outer_hole_points = std::move(hole_points);
        outer_mid_points = std::move(mid_points);
        outer_edge_type = edge_type;
      }
  }
 
  _points.resize(outer_hole_points.size());
  std::transform(outer_hole_points.begin(),
                 outer_hole_points.end(),
                 _points.begin(),
                 [](const Node * n){ return Point(*n); });
  _midpoints.resize(outer_mid_points.size());
  std::transform(outer_mid_points.begin(),
                 outer_mid_points.end(),
                 _midpoints.begin(),
                 [](const Node * n){ return Point(*n); });
 
  if (!twice_outer_area)
    report_error("Zero-area MeshedHoles are not currently supported");
 
  // We ordered ourselves counter-clockwise?  But a hole is expected
  // to be clockwise, so use the reverse order.
  if (twice_outer_area > 0)
    {
      std::reverse(_points.begin(), _points.end());
 
      // Our midpoints are numbered e.g.
      // (01a)(01b)(12a)(12b)(23a)(23b)(30a)(30b) for points 0123, but
      // if we reverse to get 3210 then we want our midpoints to be
      // (23b)(23a)(12b)(12a)(01b)(01a)(30b)(30a)
      const unsigned int n_midpoints = _midpoints.size() / _points.size();
      auto split_it = _midpoints.end() - n_midpoints;
      std::reverse(_midpoints.begin(), split_it);
      std::reverse(split_it, _midpoints.end());
    }
 
#ifdef DEBUG
  auto print_areas = [areas](){
    libMesh::out << "Found boundary areas:\n";
    static const std::vector<std::string> edgenames {"E","CW","CCW"};
    for (auto area : areas)
      libMesh::out << '(' << edgenames[area.second] << ' ' <<
        area.first << ')';
    libMesh::out << std::endl;
  };
#else
  auto print_areas = [](){};
#endif
 
  if (((twice_outer_area > 0) && outer_edge_type == 2) ||
      ((twice_outer_area < 0) && outer_edge_type == 1))
    {
      if (n_positive_areas > 1)
        {
          print_areas();
          report_error("MeshedHole found " +
                       std::to_string(n_positive_areas) +
                       " counter-clockwise boundaries and cannot choose one!");
        }
 
    }
  else if (outer_edge_type != 0)
    {
      if (n_negative_areas > 1)
        {
          print_areas();
          report_error("MeshedHole found " +
                       std::to_string(n_negative_areas) +
                       " clockwise boundaries and cannot choose one!");
        }
 
    }
 
  // Hey, no errors!  Broadcast that empty string.
  mesh.comm().broadcast(error_reported);
  mesh.comm().broadcast(_points);
  mesh.comm().broadcast(_midpoints);
}
 
 
unsigned int TriangulatorInterface::MeshedHole::n_points() const
{
  return _points.size();
}
 
 
unsigned int TriangulatorInterface::MeshedHole::n_midpoints() const
{
  libmesh_assert (!(_midpoints.size() % _points.size()));
  return _midpoints.size() / _points.size();
}
 
 
Point TriangulatorInterface::MeshedHole::point(const unsigned int n) const
{
  libmesh_assert_less (n, _points.size());
  return _points[n];
}
 
 
Point TriangulatorInterface::MeshedHole::midpoint(const unsigned int m,
                                                  const unsigned int n) const
{
  const unsigned int n_mid = this->n_midpoints();
  libmesh_assert_less (m, n_mid);
  libmesh_assert_less (n, _points.size());
  return _midpoints[n*n_mid+m];
}
 
 
Point TriangulatorInterface::MeshedHole::inside() const
{
  // This is expensive to compute, so only do it when we first need it
  if (_center(0) == std::numeric_limits<Real>::max())
    _center = this->calculate_inside_point();
 
  return _center;
}
 
 
 
} // namespace libMesh