diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..94e8fbf8562cd7a24ff69279d876b0f95e8153c2
--- /dev/null
+++ b/Tests/NumLib/TestFe.cpp
@@ -0,0 +1,436 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <gtest/gtest.h>
+
+#include <vector>
+#include <cmath>
+#ifdef OGS_USE_EIGEN
+#include <Eigen/Eigen>
+#endif
+
+#include "MeshLib/Elements/Line.h"
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "Tests/TestTools.h"
+
+using namespace NumLib;
+
+namespace
+{
+
+// copy matrix entries in upper triangle to lower triangle
+template <class T_MATRIX, typename ID_TYPE=signed>
+void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m)
+{
+ for (ID_TYPE i=0; i<dim; i++)
+ for (ID_TYPE j=0; j<i; j++)
+ m(i,j) = m(j,i);
+}
+
+// set an identity matrix
+template <class T_MATRIX, typename ID_TYPE=signed>
+void setIdentityMatrix(unsigned dim, T_MATRIX &m)
+{
+ for (unsigned i=0; i<dim; i++)
+ for (unsigned j=0; j<dim; j++)
+ m(i,j) = 0.0;
+ for (unsigned i=0; i<dim; i++)
+ m(i,i) = 1.0;
+}
+
+// Test case for LINE2
+class TestFeLINE2
+{
+public:
+ template <class T_MATRIX_TYPES>
+ using FeType = NumLib::FeLINE2<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType>;
+ typedef MeshLib::Line MeshElementType;
+ static const unsigned dim = MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+
+ MeshLib::Line* createMeshElement()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ // set an expected mass matrix
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1./3.; m(0,1) = 1./6.;
+ m(1,1) = 1./3.;
+ // make symmetric
+ copyUpperToLower(2, m);
+ }
+
+ // set an expected laplace matrix
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = -1.0;
+ m(1,1) = 1.0;
+ // make symmetric
+ copyUpperToLower(2, m);
+ m *= k;
+ }
+};
+
+// Test case for QUAD4
+class TestFeQUAD4
+{
+public:
+ template <class T_MATRIX_TYPES>
+ using FeType = NumLib::FeQUAD4<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType>;
+ typedef MeshLib::Quad MeshElementType;
+ static const unsigned dim = 2; //MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+
+ MeshLib::Quad* createMeshElement()
+ {
+ // square
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // set an expected mass matrix for 1m x 1m
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
+ m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
+ m(2,2) = 1.0; m(2,3) = 1./2;
+ m(3,3) = 1.0;
+ // make symmetric
+ copyUpperToLower(4, m);
+ m *= 1./9.;
+ }
+
+ // set an expected laplace matrix for 1m x 1m
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
+ m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
+ m(2,2) = 4.0; m(2,3) = -1.0;
+ m(3,3) = 4.0;
+ // make symmetric
+ copyUpperToLower(4, m);
+ m *= k/6.;
+ }
+};
+
+// Test case for HEX8
+class TestFeHEX8
+{
+public:
+ template <class T_MATRIX_TYPES>
+ using FeType = NumLib::FeHEX8<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType>;
+ typedef MeshLib::Hex MeshElementType;
+ static const unsigned dim = 3; //MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+
+ MeshLib::Hex* createMeshElement()
+ {
+ // cubic
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+ nodes[4] = new MeshLib::Node(0.0, 0.0, 1.0);
+ nodes[5] = new MeshLib::Node(1.0, 0.0, 1.0);
+ nodes[6] = new MeshLib::Node(1.0, 1.0, 1.0);
+ nodes[7] = new MeshLib::Node(0.0, 1.0, 1.0);
+ return new MeshLib::Hex(nodes);
+ }
+
+ // set an expected mass matrix
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2; m(0,4) = 1./2; m(0,5) = 1./4; m(0,6) = 1./8; m(0,7) = 1./4;
+ m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4; m(1,4) = 1./4; m(1,5) = 1./2; m(1,6) = 1./4; m(1,7) = 1./8;
+ m(2,2) = 1.0; m(2,3) = 1./2; m(2,4) = 1./8; m(2,5) = 1./4; m(2,6) = 1./2; m(2,7) = 1./4;
+ m(3,3) = 1.0; m(3,4) = 1./4; m(3,5) = 1./8; m(3,6) = 1./4; m(3,7) = 1./2;
+ m(4,4) = 1.0; m(4,5) = 1./2; m(4,6) = 1./4; m(4,7) = 1./2;
+ m(5,5) = 1.0; m(5,6) = 1./2; m(5,7) = 1./4;
+ m(6,6) = 1.0; m(6,7) = 1./2;
+ m(7,7) = 1.0;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= 1./27.;
+ }
+
+ // set an expected laplace matrix
+ template <class T_MATRIX, typename ID_TYPE=signed>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 2./3; m(0,1) = 0; m(0,2) = -1./6; m(0,3) = 0; m(0,4) = 0; m(0,5) = -1./6; m(0,6) = -1./6; m(0,7) = -1./6;
+ m(1,1) = 2./3; m(1,2) = 0; m(1,3) = -1./6; m(1,4) = -1./6; m(1,5) = 0; m(1,6) = -1./6; m(1,7) = -1./6;
+ m(2,2) = 2./3; m(2,3) = 0; m(2,4) = -1./6; m(2,5) = -1./6; m(2,6) = 0; m(2,7) = -1./6;
+ m(3,3) = 2./3; m(3,4) = -1./6; m(3,5) = -1./6; m(3,6) = -1./6; m(3,7) = 0;
+ m(4,4) = 2./3; m(4,5) = 0; m(4,6) = -1./6; m(4,7) = 0;
+ m(5,5) = 2./3; m(5,6) = 0; m(5,7) = -1./6;
+ m(6,6) = 2./3; m(6,7) = 0;
+ m(7,7) = 2./3;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= k/2.;
+ }
+};
+
+template <class T_MATRIX_TYPES, class T_FE>
+class NumLibFemIsoTest : public ::testing::Test, public T_FE
+{
+ public:
+ // Matrix types
+ typedef T_MATRIX_TYPES MatrixType;
+ typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix;
+ typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector;
+ typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix;
+ typedef typename T_MATRIX_TYPES::DimMatrixType DimMatrix;
+ // Finite element type
+ typedef typename T_FE::template FeType<T_MATRIX_TYPES>::type FeType;
+ // Shape matrix data type
+ typedef typename FeType::ShapeMatricesType ShapeMatricesType;
+ typedef typename T_FE::MeshElementType MeshElementType;
+
+ static const unsigned dim = T_FE::dim;
+ static const unsigned e_nnodes = T_FE::e_nnodes;
+
+ public:
+ NumLibFemIsoTest() :
+ D(dim, dim),
+ expectedM(e_nnodes,e_nnodes),
+ expectedK(e_nnodes,e_nnodes),
+ integration_method(2)
+ {
+ // create a mesh element used for testing
+ mesh_element = this->createMeshElement();
+
+ // set a conductivity tensor
+ setIdentityMatrix(dim, D);
+ D *= conductivity;
+
+ // set expected matrices
+ this->setExpectedMassMatrix(expectedM);
+ this->setExpectedLaplaceMatrix(conductivity, expectedK);
+
+ // for destructor
+ vec_eles.push_back(mesh_element);
+ for (auto e : vec_eles)
+ for (unsigned i=0; i<e->getNNodes(true); i++)
+ vec_nodes.push_back(e->getNode(i));
+ }
+
+ virtual ~NumLibFemIsoTest()
+ {
+ for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
+ delete *itr;
+ for (auto itr = vec_eles.begin(); itr!=vec_eles.end(); ++itr )
+ delete *itr;
+ }
+
+
+ static const double conductivity;
+ static const double eps;
+ DimMatrix D;
+ NodalMatrix expectedM;
+ NodalMatrix expectedK;
+ typename FeType::IntegrationMethod integration_method;
+
+ std::vector<const MeshLib::Node*> vec_nodes;
+ std::vector<const MeshElementType*> vec_eles;
+ MeshElementType* mesh_element;
+
+}; // NumLibFemIsoLineTest
+
+template <class T_MATRIX_TYPES, class T_FE>
+const double NumLibFemIsoTest<T_MATRIX_TYPES, T_FE>::conductivity = 1e-11;
+
+template <class T_MATRIX_TYPES, class T_FE>
+const double NumLibFemIsoTest<T_MATRIX_TYPES, T_FE>::eps = std::numeric_limits<double>::epsilon();
+
+template <class T_MATRIX_TYPES, class T_FE>
+const unsigned NumLibFemIsoTest<T_MATRIX_TYPES, T_FE>::dim;
+
+template <class T_MATRIX_TYPES, class T_FE>
+const unsigned NumLibFemIsoTest<T_MATRIX_TYPES, T_FE>::e_nnodes;
+
+} // namespace
+
+#ifdef OGS_USE_EIGEN
+
+template <unsigned NNODES, unsigned DIM>
+struct EigenFixedMatrixTypes
+{
+ typedef Eigen::Matrix<double, NNODES, DIM, Eigen::RowMajor> NodalMatrixType;
+ typedef Eigen::Matrix<double, NNODES, 1> NodalVectorType;
+ typedef Eigen::Matrix<double, DIM, NNODES, Eigen::RowMajor> DimNodalMatrixType;
+ typedef Eigen::Matrix<double, DIM, DIM, Eigen::RowMajor> DimMatrixType;
+};
+
+struct EigenDynamicMatrixTypes
+{
+ typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
+ typedef NodalMatrixType DimNodalMatrixType;
+ typedef NodalMatrixType DimMatrixType;
+ typedef Eigen::VectorXd NodalVectorType;
+};
+
+template <class T_FE>
+using NumLibFemIsoTestFeTypes = NumLibFemIsoTest<EigenDynamicMatrixTypes, T_FE>;
+#endif // OGS_USE_EIGEN
+
+template <class T_MATRIX_TYPES>
+using NumLibFemIsoTestMatTypes = NumLibFemIsoTest<T_MATRIX_TYPES, TestFeLINE2>;
+
+template <class T>
+struct MatrixTypes
+{
+ typedef ::testing::Types<
+#ifdef OGS_USE_EIGEN
+ EigenFixedMatrixTypes<T::e_nnodes, T::dim>
+ , EigenDynamicMatrixTypes
+#endif
+ > types;
+};
+
+
+typedef ::testing::Types<
+ TestFeLINE2,
+ TestFeQUAD4,
+ TestFeHEX8
+ > FeTypes;
+
+TYPED_TEST_CASE(NumLibFemIsoTestFeTypes, FeTypes);
+
+TYPED_TEST_CASE(NumLibFemIsoTestMatTypes, MatrixTypes<TestFeLINE2>::types);
+
+TYPED_TEST(NumLibFemIsoTestFeTypes, CheckMassMatrix)
+{
+ // Refer to typedefs in the fixture
+ typedef typename TestFixture::FeType FeType;
+ typedef typename TestFixture::NodalMatrix NodalMatrix;
+ typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+
+ // create a finite element object
+ FeType fe(*this->mesh_element);
+
+ // evaluate a mass matrix M = int{ N^T D N }dA_e
+ NodalMatrix M(this->e_nnodes, this->e_nnodes);
+ ShapeMatricesType shape(this->dim, this->e_nnodes);
+ for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
+ shape.setZero();
+ auto wp = this->integration_method.getWeightedPoint(i);
+ fe.template computeShapeFunctions<ShapeMatrixType::N_J>(wp.getCoords(), shape);
+ M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
+ }
+
+ ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
+}
+
+TYPED_TEST(NumLibFemIsoTestFeTypes, CheckLaplaceMatrix)
+{
+ // Refer to typedefs in the fixture
+ typedef typename TestFixture::FeType FeType;
+ typedef typename TestFixture::NodalMatrix NodalMatrix;
+ typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+
+ // create a finite element object
+ FeType fe(*this->mesh_element);
+
+ // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
+ NodalMatrix K(this->e_nnodes, this->e_nnodes);
+ ShapeMatricesType shape(this->dim, this->e_nnodes);
+ for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
+ shape.setZero();
+ auto wp = this->integration_method.getWeightedPoint(i);
+ fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
+ K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
+ }
+ ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
+}
+
+TYPED_TEST(NumLibFemIsoTestFeTypes, CheckMassLaplaceMatrices)
+{
+ // Refer to typedefs in the fixture
+ typedef typename TestFixture::FeType FeType;
+ typedef typename TestFixture::NodalMatrix NodalMatrix;
+ typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+
+ // create a finite element object
+ FeType fe(*this->mesh_element);
+
+ // evaluate both mass and laplace matrices at once
+ NodalMatrix M(this->e_nnodes, this->e_nnodes);
+ NodalMatrix K(this->e_nnodes, this->e_nnodes);
+ ShapeMatricesType shape(this->dim, this->e_nnodes);
+ for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
+ shape.setZero();
+ auto wp = this->integration_method.getWeightedPoint(i);
+ fe.computeShapeFunctions(wp.getCoords(), shape);
+ M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
+ K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
+ }
+ ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
+ ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
+}
+
+TYPED_TEST(NumLibFemIsoTestFeTypes, CheckGaussIntegrationLevel)
+{
+ // Refer to typedefs in the fixture
+ typedef typename TestFixture::FeType FeType;
+ typedef typename TestFixture::NodalMatrix NodalMatrix;
+ typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
+
+ // create a finite element object with gauss quadrature level 2
+ FeType fe(*this->mesh_element);
+
+ // evaluate a mass matrix
+ NodalMatrix M(this->e_nnodes, this->e_nnodes);
+ ShapeMatricesType shape(this->dim, this->e_nnodes);
+ ASSERT_EQ(std::pow(2u, this->dim), this->integration_method.getNPoints());
+ for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
+ shape.setZero();
+ auto wp = this->integration_method.getWeightedPoint(i);
+ fe.computeShapeFunctions(wp.getCoords(), shape);
+ M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
+ }
+ ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
+
+ // Change gauss quadrature level to 3
+ this->integration_method.setIntegrationOrder(3);
+ M *= .0;
+ ASSERT_EQ(std::pow(3u, this->dim), this->integration_method.getNPoints());
+ for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
+ shape.setZero();
+ auto wp = this->integration_method.getWeightedPoint(i);
+ fe.computeShapeFunctions(wp.getCoords(), shape);
+ M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
+ }
+ ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
+}
+
+
diff --git a/Tests/NumLib/TestFeHex8.cpp b/Tests/NumLib/TestFeHex8.cpp
deleted file mode 100644
index ce46bfccbb109fe9540ec982a5d2b0bc2da7f107..0000000000000000000000000000000000000000
--- a/Tests/NumLib/TestFeHex8.cpp
+++ /dev/null
@@ -1,310 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <gtest/gtest.h>
-
-#include <vector>
-#include <cmath>
-#ifdef OGS_USE_EIGEN
-#include <Eigen/Eigen>
-#endif
-
-#include "MeshLib/Elements/Hex.h"
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
-#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
-
-#include "Tests/TestTools.h"
-
-using namespace NumLib;
-
-namespace
-{
-
-static const unsigned dim = 3;
-static const unsigned e_nnodes = 8;
-
-template <class T_MATRIX_TYPES>
-class NumLibFemIsoHex8Test : public ::testing::Test
-{
- public:
- // Matrix types
- typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix;
- typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector;
- typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix;
- typedef typename T_MATRIX_TYPES::DimMatrixType DimMatrix;
- // Finite element type
- typedef typename NumLib::FeHEX8<NodalVector, DimNodalMatrix, DimMatrix>::type FeHEX8Type;
- // Shape matrix data type
- typedef typename FeHEX8Type::ShapeMatricesType ShapeMatricesType;
-
- public:
- NumLibFemIsoHex8Test() :
- D(dim, dim),
- expectedM(e_nnodes,e_nnodes),
- expectedK(e_nnodes,e_nnodes),
- integration_method(2)
- {
- // create a quad element used for testing
- mesh_element = createCubicHex(1.0);
-
- // set a conductivity tensor
- setIdentityMatrix(dim, D);
- D *= conductivity;
-
- // set expected matrices
- setExpectedMassMatrix(expectedM);
- setExpectedLaplaceMatrix(conductivity, expectedK);
-
- // for destructor
- vec_eles.push_back(mesh_element);
- for (auto e : vec_eles)
- for (unsigned i=0; i<e->getNNodes(true); i++)
- vec_nodes.push_back(e->getNode(i));
- }
-
- ~NumLibFemIsoHex8Test()
- {
- for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
- delete *itr;
- for (auto itr = vec_eles.begin(); itr!=vec_eles.end(); ++itr )
- delete *itr;
- }
-
- // Quad: square shape with a length of 1m
- MeshLib::Hex* createCubicHex(double h)
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node( h, 0.0, 0.0);
- nodes[2] = new MeshLib::Node( h, h, 0.0);
- nodes[3] = new MeshLib::Node(0.0, h, 0.0);
- nodes[4] = new MeshLib::Node(0.0, 0.0, h);
- nodes[5] = new MeshLib::Node( h, 0.0, h);
- nodes[6] = new MeshLib::Node( h, h, h);
- nodes[7] = new MeshLib::Node(0.0, h, h);
- return new MeshLib::Hex(nodes);
- }
-
- // set an identity matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setIdentityMatrix(unsigned dim, T_MATRIX &m) const
- {
- for (unsigned i=0; i<dim; i++)
- for (unsigned j=0; j<dim; j++)
- m(i,j) = 0.0;
- for (unsigned i=0; i<dim; i++)
- m(i,i) = 1.0;
- }
-
- // copy upper triangles to lower triangles in a matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m) const
- {
- for (ID_TYPE i=0; i<dim; i++)
- for (ID_TYPE j=0; j<i; j++)
- m(i,j) = m(j,i);
- }
-
- // set an expected mass matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2; m(0,4) = 1./2; m(0,5) = 1./4; m(0,6) = 1./8; m(0,7) = 1./4;
- m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4; m(1,4) = 1./4; m(1,5) = 1./2; m(1,6) = 1./4; m(1,7) = 1./8;
- m(2,2) = 1.0; m(2,3) = 1./2; m(2,4) = 1./8; m(2,5) = 1./4; m(2,6) = 1./2; m(2,7) = 1./4;
- m(3,3) = 1.0; m(3,4) = 1./4; m(3,5) = 1./8; m(3,6) = 1./4; m(3,7) = 1./2;
- m(4,4) = 1.0; m(4,5) = 1./2; m(4,6) = 1./4; m(4,7) = 1./2;
- m(5,5) = 1.0; m(5,6) = 1./2; m(5,7) = 1./4;
- m(6,6) = 1.0; m(6,7) = 1./2;
- m(7,7) = 1.0;
- // make symmetric
- copyUpperToLower(e_nnodes, m);
- m *= 1./27.;
- }
-
- // set an expected laplace matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 2./3; m(0,1) = 0; m(0,2) = -1./6; m(0,3) = 0; m(0,4) = 0; m(0,5) = -1./6; m(0,6) = -1./6; m(0,7) = -1./6;
- m(1,1) = 2./3; m(1,2) = 0; m(1,3) = -1./6; m(1,4) = -1./6; m(1,5) = 0; m(1,6) = -1./6; m(1,7) = -1./6;
- m(2,2) = 2./3; m(2,3) = 0; m(2,4) = -1./6; m(2,5) = -1./6; m(2,6) = 0; m(2,7) = -1./6;
- m(3,3) = 2./3; m(3,4) = -1./6; m(3,5) = -1./6; m(3,6) = -1./6; m(3,7) = 0;
- m(4,4) = 2./3; m(4,5) = 0; m(4,6) = -1./6; m(4,7) = 0;
- m(5,5) = 2./3; m(5,6) = 0; m(5,7) = -1./6;
- m(6,6) = 2./3; m(6,7) = 0;
- m(7,7) = 2./3;
- // make symmetric
- copyUpperToLower(e_nnodes, m);
- m *= k/2.;
- }
-
- static const double conductivity;
- static const double eps;
- DimMatrix D;
- NodalMatrix expectedM;
- NodalMatrix expectedK;
- typename FeHEX8Type::IntegrationMethod integration_method;
-
- std::vector<const MeshLib::Node*> vec_nodes;
- std::vector<const MeshLib::Hex*> vec_eles;
- MeshLib::Hex* mesh_element;
-
-}; // NumLibFemIsoQuad4Test
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoHex8Test<T_MATRIX_TYPES>::conductivity = 1e-11;
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoHex8Test<T_MATRIX_TYPES>::eps = std::numeric_limits<double>::epsilon();
-
-} // namespace
-
-#ifdef OGS_USE_EIGEN
-
-struct EigenFixedMatrixTypes
-{
- typedef Eigen::Matrix<double, e_nnodes, e_nnodes, Eigen::RowMajor> NodalMatrixType;
- typedef Eigen::Matrix<double, e_nnodes, 1> NodalVectorType;
- typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
- typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrixType;
-};
-
-struct EigenDynamicMatrixTypes
-{
- typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
- typedef NodalMatrixType DimNodalMatrixType;
- typedef NodalMatrixType DimMatrixType;
- typedef Eigen::VectorXd NodalVectorType;
-};
-
-#endif // OGS_USE_EIGEN
-
-typedef ::testing::Types<
-#ifdef OGS_USE_EIGEN
- EigenFixedMatrixTypes
- , EigenDynamicMatrixTypes
-#endif
- > MatrixTypes;
-
-TYPED_TEST_CASE(NumLibFemIsoHex8Test, MatrixTypes);
-
-TYPED_TEST(NumLibFemIsoHex8Test, CheckMassMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeHEX8Type FeHEX8Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeHEX8Type fe(*this->mesh_element);
-
- // evaluate a mass matrix M = int{ N^T D N }dA_e
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::N_J>(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
-
-// std::cout << this->expectedM << "\n M=";
-// std::cout << M;
-
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoHex8Test, CheckLaplaceMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeHEX8Type FeHEX8Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeHEX8Type fe(*this->mesh_element);
-
- // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
-// std::cout << this->expectedK << "\n M=";
-// std::cout << K;
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoHex8Test, CheckMassLaplaceMatrices)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeHEX8Type FeHEX8Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeHEX8Type fe(*this->mesh_element);
-
- // evaluate both mass and laplace matrices at once
- NodalMatrix M(e_nnodes, e_nnodes);
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoHex8Test, CheckGaussIntegrationLevel)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeHEX8Type FeHEX8Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object with gauss quadrature level 2
- FeHEX8Type fe(*this->mesh_element);
-
- // evaluate a mass matrix
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- ASSERT_EQ(8u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-
- // Change gauss quadrature level to 3
- this->integration_method.setIntegrationOrder(3);
- M *= .0;
- ASSERT_EQ(27u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-
diff --git a/Tests/NumLib/TestFeLine2.cpp b/Tests/NumLib/TestFeLine2.cpp
deleted file mode 100644
index 6f5a1312e4459bb4cc73b3da01748f2f73ba0546..0000000000000000000000000000000000000000
--- a/Tests/NumLib/TestFeLine2.cpp
+++ /dev/null
@@ -1,286 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <gtest/gtest.h>
-
-#include <vector>
-#include <cmath>
-#ifdef OGS_USE_EIGEN
-#include <Eigen/Eigen>
-#endif
-
-#include "MeshLib/Elements/Line.h"
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
-#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
-
-#include "Tests/TestTools.h"
-
-using namespace NumLib;
-
-namespace
-{
-
-static const unsigned dim = 1;
-static const unsigned e_nnodes = 2;
-
-template <class T_MATRIX_TYPES>
-class NumLibFemIsoLine2Test : public ::testing::Test
-{
- public:
- // Matrix types
- typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix;
- typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector;
- typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix;
- typedef typename T_MATRIX_TYPES::DimMatrixType DimMatrix;
- // Finite element type
- typedef typename NumLib::FeLINE2<NodalVector, DimNodalMatrix, DimMatrix>::type FeLINE2Type;
- // Shape matrix data type
- typedef typename FeLINE2Type::ShapeMatricesType ShapeMatricesType;
-
- public:
- NumLibFemIsoLine2Test() :
- D(dim, dim),
- expectedM(e_nnodes,e_nnodes),
- expectedK(e_nnodes,e_nnodes),
- integration_method(2)
- {
- // create a quad element used for testing
- unitLine = createLine(1.0);
-
- // set a conductivity tensor
- setIdentityMatrix(dim, D);
- D *= conductivity;
-
- // set expected matrices
- setExpectedMassMatrix(expectedM);
- setExpectedLaplaceMatrix(conductivity, expectedK);
-
- // for destructor
- vec_eles.push_back(unitLine);
- for (auto e : vec_eles)
- for (unsigned i=0; i<e->getNNodes(true); i++)
- vec_nodes.push_back(e->getNode(i));
- }
-
- virtual ~NumLibFemIsoLine2Test()
- {
- for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
- delete *itr;
- for (auto itr = vec_eles.begin(); itr!=vec_eles.end(); ++itr )
- delete *itr;
- }
-
- // Line: a length of 1m
- MeshLib::Line* createLine(double h)
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node( h, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-
- // set an identity matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setIdentityMatrix(unsigned dim, T_MATRIX &m) const
- {
- for (unsigned i=0; i<dim; i++)
- for (unsigned j=0; j<dim; j++)
- m(i,j) = 0.0;
- for (unsigned i=0; i<dim; i++)
- m(i,i) = 1.0;
- }
-
- // copy upper triangles to lower triangles in a matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m) const
- {
- for (ID_TYPE i=0; i<dim; i++)
- for (ID_TYPE j=0; j<i; j++)
- m(i,j) = m(j,i);
- }
-
- // set an expected mass matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1./3.; m(0,1) = 1./6.;
- m(1,1) = 1./3.;
- // make symmetric
- copyUpperToLower(2, m);
- }
-
- // set an expected laplace matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = -1.0;
- m(1,1) = 1.0;
- // make symmetric
- copyUpperToLower(2, m);
- m *= k;
- }
-
- static const double conductivity;
- static const double eps;
- DimMatrix D;
- NodalMatrix expectedM;
- NodalMatrix expectedK;
- typename FeLINE2Type::IntegrationMethod integration_method;
-
- std::vector<const MeshLib::Node*> vec_nodes;
- std::vector<const MeshLib::Line*> vec_eles;
- MeshLib::Line* unitLine;
-
-}; // NumLibFemIsoLineTest
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoLine2Test<T_MATRIX_TYPES>::conductivity = 1e-11;
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoLine2Test<T_MATRIX_TYPES>::eps = std::numeric_limits<double>::epsilon();
-
-} // namespace
-
-#ifdef OGS_USE_EIGEN
-
-struct EigenFixedMatrixTypes
-{
- typedef Eigen::Matrix<double, e_nnodes, e_nnodes, Eigen::RowMajor> NodalMatrixType;
- typedef Eigen::Matrix<double, e_nnodes, 1> NodalVectorType;
- typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
- typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrixType;
-};
-
-struct EigenDynamicMatrixTypes
-{
- typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
- typedef NodalMatrixType DimNodalMatrixType;
- typedef NodalMatrixType DimMatrixType;
- typedef Eigen::VectorXd NodalVectorType;
-};
-
-#endif // OGS_USE_EIGEN
-
-typedef ::testing::Types<
-#ifdef OGS_USE_EIGEN
- EigenFixedMatrixTypes
- , EigenDynamicMatrixTypes
-#endif
- > MatrixTypes;
-
-TYPED_TEST_CASE(NumLibFemIsoLine2Test, MatrixTypes);
-
-TYPED_TEST(NumLibFemIsoLine2Test, CheckMassMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeLINE2Type FeLINE2Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeLINE2Type fe(*this->unitLine);
-
- // evaluate a mass matrix M = int{ N^T D N }dA_e
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::N_J>(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
-
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoLine2Test, CheckLaplaceMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeLINE2Type FeLINE2Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeLINE2Type fe(*this->unitLine);
-
- // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoLine2Test, CheckMassLaplaceMatrices)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeLINE2Type FeLINE2Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeLINE2Type fe(*this->unitLine);
-
- // evaluate both mass and laplace matrices at once
- NodalMatrix M(e_nnodes, e_nnodes);
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoLine2Test, CheckGaussIntegrationLevel)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeLINE2Type FeLINE2Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object with gauss quadrature level 2
- FeLINE2Type fe(*this->unitLine);
-
- // evaluate a mass matrix
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- ASSERT_EQ(2u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-
- // Change gauss quadrature level to 3
- this->integration_method.setIntegrationOrder(3);
- M *= .0;
- ASSERT_EQ(3u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-
diff --git a/Tests/NumLib/TestFeQuad4.cpp b/Tests/NumLib/TestFeQuad4.cpp
deleted file mode 100644
index 8f47b0b7f78cd706252ffb14c71dc1804d1ecb87..0000000000000000000000000000000000000000
--- a/Tests/NumLib/TestFeQuad4.cpp
+++ /dev/null
@@ -1,297 +0,0 @@
-/**
- * \file TestFEM.cpp
- * \author Norihiro Watanabe
- * \date 2012-08-03
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <gtest/gtest.h>
-
-#include <vector>
-#include <cmath>
-#ifdef OGS_USE_EIGEN
-#include <Eigen/Eigen>
-#endif
-
-#include "MeshLib/Elements/Quad.h"
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
-#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
-
-#include "Tests/TestTools.h"
-
-using namespace NumLib;
-
-namespace
-{
-
-static const unsigned dim = 2;
-static const unsigned e_nnodes = 4;
-
-template <class T_MATRIX_TYPES>
-class NumLibFemIsoQuad4Test : public ::testing::Test
-{
- public:
- // Matrix types
- typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix;
- typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector;
- typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix;
- typedef typename T_MATRIX_TYPES::DimMatrixType DimMatrix;
- // Finite element type
- typedef typename NumLib::FeQUAD4<NodalVector, DimNodalMatrix, DimMatrix>::type FeQUAD4Type;
- // Shape matrix data type
- typedef typename FeQUAD4Type::ShapeMatricesType ShapeMatricesType;
-
- public:
- NumLibFemIsoQuad4Test() :
- D(dim, dim),
- expectedM(e_nnodes,e_nnodes),
- expectedK(e_nnodes,e_nnodes),
- integration_method(2)
- {
- // create a quad element used for testing
- unitSquareQuad = createSquareQuad(1.0);
-
- // set a conductivity tensor
- setIdentityMatrix(dim, D);
- D *= conductivity;
-
- // set expected matrices
- setExpectedMassMatrix(expectedM);
- setExpectedLaplaceMatrix(conductivity, expectedK);
-
- // for destructor
- vec_eles.push_back(unitSquareQuad);
- for (auto e : vec_eles)
- for (unsigned i=0; i<e->getNNodes(true); i++)
- vec_nodes.push_back(e->getNode(i));
- }
-
- ~NumLibFemIsoQuad4Test()
- {
- for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr )
- delete *itr;
- for (auto itr = vec_eles.begin(); itr!=vec_eles.end(); ++itr )
- delete *itr;
- }
-
- // Quad: square shape with a length of 1m
- MeshLib::Quad* createSquareQuad(double h)
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node( h, 0.0, 0.0);
- nodes[2] = new MeshLib::Node( h, h, 0.0);
- nodes[3] = new MeshLib::Node(0.0, h, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- // set an identity matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setIdentityMatrix(unsigned dim, T_MATRIX &m) const
- {
- for (unsigned i=0; i<dim; i++)
- for (unsigned j=0; j<dim; j++)
- m(i,j) = 0.0;
- for (unsigned i=0; i<dim; i++)
- m(i,i) = 1.0;
- }
-
- // copy upper triangles to lower triangles in a matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m) const
- {
- for (ID_TYPE i=0; i<dim; i++)
- for (ID_TYPE j=0; j<i; j++)
- m(i,j) = m(j,i);
- }
-
- // set an expected mass matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
- m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
- m(2,2) = 1.0; m(2,3) = 1./2;
- m(3,3) = 1.0;
- // make symmetric
- copyUpperToLower(4, m);
- m *= 1./9.;
- }
-
- // set an expected laplace matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
- m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
- m(2,2) = 4.0; m(2,3) = -1.0;
- m(3,3) = 4.0;
- // make symmetric
- copyUpperToLower(4, m);
- m *= k/6.;
- }
-
- static const double conductivity;
- static const double eps;
- DimMatrix D;
- NodalMatrix expectedM;
- NodalMatrix expectedK;
- typename FeQUAD4Type::IntegrationMethod integration_method;
-
- std::vector<const MeshLib::Node*> vec_nodes;
- std::vector<const MeshLib::Quad*> vec_eles;
- MeshLib::Quad* unitSquareQuad;
-
-}; // NumLibFemIsoQuad4Test
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoQuad4Test<T_MATRIX_TYPES>::conductivity = 1e-11;
-
-template <class T_MATRIX_TYPES>
-const double NumLibFemIsoQuad4Test<T_MATRIX_TYPES>::eps = std::numeric_limits<double>::epsilon();
-
-} // namespace
-
-#ifdef OGS_USE_EIGEN
-
-struct EigenFixedMatrixTypes
-{
- typedef Eigen::Matrix<double, e_nnodes, e_nnodes, Eigen::RowMajor> NodalMatrixType;
- typedef Eigen::Matrix<double, e_nnodes, 1> NodalVectorType;
- typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
- typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrixType;
-};
-
-struct EigenDynamicMatrixTypes
-{
- typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
- typedef NodalMatrixType DimNodalMatrixType;
- typedef NodalMatrixType DimMatrixType;
- typedef Eigen::VectorXd NodalVectorType;
-};
-
-#endif // OGS_USE_EIGEN
-
-typedef ::testing::Types<
-#ifdef OGS_USE_EIGEN
- EigenFixedMatrixTypes
- , EigenDynamicMatrixTypes
-#endif
- > MatrixTypes;
-
-TYPED_TEST_CASE(NumLibFemIsoQuad4Test, MatrixTypes);
-
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeQUAD4Type fe(*this->unitSquareQuad);
-
- // evaluate a mass matrix M = int{ N^T D N }dA_e
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::N_J>(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
-
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckLaplaceMatrix)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeQUAD4Type fe(*this->unitSquareQuad);
-
- // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckMassLaplaceMatrices)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object
- FeQUAD4Type fe(*this->unitSquareQuad);
-
- // evaluate both mass and laplace matrices at once
- NodalMatrix M(e_nnodes, e_nnodes);
- NodalMatrix K(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
- ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
-}
-
-TYPED_TEST(NumLibFemIsoQuad4Test, CheckGaussIntegrationLevel)
-{
- // Refer to typedefs in the fixture
- typedef typename TestFixture::FeQUAD4Type FeQUAD4Type;
- typedef typename TestFixture::NodalMatrix NodalMatrix;
- typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
-
- // create a finite element object with gauss quadrature level 2
- FeQUAD4Type fe(*this->unitSquareQuad);
-
- // evaluate a mass matrix
- NodalMatrix M(e_nnodes, e_nnodes);
- ShapeMatricesType shape(dim, e_nnodes);
- ASSERT_EQ(4u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-
- // Change gauss quadrature level to 3
- this->integration_method.setIntegrationOrder(3);
- M *= .0;
- ASSERT_EQ(9u, this->integration_method.getNPoints());
- for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
- shape.setZero();
- auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
- M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
- }
- ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
-}
-
-