Commit 58bf87a4 authored by JasonPries's avatar JasonPries
Browse files

Started refactorization of mesh refinement control algorithm

parent f8c716f8
......@@ -204,10 +204,17 @@ bool Mesh::refine() {
std::vector<size_t> index;
// TODO: Loop until quality is satisfied
// TODO: Iteratively decrease the minimum element size 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_ascending(quality, index);
sort_permutation_descending(radii, index);
size_t N = Triangles.size();
refine_once(index, radii, quality);
......@@ -216,12 +223,13 @@ bool Mesh::refine() {
while (M > N) {
N = M;
element_quality(radii, quality);
sort_permutation_ascending(quality, index);
//sort_permutation_ascending(quality, index);
sort_permutation_descending(radii, index);
refine_once(index, radii, quality);
M = Triangles.size();
}
return edges_are_valid();
return edges_are_valid(); // TODO: Instrument in tests
}
bool Mesh::refine_once() {
......@@ -230,11 +238,11 @@ bool Mesh::refine_once() {
std::vector<size_t> index;
element_quality(radii, quality);
sort_permutation_ascending(quality, index);
//sort_permutation_descending(radii, index);
//sort_permutation_ascending(quality, index);
sort_permutation_descending(radii, index);
refine_once(index, radii, quality);
return edges_are_valid();
return edges_are_valid(); // TODO: Instrument in tests
}
bool Mesh::in_triangle(Point const p, size_t ei) const {
......@@ -587,12 +595,12 @@ void Mesh::element_quality(std::vector<double> &radii, std::vector<double> &qual
radii.reserve(Triangles.size());
quality.reserve(Triangles.size());
for (size_t i = 0; i < Triangles.size(); ++i) {
for (size_t i = 0; i != Triangles.size(); ++i) {
double r = circumradius(Triangles[i]);
double l = shortest_edge_length(Triangles[i]);
radii.push_back(r);
quality.push_back(l / r);
quality.push_back(l / r / sqrt(3.0)); // sqrt(3.0) = (shortest edges length) / radius of equilateral triangle
}
}
......@@ -1281,8 +1289,8 @@ InsertPointResult Mesh::insert_point(Point const vc, size_t ei) {
Edge &e0 = Edges[--itr];
Edge &tri = Edges[ei];
Edge &nxt = Edges[tri.Next];//next(tri);
Edge &prv = Edges[tri.Prev];//prev(tri);
Edge &nxt = Edges[tri.Next];
Edge &prv = Edges[tri.Prev];
size_t vt = node(tri);
size_t vn = node(nxt);
......
#include "test_Mesh.hpp"
TEST(Mesh, create_triangle_domain) {
std::string test_name = "triangle_domain";
Sketch s;
auto v0 = s.new_element<Vertex>(1.0, 0.0);
......@@ -11,79 +13,85 @@ TEST(Mesh, create_triangle_domain) {
auto l1 = s.new_element<LineSegment>(v1, v2);
auto l2 = s.new_element<LineSegment>(v2, v0);
s.solve();
s.build();
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
m.create();
m.save_as(SAVE_DIR, "triangle_domain");
m.save_as(SAVE_DIR, test_name + "_mesh");
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2}, m);
{ // Test number of vertices, edges, triangles
EXPECT_TRUE(m.size_points() == 3);
EXPECT_TRUE(m.size_edges() == 3);
EXPECT_TRUE(m.size_triangles() == 1);
// Test number of vertices, edges, triangles
EXPECT_TRUE(m.size_points() == 3);
EXPECT_TRUE(m.size_edges() == 3);
EXPECT_TRUE(m.size_triangles() == 1);
EXPECT_TRUE(m.num_points() == 3);
EXPECT_TRUE(m.num_edges() == 3);
EXPECT_TRUE(m.num_triangles() == 1);
EXPECT_TRUE(m.num_points() == 3);
EXPECT_TRUE(m.num_edges() == 3);
EXPECT_TRUE(m.num_triangles() == 1);
EXPECT_TRUE(v0->x() == m.point(vmap[0]).X);
EXPECT_TRUE(v0->y() == m.point(vmap[0]).Y);
EXPECT_TRUE(v0->x() == m.point(vmap[0]).X);
EXPECT_TRUE(v0->y() == m.point(vmap[0]).Y);
EXPECT_TRUE(v1->x() == m.point(vmap[1]).X);
EXPECT_TRUE(v1->y() == m.point(vmap[1]).Y);
EXPECT_TRUE(v1->x() == m.point(vmap[1]).X);
EXPECT_TRUE(v1->y() == m.point(vmap[1]).Y);
EXPECT_TRUE(v2->x() == m.point(vmap[2]).X);
EXPECT_TRUE(v2->y() == m.point(vmap[2]).Y);
}
EXPECT_TRUE(v2->x() == m.point(vmap[2]).X);
EXPECT_TRUE(v2->y() == m.point(vmap[2]).Y);
{ // Test validity, optimality
edges_are_valid(m);
edges_are_optimal(m);
}
// Test validity, optimality
edges_are_valid(m);
edges_are_optimal(m);
{ // Test edge and node connections
const Edge e0 = m.edge(vmap[0]);
const Edge e1 = m.edge(vmap[1]);
const Edge e2 = m.edge(vmap[2]);
EXPECT_TRUE(v0->x() == m.point(e0.node()).X);
EXPECT_TRUE(v0->y() == m.point(e0.node()).Y);
EXPECT_TRUE(e1.self() == e0.next());
EXPECT_TRUE(e2.self() == e0.prev());
EXPECT_TRUE(v1->x() == m.point(e1.node()).X);
EXPECT_TRUE(v1->y() == m.point(e1.node()).Y);
EXPECT_TRUE(e2.self() == e1.next());
EXPECT_TRUE(e0.self() == e1.prev());
EXPECT_TRUE(v2->x() == m.point(e2.node()).X);
EXPECT_TRUE(v2->y() == m.point(e2.node()).Y);
EXPECT_TRUE(e0.self() == e2.next());
EXPECT_TRUE(e1.self() == e2.prev());
}
// Test edge and node connections
const Edge e0 = m.edge(vmap[0]);
const Edge e1 = m.edge(vmap[1]);
const Edge e2 = m.edge(vmap[2]);
{ // Test triangles
const Edge e = m.triangle(0);
Point cc = m.circumcenter(e.self());
EXPECT_NEAR(0.0, cc.X, TOL);
EXPECT_NEAR(sqrt(3.0) / 3.0, cc.Y, TOL);
EXPECT_TRUE(v0->x() == m.point(e0.node()).X);
EXPECT_TRUE(v0->y() == m.point(e0.node()).Y);
EXPECT_TRUE(e1.self() == e0.next());
EXPECT_TRUE(e2.self() == e0.prev());
double cr = m.circumradius(e.self());
EXPECT_NEAR(2.0 * sqrt(3.0) / 3.0, cr, TOL);
}
EXPECT_TRUE(v1->x() == m.point(e1.node()).X);
EXPECT_TRUE(v1->y() == m.point(e1.node()).Y);
EXPECT_TRUE(e2.self() == e1.next());
EXPECT_TRUE(e0.self() == e1.prev());
{ // Forced Refinement
forced_refinement(m, "triangle_domain_refine_loop", 7);
}
EXPECT_TRUE(v2->x() == m.point(e2.node()).X);
EXPECT_TRUE(v2->y() == m.point(e2.node()).Y);
EXPECT_TRUE(e0.self() == e2.next());
EXPECT_TRUE(e1.self() == e2.prev());
size_t i = m.num_points();
// Test triangles
const Edge e = m.triangle(0);
Point cc = m.circumcenter(e.self());
EXPECT_NEAR(0.0, cc.X, TOL);
EXPECT_NEAR(sqrt(3.0) / 3.0, cc.Y, TOL);
double cr = m.circumradius(e.self());
EXPECT_NEAR(2.0 * sqrt(3.0) / 3.0, cr, TOL);
// Forced Refinement
forced_refinement(m, test_name + "_mesh_refine_loop", 7);
// Test refinement algorithm
m = Mesh(s);
m.create();
m.MaximumElementSize = 0.1;
m.MinimumElementSize = 0.01;
m.MinimumElementQuality = 0.1;
m.refine();
m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
}
TEST(Mesh, create_square_domain) {
std::string test_name = "square_domain";
Sketch s;
auto v0 = s.new_element<Vertex>(0.0, 0.0);
......@@ -96,77 +104,80 @@ TEST(Mesh, create_square_domain) {
auto l2 = s.new_element<LineSegment>(v2, v3);
auto l3 = s.new_element<LineSegment>(v3, v0);
s.solve();
s.build();
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
m.create();
m.save_as(SAVE_DIR, "square_domain");
m.save_as(SAVE_DIR, test_name + "_mesh");
// Test number of verticies, edges, triangles
{
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3}, m);
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3}, m);
EXPECT_TRUE(m.size_points() == 4);
EXPECT_TRUE(m.size_edges() == 6);
EXPECT_TRUE(m.size_triangles() == 2);
EXPECT_TRUE(m.size_points() == 4);
EXPECT_TRUE(m.size_edges() == 6);
EXPECT_TRUE(m.size_triangles() == 2);
EXPECT_TRUE(m.num_points() == 4);
EXPECT_TRUE(m.num_edges() == 5);
EXPECT_TRUE(m.num_triangles() == 2);
EXPECT_TRUE(m.num_points() == 4);
EXPECT_TRUE(m.num_edges() == 5);
EXPECT_TRUE(m.num_triangles() == 2);
EXPECT_TRUE(v0->x() == m.point(vmap[0]).X);
EXPECT_TRUE(v0->y() == m.point(vmap[0]).Y);
EXPECT_TRUE(v0->x() == m.point(vmap[0]).X);
EXPECT_TRUE(v0->y() == m.point(vmap[0]).Y);
EXPECT_TRUE(v1->x() == m.point(vmap[1]).X);
EXPECT_TRUE(v1->y() == m.point(vmap[1]).Y);
EXPECT_TRUE(v1->x() == m.point(vmap[1]).X);
EXPECT_TRUE(v1->y() == m.point(vmap[1]).Y);
EXPECT_TRUE(v2->x() == m.point(vmap[2]).X);
EXPECT_TRUE(v2->y() == m.point(vmap[2]).Y);
EXPECT_TRUE(v2->x() == m.point(vmap[2]).X);
EXPECT_TRUE(v2->y() == m.point(vmap[2]).Y);
EXPECT_TRUE(v3->x() == m.point(vmap[3]).X);
EXPECT_TRUE(v3->y() == m.point(vmap[3]).Y);
}
EXPECT_TRUE(v3->x() == m.point(vmap[3]).X);
EXPECT_TRUE(v3->y() == m.point(vmap[3]).Y);
// Test validity, optimality
{
edges_are_valid(m);
edges_are_optimal(m);
}
edges_are_valid(m);
edges_are_optimal(m);
// Test edge and node connections
{
EXPECT_TRUE(m.point((size_t)0) == m.point(m.edge(0)));
EXPECT_TRUE(m.point((size_t)1) == m.point(m.edge(1)));
EXPECT_TRUE(m.point((size_t)2) == m.point(m.edge(2)));
EXPECT_TRUE(m.point((size_t)3) == m.point(m.edge(3)));
EXPECT_TRUE(m.point((size_t)1) == m.point(m.edge(4)));
EXPECT_TRUE(m.point((size_t)3) == m.point(m.edge(5)));
for (size_t i = 0; i < 5; i++) {
const Edge e = m.edge(i);
EXPECT_TRUE(e.self() == m.edge(e.next()).prev());
EXPECT_TRUE(e.self() == m.edge(e.prev()).next());
}
EXPECT_TRUE(m.point((size_t) 0) == m.point(m.edge(0)));
EXPECT_TRUE(m.point((size_t) 1) == m.point(m.edge(1)));
EXPECT_TRUE(m.point((size_t) 2) == m.point(m.edge(2)));
EXPECT_TRUE(m.point((size_t) 3) == m.point(m.edge(3)));
EXPECT_TRUE(m.point((size_t) 1) == m.point(m.edge(4)));
EXPECT_TRUE(m.point((size_t) 3) == m.point(m.edge(5)));
for (size_t i = 0; i < 5; i++) {
const Edge e = m.edge(i);
EXPECT_TRUE(e.self() == m.edge(e.next()).prev());
EXPECT_TRUE(e.self() == m.edge(e.prev()).next());
}
// Test triangles
{
for (size_t i = 0; i < m.size_triangles(); ++i) {
Point cc = m.circumcenter(m.triangle(0).self());
EXPECT_NEAR(0.5, cc.X, TOL);
EXPECT_NEAR(0.5, cc.Y, TOL);
}
for (size_t i = 0; i < m.size_triangles(); ++i) {
Point cc = m.circumcenter(m.triangle(0).self());
EXPECT_NEAR(0.5, cc.X, TOL);
EXPECT_NEAR(0.5, cc.Y, TOL);
}
// Forced Refinement
{
forced_refinement(m, "square_domain_refine_loop", 7);
}
forced_refinement(m, test_name + "_mesh_refine_loop", 7);
// Test refinement algorithm
m = Mesh(s);
m.create();
m.MaximumElementSize = 0.1;
m.MinimumElementSize = 0.01;
m.MinimumElementQuality = 0.1;
m.refine();
m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
}
TEST(Mesh, create_narrow_diamond_domain) {
std::string test_name = "narrow_diamond_domain";
Sketch s;
auto v0 = s.new_element<Vertex>(1.0, 0.0);
......@@ -179,69 +190,71 @@ TEST(Mesh, create_narrow_diamond_domain) {
auto l2 = s.new_element<LineSegment>(v2, v3);
auto l3 = s.new_element<LineSegment>(v3, v0);
s.solve();
s.build();
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
m.create();
m.save_as(SAVE_DIR, "narrow_diamond_domain");
m.save_as(SAVE_DIR, test_name + "_mesh");
// Test number of verticies, edges, triangles
{
EXPECT_TRUE(m.size_points() == 4);
EXPECT_TRUE(m.size_edges() == 6);
EXPECT_TRUE(m.size_triangles() == 2);
EXPECT_TRUE(m.num_points() == 4);
EXPECT_TRUE(m.num_edges() == 5);
EXPECT_TRUE(m.num_triangles() == 2);
}
EXPECT_TRUE(m.size_points() == 4);
EXPECT_TRUE(m.size_edges() == 6);
EXPECT_TRUE(m.size_triangles() == 2);
EXPECT_TRUE(m.num_points() == 4);
EXPECT_TRUE(m.num_edges() == 5);
EXPECT_TRUE(m.num_triangles() == 2);
// Test validity, optimality
{
edges_are_valid(m);
edges_are_optimal(m);
}
edges_are_valid(m);
edges_are_optimal(m);
// Test for proper edge swaps
{
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3}, m);
for (size_t i = 5; i < 6; i++) {
const Edge e = m.edge(5);
if (m.point(e.node()) == m.point(vmap[0])) {
EXPECT_TRUE(m.point(vmap[2]) == m.point(m.edge(e.next()).node()));
EXPECT_TRUE(m.point(vmap[2]) == m.point(m.edge(e.twin()).node()));
} else if (m.point(e.node()) == m.point(vmap[2])) {
EXPECT_TRUE(m.point(vmap[0]) == m.point(m.edge(e.next()).node()));
EXPECT_TRUE(m.point(vmap[0]) == m.point(m.edge(e.twin()).node()));
}
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3}, m);
for (size_t i = 5; i < 6; i++) {
const Edge e = m.edge(5);
if (m.point(e.node()) == m.point(vmap[0])) {
EXPECT_TRUE(m.point(vmap[2]) == m.point(m.edge(e.next()).node()));
EXPECT_TRUE(m.point(vmap[2]) == m.point(m.edge(e.twin()).node()));
} else if (m.point(e.node()) == m.point(vmap[2])) {
EXPECT_TRUE(m.point(vmap[0]) == m.point(m.edge(e.next()).node()));
EXPECT_TRUE(m.point(vmap[0]) == m.point(m.edge(e.twin()).node()));
}
}
// Test triangle circumcenters
{
Point cc0 = m.circumcenter(m.triangle(0).self());
Point cc1 = m.circumcenter(m.triangle(1).self());
Point cc0 = m.circumcenter(m.triangle(0).self());
Point cc1 = m.circumcenter(m.triangle(1).self());
EXPECT_NEAR(0.0, cc0.X, TOL);
EXPECT_NEAR(0.75, std::abs(cc0.Y), TOL);
EXPECT_NEAR(0.0, cc0.X, TOL);
EXPECT_NEAR(0.75, std::abs(cc0.Y), TOL);
EXPECT_NEAR(0.0, cc1.X, TOL);
EXPECT_NEAR(0.75, std::abs(cc1.Y), TOL);
EXPECT_NEAR(-cc0.Y, cc1.Y, TOL);
}
EXPECT_NEAR(0.0, cc1.X, TOL);
EXPECT_NEAR(0.75, std::abs(cc1.Y), TOL);
EXPECT_NEAR(-cc0.Y, cc1.Y, TOL);
// Forced Refinement
{
forced_refinement(m, "narrow_diamond_domain_refine_loop", 7);
}
forced_refinement(m, test_name + "_mesh_refine_loop", 7);
// Test refinement algorithm
m = Mesh(s);
m.create();
m.MaximumElementSize = 0.1;
m.MinimumElementSize = 0.01;
m.MinimumElementQuality = 0.1;
m.refine();
m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
}
TEST(Mesh, create_narrow_rectangle_domain) {
std::string test_name = "narrow_rectangle_domain";
Sketch s;
auto v0 = s.new_element<Vertex>(5.0, 0.0);
......@@ -256,57 +269,63 @@ TEST(Mesh, create_narrow_rectangle_domain) {
auto l3 = s.new_element<LineSegment>(v3, v4);
auto l4 = s.new_element<LineSegment>(v4, v0);
s.solve();
s.build();
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
m.create();
m.save_as(SAVE_DIR, "narrow_rectangle_domain");
m.save_as(SAVE_DIR, test_name + "_mesh");
// Test number of verticies, edges, triangles
{
EXPECT_TRUE(m.size_points() == 6);
EXPECT_TRUE(m.size_edges() == 12);
EXPECT_TRUE(m.size_triangles() == 4);
EXPECT_TRUE(m.num_points() == 6);
EXPECT_TRUE(m.num_edges() == 9);
EXPECT_TRUE(m.num_triangles() == 4);
}
EXPECT_TRUE(m.size_points() == 6);
EXPECT_TRUE(m.size_edges() == 12);
EXPECT_TRUE(m.size_triangles() == 4);
EXPECT_TRUE(m.num_points() == 6);
EXPECT_TRUE(m.num_edges() == 9);
EXPECT_TRUE(m.num_triangles() == 4);
// Test validity, optimality
{
edges_are_valid(m);
edges_are_optimal(m);
}
edges_are_valid(m);
edges_are_optimal(m);
// Test edge splits
{
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3, v4}, m);
std::vector<size_t> vmap = map_verticies_to_points({v0, v1, v2, v3, v4}, m);
Point const &v5 = m.point(5);
EXPECT_NEAR(0.0, v5.X, DBL_EPSILON);
EXPECT_NEAR(0.0, v5.Y, DBL_EPSILON);
Point const &v5 = m.point(5);
EXPECT_NEAR(0.0, v5.X, DBL_EPSILON);
EXPECT_NEAR(0.0, v5.Y, DBL_EPSILON);
for (size_t i = 0; i < 12; ++i) {
Edge const e = m.edge(i);
for (size_t i = 0; i < 12; ++i) {
Edge const e = m.edge(i);
if (m.point(e.node()) == m.point(vmap[2])) {
if (e.twin() != e.self()) {
EXPECT_TRUE(m.point(m.edge(e.twin()).node()) == v5);
}
if (m.point(e.node()) == m.point(vmap[2])) {
if (e.twin() != e.self()) {
EXPECT_TRUE(m.point(m.edge(e.twin()).node()) == v5);
}
}
}
// Forced Refinement
{
forced_refinement(m, "narrow_rectangle_domain_refine_loop", 7);
}
forced_refinement(m, test_name + "_mesh_refine_loop", 7);
// Test refinement algorithm
m = Mesh(s);
m.create();
m.MaximumElementSize = 0.1;
m.MinimumElementSize = 0.01;
m.MinimumElementQuality = 0.1;
m.refine();
m.save_as(SAVE_DIR, test_name + "_mesh_refine_algorithm");
}
TEST(Mesh, create_half_circle_domain) {
std::string test_name = "half_circle_domain";
Sketch s;
auto vc = s.new_element<Vertex>(0.0, 0.0);
......@@ -316,38 +335,45 @@ TEST(Mesh, create_half_circle_domain) {
auto arc = s.new_element<CircularArc>(v0, v1, vc, 1.0);
auto line = s.new_element<LineSegment>(v1, v0);
s.solve();
s.build();
double res_norm = s.solve();
EXPECT_LE(res_norm, FLT_EPSILON);
bool result = s.build();
ASSERT_TRUE(result);
Mesh m{s};
m.create();
m.save_as(SAVE_DIR, "half_circle_domain");
m.save_as(SAVE_DIR, test_name + "_mesh");
// Test number of verticies, edges, triangles
{
EXPECT_TRUE(m.size_points() == 6);
EXPECT_TRUE(m.size_edges() == 12);
EXPECT_TRUE(m.size_triangles() == 4);
EXPECT_TRUE(m.num_points() == 6);
EXPECT_TRUE(m.num_edges() == 9);
EXPECT_TRUE(m.num_triangles() == 4);
}
EXPECT_TRUE(m.size_points() == 6);
EXPECT_TRUE(m.size_edges() == 12);
EXPECT_TRUE(m.size_triangles() == 4);
EXPECT_TRUE(m.num_points() == 6);
EXPECT_TRUE(m.num_edges() == 9);
EXPECT_TRUE(m.num_triangles() == 4);
// Test validity, optimality
{
edges_are_valid(m);
edges_are_optimal(m);
}
edges_are_valid(m);
edges_are_optimal(m);