Triangle.h 3.81 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
#include "MatrixGroup.h"
10
#include "Node.h"
11
#include "QuadratureRule.h"
12

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

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

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

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

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

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

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

JasonPries's avatar
JasonPries committed
40
    size_t const &node(size_t const &i) const { return Node[i]; };
41
42

    template<typename... Args>
43
    Triangle(size_t id, Args &&... args) : Facet{id}, Node{std::forward<Args>(args)...} {};
44
45
46
47
48

    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>
49
    void determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes) const;
50
51

    template<size_t Q>
52
53
54
55
    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;
56
57

protected:
58
    size_t Node[NumNodes];
59
60
};

61
/*
62
63
64
65
66
67
68
69
70
71
72
73
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;
}
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

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>
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
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 {
118
    for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
119
120
121
        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];
122
123
124
125
126
    }
}

template<>
template<size_t Q>
127
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes) const {
128
129
130
    auto J = jacobian<1>(nodes);

    for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
131
132
        df.dx(i, Node[0], ID) += J(0, 0);
        df.dy(i, Node[0], ID) += J(1, 0);
133

134
135
        df.dx(i, Node[1], ID) += J(0, 1);
        df.dy(i, Node[1], ID) += J(1, 1);
136

137
138
        df.dx(i, Node[2], ID) += (-J(0, 0) - J(0, 1));
        df.dy(i, Node[2], ID) += (-J(1, 0) - J(1, 1));
139
140
141
142
    }
};

#endif //OERSTED_TRIANGLE_H