Start implementation of magnetostatic solver

src/Physics/CMakeLists.txt

@@ -5,19 +5,23 @@ set(SOURCE_FILES
 
         #./src/PhysicsCommon.h
 
-        #./src/Physics.h ./src/Physics.cpp
+        ./src/Condition.h ./src/Condition.cpp
 
-        #./src/Operator.h ./src/Operator.cpp
+        ./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
 
-        #./src/Discretization.h ./src/Discretization.cpp
+        ./src/MatrixGroup.h ./src/MatrixGroup.cpp
 
-        ./src/QuadratureRule.h ./src/QuadratureRule.cpp
+        ./src/Node.h ./src/Node.cpp
 
-        ./src/Triangle.h ./src/Triangle.cpp
+        ./src/Physics.h ./src/Physics.cpp
 
-        ./src/Node.h ./src/Node.cpp
+        #./src/Operator.h ./src/Operator.cpp
 
-        ./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
+        ./src/QuadratureRule.h ./src/QuadratureRule.cpp
+
+        ./src/Region.h ./src/Region.cpp
+
+        ./src/Triangle.h ./src/Triangle.cpp
         )
 
 add_library(physics SHARED ${SOURCE_FILES})

src/Physics/include/Physics.hpp

@@ -7,9 +7,13 @@
 //#include "../src/Discretization.h"
 //#include "../src/Operator.h"
 
+#include "../src/Condition.h"
+#include "../src/FiniteElementMesh.h"
+#include "../src/MatrixGroup.h"
 #include "../src/Node.h"
+#include "../src/Physics.h"
 #include "../src/QuadratureRule.h"
+#include "../src/Region.h"
 #include "../src/Triangle.h"
-#include "../src/FiniteElementMesh.h"
 
 #endif //OERSTED_PHYSICS_HPP

src/Physics/src/Condition.cpp

+//
+// Created by jpries on 1/2/17.
+//
+
+#include "Condition.h"

src/Physics/src/Condition.h

+#ifndef OERSTED_CONDITION_H
+#define OERSTED_CONDITION_H
+
+#include <functional>
+#include <vector>
+
+#include "Eigen"
+
+template<size_t Dimension>
+class HomogeneousForcing {
+};
+
+template<>
+class HomogeneousForcing<2> {
+public:
+    HomogeneousForcing(std::function<double(double)> ff) : F{ff} {};
+
+    double operator()(double const t) { return F(t); };
+
+protected:
+    std::function<double(double)> F;
+
+    //std::vector<size_t> Regions;
+};
+
+#endif //OERSTED_CONDITION_H

src/Physics/src/Discretization.cpp

-#include "Discretization.h"

src/Physics/src/Discretization.h

-#ifndef OERSTED_DISCRETIZATION_H
-#define OERSTED_DISCRETIZATION_H
-
-#include "PhysicsCommon.h"
-
-template<Dimension d> class Discretization {};
-
-template<Order M> class Triangle {
-public:
-    size_t Nodes[(M+2)*(M+3)/2];
-};
-
-template<Dimension d, Order M> class FiniteElementMesh : public Discretization<d> {};
-
-template<Order M> class FiniteElementMesh<2,M> : public Discretization<2> {
-public:
-    FiniteElementMesh() : Discretization() {};
-
-    std::vector<Triangle<M>> Triangles;
-};
-
-#endif //OERSTED_DISCRETIZATION_H
\ No newline at end of file

src/Physics/src/FiniteElementMesh.h

@@ -3,42 +3,10 @@
 
 #include <vector>
 
-#include "Eigen"
-
+#include "MatrixGroup.h"
 #include "Node.h"
 #include "Triangle.h"
 
