Commit 14abfb38 authored by JasonPries's avatar JasonPries

Adds DiscontinuousBoundary

This is a new route for implementation of relative motion. Renames Boundary to DiscreteBoundary (to differentiate from Sketch::ContinuousBoundary)
parent 9cdbe1be
......@@ -24,6 +24,8 @@ enum class InsertPointResult {
class Mesh { // TODO: Namespaces. Also, need to rename a bunch of things to be more consistent and clear w.r.t. how various arrays are accessed (e.g. by edge/dart index versus plain index)
public:
Mesh() {};
Mesh(Sketch &s);
friend class BoundaryConstraint;
......
......@@ -5,7 +5,7 @@ set(SOURCE_FILES
#./src/PhysicsCommon.h
./src/Boundary.h ./src/Boundary.cpp
./src/DiscreteBoundary.h ./src/DiscreteBoundary.cpp
./src/BoundaryCondition.h ./src/BoundaryCondition.cpp
......
#ifndef OERSTED_PHYSICS_HPP
#define OERSTED_PHYSICS_HPP
//#include "../src/PhysicsCommon.h"
//#include "../src/Discretization.h"
//#include "../src/Operator.h"
#include "../src/Boundary.h"
#include "../src/DiscreteBoundary.h"
#include "../src/BoundaryCondition.h"
#include "../src/FiniteElementMesh.h"
#include "../src/Forcing.h"
......@@ -14,7 +10,6 @@
#include "../src/Node.h"
#include "../src/Physics.h"
#include "../src/PhysicsConstants.h"
#include "Quadrature/src/TriangleQuadrature.h"
#include "../src/Region.h"
#include "../src/Triangle.h"
......
......@@ -7,8 +7,8 @@
#include <Eigen/Sparse>
#include "DiscreteBoundary.h"
#include "Quadrature.hpp"
#include "MatrixGroup.h"
#include "FiniteElementMesh.h"
......@@ -31,7 +31,7 @@ class ZeroDirichlet : public BoundaryCondition {
template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};
ZeroDirichlet(std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundaries) : Boundaries{boundaries} {};
ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<2, ElementOrder> const &fem) {
for (auto const &b : boundaries) {
......@@ -51,7 +51,7 @@ public:
};
protected:
std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
};
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
......@@ -129,6 +129,9 @@ public:
SlidingInterface(std::vector<size_t> first_vars, std::vector<size_t> second_vars, std::vector<double> vals, bool antiperiodic = false)
: First{first_vars}, Second{second_vars}, Value{vals}, Antiperiodic{antiperiodic} {};
SlidingInterface(std::shared_ptr<DiscontinuousBoundary<2>> db, bool antiperiodic = false)
: First{db->nodes()}, Second{db->dnodes()}, Value{db->values()}, Antiperiodic{antiperiodic} {};
SlidingInterface &operator++() {
std::rotate(Second.begin(), Second.begin() + 1, Second.end());
if (Antiperiodic) {
......
#include "DiscreteBoundary.h"
......@@ -7,35 +7,65 @@
#include <Sketch.hpp>
template<size_t Dimension>
class Boundary {
class DiscreteBoundary {
};
template<size_t Dimension>
class DiscontinuousBoundary : public DiscreteBoundary<Dimension> {
};
template<>
class Boundary<2> { // TODO: Rename DiscreteBoundary
class DiscreteBoundary<2> {
public:
Boundary() {};
DiscreteBoundary() {};
Boundary(std::vector<size_t> &&nodes) : CurvePtr{nullptr}, Nodes{nodes} {};
DiscreteBoundary(std::vector<size_t> nodes) : CurvePtr{nullptr}, Nodes{nodes} {};
Boundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> &&nodes) : CurvePtr{cptr}, Nodes{nodes} {};
DiscreteBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes) : CurvePtr{cptr}, Nodes{nodes} {};
size_t size() const { return Nodes.size(); };
auto &nodes() { return Nodes; };
std::vector<size_t> nodes() { return Nodes; };
auto const &nodes() const { return Nodes; };
std::vector<size_t> const &nodes() const { return Nodes; };
auto &node(size_t i) { return Nodes[i]; };
size_t node(size_t i) { return Nodes[i]; };
auto const &node(size_t i) const { return Nodes[i]; };
size_t const &node(size_t i) const { return Nodes[i]; };
auto curve() const { return CurvePtr; };
std::shared_ptr<Curve const> curve() const { return CurvePtr; };
protected:
std::shared_ptr<Curve const> CurvePtr;
std::vector<size_t> Nodes;
};
template<>
class DiscontinuousBoundary<2> : public DiscreteBoundary<2> {
public:
DiscontinuousBoundary() {};
DiscontinuousBoundary(std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals) : DiscreteBoundary{nodes}, DiscontinuousNodes{dnodes}, Values{vals} {};
DiscontinuousBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals) : DiscreteBoundary{cptr, nodes}, DiscontinuousNodes{dnodes}, Values{vals} {};
std::vector<size_t> dnodes() { return DiscontinuousNodes; };
std::vector<size_t> const &dnodes() const { return DiscontinuousNodes; };
std::vector<double_t> values() { return Values; };
std::vector<double_t> const &values() const { return Values; };
size_t dnode(size_t i) { return DiscontinuousNodes[i]; };
size_t const &dnode(size_t i) const {return DiscontinuousNodes[i]; };
protected:
std::vector<size_t> DiscontinuousNodes;
std::vector<double_t> Values;
};
/*
template <size_t Dimension>
class BoundaryPair {};
......
......@@ -29,7 +29,7 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
Boundaries.reserve(m.size_boundary_constraints());
for (size_t i = 0; i!= m.size_boundary_constraints(); ++i) {
Boundaries.push_back(std::make_shared<Boundary<2>>(m.curve(i), m.boundary_nodes(i)));
Boundaries.push_back(std::make_shared<DiscreteBoundary<2>>(m.curve(i), m.boundary_nodes(i)));
}
Regions.reserve(m.size_contours());
......
......@@ -6,7 +6,7 @@
#include "Mesh.hpp"
#include "Boundary.h"
#include "DiscreteBoundary.h"
#include "MatrixGroup.h"
#include "Node.h"
#include "Triangle.h"
......@@ -21,8 +21,8 @@ class FiniteElementMesh<2, Order> {
public:
FiniteElementMesh() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<Boundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {
auto comp = [](std::shared_ptr<Boundary<2>> const &x, std::shared_ptr<Boundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {
auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, std::shared_ptr<DiscreteBoundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
std::sort(Boundaries.begin(), Boundaries.end(), comp);
};
......@@ -40,15 +40,15 @@ public:
auto &boundary(size_t i) { return Boundaries[i]; };
std::shared_ptr<Boundary<2>> boundary(std::shared_ptr<Curve const> const &c) const {
auto comp = [](std::shared_ptr<Boundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
std::shared_ptr<DiscreteBoundary<2>> boundary(std::shared_ptr<Curve const> const &c) const {
auto comp = [](std::shared_ptr<DiscreteBoundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
auto result = std::lower_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), comp);
if (result != Boundaries.end()) {
return *result;
} else {
std::shared_ptr<Boundary<2>> null;
std::shared_ptr<DiscreteBoundary<2>> null;
return null;
}
};
......@@ -167,7 +167,7 @@ protected:
std::vector<Triangle<Order>> Triangles;
std::vector<std::shared_ptr<Region<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
std::vector<std::shared_ptr<DiscreteBoundary<2>>> Boundaries;
};
#endif //OERSTED_FINITEELEMENTMESH_H
......@@ -171,7 +171,7 @@ public:
ForcingCondtions.erase(ForcingCondtions.begin() + i);
}
void add_magnetic_insulation(std::vector<std::shared_ptr<Boundary<2>>> boundaries) { // TODO: should boundaries be unique_ptrs?
void add_magnetic_insulation(std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundaries) { // TODO: should boundaries be unique_ptrs?
BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
}
......@@ -189,6 +189,17 @@ public:
template <typename... Args>
auto add_sliding_interface(Args&&... args) {
/*
* TODO: I'm not sure why forwarding was used here
* Instead, accept a std::shared_ptr<Curve>
* Call Domain.make_discontinuous(std::shared_ptr<Curve>)
* Construct sliding_interface with the returned std::shared_ptr<DiscontinuousBoundary>
*
* TODO: Implement Domain.make_discontinuous(std::shared_ptr<Curve>)
* Try casting the std::shared_ptr<DiscreteBoundary> to std::shared_ptr<DiscontinuousBoundary>
* If successful, return the std::shared_ptr<DiscontinuousBoundary>
* Else, construct DiscontinuousBoundary, replace the existing std::shared_ptr<DiscreteBoundary>, and return the new std::shared_ptr<DiscontinuousBoundary>
*/
auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...);
BoundaryConditions.push_back(condition);
return condition;
......
......@@ -2,6 +2,185 @@
std::string SAVE_DIR = "./test/output/Integration/";
class Library_Integration__Quarter_Circle : public ::testing::Test {
public:
virtual void SetUp() {
SetUp_Sketch();
SetUp_Mesh();
}
void SetUp_Sketch() {
v0 = sk.new_element<Vertex>(0.0, 0.0);
v1 = sk.new_element<Vertex>(1.0, 0.0);
v2 = sk.new_element<Vertex>(2.0, 0.0);
v3 = sk.new_element<Vertex>(sqrt(3.0)/2.0, 0.5);
v4 = sk.new_element<Vertex>(sqrt(3.0), 1.0);
f0 = sk.new_element<Fixation>(v0);
l01 = sk.new_element<LineSegment>(v0, v1);
l12 = sk.new_element<LineSegment>(v1, v2);
l04 = sk.new_element<LineSegment>(v0, v4, true);
h01 = sk.new_element<Horizontal>(l01);
h12 = sk.new_element<Horizontal>(l12);
a0103 = sk.new_element<Angle>(l01, l04, 30.0);
c013 = sk.new_element<CircularArc>(v1, v3, v0, 1.0);
c024 = sk.new_element<CircularArc>(v2, v4, v0, 2.0);
r013 = sk.new_element<Radius>(c013, 1.0);
r024 = sk.new_element<Radius>(c024, 2.0);
std::vector<std::shared_ptr<Curve const>> vec{l01, l12, c013, c024};
m04 = sk.new_element<MirrorCopy>(vec, l04);
double residual = sk.solve();
EXPECT_LE(residual, FLT_EPSILON);
bool result = sk.build();
ASSERT_TRUE(result);
sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, "quarter_circle_mirror_copy_sketch");
periodic_boundary = sk.select_periodic_boundary_pairs(v0, 60.0);
{
for (auto &bp : periodic_boundary) {
if (bp.curve0().get() == l01.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l01, v0, 60.0), MatchOrientation::Reverse);
} else if (bp.curve0().get() == l12.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l12, v0, 60.0), MatchOrientation::Reverse);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found");
}
}
}
radial_boundary = sk.select_radial_boundary(v0, 2.0);
{
for (auto &c : radial_boundary) {
auto cc = std::dynamic_pointer_cast<CircularArc const>(c);
if (cc) {
EXPECT_NEAR(cc->radius(), 2.0, FLT_EPSILON);
EXPECT_EQ(cc->center().get(), v0.get());
} else {
GTEST_NONFATAL_FAILURE_("dynamic_cast of Curve to CircularArc failed");
}
}
}
v1p = sk.select_periodic_vertex(v1, v0, 60.0);
{
ASSERT_TRUE(v1p); // Make sure that returned vertex pointer is not null
double_t x0 = v0->x();
double_t y0 = v0->y();
double_t x1 = v1->x();
double_t y1 = v1->y();
double_t x2 = v1p->x();
double_t y2 = v1p->y();
double_t dx = x1 - x0;
double_t dy = y1 - y0;
double_t dr = std::hypot(dx, dy);
double_t da = std::atan2(dy, dx);
da += M_PI * 1.0 / 3.0;
x1 = x0 + dr * std::cos(da);
y1 = y0 + dr * std::sin(da);
EXPECT_NEAR(x1, x2, dr * FLT_EPSILON);
EXPECT_NEAR(y1, y2, dr * FLT_EPSILON);
}
}
void SetUp_Mesh() {
// Create Refineable Mesh
me = Mesh(sk);
me.create();
me.MaximumElementSize = 0.1;
me.MinimumElementSize = 0.01;
me.MinimumElementQuality = 0.5;
for (auto &pb : periodic_boundary) {
auto bc0 = me.boundary_constraint(pb.curve0());
if (bc0) {
bc0->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
auto bc1 = me.boundary_constraint(pb.curve1());
if (bc1) {
bc1->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
}
bool result = me.refine();
ASSERT_TRUE(result);
me.save_as(SAVE_DIR, std::string("quarter_circle_mirror_copy_mesh"));
// Create FiniteElementMesh
fem = FiniteElementMesh<2,1>(me);
for (auto & pb : periodic_boundary) {
auto const &nodes0 = fem.boundary(pb.curve0())->nodes();
auto const &nodes1 = fem.boundary(pb.curve1())->nodes();
for (size_t i = 0; i != nodes0.size(); ++i) {
auto const &n0 = fem.node(nodes0[i]);
auto const &n1 = fem.node(nodes1[nodes0.size() - 1 - i]);
EXPECT_NEAR(atan2(n0.y(), n0.x()), 0.0, FLT_EPSILON);
if (&n0 != &n1) {
EXPECT_NEAR(atan2(n1.y(), n1.x()), M_PI / 3.0, FLT_EPSILON);
}
}
}
}
Sketch sk;
Mesh me;
FiniteElementMesh<2,1> fem;
// Sketch Variables
std::shared_ptr<Vertex> v0;
std::shared_ptr<Vertex> v1;
std::shared_ptr<Vertex> v2;
std::shared_ptr<Vertex> v3;
std::shared_ptr<Vertex> v4;
std::shared_ptr<Fixation> f0;
std::shared_ptr<LineSegment> l01 ;
std::shared_ptr<LineSegment> l12;
std::shared_ptr<LineSegment> l04;
std::shared_ptr<Horizontal> h01;
std::shared_ptr<Horizontal> h12;
std::shared_ptr<Angle> a0103;
std::shared_ptr<CircularArc> c013;
std::shared_ptr<CircularArc> c024;
std::shared_ptr<Radius> r013;
std::shared_ptr<Radius> r024;
std::shared_ptr<MirrorCopy> m04;
decltype(sk.select_periodic_boundary_pairs(v0, 90.0)) periodic_boundary;
decltype(sk.select_radial_boundary(v0, 2.0)) radial_boundary;
decltype(sk.select_periodic_vertex(v1, v0, 90.0)) v1p;
};
TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
// Create Sketch
Sketch sk;
......@@ -45,7 +224,7 @@ TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
// Convert to FiniteElementMesh
FiniteElementMesh<2, 1> fem{me};
for (std::shared_ptr<Boundary<2>> const &b : fem.boundaries()) {
for (std::shared_ptr<DiscreteBoundary<2>> const &b : fem.boundaries()) {
double_t a0{-2 * M_PI};
for (size_t i : b->nodes()) {
XY const &p = fem.node(i);
......@@ -138,142 +317,7 @@ TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
}
}
TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
GTEST_NONFATAL_FAILURE_("TODO: antiperiodic boundary conditions");
// Create Sketch
Sketch sk;
auto v0 = sk.new_element<Vertex>(0.0, 0.0);
auto v1 = sk.new_element<Vertex>(1.0, 0.0);
auto v2 = sk.new_element<Vertex>(2.0, 0.0);
auto v3 = sk.new_element<Vertex>(M_SQRT1_2, M_SQRT1_2);
auto v4 = sk.new_element<Vertex>(M_SQRT2, M_SQRT2);
auto f0 = sk.new_element<Fixation>(v0);
auto l01 = sk.new_element<LineSegment>(v0, v1);
auto l12 = sk.new_element<LineSegment>(v1, v2);
auto l04 = sk.new_element<LineSegment>(v0, v4, true);
auto h01 = sk.new_element<Horizontal>(l01);
auto h12 = sk.new_element<Horizontal>(l12);
auto a0103 = sk.new_element<Angle>(l01, l04, 45.0);
auto c013 = sk.new_element<CircularArc>(v1, v3, v0, 1.0);
auto c024 = sk.new_element<CircularArc>(v2, v4, v0, 2.0);
auto r013 = sk.new_element<Radius>(c013, 1.0);
auto r024 = sk.new_element<Radius>(c024, 2.0);
std::vector<std::shared_ptr<Curve const>> vec{l01, l12, c013, c024};
auto m04 = sk.new_element<MirrorCopy>(vec, l04);
double residual = sk.solve();
EXPECT_LE(residual, FLT_EPSILON);
bool result = sk.build();
ASSERT_TRUE(result);
sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, "quarter_circle_mirror_copy_sketch");
auto periodic_boundary = sk.select_periodic_boundary_pairs(v0, 90.0);
{
for (auto &bp : periodic_boundary) {
if (bp.curve0().get() == l01.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l01, v0, 90.0), MatchOrientation::Reverse);
} else if (bp.curve0().get() == l12.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l12, v0, 90.0), MatchOrientation::Reverse);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found");
}
}
}
auto radial_boundary = sk.select_radial_boundary(v0, 2.0);
{
for (auto &c : radial_boundary) {
auto cc = std::dynamic_pointer_cast<CircularArc const>(c);
if (cc) {
EXPECT_NEAR(cc->radius(), 2.0, FLT_EPSILON);
EXPECT_EQ(cc->center().get(), v0.get());
} else {
GTEST_NONFATAL_FAILURE_("dynamic_cast of Curve to CircularArc failed");
}
}
}
auto v1p = sk.select_periodic_vertex(v1, v0, 90.0);
{
ASSERT_TRUE(v1p); // Make sure that returned vertex pointer is not null
double_t x0 = v0->x();
double_t y0 = v0->y();
double_t x1 = v1->x();
double_t y1 = v1->y();
double_t x2 = v1p->x();
double_t y2 = v1p->y();
double_t dx = x1 - x0;
double_t dy = y1 - y0;
double_t dr = std::hypot(dx, dy);
double_t da = std::atan2(dy, dx);
da += M_PI / 2.0;
x1 = x0 + dr * std::cos(da);
y1 = y0 + dr * std::sin(da);
EXPECT_NEAR(x1, x2, dr * FLT_EPSILON);
EXPECT_NEAR(y1, y2, dr * FLT_EPSILON);
}
// Create Mesh
Mesh me{sk};
me.create();
me.MaximumElementSize = 0.1;
me.MinimumElementSize = 0.01;
me.MinimumElementQuality = 0.5;
for (auto &pb : periodic_boundary) {
auto bc0 = me.boundary_constraint(pb.curve0());
if (bc0) {
bc0->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
auto bc1 = me.boundary_constraint(pb.curve1());
if (bc1) {
bc1->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
}
result = me.refine();
ASSERT_TRUE(result);
me.save_as(SAVE_DIR, std::string("quarter_circle_mirror_copy_mesh"));
// Create FiniteElementMesh
FiniteElementMesh<2,1> fem{me};
for (auto & pb : periodic_boundary) {
auto const &nodes0 = fem.boundary(pb.curve0())->nodes();
auto const &nodes1 = fem.boundary(pb.curve1())->nodes();
for (size_t i = 0; i != nodes0.size(); ++i) {
auto const &n0 = fem.node(nodes0[i]);
auto const &n1 = fem.node(nodes1[nodes0.size() - 1 - i]);
EXPECT_NEAR(n0.y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(n1.x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(n0.x(), n1.y(), FLT_EPSILON);
}
}
TEST_F(Library_Integration__Quarter_Circle, Uniform_Current_Density__Periodic) {
// Create Physics
Magnetostatic<2, 1, 1, FieldVariable::A> msph{fem};
......@@ -322,7 +366,7 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
// A = a_h * log(r) + c_h in homogeneous outer region
Eigen::VectorXd vv = msph.recover_boundary(v);
fem.write_scalar(vv, SAVE_DIR, std::string("quarter_circle_mvp"));
fem.write_scalar(vv, SAVE_DIR, std::string("uniform_current_density__periodic"));
double_t a_f = -0.5;