Skip to content

GitLab

Commit 1f394810 authored by JasonPries's avatar JasonPries
Browse files

Mesh infrastructure updates in preparation for implementation of Mesh boundary mapping constraints 

Adds a queueing mechanism to insert nodes in order to allow multiple simultaneous insertions.
Adds a BoundaryConstraint vector. Each BoundaryConstraint is associated with one curve, and contains and an index vector for referencing the associated DartConstraints.
parent 46ec5042
Hide whitespace changes
Inline Side-by-side
Showing with 425 additions and 185 deletions
src/Mesh/src/Edge.h
View file @ 1f394810
......@@ -5,7 +5,7 @@
class Mesh;
class Edge {
class Edge { // TODO: Should be renamed Dart
public:
friend class Mesh;
......
src/Mesh/src/Mesh.cpp
View file @ 1f394810
......@@ -14,7 +14,7 @@ Mesh::Mesh(Sketch &sketch) {
Contours.push_back(sketch.contour(i));
}
Constraints.push_back(DartConstraint());
new_dart_constraint(0.0, 1.0, std::make_shared<BoundaryConstraint>(nullptr));
}
bool Mesh::are_intersecting(size_t ei, size_t ej) const {
......@@ -124,32 +124,32 @@ bool Mesh::edges_are_valid() const {
}
if (is_constrained(e)) {
DartConstraint dc = Constraints[Edges[e].Constraint];
DartConstraint dc = DartConstraints[Edges[e].Constraint];
double tol = length(e) * FLT_EPSILON;
if (orientation(e)) {
Point p0 = base(e);
Point p1 = dc.ConstraintCurve->point(dc.S0);
Point p1 = dc.constraint_curve()->point(dc.S0);
if (dist(p0,p1) > tol) {
result = false;
break;
}
p0 = tip(e);
p1 = dc.ConstraintCurve->point(dc.S1);
p1 = dc.constraint_curve()->point(dc.S1);
if (dist(p0,p1) > tol) {
result = false;
break;
}
} else {
Point p0 = base(e);
Point p1 = dc.ConstraintCurve->point(dc.S1);
Point p1 = dc.constraint_curve()->point(dc.S1);
if (dist(p0,p1) > tol) {
result = false;
break;
}
p0 = tip(e);
p1 = dc.ConstraintCurve->point(dc.S0);
p1 = dc.constraint_curve()->point(dc.S0);
if (dist(p0,p1) > tol) {
result = false;
break;
......@@ -552,10 +552,12 @@ void Mesh::create_boundary_polygon() {
for (size_t i = 0; i != Boundary->size(); ++i) {
std::shared_ptr<Curve const> cc = Boundary->curve(i);
bool dir = Boundary->orientation(i);
Edge &e = new_edge(Points.size(), Constraints.size(), dir); // TODO: change to new_edge(Edge)
Edge &e = new_edge(Points.size(), DartConstraints.size(), dir); // TODO: change to new_edge(Edge)
e.Twin = e.Self;
Constraints.push_back(DartConstraint(0.0, 1.0, cc));
BoundaryConstraints.push_back(std::make_shared<BoundaryConstraint>(cc));
new_dart_constraint(0.0, 1.0, BoundaryConstraints.back());
Points.push_back(Point(dir ? cc->start() : cc->end()));
}
......@@ -658,25 +660,33 @@ void Mesh::insert_internal_boundaries() {
Point p = interior[i]->start();
LocateTriangleResult result = locate_triangle(p);
if (result == LocateTriangleResult::Interior) {
while (insert_point(p) == InsertPointResult::Midpoint);
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
// Insert end point
p = interior[i]->end();
result = locate_triangle(p);
if (result == LocateTriangleResult::Interior) {
while (insert_point(p) == InsertPointResult::Midpoint);
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
}
// Insert interior curve midpoints until constraints are naturally satisfied
std::vector<DartConstraint> queue;
// Attach edges to constraints by inserting interior curve midpoints until constraints are naturally satisfied
std::vector<size_t> queue;
for (size_t i = 0; i != interior.size(); ++i) {
queue.push_back(DartConstraint(0.0, 1.0, interior[i]));
BoundaryConstraints.push_back(std::make_shared<BoundaryConstraint>(interior[i]));
queue.push_back(DartConstraints.size());
new_dart_constraint(0.0, 1.0, BoundaryConstraints.back());
while (queue.size() != 0) {
DartConstraint &dc = queue.back();
Point p0 = dc.ConstraintCurve->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
Point p1 = dc.ConstraintCurve->point(dc.S1);
DartConstraint &dc = DartConstraints[queue.back()];
Point p0 = dc.constraint_curve()->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
Point p1 = dc.constraint_curve()->point(dc.S1);
size_t ei = Edges.size() - 1;
LocateTriangleResult result = locate_triangle(p0, ei);
......@@ -687,31 +697,38 @@ void Mesh::insert_internal_boundaries() {
if (find_attached(p1, ei)) {
Edge &e = Edges[ei];
e.Constraint = Constraints.size();
e.Constraint = dc.Self;
e.Orientation = true;
Edge &et = Edges[e.Twin];
et.Constraint = Constraints.size();
et.Constraint = dc.Self;
et.Orientation = false;
Constraints.push_back(queue.back());
queue.pop_back();
} else {
double s0 = dc.S0;
double s1 = dc.S1;
std::shared_ptr<Curve const> cc = dc.ConstraintCurve;
std::shared_ptr<Curve const> cc = dc.constraint_curve();
double sn = (s0 + s1) / 2.0;
Point const p = cc->point(sn);
dc.S1 = sn;
queue.push_back(DartConstraint(sn, s1, cc));
while (insert_point(p) == InsertPointResult::Midpoint);
queue.push_back(DartConstraints.size());
new_dart_constraint(sn, s1, dc.boundary_constraint());
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
}
}
// TODO: One final pass to insert any midpoints that have not been inserted OR length based refinement
// There are some errors that occur with half-circles becoming degenerate when the LineSegment is not encroached because this isn't done
split_encroached_edges();
}
......@@ -732,12 +749,20 @@ void Mesh::refine_once(std::vector<size_t> index, std::vector<double> radii, std
for (size_t i = 0; i != Triangles.size(); ++i) {
size_t j = index[i];
if ((triangle(j).Mark) && ((radii[j] > MaximumElementSize) || (radii[j] > MinimumElementSize && quality[j] < MinimumElementQuality))) {
insert_circumcenter(Triangles[j]);
add_circumcenter_to_queue(Triangles[j]);
insert_from_queue();
}
}
get_triangles();
}
void Mesh::insert_from_queue() {
while (!Queue.empty()) {
Queue.back()->insert(*this);
Queue.pop_back();
}
}
void Mesh::save_as(std::string path, std::string file_name) const {
/*
This is a stub for visualization
......@@ -778,58 +803,52 @@ void Mesh::split_edge(size_t ei) {
Used for initial polygon refinement.
*/
// TODO: Refactor into split_constrained_edge
size_t c{0};
if (is_constrained(ei)) {
c = Constraints.size();
DartConstraint &dc = Constraints[Edges[ei].Constraint];
double s0 = dc.S0;
double s1 = dc.S1;
std::shared_ptr<Curve const> cc = dc.ConstraintCurve;
double sn = (s0 + s1) / 2.0;
dc.S1 = sn;
Constraints.push_back(DartConstraint(sn, s1, cc));
Points.push_back(cc->point(sn));
} else { // Unconstrained Edge
Point const p0 = base(ei);
Point const p1 = tip(ei);
// TODO: Rename to split_boundary_edge
Points.push_back(Point((p0.X + p1.X) / 2.0, (p0.Y + p1.Y) / 2.0));
if (!(is_constrained(ei) && ei == twin(ei))) { // TODO: write bool Edge::is_boundary()
std::cerr << "Mesh::split_edge(size_t ei) should only be called on boundary edges" << std::endl;
return;
}
if (ei == twin(ei)) { // Boundary Edge, TODO: write bool Edge::is_boundary()
size_t itr = new_edges(1);
Edge &newe = Edges[--itr];
Edge &e = Edges[ei];
Edge &nxt = Edges[e.Next];
size_t c{DartConstraints.size()};
// Constraint Curve
newe.Orientation = e.Orientation;
if (e.Orientation) {
newe.Constraint = c;
} else {
newe.Constraint = e.Constraint;
e.Constraint = c;
}
DartConstraint &dci = DartConstraints[Edges[ei].Constraint];
double s0 = dci.S0;
double s1 = dci.S1;
std::shared_ptr<Curve const> cc = dci.constraint_curve();
// Connectivity
newe.Node = Points.size() - 1;
newe.Next = e.Next;
newe.Prev = e.Self;
newe.Twin = newe.Self;
newe.Mark = false;
double sn = (s0 + s1) / 2.0;
dci.S1 = sn;
nxt.Prev = newe.Self;
nxt.Mark = false;
new_dart_constraint(sn, s1, dci.boundary_constraint());
Points.push_back(cc->point(sn));
e.Next = newe.Self;
e.Mark = false;
size_t itr = new_edges(1);
Edge &newe = Edges[--itr];
Edge &e = Edges[ei];
Edge &nxt = Edges[e.Next];
// Constraint Curve
newe.Orientation = e.Orientation;
if (e.Orientation) {
newe.Constraint = c;
} else {
throw std::exception(); // Function should only be called on constrained edges
newe.Constraint = e.Constraint;
e.Constraint = c;
}
// Connectivity
newe.Node = Points.size() - 1;
newe.Next = e.Next;
newe.Prev = e.Self;
newe.Twin = newe.Self;
newe.Mark = false;
nxt.Prev = newe.Self;
nxt.Mark = false;
e.Next = newe.Self;
e.Mark = false;
}
void Mesh::split_encroached_edges() {
......@@ -1061,8 +1080,8 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
}
}
InsertPointResult Mesh::insert_circumcenter(size_t ei) {
return insert_point(circumcenter(ei), ei);
AddToQueueResult Mesh::add_circumcenter_to_queue(size_t dart) {
return add_point_to_queue(circumcenter(dart), dart);
}
InsertPointResult Mesh::insert_midpoint(size_t ei) {
......@@ -1072,17 +1091,17 @@ InsertPointResult Mesh::insert_midpoint(size_t ei) {
size_t c{0};
if (is_constrained(ei)) { // Constrained Edge
c = Constraints.size();
c = DartConstraints.size();
DartConstraint &dc = Constraints[Edges[ei].Constraint];
DartConstraint &dc = DartConstraints[Edges[ei].Constraint];
double s0 = dc.S0;
double s1 = dc.S1;
std::shared_ptr<Curve const> cc = dc.ConstraintCurve;
std::shared_ptr<Curve const> cc = dc.constraint_curve();
double sn = (s0 + s1) / 2.0;
dc.S1 = sn;
Constraints.push_back(DartConstraint(sn, s1, cc));
new_dart_constraint(sn, s1, dc.boundary_constraint());
Points.push_back(cc->point(sn));
} else { // Unconstrained Edge
Point const p0 = base(ei);
......@@ -1242,118 +1261,138 @@ InsertPointResult Mesh::insert_midpoint(size_t ei) {
return InsertPointResult::Midpoint;
}
InsertPointResult Mesh::insert_point(Point const vc, size_t ei) {
AddToQueueResult Mesh::add_point_to_queue(Point const vc, size_t ei) {
// Find triangle containing point
LocateTriangleResult result = locate_triangle(vc, ei);
// Do not insert point if it is duplicate
if (result == LocateTriangleResult::Point) {
return InsertPointResult::Duplicate;
return AddToQueueResult::Duplicate;
}
// Test edges in current and adjacent triangles for encroachment
// These are the only possible edges that are encroached due to empty circumcircle property
std::vector<size_t> test_edges;
test_edges.reserve(9);
std::vector<size_t> encroached_edges = get_encroached_edges(vc, ei);
bool any_encroached{encroached_edges.size() > 0};
test_edges.push_back(ei);
test_edges.push_back(next(ei));
test_edges.push_back(prev(ei));
for (size_t i = 0; i != 3; ++i) {
Edge e = Edges[test_edges[i]];
if (e.Self != e.Twin) {
test_edges.push_back(twin(e).Next);
test_edges.push_back(twin(e).Prev);
if (any_encroached) {
for (size_t e : encroached_edges) {
Queue.push_back(std::make_unique<MidpointQueuer>(e));
}
}
return AddToQueueResult::Midpoint;
// Split all encroached edges
bool encroached = false;
for (auto e : test_edges) {
if (is_encroached(vc, e)) {
insert_midpoint(e);
encroached = true;
}
}
// If none are encroached, insert circumcenter
if (encroached) {
return InsertPointResult::Midpoint;
} else if (result == LocateTriangleResult::Interior) {
size_t itr = new_edges(6);
Edge &e5 = Edges[--itr];
Edge &e4 = Edges[--itr];
Edge &e3 = Edges[--itr];
Edge &e2 = Edges[--itr];
Edge &e1 = Edges[--itr];
Edge &e0 = Edges[--itr];
Edge &tri = Edges[ei];
Edge &nxt = Edges[tri.Next];
Edge &prv = Edges[tri.Prev];
size_t vt = node(tri);
size_t vn = node(nxt);
size_t vp = node(prv);
Points.push_back(vc);
e0.Node = Points.size() - 1;
e0.Next = tri.Self;
e0.Prev = e1.Self;
e0.Twin = e5.Self;
e0.Mark = false;
e1.Node = vn;
e1.Next = e0.Self;
e1.Prev = tri.Self;
e1.Twin = e2.Self;
e1.Mark = false;
e2.Node = Points.size() - 1;
e2.Next = nxt.Self;
e2.Prev = e3.Self;
e2.Twin = e1.Self;
e2.Mark = false;
e3.Node = vp;
e3.Next = e2.Self;
e3.Prev = nxt.Self;
e3.Twin = e4.Self;
e3.Mark = false;
e4.Node = Points.size() - 1;
e4.Next = prv.Self;
e4.Prev = e5.Self;
e4.Twin = e3.Self;
e4.Mark = false;
Queue.push_back(std::make_unique<CircumcenterQueuer>(vc, ei));
return AddToQueueResult::Circumcenter;
e5.Node = vt;
e5.Next = e4.Self;
e5.Prev = prv.Self;
e5.Twin = e0.Self;
e5.Mark = false;
} else {
throw std::exception(); // TODO: No triangle could be found containing circumcenter and no boundary edge was encroached by circumcenter
}
}
nxt.Next = e3.Self;
nxt.Prev = e2.Self;
nxt.Mark = false;
InsertPointResult Mesh::insert_point(Point const vc, size_t ei) {
// TODO: Give Queuers function pointers and references to the mesh, then make these private
// ei should be the triangle containing the point vc
size_t itr = new_edges(6);
Edge &e5 = Edges[--itr];
Edge &e4 = Edges[--itr];
Edge &e3 = Edges[--itr];
Edge &e2 = Edges[--itr];
Edge &e1 = Edges[--itr];
Edge &e0 = Edges[--itr];
Edge &tri = Edges[ei];
Edge &nxt = Edges[tri.Next];
Edge &prv = Edges[tri.Prev];
size_t vt = node(tri);
size_t vn = node(nxt);
size_t vp = node(prv);
Points.push_back(vc);
e0.Node = Points.size() - 1;
e0.Next = tri.Self;
e0.Prev = e1.Self;
e0.Twin = e5.Self;
e0.Mark = false;
e1.Node = vn;
e1.Next = e0.Self;
e1.Prev = tri.Self;
e1.Twin = e2.Self;
e1.Mark = false;
e2.Node = Points.size() - 1;
e2.Next = nxt.Self;
e2.Prev = e3.Self;
e2.Twin = e1.Self;
e2.Mark = false;
e3.Node = vp;
e3.Next = e2.Self;
e3.Prev = nxt.Self;
e3.Twin = e4.Self;
e3.Mark = false;
e4.Node = Points.size() - 1;
e4.Next = prv.Self;
e4.Prev = e5.Self;
e4.Twin = e3.Self;
e4.Mark = false;
e5.Node = vt;
e5.Next = e4.Self;
e5.Prev = prv.Self;
e5.Twin = e0.Self;
e5.Mark = false;
nxt.Next = e3.Self;
nxt.Prev = e2.Self;
nxt.Mark = false;
prv.Next = e5.Self;
prv.Prev = e4.Self;
prv.Mark = false;
tri.Next = e1.Self;
tri.Prev = e0.Self;
tri.Mark = false;
recursive_swap(tri.Self);
recursive_swap(nxt.Self);
recursive_swap(prv.Self);
return InsertPointResult::Success;
}
prv.Next = e5.Self;
prv.Prev = e4.Self;
prv.Mark = false;
std::vector<size_t> Mesh::get_encroached_edges(Point const p, size_t edge) {
std::vector<size_t> encroached_edges;
encroached_edges.reserve(9);
size_t e_self = edge;
for (size_t i = 0; i != 3; ++i) {
if (is_encroached(p, e_self)) {
encroached_edges.push_back(e_self);
}
tri.Next = e1.Self;
tri.Prev = e0.Self;
tri.Mark = false;
size_t e_twin = twin(e_self);
if (e_self != e_twin) { // if not boundary edge, check adjacent triangle
size_t e_next = next(e_twin);
if (is_encroached(p, e_next)) {
encroached_edges.push_back(e_next);
}
recursive_swap(tri.Self);
recursive_swap(nxt.Self);
recursive_swap(prv.Self);
size_t e_prev = prev(e_twin);
if (is_encroached(p, e_prev)) {
encroached_edges.push_back(e_prev);
}
}
return InsertPointResult::Success;
} else {
throw std::exception(); // TODO: No triangle could be found containing circumcenter and no boundary edge was encroached by circumcenter
e_self = next(e_self);
}
return encroached_edges;
}
Point Mesh::circumcenter(size_t ei) const {
......
src/Mesh/src/Mesh.h
View file @ 1f394810
......@@ -8,6 +8,7 @@
#include <algorithm>
#include <fstream>
#include <numeric>
#include <set>
enum class LocateTriangleResult {
Interior, Exterior, Boundary, Point // TODO: Enumerate cases in function when triangle is located by a point near the boundary
......@@ -17,25 +18,122 @@ enum class InsertPointResult {
Success, Midpoint, Duplicate, Failure
};
enum class AddToQueueResult {
Circumcenter, Midpoint, Duplicate
};
class PointQueuer {
public:
virtual ~PointQueuer() {};
virtual void insert(Mesh &m) const = 0;
};
class BoundaryConstraint {
public:
BoundaryConstraint(std::shared_ptr<Curve const> cc) : ConstraintCurve(cc) {};
void add_dart(size_t d) { Darts.push_back(d); }; // TODO: Enforce uniqueness (std::set?)
std::shared_ptr<Curve const> curve() const { return ConstraintCurve; };
protected:
std::vector<size_t> Darts;
std::shared_ptr<Curve const> ConstraintCurve;
};
class CircumcenterQueuer;
class MidpointQueuer;
class DartConstraint {
public:
DartConstraint() : S0(DBL_MAX), S1(DBL_MAX), ConstraintCurve(nullptr) {};
DartConstraint(double s0, double s1, std::shared_ptr<BoundaryConstraint> bc, size_t self) : S0{s0}, S1{s1}, Constraint{bc}, Self{self} {};
DartConstraint(double s0, double s1, std::shared_ptr<Curve const> cc) : S0(s0), S1(s1), ConstraintCurve(cc) {};
std::shared_ptr<BoundaryConstraint> boundary_constraint() const { return Constraint; };
std::shared_ptr<Curve const> constraint_curve() const { return Constraint->curve(); };
double S0;
double S1;
size_t Self;
protected:
std::shared_ptr<BoundaryConstraint> Constraint;
};
/* TODO: Reverse relationship between Edge and DartConstraint
* TODO: Have DartConstraint point to the Edge it constrains
* TODO: Edge will only contain a bool IsConstrained value in order to quickly test for encroachment
* TODO:
class BoundaryConstraint {
public:
BoundaryConstraint() : ConstraintCurve(nullptr) {};
BoundaryConstraint(std::shared_ptr<Curve const> cc) : ConstraintCurve(cc) {};
std::set<size_t> DartConstraints; // TODO: Or should Darts/(Edges) be stored instead of DartConstrants? Probably not because then Darts/Edges would have to hold more data
std::shared_ptr<Curve const> ConstraintCurve;