XYIncrementalDelaunay Class Reference

Constrained Delaunay triangulation of a set of points in the plane, built one point at a time. More...

#include <XYIncrementalDelaunay.h>

Classes
struct	Point2D
	A point of the triangulation, held as plain coordinates. More...
struct	Triangle
	A triangle of the triangulation. More...

Public Types
using	Segment = std::pair< std::size_t, std::size_t >
	A constrained segment, held as a vertex id pair with the smaller id first.

Public Member Functions
void	initialize (const std::vector< Point2D > &points, const std::vector< Segment > &segments)
	Triangulates points and recovers every entry of segments as an edge of the result.
std::size_t	insertPoint (const Point2D &p)
	Inserts a point, restoring the constrained Delaunay property around it.
void	insertSegment (std::size_t v0, std::size_t v1)
	Makes the segment between two vertices an edge of the triangulation and records it as constrained, so that later insertions neither remove nor cross it.
std::size_t	numPoints () const
const Point2D &	point (std::size_t id) const
std::vector< Triangle >	getTriangles () const
bool	isConstrainedSegment (std::size_t v0, std::size_t v1) const
const std::set< Segment > &	constrainedSegments () const
std::vector< std::string >	checkInvariants () const
	Checks everything this class promises: that every triangle is counter-clockwise, that the neighbor entries agree with each other, that every constrained segment is an edge, and that every edge that is not constrained is locally Delaunay.
std::vector< std::string >	checkEmptyCircumcircle () const
	Tests every vertex against the circumcircle of every triangle.

Static Public Member Functions
static Segment	makeSegment (std::size_t v0, std::size_t v1)

Static Public Attributes
static constexpr std::size_t	invalid_index = std::numeric_limits<std::size_t>::max()
	Sentinel for a vertex, triangle or neighbor that does not exist.

Private Member Functions
const double *	xy (std::size_t v) const
std::size_t	toInternal (std::size_t id) const
std::size_t	toCaller (std::size_t v) const
bool	isConstrainedEdge (std::size_t v0, std::size_t v1) const
unsigned int	localVertexIndex (std::size_t t, std::size_t v) const
unsigned int	localEdgeIndex (std::size_t t, const Segment &edge) const
bool	containsPoint (std::size_t t, const Point2D &p) const
std::size_t	locate (const Point2D &p) const
	Finds the triangle a point falls in by walking from the triangle the last insertion produced, which keeps the search short when successive points are close together.
void	growCavity (std::size_t seed, std::size_t v_new, std::set< std::size_t > &cavity) const
	Collects the triangles of the triangulation built so far that have to make way for a new vertex, the Delaunay cavity of its insertion, starting from the triangle it falls in and spreading to every neighbor whose circumcircle contains it.
std::vector< std::size_t >	retriangulate (const std::vector< std::size_t > &removed, const std::vector< std::array< std::size_t, 3 > > &added)
	Swaps one triangulation of a region for another, reusing the slots of the triangles it removes and rebuilding the neighbor entries both inside the region and along its boundary.
void	triangulatePseudopolygon (std::size_t v_start, std::size_t v_end, const std::vector< std::size_t > &chain, std::size_t first, std::size_t last, std::vector< std::array< std::size_t, 3 > > &triangles) const
	Triangulates a polygon whose vertices are all visible from one of its edges, which is the shape a recovered segment leaves on either side of itself.
bool	isStrictlyBetween (std::size_t v_first, std::size_t v_mid, std::size_t v_last) const
std::string	vertexName (std::size_t v) const

Static Private Member Functions
static bool	isBounding (std::size_t v)

Private Attributes
std::vector< Point2D >	_vertices
	The bounding triangle vertices followed by the caller's points.
std::vector< Triangle >	_triangles
	The triangles, every one of them live and counter-clockwise.
std::vector< std::size_t >	_vertex_triangle
	One triangle touching each vertex of _vertices, which is where a walk around it starts.
std::set< Segment >	_constraints
	The constrained segments, in caller vertex ids with the smaller id first.
std::size_t	_last_triangle = invalid_index
	A triangle the last insertion produced, which is where the next point walk starts.

Static Private Attributes
static constexpr std::size_t	_num_bounding = 3
	The number of bounding triangle vertices padding the front of the vertex list.
static constexpr double	_bounding_reach = 1000.0
	How far the bounding triangle reaches beyond the points, as a multiple of their extent.

Detailed Description

Constrained Delaunay triangulation of a set of points in the plane, built one point at a time.

Points are inserted with the Bowyer-Watson construction: the triangle containing the new point is found by walking the triangulation, the Delaunay cavity (the triangles whose circumcircle contains the point) is deleted, and the point is connected to the boundary of that cavity. Constrained segments are edges that the triangulation is required to contain; the Delaunay cavity never grows across one, and a segment that is not an edge once the points are in is recovered by retriangulating the strip of triangles it crosses.

After every insertion no vertex lies strictly inside the circumcircle of any triangle, except where a constrained segment separates the two. Cocircular and collinear configurations are resolved consistently rather than arbitrarily, so the same input always gives the same triangulation. Every geometric decision is made with the exact predicates in framework/contrib/predicates, so no tolerance enters the algorithm.

The triangulation carries no mesh types: points are plain coordinate pairs and triangles are plain index triples, so it can be exercised on its own.

