Skip to content

GitLab

Commit a8054788 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit 

Elaboration of TwoRegionHex magnetostatic_forcing tests
parent 9cd34359
Hide whitespace changes
Inline Side-by-side
test/Physics/test_FiniteElementMesh.cpp
View file @ a8054788
......@@ -15,11 +15,11 @@ public:
reg.push_back(Region<2>(std::vector<size_t>{0}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0,1}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0,2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1,2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0, 1}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0, 2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1, 2}));
femesh = FiniteElementMesh<2,1>(node, tri, reg, bnd);
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
......@@ -27,7 +27,7 @@ public:
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
FiniteElementMesh<2,1> femesh;
FiniteElementMesh<2, 1> femesh;
};
TEST_F(MasterTriangle, mass_matrix) {
......@@ -111,7 +111,7 @@ public:
tri.push_back(Triangle<1>(5, 1_zu, 4_zu, 0_zu)); // x-axis reflection
tri.push_back(Triangle<1>(6, 4_zu, 3_zu, 0_zu)); // x=y reflection
reg.push_back(Region<2>(std::vector<size_t>{0,1,2,3})); // Triangles 4,5,6 are duplicates
reg.push_back(Region<2>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
......@@ -242,11 +242,11 @@ public:
reg.push_back(Region<2>(std::vector<size_t>{0}));
reg.push_back(Region<2>(std::vector<size_t>{1}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0,1}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0,2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1,2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1,3}));
bnd.push_back(Boundary<2>(std::vector<size_t>{2,3}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0, 1}));
bnd.push_back(Boundary<2>(std::vector<size_t>{0, 2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1, 2}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1, 3}));
bnd.push_back(Boundary<2>(std::vector<size_t>{2, 3}));
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
......@@ -364,26 +364,26 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
ms.add_current_density(ff, std::vector<size_t>{0});
ms.calculate_forcing(f);
EXPECT_NEAR(f(0),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(1),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(2),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(3),0.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(2), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(3), 0.0 / 6.0, FLT_EPSILON);
ms.remove_current_density(0);
ms.add_current_density(ff, std::vector<size_t>{1});
ms.calculate_forcing(f);
EXPECT_NEAR(f(0),0.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(1),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(2),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(3),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(0), 0.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(2), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(3), 1.0 / 6.0, FLT_EPSILON);
ms.remove_current_density(0);
ms.add_current_density(ff, std::vector<size_t>{0,1});
ms.add_current_density(ff, std::vector<size_t>{0, 1});
ms.calculate_forcing(f);
EXPECT_NEAR(f(0),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(1),2.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(2),2.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(3),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(1), 2.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(2), 2.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(f(3), 1.0 / 6.0, FLT_EPSILON);
ms.remove_current_density(0);
/*
......@@ -467,9 +467,9 @@ public:
tri.push_back(Triangle<1>(4, 0_zu, 5_zu, 6_zu));
tri.push_back(Triangle<1>(5, 0_zu, 6_zu, 1_zu));
reg.push_back(Region<2>(std::vector<size_t>{0,1,2,3,4,5}));
reg.push_back(Region<2>(std::vector<size_t>{0, 1, 2, 3, 4, 5}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1,2,3,4,5,6}));
bnd.push_back(Boundary<2>(std::vector<size_t>{1, 2, 3, 4, 5, 6}));
tri.push_back(Triangle<1>(6, 1_zu, 7_zu, 2_zu));
tri.push_back(Triangle<1>(7, 2_zu, 7_zu, 8_zu));
......@@ -489,9 +489,9 @@ public:
tri.push_back(Triangle<1>(16, 6_zu, 12_zu, 1_zu));
tri.push_back(Triangle<1>(17, 1_zu, 12_zu, 7_zu));
reg.push_back(Region<2>(std::vector<size_t>{6,7,8,9,10,11,12,13,14,15,16,17}));
reg.push_back(Region<2>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17}));
bnd.push_back(Boundary<2>(std::vector<size_t>{7,8,9,10,11,12}));
bnd.push_back(Boundary<2>(std::vector<size_t>{7, 8, 9, 10, 11, 12}));
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
......@@ -505,8 +505,10 @@ public:
};
TEST_F(TwoRegionHex, magnetostatic_forcing) {
Magnetostatic<2, 1, 2, Variable::A> ms{femesh};
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
// Test without boundary conditions
ms.build_matrices();
auto ff = [](double t) { return 1.0; };
......@@ -518,23 +520,21 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
double nodal_current{sqrt(3.0) / 4.0 / 3.0};
EXPECT_NEAR(f(0), nodal_current * 6.0, FLT_EPSILON);
EXPECT_NEAR(f(1), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(2), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(3), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(4), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(5), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(6), nodal_current * 2.0, FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
}
auto J = ms.init_matrix();
auto v = ms.init_vector();
auto r = ms.init_vector();
v.setZero();
ms.linearize(J,v,r,f);
ms.linearize(J, v, r, f);
//std::cout << J << std::endl;
//std::cout << r << std::endl;
// Test with boundary conditions
ms.apply_conditions();
J = ms.init_matrix();
......@@ -545,32 +545,81 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
ms.calculate_forcing(f);
EXPECT_NEAR(f(0), nodal_current * 6.0, FLT_EPSILON);
EXPECT_NEAR(f(1), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(2), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(3), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(4), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(5), nodal_current * 2.0, FLT_EPSILON);
EXPECT_NEAR(f(6), nodal_current * 2.0, FLT_EPSILON);
ms.linearize(J,v,r,f);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
}
ms.linearize(J, v, r, f);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
//std::cout << v << std::endl << std::endl;
// Verify equation is solved
ms.linearize(J, v, r, f);
for (size_t i = 0; i != r.size(); ++i) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
// Fundamental solution is 0.25*r^2 + C in forced region (0)
EXPECT_NEAR(v(0), v(1) + 0.25, FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
//std::cout << ms.derivatives().dx(0).transpose() * v << std::endl << std::endl;
//std::cout << ms.derivatives().dx(1).transpose() * v << std::endl << std::endl;
//std::cout << ms.derivatives().dx(2).transpose() * v << std::endl << std::endl;
// Test flux-density values
// Values are exact in forced region (0)
// Values are approximated in unforced region (1)
Eigen::VectorXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::VectorXd By = -ms.derivatives().dx(0).transpose() * v;
//std::cout << ms.derivatives().dy(0).transpose() * v << std::endl << std::endl;
//std::cout << ms.derivatives().dy(1).transpose() * v << std::endl << std::endl;
//std::cout << ms.derivatives().dy(2).transpose() * v << std::endl << std::endl;
Eigen::VectorXd Bmag(By.size());
Eigen::VectorXd Bang(By.size());
//std::cout << v << std::endl;
//std::cout << J * v << std::endl;
//std::cout << f << std::endl;
//std::cout << J * v - f << std::endl;
for (size_t i = 0; i != By.size(); ++i) {
Bmag(i) = hypot(Bx(i), By(i));
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
}
double x_q = (1.0 + cos(M_PI / 3.0)) / 3.0;
double y_q = sin(M_PI / 3.0) / 3.0;
double r_q = hypot(x_q, y_q);
double B_q = r_q / 2.0;
for (size_t i = 0; i != 6; ++i) {
EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
double ang = 120.0 + 60.0 * i;
ang -= 360.0 * (ang > 180.0);
EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
}
x_q = (1.0 + cos(M_PI / 3.0)) / 2.0;
y_q = sin(M_PI / 3.0) / 2.0;
r_q = 2.0 - hypot(x_q, y_q);
B_q = v(1) / r_q;
for (size_t i = 6; i < 18; i += 2) {
EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
double ang = 120.0 + 60.0 * (i - 6) / 2.0;
ang -= 360.0 * (ang > 180.0);
EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
}
x_q = (2.0 * cos(M_PI / 6.0) + 2.0 * cos(M_PI / 6.0 * 3.0)) / 2.0;
y_q = (2.0 * sin(M_PI / 6.0) + 2.0 * sin(M_PI / 6.0 * 3.0)) / 2.0;
r_q = hypot(x_q, y_q) - 1.0;
B_q = v(1) / r_q;
for (size_t i = 7; i < 18; i += 2) {
EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
double ang = 120.0 + 60.0 * (i - 6) / 2.0;
ang -= 360.0 * (ang > 180.0);
EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
}
}
\ No newline at end of file
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