Commit b3c33958 authored by JasonPries's avatar JasonPries
Browse files

Introduce DerivativeMatrix Group.

Refactor tests and matrix construction to utilize interfaces for XMatrixGroup classes
parent 4d29737e
......@@ -38,18 +38,42 @@ private:
XYZ d[Q];
};
template<size_t M, size_t N>
class FiniteElementMesh {
};
template<size_t D, size_t P>
class FiniteElementMesh {};
template<size_t P>
class FiniteElementMesh<2, P> {
class FiniteElementMesh<2,P> {
public:
FiniteElementMesh() {};
FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<P>> tris) : Nodes(nodes), Triangles(tris) {};
std::vector<XY> const nodes() { return Nodes; };
std::vector<XY> const nodes() const { return Nodes; };
std::vector<Triangle<P>> const triangles() const { return Triangles; };
template<size_t Q>
SparseMatrixGroup<TriangleQuadratureRule<Q>::size> basis() const {
SparseMatrixGroup<TriangleQuadratureRule<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<P>::N);
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].basis<Q>(mat, Nodes, i);
}
return mat;
};
template<size_t Q>
DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> derivative() const {
DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<P>::N);
for (size_t i = 0; i!= Triangles.size();++i) {
Triangles[i].derivative<Q>(df, Nodes, i);
}
std::vector<Triangle<P>> const triangles() { return Triangles; };
return df;
}
protected:
std::vector<XY> Nodes;
......
......@@ -10,17 +10,36 @@
template<size_t Q>
class SparseMatrixGroup {
public:
SparseMatrixGroup(size_t const rows, size_t const cols, Eigen::VectorXi const &&nnz) {
SparseMatrixGroup(size_t const rows, size_t const cols, size_t const nnz) {
for (size_t i = 0; i != Q; ++i) {
Mat[i].resize(rows, cols);
Mat[i].reserve(nnz);
Matrices[i].resize(rows, cols);
Matrices[i].reserve(Eigen::VectorXi::Constant(cols, nnz));
}
};
auto &operator[](size_t const i) { return Mat[i]; };
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); };
protected:
std::array<Eigen::SparseMatrix<double>, Q> Mat;
SparseMatrixGroup<Q> Dx;
SparseMatrixGroup<Q> Dy;
};
// TODO: Curved elements
......@@ -32,7 +51,7 @@ public:
Triangle() : Node{} {};
size_t const &node(size_t const &i) { return Node[i]; };
size_t const &node(size_t const &i) const { return Node[i]; };
template<typename... Args>
Triangle(Args &&... args) : Node{std::forward<Args>(args)...} {};
......@@ -41,10 +60,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(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, size_t const &tri) const;
void basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const;
template<size_t Q>
void derivative(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dx, SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dy, std::vector<XY> const &nodes, size_t const &tri) const;
void derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, std::vector<XY> const &nodes, size_t const &tri) const;
protected:
size_t Node[N];
......@@ -88,28 +107,28 @@ Eigen::Matrix<double, 2, 2> Triangle<1>::jacobian<1>(std::vector<XY> const &node
template<>
template<size_t Q>
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, size_t const &tri) const {
void Triangle<1>::basis(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &mat, std::vector<XY> const &nodes, size_t const &tri) const {
for (size_t i = 0; i != TriangleQuadratureRule<Q>::size; ++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];
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];
}
}
template<>
template<size_t Q>
void Triangle<1>::derivative(SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dx, SparseMatrixGroup<TriangleQuadratureRule<Q>::size> &dy, std::vector<XY> const &nodes, size_t const &tri) const {
void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadratureRule<Q>::size> &df, 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) {
dx[i].coeffRef(Node[0], tri) += J(0, 0);
dy[i].coeffRef(Node[0], tri) += J(1, 0);
df.dx(i, Node[0], tri) += J(0, 0);
df.dy(i, Node[0], tri) += J(1, 0);
dx[i].coeffRef(Node[1], tri) += J(0, 1);
dy[i].coeffRef(Node[1], tri) += J(1, 1);
df.dx(i, Node[1], tri) += J(0, 1);
df.dy(i, Node[1], tri) += 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));
df.dx(i, Node[2], tri) += (-J(0, 0) - J(0, 1));
df.dy(i, Node[2], tri) += (-J(1, 0) - J(1, 1));
}
};
......
......@@ -4,7 +4,7 @@ inline constexpr size_t operator "" _zu(unsigned long long int i) {
return (size_t) i;
}
class Master_Triangle_1 : public ::testing::Test {
class Master_Triangle : public ::testing::Test {
public:
virtual void SetUp() {
node.push_back(XY{0.0, 0.0});
......@@ -20,13 +20,17 @@ public:
tri.push_back(Triangle<1>(3_zu, 2_zu, 0_zu)); // y-axis reflection
tri.push_back(Triangle<1>(1_zu, 4_zu, 0_zu)); // x-axis reflection
tri.push_back(Triangle<1>(4_zu, 3_zu, 0_zu)); // x=y reflection
femesh = FiniteElementMesh<2, 1>(node, tri);
}
std::vector<XY> node;
std::vector<Triangle<1>> tri;
FiniteElementMesh<2, 1> femesh;
};
TEST_F(Master_Triangle_1, jacobian_1) {
TEST_F(Master_Triangle, jacobian_1) {
auto J = tri[0].jacobian<1>(node);
EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
......@@ -77,72 +81,162 @@ TEST_F(Master_Triangle_1, jacobian_1) {
EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
}
TEST_F(Master_Triangle_1, basis_1) {
SparseMatrixGroup<1> mat(node.size(), tri.size(), Eigen::VectorXi::Constant(tri.size(), 6));
TEST_F(Master_Triangle, basis_1) {
auto mat = femesh.basis<1>();
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);
EXPECT_NEAR(mat(0, tri[i].node(0), i), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat(0, tri[i].node(2), i), 1.0 / 3.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_1) {
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));
TEST_F(Master_Triangle, derivative_1) {
auto df = femesh.derivative<1>();
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);
EXPECT_NEAR(df.dx(0, tri[0].node(0), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(0, tri[0].node(0), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(df.dx(0, tri[0].node(1), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(0, tri[0].node(1), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(df.dx(0, tri[0].node(2), 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(0, tri[0].node(2), 0), -1.0, FLT_EPSILON);
}
TEST_F(Master_Triangle_1, basis_2) {
SparseMatrixGroup<3> mat(node.size(),tri.size(),Eigen::VectorXi::Constant(tri.size(),6));
TEST_F(Master_Triangle, basis_2) {
auto mat = femesh.basis<2>();
for (size_t i = 0; i != tri.size(); ++i) {
tri[i].basis<2>(mat, i);
EXPECT_NEAR(mat(0, tri[i].node(0), i), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat(1, tri[i].node(0), i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat(1, tri[i].node(1), i), 2.0 / 3.0, FLT_EPSILON);
EXPECT_NEAR(mat(1, tri[i].node(2), i), 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(mat(2, tri[i].node(0), i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat(2, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);
EXPECT_NEAR(mat(2, tri[i].node(2), i), 2.0 / 3.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);
TEST_F(Master_Triangle, derivative_2) {
auto df = femesh.derivative<2>();
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);
for (size_t q = 0; q != 3; ++q) {
EXPECT_NEAR(df.dx(q, tri[0].node(0), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(q, tri[0].node(0), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(df.dx(q, tri[0].node(1), 0), 0.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(q, tri[0].node(1), 0), 1.0, FLT_EPSILON);
EXPECT_NEAR(df.dx(q, tri[0].node(2), 0), -1.0, FLT_EPSILON);
EXPECT_NEAR(df.dy(q, tri[0].node(2), 0), -1.0, FLT_EPSILON);
}
}
TEST_F(Master_Triangle_1, derivative_2) {
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));
class Simple_Square : public ::testing::Test {
public:
virtual void SetUp() {
node.push_back(XY{0.0, 0.0});
node.push_back(XY{1.0, 0.0});
node.push_back(XY{0.0, 1.0});
node.push_back(XY{1.0, 1.0});
tri.push_back(Triangle<1>(0_zu, 1_zu, 2_zu));
tri.push_back(Triangle<1>(1_zu, 3_zu, 2_zu));
for(size_t i =0;i!=tri.size();++i) {
tri[i].derivative<2>(dx, dy, node, i);
femesh = FiniteElementMesh<2, 1>(node, tri);
vx_m << 1.0, 0.0, 1.0, 0.0;
vx_p << 0.0, 1.0, 0.0, 1.0;
vy_m << 1.0, 1.0, 0.0, 0.0;
vy_p << 0.0, 0.0, 1.0, 1.0;
}
std::cout << dx[0] << std::endl;
std::cout << dy[0] << std::endl;
std::vector<XY> node;
std::vector<Triangle<1>> tri;
Eigen::Vector4d vx_m;
Eigen::Vector4d vx_p;
Eigen::Vector4d vy_m;
Eigen::Vector4d vy_p;
FiniteElementMesh<2, 1> femesh;
};
TEST_F(Simple_Square, derivative_1_temp) {
auto df = femesh.derivative<1>();
Eigen::Vector2d tmp;
tmp = df.dx(0).transpose() * vx_p;
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
tmp = df.dx(0).transpose() * vx_m;
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
tmp = df.dx(0).transpose() * vy_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dx(0).transpose() * vy_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dy(0).transpose() * vy_p;
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
tmp = df.dy(0).transpose() * vy_m;
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
tmp = df.dy(0).transpose() * vx_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dy(0).transpose() * vx_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
}
TEST_F(Simple_Square, derivative_2_temp) {
auto df = femesh.derivative<2>();
Eigen::Vector2d tmp;
// 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);
tmp = df.dx(i).transpose() * vx_p;
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
tmp = df.dx(i).transpose() * vx_m;
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
tmp = df.dx(i).transpose() * vy_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dx(i).transpose() * vy_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dy(i).transpose() * vy_p;
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
tmp = df.dy(i).transpose() * vy_m;
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
tmp = df.dy(i).transpose() * vx_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
tmp = df.dy(i).transpose() * vx_m;
EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
EXPECT_NEAR(tmp(0), 0.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