Internally the point set is padded with the three vertices of a bounding triangle that encloses everything, which is what lets the Delaunay cavity of every insertion be a closed polygon. Those three vertices are not part of the caller's point set and never appear in the public interface: the vertex ids used here run from 0 to numPoints() - 1 in the order the points were added, and getTriangles() drops the triangles that touch the bounding triangle.

Definition at line 45 of file XYIncrementalDelaunay.h.

Member Typedef Documentation

Segment

using XYIncrementalDelaunay::Segment = std::pair<std::size_t, std::size_t>

A constrained segment, held as a vertex id pair with the smaller id first.

Definition at line 72 of file XYIncrementalDelaunay.h.

Member Function Documentation

checkEmptyCircumcircle()

std::vector< std::string > XYIncrementalDelaunay::checkEmptyCircumcircle

(

)

const

Tests every vertex against the circumcircle of every triangle.

This ignores constrained segments, which a vertex is allowed to sit behind, so it is the right criterion only for a triangulation with no constrained segments; use checkInvariants() otherwise.

Returns

One entry describing each vertex that lies strictly inside a circumcircle, empty if there are none

Definition at line 710 of file XYIncrementalDelaunay.C.

{
  Moose::initPredicates();
 
  std::vector<std::string> violations;
  for (const auto t : index_range(_triangles))
  {
    const auto & triangle = _triangles[t];
    for (const auto v : index_range(_vertices))
    {
      if (v == triangle.vertices[0] || v == triangle.vertices[1] || v == triangle.vertices[2])
        continue;
      if (moose_incircle(
              xy(triangle.vertices[0]), xy(triangle.vertices[1]), xy(triangle.vertices[2]), xy(v)) >
          0.0)
        violations.push_back(vertexName(v) + " is inside the circumcircle of triangle " +
                             std::to_string(t));
    }
  }
  return violations;
}

checkInvariants()

std::vector< std::string > XYIncrementalDelaunay::checkInvariants

(

)

const

Checks everything this class promises: that every triangle is counter-clockwise, that the neighbor entries agree with each other, that every constrained segment is an edge, and that every edge that is not constrained is locally Delaunay.

A triangulation whose unconstrained edges are all locally Delaunay is a constrained Delaunay triangulation.

Returns

One entry describing each violation found, empty if there are none

Definition at line 643 of file XYIncrementalDelaunay.C.

{
  Moose::initPredicates();
 
  std::vector<std::string> violations;
  std::set<Segment> edges;
 
  for (const auto t : index_range(_triangles))
  {
    const auto & triangle = _triangles[t];
    if (moose_orient2d(
            xy(triangle.vertices[0]), xy(triangle.vertices[1]), xy(triangle.vertices[2])) <= 0.0)
      violations.push_back("triangle " + std::to_string(t) +
                           " is not counter-clockwise, so its area is zero or negative");
 
    for (const auto i : make_range(3u))
    {
      const auto v0 = triangle.vertices[(i + 1) % 3];
      const auto v1 = triangle.vertices[(i + 2) % 3];
      edges.insert(makeSegment(v0, v1));
 
      const auto n = triangle.neighbors[i];
      if (n == invalid_index)
        continue;
      if (n >= _triangles.size())
      {
        violations.push_back("triangle " + std::to_string(t) + " has neighbor " +
                             std::to_string(n) + ", which is not a triangle");
        continue;
      }
 
      unsigned int j = 3;
      for (const auto k : make_range(3u))
        if (makeSegment(_triangles[n].vertices[(k + 1) % 3], _triangles[n].vertices[(k + 2) % 3]) ==
            makeSegment(v0, v1))
          j = k;
 
      if (j == 3)
      {
        violations.push_back("triangles " + std::to_string(t) + " and " + std::to_string(n) +
                             " are neighbors but share no edge");
        continue;
      }
      if (_triangles[n].neighbors[j] != t)
        violations.push_back("triangle " + std::to_string(t) + " has neighbor " +
                             std::to_string(n) + ", which does not have it back");
 
      if (!isConstrainedEdge(v0, v1) && moose_incircle(xy(triangle.vertices[0]),
                                                       xy(triangle.vertices[1]),
                                                       xy(triangle.vertices[2]),
                                                       xy(_triangles[n].vertices[j])) > 0.0)
        violations.push_back("the edge between " + vertexName(v0) + " and " + vertexName(v1) +
                             " is not constrained and is not locally Delaunay, because " +
                             vertexName(_triangles[n].vertices[j]) +
                             " is inside the circumcircle of triangle " + std::to_string(t));
    }
  }
 
  for (const auto & [v0, v1] : _constraints)
    if (edges.count(makeSegment(toInternal(v0), toInternal(v1))) == 0)
      violations.push_back("constrained segment (" + std::to_string(v0) + ", " +
                           std::to_string(v1) + ") is not an edge of the triangulation");
 
  return violations;
}

constrainedSegments()

const std::set< Segment > & XYIncrementalDelaunay::constrainedSegments ( ) const

inline

Returns

The constrained segments, each with the smaller vertex id first

Definition at line 127 of file XYIncrementalDelaunay.h.

{ return _constraints; }

Referenced by XYFrontalDelaunayGenerator::advanceFront(), and XYFrontalDelaunayGenerator::recordSplitBoundaryIds().

containsPoint()

bool XYIncrementalDelaunay::containsPoint ( std::size_t  t, const Point2D &  p  ) const

private

Returns

Whether a point lies inside or on the boundary of a triangle

Definition at line 110 of file XYIncrementalDelaunay.C.

{
  for (const auto i : make_range(3u))
    if (moose_orient2d(xy(_triangles[t].vertices[(i + 1) % 3]),
                       xy(_triangles[t].vertices[(i + 2) % 3]),
                       coords(p)) < 0.0)
      return false;
  return true;
}

Referenced by locate().

getTriangles()

std::vector< XYIncrementalDelaunay::Triangle > XYIncrementalDelaunay::getTriangles

(

)

const

Returns

The triangles covering the convex hull of the points, counter-clockwise, ordered so that the same input always gives the same list. The neighbor entries index into this same list and are invalid_index on the hull, where no triangle lies across the edge.

Definition at line 616 of file XYIncrementalDelaunay.C.

{
  std::vector<std::size_t> position(_triangles.size(), invalid_index);
  std::size_t count = 0;
  for (const auto t : index_range(_triangles))
    if (!isBounding(_triangles[t].vertices[0]) && !isBounding(_triangles[t].vertices[1]) &&
        !isBounding(_triangles[t].vertices[2]))
      position[t] = count++;
 
  std::vector<Triangle> triangles(count);
  for (const auto t : index_range(_triangles))
  {
    if (position[t] == invalid_index)
      continue;
 
    auto & out = triangles[position[t]];
    for (const auto i : make_range(3u))
    {
      out.vertices[i] = toCaller(_triangles[t].vertices[i]);
      const auto n = _triangles[t].neighbors[i];
      out.neighbors[i] = n == invalid_index ? invalid_index : position[n];
    }
  }
  return triangles;
}

Referenced by XYFrontalDelaunayGenerator::advanceFront(), and XYFrontalDelaunayGenerator::buildTriangleMesh().

growCavity()

void XYIncrementalDelaunay::growCavity ( std::size_t  seed, std::size_t  v_new, std::set< std::size_t > &  cavity  ) const

private

Collects the triangles of the triangulation built so far that have to make way for a new vertex, the Delaunay cavity of its insertion, starting from the triangle it falls in and spreading to every neighbor whose circumcircle contains it.

Constrained segments stop the spread. A neighbor is also taken in when the shared edge would otherwise produce a triangle of zero or negative area, which is what keeps the cavity star shaped about the vertex.

Parameters

seed	The triangle the new vertex falls in
v_new	The internal vertex id of the new vertex
cavity	Filled with the triangles to remove

Definition at line 185 of file XYIncrementalDelaunay.C.

{
  cavity.insert(seed);
  std::vector<std::size_t> pending{seed};
 
  while (!pending.empty())
  {
    const auto t = pending.back();
    pending.pop_back();
 
    for (const auto i : make_range(3u))
    {
      const auto v0 = _triangles[t].vertices[(i + 1) % 3];
      const auto v1 = _triangles[t].vertices[(i + 2) % 3];
      const auto side = moose_orient2d(xy(v0), xy(v1), xy(v_new));
 
      if (isConstrainedEdge(v0, v1))
      {
        if (side <= 0.0)
          mooseError("XYIncrementalDelaunay: the point (",
                     _vertices[v_new].x,
                     ", ",
                     _vertices[v_new].y,
                     ") is on or beyond constrained segment (",
                     toCaller(v0),
                     ", ",
                     toCaller(v1),
                     "), so it lies outside the region that segment bounds.");
        continue;
      }
 
      const auto n = _triangles[t].neighbors[i];
      if (n == invalid_index || cavity.count(n) > 0)
        continue;
 
      // Taking the neighbor in when the shared edge would give a triangle of zero or negative area
      // is what keeps the cavity star shaped about the new vertex; with exact predicates that only
      // happens when the new vertex falls exactly on the edge.
      if (side <= 0.0 || moose_incircle(xy(_triangles[n].vertices[0]),
                                        xy(_triangles[n].vertices[1]),
                                        xy(_triangles[n].vertices[2]),
                                        xy(v_new)) > 0.0)
      {
        cavity.insert(n);
        pending.push_back(n);
      }
    }
  }
}

Referenced by insertPoint().

initialize()

void XYIncrementalDelaunay::initialize	(	const std::vector< Point2D > &	points,
		const std::vector< Segment > &	segments
	)

Triangulates points and recovers every entry of segments as an edge of the result.

Points are inserted in the order given, so vertex id i refers to points[i]. Any later insertPoint() must lie within the region this call encloses, which is the bounding box of points grown by a large multiple of its own size; the interior points a frontal advance adds always do.

Parameters

points	The points to triangulate, which must not contain two identical points
segments	The segments the triangulation is required to contain, as vertex id pairs in either order

Definition at line 346 of file XYIncrementalDelaunay.C.

{
  Moose::initPredicates();
 
  if (points.empty())
    mooseError("XYIncrementalDelaunay: initialize() needs at least one point.");
 
  _vertices.clear();
  _triangles.clear();
  _vertex_triangle.clear();
  _constraints.clear();
 
  auto x_min = points.front().x;
  auto x_max = x_min;
  auto y_min = points.front().y;
  auto y_max = y_min;
  for (const auto & p : points)
  {
    x_min = std::min(x_min, p.x);
    x_max = std::max(x_max, p.x);
    y_min = std::min(y_min, p.y);
    y_max = std::max(y_max, p.y);
  }
 
  // A triangle far enough outside the points that all of them are strictly inside it. Inserting a
  // point deletes the Delaunay cavity of that insertion, the triangles whose circumcircle contains
  // the point (see growCavity()), and fans the hole they leave out from the point. With every
  // point inside the bounding triangle that hole is always a closed polygon, so no insertion ever
  // has to extend the convex hull of the triangulation.
  const auto reach = _bounding_reach * std::max({x_max - x_min, y_max - y_min, 1.0});
  const auto x_mid = 0.5 * (x_min + x_max);
  const auto y_mid = 0.5 * (y_min + y_max);
  _vertices.push_back({x_mid - reach, y_mid - reach});
  _vertices.push_back({x_mid + reach, y_mid - reach});
  _vertices.push_back({x_mid, y_mid + reach});
  _vertex_triangle.assign(_num_bounding, std::size_t(0));
 
  Triangle bounding;
  bounding.vertices = {0, 1, 2};
  bounding.neighbors = {invalid_index, invalid_index, invalid_index};
  _triangles.push_back(bounding);
  _last_triangle = 0;
 
  for (const auto i : index_range(points))
  {
    const auto id = insertPoint(points[i]);
    if (id != i)
      mooseError("XYIncrementalDelaunay: points ",
                 id,
                 " and ",
                 i,
                 " are the same point; initialize() needs the points to be distinct, because the "
                 "constrained segments refer to them by position.");
  }
 
  for (const auto & [v0, v1] : segments)
    insertSegment(v0, v1);
}

insertPoint()

std::size_t XYIncrementalDelaunay::insertPoint

(

const Point2D &

p

)

Inserts a point, restoring the constrained Delaunay property around it.

Parameters

p	The point to insert

Returns

The vertex id of p. A point identical to an existing vertex is not inserted a second time; the id of that vertex is returned and the triangulation is left alone. A point that falls exactly on a constrained segment splits it, replacing that segment by the two halves the point divides it into.

Definition at line 407 of file XYIncrementalDelaunay.C.

{
  Moose::initPredicates();
 
  if (_triangles.empty())
    mooseError("XYIncrementalDelaunay: initialize() has to run before a point can be inserted.");
 
  const auto seed = locate(p);
  if (seed == invalid_index)
    mooseError("XYIncrementalDelaunay: the point (",
               p.x,
               ", ",
               p.y,
               ") is outside the region initialize() enclosed, so it cannot be inserted.");
 
  // Two points count as the same point only when their coordinates agree exactly, which is the
  // criterion the exact predicates use as well.
  for (const auto v : _triangles[seed].vertices)
    if (_vertices[v].x == p.x && _vertices[v].y == p.y)
    {
      if (isBounding(v))
        mooseError("XYIncrementalDelaunay: the point (",
                   p.x,
                   ", ",
                   p.y,
                   ") is a vertex of the bounding triangle, so it is far outside the region "
                   "initialize() enclosed.");
      return toCaller(v);
    }
 
  // A point landing on a constrained segment divides it, because the Delaunay cavity may not cross
  // a constrained segment and a segment with a vertex on it can no longer be a single edge.
  Segment split(invalid_index, invalid_index);
  for (const auto i : make_range(3u))
  {
    const auto v0 = _triangles[seed].vertices[(i + 1) % 3];
    const auto v1 = _triangles[seed].vertices[(i + 2) % 3];
    if (moose_orient2d(xy(v0), xy(v1), coords(p)) == 0.0 && isConstrainedEdge(v0, v1))
    {
      split = makeSegment(toCaller(v0), toCaller(v1));
      _constraints.erase(split);
      break;
    }
  }
 
  _vertices.push_back(p);
  _vertex_triangle.push_back(invalid_index);
  const auto v_new = _vertices.size() - 1;
 
  std::set<std::size_t> cavity;
  growCavity(seed, v_new, cavity);
 
  std::vector<std::array<std::size_t, 3>> added;
  for (const auto t : cavity)
    for (const auto i : make_range(3u))
    {
      const auto n = _triangles[t].neighbors[i];
      if (n == invalid_index || cavity.count(n) == 0)
        added.push_back(
            {_triangles[t].vertices[(i + 1) % 3], _triangles[t].vertices[(i + 2) % 3], v_new});
    }
 
  mooseAssert(added.size() == cavity.size() + 2,
              "The Delaunay cavity of an insertion is a disc, whose boundary has two more edges "
              "than the disc has triangles");
 
  const std::vector<std::size_t> removed(cavity.begin(), cavity.end());
  _last_triangle = retriangulate(removed, added).front();
 
  const auto id = toCaller(v_new);
  if (split.first != invalid_index)
  {
    insertSegment(split.first, id);
    insertSegment(id, split.second);
  }
  return id;
}

Referenced by XYFrontalDelaunayGenerator::advanceFront(), and initialize().

insertSegment()

void XYIncrementalDelaunay::insertSegment	(	std::size_t	v0,
		std::size_t	v1
	)

Makes the segment between two vertices an edge of the triangulation and records it as constrained, so that later insertions neither remove nor cross it.

Doing this to a segment that is already an edge only records the constraint.

A segment that runs exactly through a third vertex cannot be a single edge, and is recorded as the pieces that vertex divides it into rather than as itself.

Parameters

v0	The vertex id of one end of the segment
v1	The vertex id of the other end of the segment

Definition at line 486 of file XYIncrementalDelaunay.C.

