-
The LocalIntegralStiffness assembles a tangent matrix by suitably combining derivatives of the energy density with derivatives of the geometric FE interpolation. See the detailed description in dune-gfe-manual.pdf. This patch also includes a new test in localintegralstiffnesstest.cc. It checks if the assembled matrix is the same as the one computed by LocalGeodesicFEADOLCStiffness. For a 3d Cosserat material I get a speedup of about 3x.
The LocalIntegralStiffness assembles a tangent matrix by suitably combining derivatives of the energy density with derivatives of the geometric FE interpolation. See the detailed description in dune-gfe-manual.pdf. This patch also includes a new test in localintegralstiffnesstest.cc. It checks if the assembled matrix is the same as the one computed by LocalGeodesicFEADOLCStiffness. For a 3d Cosserat material I get a speedup of about 3x.
localintegralstiffnesstest.cc 25.53 KiB
#include <config.h>
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/common/timer.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#include <dune/gfe/assemblers/sumenergy.hh>
#include <dune/gfe/assemblers/geodesicfeassembler.hh>
#include <dune/gfe/assemblers/localintegralenergy.hh>
#include <dune/gfe/assemblers/localintegralstiffness.hh>
#include <dune/gfe/densities/localdensity.hh>
#include <dune/gfe/densities/bulkcosseratdensity.hh>
#include <dune/gfe/densities/harmonicdensity.hh>
#include <dune/gfe/spaces/realtuple.hh>
// grid dimension
const int dim = 3;
using namespace Dune;
using namespace Dune::Indices;
using namespace Functions::BasisFactory;
enum InterpolationType {Geodesic, ProjectionBased};
template <class GridView, InterpolationType interpolationType>
int testHarmonicMapIntoSphere(TestSuite& test, const GridView& gridView)
{
using TargetSpace = UnitVector<double,dim>;
// Finite element order
const int order = 1;
const static int blocksize = TargetSpace::TangentVector::dimension;
using Correction = BlockVector<TargetSpace::TangentVector>;
using Matrix = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
//////////////////////////////////////////////////////////////////////////////////
// Construct all needed function space bases
//////////////////////////////////////////////////////////////////////////////////
auto tangentBasis = makeBasis(
gridView,
power<blocksize>(
lagrange<order>()
)
);
auto valueBasis = makeBasis(gridView, power<TargetSpace::CoordinateType::dimension>(lagrange<order>()));
/////////////////////////////////////////////////////////////////////////
// Construct the configuration where to assemble the tangent problem
/////////////////////////////////////////////////////////////////////////
using CoefficientVector = std::vector<TargetSpace>;
CoefficientVector x(tangentBasis.size());
BlockVector<FieldVector<double,dim> > initialIterate;
Functions::interpolate(valueBasis, initialIterate, [](FieldVector<double,dim> x)
-> FieldVector<double,dim> {
// Interpret x[0] and x[1] as spherical coordinates
auto phi = x[0];
auto theta = x[1];
return {sin(phi)*cos(theta), sin(phi)*sin(theta), cos(theta)};
});
for (size_t i = 0; i < x.size(); i++)
x[i] = initialIterate[i];
using ATargetSpace = typename TargetSpace::rebind<adouble>::other;
// Select which type of geometric interpolation to use
// TODO: Use the tangent basis here!
using FEBasis = Functions::LagrangeBasis<GridView,order>;
Functions::LagrangeBasis<GridView,order> feBasis(gridView);
//////////////////////////////////////////////////////////////
// Set up the two assemblers
//////////////////////////////////////////////////////////////
std::shared_ptr<GeodesicFEAssembler<FEBasis,TargetSpace> > assemblerADOLC;
std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,TargetSpace> > localIntegralStiffness;
using LocalCoordinate = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
auto harmonicDensity = std::make_shared<GFE::HarmonicDensity<LocalCoordinate,ATargetSpace> >();
if constexpr (interpolationType==Geodesic)
{
std::cout << "Using geodesic interpolation" << std::endl;
using LocalInterpolationRule = LocalGeodesicFEFunction<dim,
typename GridView::ctype,
decltype(feBasis.localView().tree().finiteElement()),
TargetSpace>;
using ALocalInterpolationRule = LocalGeodesicFEFunction<dim,
typename GridView::ctype,
decltype(feBasis.localView().tree().finiteElement()),
ATargetSpace>;
// Assemble using the old assembler
auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensity);
auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
// Assemble using the new assembler
localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<FEBasis,LocalInterpolationRule,TargetSpace> >(harmonicDensity);
}
else
{
std::cout << "Using projection-based interpolation" << std::endl;
using LocalInterpolationRule = GFE::LocalProjectedFEFunction<dim,
typename GridView::ctype,
decltype(feBasis.localView().tree().finiteElement()),
TargetSpace>;
using ALocalInterpolationRule = GFE::LocalProjectedFEFunction<dim,
typename GridView::ctype,
decltype(feBasis.localView().tree().finiteElement()),
ATargetSpace>;
// Assemble using the old assembler
auto energy = std::make_shared<GFE::LocalIntegralEnergy<FEBasis, ALocalInterpolationRule,ATargetSpace> >(harmonicDensity);
auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> >(energy);
assemblerADOLC = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localGFEADOLCStiffness);
// Assemble using the new assembler
localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<FEBasis,LocalInterpolationRule,TargetSpace> >(harmonicDensity);
}
// TODO: The assembler should really get the tangent basis.
auto assemblerIntegralStiffness = std::make_shared<GeodesicFEAssembler<FEBasis,TargetSpace> >(feBasis, localIntegralStiffness);
//////////////////////////////////////////////////////////////
// Assemble
//////////////////////////////////////////////////////////////
Correction gradient;
Matrix hesseMatrix;
Correction gradientSmart;
Matrix hesseMatrixSmart;
Timer assemblerTimer;
assemblerADOLC->assembleGradientAndHessian(x, gradient, hesseMatrix,true);
std::cout << "assemblerTimer took: " << assemblerTimer.elapsed() << std::endl;
Timer assemblerSmartTimer;
assemblerIntegralStiffness->assembleGradientAndHessian(x, gradientSmart, hesseMatrixSmart,true);
std::cout << "assemblerSmartTimer took: " << assemblerSmartTimer.elapsed() << std::endl;
//////////////////////////////////////////////////////////////
// Compare the results
//////////////////////////////////////////////////////////////
double gradientInfinityNorm = gradient.infinity_norm();
double matrixFrobeniusNorm = hesseMatrix.frobenius_norm();
double gradientSmartInfinityNorm = gradientSmart.infinity_norm();
double matrixSmartFrobeniusNorm = hesseMatrixSmart.frobenius_norm();
if (std::isnan(gradientInfinityNorm) || std::isnan(gradientSmartInfinityNorm) || std::abs(gradientInfinityNorm - gradientSmartInfinityNorm)/gradientInfinityNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Gradient max norm from localADOLCStiffness is " << gradientInfinityNorm
<< ", from LocalIntegralStiffness is " << gradientSmartInfinityNorm << " :(" << std::endl;
return 1;
}
if (std::isnan(matrixFrobeniusNorm) || std::isnan(matrixSmartFrobeniusNorm) || std::abs(matrixFrobeniusNorm - matrixSmartFrobeniusNorm)/matrixFrobeniusNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Matrix norm from localADOLCStiffness is " << matrixFrobeniusNorm
<< ", from LocalIntegralStiffness is " << matrixSmartFrobeniusNorm << " :(" << std::endl;
return 1;
}
// TODO: Use rangeFor
for (auto rowIt = hesseMatrixSmart.begin(); rowIt != hesseMatrixSmart.end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for (size_t k = 0; k < blocksize; k++)
for (size_t l = 0; l < blocksize; l++)
if (std::abs((*colIt)[k][l]) > 1e-6 and std::abs(hesseMatrix[rowIt.index()][colIt.index()][k][l]) > 1e-6
and std::abs((*colIt)[k][l] - hesseMatrix[rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-4) {
std::cerr << "Matrix: Index " << " " << rowIt.index() << " " << colIt.index() << " " << k << " " << l << " is wrong: " << hesseMatrix[rowIt.index()][colIt.index()][k][l] << " - " << (*colIt)[k][l] << std::endl;
return 1;
}
return 0;
}
template <class GridView, InterpolationType interpolationType>
int testCosseratBulkModel(TestSuite& test, const GridView& gridView)
{
// Order of the approximation spaces
// Deliberately test with different orders, to find more bugs
const int deformationOrder = 2;
const int rotationOrder = 1;
const static int deformationBlocksize = RealTuple<double,dim>::TangentVector::dimension;
const static int orientationBlocksize = Rotation<double,dim>::TangentVector::dimension;
using CorrectionType0 = BlockVector<FieldVector<double, deformationBlocksize> >;
using CorrectionType1 = BlockVector<FieldVector<double, orientationBlocksize> >;
using MatrixType00 = BCRSMatrix<FieldMatrix<double, deformationBlocksize, deformationBlocksize> >;
using MatrixType01 = BCRSMatrix<FieldMatrix<double, deformationBlocksize, orientationBlocksize> >;
using MatrixType10 = BCRSMatrix<FieldMatrix<double, orientationBlocksize, deformationBlocksize> >;
using MatrixType11 = BCRSMatrix<FieldMatrix<double, orientationBlocksize, orientationBlocksize> >;
using MatrixType = MultiTypeBlockMatrix<MultiTypeBlockVector<MatrixType00,MatrixType01>,
MultiTypeBlockVector<MatrixType10,MatrixType11> >;
//////////////////////////////////////////////////////////////////////////////////
// Construct all needed function space bases
//////////////////////////////////////////////////////////////////////////////////
const int dimRotation = Rotation<double,dim>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
composite(
power<dim>(
lagrange<deformationOrder>()
),
power<dimRotation>(
lagrange<rotationOrder>()
)
));
using CompositeBasis = decltype(compositeBasis);
/////////////////////////////////////////////////////////////////////////
// Construct the configuration where to assemble the tangent problem
/////////////////////////////////////////////////////////////////////////
using CoefficientVector = TupleVector<std::vector<RealTuple<double,dim> >,std::vector<Rotation<double,dim> > >;
CoefficientVector x;
x[_0].resize(compositeBasis.size({0}));
x[_1].resize(compositeBasis.size({1}));
MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >,BlockVector<FieldVector<double,dimRotation> > > initialIterate;
Functions::interpolate(Functions::subspaceBasis(compositeBasis, _0), initialIterate, [](FieldVector<double,dim> x){
FieldVector<double,dim> y = {0.5*x[0]*x[0], 0.5*x[1] + 0.2*x[0], 4*std::abs(std::sin(x[2]))};
return y;
});
for (size_t i = 0; i < x[_0].size(); i++)
x[_0][i] = initialIterate[_0][i];
Functions::interpolate(Functions::subspaceBasis(compositeBasis, _1), initialIterate, [](FieldVector<double,dim> x){
// Just any rotation field, for testing
FieldVector<double,4> y = {0.5*x[0], 0.7*x[1], 1.0*std::abs(std::sin(x[2])), 1.0 - 0.2*x[0]*x[1]};
return y;
});
for (size_t i = 0; i < x[_1].size(); i++)
x[_1][i] = Rotation<double,dim>(initialIterate[_1][i]);
// Set of material parameters
ParameterTree parameters; parameters["thickness"] = "0.1";
parameters["mu"] = "1";
parameters["lambda"] = "1";
parameters["mu_c"] = "1";
parameters["L_c"] = "0.1";
parameters["q"] = "2";
parameters["kappa"] = "1";
parameters["b1"] = "1";
parameters["b2"] = "1";
parameters["b3"] = "1";
using RigidBodyMotion = GFE::ProductManifold<RealTuple<double,dim>, Rotation<double,dim> >;
using ARigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
using LocalCoordinate = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<LocalCoordinate,double> >(parameters);
auto aBulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<LocalCoordinate,adouble> >(parameters);
// Select which type of geometric interpolation to use
using DeformationFEBasis = Functions::LagrangeBasis<GridView,deformationOrder>;
using OrientationFEBasis = Functions::LagrangeBasis<GridView,rotationOrder>;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
//////////////////////////////////////////////////////////////
// Set up the two assemblers
//////////////////////////////////////////////////////////////
std::shared_ptr<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> > assemblerADOLC;
std::shared_ptr<LocalGeodesicFEStiffness<CompositeBasis,RigidBodyMotion> > localIntegralStiffness;
if constexpr (interpolationType==Geodesic)
{
std::cout << "Using geodesic interpolation" << std::endl;
using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<double,dim> >;
using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<double,dim> >;
using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
using ALocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
using ALocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
// Assemble using the ADOL-C assembler
auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
// Assemble using the new assembler
localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(aBulkCosseratDensity);
}
else
{
std::cout << "Using projection-based interpolation" << std::endl;
using LocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<double,dim> >;
using LocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<double,dim> >;
using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
using ALocalDeformationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
using ALocalOrientationInterpolationRule = GFE::LocalProjectedFEFunction<dim, typename GridView::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
using ALocalInterpolationRule = std::tuple<ALocalDeformationInterpolationRule,ALocalOrientationInterpolationRule>;
// Assemble using the ADOL-C assembler
auto energy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ALocalInterpolationRule,ARigidBodyMotion> >(aBulkCosseratDensity);
auto localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<CompositeBasis,RigidBodyMotion> >(energy);
assemblerADOLC = std::make_shared<MixedGFEAssembler<CompositeBasis,RigidBodyMotion> >(compositeBasis, localGFEADOLCStiffness);
// Assemble using the new assembler
localIntegralStiffness = std::make_shared<GFE::LocalIntegralStiffness<CompositeBasis,LocalInterpolationRule,RigidBodyMotion> >(aBulkCosseratDensity);
}
MixedGFEAssembler<CompositeBasis,RigidBodyMotion> mixedAssemblerSmart(compositeBasis,
localIntegralStiffness);
//////////////////////////////////////////////////////////////
// Assemble
//////////////////////////////////////////////////////////////
CorrectionType0 gradient0;
CorrectionType1 gradient1;
MatrixType hesseMatrix;
CorrectionType0 gradientSmart0;
CorrectionType1 gradientSmart1;
MatrixType hesseMatrixSmart;
Timer assemblerTimer;
assemblerADOLC->assembleGradientAndHessian(x[_0], x[_1], gradient0, gradient1, hesseMatrix,true);
std::cout << "assemblerTimer took: " << assemblerTimer.elapsed() << std::endl;
Timer assemblerSmartTimer;
mixedAssemblerSmart.assembleGradientAndHessian(x[_0], x[_1], gradientSmart0, gradientSmart1, hesseMatrixSmart,true);
std::cout << "assemblerSmartTimer took: " << assemblerSmartTimer.elapsed() << std::endl;
//////////////////////////////////////////////////////////////
// Compare the results
//////////////////////////////////////////////////////////////
double gradient0InfinityNorm = gradient0.infinity_norm();
double gradient1InfinityNorm = gradient1.infinity_norm();
double matrix00FrobeniusNorm = hesseMatrix[_0][_0].frobenius_norm();
double matrix01FrobeniusNorm = hesseMatrix[_0][_1].frobenius_norm();
double matrix10FrobeniusNorm = hesseMatrix[_1][_0].frobenius_norm();
double matrix11FrobeniusNorm = hesseMatrix[_1][_1].frobenius_norm();
double gradientSmart0InfinityNorm = gradientSmart0.infinity_norm();
double gradientSmart1InfinityNorm = gradientSmart1.infinity_norm();
double matrixSmart00FrobeniusNorm = hesseMatrixSmart[_0][_0].frobenius_norm();
double matrixSmart01FrobeniusNorm = hesseMatrixSmart[_0][_1].frobenius_norm();
double matrixSmart10FrobeniusNorm = hesseMatrixSmart[_1][_0].frobenius_norm();
double matrixSmart11FrobeniusNorm = hesseMatrixSmart[_1][_1].frobenius_norm();
if (std::isnan(gradient0InfinityNorm) || std::isnan(gradientSmart0InfinityNorm) || std::abs(gradient0InfinityNorm - gradientSmart0InfinityNorm)/gradient0InfinityNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Gradient (deformation part) max norm from localADOLCStiffness is " << gradient0InfinityNorm << ", from localADOLCIntegralStiffness is " << gradientSmart0InfinityNorm << " :(" << std::endl;
return 1;
}
if (std::isnan(gradient1InfinityNorm) || std::isnan(gradientSmart1InfinityNorm) || std::abs(gradient1InfinityNorm - gradientSmart1InfinityNorm)/gradient1InfinityNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Gradient (orientation part) max norm from localADOLCStiffness is " << gradient1InfinityNorm << ", from localADOLCIntegralStiffness is " << gradientSmart1InfinityNorm << " :(" << std::endl;
return 1;
}
if (std::isnan(matrix00FrobeniusNorm) || std::isnan(matrixSmart00FrobeniusNorm) || std::abs(matrix00FrobeniusNorm - matrixSmart00FrobeniusNorm)/matrix00FrobeniusNorm > 1e-3)
{
std::cerr << std::setprecision(9); std::cerr << "Matrix00 norm from localADOLCStiffness is " << matrix00FrobeniusNorm << ", from localADOLCIntegralStiffness is " << matrixSmart00FrobeniusNorm << " :(" << std::endl;
return 1;
}
for (auto rowIt = hesseMatrixSmart[_0][_0].begin(); rowIt != hesseMatrixSmart[_0][_0].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for (size_t k = 0; k < deformationBlocksize; k++)
for (size_t l = 0; l < deformationBlocksize; l++)
if (std::abs((*colIt)[k][l]) > 1e-6 and std::abs(hesseMatrix[_0][_0][rowIt.index()][colIt.index()][k][l]) > 1e-6
and std::abs((*colIt)[k][l] - hesseMatrix[_0][_0][rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-4) {
std::cerr << "Matrix00: Index " << " " << rowIt.index() << " " << colIt.index() << " " << k << " " << l << " is wrong: " << hesseMatrix[_0][_0][rowIt.index()][colIt.index()][k][l] << " - " << (*colIt)[k][l] << std::endl;
return 1;
}
if (std::isnan(matrix01FrobeniusNorm) || std::isnan(matrixSmart01FrobeniusNorm) || std::abs(matrix01FrobeniusNorm - matrixSmart01FrobeniusNorm)/matrix01FrobeniusNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Matrix01 norm from localADOLCStiffness is " << matrix01FrobeniusNorm << ", from localADOLCIntegralStiffness is " << matrixSmart01FrobeniusNorm << " :(" << std::endl;
return 1;
}
for (auto rowIt = hesseMatrixSmart[_0][_1].begin(); rowIt != hesseMatrixSmart[_0][_1].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for (size_t k = 0; k < deformationBlocksize; k++)
for (size_t l = 0; l < deformationBlocksize; l++)
if (std::abs((*colIt)[k][l]) > 1e-6 and std::abs(hesseMatrix[_0][_1][rowIt.index()][colIt.index()][k][l]) > 1e-6
and std::abs((*colIt)[k][l] - hesseMatrix[_0][_1][rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-4) {
std::cerr << "Matrix01 Index " << " " << rowIt.index() << " " << colIt.index() << " " << k << " " << l << " is wrong: " << hesseMatrix[_0][_1][rowIt.index()][colIt.index()][k][l] << " - " << (*colIt)[k][l] << std::endl;
return 1;
}
if (std::isnan(matrix10FrobeniusNorm) || std::isnan(matrixSmart10FrobeniusNorm) || std::abs(matrix10FrobeniusNorm - matrixSmart10FrobeniusNorm)/matrix10FrobeniusNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Matrix10 norm from localADOLCStiffness is " << matrix10FrobeniusNorm << ", from localADOLCIntegralStiffness is " << matrixSmart10FrobeniusNorm << " :(" << std::endl;
return 1;
}
for (auto rowIt = hesseMatrixSmart[_1][_0].begin(); rowIt != hesseMatrixSmart[_1][_0].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for (size_t k = 0; k < deformationBlocksize; k++)
for (size_t l = 0; l < deformationBlocksize; l++)
if (std::abs((*colIt)[k][l]) > 1e-6 and std::abs(hesseMatrix[_1][_0][rowIt.index()][colIt.index()][k][l]) > 1e-6
and std::abs((*colIt)[k][l] - hesseMatrix[_1][_0][rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-4) {
std::cerr << "Matrix10 Index " << " " << rowIt.index() << " " << colIt.index() << " " << k << " " << l << " is wrong: " << hesseMatrix[_1][_0][rowIt.index()][colIt.index()][k][l] << " - " << (*colIt)[k][l] << std::endl;
return 1;
}
if (std::isnan(matrix11FrobeniusNorm) || std::isnan(matrixSmart00FrobeniusNorm) || std::abs(matrix11FrobeniusNorm - matrixSmart11FrobeniusNorm)/matrix11FrobeniusNorm > 1e-3)
{
std::cerr << std::setprecision(9);
std::cerr << "Matrix11 norm from localADOLCStiffness is " << matrix11FrobeniusNorm << ", from localADOLCIntegralStiffness is " << matrixSmart11FrobeniusNorm << " :(" << std::endl;
return 1;
}
for (auto rowIt = hesseMatrixSmart[_1][_1].begin(); rowIt != hesseMatrixSmart[_1][_1].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for (size_t k = 0; k < deformationBlocksize; k++)
for (size_t l = 0; l < deformationBlocksize; l++)
if (std::abs((*colIt)[k][l]) > 1e-6 and std::abs(hesseMatrix[_1][_1][rowIt.index()][colIt.index()][k][l]) > 1e-6
and std::abs((*colIt)[k][l] - hesseMatrix[_1][_1][rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-4) {
std::cerr << "Matrix11 Index " << " " << rowIt.index() << " " << colIt.index() << " " << k << " " << l << " is wrong: " << hesseMatrix[_1][_1][rowIt.index()][colIt.index()][k][l] << " - " << (*colIt)[k][l] << std::endl;
return 1;
}
return 0;
}
int main (int argc, char *argv[]){
MPIHelper::instance(argc, argv);
TestSuite test;
/////////////////////////////////////////
// Create the grid
/////////////////////////////////////////
// Create a grid with 2 elements, one element will get refined to create different element types, the other one stays a cube
using GridType = UGGrid<dim>;
auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {2,1,1});
//Refine once
for (auto&& e : elements(grid->leafGridView())) {
bool refineThisElement = false;
for (int i = 0; i < e.geometry().corners(); i++) {
refineThisElement = refineThisElement || (e.geometry().corner(i)[0] > 0.99);
}
grid->mark(refineThisElement ? 1 : 0, e);
}
grid->adapt();
grid->loadBalance();
using GridView = GridType::LeafGridView;
auto gridView = grid->leafGridView();
/////////////////////////////////////////////
// Test different energies
/////////////////////////////////////////////
// TODO: Use test framework
testHarmonicMapIntoSphere<GridView,Geodesic>(test, gridView);
testHarmonicMapIntoSphere<GridView,ProjectionBased>(test, gridView);
testCosseratBulkModel<GridView,Geodesic>(test, gridView);
testCosseratBulkModel<GridView,ProjectionBased>(test, gridView);
return test.exit();
}