Commit ee1ea055 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit

Work on HomogeneousForcing for Poisson-like equations
Started implementation of Physics template interface
parent e8ee275c
......@@ -6,21 +6,37 @@
#include "Eigen"
template<size_t Dimension>
#include "FiniteElementMesh.h"
#include "Physics.h"
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class HomogeneousForcing {
};
template<>
class HomogeneousForcing<2> {
template<size_t ElementOrder, size_t QuadratureOrder>
class HomogeneousForcing<2, ElementOrder, QuadratureOrder> {
public:
HomogeneousForcing(std::function<double(double)> ff) : F{ff} {};
HomogeneousForcing(Physics<2, ElementOrder, QuadratureOrder> &op, std::vector<size_t> r, std::function<double(double)> f) : Operator{op}, Regions{r}, Function{f} {};
Eigen::VectorXd operator()(double const t) {
// TODO: Regions only
Eigen::VectorXd output = Operator.basis()[0] * Operator.weights()(0);
double operator()(double const t) { return F(t); };
for (size_t i = 1; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
output += Operator.basis()[i] * Operator.weights()(i);
}
output *= Function(t);
return output;
};
protected:
std::function<double(double)> F;
Physics<2, ElementOrder, QuadratureOrder> &Operator;
std::vector<size_t> Regions;
//std::vector<size_t> Regions;
std::function<double(double)> Function;
};
#endif //OERSTED_CONDITION_H
......@@ -6,6 +6,7 @@
#include "MatrixGroup.h"
#include "Node.h"
#include "Triangle.h"
#include "Region.h"
template<size_t D, size_t P>
class FiniteElementMesh {};
......@@ -15,7 +16,7 @@ class FiniteElementMesh<2,P> {
public:
FiniteElementMesh() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<P>> tris) : Nodes(nodes), Triangles(tris) {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<P>> tris, std::vector<Region<2>> regs) : Nodes(nodes), Triangles(tris), Regions(regs) {};
std::vector<XY> const nodes() const { return Nodes; };
......@@ -57,6 +58,7 @@ public:
protected:
std::vector<XY> Nodes;
std::vector<Triangle<P>> Triangles;
std::vector<Region<2>> Regions;
};
#endif //OERSTED_FINITEELEMENTMESH_H
......@@ -37,7 +37,8 @@ public:
}
}
auto &operator()(size_t q) { return Matrices[q]; };
auto &operator()(size_t const q) const { return Matrices[q]; };
auto &operator()(size_t const q) { return Matrices[q]; };
protected:
std::array<Eigen::DiagonalMatrix<double, Eigen::Dynamic>, Q> Matrices;
......
......@@ -27,17 +27,35 @@ enum class Variable {
Psi = 26 + 23
};
template<size_t Dimension, size_t ElementOrder, Variable V>
class Magnetostatic {
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class Physics {
public:
static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
/*
virtual Eigen::VectorXd init_residual() = 0;
virtual Eigen::SparseMatrix<double> init_jacobian() = 0;
virtual void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) = 0;
virtual void build_matrices() = 0;
virtual void apply_conditions() = 0;
*/
virtual DiagonalMatrixGroup<QuadraturePoints> const & weights() = 0;
virtual SparseMatrixGroup<QuadraturePoints> const & basis() = 0;
virtual DerivativeMatrixGroup<QuadraturePoints> const & derivatives() = 0;
};
template<size_t ElementOrder>
class Magnetostatic<2, ElementOrder, Variable::MagneticVectorPotential> {
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, Variable V>
class Magnetostatic : public Physics<Dimension, ElementOrder, QuadratureOrder> {
};
template<size_t ElementOrder, size_t QuadratureOrder>
class Magnetostatic<2, ElementOrder, QuadratureOrder, Variable::MagneticVectorPotential> : public Physics<2, ElementOrder, QuadratureOrder> {
public:
static constexpr size_t QuadratureOrder = 2 * ElementOrder - 1;
static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
Magnetostatic(FiniteElementMesh<2, ElementOrder> d) : Domain{d},
Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
ElementWeights{d.triangles().size()},
Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
......@@ -68,13 +86,16 @@ public:
void apply_conditions() {};
DiagonalMatrixGroup<QuadraturePoints> const & weights() override {return ElementWeights;};
SparseMatrixGroup<QuadraturePoints> const & basis() override {return Basis;};
DerivativeMatrixGroup<QuadraturePoints> const & derivatives() override {return Derivatives;};
public:
FiniteElementMesh<2, ElementOrder> Domain;
FiniteElementMesh<2, ElementOrder> &Domain;
DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
SparseMatrixGroup<QuadraturePoints> Basis;
DerivativeMatrixGroup<QuadraturePoints> Derivatives;
};
#endif //OERSTED_PHYSICS_H
......@@ -7,17 +7,20 @@ inline constexpr size_t operator "" _zu(unsigned long long int i) {
class MasterTriangle : public ::testing::Test {
public:
virtual void SetUp() {
node.push_back(XY{0.0,0.0});
node.push_back(XY{1.0,0.0});
node.push_back(XY{0.0,1.0});
node.push_back(XY{0.0, 0.0});
node.push_back(XY{1.0, 0.0});
node.push_back(XY{0.0, 1.0});
tri.push_back(Triangle<1>(0, 1_zu, 2_zu, 0_zu));
femesh = FiniteElementMesh<2,1>(node,tri);
reg.push_back(Region<2>(std::vector<size_t>{0}));
femesh = FiniteElementMesh<2,1>(node, tri, reg);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
FiniteElementMesh<2,1> femesh;
};
......@@ -58,25 +61,27 @@ TEST_F(MasterTriangle, stiffness_matrix) {
EXPECT_NEAR(K.coeffRef(2, 2), 1.0 / 2.0, FLT_EPSILON);
}
TEST_F(MasterTriangle, magnetostatic) {
Magnetostatic<2,1,Variable::A> ms(femesh);
static constexpr size_t Q = 1;
Magnetostatic<2, 1, Q, Variable::A> ms(femesh);
Magnetostatic<2, 1, Q, Variable::A> *ms_ptr{&ms};
Physics<2,1,Q> *ph_ptr = dynamic_cast<Physics<2,1,Q>*>(ms_ptr);
ms.build_matrices();
auto J = ms.init_jacobian();
auto r = ms.init_residual();
ms.residual(J,r);
std::cout << J << std::endl;
std::cout << r << std::endl;
ms.residual(J, r);
HomogeneousForcing<2> hf([](double t){return 2.0 * t;});
auto ff = [](double t) { return 1.0 * t; };
std::cout << hf(M_PI) << std::endl;
HomogeneousForcing<2,1,Q> hf{ms, std::vector<size_t>{0}, ff};
std::cout << sizeof(size_t) << std::endl;
std::cout << hf(1.0) << std::endl;
}
class MasterTriangleRotationReflection : public ::testing::Test {
......@@ -96,11 +101,14 @@ 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
femesh = FiniteElementMesh<2, 1>(node, tri);
reg.push_back(Region<2>(std::vector<size_t>{0,1,2,3})); // Triangles 4,5,6 are duplicates
femesh = FiniteElementMesh<2, 1>(node, tri, reg);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
FiniteElementMesh<2, 1> femesh;
};
......@@ -219,7 +227,10 @@ 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));
femesh = FiniteElementMesh<2, 1>(node, tri);
reg.push_back(Region<2>(std::vector<size_t>{0}));
reg.push_back(Region<2>(std::vector<size_t>{1}));
femesh = FiniteElementMesh<2, 1>(node, tri, reg);
vx_m << 1.0, 0.0, 1.0, 0.0;
vx_p << 0.0, 1.0, 0.0, 1.0;
......@@ -229,6 +240,8 @@ public:
std::vector<XY> node;
std::vector<Triangle<1>> tri;
std::vector<Region<2>> reg;
Eigen::Vector4d vx_m;
Eigen::Vector4d vx_p;
......
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