Start of mesh conversion rewrite to enable higher order elements

Reorganize some code into separate shared libraries

CMakeLists.txt

@@ -17,6 +17,9 @@ include_directories(./src/)
 include_directories(./src/Sketch/include/)
 include_directories(./src/Mesh/include/)
 include_directories(./src/Quadrature/include/)
+include_directories(./src/Elements/include/)
+include_directories(./src/Matrix/include/)
+include_directories(./src/Materials/include/)
 include_directories(./src/Physics/include/)
 
 add_subdirectory(./src/)

src/CMakeLists.txt

 include_directories(./Sketch/)
 include_directories(./Mesh/)
 include_directories(./Quadrature/)
+include_directories(./Elements/)
+include_directories(./Matrix/)
+include_directories(./Materials/)
 include_directories(./Physics/)
 
+
 add_subdirectory(./Sketch/)
 add_subdirectory(./Mesh/)
 add_subdirectory(./Quadrature/)
+add_subdirectory(./Elements/)
+add_subdirectory(./Matrix/)
+add_subdirectory(./Materials/)
 add_subdirectory(./Physics/)
\ No newline at end of file

src/Elements/CMakeLists.txt

+PROJECT(Oersted_Elements)
+
+set(SOURCE_FILES
+        ./include/Elements.hpp
+
+        ./src/Triangle.cpp  ./src/Triangle.h
+        ./src/Node.cpp      ./src/Node.h
+        ./src/Element.cpp   ./src/Element.h
+        ./src/Facet.cpp     ./src/Facet.h
+        ../Materials/include/Materials.hpp)
+
+add_library(elements SHARED ${SOURCE_FILES})
+
+target_link_libraries(elements)
\ No newline at end of file

src/Elements/include/Elements.hpp

+#ifndef OERSTED_ELEMENTS_HP
+#define OERSTED_ELEMENTS_HP
+
+#include "../src/Node.h"
+#include "../src/Triangle.h"
+
+#endif //OERSTED_ELEMENTS_HPP

src/Elements/src/Element.cpp

+//
+// Created by p7k on 8/23/17.
+//
+
+#include "Element.h"

src/Elements/src/Element.h

+#ifndef OERSTED_ELEMENT_H
+#define OERSTED_ELEMENT_H
+
+#include <cstddef>
+
+template<size_t Dimension>
+class Element {
+public:
+    Element() : ID{0} {};
+
+    Element(size_t id) : ID{id} {};
+
+    size_t const id() const { return ID; };
+
+    void id(size_t i) { ID = i; };
+
+protected:
+    size_t ID;
+};
+
+#endif //OERSTED_ELEMENT_H

src/Elements/src/Facet.cpp

+//
+// Created by p7k on 8/23/17.
+//
+
+#include "Facet.h"

src/Elements/src/Facet.h

+#ifndef OERSTED_FACET_H
+#define OERSTED_FACET_H
+
+#include "Element.h"
+
+class Facet : public Element<2> {
+public:
+    Facet() : Element{} {};
+
+    Facet(size_t id) : Element{id} {};
+};
+
+#endif //OERSTED_FACET_H

src/Physics/src/Node.cpp → src/Elements/src/Node.cpp

+//
+// Created by p7k on 8/23/17.
+//
+
 #include "Node.h"

src/Elements/src/Node.h

+#ifndef OERSTED_NODE_H
+#define OERSTED_NODE_H
+
+// TODO: Const
+
+#include <cstddef>
+
+#include "Mesh.hpp"
+
+// TODO: Combine different node/point classes from Mesh and Elements
+
+class XY {
+public:
+    XY() : Value{} {};
+
+    XY(double_t x, double_t y) : Value{x, y} {};
+
+    XY(XY const &v) : Value{v.x(), v.y()} {};
+
+    XY(Point const& p) : Value{p.X, p.Y} {};
+
+    double_t value(size_t i) const { return Value[i]; };
+
+    double_t x() const { return Value[0]; };
+
+    void x(double_t val) { Value[0] = val; };
+
+    double_t y() const { return Value[1]; };
+
+    void y(double_t val) { Value[1] = val; };
+
+protected:
+    double_t Value[2];
+};
+
+#endif //OERSTED_NODE_H
\ No newline at end of file

src/Elements/src/Triangle.cpp

+#include "Triangle.h"

src/Physics/src/Triangle.h → src/Elements/src/Triangle.h

@@ -3,36 +3,14 @@
 
 #include <vector>
 
-#include "Eigen"
-#include "Eigen/Sparse"
-
+#include "Matrix.hpp"
 #include "Quadrature.hpp"
 
-#include "MatrixGroup.h"
 #include "Node.h"
+#include "Facet.h"
 
 // TODO: Curved elements
 
-template<size_t Dimension>
-class Element {
-public:
-    Element() : ID{0} {};
-    Element(size_t id) : ID{id} {};
-
-    size_t const id() const { return ID; };
-
-    void id(size_t i) { ID = i; };
-
-protected:
-    size_t ID;
-};
-
-class Facet : public Element<2> {
-public:
-    Facet() : Element{} {};
-    Facet(size_t id) : Element{id} {};
-};
-
 template<size_t P>
 class Triangle : public Facet {
 public:
@@ -40,18 +18,14 @@ public:
 
     Triangle() : Facet{}, Node{} {};
 
-    Triangle(size_t id, std::array<size_t, NumNodes> const &n) : Facet{id} {
-        for (size_t i = 0; i != NumNodes; ++i) {
-            Node[i] = n[i];
-        }
-    };
+    Triangle(size_t id, std::array<size_t, NumNodes> const &n) : Facet{id}, Node{n} {};
 
     void node(size_t const &i, size_t n) { Node[i] = n; };
 
     size_t const &node(size_t const &i) const { return Node[i]; };
 
     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
+    Oe::Matrix<double_t, 2, 2> jacobian(std::vector<XY> const &nodes) const; // TODO: Accept Eigen::Matrix<double,2,2> as input
 
     template<size_t Q>
     void determinant(DiagonalMatrixGroup &mat, std::vector<XY> const &nodes) const;
@@ -66,38 +40,23 @@ public:
     std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const;
 
 protected:
-    size_t Node[NumNodes]; // TODO: std::array instead?
+    std::array<size_t, NumNodes> Node;
 };
 
-/*
 template<>
 template<>
-inline Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &nodes) const {
-    Eigen::Matrix<double, 2, 2> value;
-
-    value(0, 0) = 1.0;
-    value(0, 1) = 0.0;
-    value(1, 0) = 0.0;
-    value(1, 1) = 1.0;
-
-    return value;
-}
-*/
-
-template<>
-template<>
-inline Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
-    Eigen::Matrix<double, 2, 2> value;
+inline Oe::Matrix<double_t, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
+    Oe::Matrix<double_t, 2, 2> value;
 
     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();
-    double det = xx * yy - xy * yx;
+    double_t xx = p0.x() - p2.x();
+    double_t xy = p0.y() - p2.y();
+    double_t yx = p1.x() - p2.x();
+    double_t yy = p1.y() - p2.y();
+    double_t det = xx * yy - xy * yx;
 
     value(0, 0) = yy / det;
     value(0, 1) = -xy / det;
@@ -115,18 +74,19 @@ void Triangle<1>::determinant(DiagonalMatrixGroup &mat, std::vector<XY> const &n
         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();
+        double_t xx = p0.x() - p2.x();
+        double_t xy = p0.y() - p2.y();
+        double_t yx = p1.x() - p2.x();
+        double_t yy = p1.y() - p2.y();
         mat(i)(ID) = xx * yy - xy * yx;
     }
 }
 
 
+// TODO? Return Triplets instead of modifying matrix directly
+
 template<>
 template<size_t Q>