{
  Moose::initPredicates();
 
  if (v0 == v1)
    mooseError(
        "XYIncrementalDelaunay: a constrained segment needs two different vertices, but both "
        "ends of this one are vertex ",
        v0,
        ".");
 
  const auto v_from = toInternal(v0);
  const auto v_to = toInternal(v1);
 
  // Turn around v_from until the triangle the segment leaves through is found. The segment is
  // already an edge if v_to is met on the way, and needs splitting if some other vertex is.
  auto entered = invalid_index;
  auto right = invalid_index;
  auto left = invalid_index;
  const auto start = _vertex_triangle[v_from];
  mooseAssert(start < _triangles.size(), "Every vertex records a triangle it belongs to");
 
  auto current = start;
  do
  {
    const auto i = localVertexIndex(current, v_from);
    const auto next_v = _triangles[current].vertices[(i + 1) % 3];
    const auto far_v = _triangles[current].vertices[(i + 2) % 3];
 
    if (next_v == v_to)
    {
      _constraints.insert(makeSegment(v0, v1));
      return;
    }
 
    const auto next_side = moose_orient2d(xy(v_from), xy(next_v), xy(v_to));
 
    // Nothing has been changed yet, so recovering the two halves from scratch is safe.
    if (next_side == 0.0 && isStrictlyBetween(v_from, next_v, v_to))
    {
      insertSegment(v0, toCaller(next_v));
      insertSegment(toCaller(next_v), v1);
      return;
    }
 
    if (next_side > 0.0 && moose_orient2d(xy(v_from), xy(far_v), xy(v_to)) < 0.0)
    {
      entered = current;
      right = next_v;
      left = far_v;
      break;
    }
 
    current = _triangles[current].neighbors[(i + 1) % 3];
  } while (current != start && current != invalid_index);
 
  if (entered == invalid_index)
    mooseError("XYIncrementalDelaunay: no triangle around vertex ",
               v0,
               " is entered by the segment to vertex ",
               v1,
               ".");
 
  // Follow the segment across the triangulation, keeping the vertices it passes on either side.
  std::vector<std::size_t> crossed{entered};
  std::vector<std::size_t> right_chain{right};
  std::vector<std::size_t> left_chain{left};
 
  while (true)
  {
    if (isConstrainedEdge(right, left))
      mooseError("XYIncrementalDelaunay: the segment from vertex ",
                 v0,
                 " to vertex ",
                 v1,
                 " crosses constrained segment (",
                 toCaller(right),
                 ", ",
                 toCaller(left),
                 "); constrained segments may only meet at their ends.");
 
    const auto edge = makeSegment(right, left);
    const auto next_t = _triangles[crossed.back()].neighbors[localEdgeIndex(crossed.back(), edge)];
    if (next_t == invalid_index)
      mooseError("XYIncrementalDelaunay: the segment from vertex ",
                 v0,
                 " to vertex ",
                 v1,
                 " leaves the triangulation before it reaches its end.");
 
    const auto apex = _triangles[next_t].vertices[localEdgeIndex(next_t, edge)];
    if (apex == v_to)
    {
      crossed.push_back(next_t);
      break;
    }
 
    const auto side = moose_orient2d(xy(v_from), xy(v_to), xy(apex));
    if (side == 0.0)
    {
      insertSegment(v0, toCaller(apex));
      insertSegment(toCaller(apex), v1);
      return;
    }
 
    crossed.push_back(next_t);
    if (side > 0.0)
    {
      left_chain.push_back(apex);
      left = apex;
    }
    else
    {
      right_chain.push_back(apex);
      right = apex;
    }
  }
 
  // The segment leaves a polygon on either side of itself, every vertex of which the segment sees,
  // so triangulating each on its own restores the Delaunay property everywhere.
  std::vector<std::array<std::size_t, 3>> added;
  const std::vector<std::size_t> left_polygon(left_chain.rbegin(), left_chain.rend());
  triangulatePseudopolygon(v_from, v_to, left_polygon, 0, left_polygon.size(), added);
  triangulatePseudopolygon(v_to, v_from, right_chain, 0, right_chain.size(), added);
 
  retriangulate(crossed, added);
  _constraints.insert(makeSegment(v0, v1));
}

Referenced by initialize(), and insertSegment().

isBounding()

static bool XYIncrementalDelaunay::isBounding ( std::size_t  v )

inlinestaticprivate

Returns

Whether an internal vertex belongs to the bounding triangle rather than to the caller

Definition at line 158 of file XYIncrementalDelaunay.h.

{ return v < _num_bounding; }

Referenced by getTriangles(), insertPoint(), isConstrainedEdge(), toCaller(), and vertexName().

isConstrainedEdge()

bool XYIncrementalDelaunay::isConstrainedEdge ( std::size_t  v0, std::size_t  v1  ) const

private

Returns

Whether the edge between two internal vertices is a constrained segment

Definition at line 77 of file XYIncrementalDelaunay.C.

{
  if (isBounding(v0) || isBounding(v1))
    return false;
  return _constraints.count(makeSegment(toCaller(v0), toCaller(v1))) > 0;
}

Referenced by checkInvariants(), growCavity(), insertPoint(), and insertSegment().

isConstrainedSegment()

bool XYIncrementalDelaunay::isConstrainedSegment	(	std::size_t	v0,
		std::size_t	v1
	)		const

Returns

