Triangle.h 2.89 KB
 JasonPries committed Dec 27, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 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 ``````#ifndef OERSTED_TRIANGLE_H #define OERSTED_TRIANGLE_H #include "Eigen" #include "Eigen/Sparse" #include "QuadratureRule.h" #include "Node.h" // TODO: Curved elements template class Triangle { public: static constexpr size_t N{(P + 1) * (P + 2) / 2}; Triangle() : Node{} {}; size_t const &node(size_t const &i) { return Node[i]; }; template Triangle(Args &&... args) : Node{std::forward(args)...} {}; template Eigen::Matrix jacobian(std::vector const &nodes) const; // TODO: Accept Eigen::Matrix as input template void basis(std::array, TriangleQuadratureRule::size> &sp_arr, size_t &row); template std::array, TriangleQuadratureRule::size> derivative(std::vector const &nodes) const; // TODO: Rewrite ala ::basis(...); protected: size_t Node[N]; }; template<> template<> Eigen::Matrix Triangle<1>::jacobian<0>(std::vector const &nodes) const { Eigen::Matrix 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 Triangle<1>::jacobian<1>(std::vector const &nodes) const { Eigen::Matrix 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 void Triangle<1>::basis(std::array, TriangleQuadratureRule::size> &sp_arr, size_t &row) { for (size_t i = 0; i != TriangleQuadratureRule::size; ++i) { sp_arr[i].coeffRef(row, Node[0]) += TriangleQuadratureRule::a[i]; sp_arr[i].coeffRef(row, Node[1]) += TriangleQuadratureRule::b[i]; sp_arr[i].coeffRef(row, Node[2]) += 1.0 - TriangleQuadratureRule::a[i] - TriangleQuadratureRule::b[i]; } } template<> template std::array::N>, TriangleQuadratureRule::size> Triangle<1>::derivative(std::vector const &nodes) const { std::array::N>, TriangleQuadratureRule::size> value{}; auto J = jacobian<1>(nodes); for (size_t i = 0; i != TriangleQuadratureRule::size; ++i) { value[i][0].x(J(0, 0)); value[i][0].y(J(1, 0)); value[i][1].x(J(0, 1)); value[i][1].y(J(1, 1)); value[i][2].x(-J(0, 0) - J(0, 1)); value[i][2].y(-J(1, 0) - J(1, 1)); } return value; }; #endif //OERSTED_TRIANGLE_H``````