-template<size_t M, size_t Q>
-class QuadratureDerivatives {
-private:
-    SmallVector<M> Value[Q];
-};
-
-template<size_t Q>
-class QuadratureDerivatives<2, Q> {
-public:
-    double dx(size_t i) { return d[i].x(); };
-
-    double dy(size_t i) { return d[i].y(); };
-
-private:
-    XY d[Q];
-};
-
-template<size_t Q>
-class QuadratureDerivatives<3, Q> {
-public:
-    double dx(size_t i) { return d[i].x(); };
-
-    double dy(size_t i) { return d[i].y(); };
-
-    double dz(size_t i) { return d[i].z(); };
-
-private:
-    XYZ d[Q];
-};
-
-
 template<size_t D, size_t P>
 class FiniteElementMesh {};
 
@@ -53,23 +21,34 @@ public:
 
     std::vector<Triangle<P>> const triangles() const { return Triangles; };
 
+    template<size_t Q>
+    DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> determinant() const {
+        DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Triangles.size());
+
+        for (size_t i = 0; i != Triangles.size(); ++i) {
+            Triangles[i].determinant<Q>(mat, Nodes);
+        }
+
+        return mat;
+    }
+
     template<size_t Q>
     SparseMatrixGroup<TriangleQuadratureRule<Q>::size> basis() const {
-        SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<P>::N);
+        SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
 
         for (size_t i = 0; i != Triangles.size(); ++i) {
-            Triangles[i].basis<Q>(mat, Nodes, i);
+            Triangles[i].basis<Q>(mat, Nodes);
         }
 
         return mat;
-    };
+    }
 
     template<size_t Q>
     DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> derivative() const {
-        DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<P>::N);
+        DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<P>::NumNodes);
 
         for (size_t i = 0; i!= Triangles.size();++i) {
-            Triangles[i].derivative<Q>(df, Nodes, i);
+            Triangles[i].derivative<Q>(df, Nodes);
         }
 
         return df;

src/Physics/src/MatrixGroup.cpp

+//
+// Created by jpries on 1/2/17.
+//
+
+#include "MatrixGroup.h"

src/Physics/src/MatrixGroup.h

+#ifndef OERSTED_MATRIXGROUP_H
+#define OERSTED_MATRIXGROUP_H
+
+#include "Eigen"
+#include "Eigen/Sparse"
+
+// TODO: template<size_t Q, size_t D> MatrixGroup {};
+// TODO: using template<size_t Q> BasisMatrixGroup = MatrixGroup<Q,1>;
+// TODO: using template<size_t Q> DerivativeMatrixGroup = MatrixGroup<Q,2>;
+
+template<size_t Q>
+class SparseMatrixGroup {
+public:
+    SparseMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) {
+        for (size_t i = 0; i != Q; ++i) {
+            Matrices[i].resize(rows, cols);
+            Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
+        }
+    };
+
+    double &operator()(size_t const q, size_t const i, size_t const j) { return Matrices[q].coeffRef(i, j); };
+
+    auto &operator[](size_t const q) { return Matrices[q]; };
+
+    auto &operator[](size_t const q) const { return Matrices[q]; };
+
+protected:
+    std::array<Eigen::SparseMatrix<double>, Q> Matrices;
+};
+
+template<size_t Q>
+class DiagonalMatrixGroup {
+public:
+    DiagonalMatrixGroup(size_t const dim) {
+        for (size_t i = 0; i != Q; ++i) {
+            Matrices[i].resize(dim);
+        }
+    }
+
+    auto &operator()(size_t q) { return Matrices[q]; };
+
+protected:
+    std::array<Eigen::DiagonalMatrix<double, Eigen::Dynamic>, Q> Matrices;
+};
+
+template<size_t Q>
+class DerivativeMatrixGroup {
+public:
+    DerivativeMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) : Dx{rows, cols, nnz}, Dy{rows, cols, nnz} {};
+
+    auto &dx(size_t const q) { return Dx[q]; };
+
+    auto &dy(size_t const q) { return Dy[q]; };
+
+    auto &dx(size_t const q, size_t const i, size_t const j) { return Dx(q, i, j); };
+
+    auto &dy(size_t const q, size_t const i, size_t const j) { return Dy(q, i, j); };
+
+protected:
+    SparseMatrixGroup<Q> Dx;
+    SparseMatrixGroup<Q> Dy;
+};
+
+#endif //OERSTED_MATRIXGROUP_H