Whether the segment between two vertices is constrained, given its ends in either order

Definition at line 71 of file XYIncrementalDelaunay.C.

{
  return _constraints.count(makeSegment(v0, v1)) > 0;
}

isStrictlyBetween()

bool XYIncrementalDelaunay::isStrictlyBetween ( std::size_t  v_first, std::size_t  v_mid, std::size_t  v_last  ) const

private

Returns

Whether v_mid lies strictly between v_first and v_last, which must all be collinear

Definition at line 121 of file XYIncrementalDelaunay.C.

{
  mooseAssert(moose_orient2d(xy(v_first), xy(v_mid), xy(v_last)) == 0.0,
              "The three vertices have to be collinear for a coordinate comparison to order them");
 
  // The three are collinear, so whichever coordinate separates the two ends orders all three.
  const auto & first = _vertices[v_first];
  const auto & mid = _vertices[v_mid];
  const auto & last = _vertices[v_last];
  if (first.x != last.x)
    return (first.x < mid.x && mid.x < last.x) || (last.x < mid.x && mid.x < first.x);
  return (first.y < mid.y && mid.y < last.y) || (last.y < mid.y && mid.y < first.y);
}

Referenced by insertSegment().

localEdgeIndex()

unsigned int XYIncrementalDelaunay::localEdgeIndex ( std::size_t  t, const Segment &  edge  ) const

private

Returns

The position within a triangle of the edge opposite the given edge's vertex

Definition at line 94 of file XYIncrementalDelaunay.C.

{
  for (const auto i : make_range(3u))
    if (makeSegment(_triangles[t].vertices[(i + 1) % 3], _triangles[t].vertices[(i + 2) % 3]) ==
        edge)
      return i;
  mooseError("XYIncrementalDelaunay: triangle ",
             t,
             " does not have an edge between vertices ",
             edge.first,
             " and ",
             edge.second,
             ".");
}

Referenced by insertSegment(), and retriangulate().

localVertexIndex()

unsigned int XYIncrementalDelaunay::localVertexIndex ( std::size_t  t, std::size_t  v  ) const

private

Returns

The position of an internal vertex within a triangle

Definition at line 85 of file XYIncrementalDelaunay.C.

{
  for (const auto i : make_range(3u))
    if (_triangles[t].vertices[i] == v)
      return i;
  mooseError("XYIncrementalDelaunay: triangle ", t, " does not have vertex ", v, ".");
}

Referenced by insertSegment().

locate()

std::size_t XYIncrementalDelaunay::locate ( const Point2D &  p ) const

private

Finds the triangle a point falls in by walking from the triangle the last insertion produced, which keeps the search short when successive points are close together.

Falls back to a scan of every triangle if the walk leaves the triangulation or fails to settle.

Returns

The triangle containing p, or invalid_index if p is outside the triangulation

Definition at line 146 of file XYIncrementalDelaunay.C.

{
  // start from the last insertion if existing
  auto current = _last_triangle < _triangles.size() ? _last_triangle : std::size_t(0);
 
  // The walk always arrives in a Delaunay triangulation, but a constrained one can in principle
  // send it round in a circle, so it gives up and scans rather than spin.
  auto steps_left = 2 * _triangles.size() + 8;
  while (steps_left > 0)
  {
    --steps_left;
 
    auto next = invalid_index;
    bool inside = true;
    for (const auto i : make_range(3u))
      if (moose_orient2d(xy(_triangles[current].vertices[(i + 1) % 3]),
                         xy(_triangles[current].vertices[(i + 2) % 3]),
                         coords(p)) < 0.0)
      {
        inside = false;
        next = _triangles[current].neighbors[i];
        break;
      }
 
    if (inside)
      return current;
    if (next == invalid_index)
      break;
    current = next;
  }
 
  // search all triangles if starting from the last inserted one did not succeed
  for (const auto t : index_range(_triangles))
    if (containsPoint(t, p))
      return t;
  return invalid_index;
}

Referenced by insertPoint().

makeSegment()

XYIncrementalDelaunay::Segment XYIncrementalDelaunay::makeSegment ( std::size_t  v0, std::size_t  v1  )

static

Returns

A vertex id pair with the smaller id first, the form constrained segments are held in

Definition at line 65 of file XYIncrementalDelaunay.C.

{
  return {std::min(v0, v1), std::max(v0, v1)};
}

Referenced by XYFrontalDelaunayGenerator::appendLoop(), XYFrontalDelaunayGenerator::buildTriangleMesh(), checkInvariants(), insertPoint(), insertSegment(), isConstrainedEdge(), isConstrainedSegment(), localEdgeIndex(), XYFrontalDelaunayGenerator::recordSplitBoundaryIds(), and retriangulate().

numPoints()

std::size_t XYIncrementalDelaunay::numPoints

(

)

const

Returns

The number of points in the triangulation, whose vertex ids run from 0 to this minus 1

Definition at line 34 of file XYIncrementalDelaunay.C.

{
  return _vertices.empty() ? 0 : _vertices.size() - _num_bounding;
}

Referenced by XYFrontalDelaunayGenerator::advanceFront(), and toInternal().

point()

const XYIncrementalDelaunay::Point2D & XYIncrementalDelaunay::point

(

std::size_t

id

)

const

Returns

The coordinates of the vertex with the given id

Definition at line 59 of file XYIncrementalDelaunay.C.

{
  return _vertices[toInternal(id)];
}

Referenced by XYFrontalDelaunayGenerator::buildTriangleMesh(), XYFrontalDelaunayGenerator::collectFront(), XYFrontalDelaunayGenerator::hasVertexWithin(), and XYFrontalDelaunayGenerator::placePoint().

retriangulate()

std::vector< std::size_t > XYIncrementalDelaunay::retriangulate ( const std::vector< std::size_t > &  removed, const std::vector< std::array< std::size_t, 3 > > &  added  )

private

Swaps one triangulation of a region for another, reusing the slots of the triangles it removes and rebuilding the neighbor entries both inside the region and along its boundary.

Parameters

removed	The triangles covering the region, in any order
added	The vertices of the triangles to cover it with, counter-clockwise

Returns

The slots the added triangles were placed in

Definition at line 238 of file XYIncrementalDelaunay.C.

{
  mooseAssert(added.size() >= removed.size(),
              "A region takes at least as many triangles to cover as it is emptied of");
 
  const std::set<std::size_t> emptied(removed.begin(), removed.end());
 
  // What lies outside the region being replaced, keyed on the edge it is joined along.
  std::map<Segment, std::size_t> outside;
  for (const auto t : emptied)
    for (const auto i : make_range(3u))
    {
      const auto n = _triangles[t].neighbors[i];
      if (n != invalid_index && emptied.count(n) == 0)
        outside.emplace(
            makeSegment(_triangles[t].vertices[(i + 1) % 3], _triangles[t].vertices[(i + 2) % 3]),
            n);
    }
 
  // Filling the emptied slots in id order, rather than in the order the region happened to be
  // walked, is what makes the triangle numbering repeatable from one run to the next.
  std::vector<std::size_t> slots(emptied.begin(), emptied.end());
  while (slots.size() < added.size())
  {
    slots.push_back(_triangles.size());
    _triangles.emplace_back();
  }
 
  for (const auto i : index_range(added))
  {
    auto & triangle = _triangles[slots[i]];
    triangle.vertices = added[i];
    triangle.neighbors = {invalid_index, invalid_index, invalid_index};
    for (const auto v : triangle.vertices)
      _vertex_triangle[v] = slots[i];
  }
 
  // An edge two new triangles share joins them to each other; one only a single new triangle has is
  // on the boundary of the region and joins it to what was already outside.
  std::map<Segment, std::pair<std::size_t, unsigned int>> open_edges;
  for (const auto i : index_range(added))
  {
    const auto t = slots[i];
    for (const auto e : make_range(3u))
    {
      const auto edge =
          makeSegment(_triangles[t].vertices[(e + 1) % 3], _triangles[t].vertices[(e + 2) % 3]);
      const auto it = open_edges.find(edge);
      if (it == open_edges.end())
        open_edges.emplace(edge, std::make_pair(t, e));
      else
      {
        const auto [other_t, other_e] = it->second;
        _triangles[t].neighbors[e] = other_t;
        _triangles[other_t].neighbors[other_e] = t;
        open_edges.erase(it);
      }
    }
  }
 
  for (const auto & [edge, entry] : open_edges)
  {
    const auto [t, e] = entry;
    const auto it = outside.find(edge);
    // true outer boundary (for now)
    if (it == outside.end())
      continue;
    // boundary with 'outside' the retriangulated area
    _triangles[t].neighbors[e] = it->second;
    _triangles[it->second].neighbors[localEdgeIndex(it->second, edge)] = t;
  }
 
  return slots;
}

Referenced by insertSegment().

toCaller()

std::size_t XYIncrementalDelaunay::toCaller ( std::size_t  v ) const

private

Returns

The caller vertex id of an internal vertex id, which must not be a bounding vertex

Definition at line 52 of file XYIncrementalDelaunay.C.

{
  mooseAssert(!isBounding(v), "A bounding triangle vertex has no caller vertex id");
  return v - _num_bounding;
}

Referenced by getTriangles(), growCavity(), insertPoint(), insertSegment(), isConstrainedEdge(), and vertexName().

toInternal()

std::size_t XYIncrementalDelaunay::toInternal ( std::size_t  id ) const

private

Returns

The internal vertex id of a caller vertex id, which must be in range

Definition at line 40 of file XYIncrementalDelaunay.C.

{
  if (id >= numPoints())
    mooseError("XYIncrementalDelaunay: vertex id ",
               id,
               " does not exist; the triangulation has ",
               numPoints(),
               " points.");
  return id + _num_bounding;
}

Referenced by checkInvariants(), insertSegment(), and point().

triangulatePseudopolygon()

void XYIncrementalDelaunay::triangulatePseudopolygon ( std::size_t  v_start, std::size_t  v_end, const std::vector< std::size_t > &  chain, std::size_t  first, std::size_t  last, std::vector< std::array< std::size_t, 3 > > &  triangles  ) const

private

Triangulates a polygon whose vertices are all visible from one of its edges, which is the shape a recovered segment leaves on either side of itself.

The polygon runs counter-clockwise from v_start to v_end and back along chain[first] to chain[last - 1]. Splitting it at the vertex whose circumcircle holds no other and recursing on the two halves makes it Delaunay.

Parameters

v_start	The internal vertex at the start of the edge every vertex sees
v_end	The internal vertex at the end of that edge
chain	The internal vertices of the rest of the polygon, in counter-clockwise order
first	The position in chain where this polygon begins
last	One past the position in chain where this polygon ends
triangles	Appended with the vertices of the triangles covering the polygon