-//void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const {
 void Triangle<1>::basis(SparseMatrixGroup &mat, std::vector<XY> const &nodes) const {
     for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
         mat(i, Node[0], ID) += TriangleQuadrature<Q>::a[i];

src/Materials/CMakeLists.txt

+PROJECT(Oersted_Materials)
+
+set(SOURCE_FILES
+        ./include/Materials.hpp
+
+        ./src/PhysicsConstants.h
+
+        ./src/MaterialProperties.cpp ./src/MaterialProperties.h
+        )
+
+add_library(materials SHARED ${SOURCE_FILES})
+
+target_link_libraries(materials)
\ No newline at end of file

src/Materials/include/Materials.hpp

+#ifndef OERSTED_MATERIALS_HPP
+#define OERSTED_MATERIALS_HPP
+
+#include "../src/PhysicsConstants.h"
+#include "../src/MaterialProperties.h"
+
+#endif //OERSTED_MATERIALS_HPP

src/Physics/src/MaterialProperties.h → src/Materials/src/MaterialProperties.h

@@ -7,7 +7,7 @@
 #include <utility>
 #include <vector>
 
-#include "Eigen"
+#include "Matrix.hpp"
 
 #include "PhysicsConstants.h"
 
@@ -15,14 +15,14 @@ class MagneticMaterialInterface {
 public:
     virtual ~MagneticMaterialInterface(){};
 
-    virtual void h_from_b(std::vector<size_t> const &index, Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) = 0;
+    virtual void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) = 0;
 };
 
 class LinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
 public:
     LinearIsotropicMagneticMaterial(double relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};
 
