Commit 4d29737e authored by JasonPries's avatar JasonPries
Browse files

Introduce SparseMatrixGroup class for handling numerical quadrature

parent e7e06b95
......@@ -7,6 +7,22 @@
#include "QuadratureRule.h"
#include "Node.h"
template<size_t Q>
class SparseMatrixGroup {
public:
SparseMatrixGroup(size_t const rows, size_t const cols, Eigen::VectorXi const &&nnz) {
for (size_t i = 0; i != Q; ++i) {
Mat[i].resize(rows, cols);
Mat[i].reserve(nnz);
}
};
auto &operator[](size_t const i) { return Mat[i]; };
protected:
std::array<Eigen::SparseMatrix<double>, Q> Mat;
};
// TODO: Curved elements
template<size_t P>
......@@ -25,10 +41,10 @@ public:
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>
void basis(std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, TriangleQuadratureRule<Q>::size> &sp_arr, size_t &row);
void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, size_t const &tri) const;
template<size_t Q>
std::array<std::array<XY, N>, TriangleQuadratureRule<Q>::size> derivative(std::vector<XY> const &nodes) const; // TODO: Rewrite ala ::basis(...);
void derivative(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dx, SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dy, std::vector<XY> const &nodes, size_t const &tri) const;
protected:
size_t Node[N];
......@@ -72,32 +88,29 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &node
template<>
template<size_t Q>
void Triangle<1>::basis(std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, TriangleQuadratureRule<Q>::size> &sp_arr, size_t &row) {
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, size_t const &tri) const {
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
sp_arr[i].coeffRef(row, Node[0]) += TriangleQuadratureRule<Q>::a[i];
sp_arr[i].coeffRef(row, Node[1]) += TriangleQuadratureRule<Q>::b[i];
sp_arr[i].coeffRef(row, Node[2]) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
mat[i].coeffRef(Node[0], tri) += TriangleQuadratureRule<Q>::a[i];
mat[i].coeffRef(Node[1], tri) += TriangleQuadratureRule<Q>::b[i];
mat[i].coeffRef(Node[2], tri) += 1.0 - TriangleQuadratureRule<Q>::a[i] - TriangleQuadratureRule<Q>::b[i];
}
}
template<>
template<size_t Q>
std::array<std::array<XY, Triangle<1>::N>, TriangleQuadratureRule<Q>::size> Triangle<1>::derivative(std::vector<XY> const &nodes) const {
std::array<std::array<XY, Triangle<1>::N>, TriangleQuadratureRule<Q>::size> value{};
void Triangle<1>::derivative(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dx, SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dy, std::vector<XY> const &nodes, size_t const &tri) const {
auto J = jacobian<1>(nodes);
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++i) {
value[i][0].x(J(0, 0));
value[i][0].y(J(1, 0));
dx[i].coeffRef(Node[0], tri) += J(0, 0);
dy[i].coeffRef(Node[0], tri) += J(1, 0);
value[i][1].x(J(0, 1));
value[i][1].y(J(1, 1));
dx[i].coeffRef(Node[1], tri) += J(0, 1);
dy[i].coeffRef(Node[1], tri) += J(1, 1);
value[i][2].x(-J(0, 0) - J(0, 1));
value[i][2].y(-J(1, 0) - J(1, 1));
dx[i].coeffRef(Node[2], tri) += (-J(0, 0) - J(0, 1));
dy[i].coeffRef(Node[2], tri) += (-J(1, 0) - J(1, 1));
}
return value;
};
#endif //OERSTED_TRIANGLE_H
\ No newline at end of file
......@@ -78,63 +78,71 @@ TEST_F(Master_Triangle_1, jacobian_1) {
}
TEST_F(Master_Triangle_1, basis_1) {
std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, 1> sp_arr;
sp_arr[0].resize(tri.size(), node.size());
sp_arr[0].reserve(Eigen::VectorXi::Constant(tri.size(),6));
for (size_t i = 0;i != tri.size();++i) {
tri[i].basis<1>(sp_arr, i);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(0)), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(1)), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(2)), 1.0 / 3.0, FLT_EPSILON);
SparseMatrixGroup<1> mat(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
for (size_t i = 0; i != tri.size(); ++i) {
tri[i].basis<1>(mat, i);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(0),i), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(1),i), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(2),i), 1.0 / 3.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_1) {
auto d = tri[0].derivative<1>(node);
EXPECT_NEAR(d[0][0].x(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][0].y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[0][1].x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[0][1].y(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][2].x(), -1.0, FLT_EPSILON);
EXPECT_NEAR(d[0][2].y(), -1.0, FLT_EPSILON);
SparseMatrixGroup<1> dx(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
SparseMatrixGroup<1> dy(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
for(size_t i = 0;i!=tri.size();++i) {
tri[i].derivative<1>(dx, dy, node, i);
}
// TODO: Test all elements
EXPECT_NEAR(dx[0].coeffRef(tri[0].node(0), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(dy[0].coeffRef(tri[0].node(0), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(dx[0].coeffRef(tri[0].node(1), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(dy[0].coeffRef(tri[0].node(1), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(dx[0].coeffRef(tri[0].node(2), 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(dy[0].coeffRef(tri[0].node(2), 0), -1.0, FLT_EPSILON);
}
TEST_F(Master_Triangle_1, basis_2) {
std::array<Eigen::SparseMatrix<double, Eigen::RowMajor>, 3> sp_arr;
for (size_t i = 0;i!=3;++i) {
sp_arr[i].resize(tri.size(), node.size());
sp_arr[i].reserve(Eigen::VectorXi::Constant(tri.size(),6));
}
SparseMatrixGroup<3> mat(node.size(),tri.size(),Eigen::VectorXi::Constant(tri.size(),6));
for (size_t i = 0; i != tri.size(); ++i) {
tri[i].basis<2>(sp_arr, i);
tri[i].basis<2>(mat, i);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(0)), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(1)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[0].coeffRef(i,tri[i].node(2)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(0),i), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(1),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[0].coeffRef(tri[i].node(2),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(0)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(1)), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[1].coeffRef(i,tri[i].node(2)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[1].coeffRef(tri[i].node(0),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[1].coeffRef(tri[i].node(1),i), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat[1].coeffRef(tri[i].node(2),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(0)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(1)), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(sp_arr[2].coeffRef(i,tri[i].node(2)), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat[2].coeffRef(tri[i].node(0),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[2].coeffRef(tri[i].node(1),i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat[2].coeffRef(tri[i].node(2),i), 2.0 / 3.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_2) {
auto d = tri[0].derivative<2>(node);
SparseMatrixGroup<3> dx(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
SparseMatrixGroup<3> dy(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
for (size_t i = 0; i != 3; ++i) {
EXPECT_NEAR(d[i][0].x(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][0].y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[i][1].x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(d[i][1].y(), 1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][2].x(), -1.0, FLT_EPSILON);
EXPECT_NEAR(d[i][2].y(), -1.0, FLT_EPSILON);
for(size_t i =0;i!=tri.size();++i) {
tri[i].derivative<2>(dx, dy, node, i);
}
std::cout << dx[0] << std::endl;
std::cout << dy[0] << std::endl;
// TODO: Test all elements
for (size_t i = 0; i != 3; ++i) {
EXPECT_NEAR(dx[i].coeffRef(tri[0].node(0),0), 1.0, FLT_EPSILON);
EXPECT_NEAR(dy[i].coeffRef(tri[0].node(0),0), 0.0, FLT_EPSILON);
EXPECT_NEAR(dx[i].coeffRef(tri[0].node(1),0), 0.0, FLT_EPSILON);
EXPECT_NEAR(dy[i].coeffRef(tri[0].node(1),0), 1.0, FLT_EPSILON);
EXPECT_NEAR(dx[i].coeffRef(tri[0].node(2),0), -1.0, FLT_EPSILON);
EXPECT_NEAR(dy[i].coeffRef(tri[0].node(2),0), -1.0, FLT_EPSILON);
}
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment