Revision 3b5cec3d6ec5aafa79f45902e999e53649f3084d authored by Lars Bilke on 06 February 2023, 08:02:08 UTC, committed by Lars Bilke on 06 February 2023, 08:02:08 UTC
Two minor updates See merge request ogs/ogs!4472
TestFe.cpp
/**
* \copyright
* Copyright (c) 2012-2023, 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 <Eigen/Core>
#include <cmath>
#include <vector>
#include "FeTestData/TestFeHEX8.h"
#include "FeTestData/TestFeLINE2.h"
#include "FeTestData/TestFeLINE2Y.h"
#include "FeTestData/TestFeLINE3.h"
#include "FeTestData/TestFePRISM6.h"
#include "FeTestData/TestFePYRA5.h"
#include "FeTestData/TestFeQUAD4.h"
#include "FeTestData/TestFeTET4.h"
#include "FeTestData/TestFeTRI3.h"
#include "MeshLib/CoordinateSystem.h"
#include "MeshLib/ElementCoordinatesMappingLocal.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "Tests/TestTools.h"
using namespace NumLib;
using namespace FeTestData;
namespace
{
// test cases
template <class TestFeType_, template <typename, int> class ShapeMatrixPolicy_>
struct TestCase
{
using TestFeType = TestFeType_;
static const int GlobalDim = TestFeType::global_dim;
using ShapeMatrixTypes =
ShapeMatrixPolicy_<typename TestFeType::ShapeFunction, GlobalDim>;
template <typename X>
using ShapeMatrixPolicy = ShapeMatrixPolicy_<X, GlobalDim>;
using ShapeFunction = typename TestFeType::ShapeFunction;
};
using TestTypes =
::testing::Types<TestCase<TestFeHEX8, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeLINE2, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeLINE2Y, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeLINE3, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFePRISM6, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFePYRA5, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeQUAD4, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeTET4, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeTRI3, EigenDynamicShapeMatrixPolicy>,
TestCase<TestFeHEX8, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeLINE2, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeLINE2Y, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeLINE3, EigenFixedShapeMatrixPolicy>,
TestCase<TestFePRISM6, EigenFixedShapeMatrixPolicy>,
TestCase<TestFePYRA5, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeQUAD4, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeTET4, EigenFixedShapeMatrixPolicy>,
TestCase<TestFeTRI3, EigenFixedShapeMatrixPolicy>>;
} // namespace
template <class T>
class NumLibFemIsoTest : public ::testing::Test, public T::TestFeType
{
public:
using ShapeMatrixTypes = typename T::ShapeMatrixTypes;
using ShapeFunction = typename T::ShapeFunction;
using TestFeType = typename T::TestFeType;
// Matrix types
using NodalMatrix = typename ShapeMatrixTypes::NodalMatrixType;
using NodalVector = typename ShapeMatrixTypes::NodalVectorType;
using DimNodalMatrix = typename ShapeMatrixTypes::DimNodalMatrixType;
using DimMatrix = typename ShapeMatrixTypes::DimMatrixType;
using GlobalDimMatrixType = typename ShapeMatrixTypes::GlobalDimMatrixType;
// Finite element type
template <typename X>
using ShapeMatrixPolicy = typename T::template ShapeMatrixPolicy<X>;
using FeType =
typename TestFeType::template FeType<ShapeMatrixPolicy>::type;
// Shape matrix data type
using ShapeMatricesType = typename ShapeMatrixTypes::ShapeMatrices;
using MeshElementType = typename TestFeType::MeshElementType;
static const unsigned dim = TestFeType::dim;
static const unsigned e_nnodes = TestFeType::e_nnodes;
static const unsigned n_sample_pt_order2 = TestFeType::n_sample_pt_order2;
static const unsigned n_sample_pt_order3 = TestFeType::n_sample_pt_order3;
using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
MeshElementType>::IntegrationMethod;
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;
MeshLib::ElementCoordinatesMappingLocal ele_local_coord(
*mesh_element,
MeshLib::CoordinateSystem(*mesh_element).getDimension());
auto R = ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
TestFeType::dim, TestFeType::global_dim);
globalD.noalias() = R.transpose() * (D * R);
// 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->getNumberOfNodes(); i++)
{
vec_nodes.push_back(e->getNode(i));
}
}
}
~NumLibFemIsoTest() override
{
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;
IntegrationMethod integration_method;
GlobalDimMatrixType globalD;
std::vector<const MeshLib::Node*> vec_nodes;
std::vector<const MeshElementType*> vec_eles;
MeshElementType* mesh_element;
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
}; // NumLibFemIsoTest
template <class T>
const double NumLibFemIsoTest<T>::conductivity = 1e-11;
template <class T>
const double NumLibFemIsoTest<T>::eps = 2 *
std::numeric_limits<double>::epsilon();
template <class T>
const unsigned NumLibFemIsoTest<T>::dim;
template <class T>
const unsigned NumLibFemIsoTest<T>::e_nnodes;
template <class T>
const unsigned NumLibFemIsoTest<T>::n_sample_pt_order2;
template <class T>
const unsigned NumLibFemIsoTest<T>::n_sample_pt_order3;
TYPED_TEST_SUITE(NumLibFemIsoTest, TestTypes);
TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
{
// Refer to typedefs in the fixture
using NodalMatrix = typename TestFixture::NodalMatrix;
auto const shape_matrices =
NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
typename TestFixture::ShapeMatrixTypes,
TestFixture::global_dim,
ShapeMatrixType::N_J>(
*this->mesh_element, false /*is_axially_symmetric*/,
this->integration_method);
// evaluate a mass matrix M = int{ N^T D N }dA_e
NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
for (std::size_t i = 0; i < this->integration_method.getNumberOfPoints();
i++)
{
auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
M.noalias() +=
shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
}
// std::cout << "M=\n" << M;
ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
}
TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
{
// Refer to typedefs in the fixture
using NodalMatrix = typename TestFixture::NodalMatrix;
auto const shape_matrices =
NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
typename TestFixture::ShapeMatrixTypes,
TestFixture::global_dim,
ShapeMatrixType::DNDX>(
*this->mesh_element, false /*is_axially_symmetric*/,
this->integration_method);
// evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
for (std::size_t i = 0; i < this->integration_method.getNumberOfPoints();
i++)
{
auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
K.noalias() += shape.dNdx.transpose() * this->globalD * shape.dNdx *
shape.detJ * wp.getWeight();
}
// std::cout << "K=\n" << K << std::endl;
ASSERT_ARRAY_NEAR(this->expectedK.data(), K.data(), K.size(), this->eps);
}
TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
{
// Refer to typedefs in the fixture
using NodalMatrix = typename TestFixture::NodalMatrix;
auto const shape_matrices =
NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
typename TestFixture::ShapeMatrixTypes,
TestFixture::global_dim>(
*this->mesh_element, false /*is_axially_symmetric*/,
this->integration_method);
// evaluate both mass and laplace matrices at once
NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
for (std::size_t i = 0; i < this->integration_method.getNumberOfPoints();
i++)
{
auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
M.noalias() +=
shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
K.noalias() += shape.dNdx.transpose() * (this->globalD * 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);
}
#if 0
TYPED_TEST(NumLibFemIsoTest, CheckGaussLegendreIntegrationLevel)
{
// Refer to typedefs in the fixture
using NodalMatrix = typename TestFixture::NodalMatrix;
auto shape_matrices =
NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
typename TestFixture::ShapeMatrixTypes,
TestFixture::global_dim>(
*this->mesh_element, false /*is_axially_symmetric*/,
this->integration_method);
// evaluate a mass matrix
NodalMatrix M(this->e_nnodes, this->e_nnodes);
M.setZero();
ASSERT_EQ(TestFixture::n_sample_pt_order2, this->integration_method.getNumberOfPoints());
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M << std::endl;
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(TestFixture::n_sample_pt_order3, this->integration_method.getNumberOfPoints());
shape_matrices =
NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
typename TestFixture::ShapeMatrixTypes,
TestFixture::global_dim>(
*this->mesh_element, false /*is_axially_symmetric*/,
this->integration_method);
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M << std::endl;
ASSERT_ARRAY_NEAR(this->expectedM.data(), M.data(), M.size(), this->eps);
}
#endif
Computing file changes ...