Triangle.h 4.73 KB
Newer Older
1 2 3
#ifndef OERSTED_TRIANGLE_H
#define OERSTED_TRIANGLE_H

JasonPries's avatar
JasonPries committed
4 5
#include <vector>

6 7 8
#include "Eigen"
#include "Eigen/Sparse"

9 10
#include "Quadrature.hpp"

11
#include "MatrixGroup.h"
12 13
#include "Node.h"

14
// TODO: Curved elements
JasonPries's avatar
JasonPries committed
15

16 17 18 19 20
template<size_t Dimension>
class Element {
public:
    Element() : ID{0} {};
    Element(size_t id) : ID{id} {};
JasonPries's avatar
JasonPries committed
21

JasonPries's avatar
JasonPries committed
22 23
    size_t const id() const { return ID; };

JasonPries's avatar
JasonPries committed
24
protected:
25
    size_t ID;
JasonPries's avatar
JasonPries committed
26 27
};

28
class Facet : public Element<2> {
JasonPries's avatar
JasonPries committed
29
public:
30 31
    Facet() : Element{} {};
    Facet(size_t id) : Element{id} {};
32 33
};

34
template<size_t P>
35
class Triangle : public Facet {
36
public:
37
    static constexpr size_t NumNodes{(P + 1) * (P + 2) / 2};
38

39
    Triangle() : Facet{}, Node{} {};
40

41 42 43 44 45
    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];
        }
    };
46

47 48
    void node(size_t const &i, size_t n) { Node[i] = n; };

49
    size_t const &node(size_t const &i) const { return Node[i]; };
50 51 52 53 54

    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>
55
    void determinant(DiagonalMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const;
56 57

    template<size_t Q>
58
    void basis(SparseMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const;
59 60

    template<size_t Q>
61
    void derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size> &df, std::vector<XY> const &nodes) const;
62

63 64 65
    template<size_t Q>
    std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const;

66
protected:
67
    size_t Node[NumNodes]; // TODO: std::array instead?
68 69
};

70
/*
71 72
template<>
template<>
73
inline Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<0>(std::vector<XY> const &nodes) const {
74 75 76 77 78 79 80 81 82
    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;
}
83
*/
84 85 86

template<>
template<>
87
inline Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
    Eigen::Matrix<double, 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;

    value(0, 0) = yy / det;
    value(0, 1) = -xy / det;
    value(1, 0) = -yx / det;
    value(1, 1) = xx / det;

    return value;
};

template<>
template<size_t Q>
110 111
void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadrature<Q>::size> &mat, std::vector<XY> const &nodes) const {
    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
112 113 114 115 116 117 118 119
        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();
JasonPries's avatar
JasonPries committed
120
        mat(i)(ID) = xx * yy - xy * yx;
121 122 123
    }
}

124

125 126
template<>
template<size_t Q>
127 128 129 130 131
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadrature<Q>::size> &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];
132 133 134 135 136
    }
}

template<>
template<size_t Q>
137
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size> &df, std::vector<XY> const &nodes) const {
138 139
    auto J = jacobian<1>(nodes);

140
    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
141 142
        df.dx(i, Node[0], ID) += J(0, 0);
        df.dy(i, Node[0], ID) += J(1, 0);
143

144 145
        df.dx(i, Node[1], ID) += J(0, 1);
        df.dy(i, Node[1], ID) += J(1, 1);
146

147 148
        df.dx(i, Node[2], ID) += (-J(0, 0) - J(0, 1));
        df.dy(i, Node[2], ID) += (-J(1, 0) - J(1, 1));
149 150 151
    }
};

152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
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);
    }

    return qp;
}

171
#endif //OERSTED_TRIANGLE_H