Skip to content

GitLab

Commit d093a288 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit 

More infrastructure work for periodic b.c.s
parent 2c213373
Hide whitespace changes
Inline Side-by-side
src/Physics/src/Boundary.h
View file @ d093a288
......@@ -10,12 +10,24 @@ class Boundary<2> {
public:
Boundary(std::vector<size_t> nodes) : Nodes{nodes} {};
size_t size() const { return Nodes.size(); };
std::vector<size_t> &nodes() { return Nodes; };
std::vector<size_t> const &nodes() const { return Nodes; };
size_t node(size_t i) { return Nodes[i]; };
size_t const &node(size_t i) const { return Nodes[i]; };
protected:
std::vector<size_t> Nodes;
};
template <size_t Dimension>
class BoundaryPair {
};
#endif //OERSTED_BOUNDARY_H
src/Physics/src/BoundaryCondition.h
View file @ d093a288
......@@ -28,18 +28,20 @@ class ZeroDirichlet : public BoundaryCondition{};
template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
public:
ZeroDirichlet(std::vector<size_t> nodes) : Nodes{nodes} {};
ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
void reduce(std::set<size_t,std::greater<size_t>> &index) const override {
for (size_t i = 0; i != Nodes.size(); ++i) {
index.insert(Nodes[i]);
for (auto const& b : Boundaries) {
for (auto const& i : b->nodes()) {
index.insert(i);
}
}
};
protected:
std::vector<size_t> Nodes;
std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
};
......@@ -49,7 +51,7 @@ class PeriodicBoundaryCondition : public BoundaryCondition{};
template<size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
public:
PeriodicBoundaryCondition(std::vector<size_t> first_nodes, std::vector<size_t> second_nodes, bool antiperiodic) : FirstNodes{first_nodes}, SecondNodes{second_nodes}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<std::shared_ptr<BoundaryPair<2>>> boundary_pairs, bool antiperiodic) : BoundaryPairs{boundary_pairs}, Antiperiodic{antiperiodic} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double sign{1.0 - 2.0 * Antiperiodic};
......@@ -64,14 +66,17 @@ public:
}
void reduce(std::set<size_t, std::greater<size_t>> &index) const override {
for (size_t i = 0; i != SecondNodes.size(); ++i) {
index.insert(SecondNodes[i]);
/*
for (auto const &bp : BoundaryPairs) {
for(auto const &ij : bp->nodes()) {
index.insert(ij.second);
}
}
*/
}
protected:
std::vector<size_t> FirstNodes;
std::vector<size_t> SecondNodes;
std::vector<std::shared_ptr<BoundaryPair<2>>> BoundaryPairs;
bool Antiperiodic;
};
......
src/Physics/src/FiniteElementMesh.h
View file @ d093a288
......@@ -17,29 +17,37 @@ class FiniteElementMesh<2,Order> {
public:
FiniteElementMesh() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<Region<2>> r, std::vector<Boundary<2>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris,
std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<Boundary<2>>> b)
: Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
size_t size_nodes() const { return Nodes.size(); };
size_t size_elements() const { return Triangles.size(); };
std::vector<XY> const &nodes() const { return Nodes; };
XY const &node(size_t i) const { return Nodes[i]; };
size_t size_elements() const { return Triangles.size(); };
std::vector<Triangle<Order>> const &triangles() const { return Triangles; };
std::vector<Region<2>> &regions() { return Regions; };
std::vector<Region<2>> const &regions() const { return Regions; };
Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
std::vector<Boundary<2>> const &boundaries() const { return Boundaries; };
auto &regions() { return Regions; };
XY const &node(size_t i) const { return Nodes[i]; };
auto &region(size_t i) { return Regions[i]; };
Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
auto const &regions() const { return Regions; };
auto const &region(size_t i) const { return Regions[i]; };
auto &boundaries() {return Boundaries;};
auto &boundary(size_t i) {return Boundaries[i];};
Region<2> &region(size_t i) { return Regions[i]; };
Region<2> const &region(size_t i) const { return Regions[i]; };
auto const &boundaries() const { return Boundaries; };
Boundary<2> const &boundary(size_t i) const { return Boundaries[i]; };
auto const &boundary(size_t i) const { return Boundaries[i]; };
template<size_t Q>
DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> determinant() const {
......@@ -78,8 +86,8 @@ protected:
std::vector<XY> Nodes;
std::vector<Triangle<Order>> Triangles;
std::vector<Region<2>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std::vector<Boundary<2>> Boundaries;
std::vector<std::shared_ptr<Region<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
};
#endif //OERSTED_FINITEELEMENTMESH_H
src/Physics/src/Forcing.h
View file @ d093a288
......@@ -25,16 +25,18 @@ template<size_t ElementOrder, size_t QuadratureOrder>
class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Forcing {
public:
HomogeneousForcing(std::function<double(double)> f,
std::vector<size_t> elements,
std::vector<std::shared_ptr<Region<2>>> regions,
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &basis,
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &weights)
: Function{f}, Elements{elements}, Basis(basis), Weights(weights) {
: Function{f}, Regions{regions}, Basis(basis), Weights(weights) {
// TODO: May need to adjust implementation once quadrilaterals are added
Indicator.resize(Weights(0).size());
Indicator.setZero();
for (size_t const &i : Elements) {
Indicator(i) = 1.0;
for (auto const &r : Regions) {
for(auto const &i : r->elements()) {
Indicator(i) = 1.0;
}
}
};
......@@ -53,7 +55,7 @@ public:
protected:
std::function<double(double)> Function;
std::vector<size_t> Elements;
std::vector<std::shared_ptr<Region<2>>> Regions;
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Basis;
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Weights;
......
src/Physics/src/Physics.h
View file @ d093a288
......@@ -121,7 +121,7 @@ public:
// Calculate polarization
for(auto &dr : Domain.regions()) {
dr.material().h_from_b(dr.triangles(),Fx,Fy,dFxdx,dFydy,dFxdy);
dr->material().h_from_b(dr->elements(),Fx,Fy,dFxdx,dFydy,dFxdy);
}
// Calculate residual
......@@ -154,16 +154,16 @@ public:
}
}
void add_current_density(std::function<double(double)> func, std::vector<size_t> elements) {
ForcingCondtions.push_back(std::make_unique<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, elements, Basis, ElementWeights));
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));
}
void remove_current_density(size_t i) {
ForcingCondtions.erase(ForcingCondtions.begin() + i);
}
void add_magnetic_insulation(std::vector<size_t> nodes) {
BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(nodes));
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));
}
void apply_conditions() {
......@@ -183,16 +183,28 @@ public:
BoundaryConditions[i]->reduce(index);
}
// TODO: Rewrite following section:
// (Need to handle arbitrary number of entries in the b.c. "permutation" matrix)
// First, search for triplets with row() == reduce and erase()
// Then, search for skipped rows and cumulatively renumber
// BEGIN
size_t rows{Domain.nodes().size()};
size_t cols{Domain.nodes().size()};
for (auto &i : index) {
triplets.erase(triplets.begin()+i);
for (auto &i : index) { // index is sorted from greatest to least, so triplets are erased in reverse order
auto iter = (triplets.begin() + i);
if(iter->row() == i) {
triplets.erase(iter);
}
--rows;
}
for (size_t i = 0; i != triplets.size(); ++i) {
triplets[i] = Eigen::Triplet<double>(i, triplets[i].col(), triplets[i].value());
}
// END
Eigen::SparseMatrix<double> bc_matrix(rows,cols);
bc_matrix.setFromTriplets(triplets.begin(),triplets.end());
......
src/Physics/src/Region.h
View file @ d093a288
......@@ -12,19 +12,19 @@ class Region {
template<>
class Region<2> { // TODO: Rename Triangles to Elements?
public:
Region(std::vector<size_t> tris) : Triangles{tris}, Material{Air()} {};
Region(std::vector<size_t> tris, MaterialProperties mat) : Triangles{tris}, Material{mat} {};
Region(std::vector<size_t> elements) : Elements{elements}, Material{Air()} {};
Region(std::vector<size_t> elements, MaterialProperties material) : Elements{elements}, Material{material} {};
std::vector<size_t> const &triangles() const { return Triangles; };
std::vector<size_t> const &elements() const { return Elements; };
size_t const &triangle(size_t i) const { return Triangles[i]; };
size_t const &element(size_t i) const { return Elements[i]; };
MaterialProperties &material() { return Material; }; // non-const because internal state of the material may be altered
void material(MaterialProperties mat) { Material = mat; };
protected:
std::vector<size_t> Triangles;
std::vector<size_t> Elements;
MaterialProperties Material;
};
......
test/Physics/test_FiniteElementMesh.cpp
View file @ d093a288
......@@ -13,19 +13,20 @@ public:
tri.push_back(Triangle<1>(0, 1_zu, 2_zu, 0_zu));
reg.push_back(Region<2>(std::vector<size_t>{0}));
reg.push_back(std::make_shared<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(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 1}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 2}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2}));
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
FiniteElementMesh<2, 1> femesh;
};
......@@ -76,7 +77,7 @@ TEST_F(MasterTriangle, magnetostatic) {
auto ff = [](double t) { return 1.0 * t; };
ms.add_current_density(ff, std::vector<size_t>{0});
ms.add_current_density(ff, reg);
ms.build_matrices();
......@@ -117,16 +118,16 @@ 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(std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
FiniteElementMesh<2, 1> femesh;
};
......@@ -245,14 +246,14 @@ public:
tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
tri.push_back(Triangle<1>(1, 1_zu, 3_zu, 2_zu));
reg.push_back(Region<2>(std::vector<size_t>{0}));
reg.push_back(Region<2>(std::vector<size_t>{1}));
reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0}));
reg.push_back(std::make_shared<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(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 1}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 2}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 3}));
bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{2, 3}));
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
......@@ -264,8 +265,9 @@ public:
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
Eigen::Vector4d vx_m;
Eigen::Vector4d vx_p;
......@@ -368,7 +370,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
auto ff = [](double t) { return 1.0 * t; };
auto f = ms.init_unknown_vector();
ms.add_current_density(ff, std::vector<size_t>{0});
ms.add_current_density(ff, {reg[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);
......@@ -376,7 +378,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
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.add_current_density(ff, {reg[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);
......@@ -384,7 +386,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
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, {reg[0], reg[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);
......@@ -463,7 +465,6 @@ public:
r = 2.0;
for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0 + M_PI / 6.0;
outer_boundary.push_back(node.size());
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
......@@ -474,11 +475,11 @@ 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));
inner_region = {0,1,2,3,4,5};
inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3, 4, 5});
inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6});
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}));
reg.push_back(inner_region);
bnd.push_back(inner_boundary);
tri.push_back(Triangle<1>(6, 1_zu, 7_zu, 2_zu));
tri.push_back(Triangle<1>(7, 2_zu, 7_zu, 8_zu));
......@@ -498,24 +499,25 @@ 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));
outer_region = {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}));
outer_region = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
outer_boundary = std::make_shared<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}));
reg.push_back(outer_region);
bnd.push_back(outer_boundary);
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
std::vector<std::shared_ptr<Region<2>>> reg;
std::vector<std::shared_ptr<Boundary<2>>> bnd;
std::vector<size_t> inner_region;
std::vector<size_t> outer_region;
std::shared_ptr<Region<2>> inner_region;
std::shared_ptr<Region<2>> outer_region;
std::vector<size_t> outer_boundary;
std::shared_ptr<Boundary<2>> inner_boundary;
std::shared_ptr<Boundary<2>> outer_boundary;
FiniteElementMesh<2, 1> femesh;
......@@ -576,8 +578,8 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
auto ff = [](double t) { return 1.0; };
auto f = ms.init_unknown_vector();
ms.add_current_density(ff, inner_region);
ms.add_magnetic_insulation(outer_boundary);
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
ms.calculate_forcing(f);
double nodal_current{sqrt(3.0) / 4.0 / 3.0};
......@@ -655,7 +657,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
}
TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
femesh.region(1).material(Linear1000);
femesh.region(1)->material(Linear1000);
Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
ms.time(0.0);
......@@ -663,8 +665,8 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, inner_region);
ms.add_magnetic_insulation(outer_boundary);
ms.add_current_density(ff, {inner_region});
ms.add_magnetic_insulation({outer_boundary});
ms.apply_conditions();
......
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