Skip to content
Snippets Groups Projects
Forked from agnumpde / dune-elasticity
168 commits behind the upstream repository.
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;
}