Skip to content
GitLab
Commit a49cc27a authored by Pries, Jason's avatar Pries, Jason
Browse files

Fixed some issues with meshing by adjust delaunay optimality test and reverted... 

Fixed some issues with meshing by adjust delaunay optimality test and reverted back to old Mesh::insert_internal_boundaries(); implementation
parent 87085af2
Hide whitespace changes
Inline Side-by-side
src/Mesh/src/Mesh.cpp
View file @ a49cc27a
......@@ -448,12 +448,14 @@ bool Mesh::is_valid(size_t ei) const {
bool Mesh::recursive_swap(size_t ei) {
// TODO, May need to have two different recursive swap methods, one for midpoint insertion and one for circumcenter insertion
if (!is_locally_optimal(ei) && swap(ei)) {
if (!is_locally_optimal(ei)){//} && swap(ei)) {
size_t enext = next(ei);
size_t eprev = prev(ei);
size_t tnext = next(twin(ei));
size_t tprev = prev(twin(ei));
swap(ei);
recursive_swap(enext);
recursive_swap(eprev);
recursive_swap(tnext);
......@@ -595,7 +597,7 @@ void Mesh::create() {
insert_internal_boundaries();
//get_triangles();
get_triangles();
}
void Mesh::create_boundary_polygon() {
......@@ -765,7 +767,8 @@ void Mesh::insert_internal_boundaries() {
// TODO: Then repeat for internal boundaries with one existing end point
// TODO: Or just add a generic bounding box + insert
sort_constraints_by_length(); // TODO: Method is sensitive to sorting of constraints
sort_constraints_by_length();
// TODO: Method is sometimes sensitive to sorting of constraints
// TODO: Method is also sensitive to lines connecting boundary verticies intersecting with internal features
// TODO: I.E., straight line approximations of the geometry hull must not intersect internal boundaries
......@@ -786,7 +789,6 @@ void Mesh::insert_internal_boundaries() {
}
}
/*
// Attach edges to constraints by inserting interior curve midpoints until constraints are naturally satisfied
std::list<size_t> queue;
for (size_t i : interior_index) {
......@@ -812,153 +814,53 @@ void Mesh::insert_internal_boundaries() {
queue.push_back(DartConstraints.size());
new_dart_constraint(0.0, 1.0, BoundaryConstraints[i]);
while (queue.size() != 0) {
DartConstraint &dc = DartConstraints[queue.back()];
Point p0 = dc.curve()->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
Point p1 = dc.curve()->point(dc.S1);
size_t ei = Edges.size() - 1;
LocateTriangleResult result = locate_triangle(p0, ei);
if (result != LocateTriangleResult::Point) {
std::cerr << "Could not locate edge referencing point, aborting mesh" << std::endl;
return;
}
if (find_attached(p1, ei)) {
// if ei is an edge attaching p0 and p1, ei coincides with dc and is a constrained edge for BoundaryConstraint[i]
Edge &e = Edges[ei];
e.Constraint = dc.Self;
e.Orientation = true;
dc.forward_dart(e.self());
Edge &et = Edges[e.Twin];
et.Constraint = dc.Self;
et.Orientation = false;
dc.reverse_dart(et.self());
queue.pop_back();
} else {
// if p0 and p1 are not connected by ei, insert the midpoint of dc
double_t s0 = dc.S0;
double_t s1 = dc.S1;
double_t sn = (s0 + s1) / 2.0;
Point const p = dc.curve()->point(sn);
//dc.S1 = sn;
dc.S0 = sn;
queue.push_back(DartConstraints.size());
//new_dart_constraint(sn, s1, dc.boundary_constraint());
new_dart_constraint(s0, sn, dc.boundary_constraint());
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
}
}
*/
// Insert all interior start and end points
double_t len_min = std::min(BoundaryConstraints[interior_index[0]]->curve()->length(),BoundaryConstraints[interior_index.back()]->curve()->length());
for (size_t i : interior_index) {
size_t N = (size_t)std::pow(2,std::ceil(std::log2(BoundaryConstraints[i]->curve()->length() / len_min / 2.0 + 1) + 1));
std::cout << N << std::endl;
for (size_t j = 0; j < N; ++j) {
double_t s = j * 1.0 / (N - 1);
Point p = BoundaryConstraints[i]->curve()->point(s);
LocateTriangleResult result = locate_triangle(p);
if (result == LocateTriangleResult::Interior) {
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
}
make_edges_optimal();
for (size_t j = 0; j < N - 1; ++j) {
double_t s0 = j * 1.0 / (N - 1);
Point p0 = BoundaryConstraints[i]->curve()->point(s0);
double_t s1 = (j+1) * 1.0 / (N - 1);
Point p1 = BoundaryConstraints[i]->curve()->point(s1);
size_t e{0};
locate_triangle(p1,e);
if (!find_attached(p0,e)) {
std::cout << "SOMETHING IS STILL WRONG" << std::endl;
}
}
}
return;
/*
while (queue.size() > 0) {
DartConstraint &dc = queue.back();
while (queue.size() != 0) {
DartConstraint &dc = DartConstraints[queue.back()];
Point p0 = dc.curve()->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
Point p1 = dc.curve()->point(dc.S1);
size_t ei = Edges.size() - 1;
LocateTriangleResult result = locate_triangle(p0, ei);
if (result == LocateTriangleResult::Exterior) {
std::cout << "Point is exterior to mesh, aborting" << std::endl;
std::cout << p0.X << "," << p0.Y << std::endl;
return;
} else if (result == LocateTriangleResult::Interior) {
std::cout << "Point is interior to mesh element, aborting" << std::endl;
std::cout << p0.X << "," << p0.Y << std::endl;
return;
if (result != LocateTriangleResult::Point) {
std::runtime_error(std::string("Could not locate edge referencing point which should exist in mesh"));
}
if ((result == LocateTriangleResult::Point) && find_attached(p1, ei)) {
if (find_attached(p1, ei)) {
// if ei is an edge attaching p0 and p1, ei coincides with dc and is a constrained edge for BoundaryConstraint[i]
new_dart_constraint(dc.S0, dc.S1, dc.boundary_constraint());
DartConstraint &dcnew = DartConstraints.back();
Edge &e = Edges[ei];
e.Constraint = dcnew.Self;
e.Constraint = dc.Self;
e.Orientation = true;
dcnew.forward_dart(e.self());
dc.forward_dart(e.self());
Edge &et = Edges[e.Twin];
et.Constraint = dcnew.Self;
et.Constraint = dc.Self;
et.Orientation = false;
dcnew.reverse_dart(et.self());
dc.reverse_dart(et.self());
queue.pop_back();
split_encroached_edges();
} else {
// if p0 and p1 are not connected by ei, insert the midpoint of dc
double_t s0 = dc.S0;
double_t s1 = dc.S1;
double_t sn = (s0 + s1) / 2.0;
Point const p = dc.curve()->point(sn);
AddToQueueResult atqr = add_point_to_queue(p);
if (atqr != AddToQueueResult::Duplicate) {
while (atqr == AddToQueueResult::Midpoint) {
insert_from_queue();
atqr = add_point_to_queue(p);
}
insert_from_queue();
} else {
std::cout << "Duplicate detected" << std::endl;
}
dc.S0 = sn;
queue.emplace_back(s0, sn, dc.boundary_constraint(), dc.self());
queue.push_back(DartConstraints.size());
new_dart_constraint(s0, sn, dc.boundary_constraint());
while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
insert_from_queue();
}
insert_from_queue();
}
}
*/
return;
split_encroached_edges();
for (auto bc : BoundaryConstraints) { // Length based refinement
......@@ -1148,6 +1050,7 @@ void Mesh::sort_constraints() {
* Sorts BoundaryConstraints by pointer address for fast lookup
*/
std::cerr << "DO NOT SORT BY POINTER ADDRESS. THIS MAKES THINGS NONDETERMINISTIC" << std::endl;
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);
......@@ -1158,8 +1061,10 @@ 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()); }; // smallest to largest
//auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); }; // largest to smallest
// TODO: sorting order should not matter for mesh construction, but sometimes it might
//auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); }; // smallest to largest
auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); }; // largest to smallest
std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(), comp);
};
......
test/LibraryIntegration/test_Circle.cpp
View file @ a49cc27a
......@@ -269,6 +269,6 @@ TEST_F(Circle, magnetostatic_sinusoidal_forcing) {
test_msf_solution(magnetostatic_sinusoidal_forcing<1,1>(), 0.025, 0.05, FLT_MAX, FLT_MAX);
test_msf_solution(magnetostatic_sinusoidal_forcing<1,2>(), 0.025, 0.1, FLT_MAX, FLT_MAX);
test_msf_solution(magnetostatic_sinusoidal_forcing<2,2>(), 0.005, 0.01, 0.75, 0.001);
test_msf_solution(magnetostatic_sinusoidal_forcing<2,2>(), 0.01, 0.013, 0.75, 0.001);
test_msf_solution(magnetostatic_sinusoidal_forcing<2,4>(), 0.005, 0.03, 1.0, 0.001);
}
\ No newline at end of file
test/Mesh/Mesh_Test.cpp
View file @ a49cc27a
......@@ -1039,17 +1039,17 @@ TEST(Mesh, locate_triangle__triangular_domain) {
size_t loc = i;
vp = ve0;
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Interior);
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[0])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[1])));
vp = ve1;
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Interior);
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[1])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[2])));
vp = ve2;
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Interior);
EXPECT_TRUE(mesh.locate_triangle(vp, loc) == LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[2])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[0])));
}
......@@ -1138,25 +1138,25 @@ TEST(Mesh, locate_triange__square_domain) {
vp = ve0;
result = mesh.locate_triangle(vp, loc);
EXPECT_EQ(result, LocateTriangleResult::Interior);
EXPECT_EQ(result, LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[0])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[1])));
vp = ve1;
result = mesh.locate_triangle(vp, loc);
EXPECT_EQ(result, LocateTriangleResult::Interior);
EXPECT_EQ(result, LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[1])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[2])));
vp = ve2;
result = mesh.locate_triangle(vp, loc);
EXPECT_EQ(result, LocateTriangleResult::Interior);
EXPECT_EQ(result, LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[2])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[3])));
vp = ve3;
result = mesh.locate_triangle(vp, loc);
EXPECT_EQ(result, LocateTriangleResult::Interior);
EXPECT_EQ(result, LocateTriangleResult::Exterior);
e = mesh.edge(loc);
EXPECT_TRUE((mesh.point(e.node()) == mesh.point(vmap[3])) && (mesh.point(mesh.edge(e.next()).node()) == mesh.point(vmap[0])));
}
......
test/Model Templates/test_SynchronousReluctanceMachine.cpp
View file @ a49cc27a
......@@ -58,9 +58,12 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
EXPECT_TRUE(mesh.edges_are_valid());
EXPECT_TRUE(mesh.edges_are_optimal());
//mesh.MinimumElementQuality = 0.5;
mesh.MinimumElementQuality = 0.5;
//EXPECT_TRUE(mesh.refine());
EXPECT_TRUE(mesh.refine());
EXPECT_TRUE(mesh.edges_are_valid());
EXPECT_TRUE(mesh.edges_are_optimal());
mesh.save_as(SAVE_DIR,"Synchronous_Reluctance_Machine");
}
\ No newline at end of file
test/Model Templates/test_SynchronousReluctanceRotor.cpp
View file @ a49cc27a
......@@ -38,9 +38,15 @@ TEST_F(Synchronous_Reluctance_Rotor, basic_test) {
mesh.create();
EXPECT_TRUE(mesh.edges_are_valid());
EXPECT_TRUE(mesh.edges_are_optimal());
mesh.MinimumElementQuality = 0.5;
EXPECT_TRUE(mesh.refine());
EXPECT_TRUE(mesh.edges_are_valid());
EXPECT_TRUE(mesh.edges_are_optimal());
mesh.save_as(SAVE_DIR,"Synchronous_Reluctance_Rotor");
}
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment