Implement basic functions for second order triangular elements

src/Elements/src/Triangle.h

@@ -24,111 +24,212 @@ public:
 
     size_t const &node(size_t const &i) const { return Node[i]; };
 
-    template<size_t D>
-    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>
+    std::array<double_t, NumNodes> basis(size_t q) const;
+
+    template<size_t Q>
+    std::array<double_t, NumNodes> da(size_t q) const;
 
     template<size_t Q>
-    void determinant(DiagonalMatrixGroup &mat, std::vector<XY> const &nodes) const;
+    std::array<double_t, NumNodes> db(size_t q) const;
 
     template<size_t Q>
-    void basis(SparseMatrixGroup &mat, std::vector<XY> const &nodes) const;
+    Oe::Matrix<double_t, 2, 2> jacobian(std::vector<XY> const &nodes, size_t q) const {
+        Oe::Matrix<double_t, 2, 2> value;
+
+        auto const &daq = da<Q>(q);
+        auto const &dbq = db<Q>(q);
+
+        double_t dxda{0.0};
+        double_t dyda{0.0};
+        double_t dxdb{0.0};
+        double_t dydb{0.0};
+        for (size_t i = 0; i != NumNodes; ++i) {
+            XY const &p = nodes[Node[i]];
+
+            dxda += daq[i] * p.x();
+            dxdb += dbq[i] * p.x();
+
+            dydb += dbq[i] * p.y();
+            dyda += daq[i] * p.y();
+        }
+
+        double_t det = dxda * dydb - dyda * dxdb;
+
+        value(0, 0) = dydb / det;
+        value(0, 1) = -dyda / det;
+        value(1, 0) = -dxdb / det;
+        value(1, 1) = dxda / det;
+
+        return value;
+    }
 
     template<size_t Q>
-    void derivative(DerivativeMatrixGroup &df, std::vector<XY> const &nodes) const;
+    void determinant(DiagonalMatrixGroup &mat, std::vector<XY> const &nodes) const {
+        for (size_t q = 0; q != TriangleQuadrature<Q>::size; ++q) {
+            auto const &daq = da<Q>(q);
+            auto const &dbq = db<Q>(q);
+
+            double_t dxda{0.0};
+            double_t dyda{0.0};
+            double_t dxdb{0.0};
+            double_t dydb{0.0};
+            for (size_t i = 0; i != NumNodes; ++i) {
+                XY const &p = nodes[Node[i]];
+
+                dxda += daq[i] * p.x();
+                dxdb += dbq[i] * p.x();
+
+                dydb += dbq[i] * p.y();
+                dyda += daq[i] * p.y();
+            }
+
+            mat(q)(ID) = dxda * dydb - dyda * dxdb;
+        }
+    }
 
     template<size_t Q>
-    std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const;
+    void basis(SparseMatrixGroup &mat, std::vector<XY> const &nodes) const {
+        for (size_t q = 0; q != TriangleQuadrature<Q>::size; ++q) {
+            auto const &value = basis<Q>(q);
+            for (size_t i = 0; i != NumNodes; ++i) {
+                mat(q, Node[i], ID) += value[i];
+            }
+        }
+    }
+
+    template<size_t Q>
+    void derivative(DerivativeMatrixGroup &df, std::vector<XY> const &nodes) const {
+        for (size_t q = 0; q != TriangleQuadrature<Q>::size; ++q) {
+            auto J = jacobian<Q>(nodes, q);
+
+            auto const &daq = da<Q>(q);
+            auto const &dbq = db<Q>(q);
+
+            for(size_t i = 0; i != NumNodes; ++i) {
+                df.dx(q, Node[i], ID) += J(0,0) * daq[i] + J(0,1) * dbq[i];
+                df.dy(q, Node[i], ID) += J(1,0) * daq[i] + J(1,1) * dbq[i];
+            }
+        }
+    }
+
+    template<size_t Q>
+    std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const {
+        std::vector<XY> qp;
+        qp.reserve(TriangleQuadrature<Q>::size);
+
+        for (size_t q = 0; q != TriangleQuadrature<Q>::size; ++q) {
+            auto const &b = basis<Q>(q);
+            double_t x{0.0};
+            double_t y{0.0};
+
+            for (size_t i = 0; i != NumNodes; ++i) {
+                auto const &p = nodes[Node[i]];
+                x += b[i] * p.x();
+                y += b[i] * p.y();
+            }
+
+            qp.emplace_back(x,y);
+        }
+
+        return qp;
+    }
 
 protected:
     std::array<size_t, NumNodes> Node;
 };
 
 template<>