src/Physics/src/Physics.h

@@ -3,20 +3,78 @@
 
 #include <memory>
 
-#include "PhysicsCommon.h"
-#include "Operator.h"
-#include "Discretization.h"
+#include "Node.h"
+#include "Triangle.h"
+#include "QuadratureRule.h"
+#include "FiniteElementMesh.h"
+
+enum class Variable {
+    MagneticVectorPotential = 1,
+    A = 1,
+    MagneticFluxDensity = 2,
+    B = 2,
+    ElectricFluxDensity = 4,
+    D = 4,
+    ElectricFieldIntensity = 5,
+    E = 5,
+    MagneticFieldIntensity = 8,
+    H = 8,
+    ElectricCurrentDensity = 10,
+    J = 10,
+    ElectricScalarPotential = 26 + 21,
+    Phi = 26 + 21,
+    MagneticScalarPotential = 26 + 23,
+    Psi = 26 + 23
+};
+
+template<size_t Dimension, size_t ElementOrder, Variable V>
+class Magnetostatic {
+};
+
+template<size_t ElementOrder>
+class Magnetostatic<2, ElementOrder, Variable::MagneticVectorPotential> {
+public:
+    static constexpr size_t QuadratureOrder = 2 * ElementOrder - 1;
+    static constexpr size_t QuadraturePoints = TriangleQuadratureRule<QuadratureOrder>::size;
+
+    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} {};
+
+    Eigen::VectorXd init_residual() { return Eigen::VectorXd(Domain.nodes().size()); };
+
+    Eigen::SparseMatrix<double> init_jacobian() { return Eigen::SparseMatrix<double>(Domain.nodes().size(), Domain.nodes().size()); };
+
+    void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) {
+        J.setZero();
+        r.setZero(); // TODO: Unecessary for r?
+
+        for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
+            J += (Derivatives.dx(i) * ElementWeights(i) * Derivatives.dx(i).transpose())
+                 + (Derivatives.dy(i) * ElementWeights(i) * Derivatives.dy(i).transpose());
+        }
+    };
+
+    void build_matrices() {
+        ElementWeights = Domain.template determinant<QuadratureOrder>();
+        Basis = Domain.template basis<QuadratureOrder>();
+        Derivatives = Domain.template derivative<QuadratureOrder>();
+
+        for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
+            ElementWeights(i) = ElementWeights(i) * TriangleQuadratureRule<QuadratureOrder>::w[i];
+        }
+    };
+
+    void apply_conditions() {};
 
-template<class op_p, class op_d, Dimension d>
-class Physics {
 public:
-    Physics(std::shared_ptr<Discretization<d>> dsc) : PrimalOperator(std::make_shared<op_p>(dsc)), DualOperator(std::make_shared<op_d>(dsc)), Domain(dsc) {};
+    FiniteElementMesh<2, ElementOrder> Domain;
 
-protected:
-    std::shared_ptr<Operator<Primal, d>> PrimalOperator;
-    std::shared_ptr<Operator<Dual, d>> DualOperator;
+    DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
+    SparseMatrixGroup<QuadraturePoints> Basis;
+    DerivativeMatrixGroup<QuadraturePoints> Derivatives;
 
-    std::shared_ptr<Discretization<d>> Domain;
 };
 
 #endif //OERSTED_PHYSICS_H

src/Physics/src/QuadratureRule.h

@@ -10,7 +10,7 @@ struct TriangleQuadratureRule<1> {
     static constexpr size_t size{1};
     static constexpr double a[1]{1.0 / 3.0};
     static constexpr double b[1]{1.0 / 3.0};
-    static constexpr double w[1]{1.0};
+    static constexpr double w[1]{1.0 / 2.0};
 };
 constexpr double TriangleQuadratureRule<1>::a[];
 constexpr double TriangleQuadratureRule<1>::b[];
