Commit 3f264747 authored by JasonPries's avatar JasonPries
Browse files

Nightly Commit

Work on HomogeneousForcing for Poisson-like equations
Started implementation of Physics template interface
parent ee1ea055
......@@ -16,15 +16,27 @@ class HomogeneousForcing {
template<size_t ElementOrder, size_t QuadratureOrder>
class HomogeneousForcing<2, ElementOrder, QuadratureOrder> {
public:
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
HomogeneousForcing(Physics<2, ElementOrder, QuadratureOrder> &op, std::vector<size_t> r, std::function<double(double)> f) : Operator{op}, Regions{r}, Function{f} {
auto const& dom = op.domain();
auto const& tris = dom.triangles();
Indicator.resize(tris.size());
Indicator.setZero();
for(size_t i = 0;i!=r.size();++i){
auto const& t = dom.region(r[i]).triangles();
for(size_t j = 0;j!=t.size();++j){
Indicator(tris[t[j]].id()) = 1.0;
}
}
};
Eigen::VectorXd output = Operator.basis()[0] * Operator.weights()(0);
auto const &indicator() const { return Indicator; };
Eigen::VectorXd operator()(double const t) const {
Eigen::VectorXd output = (Operator.basis()[0] * (Operator.weights()(0) * Indicator)); // TODO: ?SparseVector?
for (size_t i = 1; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
output += Operator.basis()[i] * Operator.weights()(i);
output += (Operator.basis()[i] * (Operator.weights()(i) * Indicator));
}
output *= Function(t);
......@@ -36,6 +48,8 @@ protected:
std::vector<size_t> Regions;
Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
std::function<double(double)> Function;
};
......
......@@ -18,9 +18,17 @@ public:
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; };
std::vector<XY> const &nodes() const { return Nodes; };
std::vector<Triangle<P>> const triangles() const { return Triangles; };
std::vector<Triangle<P>> const &triangles() const { return Triangles; };
std::vector<Region<2>> const &regions() const { return Regions; };
XY const &node(size_t i) const { return Nodes[i]; };
Triangle<P> const &triangle(size_t i) const { return Triangles[i]; };
Region<2> const &region(size_t i) const { return Regions[i]; };
template<size_t Q>
DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> determinant() const {
......
......@@ -7,29 +7,32 @@
#include "PhysicsCommon.h"
template<Dimension d> class Discretization;
template<Dimension d>
class Discretization;
template<VectorSpace vs, Dimension d>
class Operator {
public:
Eigen::VectorXd operator()(Eigen::VectorXd vec) { return Matrix * vec;};
static const VectorSpace Space = vs;
static const Dimension Dim = d;
Operator(std::shared_ptr <Discretization<d>> dsc) : Domain(dsc) {};
Eigen::VectorXd operator()(Eigen::VectorXd vec) { return Matrix * vec; };
Eigen::MatrixXd operator()(Eigen::MatrixXd mat) { return Matrix * mat; };
static const VectorSpace Space = vs;
static const Dimension Dim = d;
Discretization<d> const &domain() { return Domain->get(); };
protected:
Operator(std::shared_ptr<Discretization<d>> dsc) : Domain(dsc) {};
std::shared_ptr<Discretization<d>> Domain;
std::shared_ptr <Discretization<d>> Domain;
Eigen::MatrixXd Matrix;
};
template<VectorSpace vs, Dimension d> class Curl : public Operator<vs,d> {
template<VectorSpace vs, Dimension d>
class Curl : public Operator<vs, d> {
public:
Curl(std::shared_ptr<Discretization<d>> dsc) : Operator<vs,d>(dsc) {};
Curl(std::shared_ptr <Discretization<d>> dsc) : Operator<vs, d>(dsc) {};
};
#endif //OERSTED_OPERATOR_H
......@@ -41,9 +41,25 @@ public:
virtual void apply_conditions() = 0;
*/
virtual DiagonalMatrixGroup<QuadraturePoints> const & weights() = 0;
virtual SparseMatrixGroup<QuadraturePoints> const & basis() = 0;
virtual DerivativeMatrixGroup<QuadraturePoints> const & derivatives() = 0;
Physics(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} {};
FiniteElementMesh<Dimension, ElementOrder> &domain() const { return Domain; };
DiagonalMatrixGroup<QuadraturePoints> const &weights() { return ElementWeights; };
SparseMatrixGroup<QuadraturePoints> const &basis() { return Basis; };
DerivativeMatrixGroup<QuadraturePoints> const &derivatives() { return Derivatives; };
protected:
FiniteElementMesh<Dimension, ElementOrder> &Domain;
DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
SparseMatrixGroup<QuadraturePoints> Basis;
DerivativeMatrixGroup<QuadraturePoints> Derivatives;
};
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, Variable V>
......@@ -53,12 +69,21 @@ 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 QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
using Physics<2, ElementOrder, QuadratureOrder>::QuadraturePoints;
using Physics<2, ElementOrder, QuadratureOrder>::domain;
using Physics<2, ElementOrder, QuadratureOrder>::derivatives;
using Physics<2, ElementOrder, QuadratureOrder>::weights;
using Physics<2, ElementOrder, QuadratureOrder>::basis;
/*
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} {};
*/
Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : Physics<2, ElementOrder, QuadratureOrder>{d} {};
Eigen::VectorXd init_residual() { return Eigen::VectorXd(Domain.nodes().size()); };
......@@ -86,16 +111,11 @@ 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;
DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
SparseMatrixGroup<QuadraturePoints> Basis;
DerivativeMatrixGroup<QuadraturePoints> Derivatives;
protected:
using Physics<2, ElementOrder, QuadratureOrder>::Domain;
using Physics<2, ElementOrder, QuadratureOrder>::Derivatives;
using Physics<2, ElementOrder, QuadratureOrder>::ElementWeights;
using Physics<2, ElementOrder, QuadratureOrder>::Basis;
};
#endif //OERSTED_PHYSICS_H
#ifndef OERSTED_REGION_H
#define OERSTED_REGION_H
#include <vector>
template<size_t Dimension>
class Region {};
......@@ -9,6 +11,10 @@ class Region<2> {
public:
Region(std::vector<size_t> tris) : Triangles(tris) {};
size_t const &triangle(size_t i) const { return Triangles[i]; };
std::vector<size_t> const& triangles() const { return Triangles; };
protected:
std::vector<size_t> Triangles;
};
......
#ifndef OERSTED_TRIANGLE_H
#define OERSTED_TRIANGLE_H
#include <vector>
#include "Eigen"
#include "Eigen/Sparse"
......@@ -16,6 +18,8 @@ public:
Element() : ID{0} {};
Element(size_t id) : ID{id} {};
size_t const id() const { return ID; };
protected:
size_t ID;
};
......
......@@ -62,6 +62,7 @@ TEST_F(MasterTriangle, stiffness_matrix) {
}
TEST_F(MasterTriangle, magnetostatic) {
//TODO: This isn't really testing anything...
static constexpr size_t Q = 1;
Magnetostatic<2, 1, Q, Variable::A> ms(femesh);
......@@ -327,4 +328,43 @@ TEST_F(SimpleSquare, derivative_2_temp) {
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
}
}
TEST_F(SimpleSquare, magnetostatic_forcing) {
// TODO: Write Physics::add_condition() (templated?)
// TODO: BoundaryCondition, ForcingCondition, MotionCondition?
// TODO: ?then template parameters of conditions can be hidden?
// TODO: Hierarchy of conditions: Those altering RHS only (forcing, some bcs), those altering matrices in a fixed way (some bcs), those altering the mesh (motion, some bcs)
// TODO: Applied in reverse order, Motion->BoundaryCondition->Forcing,
// TODO: Or, TimeVarying->TimeInvariant->Forcing?
Magnetostatic<2, 1, 2, Variable::A> ms{femesh};
ms.build_matrices();
auto ff = [](double t) { return 1.0 * t; };
HomogeneousForcing<2,1,2> hf0{ms, std::vector<size_t>{0}, ff};
//ms.add_current_density(std::vector<size_t{0}, ff};
HomogeneousForcing<2,1,2> hf1{ms, std::vector<size_t>{1}, ff};
HomogeneousForcing<2,1,2> hf01{ms, std::vector<size_t>{0,1}, ff};
auto v = hf0(1.0);
EXPECT_NEAR(v(0),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(1),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(2),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(3),0.0/6.0,FLT_EPSILON);
v = hf1(1.0);
EXPECT_NEAR(v(0),0.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(1),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(2),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(3),1.0/6.0,FLT_EPSILON);
v = hf01(1.0);
EXPECT_NEAR(v(0),1.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(1),2.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(2),2.0/6.0,FLT_EPSILON);
EXPECT_NEAR(v(3),1.0/6.0,FLT_EPSILON);
}
\ 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