Commit 20229513 authored by JasonPries's avatar JasonPries
Browse files

Basic implementation of optional uniform boundary discretizations

This is groundwork for the initial implementation strategy for virtual motion with node renumbering and periodic boundary constraints in the FEA solvers
parent 1865f2b1
...@@ -2,61 +2,83 @@ ...@@ -2,61 +2,83 @@
#include "DartConstraint.h" #include "DartConstraint.h"
#include "Mesh.h" #include "Mesh.h"
void BoundaryConstraint::add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints, size_t dart_constraint) { bool BoundaryConstraint::is_uniform(Mesh const &m) const {
if (UniformDiscretization) { DartConstraint const &dc = m.dart_constraint(DartConstraints[0]);
for (size_t dc : DartConstraints) { double_t delta_0{dc.delta()};
size_t dart = constraints[dc].dart();
queue.push_back(std::make_unique<MidpointQueuer>(dart)); for(size_t i : DartConstraints) {
double delta_i{m.dart_constraint(i).delta()};
if (std::abs(delta_0 - delta_i) > FLT_EPSILON) {
return false;
} }
} else {
size_t dart = constraints[dart_constraint].dart();
queue.push_back(std::make_unique<MidpointQueuer>(dart));
} }
return true;
} }
void BoundaryConstraint::refine(Mesh &m, double_t tol) { double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
double_t err_max{DBL_MAX}; double_t delta_min{DBL_MAX};
while (err_max > tol) {
err_max = 0.0;
for (size_t i : DartConstraints) { for (size_t i : DartConstraints) {
DartConstraint const &dc = m.dart_constraint(i); DartConstraint const &dc = m.dart_constraint(i);
delta_min = std::min(delta_min, std::abs(dc.S1 - dc.S0));
}
double_t len_c = ConstraintCurve->length(dc.S0, dc.S1); return delta_min;
double_t len_e = m.edge(dc.dart()).length(m); }
double_t err_i = std::abs(len_c - len_e) / len_c;
if (err_i > tol) { double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
double delta_max{delta_min};
for(size_t i : DartConstraints) {
DartConstraint const &dc = m.dart_constraint(i);
double_t delta = std::abs(dc.S1 - dc.S0);
if (delta > delta_min) {
m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart())); m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
err_max = std::max(err_max, err_i); delta_max = std::max(delta_max, delta);
} }
} }
m.insert_from_queue(); return delta_max;
}
void BoundaryConstraint::add_to_queue(Mesh &m, size_t dart) const {
if (UniformDiscretization) {
queue_uniform(m, smallest_parametric_edge_length(m) / 2.0);
} else {
m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
} }
} }
void BoundaryConstraint::make_uniform(Mesh &m) {
double_t delta_min{DBL_MAX};
for (size_t i : DartConstraints) { void BoundaryConstraint::make_uniform(Mesh &m) { make_uniform(m, smallest_parametric_edge_length(m)); };
DartConstraint const &dc = m.dart_constraint(i);
delta_min = std::min(delta_min, std::abs(dc.S1 - dc.S0));
}
delta_min *= (1 + FLT_EPSILON);
void BoundaryConstraint::make_uniform(Mesh &m, double delta_min) {
delta_min *= (1 + FLT_EPSILON);
double_t delta_max{DBL_MAX}; double_t delta_max{DBL_MAX};
while (delta_max > delta_min) { while (delta_max > delta_min) {
delta_max = delta_min; delta_max = queue_uniform(m, delta_min);
m.insert_from_queue();
}
}
for(size_t i : DartConstraints) { void BoundaryConstraint::refine(Mesh &m, double_t tol) {
double_t err_max{DBL_MAX};
while (err_max > tol) {
err_max = 0.0;
for (size_t i : DartConstraints) {
DartConstraint const &dc = m.dart_constraint(i); DartConstraint const &dc = m.dart_constraint(i);
double_t delta = std::abs(dc.S1 - dc.S0);
if (delta > delta_min) { double_t len_c = ConstraintCurve->length(dc.S0, dc.S1);
double_t len_e = m.edge(dc.dart()).length(m);
double_t err_i = std::abs(len_c - len_e) / len_c;
if (err_i > tol) {
m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart())); m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
delta_max = std::max(delta_max, delta); err_max = std::max(err_max, err_i);
} }
} }
......
...@@ -14,23 +14,31 @@ class BoundaryConstraint { ...@@ -14,23 +14,31 @@ class BoundaryConstraint {
public: public:
BoundaryConstraint(std::shared_ptr<Curve const> cc) : ConstraintCurve(cc) {}; BoundaryConstraint(std::shared_ptr<Curve const> cc) : ConstraintCurve(cc) {};
void add_dart_constraint(size_t d) { DartConstraints.push_back(d); }; // TODO: Enforce uniqueness (std::set?) size_t size_dart_constraints() const { return DartConstraints.size(); };
std::shared_ptr<Curve const> curve() const { return ConstraintCurve; }; size_t dart(size_t i) const { return DartConstraints[i]; };
bool is_uniform(Mesh const &m) const;
bool uniform_discretizaiton() const { return UniformDiscretization; }; bool uniform_discretizaiton() const { return UniformDiscretization; };
double smallest_parametric_edge_length(Mesh &m) const;
double queue_uniform(Mesh &m, double delta_min) const;
void add_dart_constraint(size_t d) { DartConstraints.push_back(d); }; // TODO: Enforce uniqueness (std::set?)
void uniform_discretization(bool flag) { UniformDiscretization = flag; }; void uniform_discretization(bool flag) { UniformDiscretization = flag; };
void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints, size_t dart); void add_to_queue(Mesh &m, size_t dart) const;
void refine(Mesh &m, double_t tol); void refine(Mesh &m, double_t tol);
void make_uniform(Mesh &m); void make_uniform(Mesh &m);
size_t size_dart_constraints() const { return DartConstraints.size(); }; void make_uniform(Mesh &m, double delta_min);
size_t dart(size_t i) const { return DartConstraints[i]; }; std::shared_ptr<Curve const> curve() const { return ConstraintCurve; };
protected: protected:
std::vector<size_t> DartConstraints; std::vector<size_t> DartConstraints;
......
...@@ -24,9 +24,9 @@ public: ...@@ -24,9 +24,9 @@ public:
size_t reverse_dart() const { return ReverseDart; }; size_t reverse_dart() const { return ReverseDart; };
void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints) { double delta() const { return std::abs(S1 - S0); };
Constraint->add_to_queue(queue, constraints, Self);
} void add_to_queue(Mesh &m) const { Constraint->add_to_queue(m, dart()); };
void forward_dart(size_t d) { ForwardDart = d; }; void forward_dart(size_t d) { ForwardDart = d; };
......
...@@ -4,3 +4,7 @@ ...@@ -4,3 +4,7 @@
double_t Edge::length(Mesh const &m) const { double_t Edge::length(Mesh const &m) const {
return dist(m.point(node()), m.point(m.edge(next()).node())); return dist(m.point(node()), m.point(m.edge(next()).node()));
} }
void Edge::add_to_queue(Mesh &m) {
m.constraint(Self).add_to_queue(m);
}
\ No newline at end of file
...@@ -47,9 +47,7 @@ public: ...@@ -47,9 +47,7 @@ public:
double_t length(Mesh const &m) const; double_t length(Mesh const &m) const;
void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints) { void add_to_queue(Mesh &m);
constraints[Constraint].add_to_queue(queue, constraints);
}
protected: protected:
size_t Node; //Point at start of this edge size_t Node; //Point at start of this edge
......
...@@ -7,8 +7,10 @@ Mesh::Mesh(Sketch &sketch) { ...@@ -7,8 +7,10 @@ Mesh::Mesh(Sketch &sketch) {
auto c = sketch.curve(i); auto c = sketch.curve(i);
if (!(c->for_construction())) { if (!(c->for_construction())) {
Curves.push_back(c); Curves.push_back(c);
BoundaryConstraints.push_back(std::make_shared<BoundaryConstraint>(c));
} }
} }
sort_constraints();
for (size_t i = 0; i != sketch.size_contours(); ++i) { for (size_t i = 0; i != sketch.size_contours(); ++i) {
Contours.push_back(sketch.contour(i)); Contours.push_back(sketch.contour(i));
...@@ -555,8 +557,7 @@ void Mesh::create_boundary_polygon() { ...@@ -555,8 +557,7 @@ void Mesh::create_boundary_polygon() {
Edge &e = new_edge(Points.size(), DartConstraints.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; e.Twin = e.Self;
BoundaryConstraints.push_back(std::make_shared<BoundaryConstraint>(cc)); DartConstraint &dc = new_dart_constraint(0.0, 1.0, boundary_constraint(cc));
DartConstraint &dc = new_dart_constraint(0.0, 1.0, BoundaryConstraints.back());
if (dir) { if (dir) {
dc.forward_dart(e.self()); dc.forward_dart(e.self());
...@@ -575,8 +576,6 @@ void Mesh::create_boundary_polygon() { ...@@ -575,8 +576,6 @@ void Mesh::create_boundary_polygon() {
Edges[j].Prev = Edges[i].Self; Edges[j].Prev = Edges[i].Self;
} }
// TODO: Length Based Refinement of Boundary curves
// Some edges may intersect due to discretization error // Some edges may intersect due to discretization error
// If two edges intersect, split the longest edge // If two edges intersect, split the longest edge
bool any_split = true; bool any_split = true;
...@@ -686,9 +685,8 @@ void Mesh::insert_internal_boundaries() { ...@@ -686,9 +685,8 @@ void Mesh::insert_internal_boundaries() {
// Attach edges to constraints by inserting interior curve midpoints until constraints are naturally satisfied // Attach edges to constraints by inserting interior curve midpoints until constraints are naturally satisfied
std::vector<size_t> queue; std::vector<size_t> queue;
for (size_t i = 0; i != interior.size(); ++i) { for (size_t i = 0; i != interior.size(); ++i) {
BoundaryConstraints.push_back(std::make_shared<BoundaryConstraint>(interior[i]));
queue.push_back(DartConstraints.size()); queue.push_back(DartConstraints.size());
new_dart_constraint(0.0, 1.0, BoundaryConstraints.back()); new_dart_constraint(0.0, 1.0, boundary_constraint(interior[i]));
while (queue.size() != 0) { while (queue.size() != 0) {
DartConstraint &dc = DartConstraints[queue.back()]; 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 p0 = dc.constraint_curve()->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
...@@ -735,26 +733,17 @@ void Mesh::insert_internal_boundaries() { ...@@ -735,26 +733,17 @@ void Mesh::insert_internal_boundaries() {
} }
} }
// TODO: Do length based refinement instead for (auto bc : BoundaryConstraints) { // 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
for (auto bc : BoundaryConstraints) {
bc->refine(*this, 0.01); bc->refine(*this, 0.01);
/*
auto &temp = *(bc.get());
if(bc->size_dart_constraints() == 1) {
bc->add_to_queue(Queue, DartConstraints, bc->dart(0));
insert_from_queue();
}
*/
} }
/* split_encroached_edges();
for (auto bc : BoundaryConstraints) {
for (auto bc: BoundaryConstraints) {
if (bc->uniform_discretizaiton()) {
bc->make_uniform(*this); bc->make_uniform(*this);
} }
*/ }
split_encroached_edges();
// TODO: Enforce BoundaryConstraint conditions (e.g. uniform discretization, boundary maps) // TODO: Enforce BoundaryConstraint conditions (e.g. uniform discretization, boundary maps)
} }
...@@ -812,6 +801,14 @@ void Mesh::save_as(std::string path, std::string file_name) const { ...@@ -812,6 +801,14 @@ void Mesh::save_as(std::string path, std::string file_name) const {
fs.close(); fs.close();
} }
void Mesh::sort_constraints() {
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(),
[](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) {
return (size_t)(x->curve().get()) < (size_t)(y->curve().get());
}
);
}
void Mesh::sort_permutation_ascending(std::vector<double> &values, std::vector<size_t> &index) const { void Mesh::sort_permutation_ascending(std::vector<double> &values, std::vector<size_t> &index) const {
index.resize(values.size()); index.resize(values.size());
std::iota(index.begin(), index.end(), 0); std::iota(index.begin(), index.end(), 0);
...@@ -947,6 +944,15 @@ void Mesh::triangulate_boundary_polygon() { ...@@ -947,6 +944,15 @@ void Mesh::triangulate_boundary_polygon() {
} }
split_encroached_edges(); split_encroached_edges();
for (size_t i = 0; i != Boundary->size(); ++i) { // Length based-refinement of boundary curves
std::shared_ptr<Curve const> cc = Boundary->curve(i);
std::shared_ptr<BoundaryConstraint> bc = boundary_constraint(cc);
bc->refine(*this, 0.01);
if (bc->uniform_discretizaiton()) {
bc->make_uniform(*this);
}
}
} }
LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const { LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
...@@ -1344,8 +1350,7 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) { ...@@ -1344,8 +1350,7 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) {
if (any_encroached) { if (any_encroached) {
for (size_t e : encroached_edges) { for (size_t e : encroached_edges) {
//Queue.push_back(std::make_unique<MidpointQueuer>(e)); Edges[e].add_to_queue(*this);
Edges[e].add_to_queue(Queue, DartConstraints);
} }
return AddToQueueResult::Midpoint; return AddToQueueResult::Midpoint;
...@@ -1483,3 +1488,15 @@ Point Mesh::circumcenter(size_t ei) const { ...@@ -1483,3 +1488,15 @@ Point Mesh::circumcenter(size_t ei) const {
return Point{x, y}; return Point{x, y};
} }
std::shared_ptr<BoundaryConstraint> Mesh::boundary_constraint(std::shared_ptr<Curve const> const &c) const {
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, size_t const y) { return (size_t)(x->curve().get()) < y; };
std::shared_ptr<BoundaryConstraint> result = *std::lower_bound(BoundaryConstraints.begin(), BoundaryConstraints.end(), (size_t)c.get(), comp);
if((size_t)(result->curve().get()) == (size_t)(c.get())) {
return result;
} else {
return nullptr;
}
}
\ No newline at end of file
...@@ -21,68 +21,6 @@ enum class InsertPointResult { ...@@ -21,68 +21,6 @@ enum class InsertPointResult {
Success, Midpoint, Duplicate, Failure Success, Midpoint, Duplicate, Failure
}; };
/* 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;
std::set<size_t> RefinementQueue;
void add_to_queue(size_t d) {
// TODO: Add dart constraints to queue for any mapped boundaries
if (DartConstraint.count(d) == 1) {
if (Uniform) {
RefinementQueue = DartConstraints;
} else {
RefinementQueue.insert(d);
}
} else { // Don't do that, give warning
throw;
}
}
bool insert(Mesh &m) {
for(size_t d : RefinementQueue) {
Mesh.insert_midpoint(Mesh.DartConstraints[d].DartForward);
}
RefinementQueue = std::set<size_t>(); // empty
};
// TODO: Have two seperate methods
// TODO: One to add darts to a refinement queue (set with unique edges)
// TODO: An insert method which empties the refinement queue
// TODO: This is necessary because there is a possibility with encroached edges and boundary mappings that too many refinements will occur
//something like std::vector<BoundaryMapping> Map to control, e.g. periodic boundary mapping
//bool Uniform; to control boundary discretizations which are made of evenly spaced points in the parameter space
};
class DartConstraint {
public:
DartConstraint() : S0(DBL_MAX), S1(DBL_MAX), BoundaryConstraints(SIZE_MAX), DartForward(SIZE_MAX), DartReverse(SIZE_MAX) {};
DartConstraint(double s0, double s1) : DartForward(d), S0(s0), S1(s1) {};
void add_to_queue(Mesh &m) { m.boundary_constraint(BoundaryConstraint).add_to_queue(m, Self)}; // TODO: Have dart refinement queue in Mesh class???
double S0;
double S1;
size_t Self;
size_t BoundaryConstraint;
size_t DartForward; // TODO: Register Edge with DartConstraint on assignment
size_t DartReverse; // TODO: Register Edge with DartConstraint on assignment
}
*/
class Mesh { class Mesh {
public: public:
Mesh(Sketch &s); Mesh(Sketch &s);
...@@ -153,13 +91,16 @@ public: ...@@ -153,13 +91,16 @@ public:
void save_as(std::string path, std::string file_name) const; void save_as(std::string path, std::string file_name) const;
// TODO: rename these methods something like, dart_constraint_from_edge_index, curve_from_edge_index
DartConstraint const constraint(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; }; DartConstraint const constraint(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; };
std::shared_ptr<Curve const> constraint_curve(size_t ei) const { return DartConstraints[Edges[ei].Constraint].constraint_curve(); }; std::shared_ptr<Curve const> constraint_curve(size_t ei) const { return DartConstraints[Edges[ei].Constraint].constraint_curve(); };
std::shared_ptr<BoundaryConstraint> boundary_constraint(std::shared_ptr<Curve const> const &c) const;
Point circumcenter(size_t ei) const; Point circumcenter(size_t ei) const;
DartConstraint const dart_constraint(size_t i) { return DartConstraints[i]; }; DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
Point const &base(Edge const &e) const { return Points[e.Node]; }; Point const &base(Edge const &e) const { return Points[e.Node]; };
...@@ -206,6 +147,7 @@ public: ...@@ -206,6 +147,7 @@ public:
AddToQueueResult add_point_to_queue(Point const p) { return add_point_to_queue(p, Edges.size() - 1); }; AddToQueueResult add_point_to_queue(Point const p) { return add_point_to_queue(p, Edges.size() - 1); };
protected: protected:
// TODO: There is some amount of redundancy here (e.g. There is one Curve for each BoundaryConstraint, Countours and Boundary also contain those curves
std::shared_ptr<Contour const> Boundary; std::shared_ptr<Contour const> Boundary;
std::vector<std::shared_ptr<Curve const>> Curves; std::vector<std::shared_ptr<Curve const>> Curves;
std::vector<std::shared_ptr<Contour const>> Contours; std::vector<std::shared_ptr<Contour const>> Contours;
...@@ -264,6 +206,8 @@ private: ...@@ -264,6 +206,8 @@ private:
void refine_once(std::vector<size_t> index, std::vector<double> circumradius, std::vector<double> quality); void refine_once(std::vector<size_t> index, std::vector<double> circumradius, std::vector<double> quality);
void sort_constraints();
void sort_permutation_ascending(std::vector<double> &value, std::vector<size_t> &index) const; void sort_permutation_ascending(std::vector<double> &value, std::vector<size_t> &index) const;
void sort_permutation_descending(std::vector<double> &values, std::vector<size_t> &index) const; void sort_permutation_descending(std::vector<double> &values, std::vector<size_t> &index) const;
......
...@@ -627,7 +627,8 @@ public: ...@@ -627,7 +627,8 @@ public:
auto l4 = s.new_element<LineSegment>(v4, v5); auto l4 = s.new_element<LineSegment>(v4, v5);
auto l5 = s.new_element<LineSegment>(v4, v6); auto l5 = s.new_element<LineSegment>(v4, v6);
auto c0 = s.new_element<CircularArc>(v6, v5, v4, 0.28); c0 = s.new_element<CircularArc>(v6, v5, v4, 0.28);
cc = c0;
auto fixation0 = s.new_element<Fixation>(v4); auto fixation0 = s.new_element<Fixation>(v4);
auto radius0 = s.new_element<Radius>(c0, 0.28); auto radius0 = s.new_element<Radius>(c0, 0.28);
...@@ -636,6 +637,9 @@ public: ...@@ -636,6 +637,9 @@ public:
auto vertical1 = s.new_element<Vertical>(l5); auto vertical1 = s.new_element<Vertical>(l5);
} }
std::shared_ptr<Curve> cc;
std::shared_ptr<CircularArc> c0;
Sketch s; Sketch s;
}; };
...@@ -670,6 +674,44 @@ TEST_F(HalfCircleInSquare, create) { ...@@ -670,6 +674,44 @@ TEST_F(HalfCircleInSquare, create) {
m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm"); m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
} }
TEST_F(HalfCircleInSquare, uniform_circle) {
std::string test_name{"half_circle_in_square_uniform_circle"};
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
s.save_as<SaveMethod::Rasterize>(SAVE_DIR, test_name);
m.MaximumElementSize = 0.1;
m.MinimumElementSize = 0.001;
m.MinimumElementQuality = 0.5;
auto bc = m.boundary_constraint(c0);
bc->uniform_discretization(true);
m.create();
m.save_as(SAVE_DIR, test_name + "_mesh");
edges_are_valid(m);
edges_are_optimal(m);
EXPECT_TRUE(bc->is_uniform(m));
m.refine();
edges_are_valid(m);
edges_are_optimal(m);
EXPECT_TRUE(bc->is_uniform(m));