Commit 27c32224 authored by JasonPries's avatar JasonPries

Implements Mesh smoothing and MappedBoundaryPair

Mesh::smooth() is based on satisfaction of a Laplacian condition

MappedBoundaryPair guarantees correct behavior of PeriodicBoundaryConditions when the UniformDiscretization option is set to false
parent 8b95e011
......@@ -7,9 +7,10 @@ set(SOURCE_FILES
./src/DartConstraint.h ./src/DartConstraint.cpp
./src/Edge.h ./src/Edge.cpp
./src/InsertionQueuer.h ./src/InsertionQueuer.cpp
./src/MappedBoundaryPair.h ./src/MappedBoundaryPair.cpp
./src/Mesh.h ./src/Mesh.cpp
./src/Point.h ./src/Point.cpp
./src/DartTriangle.h ./src/DartTriangle.cpp
./src/DartTriangle.h ./src/DartTriangle.cpp
)
add_library(mesh SHARED ${SOURCE_FILES})
......
......@@ -6,6 +6,7 @@
#include "../src/DartTriangle.h"
#include "../src/Edge.h"
#include "../src/InsertionQueuer.h"
#include "../src/MappedBoundaryPair.h"
#include "../src/Mesh.h"
#include "../src/Point.h"
......
......@@ -6,7 +6,7 @@ bool BoundaryConstraint::is_uniform(Mesh const &m) const {
DartConstraint const &dc = m.dart_constraint(DartConstraints[0]);
double_t delta_0{dc.delta()};
for(size_t i : DartConstraints) {
for (size_t i : DartConstraints) {
double delta_i{m.dart_constraint(i).delta()};
if (std::abs(delta_0 - delta_i) > FLT_EPSILON) {
return false;
......@@ -28,14 +28,15 @@ double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
}
double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
std::cout << "//TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued" << std::endl;
double delta_max{delta_min};
for(size_t i : DartConstraints) {
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_dart_to_queue(dc.dart());
delta_max = std::max(delta_max, delta);
}
}
......@@ -47,10 +48,17 @@ 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));
m.add_dart_to_queue(dart);
//m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
}
}
void BoundaryConstraint::insert_parameters(Mesh &mesh, std::vector<double_t> s) {
while (s.size() > 0) {
queue_parameters(mesh, s);
mesh.insert_from_queue();
}
}
void BoundaryConstraint::make_uniform(Mesh &m) { make_uniform(m, smallest_parametric_edge_length(m)); };
......@@ -63,6 +71,43 @@ void BoundaryConstraint::make_uniform(Mesh &m, double delta_min) {
}
}
void BoundaryConstraint::queue_parameters(Mesh &mesh, std::vector<double_t> &s) {
/*
Enqueues all darts having midpoints coincident with the parameter vector s
The parameter vector s is consumed, with s[i] being erased if an associated dart is enqueued
TODO: This method assumes all values in s are powers of 2 and should check that assumption
*/
sort_dart_constraints(mesh);
std::sort(s.begin(), s.end());
size_t j{0};
for (size_t i : DartConstraints) {
auto const &dc = mesh.dart_constraint(i);
while (j != s.size() && s[j] < dc.S1 + FLT_EPSILON) {
if (std::abs(dc.S0 + dc.S1 - 2.0 * s[j]) < FLT_EPSILON) { // if s[j] is the dart midpoint, split the dart in half
mesh.add_dart_to_queue(dc.dart());
s.erase(s.begin() + j);
break;
}
if (std::abs(dc.S0 - s[j]) < FLT_EPSILON || std::abs(dc.S1 - s[j]) < FLT_EPSILON) { // if s[j] is either dart endpoint, remove s from the queue
s.erase(s.begin() + j);
break;
}
++j;
}
if (s.size() == 0) {
break;
}
}
}
void BoundaryConstraint::refine(Mesh &m, double_t tol) {
double_t err_max{DBL_MAX};
......@@ -87,9 +132,8 @@ void BoundaryConstraint::refine(Mesh &m, double_t tol) {
}
void BoundaryConstraint::sort_dart_constraints(Mesh const &m) {
std::sort(DartConstraints.begin(), DartConstraints.end(),
[&](size_t i, size_t j) { return m.dart_constraint(i).S0 < m.dart_constraint(j).S0; }
);
auto comp = [&](size_t i, size_t j) { return m.dart_constraint(i).S0 < m.dart_constraint(j).S0; };
std::sort(DartConstraints.begin(), DartConstraints.end(), comp);
}
std::vector<size_t> BoundaryConstraint::nodes(Mesh const &m) {
......@@ -102,4 +146,14 @@ std::vector<size_t> BoundaryConstraint::nodes(Mesh const &m) {
n.push_back(m.dart_constraint(DartConstraints.back()).tip(m));
return n;
}
std::vector<double_t> BoundaryConstraint::parameters(Mesh const &m) {
std::vector<double_t> s;
for (size_t d : DartConstraints) {
s.emplace_back(m.dart_constraint(d).S0);
}
s.emplace_back(1.0);
return s;
}
\ No newline at end of file
......@@ -28,27 +28,43 @@ public:
void add_dart_constraint(size_t d) { DartConstraints.push_back(d); }; // TODO: Enforce uniqueness (std::set?)
void uniform_discretization(bool flag) { UniformDiscretization = flag; };
void add_to_queue(Mesh &m, size_t dart) const;
void refine(Mesh &m, double_t tol);
void insert_parameters(Mesh &mesh, std::vector<double_t> s);
void make_uniform(Mesh &m);
void make_uniform(Mesh &m, double delta_min);
void queue_parameters(Mesh &mesh, std::vector<double_t> &s);
void refine(Mesh &m, double_t tol);
void sort_dart_constraints(Mesh const &m);
void uniform_discretization(bool flag) { UniformDiscretization = flag; };
std::shared_ptr<Curve const> curve() const { return ConstraintCurve; };
std::vector<double_t> parameters(Mesh const &m);
std::vector<size_t> nodes(Mesh const &m);
protected:
std::vector<size_t> DartConstraints;
std::shared_ptr<Curve const> ConstraintCurve;
bool UniformDiscretization{false};
bool UniformDiscretization{false}; // TODO: Revist this strategy (think about making a UniformBoundaryConstraint object)
// TODO: Mapped Boundaries
// TODO: Accept PeriodicBoundaryPair (from Sketch)
// TODO: construct a DiscreteBoundaryMap
// TODO: When a Dart is queued, pass the parametric point and this pointer to DiscreteBoundaryMap
// TODO: DiscreteBoundaryMap then calls BoundaryConstraint->insert(double_t) on all BoundaryConstraint objects except the one passed to the DiscreteBoundaryMap
// TODO: Must implement BoundaryConstraint->insert(double_t) which searches through the DartConstraints to find the DartConstraint representing the correct range (param min < param insert < param max)
// TODO: Could just implement this first and take care of boundary mapping after mesh refinement is complete
// TODO: OR: Enforce boundary mapping after each refinement iteration <- this is better because there is no need to prevent infinite recursion <<< DO THIS ONE
};
......
#include "MappedBoundaryPair.h"
#include "Mesh.h"
MappedBoundaryPair::MappedBoundaryPair(Mesh &mesh, ContinuousBoundaryPair const &cbp) {
Boundary0 = mesh.boundary_constraint(cbp.curve0());
Boundary1 = mesh.boundary_constraint(cbp.curve1());
Orientation = cbp.orientation();
}
bool MappedBoundaryPair::enforce_mapping(Mesh &mesh) const {
/*
*
*/
// Get parameters
std::vector<double_t> p0 = Boundary0->parameters(mesh);
std::vector<double_t> p1 = Boundary1->parameters(mesh);
if (!Orientation) { // reverse direction of p0
for (double_t &val : p0) {
val = 1.0 - val;
}
}
// Determine unique parameter mappings
std::sort(p0.begin(), p0.end());
std::sort(p1.begin(), p1.end());
std::vector<double_t> p0_diff;
std::set_difference(p0.begin(), p0.end(), p1.begin(), p1.end(), back_inserter(p0_diff));
std::vector<double_t> p1_diff;
std::set_difference(p1.begin(), p1.end(), p0.begin(), p0.end(), back_inserter(p1_diff));
if (p0_diff.size() == 0 && p1_diff.size() == 0) { // If boundary mapping is already satisfied, return
return true;
}
if (!Orientation) { // reverse direction of p1
for (double_t &val : p1_diff) {
val = 1.0 - val;
}
}
Boundary0->insert_parameters(mesh, p1_diff);
Boundary1->insert_parameters(mesh, p0_diff);
return false;
}
\ No newline at end of file
#ifndef OERSTED_MAPPEDBOUNDARYPAIR_H
#define OERSTED_MAPPEDBOUNDARYPAIR_H
#include "Sketch.hpp"
class Mesh;
#include "BoundaryConstraint.h"
class MappedBoundaryPair {
public:
MappedBoundaryPair(Mesh &mesh, ContinuousBoundaryPair const &cbp);
bool enforce_mapping(Mesh &mesh) const;
protected:
std::shared_ptr<BoundaryConstraint> Boundary0;
std::shared_ptr<BoundaryConstraint> Boundary1;
bool Orientation;
};
#endif //OERSTED_MAPPEDBOUNDARYPAIR_H
......@@ -85,6 +85,20 @@ bool Mesh::are_intersecting(size_t ei, size_t ej) const {
}
}
bool Mesh::edges_are_optimal() const {
/*
* Returns true if all unconstrained edges satisfy the empty circumcircle property
*/
for (size_t e = 0; e != Edges.size(); ++e) {
if (!is_locally_optimal(e)) {
return false;
}
}
return true;
}
bool Mesh::edges_are_valid() const {
bool result = true;
......@@ -200,40 +214,51 @@ bool Mesh::find_attached(Point const p, size_t &e_out) {
}
bool Mesh::refine() {
/*
* Iteratively refines the mesh based on user defined size and quality criteria
*
* TODO: Loop until quality is satisfied
* TODO: Iteratively decrease the min and max element size until quality is satisfied
* TODO: First: refine until maximum element size criteria is satisfied
* TODO: plan() out iterative maximum element size refinement
* TODO: Then: refine until element quality criteria is satisfied
* TODO: ?somehow iterate?
*
* TODO: Implement sort permutation option
*/
std::vector<double> radii;
std::vector<double> quality;
std::vector<size_t> index;
// TODO: Loop until quality is satisfied
// TODO: Iteratively decrease the min and max element size until quality is satisfied
// TODO: plan() bounds on maximum element quality
// TODO: First: refine until maximum element size criteria is satisfied
// TODO: plan() out iterative maximum element size refinement
// TODO: Then: refine until element quality cirteria is satisfied
// TODO: ?somehow iterate?
// TODO: SMOOTHING!
element_quality(radii, quality);
//sort_permutation_ascending(quality, index);
sort_permutation_descending(radii, index);
size_t N = Triangles.size();
refine_once(index, radii, quality);
size_t M = Triangles.size();
size_t prev_points_size = 0;
size_t smooth_threshold = Points.size() + Edges.size() / 3; // Heuristic, 3 Darts per triangle
while (Points.size() > prev_points_size) {
prev_points_size = Points.size();
while (M > N) {
N = M;
element_quality(radii, quality);
//sort_permutation_ascending(quality, index);
// TODO: Implement sort permutation option
//sort_permutation_ascending(quality, index);s
sort_permutation_descending(radii, index);
refine_once(index, radii, quality);
M = Triangles.size();
if (Points.size() > smooth_threshold || Points.size() <= prev_points_size) { // if the loop will exit, performing a smoothing iteration
enforce_boundary_mappings();
smooth();
smooth_threshold = Points.size() + Edges.size() / 3;
}
get_triangles();
}
return edges_are_valid(); // TODO: Instrument in tests
}
bool Mesh::refine_once() {
/*
* Performs one iteration of refinement without smoothing or enforcing boundary mapping constraints
*/
std::vector<double> radii;
std::vector<double> quality;
std::vector<size_t> index;
......@@ -242,6 +267,7 @@ bool Mesh::refine_once() {
//sort_permutation_ascending(quality, index);
sort_permutation_descending(radii, index);
refine_once(index, radii, quality);
get_triangles();
return edges_are_valid(); // TODO: Instrument in tests
}
......@@ -537,11 +563,29 @@ size_t Mesh::num_edges() const {
return count;
}
void Mesh::add_dart_to_queue(size_t d) {
/*
* Adds a Dart and all other darts encroached by its midpoint to the refinement queue
*/
std::vector<size_t> encroached_darts = get_encroached_edges(midpoint(d), d);
for (size_t ed : encroached_darts) {
add_to_queue(std::make_unique<MidpointQueuer>(ed));
}
}
void Mesh::create() {
/*
* Creates the initial mesh by inserting nodes only on the mesh constraint boundaries
*/
create_boundary_polygon();
triangulate_boundary_polygon();
insert_internal_boundaries();
enforce_boundary_mappings();
get_triangles();
}
......@@ -609,7 +653,29 @@ void Mesh::element_quality(std::vector<double> &radii, std::vector<double> &qual
}
}
void Mesh::enforce_boundary_mappings() {
/*
* Inserts nodes on boundaries until all boundary mapping constraints are satisfied
*
* TODO: Will this strategy converge for cyclical mappings?
*/
bool all_enforced{false};
while (!all_enforced) {
all_enforced = true;
for (auto &mb : MappedBoundaries) {
all_enforced &= mb->enforce_mapping(*this);
}
}
}
void Mesh::get_triangles() {
/*
* Parses the mesh for darts representing unique triangles
*
* // TODO: Instead of starting from scratch for each iteration, sweep and mark existing triangles, then parse the remainder
*/
Triangles.resize(0);
Triangles.reserve(2 * num_points());
......@@ -674,6 +740,8 @@ void Mesh::insert_internal_boundaries() {
If encroached, split edge
*/
sort_constraints_by_length();
// Find interior curves
std::vector<size_t> interior_index;
for (size_t i = 0; i != BoundaryConstraints.size(); ++i) {
......@@ -775,9 +843,20 @@ void Mesh::insert_internal_boundaries() {
}
}
sort_constraints();
// TODO: Enforce BoundaryConstraint conditions (e.g. uniform discretization, boundary maps)
}
void Mesh::make_edges_optimal() {
/*
* Swaps edges in the triangulation until the circumcircle property is satisfied for all unconstrained edges
*/
for (size_t i = 0; i != Edges.size(); ++i) {
recursive_swap(i);
}
}
void Mesh::mark_triangles() {
for (size_t i = 0; i != Edges.size(); ++i) {
Edges[i].Mark = true;
......@@ -799,7 +878,6 @@ void Mesh::refine_once(std::vector<size_t> index, std::vector<double> radii, std
insert_from_queue();
}
}
get_triangles();
}
void Mesh::insert_from_queue() {
......@@ -831,14 +909,108 @@ void Mesh::save_as(std::string path, std::string file_name) const {
fs.close();
}
void Mesh::smooth() {
/*
* Moves unconstrained nodes in the triangulation so that the mesh satisfies a laplacian optimality condition
*
* TODO: Apply boundary conditions such that the A matrix is symmetric => can use Cholesky factorization
* TODO: Could also use iterative method since current node locations give good initial guess
* TODO: Could also think of this a discrete system, where the present implementation is an implicit time step
* TODO: The question then is, could we use an explicit time step and still have the smoothing converge?
* TODO: The question is equivalent to, what are the eigenvalues of the matrix?
*/
size_t num_points = Points.size();
std::vector<Eigen::Triplet<double_t>> triplets;
Eigen::VectorXd b = Eigen::VectorXd::Zero(2 * num_points);
Eigen::VectorXd q;
for (auto &e_i : Edges) {
bool node_is_constrained = e_i.is_constrained();
auto e_j = e_i;
while (!node_is_constrained) {
e_j = next(twin(e_j));
if (e_j == e_i) {
break;
}
node_is_constrained = node_is_constrained || e_j.is_constrained();
}
if (node_is_constrained) {
size_t n = node(e_i);
Point p = Points[n];
n *= 2;
triplets.emplace_back(n + 0, n + 0, 1.0);
triplets.emplace_back(n + 1, n + 1, 1.0);
b(n + 0) += p.X;
b(n + 1) += p.Y;
} else {
size_t n0 = 2 * node(e_i);
size_t n1 = 2 * node(next(e_i));
triplets.emplace_back(n0 + 0, n0 + 0, 1.0);
triplets.emplace_back(n0 + 0, n1 + 0, -1.0);
triplets.emplace_back(n0 + 1, n0 + 1, 1.0);
triplets.emplace_back(n0 + 1, n1 + 1, -1.0);
}
}
Eigen::SparseMatrix<double> A(2 * num_points, 2 * num_points);
A.setFromTriplets(triplets.begin(), triplets.end());
Eigen::SparseLU<Eigen::SparseMatrix<double>> LU;
LU.compute(A);
if(LU.info()!=Eigen::Success) {
std::cerr << "Factorization of mesh smoothing matrix was not successful";
return;
}
q = LU.solve(b);
if(LU.info()!=Eigen::Success) {
std::cerr << "Mesh smoothing equation could not be solved";
return;
}
for (size_t i = 0; i != num_points; ++i) {
size_t n = 2 * i;
Points[i].X = q(n);
Points[i].Y = q(n + 1);
}
make_edges_optimal();
add_encroached_edges_to_queue_and_insert();
if(!edges_are_optimal()) {
std::cerr << "Edges are not optimal" << std::endl;
throw;
}
};
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());
}
);
/*
* Sorts BoundaryConstraints by pointer address for fast lookup
*/
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) {return (size_t) (x->curve().get()) < (size_t) (y->curve().get());};
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(), comp);
}
void Mesh::sort_constraints_by_length() {
/*
* Sorts BoundaryConstraints by length for internal boundary insertion
*/
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); };
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(), comp);
};
void Mesh::sort_permutation_ascending(std::vector<double> &values, std::vector<size_t> &index) const {
index.resize(values.size());
std::iota(index.begin(), index.end(), 0);
......@@ -855,9 +1027,9 @@ void Mesh::split_edge(size_t ei) {
/*
Splits edge into two edges at the midpoint without creating any new triangles.
Used for initial polygon refinement.
*/
// TODO: Rename to split_boundary_edge
// TODO: Rename to split_boundary_edge
*/
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;
......@@ -913,6 +1085,12 @@ void Mesh::split_edge(size_t ei) {
}
void Mesh::split_encroached_edges() {
/*
* Inserts the midpoints of all edges until they are no longer encroached
*
* // TODO: Review where this should be used versus add_encroached_edges_to_queue_and_insert
*/
bool any_split = true;
while (any_split) {
any_split = false;
......@@ -927,7 +1105,38 @@ void Mesh::split_encroached_edges() {
}
}
void Mesh::add_encroached_edges_to_queue_and_insert() {
/*
* Inserts the midpoints of all edges until they are no longer encroached
*/
bool any_encroached = true;
while (any_encroached) {
any_encroached = false;
size_t num_edges_at_start = Edges.size();
for (size_t e = 0; e != num_edges_at_start; ++e) {
if (is_constrained(e)) {
if (is_encroached(base(prev(e)), e)) { // base(prev(e)) == node in triangle opposite of edge e
any_encroached = true;
std::vector<size_t> encroached_edges = get_encroached_edges(midpoint(e), e);
for (size_t ee : encroached_edges) {
Edges[ee].add_to_queue(*this);
}
insert_from_queue();