Definition at line 315 of file XYIncrementalDelaunay.C.

{
  // An empty chain leaves nothing but the edge from v_start to v_end itself, so there is no
  // polygon to cover
  if (first == last)
    return;
 
  // The triangle the edge (v_start, v_end) belongs to in a Delaunay triangulation of the polygon
  // is the one whose apex leaves no other chain vertex inside its circumcircle. Every chain vertex
  // sees the edge, so whenever a vertex lies inside the circumcircle of the current apex, that
  // vertex takes over as the apex; the one left standing has an empty circumcircle
  auto best = first;
  for (const auto i : make_range(first + 1, last))
    if (moose_incircle(xy(v_start), xy(v_end), xy(chain[best]), xy(chain[i])) > 0.0)
      best = i;
 
  // That triangle splits the polygon in two: the chain vertices before the apex form a polygon
  // seen from the new edge (apex, v_end), and those after it a polygon seen from (v_start, apex),
  // each of which is covered the same way
  triangles.push_back({v_start, v_end, chain[best]});
  triangulatePseudopolygon(chain[best], v_end, chain, first, best, triangles);
  triangulatePseudopolygon(v_start, chain[best], chain, best + 1, last, triangles);
}

Referenced by insertSegment(), and triangulatePseudopolygon().

vertexName()

std::string XYIncrementalDelaunay::vertexName ( std::size_t  v ) const

private

Returns

A name for an internal vertex, for the messages the check methods return

Definition at line 138 of file XYIncrementalDelaunay.C.

{
  if (isBounding(v))
    return "bounding vertex " + std::to_string(v);
  return "point " + std::to_string(toCaller(v));
}

Referenced by checkEmptyCircumcircle(), and checkInvariants().

xy()

const double * XYIncrementalDelaunay::xy ( std::size_t  v ) const

inlineprivate

Returns

The two coordinates of an internal vertex, in the layout the exact predicates take

Definition at line 155 of file XYIncrementalDelaunay.h.

{ return &_vertices[v].x; }

Referenced by checkEmptyCircumcircle(), checkInvariants(), containsPoint(), growCavity(), insertPoint(), insertSegment(), isStrictlyBetween(), locate(), and triangulatePseudopolygon().

Member Data Documentation

_bounding_reach

constexpr double XYIncrementalDelaunay::_bounding_reach = 1000.0

staticconstexprprivate

How far the bounding triangle reaches beyond the points, as a multiple of their extent.

Definition at line 152 of file XYIncrementalDelaunay.h.

Referenced by initialize().

_constraints

std::set<Segment> XYIncrementalDelaunay::_constraints

private

The constrained segments, in caller vertex ids with the smaller id first.

Definition at line 243 of file XYIncrementalDelaunay.h.

Referenced by checkInvariants(), constrainedSegments(), initialize(), insertPoint(), insertSegment(), isConstrainedEdge(), and isConstrainedSegment().

_last_triangle

std::size_t XYIncrementalDelaunay::_last_triangle = invalid_index

private

A triangle the last insertion produced, which is where the next point walk starts.

Definition at line 246 of file XYIncrementalDelaunay.h.

Referenced by initialize(), and locate().

_num_bounding

constexpr std::size_t XYIncrementalDelaunay::_num_bounding = 3

staticconstexprprivate

The number of bounding triangle vertices padding the front of the vertex list.

Definition at line 149 of file XYIncrementalDelaunay.h.

Referenced by initialize(), isBounding(), numPoints(), toCaller(), and toInternal().

_triangles

std::vector<Triangle> XYIncrementalDelaunay::_triangles

private

The triangles, every one of them live and counter-clockwise.

Definition at line 237 of file XYIncrementalDelaunay.h.

Referenced by checkEmptyCircumcircle(), checkInvariants(), containsPoint(), getTriangles(), growCavity(), initialize(), insertPoint(), insertSegment(), localEdgeIndex(), localVertexIndex(), locate(), and retriangulate().

_vertex_triangle

std::vector<std::size_t> XYIncrementalDelaunay::_vertex_triangle

private

One triangle touching each vertex of _vertices, which is where a walk around it starts.

Definition at line 240 of file XYIncrementalDelaunay.h.

Referenced by initialize(), insertPoint(), insertSegment(), and retriangulate().

_vertices

std::vector<Point2D> XYIncrementalDelaunay::_vertices

private

The bounding triangle vertices followed by the caller's points.

Definition at line 234 of file XYIncrementalDelaunay.h.

Referenced by checkEmptyCircumcircle(), growCavity(), initialize(), insertPoint(), isStrictlyBetween(), numPoints(), point(), and xy().

invalid_index

constexpr std::size_t XYIncrementalDelaunay::invalid_index = std::numeric_limits<std::size_t>::max()

staticconstexpr

Sentinel for a vertex, triangle or neighbor that does not exist.

Definition at line 49 of file XYIncrementalDelaunay.h.

Referenced by XYFrontalDelaunayGenerator::buildTriangleMesh(), checkInvariants(), XYFrontalDelaunayGenerator::collectFront(), getTriangles(), growCavity(), initialize(), insertPoint(), insertSegment(), locate(), and retriangulate().

---

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

• include/utils/XYIncrementalDelaunay.h
• src/utils/XYIncrementalDelaunay.C