-template<>
-inline Oe::Matrix<double_t, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
-    Oe::Matrix<double_t, 2, 2> value;
+template<size_t Q>
+inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::basis(size_t q) const {
+    std::array<double_t, NumNodes> value;
+
+    double_t const &a = TriangleQuadrature<Q>::a[q];
+    double_t const &b = TriangleQuadrature<Q>::b[q];
+
+    value[0] = a;
+    value[1] = b;
+    value[2] = 1.0 - a - b;
+
+    return value;
+}
 
-    XY const &p0 = nodes[Node[0]];
-    XY const &p1 = nodes[Node[1]];
-    XY const &p2 = nodes[Node[2]];
+template<>
+template<size_t Q>
+inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::basis(size_t q) const {
+    std::array<double_t, NumNodes> value;
 
-    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;
+    double_t const &a = TriangleQuadrature<Q>::a[q];
+    double_t const &b = TriangleQuadrature<Q>::b[q];
 
-    value(0, 0) = yy / det;
-    value(0, 1) = -xy / det;
-    value(1, 0) = -yx / det;
-    value(1, 1) = xx / det;
+    value[0] = 2.0 * (1.0 - a - b) * (0.5 - a - b);
+    value[1] = 4.0 * (1.0 - a - b) * a;
+    value[2] = 2.0 * (a - 0.5) * a;
+    value[3] = 4.0 * a * b;
+    value[4] = 2.0 * (b - 0.5) * b;
+    value[5] = 4.0 * (1.0 - a - b) * b;
 
     return value;
 };
 
 template<>
 template<size_t Q>
-void Triangle<1>::determinant(DiagonalMatrixGroup &mat, std::vector<XY> const &nodes) const {
-    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
-        XY const &p0 = nodes[Node[0]];
-        XY const &p1 = nodes[Node[1]];
-        XY const &p2 = nodes[Node[2]];
-
-        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;
-    }
-}
+inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::da(size_t q) const {
+    std::array<double_t, NumNodes> value;
 
+    value[0] = 1.0;
+    value[1] = 0.0;
+    value[2] = -1.0;
 
-// TODO? Return Triplets instead of modifying matrix directly
+    return value;
+}
 
 template<>
 template<size_t Q>
-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];
-        mat(i, Node[1], ID) += TriangleQuadrature<Q>::b[i];
-        mat(i, Node[2], ID) += 1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i];
-    }
+inline std::array<double_t, Triangle<1>::NumNodes> Triangle<1>::db(size_t q) const {
+    std::array<double_t, NumNodes> value;
+
+    value[0] = 0.0;
+    value[1] = 1.0;
+    value[2] = -1.0;
+
+    return value;
 }
 
 template<>
 template<size_t Q>
