Commit 32dddf26 authored by Pries, Jason's avatar Pries, Jason

Implements second order curved elements and tests on circular domain

parent 27535235
......@@ -22,9 +22,9 @@ public:
bool uniform_discretizaiton() const { return UniformDiscretization; };
double smallest_parametric_edge_length(Mesh &m) const;
double_t smallest_parametric_edge_length(Mesh &m) const;
double queue_uniform(Mesh &m, double delta_min) const;
double_t queue_uniform(Mesh &m, double delta_min) const;
void add_dart_constraint(size_t d) { DartConstraints.push_back(d); }; // TODO: Enforce uniqueness (std::set?)
......
......@@ -28,7 +28,7 @@ public:
size_t tip(Mesh const &m) const;
double delta() const { return std::abs(S1 - S0); };
double_t delta() const { return std::abs(S1 - S0); };
void add_to_queue(Mesh &m) const { Constraint->add_to_queue(m, dart()); };
......@@ -36,8 +36,8 @@ public:
void reverse_dart(size_t d) { ReverseDart = d; };
double S0;
double S1;
double_t S0;
double_t S1;
size_t Self;
......
......@@ -250,6 +250,8 @@ public:
for (size_t k = 0; k != 3; ++k) {
DartConstraint dc = dart_constraint_from_edge(e.Self);
std::shared_ptr<const Curve> cc = dc.curve();
Point point0 = base(e);
Point point1 = tip(e);
......@@ -264,10 +266,16 @@ public:
for (size_t i = 0; i != P + 1; i++) {
double_t s = i / (double_t) (P);
Point pointi = point1 * s + point0 * (1.0 - s);
double_t parami = param1 * s + param0 * (1.0 - s);
data.emplace_back(t, k, i, dc.curve(), parami, pointi);
Point pointi;
if (i == 0 || i == P || cc == nullptr) { // Strict requirement on equality of nodes at the beginning and ending of curves
pointi = point1 * s + point0 * (1.0 - s);
} else {
pointi = cc->point(parami);
}
data.emplace_back(t, k, i, cc, parami, pointi);
}
e = next(e);
......
......@@ -729,7 +729,7 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Periodic_Rotation) {
double_t a_f = -0.5;
double_t c_f = -2.0 * a_f * std::log(2.0) - a_f;
double_t a_h = 2 * a_f;
double_t c_h = -2.0 * a_f * log(2.0);
double_t c_h = -2.0 * a_f * std::log(2.0);
for (size_t i = 0; i != fem->nodes_size(); ++i) {
XY const &n = fem->node(i);
double_t r = std::hypot(n.x(), n.y());
......
......@@ -3,21 +3,21 @@
class Circle : public ::testing::Test {
public:
virtual void SetUp() {
// SetUp_Sketch();
// SetUp_Mesh();
SetUp_Sketch();
SetUp_Mesh();
}
/*
void SetUp_Sketch() {
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>(1.0, 1.0);
auto v3 = sk.new_element<Vertex>(0.0, 1.0);
auto v1 = sk.new_element<Vertex>(+1.0, 0.0);
auto v2 = sk.new_element<Vertex>(0.0, +1.0);
auto v3 = sk.new_element<Vertex>(-1.0, 0.0);
auto v4 = sk.new_element<Vertex>(0.0, -1.0);
boundaries.push_back(sk.new_element<LineSegment>(v0,v1));
boundaries.push_back(sk.new_element<LineSegment>(v1,v2));
boundaries.push_back(sk.new_element<LineSegment>(v2,v3));
boundaries.push_back(sk.new_element<LineSegment>(v3,v0));
boundaries.push_back(sk.new_element<CircularArc>(v1,v2,v0,1.0));
boundaries.push_back(sk.new_element<CircularArc>(v2,v3,v0,1.0));
boundaries.push_back(sk.new_element<CircularArc>(v3,v4,v0,1.0));
boundaries.push_back(sk.new_element<CircularArc>(v4,v1,v0,1.0));
double residual = sk.solve();
EXPECT_LE(residual, FLT_EPSILON);
......@@ -29,7 +29,6 @@ public:
}
void SetUp_Mesh() {
// Create Refineable Mesh
me = Mesh(sk);
me.create();
......@@ -43,14 +42,14 @@ public:
me.save_as(SAVE_DIR, std::string("mesh"));
square_contour = me.select_contour(Point{0.5,0.5});
circle_contour = me.select_contour(Point{0.0,0.0});
}
std::string SAVE_DIR = "./test/output/Integration/test_Square";
std::string SAVE_DIR = "./test/output/Integration/test_Circle/";
std::vector<std::shared_ptr<Curve>> boundaries;
std::shared_ptr<const Contour> square_contour;
std::shared_ptr<const Contour> circle_contour;
Sketch sk;
Mesh me;
......@@ -58,15 +57,15 @@ public:
std::shared_ptr<FiniteElementMeshInterface> fem;
template<size_t P, size_t Q>
auto magnetostatic(size_t b) {
auto magnetostatic() {
fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
Magnetostatic ms{fem};
FunctionArguments fargs;
ms.add_current_density([](FunctionArguments args) { return 1.0; }, {square_contour});
ms.add_current_density([](FunctionArguments args) { return 1.0; }, {circle_contour});
ms.add_magnetic_insulation({boundaries[b]});
ms.add_magnetic_insulation({boundaries[0],boundaries[1],boundaries[2],boundaries[3]});
ms.assemble();
......@@ -79,27 +78,23 @@ public:
return v;
};
void test_boundary0_solution(Oe::VectorXd const &v, double_t tol) {
// v(x,y) = -0.5 * mu_o * y^2 + mu_o * y
for (size_t i = 0; i != fem->nodes().size(); ++i) {
double_t y = fem->node(i).y();
double_t v_expected = -0.5 * mu_0 * y * y + mu_0 * y;
EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
}
}
void test_solution(Oe::VectorXd const &v, double_t tol) {
// v(r) = a * (1 - r^2), a = mu_0 / 4.0;
void test_boundary3_solution(Oe::VectorXd const &v, double_t tol) {
// v(x,y) = -0.5 * mu_o * x^2 + mu_o * x
double_t a = mu_0 / 4.0;
for (size_t i = 0; i != fem->nodes().size(); ++i) {
double_t x = fem->node(i).x();
double_t v_expected = -0.5 * mu_0 * x * x + mu_0 * x;
EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
double_t r = std::hypot(fem->node(i).x(),fem->node(i).y());
double_t v_expected = a * (1.0 - r * r);
EXPECT_NEAR(v(i), v_expected, tol * a);
}
}
*/
};
TEST_F(Circle, magnetostatic) {
test_solution(magnetostatic<1,1>(), 0.01);
test_solution(magnetostatic<1,2>(), 0.01);
test_solution(magnetostatic<2,2>(), 0.0001); // Does not obtain FLT_EPSILON accuracy because of boundary discretization error
FAIL();
// TODO: Implement spatially dependent current distribution forcing function
// TODO: Implement tests on circle with sinusoidal current distribution (cylindrical coordinates)
......
......@@ -45,7 +45,7 @@ public:
square_contour = me.select_contour(Point{0.5,0.5});
}
std::string SAVE_DIR = "./test/output/Integration/test_Square";
std::string SAVE_DIR = "./test/output/Integration/test_Square/";
std::vector<std::shared_ptr<Curve>> boundaries;
......
Markdown is supported
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