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

Nightly Commit

Simplification and clarification of (discrete) boundary conditions effecting the implementation of FiniteElementMesh classes
parent 91622842
......@@ -11,34 +11,24 @@
#include "QuadratureRule.h"
#include "FiniteElementMesh.h"
// TODO: DiscreteBoundaryCondition versus ContinuousBoundaryCondition? Somehow need to apply to top level model and have that be translated to the discrete model
// TODO: Boundary conditions should represent simple containers that know how to perform actions on matrices, parsing of boundaries to extract nodes should take place at a higher level
class BoundaryCondition {
public:
BoundaryCondition(std::vector<size_t> boundaries) : Boundaries{boundaries} {};
virtual ~BoundaryCondition(){};
virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;
virtual void reduce(std::set<size_t,std::greater<size_t>> &index) const = 0;
protected:
std::vector<size_t> Boundaries;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet : public BoundaryCondition{};
template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
public:
ZeroDirichlet(std::vector<size_t> boundaries,
FiniteElementMesh<2,ElementOrder> const &domain) : BoundaryCondition{boundaries} {
for(size_t i : boundaries) {
for(size_t j : domain.boundary(i).nodes()) {
Nodes.push_back(j);
}
}
};
ZeroDirichlet(std::vector<size_t> nodes) : Nodes{nodes} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};
......@@ -51,4 +41,40 @@ public:
protected:
std::vector<size_t> Nodes;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
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} {};
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double sign{1.0 - 2.0 * Antiperiodic};
// TODO:
/*
for(size_t i = 0;i!=FirstNodes.size();++i){
triplets[FirstNodes[i]].row() = SecondNodes[i];
triplest[FirstNodes[i]].value() = sign;
}
*/
}
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]);
}
}
protected:
std::vector<size_t> FirstNodes;
std::vector<size_t> SecondNodes;
bool Antiperiodic;
};
#endif //OERSTED_BOUNDARYCONDITION_H
......@@ -25,18 +25,16 @@ 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> reg,
std::vector<size_t> elements,
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &basis,
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &weights,
FiniteElementMesh<2,ElementOrder> const &domain) : Regions{reg}, Basis(basis), Weights(weights), Function{f} {
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &weights)
: Function{f}, Elements{elements}, Basis(basis), Weights(weights) {
Indicator.resize(domain.triangles().size()); // TODO: num_elements (once quadrilateral elements added)
// TODO: May need to adjust implementation once quadrilaterals are added
Indicator.resize(Weights(0).size());
Indicator.setZero();
for (size_t const &i : reg) {
for (size_t const &j : domain.region(i).triangles()) {
Indicator(j) = 1.0;
}
for (size_t const &i : Elements) {
Indicator(i) = 1.0;
}
};
......@@ -55,7 +53,7 @@ public:
protected:
std::function<double(double)> Function;
std::vector<size_t> Regions;
std::vector<size_t> Elements;
SparseMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Basis;
DiagonalMatrixGroup<TriangleQuadratureRule<QuadratureOrder>::size> const &Weights;
......
......@@ -154,16 +154,16 @@ public:
}
}
void add_current_density(std::function<double(double)> func, std::vector<size_t> regions) {
ForcingCondtions.push_back(std::make_unique<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights, Domain));
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 remove_current_density(size_t i) {
ForcingCondtions.erase(ForcingCondtions.begin() + i);
}
void add_magnetic_insulation(std::vector<size_t> boundaries) {
BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries, Domain));
void add_magnetic_insulation(std::vector<size_t> nodes) {
BoundaryConditions.push_back(std::make_unique<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(nodes));
}
void apply_conditions() {
......
......@@ -53,9 +53,8 @@ TEST_F(MasterTriangle, stiffness_matrix) {
auto df = femesh.derivative<1>();
auto det = femesh.determinant<1>();
Eigen::SparseMatrix<double> K =
(df.dx(0) * det(0).matrix().asDiagonal() * df.dx(0).transpose() + df.dy(0) * det(0).matrix().asDiagonal() * df.dy(0).transpose()) *
TriangleQuadratureRule<1>::w[0];
Eigen::SparseMatrix<double> K = (df.dx(0) * det(0).matrix().asDiagonal() * df.dx(0).transpose()
+ df.dy(0) * det(0).matrix().asDiagonal() * df.dy(0).transpose()) * TriangleQuadratureRule<1>::w[0];
EXPECT_NEAR(K.coeffRef(0, 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(K.coeffRef(1, 0), -1.0 / 2.0, FLT_EPSILON);
......@@ -464,6 +463,7 @@ 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,6 +474,8 @@ 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};
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}));
......@@ -496,6 +498,8 @@ 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}));
bnd.push_back(Boundary<2>(std::vector<size_t>{7, 8, 9, 10, 11, 12}));
......@@ -508,6 +512,11 @@ public:
std::vector<Region<2>> reg;
std::vector<Boundary<2>> bnd;
std::vector<size_t> inner_region;
std::vector<size_t> outer_region;
std::vector<size_t> outer_boundary;
FiniteElementMesh<2, 1> femesh;
MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
......@@ -567,8 +576,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, std::vector<size_t>{0});
ms.add_magnetic_insulation(std::vector<size_t>{1});
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};
......@@ -654,8 +663,8 @@ TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
auto ff = [](double t) { return 1.0; };
ms.add_current_density(ff, std::vector<size_t>{0});
ms.add_magnetic_insulation(std::vector<size_t>{1});
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