-void Triangle<1>::derivative(DerivativeMatrixGroup &df, std::vector<XY> const &nodes) const {
-    auto J = jacobian<1>(nodes);
+inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::da(size_t q) const {
+    std::array<double_t, NumNodes> value;
 
-    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
-        df.dx(i, Node[0], ID) += J(0, 0);
-        df.dy(i, Node[0], ID) += J(1, 0);
+    double_t const &a = TriangleQuadrature<Q>::a[q];
+    double_t const &b = TriangleQuadrature<Q>::b[q];
 
-        df.dx(i, Node[1], ID) += J(0, 1);
-        df.dy(i, Node[1], ID) += J(1, 1);
+    value[0] = +4.0 * a + 4.0 * b - 3.0;
+    value[1] = +4.0 - 4.0 * b - 8.0 * a;
+    value[2] = +4.0 * a - 1.0;
+    value[3] = +4.0 * b;
+    value[4] = +0.0;
+    value[5] = -4.0 * b;
 
-        df.dx(i, Node[2], ID) += (-J(0, 0) - J(0, 1));
-        df.dy(i, Node[2], ID) += (-J(1, 0) - J(1, 1));
-    }
+    return value;
 }
 
 template<>
 template<size_t Q>
-std::vector<XY> Triangle<1>::quadrature_points(std::vector<XY> const& nodes) const {
-    std::vector<XY> qp;
-    qp.reserve(TriangleQuadrature<Q>::size);
-
-    XY const &p0 = nodes[Node[0]];
-    XY const &p1 = nodes[Node[1]];
-    XY const &p2 = nodes[Node[2]];
-
-    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
-        double_t x = p0.x() * TriangleQuadrature<Q>::a[i] + p1.x() * TriangleQuadrature<Q>::b[i] + p2.x() * (1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i]);
-        double_t y = p0.y() * TriangleQuadrature<Q>::a[i] + p1.y() * TriangleQuadrature<Q>::b[i] + p2.y() * (1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i]);
-        qp.emplace_back(x,y);
-    }
+inline std::array<double_t, Triangle<2>::NumNodes> Triangle<2>::db(size_t q) const {
+    std::array<double_t, NumNodes> value;
+
+    double_t const &a = TriangleQuadrature<Q>::a[q];
+    double_t const &b = TriangleQuadrature<Q>::b[q];
 
-    return qp;
+    value[0] = +4.0 * a + 4.0 * b - 3.0;
+    value[1] = -4.0 * a;
+    value[2] = +0.0;
+    value[3] = +4.0 * a;
+    value[4] = +4.0 * b - 1.0;
+    value[5] = +4.0 - 8.0 * b - 4.0 * a;
+
+    return value;
 }
 
 #endif //OERSTED_TRIANGLE_H
\ No newline at end of file

src/Matrix/src/DerivativeMatrixGroup.h

@@ -8,7 +8,7 @@ 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} {};
+    DerivativeMatrixGroup(size_t const qpts, size_t const nodes, size_t const nels, size_t const el_nnz) : Dx{qpts, nodes, nels, el_nnz}, Dy{qpts, nodes, nels, el_nnz} {};
 
     auto const &dx(size_t const q) const { return Dx(q); };
 

src/Matrix/src/DiagonalMatrixGroup.h

@@ -7,10 +7,10 @@ class DiagonalMatrixGroup {
 public:
     DiagonalMatrixGroup() {};
 
-    DiagonalMatrixGroup(size_t Q, size_t const dim) {
-        Matrices.resize(Q);
+    DiagonalMatrixGroup(size_t const qpts, size_t const nels) {
+        Matrices.resize(qpts);
         for (size_t i = 0; i != Matrices.size(); ++i) {
-            Matrices[i].resize(dim);
+            Matrices[i].resize(nels);
         }
     }
 

src/Matrix/src/SparseMatrixGroup.h

@@ -7,11 +7,11 @@ class SparseMatrixGroup {
 public:
     SparseMatrixGroup() {};
 
-    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);
+    SparseMatrixGroup(size_t const qpts, size_t const nodes, size_t const nels, size_t const el_nnz) { // Quadrature Points, Nodes, Elements, NNZ per Element
+        Matrices.resize(qpts);
         for (size_t i = 0; i != Matrices.size(); ++i) {
-            Matrices[i].resize(rows, cols);
-            Matrices[i].reserve(Oe::VectorXd::Constant(cols, nnz));
+            Matrices[i].resize(nodes, nels);
+            Matrices[i].reserve(Oe::VectorXd::Constant(nels, el_nnz));
         }
     };
 

src/Mesh/src/Mesh.h

@@ -234,6 +234,8 @@ public:
     };
 
     MeshData mesh_data(size_t p) const {
+        // TODO: Modify this to return triangles and nodes instead of std::vector<std::vector<std::vector<size_t>>>
+
         std::vector<NodeData> data;
         std::vector<Point> points;
         std::vector<std::vector<std::vector<size_t>>> edges;

src/Physics/src/FiniteElementMesh.cpp

@@ -18,7 +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)
      */
 
-    //TODO: CLEAN ME and MESH
+    // TODO: CLEAN ME and MESH
 
     MeshData md = m.mesh_data(P);
 

test/Elements/test_Triangle.cpp

@@ -14,13 +14,22 @@ public:
     Triangle<1> tri;
 };
 