@@ -21,7 +21,7 @@ struct TriangleQuadratureRule<2> {
     static constexpr size_t size{3};
     static constexpr double a[3]{2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0};
     static constexpr double b[3]{1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0};
-    static constexpr double w[3]{1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+    static constexpr double w[3]{1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0};
 };
 constexpr double TriangleQuadratureRule<2>::a[];
 constexpr double TriangleQuadratureRule<2>::b[];

src/Physics/src/Region.cpp

+#include "Region.h"

src/Physics/src/Region.h

+#ifndef OERSTED_REGION_H
+#define OERSTED_REGION_H
+
+template<size_t Dimension>
+class Region {};
+
+template<>
+class Region<2> {
+public:
+    Region(std::vector<size_t> tris) : Triangles(tris) {};
+
+protected:
+    std::vector<size_t> Triangles;
+};
+
+#endif //OERSTED_REGION_H

src/Physics/src/Triangle.h

@@ -4,71 +4,57 @@
 #include "Eigen"
 #include "Eigen/Sparse"
 
-#include "QuadratureRule.h"
+#include "MatrixGroup.h"
 #include "Node.h"
+#include "QuadratureRule.h"
 
-template<size_t Q>
-class SparseMatrixGroup {
-public:
-    SparseMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) {
-        for (size_t i = 0; i != Q; ++i) {
-            Matrices[i].resize(rows, cols);
-            Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
-        }
-    };
-
-    double &operator()(size_t const q, size_t const i, size_t const j) { return Matrices[q].coeffRef(i, j); };
+// TODO: Curved elements
 
-    auto &operator[](size_t const q) { return Matrices[q]; };
-    auto &operator[](size_t const q) const { return Matrices[q]; };
+template<size_t Dimension>
+class Element {
+public:
+    Element() : ID{0} {};
+    Element(size_t id) : ID{id} {};
 
 protected:
-    std::array<Eigen::SparseMatrix<double>, Q> Matrices;
+    size_t ID;
 };
 
-template<size_t Q>
-class DerivativeMatrixGroup {
+class Facet : public Element<2> {
 public:
-    DerivativeMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) : Dx{rows, cols, nnz}, Dy{rows, cols, nnz} {};
-
-    auto &dx(size_t const q) { return Dx[q]; };
-    auto &dy(size_t const q) { return Dy[q]; };
-
-    auto &dx(size_t const q, size_t const i, size_t const j) { return Dx(q,i,j); };
-    auto &dy(size_t const q, size_t const i, size_t const j) { return Dy(q,i,j); };
-
-protected:
-    SparseMatrixGroup<Q> Dx;
-    SparseMatrixGroup<Q> Dy;
+    Facet() : Element{} {};
+    Facet(size_t id) : Element{id} {};
 };
 
-// TODO: Curved elements
-
 template<size_t P>
-class Triangle {
+class Triangle : public Facet {
 public:
-    static constexpr size_t N{(P + 1) * (P + 2) / 2};
+    static constexpr size_t NumNodes{(P + 1) * (P + 2) / 2};
 
-    Triangle() : Node{} {};
+    Triangle() : Facet{}, Node{} {};
 
     size_t const &node(size_t const &i) const { return Node[i]; };
 
     template<typename... Args>
-    Triangle(Args &&... args) : Node{std::forward<Args>(args)...} {};
+    Triangle(size_t id, Args &&... args) : Facet{id}, Node{std::forward<Args>(args)...} {};
 
     template<size_t D>
     Eigen::Matrix<double, 2, 2> jacobian(std::vector<XY> const &nodes) const; // TODO: Accept Eigen::Matrix<double,2,2> as input
 
     template<size_t Q>
-    void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const;
+    void determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const;
 
     template<size_t Q>
-    void derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes, size_t const &tri) const;
+    void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const;
+
+    template<size_t Q>
+    void derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes) const;
 
 protected:
-    size_t Node[N];
+    size_t Node[NumNodes];
 };
 
