Triangle.h 3.94 KB
Newer Older
1
2
3
4
5
6
7
8
9
#ifndef OERSTED_TRIANGLE_H
#define OERSTED_TRIANGLE_H

#include "Eigen"
#include "Eigen/Sparse"

#include "QuadratureRule.h"
#include "Node.h"

10
11
12
template<size_t Q>
class SparseMatrixGroup {
public:
JasonPries's avatar
JasonPries committed
13
    SparseMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) {
14
        for (size_t i = 0; i != Q; ++i) {
JasonPries's avatar
JasonPries committed
15
16
            Matrices[i].resize(rows, cols);
            Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
17
18
19
        }
    };

JasonPries's avatar
JasonPries committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
    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 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); };
39
40

protected:
JasonPries's avatar
JasonPries committed
41
42
    SparseMatrixGroup<Q> Dx;
    SparseMatrixGroup<Q> Dy;
43
44
};

45
46
47
48
49
50
51
52
53
// TODO: Curved elements

template<size_t P>
class Triangle {
public:
    static constexpr size_t N{(P + 1) * (P + 2) / 2};

    Triangle() : Node{} {};

JasonPries's avatar
JasonPries committed
54
    size_t const &node(size_t const &i) const { return Node[i]; };
55
56
57
58
59
60
61
62

    template<typename... Args>
    Triangle(Args &&... args) : 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>
JasonPries's avatar
JasonPries committed
63
    void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const;
64
65

    template<size_t Q>
JasonPries's avatar
JasonPries committed
66
    void derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes, size_t const &tri) const;
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

protected:
    size_t Node[N];
};

template<>
template<>
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<>
Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &nodes) const {
    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>
JasonPries's avatar
JasonPries committed
110
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const {
111
    for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
JasonPries's avatar
JasonPries committed
112
113
114
        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];
115
116
117
118
119
    }
}

template<>
template<size_t Q>
JasonPries's avatar
JasonPries committed
120
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes, size_t const &tri) const {
121
122
123
    auto J = jacobian<1>(nodes);

    for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
JasonPries's avatar
JasonPries committed
124
125
        df.dx(i, Node[0], tri) += J(0, 0);
        df.dy(i, Node[0], tri) += J(1, 0);
126

JasonPries's avatar
JasonPries committed
127
128
        df.dx(i, Node[1], tri) += J(0, 1);
        df.dy(i, Node[1], tri) += J(1, 1);
129

JasonPries's avatar
JasonPries committed
130
131
        df.dx(i, Node[2], tri) += (-J(0, 0) - J(0, 1));
        df.dy(i, Node[2], tri) += (-J(1, 0) - J(1, 1));
132
133
134
135
    }
};

#endif //OERSTED_TRIANGLE_H