Forked from
agnumpde / dune-elasticity
168 commits behind the upstream repository.
-
Jonathan Youett authoredJonathan Youett authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
materialtest.cc 6.77 KiB
#include <config.h>
#include <iomanip>
#include <random>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/grid/geometrygrid/grid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/deformationfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/makering.hh>
#include <dune/fufem/makesphere.hh>
#define LAURSEN
#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include <dune/elasticity/materials/neohookeanmaterial.hh>
#include <dune/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/materials/ogdenmaterial.hh>
template <class Basis, class Material>
bool checkHessian(const Basis& basis, Material& material, double eps) {
constexpr int dim = Basis::GridView::dimension;
using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
double lower_bound{0}, upper_bound{0.005};
std::uniform_real_distribution<double> unif(lower_bound, upper_bound);
std::default_random_engine re;
Vector conf(basis.getGridView().size(dim));
for (size_t i = 0; i < conf.size(); ++i)
for (int j = 0; j < dim; ++j)
conf[i][j] = unif(re);
using GridFunc = BasisGridFunction<Basis, Vector>;
auto confG = std::make_shared<GridFunc>(basis, conf);
// assemble Hessian using the material assemblers
auto& locHessian = material.secondDerivative(confG);
OperatorAssembler<Basis, Basis> opAss(basis, basis);
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>>;
Matrix exactHessian;
opAss.assemble(locHessian, exactHessian);
// compute finite difference approximation
// dense matrix for the FD approximation
Dune::MatrixIndexSet fdIndex(conf.size(), conf.size());
for (size_t i = 0; i < conf.size(); i++)
for (size_t j = 0; j < conf.size(); ++j)
fdIndex.add(i, j);
Matrix fdMat;
fdIndex.exportIdx(fdMat);
fdMat = 0;
FunctionalAssembler<Basis> funcAss(basis);
for (size_t i = 0; i < conf.size(); ++i)
for (int j = 0; j < dim; ++j) {
auto forward = conf;
auto backward = conf;
forward[i][j] += eps;
backward[i][j] -= eps;
auto fdConf = std::make_shared<GridFunc>(basis, forward);
Vector forw;
funcAss.assemble(material.firstDerivative(fdConf), forw);
fdConf = std::make_shared<GridFunc>(basis, backward);
Vector backw;
funcAss.assemble(material.firstDerivative(fdConf), backw);
for (size_t k = 0; k < forw.size(); ++k)
for (int l = 0; l < dim; ++l)
fdMat[k][i][l][j] = (forw[k][l] - backw[k][l]) / (2 * eps);
}
auto diff = fdMat;
diff -= exactHessian;
std::cout << "Maximal difference " << std::setprecision(12) << diff.infinity_norm() << "\n";
return (diff.infinity_norm() < 1e-3);
}
template <class Basis, class Material>
bool checkLinearization(const Basis& basis, Material& material, double eps) {
constexpr int dim = Basis::GridView::dimension;
using Vector = Dune::BlockVector<Dune::FieldVector<double, dim> >;
double lower_bound{0}, upper_bound{0.05};
std::uniform_real_distribution<double> unif(lower_bound, upper_bound);
std::default_random_engine re;
Vector conf(basis.getGridView().size(dim));
for (size_t i = 0; i < conf.size(); ++i)
for (int j = 0; j < dim; ++j)
conf[i][j] = unif(re);
using GridFunc = BasisGridFunction<Basis, Vector>;
auto confG = std::make_shared<GridFunc>(basis, conf);
// assemble Hessian using the material assemblers
auto& locLinear = material.firstDerivative(confG);
FunctionalAssembler<Basis> funcAss(basis);
Vector exactLinear;
funcAss.assemble(locLinear, exactLinear);
// compute finite difference approximation
auto fdLin = conf;
auto forward = conf;
auto backward = conf;
for (size_t i = 0; i < conf.size(); ++i)
for (int j = 0; j < dim; ++j) {
forward[i][j] += eps;
backward[i][j] -= eps;
auto forE = material.energy(std::make_shared<GridFunc>(basis, forward));
auto backE = material.energy(std::make_shared<GridFunc>(basis, backward));
fdLin[i][j] = (forE - backE) / (2 * eps);
forward[i][j] -= eps;
backward[i][j] += eps;
}
auto diff = fdLin;
diff -= exactLinear;
std::cout << "Two norm " << std::setprecision(12) << diff.two_norm() << "\n";
return (diff.two_norm() < 1e-4);
}
using namespace Dune;
int main(int argc, char* argv[]) try {
Dune::MPIHelper::instance(argc, argv);
// parse eps
double eps = 1e-6;
if (argc == 2) eps = atof(argv[1]);
std::cout << "eps " << eps << std::endl;
using Grid = UGGrid<3>;
// make sphere grid
FieldVector<double, 3> center(0);
double radius = 1;
auto grid = makeSphere<Grid>(center, radius);
grid->setRefinementType(Grid::COPY);
for (int i = 0; i < 1; i++)
grid->globalRefine(1);
// Create the materials
using P1Basis = P1NodalBasis<Grid::LeafGridView>;
P1Basis p1Basis(grid->leafGridView());
// make a material
/*
using MoonRivlinMaterial = MooneyRivlinMaterial<P1Basis>;
using MaterialType = AdolcMaterial<P1Basis>;
MoonRivlinMaterial materialM;
materialM.setup(p1Basis, 1e5, 0.3);
MaterialType material(p1Basis, materialM);
*/
using MaterialType = NeoHookeanMaterial<P1Basis>;
// typedef OgdenMaterial<P1Basis> MaterialType;
// typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
MaterialType material(p1Basis, 1e5, 0.3);
//materialM.setup(p1Basis, 1e5, 0.3);
bool passed = checkLinearization(p1Basis, material, eps);
passed = checkHessian(p1Basis, material, eps) and passed;
std::cout << "Checking 2D case!\n";
using Grid2 = UGGrid<2>;
FieldVector<double, 2> cent2(0);
auto grid2 = makeRingSegment2D<Grid2>(cent2, 0.2, 1, 0, 1.5, 10, false);
grid2->setRefinementType(Grid2::COPY);
for (int i = 0; i < 1; i++)
grid2->globalRefine(1);
// Create the materials
using P1Basis2 = P1NodalBasis<Grid2::LeafGridView>;
P1Basis2 p1Basis2(grid2->leafGridView());
// make a material
//using MaterialType2 = MooneyRivlinMaterial<P1Basis2>;
using MaterialType2 = NeoHookeanMaterial<P1Basis2>;
// typedef OgdenMaterial<P1Basis2> MaterialType2;
// typedef GeomExactStVenantMaterial<P1Basis2> MaterialType2;
//MaterialType2 material2(p1Basis2, 1e5, 0.3);
MaterialType2 material2(p1Basis2, 1e5, 0.3);
passed = checkLinearization(p1Basis2, material2, eps) and passed;
passed = checkHessian(p1Basis2, material2, eps) and passed;
return !passed;
} catch (Exception e) {
std::cout << e << std::endl;
return 1;
}