#ifndef OERSTED_TRIANGLE_H #define OERSTED_TRIANGLE_H #include #include "Eigen" #include "Eigen/Sparse" #include "Quadrature.hpp" #include "MatrixGroup.h" #include "Node.h" // TODO: Curved elements template class Element { public: Element() : ID{0} {}; Element(size_t id) : ID{id} {}; size_t const id() const { return ID; }; protected: size_t ID; }; class Facet : public Element<2> { public: Facet() : Element{} {}; Facet(size_t id) : Element{id} {}; }; template class Triangle : public Facet { public: static constexpr size_t NumNodes{(P + 1) * (P + 2) / 2}; Triangle() : Facet{}, Node{} {}; Triangle(size_t id, std::array const &n) : Facet{id} { for (size_t i = 0; i != NumNodes; ++i) { Node[i] = n[i]; } }; size_t const &node(size_t const &i) const { return Node[i]; }; template Eigen::Matrix jacobian(std::vector const &nodes) const; // TODO: Accept Eigen::Matrix as input template void determinant(DiagonalMatrixGroup::size> &mat, std::vector const &nodes) const; template void basis(SparseMatrixGroup::size> &mat, std::vector const &nodes) const; template void derivative(DerivativeMatrixGroup::size> &df, std::vector const &nodes) const; protected: size_t Node[NumNodes]; }; /* template<> template<> inline 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<> inline 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>::determinant(DiagonalMatrixGroup::size> &mat, std::vector const &nodes) const { for (size_t i = 0; i != TriangleQuadrature::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)(ID) = xx * yy - xy * yx; } } template<> template void Triangle<1>::basis(SparseMatrixGroup::size> &mat, std::vector const &nodes) const { for (size_t i = 0; i != TriangleQuadrature::size; ++i) { mat(i, Node[0], ID) += TriangleQuadrature::a[i]; mat(i, Node[1], ID) += TriangleQuadrature::b[i]; mat(i, Node[2], ID) += 1.0 - TriangleQuadrature::a[i] - TriangleQuadrature::b[i]; } } template<> template void Triangle<1>::derivative(DerivativeMatrixGroup::size> &df, std::vector const &nodes) const { auto J = jacobian<1>(nodes); for (size_t i = 0; i != TriangleQuadrature::size; ++i) { df.dx(i, Node[0], ID) += J(0, 0); df.dy(i, Node[0], ID) += J(1, 0); df.dx(i, Node[1], ID) += J(0, 1); df.dy(i, Node[1], ID) += J(1, 1); df.dx(i, Node[2], ID) += (-J(0, 0) - J(0, 1)); df.dy(i, Node[2], ID) += (-J(1, 0) - J(1, 1)); } }; #endif //OERSTED_TRIANGLE_H