-    void h_from_b(std::vector<size_t> const &index, Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) override {
+    void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) override {
         for(size_t const &i : index) {
             Fx(i) *= Reluctivity;
             Fy(i) *= Reluctivity;
@@ -74,7 +74,7 @@ public:
         MSaturation = (s[0] * BSaturation + s[1]) * BSaturation + s[2];
     }
 
-    void h_from_b(std::vector<size_t> const &index, Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) override {
+    void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) override {
         for(size_t const &i : index) {
             const double_t &Bx = Fx(i);
             const double_t &By = Fy(i);

src/Matrix/CMakeLists.txt

+PROJECT(Oersted_Matrix)
+
+set(SOURCE_FILES
+        ./include/Matrix.hpp
+
+        ./src/Alias.h
+
+        ./src/MatrixGroup.cpp               ./src/MatrixGroup.h
+        ./src/SparseMatrixGroup.cpp         ./src/SparseMatrixGroup.h
+        ./src/DiagonalMatrixGroup.cpp       ./src/DiagonalMatrixGroup.h
+        ./src/DerivativeMatrixGroup.cpp     ./src/DerivativeMatrixGroup.h
+
+        )
+
+add_library(matrix SHARED ${SOURCE_FILES})
+
+target_link_libraries(matrix)
\ No newline at end of file

src/Matrix/src/Alias.h

+#ifndef OERSTED_ALIAS_H
+#define OERSTED_ALIAS_H
+
+#include "./Alias/EigenAlias.h"
+
+#endif //OERSTED_ALIAS_H

src/Matrix/src/DerivativeMatrixGroup.cpp

+#include "DerivativeMatrixGroup.h"

src/Matrix/src/DerivativeMatrixGroup.h

+#ifndef OERSTED_DERIVATIVEMATRIXGROUP_H
+#define OERSTED_DERIVATIVEMATRIXGROUP_H
+
+#include "MatrixGroup.h"
+#include "SparseMatrixGroup.h"
+
+class DerivativeMatrixGroup {
+public:
+    DerivativeMatrixGroup() {};
+
+    DerivativeMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) : Dx{Q, rows, cols, nnz}, Dy{Q, rows, cols, nnz} {};
+
+    auto const &dx(size_t const q) const { return Dx(q); };
+
+    auto const &dy(size_t const q) const { 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); };
+
+    size_t size() const { return Dx.size(); };
+
+    void transform(Oe::SparseMatrix<double_t> &A) {
+        Dx.transform(A);
+        Dy.transform(A);
+    }
+
+protected:
+    SparseMatrixGroup Dx;
+    SparseMatrixGroup Dy;
+};
+
+#endif //OERSTED_DERIVATIVEMATRIXGROUP_H

src/Matrix/src/DiagonalMatrixGroup.cpp

+#include "DiagonalMatrixGroup.h"

src/Matrix/src/DiagonalMatrixGroup.h

+#ifndef OERSTED_DIAGONALMATRIXGROUP_H
+#define OERSTED_DIAGONALMATRIXGROUP_H
+
+#include "MatrixGroup.h"
+
+class DiagonalMatrixGroup {
+public:
+    DiagonalMatrixGroup() {};
+
+    DiagonalMatrixGroup(size_t Q, size_t const dim) {
+        Matrices.resize(Q);
+        for (size_t i = 0; i != Matrices.size(); ++i) {
+            Matrices[i].resize(dim);
+        }
+    }
+
+    auto &operator()(size_t const q) const { return Matrices[q]; };
+
+    auto &operator()(size_t const q) { return Matrices[q]; };
+
+    size_t size() const { return (size_t)(Matrices[0].size()); };
+
+protected:
+    std::vector<Oe::ArrayXd> Matrices;
+};
+
+#endif //OERSTED_DIAGONALMATRIXGROUP_H

src/Matrix/src/MatrixGroup.h

+#ifndef OERSTED_MATRIXGROUP_H
+#define OERSTED_MATRIXGROUP_H
+
+#include <cstddef>
+#include <cmath>
+
+#include "Alias.h"
+
+// TODO: MatrixGroup interface class
+
+#endif //OERSTED_MATRIXGROUP_H

src/Matrix/src/SparseMatrixGroup.cpp

+#include "SparseMatrixGroup.h"

src/Physics/src/MatrixGroup.h → src/Matrix/src/SparseMatrixGroup.h

-#ifndef OERSTED_MATRIXGROUP_H
-#define OERSTED_MATRIXGROUP_H
+#ifndef OERSTED_SPARSEMATRIXGROUP_H
+#define OERSTED_SPARSEMATRIXGROUP_H
 
-#include "Eigen"
-#include "Eigen/Sparse"
-
-// TODO: MatrixGroup interface class
+#include "MatrixGroup.h"
 
 class SparseMatrixGroup {
 public:
     SparseMatrixGroup() {};
 
-    SparseMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) {
+    SparseMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) { // Quadrature Points, Nodes, Elements, NNZ per Element
         Matrices.resize(Q);
         for (size_t i = 0; i != Matrices.size(); ++i) {
             Matrices[i].resize(rows, cols);
-            Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
+            Matrices[i].reserve(Oe::VectorXd::Constant(cols, nnz));
         }
     };
 
@@ -26,61 +23,14 @@ public:
 
     size_t size() const { return Matrices.size(); };
 
-    void transform(Eigen::SparseMatrix<double> &A) {
+    void transform(Oe::SparseMatrix<double_t> &A) {
         for (size_t i = 0; i != Matrices.size(); ++i) {
             Matrices[i] = A * Matrices[i];
         }
     }
 
 protected:
-    std::vector<Eigen::SparseMatrix<double_t>> Matrices;
-};
-
-class DiagonalMatrixGroup {
-public:
-    DiagonalMatrixGroup() {};
-
-    DiagonalMatrixGroup(size_t Q, size_t const dim) {
-        Matrices.resize(Q);
-        for (size_t i = 0; i != Matrices.size(); ++i) {
-            Matrices[i].resize(dim);
-        }
-    }
-
-    auto &operator()(size_t const q) const { return Matrices[q]; };
-
-    auto &operator()(size_t const q) { return Matrices[q]; };
-
-    size_t size() const { return (size_t)(Matrices[0].size()); };
-
-protected:
-    std::vector<Eigen::ArrayXd> Matrices;
-};
-
-class DerivativeMatrixGroup {
-public:
-    DerivativeMatrixGroup() {};
-
-    DerivativeMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) : Dx{Q, rows, cols, nnz}, Dy{Q, rows, cols, nnz} {};
-
-    auto const &dx(size_t const q) const { return Dx(q); };
-
-    auto const &dy(size_t const q) const { 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); };
-
-    size_t size() const { return Dx.size(); };
-
-    void transform(Eigen::SparseMatrix<double> &A) {
-        Dx.transform(A);
-        Dy.transform(A);
-    }
-
-protected:
-    SparseMatrixGroup Dx;
-    SparseMatrixGroup Dy;
+    std::vector<Oe::SparseMatrix<double_t>> Matrices;
 };
 
