Nightly Commit

Work on HomogeneousForcing for Poisson-like equations
Started implementation of Magnetostatic solver

src/Physics/src/Condition.h

@@ -6,34 +6,35 @@
 
 #include "Eigen"
 
-#include "FiniteElementMesh.h"
+#include "Region.h"
 #include "Physics.h"
 
+class Condition {
+public:
+    virtual ~Condition() {};
+
+    virtual Eigen::VectorXd operator()(double const t) const = 0;
+};
+
 template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
-class HomogeneousForcing {
+class HomogeneousForcing : public Condition {
 };
 
 template<size_t ElementOrder, size_t QuadratureOrder>
-class HomogeneousForcing<2, ElementOrder, QuadratureOrder> {
+class HomogeneousForcing<2, ElementOrder, QuadratureOrder> : public Condition {
 public:
-    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());
+    HomogeneousForcing(Physics<2, ElementOrder, QuadratureOrder> &op, std::vector<size_t> reg, std::function<double(double)> f) : Operator{op}, Regions{reg}, Function{f} {
+        Indicator.resize(op.domain().triangles().size());  // TODO: num_elements (once quadrilateral elements added)
         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;
+        for (size_t const &i : reg) {
+            for (size_t const &j : op.domain().region(i).triangles()) {
+                Indicator(j) = 1.0;
             }
-        }
+    }
     };
 
-    auto const &indicator() const { return Indicator; };
-
-    Eigen::VectorXd operator()(double const t) const {
+    Eigen::VectorXd operator()(double const t) const override {
         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) * Indicator));
@@ -43,14 +44,16 @@ public:
         return output;
     };
 
+    auto const &indicator() const { return Indicator; };
+
 protected:
     Physics<2, ElementOrder, QuadratureOrder> &Operator;
 
     std::vector<size_t> Regions;
 
-    Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
-
     std::function<double(double)> Function;
+
+    Eigen::VectorXd Indicator; // TODO: ?Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator?
 };
 
 #endif //OERSTED_CONDITION_H

src/Physics/src/FiniteElementMesh.h

 #ifndef OERSTED_FINITEELEMENTMESH_H
 #define OERSTED_FINITEELEMENTMESH_H
 
+#include <memory>
 #include <vector>
 
 #include "MatrixGroup.h"
@@ -8,25 +9,25 @@
 #include "Triangle.h"
 #include "Region.h"
 
-template<size_t D, size_t P>
+template<size_t Dimension, size_t Order>
 class FiniteElementMesh {};
 
-template<size_t P>
-class FiniteElementMesh<2,P> {
+template<size_t Order>
+class FiniteElementMesh<2,Order> {
 public:
     FiniteElementMesh() {};
 
-    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<P>> tris, std::vector<Region<2>> regs) : Nodes(nodes), Triangles(tris), Regions(regs) {};
+    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<Region<2>> r) : Nodes(nodes), Triangles(tris), Regions(r) {};
 
     std::vector<XY> const &nodes() const { return Nodes; };
 
-    std::vector<Triangle<P>> const &triangles() const { return Triangles; };
+    std::vector<Triangle<Order>> 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]; };
+    Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
 
     Region<2> const &region(size_t i) const { return Regions[i]; };
 
@@ -43,7 +44,7 @@ public:
 
     template<size_t Q>
     SparseMatrixGroup<TriangleQuadratureRule<Q>::size> basis() const {
-        SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
+        SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
 
         for (size_t i = 0; i != Triangles.size(); ++i) {
             Triangles[i].basis<Q>(mat, Nodes);
@@ -54,7 +55,7 @@ public:
 
     template<size_t Q>
     DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> derivative() const {
-        DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
+        DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
 
         for (size_t i = 0; i!= Triangles.size();++i) {
             Triangles[i].derivative<Q>(df, Nodes);
@@ -65,8 +66,9 @@ public:
 
 protected:
     std::vector<XY> Nodes;
-    std::vector<Triangle<P>> Triangles;
-    std::vector<Region<2>> Regions;
+    std::vector<Triangle<Order>> Triangles;
+
+    std::vector<Region<2>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
 };
 
 #endif //OERSTED_FINITEELEMENTMESH_H

src/Physics/src/Physics.h

@@ -27,13 +27,12 @@ enum class Variable {
     Psi = 26 + 23
 };
 
-
 template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
 class Physics {
 public:
     static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
 
-    /*
+    /* TODO: template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder> class PhysicsData : public PhysicsInterface {...};
     virtual Eigen::VectorXd init_residual() = 0;
     virtual Eigen::SparseMatrix<double> init_jacobian() = 0;
     virtual void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) = 0;
@@ -76,13 +75,6 @@ public:
     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()); };

src/Physics/src/Region.h

@@ -2,18 +2,20 @@
 #define OERSTED_REGION_H
 
 #include <vector>
+#include <cstddef>
 
 template<size_t Dimension>
-class Region {};
+class Region {
+};
 
 template<>
 class Region<2> {
 public:
-    Region(std::vector<size_t> tris) : Triangles(tris) {};
+    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; };
 
-    std::vector<size_t> const& triangles() const { return Triangles; };
+    size_t const &triangle(size_t i) const { return Triangles[i]; };
 
 protected:
     std::vector<size_t> Triangles;