Revision e75c84b4d4c8a26c254c972e0ad3b35075980133 authored by Thomas Fischer on 18 June 2021, 06:22:59 UTC, committed by Thomas Fischer on 18 June 2021, 15:56:34 UTC
1 parent 3908c72
TestShapeFunctions.cpp
/**
* \file
* \author Norihiro Watanabe
* \date 2013-09-03
*
* \copyright
* Copyright (c) 2012-2021, 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 <algorithm>
#include <limits>
#include <numeric>
#include <type_traits>
#include <valarray>
#include "MeshLib/Elements/Elements.h"
#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
#include "Tests/AutoCheckTools.h"
#include "Tests/TestTools.h"
namespace autocheck
{
template <typename ShapeFunction,
typename Gen = randomTupleGenerator<double, ShapeFunction::DIM>>
struct NaturalPointGenerator
{
Gen generator; // The default interval is fixed to [-1, 1].
using result_type = std::array<double, ShapeFunction::DIM>;
result_type intervalMap(result_type const& tuple) const
{
if (std::is_same_v<ShapeFunction, NumLib::ShapeLine2> ||
std::is_same_v<ShapeFunction, NumLib::ShapeLine3> ||
std::is_same_v<ShapeFunction, NumLib::ShapeQuad4> ||
std::is_same_v<ShapeFunction, NumLib::ShapeQuad8> ||
std::is_same_v<ShapeFunction, NumLib::ShapeQuad9> ||
std::is_same_v<ShapeFunction, NumLib::ShapeHex8> ||
std::is_same_v<ShapeFunction, NumLib::ShapeHex20>)
{
return tuple;
}
if (std::is_same_v<ShapeFunction, NumLib::ShapeTri3> ||
std::is_same_v<ShapeFunction, NumLib::ShapeTri6> ||
std::is_same_v<ShapeFunction, NumLib::ShapeTet4> ||
std::is_same_v<ShapeFunction, NumLib::ShapeTet10>)
{
// Map square (x, y) \in [-1, 1]^2 to a triangle such that x,y
// \in [0, 1], and x+y \in [0, 2/2].
// Same for three-d case with x + y + z \in [0, 3/2].
result_type mapped_tuple;
std::transform(std::begin(tuple), std::end(tuple),
std::begin(mapped_tuple),
[](double const& v) { return (v + 1.) / 2.; });
if (std::accumulate(std::begin(mapped_tuple),
std::end(mapped_tuple),
0.) > ShapeFunction::DIM / 2.)
{
std::transform(std::begin(mapped_tuple), std::end(mapped_tuple),
std::begin(mapped_tuple),
[](double const& v) { return 1. - v; });
}
return mapped_tuple;
}
OGS_FATAL("NaturalPointGenerator: Unknown shape function type.");
}
result_type operator()(std::size_t /*size*/ = 0)
{
return intervalMap(fix(1, generator)());
}
};
} // namespace autocheck
namespace ac = autocheck;
using namespace NumLib;
template <typename ShapeFunction>
struct ShapeFunctionTest : public ::testing::Test
{
ac::NaturalPointGenerator<ShapeFunction> natural_point_generator;
ac::gtest_reporter gtest_reporter;
};
using ShapeFunctionTestTypes =
::testing::Types<ShapeLine2, ShapeLine3, ShapeTri3, ShapeTri6, ShapeTet4,
ShapeTet10, ShapeQuad4, ShapeQuad8, ShapeQuad9, ShapeHex8,
ShapeHex20>;
TYPED_TEST_SUITE(ShapeFunctionTest, ShapeFunctionTestTypes);
// TypeParam is the type of the ShapeFunction.
// Access private members via this pointer or TestFixture:: for types
TYPED_TEST(ShapeFunctionTest, PartitionOfUnity)
{
auto isPartitionOfUnity =
[](std::array<double, TypeParam::DIM>& natural_coordinates_point)
-> bool {
// compute shape functions
std::array<double, TypeParam::NPOINTS> N;
TypeParam::computeShapeFunction(natural_coordinates_point, N);
double const sum = std::accumulate(std::begin(N), std::end(N), 0.);
return std::abs(1. - sum) < 1e-15;
};
ac::check<std::array<double, TypeParam::DIM>>(
isPartitionOfUnity,
1000,
ac::make_arbitrary(this->natural_point_generator),
this->gtest_reporter);
}
TYPED_TEST(ShapeFunctionTest, SumOfGradientsIsZero)
{
auto isSumOfGradientsZero =
[](std::array<double, TypeParam::DIM>& natural_coordinates_point)
-> bool {
// compute shape functions
std::array<double, TypeParam::DIM * TypeParam::NPOINTS> dNdr;
TypeParam::computeGradShapeFunction(natural_coordinates_point, dNdr);
double const sum =
std::accumulate(std::begin(dNdr), std::end(dNdr), 0.);
return std::abs(sum) <= 5e-15;
};
ac::check<std::array<double, TypeParam::DIM>>(
isSumOfGradientsZero,
1000,
ac::make_arbitrary(this->natural_point_generator),
this->gtest_reporter);
}
TEST(NumLib, FemShapeQuad4)
{
static const double eps = std::numeric_limits<double>::epsilon();
static const unsigned NNodes = 4;
static const unsigned dim = 2;
std::valarray<double> r(dim);
std::valarray<double> N(NNodes);
std::valarray<double> dN(NNodes * dim);
// check N, dN at specific location
{
r = .5; // r = (0,5, 0.5)
ShapeQuad4::computeShapeFunction(r, N);
ShapeQuad4::computeGradShapeFunction(r, dN);
double exp_N[] = {0.5625, 0.1875, 0.0625, 0.1875};
double exp_dN[] = {0.375, -0.375, -0.125, 0.125,
0.375, 0.125, -0.125, -0.375};
ASSERT_ARRAY_NEAR(exp_N, N, N.size(), eps);
ASSERT_ARRAY_NEAR(exp_dN, dN, dN.size(), eps);
}
std::valarray<double> exp_N(NNodes);
// check N_i(r_j)= {i==j: 1, i!=j: 0}
for (unsigned i = 0; i < NNodes; i++)
{
r[0] = (i == 0 || i == 3) ? 1 : -1;
r[1] = i < 2 ? 1 : -1;
exp_N = .0;
exp_N[i] = 1.0;
ShapeQuad4::computeShapeFunction(r, N);
ASSERT_ARRAY_NEAR(exp_N, N, NNodes, eps);
}
// check sum_i[N_i(r)] = 1
const double dist = 0.5;
for (unsigned i = 0; i < NNodes; i++)
{
r[0] = (i == 0 || i == 3) ? dist : -dist;
r[1] = i < 2 ? dist : -dist;
ShapeQuad4::computeShapeFunction(r, N);
ASSERT_NEAR(1.0, N.sum(), eps);
// for(auto v: N)
// std::cout << v << " ";
// std::cout << "\n";
}
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...