Skip to content
Snippets Groups Projects
Commit d5de01dc authored by Norihiro Watanabe's avatar Norihiro Watanabe
Browse files

merge tests for Fe types to one

parent 80d9feb9
No related branches found
No related tags found
No related merge requests found
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
#include <Eigen/Eigen> #include <Eigen/Eigen>
#endif #endif
#include "MeshLib/Elements/Hex.h" #include "MeshLib/Elements/Line.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h" #include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h" #include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
...@@ -26,92 +26,144 @@ using namespace NumLib; ...@@ -26,92 +26,144 @@ using namespace NumLib;
namespace namespace
{ {
static const unsigned dim = 3; // copy matrix entries in upper triangle to lower triangle
static const unsigned e_nnodes = 8; 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);
}
template <class T_MATRIX_TYPES> // set an identity matrix
class NumLibFemIsoHex8Test : public ::testing::Test template <class T_MATRIX, typename ID_TYPE=signed>
void setIdentityMatrix(unsigned dim, T_MATRIX &m)
{ {
public: for (unsigned i=0; i<dim; i++)
// Matrix types for (unsigned j=0; j<dim; j++)
typedef typename T_MATRIX_TYPES::NodalMatrixType NodalMatrix; m(i,j) = 0.0;
typedef typename T_MATRIX_TYPES::NodalVectorType NodalVector; for (unsigned i=0; i<dim; i++)
typedef typename T_MATRIX_TYPES::DimNodalMatrixType DimNodalMatrix; m(i,i) = 1.0;
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: // Test case for LINE2
NumLibFemIsoHex8Test() : class TestFeLINE2
D(dim, dim), {
expectedM(e_nnodes,e_nnodes), public:
expectedK(e_nnodes,e_nnodes), template <class T_MATRIX_TYPES>
integration_method(2) 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()
{ {
// create a quad element used for testing MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
mesh_element = createCubicHex(1.0); nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
// set a conductivity tensor return new MeshLib::Line(nodes);
setIdentityMatrix(dim, D); }
D *= conductivity;
// set expected matrices
setExpectedMassMatrix(expectedM);
setExpectedLaplaceMatrix(conductivity, expectedK);
// for destructor // set an expected mass matrix
vec_eles.push_back(mesh_element); template <class T_MATRIX, typename ID_TYPE=signed>
for (auto e : vec_eles) void setExpectedMassMatrix(T_MATRIX &m)
for (unsigned i=0; i<e->getNNodes(true); i++) {
vec_nodes.push_back(e->getNode(i)); // set upper triangle entries
m(0,0) = 1./3.; m(0,1) = 1./6.;
m(1,1) = 1./3.;
// make symmetric
copyUpperToLower(2, m);
} }
~NumLibFemIsoHex8Test() // set an expected laplace matrix
template <class T_MATRIX, typename ID_TYPE=signed>
void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
{ {
for (auto itr = vec_nodes.begin(); itr!=vec_nodes.end(); ++itr ) // set upper triangle entries
delete *itr; m(0,0) = 1.0; m(0,1) = -1.0;
for (auto itr = vec_eles.begin(); itr!=vec_eles.end(); ++itr ) m(1,1) = 1.0;
delete *itr; // make symmetric
copyUpperToLower(2, m);
m *= k;
} }
};
// Quad: square shape with a length of 1m // Test case for QUAD4
MeshLib::Hex* createCubicHex(double h) 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]; MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0); nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
nodes[1] = new MeshLib::Node( h, 0.0, 0.0); nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
nodes[2] = new MeshLib::Node( h, h, 0.0); nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
nodes[3] = new MeshLib::Node(0.0, h, 0.0); nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
nodes[4] = new MeshLib::Node(0.0, 0.0, h); return new MeshLib::Quad(nodes);
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 // set an expected mass matrix for 1m x 1m
template <class T_MATRIX, typename ID_TYPE=signed> template <class T_MATRIX, typename ID_TYPE=signed>
void setIdentityMatrix(unsigned dim, T_MATRIX &m) const void setExpectedMassMatrix(T_MATRIX &m)
{ {
for (unsigned i=0; i<dim; i++) // set upper triangle entries
for (unsigned j=0; j<dim; j++) m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
m(i,j) = 0.0; m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
for (unsigned i=0; i<dim; i++) m(2,2) = 1.0; m(2,3) = 1./2;
m(i,i) = 1.0; m(3,3) = 1.0;
// make symmetric
copyUpperToLower(4, m);
m *= 1./9.;
} }
// copy upper triangles to lower triangles in a matrix // set an expected laplace matrix for 1m x 1m
template <class T_MATRIX, typename ID_TYPE=signed> template <class T_MATRIX, typename ID_TYPE=signed>
void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m) const void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
{ {
for (ID_TYPE i=0; i<dim; i++) // set upper triangle entries
for (ID_TYPE j=0; j<i; j++) m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
m(i,j) = m(j,i); 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.;
} }
};
// set an expected mass matrix for 1m x 1m // 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> template <class T_MATRIX, typename ID_TYPE=signed>
void setExpectedMassMatrix(T_MATRIX &m) void setExpectedMassMatrix(T_MATRIX &m)
{ {
...@@ -129,7 +181,7 @@ class NumLibFemIsoHex8Test : public ::testing::Test ...@@ -129,7 +181,7 @@ class NumLibFemIsoHex8Test : public ::testing::Test
m *= 1./27.; m *= 1./27.;
} }
// set an expected laplace matrix for 1m x 1m // set an expected laplace matrix
template <class T_MATRIX, typename ID_TYPE=signed> template <class T_MATRIX, typename ID_TYPE=signed>
void setExpectedLaplaceMatrix(double k, T_MATRIX &m) void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
{ {
...@@ -146,36 +198,97 @@ class NumLibFemIsoHex8Test : public ::testing::Test ...@@ -146,36 +198,97 @@ class NumLibFemIsoHex8Test : public ::testing::Test
copyUpperToLower(e_nnodes, m); copyUpperToLower(e_nnodes, m);
m *= k/2.; 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 conductivity;
static const double eps; static const double eps;
DimMatrix D; DimMatrix D;
NodalMatrix expectedM; NodalMatrix expectedM;
NodalMatrix expectedK; NodalMatrix expectedK;
typename FeHEX8Type::IntegrationMethod integration_method; typename FeType::IntegrationMethod integration_method;
std::vector<const MeshLib::Node*> vec_nodes; std::vector<const MeshLib::Node*> vec_nodes;
std::vector<const MeshLib::Hex*> vec_eles; std::vector<const MeshElementType*> vec_eles;
MeshLib::Hex* mesh_element; MeshElementType* mesh_element;
}; // NumLibFemIsoQuad4Test }; // NumLibFemIsoLineTest
template <class T_MATRIX_TYPES> template <class T_MATRIX_TYPES, class T_FE>
const double NumLibFemIsoHex8Test<T_MATRIX_TYPES>::conductivity = 1e-11; const double NumLibFemIsoTest<T_MATRIX_TYPES, T_FE>::conductivity = 1e-11;
template <class T_MATRIX_TYPES> template <class T_MATRIX_TYPES, class T_FE>
const double NumLibFemIsoHex8Test<T_MATRIX_TYPES>::eps = std::numeric_limits<double>::epsilon(); 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 } // namespace
#ifdef OGS_USE_EIGEN #ifdef OGS_USE_EIGEN
template <unsigned NNODES, unsigned DIM>
struct EigenFixedMatrixTypes struct EigenFixedMatrixTypes
{ {
typedef Eigen::Matrix<double, e_nnodes, e_nnodes, Eigen::RowMajor> NodalMatrixType; typedef Eigen::Matrix<double, NNODES, DIM, Eigen::RowMajor> NodalMatrixType;
typedef Eigen::Matrix<double, e_nnodes, 1> NodalVectorType; typedef Eigen::Matrix<double, NNODES, 1> NodalVectorType;
typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrixType; typedef Eigen::Matrix<double, DIM, NNODES, Eigen::RowMajor> DimNodalMatrixType;
typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrixType; typedef Eigen::Matrix<double, DIM, DIM, Eigen::RowMajor> DimMatrixType;
}; };
struct EigenDynamicMatrixTypes struct EigenDynamicMatrixTypes
...@@ -186,30 +299,48 @@ struct EigenDynamicMatrixTypes ...@@ -186,30 +299,48 @@ struct EigenDynamicMatrixTypes
typedef Eigen::VectorXd NodalVectorType; typedef Eigen::VectorXd NodalVectorType;
}; };
template <class T_FE>
using NumLibFemIsoTestFeTypes = NumLibFemIsoTest<EigenDynamicMatrixTypes, T_FE>;
#endif // OGS_USE_EIGEN #endif // OGS_USE_EIGEN
typedef ::testing::Types< template <class T_MATRIX_TYPES>
using NumLibFemIsoTestMatTypes = NumLibFemIsoTest<T_MATRIX_TYPES, TestFeLINE2>;
template <class T>
struct MatrixTypes
{
typedef ::testing::Types<
#ifdef OGS_USE_EIGEN #ifdef OGS_USE_EIGEN
EigenFixedMatrixTypes EigenFixedMatrixTypes<T::e_nnodes, T::dim>
, EigenDynamicMatrixTypes , EigenDynamicMatrixTypes
#endif #endif
> MatrixTypes; > types;
};
TYPED_TEST_CASE(NumLibFemIsoHex8Test, MatrixTypes);
TYPED_TEST(NumLibFemIsoHex8Test, CheckMassMatrix) 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 // Refer to typedefs in the fixture
typedef typename TestFixture::FeHEX8Type FeHEX8Type; typedef typename TestFixture::FeType FeType;
typedef typename TestFixture::NodalMatrix NodalMatrix; typedef typename TestFixture::NodalMatrix NodalMatrix;
typedef typename TestFixture::ShapeMatricesType ShapeMatricesType; typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
// create a finite element object // create a finite element object
FeHEX8Type fe(*this->mesh_element); FeType fe(*this->mesh_element);
// evaluate a mass matrix M = int{ N^T D N }dA_e // evaluate a mass matrix M = int{ N^T D N }dA_e
NodalMatrix M(e_nnodes, e_nnodes); NodalMatrix M(this->e_nnodes, this->e_nnodes);
ShapeMatricesType shape(dim, e_nnodes); ShapeMatricesType shape(this->dim, this->e_nnodes);
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) { for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero(); shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i); auto wp = this->integration_method.getWeightedPoint(i);
...@@ -217,50 +348,45 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckMassMatrix) ...@@ -217,50 +348,45 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckMassMatrix)
M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight(); 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); ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
} }
TYPED_TEST(NumLibFemIsoHex8Test, CheckLaplaceMatrix) TYPED_TEST(NumLibFemIsoTestFeTypes, CheckLaplaceMatrix)
{ {
// Refer to typedefs in the fixture // Refer to typedefs in the fixture
typedef typename TestFixture::FeHEX8Type FeHEX8Type; typedef typename TestFixture::FeType FeType;
typedef typename TestFixture::NodalMatrix NodalMatrix; typedef typename TestFixture::NodalMatrix NodalMatrix;
typedef typename TestFixture::ShapeMatricesType ShapeMatricesType; typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
// create a finite element object // create a finite element object
FeHEX8Type fe(*this->mesh_element); FeType fe(*this->mesh_element);
// evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
NodalMatrix K(e_nnodes, e_nnodes); NodalMatrix K(this->e_nnodes, this->e_nnodes);
ShapeMatricesType shape(dim, e_nnodes); ShapeMatricesType shape(this->dim, this->e_nnodes);
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) { for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero(); shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i); auto wp = this->integration_method.getWeightedPoint(i);
fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape); fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
K.noalias() += shape.dNdx.transpose() * this->D * shape.dNdx * shape.detJ * wp.getWeight(); 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); ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
} }
TYPED_TEST(NumLibFemIsoHex8Test, CheckMassLaplaceMatrices) TYPED_TEST(NumLibFemIsoTestFeTypes, CheckMassLaplaceMatrices)
{ {
// Refer to typedefs in the fixture // Refer to typedefs in the fixture
typedef typename TestFixture::FeHEX8Type FeHEX8Type; typedef typename TestFixture::FeType FeType;
typedef typename TestFixture::NodalMatrix NodalMatrix; typedef typename TestFixture::NodalMatrix NodalMatrix;
typedef typename TestFixture::ShapeMatricesType ShapeMatricesType; typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
// create a finite element object // create a finite element object
FeHEX8Type fe(*this->mesh_element); FeType fe(*this->mesh_element);
// evaluate both mass and laplace matrices at once // evaluate both mass and laplace matrices at once
NodalMatrix M(e_nnodes, e_nnodes); NodalMatrix M(this->e_nnodes, this->e_nnodes);
NodalMatrix K(e_nnodes, e_nnodes); NodalMatrix K(this->e_nnodes, this->e_nnodes);
ShapeMatricesType shape(dim, e_nnodes); ShapeMatricesType shape(this->dim, this->e_nnodes);
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) { for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero(); shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i); auto wp = this->integration_method.getWeightedPoint(i);
...@@ -272,20 +398,20 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckMassLaplaceMatrices) ...@@ -272,20 +398,20 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckMassLaplaceMatrices)
ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps); ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
} }
TYPED_TEST(NumLibFemIsoHex8Test, CheckGaussIntegrationLevel) TYPED_TEST(NumLibFemIsoTestFeTypes, CheckGaussIntegrationLevel)
{ {
// Refer to typedefs in the fixture // Refer to typedefs in the fixture
typedef typename TestFixture::FeHEX8Type FeHEX8Type; typedef typename TestFixture::FeType FeType;
typedef typename TestFixture::NodalMatrix NodalMatrix; typedef typename TestFixture::NodalMatrix NodalMatrix;
typedef typename TestFixture::ShapeMatricesType ShapeMatricesType; typedef typename TestFixture::ShapeMatricesType ShapeMatricesType;
// create a finite element object with gauss quadrature level 2 // create a finite element object with gauss quadrature level 2
FeHEX8Type fe(*this->mesh_element); FeType fe(*this->mesh_element);
// evaluate a mass matrix // evaluate a mass matrix
NodalMatrix M(e_nnodes, e_nnodes); NodalMatrix M(this->e_nnodes, this->e_nnodes);
ShapeMatricesType shape(dim, e_nnodes); ShapeMatricesType shape(this->dim, this->e_nnodes);
ASSERT_EQ(8u, this->integration_method.getNPoints()); ASSERT_EQ(std::pow(2u, this->dim), this->integration_method.getNPoints());
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) { for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero(); shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i); auto wp = this->integration_method.getWeightedPoint(i);
...@@ -297,7 +423,7 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckGaussIntegrationLevel) ...@@ -297,7 +423,7 @@ TYPED_TEST(NumLibFemIsoHex8Test, CheckGaussIntegrationLevel)
// Change gauss quadrature level to 3 // Change gauss quadrature level to 3
this->integration_method.setIntegrationOrder(3); this->integration_method.setIntegrationOrder(3);
M *= .0; M *= .0;
ASSERT_EQ(27u, this->integration_method.getNPoints()); ASSERT_EQ(std::pow(3u, this->dim), this->integration_method.getNPoints());
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) { for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero(); shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i); auto wp = this->integration_method.getWeightedPoint(i);
......
/**
* \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);
}
/**
* \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);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment