Commit 2aee9a7b authored by JasonPries's avatar JasonPries
Browse files

Modify PeriodicBoundaryCondition to account for duplicate mapped nodes and mapping of the origin

parent 08f11fbf
......@@ -57,7 +57,7 @@ protected:
class VariableMap {
public:
VariableMap(size_t first, size_t second, double value) : First{first}, Second{second}, Value{value} {};
VariableMap(size_t first, size_t second, double value) : First{std::min(first,second)}, Second{std::max(first,second)}, Value{value} {};
size_t first() const { return First; };
......@@ -65,6 +65,14 @@ public:
double value() const { return Value; };
inline bool operator==(VariableMap const &rhs) {
return (First == rhs.First && Second == rhs.Second && Value == rhs.Value);
}
inline bool operator<(VariableMap const &rhs) {
return (First < rhs.First) || (First == rhs.First && Second < rhs.Second);
}
protected:
size_t First;
size_t Second;
......
......@@ -25,13 +25,21 @@ 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 {
public:
ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};
ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<2, ElementOrder> const &fem) {
for (auto const &b : boundaries) {
auto fb = fem.boundary(b);
if (fb) { Boundaries.push_back(fb); }
}
}
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
......@@ -49,73 +57,79 @@ protected:
// 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 {
public:
PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<VariableMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2,ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
Map.reserve(periodic_boundary_pairs.size());
PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2, ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
double_t sgn = sign();
double_t sign = (antiperiodic ? -1.0 : 1.0);
for(auto const &pbp : periodic_boundary_pairs) {
// TODO: Change Map type to std::set<VariableMap>, add comparison operator to VariableMap, and make sure node indicies are lexigraphically ordered in VariableMap
for (auto const &pbp : periodic_boundary_pairs) {
auto const &n0 = fem.boundary(pbp.curve0())->nodes();
auto const &n1 = fem.boundary(pbp.curve1())->nodes();
std::vector<VariableMap> vmap;
vmap.reserve(n0.size());
if (pbp.orientation()) {
for (size_t i = 0; i != n0.size(); ++i) {
vmap.emplace_back(n0[i], n1[i], sign);
Map.emplace_back(n0[i], n1[i], sgn);
}
} else {
size_t j = n1.size();
for (size_t i = 0; i != n0.size(); ++i) {
vmap.emplace_back(n0[i], n1[--j], sign);
Map.emplace_back(n0[i], n1[--j], sgn);
}
}
}
Map.emplace_back(vmap);// TODO: needs testing, allows duplicate nodes which may not be handle correctly by constraint algorithm
// Remove duplicates
std::sort(Map.begin(), Map.end());
auto last = std::unique(Map.begin(), Map.end());
Map.erase(last, Map.end());
if (!antiperiodic) { // Remove center node if it exists
auto eq = [](VariableMap const &x) { return (x.first() == x.second()); };
auto center = std::find_if(Map.begin(), Map.end(), eq);
if (center != Map.end()) { Map.erase(center); };
}
// TODO: needs testing
}
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double sign{1.0 - 2.0 * Antiperiodic};
double_t sgn = sign();
for (BoundaryMap const &b : Map) {
for (VariableMap const &v : b.map()) {
triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sign * v.value()));
}
for (VariableMap const &v : Map) {
triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sgn * v.value()));
}
}
void reduce(std::set<size_t, std::less<size_t>> &index) const override {
for (BoundaryMap const &b : Map) {
for (VariableMap const &v : b.map()) {
for (VariableMap const &v : Map) {
index.insert(v.second());
}
}
}
double_t sign() const { return (Antiperiodic ? -1.0 : 1.0); };
protected:
std::vector<BoundaryMap> Map;
std::vector<VariableMap> Map;
bool Antiperiodic;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface : public BoundaryCondition {};
class SlidingInterface : public BoundaryCondition {
};
template<size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
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++() {
SlidingInterface &operator++() {
std::rotate(Second.begin(), Second.begin() + 1, Second.end());
if (Antiperiodic) {
std::rotate(Value.begin(), Value.begin() + 1, Value.end());
......@@ -126,13 +140,13 @@ public:
// TODO: SlidingInterface& operator--() {...};
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]));
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){
for (size_t i = 0; i != Second.size(); ++i) {
index.insert(Second[i]);
}
}
......@@ -144,4 +158,5 @@ protected:
bool Antiperiodic;
};
#endif //OERSTED_BOUNDARYCONDITION_H
......@@ -169,7 +169,11 @@ public:
BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
}
void add_periodic_boundary(std::vector<BoundaryMap> map, bool antiperiodic) {
void add_magnetic_insulation(std::vector<std::shared_ptr<Curve const>> boundaries) {
BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries, Domain));
}
void add_periodic_boundary(std::vector<VariableMap> map, bool antiperiodic) {
BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
}
......
......@@ -10,7 +10,7 @@ class Region {
};
template<>
class Region<2> { // TODO: Rename Triangles to Elements?
class Region<2> {
public:
Region() : Material{Air()} {};
Region(std::vector<size_t> elements) : Elements{elements}, Material{Air()} {};
......
......@@ -275,15 +275,64 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
// Create Physics
Magnetostatic<2, 1, 1, FieldVariable::A> msph{fem};
//msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(0)});
// 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>
//msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(1)});
//msph.add_magnetic_insulation({fem.boundary(0), fem.boundary(1), fem.boundary(2)});
msph.add_magnetic_insulation(radial_boundary);
msph.add_periodic_boundary(periodic_boundary, false); GTEST_FATAL_FAILURE_("Needs testing, allows duplicate nodes which may not be handle correctly by constraint algorithm");
msph.add_periodic_boundary(periodic_boundary, false);
//msph.build_matrices();
msph.build_matrices();
msph.apply_conditions();
auto J = msph.init_unknown_matrix();
auto v = msph.init_unknown_vector();
auto r = msph.init_unknown_vector();
auto f = msph.init_unknown_vector();
auto Fx = msph.init_element_array();
auto Fy = msph.init_element_array();
auto dFxdx = msph.init_element_array();
auto dFydy = msph.init_element_array();
auto dFxdy = msph.init_element_array();
v.setZero();
msph.calculate_forcing(f);
msph.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
msph.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);
}
// Test flux-density values
Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -msph.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;
}
std::cout << Bmag << std::endl << std::endl;
std::cout << Bang << std::endl << std::endl;
std::cout << f << std::endl << std::endl;
//msph.apply_conditions();
std::cout << v << std::endl << std::endl;
std::cout << std::endl;
GTEST_NONFATAL_FAILURE_("Solution needs testing");
}
\ No newline at end of file
......@@ -795,7 +795,7 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_periodic) {
ms.add_magnetic_insulation({outer_boundary});
ms.add_periodic_boundary({variable_map}, false);
ms.add_periodic_boundary(variable_map, false);
ms.build_matrices();
......@@ -873,7 +873,7 @@ TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
ms.add_magnetic_insulation({outer_boundary});
variable_map.push_back(VariableMap(0, 0, -1.0));
ms.add_periodic_boundary({variable_map}, true);
ms.add_periodic_boundary(variable_map, true);
ms.build_matrices();
......@@ -1166,3 +1166,5 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
++*position;
}
}
// TODO: Test rotation + periodic boundaries
\ 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