+/*
 template<>
 template<>
 Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &nodes) const {
@@ -81,6 +67,7 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &node
 
     return value;
 }
+*/
 
 template<>
 template<>
@@ -107,28 +94,44 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &node
 
 template<>
 template<size_t Q>
-void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const {
+void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const {
+    for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
+        XY const &p0 = nodes[Node[0]];
+        XY const &p1 = nodes[Node[1]];
+        XY const &p2 = nodes[Node[2]];
+
+        double xx = p0.x() - p2.x();
+        double xy = p0.y() - p2.y();
+        double yx = p1.x() - p2.x();
+        double yy = p1.y() - p2.y();
+        mat(i).diagonal()(ID) = xx * yy - xy * yx;
+    }
+}
+
+template<>
+template<size_t Q>
+void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const {
     for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
-        mat(i, Node[0], tri) += TriangleQuadratureRule<Q>::a[i];
-        mat(i, Node[1], tri) += TriangleQuadratureRule<Q>::b[i];
-        mat(i, Node[2], tri) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
+        mat(i, Node[0], ID) += TriangleQuadratureRule<Q>::a[i];
+        mat(i, Node[1], ID) += TriangleQuadratureRule<Q>::b[i];
+        mat(i, Node[2], ID) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
     }
 }
 
 template<>
 template<size_t Q>
-void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes, size_t const &tri) const {
+void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes) const {
     auto J = jacobian<1>(nodes);
 
     for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
-        df.dx(i, Node[0], tri) += J(0, 0);
-        df.dy(i, Node[0], tri) += J(1, 0);
+        df.dx(i, Node[0], ID) += J(0, 0);
+        df.dy(i, Node[0], ID) += J(1, 0);
 
-        df.dx(i, Node[1], tri) += J(0, 1);
-        df.dy(i, Node[1], tri) += J(1, 1);
+        df.dx(i, Node[1], ID) += J(0, 1);
+        df.dy(i, Node[1], ID) += J(1, 1);
 
-        df.dx(i, Node[2], tri) += (-J(0, 0) - J(0, 1));
-        df.dy(i, Node[2], tri) += (-J(1, 0) - J(1, 1));
+        df.dx(i, Node[2], ID) += (-J(0, 0) - J(0, 1));
+        df.dy(i, Node[2], ID) += (-J(1, 0) - J(1, 1));
     }
 };
 

test/Physics/test_FiniteElementMesh.cpp

@@ -4,7 +4,82 @@ inline constexpr size_t operator "" _zu(unsigned long long int i) {
     return (size_t) i;
 }
 
-class Master_Triangle : public ::testing::Test {
+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});
+
+        tri.push_back(Triangle<1>(0, 1_zu, 2_zu, 0_zu));
+
+        femesh = FiniteElementMesh<2,1>(node,tri);
+    }
+
+    std::vector<XY> node;
+    std::vector<Triangle<1>> tri;
+
+    FiniteElementMesh<2,1> femesh;
+};
+
+TEST_F(MasterTriangle, mass_matrix) {
+    auto mat = femesh.basis<2>();
+    auto det = femesh.determinant<2>();
+
+    Eigen::SparseMatrix<double> M = (mat[0] * det(0) * mat[0].transpose()) * TriangleQuadratureRule<2>::w[0]
+                                    + (mat[1] * det(1) * mat[1].transpose()) * TriangleQuadratureRule<2>::w[1]
+                                    + (mat[2] * det(2) * mat[2].transpose()) * TriangleQuadratureRule<2>::w[2];
+
+    EXPECT_NEAR(M.coeffRef(0, 0), 1.0 / 12.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(1, 0), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(2, 0), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(0, 1), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(1, 1), 1.0 / 12.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(2, 1), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(0, 2), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(1, 2), 1.0 / 24.0, FLT_EPSILON);
+    EXPECT_NEAR(M.coeffRef(2, 2), 1.0 / 12.0, FLT_EPSILON);
+}
+
+TEST_F(MasterTriangle, stiffness_matrix) {