Commit 0c9df1d0 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit

Initial implementation and testing of relative motion with a SlidingInterface
parent 5d7843f2
......@@ -24,8 +24,7 @@ public:
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet : public BoundaryCondition {
};
class ZeroDirichlet : public BoundaryCondition {};
template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
......@@ -46,10 +45,10 @@ protected:
std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
};
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition : public BoundaryCondition {
};
class PeriodicBoundaryCondition : public BoundaryCondition {};
template<size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
......@@ -80,5 +79,42 @@ protected:
bool Antiperiodic;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface : public BoundaryCondition {};
template<size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
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& operator++() {
std::rotate(Second.begin(), Second.begin() + 1, Second.end());
if (Antiperiodic) {
std::rotate(Value.begin(), Value.begin() + 1, Value.end());
*(Value.end()) = -*(Value.end());
}
}
// TODO: SlidingInterface& opertor--() {...};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
for(size_t i = 0; i!=Value.size();++i){
triplets.push_back(Eigen::Triplet<double>(First[i],Second[i],Value[i]));
}
}
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (size_t i = 0; i!=Second.size();++i){
index.insert(Second[i]);
}
}
protected:
std::vector<size_t> First;
std::vector<size_t> Second;
std::vector<double> Value;
bool Antiperiodic;
};
#endif //OERSTED_BOUNDARYCONDITION_H
......@@ -15,6 +15,10 @@
#include "QuadratureRule.h"
#include "Triangle.h"
// TODO: These classes should be generic (e.g. CurlCurl, ScalarLaplacian, VectorLaplacian)
// TODO: Higher level classes should implement interfaces with appropriate physics grammar
// TODO: Implement material property dependencies using pointers to virtual member functions
enum class Variable {
MagneticVectorPotential = 1,
A = 1,
......@@ -43,8 +47,8 @@ public:
protected:
double Time;
std::vector<std::unique_ptr<Forcing>> ForcingCondtions;
std::vector<std::unique_ptr<BoundaryCondition>> BoundaryConditions;
std::vector<std::shared_ptr<Forcing>> ForcingCondtions;
std::vector<std::shared_ptr<BoundaryCondition>> BoundaryConditions;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
......@@ -155,7 +159,7 @@ public:
}
void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<Region<2>>> regions) {
ForcingCondtions.push_back(std::make_unique<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
ForcingCondtions.push_back(std::make_shared<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
}
void remove_current_density(size_t i) {
......@@ -163,11 +167,18 @@ public:
}
void add_magnetic_insulation(std::vector<std::shared_ptr<Boundary<2>>> boundaries) { // TODO: should boundaries be unique_ptrs?
BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
}
void add_periodic_boundary(std::vector<BoundaryMap> map, bool antiperiodic) {
BoundaryConditions.push_back(std::make_unique<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
}
template <typename... Args>
auto add_sliding_interface(Args&&... args) {
auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...);
BoundaryConditions.push_back(condition);
return condition;
}
void apply_conditions() {
......
......@@ -859,7 +859,7 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_periodic) {
std::vector<double> b_mag_expected(6, 3.86994e-7);
std::vector<double> b_ang_expected{112.5, 157.5, 112.5, 112.5, 157.5, 157.5};
for(size_t i = 0;i!=6;++i) {
for (size_t i = 0; i != 6; ++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-12);
EXPECT_NEAR(Bang[i], b_ang_expected[i], 1e-1);
}
......@@ -914,7 +914,7 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
}
std::vector<double> v_expected{0.0,6.60365e-8};
std::vector<double> v_expected{0.0, 6.60365e-8};
for (size_t i = 0; i != 2; ++i) {
EXPECT_NEAR(v[i], v_expected[i], 1e-13);
......@@ -935,8 +935,220 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
std::vector<double> b_mag_expected{9.33897e-8, 9.33897e-8, 9.33897e-8, 7.14774e-8, 7.14774e-8, 9.33897e-8};
std::vector<double> b_ang_expected{0.0, -90, 0.0, 112.5, 157.5, -90.0};
for(size_t i = 0; i!= 6;++i) {
for (size_t i = 0; i != 6; ++i) {
EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-13);
EXPECT_NEAR(Bang[i], b_ang_expected[i], FLT_EPSILON);
}
}
class TwoRegionHexagonRotation : public ::testing::Test {
public:
virtual void SetUp() {
// Inner Region
node.push_back(XY{0.0, 0.0});
double r{1.0};
for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
tri.push_back(Triangle<1>(1, 0_zu, 2_zu, 3_zu));
tri.push_back(Triangle<1>(2, 0_zu, 3_zu, 4_zu));
tri.push_back(Triangle<1>(3, 0_zu, 4_zu, 5_zu));
tri.push_back(Triangle<1>(4, 0_zu, 5_zu, 6_zu));
tri.push_back(Triangle<1>(5, 0_zu, 6_zu, 1_zu));
region_i = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3, 4, 5});
boundary_io = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6});
reg.push_back(region_i);
bnd.push_back(boundary_io);
// Outer Region
r = 1.0;
for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
r = 2.0;
for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0 + M_PI / 6.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
tri.push_back(Triangle<1>(6, 7_zu, 13_zu, 8_zu));
tri.push_back(Triangle<1>(7, 8_zu, 13_zu, 14_zu));
tri.push_back(Triangle<1>(8, 8_zu, 14_zu, 9_zu));
tri.push_back(Triangle<1>(9, 9_zu, 14_zu, 15_zu));
tri.push_back(Triangle<1>(10, 9_zu, 15_zu, 10_zu));
tri.push_back(Triangle<1>(11, 10_zu, 15_zu, 16_zu));
tri.push_back(Triangle<1>(12, 10_zu, 16_zu, 11_zu));
tri.push_back(Triangle<1>(13, 11_zu, 16_zu, 17_zu));
tri.push_back(Triangle<1>(14, 11_zu, 17_zu, 12_zu));
tri.push_back(Triangle<1>(15, 12_zu, 17_zu, 18_zu));
tri.push_back(Triangle<1>(16, 12_zu, 18_zu, 7_zu));
tri.push_back(Triangle<1>(17, 7_zu, 18_zu, 13_zu));
region_o = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
boundary_oi = std::make_shared<Boundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12});
boundary_oo = std::make_shared<Boundary<2>>(std::vector<size_t>{13, 14, 15, 16, 17, 18});
reg.push_back(region_o);
bnd.push_back(boundary_oi);
bnd.push_back(boundary_oo);
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
std::shared_ptr<Region<2>> region_i;
std::shared_ptr<Region<2>> region_o;
std::shared_ptr<Boundary<2>> boundary_io;
std::shared_ptr<Boundary<2>> boundary_oi;
std::shared_ptr<Boundary<2>> boundary_oo;
FiniteElementMesh<2, 1> femesh;
MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
void test_flux_density_element_order_1(Eigen::ArrayXd v, Eigen::ArrayXd Bx, Eigen::ArrayXd By, Eigen::ArrayXd Bmag, Eigen::ArrayXd Bang) {
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 = mu_0 * 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);
if (std::abs(180.0 + Bang(i)) < 0.1) {
Bang(i) = -Bang(i);
}
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);
if (std::abs(180.0 + Bang(i)) < 0.1) {
Bang(i) = -Bang(i);
}
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);
if (std::abs(180.0 + Bang(i)) < 0.1) {
Bang(i) = -Bang(i);
}
EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
}
}
};
TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
// TODO: Test with spatially varying forcing condition instead
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, {region_i});
ms.add_magnetic_insulation({boundary_oo});
auto position = ms.add_sliding_interface(boundary_io->nodes(), boundary_oi->nodes(), std::vector<double>(6, 1.0), false);
for (size_t iter = 0; iter != 6; ++iter) {
ms.build_matrices();
ms.apply_conditions();
auto J = ms.init_unknown_matrix();
auto v = ms.init_unknown_vector();
auto r = ms.init_unknown_vector();
auto f = ms.init_unknown_vector();
auto Fx = ms.init_element_array();
auto Fy = ms.init_element_array();
auto dFxdx = ms.init_element_array();
auto dFydy = ms.init_element_array();
auto dFxdy = ms.init_element_array();
v.setZero();
ms.calculate_forcing(f);
double nodal_current{sqrt(3.0) / 4.0 / 3.0};
EXPECT_NEAR(f(0), nodal_current * 6.0, FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
}
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
ldlt.compute(J);
ASSERT_EQ(ldlt.info(), Eigen::Success);
v -= ldlt.solve(r);
// Verify equation is solved
ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
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 * mu_0, mu_0 * FLT_EPSILON);
for (size_t i = 1; i != 7; ++i) {
EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
}
// Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
Eigen::ArrayXd Bmag(By.size());
Eigen::ArrayXd Bang(By.size());
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;
}
test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
++*position;
}
}
\ 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