-TEST_F(MasterTriangleOrder1, determinant) {
-    DiagonalMatrixGroup dmg{1,1};
+class MasterTriangleOrder2 : public ::testing::Test {
+public:
+    virtual void SetUp() {
+        nodes.push_back(XY{0.0, 0.0});
+        nodes.push_back(XY{0.5, 0.0});
+        nodes.push_back(XY{1.0, 0.0});
+        nodes.push_back(XY{0.5, 0.5});
+        nodes.push_back(XY{0.0, 1.0});
+        nodes.push_back(XY{0.0, 0.5});
 
-    tri.determinant<1>(dmg,nodes);
+        tri = Triangle<2>(0, {0,1,2,3,4,5});
+    }
 
-    EXPECT_NEAR(dmg(0)(0), 1.0, FLT_EPSILON);
-}
+    std::vector<XY> nodes;
+    Triangle<2> tri;
+};
 
 TEST_F(MasterTriangleOrder1, basis) {
     SparseMatrixGroup smg{1,3,1,3};
@@ -43,4 +52,123 @@ TEST_F(MasterTriangleOrder1, derivative) {
         EXPECT_NEAR(dmg.dx(0,i,0), dx[i], FLT_EPSILON);
         EXPECT_NEAR(dmg.dy(0,i,0), dy[i], FLT_EPSILON);
     }
+}
+
+TEST_F(MasterTriangleOrder1, determinant) {
+    DiagonalMatrixGroup dmg{1,1};
+
+    tri.determinant<1>(dmg,nodes);
+
+    EXPECT_NEAR(dmg(0)(0), 1.0, FLT_EPSILON);
+}
+
+TEST_F(MasterTriangleOrder2, basis) {
+    SparseMatrixGroup smg{3,6,1,6};
+
+    tri.basis<2>(smg, nodes);
+
+    std::vector<double_t> vals{-1.0,+4.0,+2.0,+4.0,-1.0,+1.0};
+    for (double_t &v : vals) {
+        v /= 9.0;
+    }
+
+    for (size_t q = 0; q != TriangleQuadrature<2>::size; ++q) {
+        for (size_t i = 0; i != 6; ++i) {
+            EXPECT_NEAR(smg(q,i,0), vals[i], FLT_EPSILON);
+        }
+        std::rotate(vals.begin(), vals.end() - 2, vals.end());
+    }
+}
+
+TEST_F(MasterTriangleOrder2, jacobian) {
+
+    auto J = tri.jacobian<1>(nodes, 0);
+
+    EXPECT_NEAR(J(0,0),1.0, FLT_EPSILON);
+    EXPECT_NEAR(J(0,1),0.0, FLT_EPSILON);
+    EXPECT_NEAR(J(1,0),0.0, FLT_EPSILON);
+    EXPECT_NEAR(J(1,1),1.0, FLT_EPSILON);
+
+    for (size_t q = 0; q != TriangleQuadrature<2>::size; q++) {
+        auto J = tri.jacobian<2>(nodes, q);
+
+        EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
+        EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
+        EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
+        EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);
+    }
+}
+
+TEST_F(MasterTriangleOrder2, determinant) {
+    DiagonalMatrixGroup dmg{3,1};
+
+    tri.determinant<2>(dmg, nodes);
+
+    for (size_t q = 0; q != 3; ++q) {
+        EXPECT_NEAR(dmg(q)(0), 1.0, FLT_EPSILON);
+    }
+
+    nodes[3].x(M_SQRT1_2);
+    nodes[3].y(M_SQRT1_2);
+
+    tri.determinant<2>(dmg, nodes);
+
+    EXPECT_NEAR(dmg(0)(0), dmg(1)(0), FLT_EPSILON);
+    EXPECT_NEAR((dmg(0)(0) + dmg(1)(0) + dmg(1)(0)) / 6.0, M_PI_4, 0.1 * M_PI_4);
+}
+
+TEST_F(MasterTriangleOrder2, derivative) {
+    DerivativeMatrixGroup dmg{3,6,1,6};
+
+    tri.derivative<2>(dmg, nodes);
+
+    for (size_t q = 0; q != TriangleQuadrature<2>::size; ++q) {
+        double_t const &a = TriangleQuadrature<2>::a[q];
+        double_t const &b = TriangleQuadrature<2>::b[q];
+
+        double_t dx;
+        double_t dy;
+
+        // l0
+        dx = 4 * a + 4 * b - 3;
+        dy = 4 * a + 4 * b - 3;
+
+        EXPECT_NEAR(dmg.dx(q,0,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,0,0), dy, FLT_EPSILON);
+
+        // l1
+        dx = 4 - 4 * b - 8 * a;
+        dy = -4 * a;
+
+        EXPECT_NEAR(dmg.dx(q,1,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,1,0), dy, FLT_EPSILON);
+
+        // l2
+        dx = 4 * a - 1;
+        dy = 0.0;
+
+        EXPECT_NEAR(dmg.dx(q,2,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,2,0), dy, FLT_EPSILON);
+
+        // l3
+        dx = 4 * b;
+        dy = 4 * a;
+
+        EXPECT_NEAR(dmg.dx(q,3,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,3,0), dy, FLT_EPSILON);
+
+        // l4
+        dx = 0;
+        dy = 4 * b - 1;
+
+        EXPECT_NEAR(dmg.dx(q,4,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,4,0), dy, FLT_EPSILON);
+
+        // l5
+        dx = - 4 * b;
+        dy = 4 - 8 * b - 4 * a;
+
+        EXPECT_NEAR(dmg.dx(q,5,0), dx, FLT_EPSILON);
+        EXPECT_NEAR(dmg.dy(q,5,0), dy, FLT_EPSILON);
+    }
 }
\ No newline at end of file

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -137,49 +137,49 @@ public:
 };
 
 TEST_F(MasterTriangleRotationReflection, jacobian_1) {
-    auto J = tri[0].jacobian<1>(node);
+    auto J = tri[0].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);
 
-    J = tri[1].jacobian<1>(node);
+    J = tri[1].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), 1.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
 
-    J = tri[2].jacobian<1>(node);
+    J = tri[2].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);
 
-    J = tri[3].jacobian<1>(node);
+    J = tri[3].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), 1.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), -1.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
 
-    J = tri[4].jacobian<1>(node);
+    J = tri[4].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);
 
-    J = tri[5].jacobian<1>(node);
+    J = tri[5].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);
 
-    J = tri[6].jacobian<1>(node);
+    J = tri[6].jacobian<1>(node, 0);
 
     EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
     EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
@@ -462,7 +462,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     std::cout << r << std::endl;
     std::cout << f << std::endl;
 
-    std::cout << ms.derivatives().dx(0).transpose() * v << std::endl;
+    std::cout << ms.derivatives().da(0).transpose() * v << std::endl;
     std::cout << ms.derivatives().dy(0).transpose() * v << std::endl;
     */
 }