-#endif //OERSTED_MATRIXGROUP_H
+#endif //OERSTED_SPARSEMATRIXGROUP_H

src/Mesh/CMakeLists.txt

@@ -10,9 +10,8 @@ set(SOURCE_FILES
         ./src/MappedBoundaryPair.h  ./src/MappedBoundaryPair.cpp
         ./src/Mesh.h                ./src/Mesh.cpp
         ./src/Point.h               ./src/Point.cpp
-        ./src/DartTriangle.h        ./src/DartTriangle.cpp
-        )
+        ./src/DartTriangle.h        ./src/DartTriangle.cpp)
 
 add_library(mesh SHARED ${SOURCE_FILES})
 
-target_link_libraries(mesh ${Boost_LIBRARIES})
\ No newline at end of file
+target_link_libraries(mesh)
\ No newline at end of file

src/Physics/CMakeLists.txt

@@ -8,13 +8,8 @@ set(SOURCE_FILES
         ./src/DiscreteRegion.h ./src/DiscreteRegion.cpp
         ./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
         ./src/Forcing.h ./src/Forcing.cpp
-        ./src/MaterialProperties.h ./src/MaterialProperties.cpp
-        ./src/MatrixGroup.h ./src/MatrixGroup.cpp
         ./src/PostProcessorInterface.h ./src/PostProcessorInterface.cpp
-        ./src/Node.h ./src/Node.cpp
         ./src/Physics.h ./src/Physics.cpp
-        ./src/PhysicsConstants.h
-        ./src/Triangle.h ./src/Triangle.cpp
         ./src/Solution.cpp ./src/Solution.h
         ./src/ZeroDirichlet.cpp ./src/ZeroDirichlet.h
         ./src/PeriodicBoundaryCondition.cpp ./src/PeriodicBoundaryCondition.h
@@ -25,9 +20,8 @@ set(SOURCE_FILES
         ./src/HomogeneousForcing.cpp ./src/HomogeneousForcing.h
         ./src/FunctionArguments.cpp ./src/FunctionArguments.h
         ./src/FiniteElementMeshInterface.cpp ./src/FiniteElementMeshInterface.h
-
         )
 
 add_library(physics SHARED ${SOURCE_FILES})
 
-target_link_libraries(physics ${Boost_LIBRARIES})
\ No newline at end of file
+target_link_libraries(physics)
\ No newline at end of file

src/Physics/src/BoundaryCondition.h

@@ -5,11 +5,10 @@
 #include <set>
 #include <algorithm>
 
-#include <Eigen/Sparse>
+#include "Matrix.hpp"
+#include "Quadrature.hpp"
 
 #include "DiscreteBoundary.h"
-#include "Quadrature.hpp"
-#include "MatrixGroup.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
@@ -17,7 +16,7 @@
 
 class BoundaryCondition {
 public:
-    virtual void apply(Eigen::SparseMatrix<double_t> &bc_matrix) const = 0;
+    virtual void apply(Oe::SparseMatrix<double_t> &bc_matrix) const = 0;
 
     virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
 };

src/Physics/src/DiscreteRegion.h

@@ -4,8 +4,7 @@
 #include <vector>
 
 #include "Mesh.hpp"
-
-#include "MaterialProperties.h"
+#include "Materials.hpp"
 
 class DiscreteRegion {
 public:

src/Physics/src/FiniteElementMesh.cpp

@@ -18,38 +18,7 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
             Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
      */
 
-    // std::vector<XY> nodes
-    // std::vector<Triangle<Order>> tris
-    // std::vector<std::shared_ptr<Region<2>>> r
-    // std::vector<std::shared_ptr<Boundary<2>>> b
-
-    /*
-    Nodes.reserve(m.size_points()); // TODO: This will be larger when ElementOrder > 1
-    for (size_t i = 0; i != m.size_points(); ++i) {