Commit e21f3e59 authored by JasonPries's avatar JasonPries
Browse files

Adds tests to solution and flux density of Quarter_Circle_Mirror_Copy_Uniform_Current_Density

parent 2aee9a7b
...@@ -76,12 +76,16 @@ public: ...@@ -76,12 +76,16 @@ public:
DerivativeMatrixGroup<QuadraturePoints> const &derivatives() { return Derivatives; }; DerivativeMatrixGroup<QuadraturePoints> const &derivatives() { return Derivatives; };
Eigen::VectorXd recover_boundary(Eigen::VectorXd const &v) const { return BoundaryConditionMatrix.transpose() * v; };
protected: protected:
FiniteElementMesh<Dimension, ElementOrder> &Domain; FiniteElementMesh<Dimension, ElementOrder> &Domain;
DiagonalMatrixGroup<QuadraturePoints> ElementWeights; DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
SparseMatrixGroup<QuadraturePoints> Basis; SparseMatrixGroup<QuadraturePoints> Basis;
DerivativeMatrixGroup<QuadraturePoints> Derivatives; DerivativeMatrixGroup<QuadraturePoints> Derivatives;
Eigen::SparseMatrix<double> BoundaryConditionMatrix;
}; };
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, FieldVariable V> template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, FieldVariable V>
...@@ -98,6 +102,8 @@ public: ...@@ -98,6 +102,8 @@ public:
using PhysicsData<2, ElementOrder, QuadratureOrder>::weights; using PhysicsData<2, ElementOrder, QuadratureOrder>::weights;
using PhysicsData<2, ElementOrder, QuadratureOrder>::basis; using PhysicsData<2, ElementOrder, QuadratureOrder>::basis;
using PhysicsData<2, ElementOrder, QuadratureOrder>::recover_boundary;
Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.size_nodes()}, Elements{d.size_elements()} {}; Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.size_nodes()}, Elements{d.size_elements()} {};
Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); }; Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); };
...@@ -238,11 +244,11 @@ public: ...@@ -238,11 +244,11 @@ public:
size_t rows{cols - num_erased}; size_t rows{cols - num_erased};
// TODO: END // TODO: END
Eigen::SparseMatrix<double> bc_matrix(rows, cols); BoundaryConditionMatrix.resize(rows, cols);
bc_matrix.setFromTriplets(triplets.begin(), triplets.end()); BoundaryConditionMatrix.setFromTriplets(triplets.begin(), triplets.end());
Basis.transform(bc_matrix); Basis.transform(BoundaryConditionMatrix);
Derivatives.transform(bc_matrix); Derivatives.transform(BoundaryConditionMatrix);
Unknowns = rows; Unknowns = rows;
}; };
...@@ -251,6 +257,7 @@ protected: ...@@ -251,6 +257,7 @@ protected:
using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives; using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives;
using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights; using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights;
using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis; using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis;
using PhysicsData<2, ElementOrder, QuadratureOrder>::BoundaryConditionMatrix;
private: private:
size_t Unknowns; size_t Unknowns;
......
...@@ -277,7 +277,8 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) { ...@@ -277,7 +277,8 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
// TODO: Need a method for selecting Contour from Sketch and selecting FiniteElementMesh Region<2> given a contour // TODO: Need a method for selecting Contour from Sketch and selecting FiniteElementMesh Region<2> given a contour
// TODO: Add std::shared_ptr<Contour const> to Region<2> // TODO: Add std::shared_ptr<Contour const> to Region<2>
//msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(1)}); GTEST_NONFATAL_FAILURE_("TODO: region selection");
msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(1)});
msph.add_magnetic_insulation(radial_boundary); msph.add_magnetic_insulation(radial_boundary);
...@@ -314,7 +315,31 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) { ...@@ -314,7 +315,31 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON); EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
} }
// Test solution:
// A = a_f * r^2 + c_f in forced inner region
// A = a_h * log(r) + c_h in homogeneous outer region
Eigen::VectorXd vv = msph.recover_boundary(v);
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);
for(size_t i = 0; i!= fem.size_nodes(); ++i) {
XY const &n = fem.node(i);
double_t r = std::hypot(n.x(), n.y());
if (r <= 1.0) {
double_t val = a_f * r * r + c_f;
EXPECT_NEAR(val, vv(i), c_f * 0.01);
} else {
double_t val = a_h * log(r) + c_h;
EXPECT_NEAR(val, vv(i), c_h * 0.01);
}
}
// Test flux-density values // Test flux-density values
// B_r = 2 * a_f * r in forced inner region
// B_r = a_h / r in homogenous outer region
// TODO: The correct way to test may be to calculate average flux in the element from analytical solution
Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * v; Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * v; Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * v;
...@@ -326,13 +351,22 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) { ...@@ -326,13 +351,22 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI; Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
} }
std::cout << Bmag << std::endl << std::endl; std::vector<std::vector<XY>> qp = fem.quadrature_points<1>();
std::cout << Bang << std::endl << std::endl; for (size_t i = 0; i != qp.size(); ++i) {
double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());
std::cout << f << std::endl << std::endl; if (r <= 1) {
double_t val = std::abs(2.0 * a_f * r);
EXPECT_NEAR(Bmag(i), val, std::abs(2 * a_f * 0.05));
} else {
double_t val = std::abs(a_h / r);
EXPECT_NEAR(Bmag(i), val, std::abs(a_h* 0.05));
}
std::cout << v << std::endl << std::endl; double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;
GTEST_NONFATAL_FAILURE_("Solution needs testing"); if (Bang(i) < 0.0) { Bang(i) += 360.0; };
EXPECT_NEAR(Bang(i), a, 90 * 0.05);
}
} }
\ No newline at end of file
...@@ -1128,7 +1128,6 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) { ...@@ -1128,7 +1128,6 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON); EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
} }
// Fundamental solution is 0.25*r^2 + C in forced region (0)
for (size_t i = 0; i != 7; ++i) { for (size_t i = 0; i != 7; ++i) {
EXPECT_NEAR(v(i), v_expected[i], FLT_EPSILON); EXPECT_NEAR(v(i), v_expected[i], FLT_EPSILON);
} }
......
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