Skip to content
Snippets Groups Projects
Commit 7a801766 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Add test that checks the linearization/hessian of a material class

against what Adol-C computes.

Turns out that there is a bug in the NeoHook-hessian assembler class.
parent 84853f79
No related branches found
No related tags found
No related merge requests found
#include<config.h>
#include <fstream>
#include <dune/fufem/assemblers/localassemblers/adolclocalenergy.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/timer.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/geometrygrid/grid.hh>
#include <dune/fufem/makesphere.hh>
#include <dune/fufem/makering.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/deformationfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/adolclinearizationassembler.hh>
#include <dune/fufem/assemblers/localassemblers/adolchessianassembler.hh>
//#define LAURSEN
#include <dune/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/materials/neohookeanmaterial.hh>
//#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
// grid dimension
const int dim = 3;
using namespace Dune;
int main (int argc, char *argv[]) try
{
double E = 1e5;
if (argc==2)
E = atof(argv[1]);
typedef FieldVector<double,dim> FVector;
typedef FieldMatrix<double,dim,dim> FMatrix;
typedef BlockVector<FVector> VectorType;
typedef BCRSMatrix<FMatrix> MatrixType;
typedef std::vector<FVector> CoeffType;
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
typedef UGGrid<dim> GridType;
GridType* grid = new GridType;
makeSphere(*grid,FVector(0),0.2);
grid->globalRefine(3);
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
typedef P1NodalBasis<GridView> FEBasis;
FEBasis feBasis(gridView);
// //////////////////////////
// Initial iterate
// //////////////////////////
std::ifstream infile("adolctest.sol");
std::string line;
VectorType x(feBasis.size());
int i=0;
while (std::getline(infile, line))
{
std::istringstream iss(line);
for (int j=0; j<dim;j++)
iss >> x[i][j];
i++;
}
auto xGridFunc = make_shared<BasisGridFunction<FEBasis,VectorType> >(feBasis,x);
// ////////////////////////////////////////////////////////////
// Create an assembler for the energy functional
// ////////////////////////////////////////////////////////////
double nu = 0.4;
typedef FEBasis::LocalFiniteElement Lfe;
// make a material
typedef NeoHookeanMaterial<FEBasis> MaterialType;
MaterialType material(feBasis,E,nu);
// check against the Adol-C material
typedef NeoHookeMaterial<FEBasis> AdolcMaterialType;
AdolcMaterialType adolcMaterial(feBasis,E,nu);
auto& linPaperAssembler = material.firstDerivative(xGridFunc);
auto& linAdolcAssembler = adolcMaterial.firstDerivative(xGridFunc);
FunctionalAssembler<FEBasis> functionalAssembler(feBasis);
VectorType adolcGradient(feBasis.size());
VectorType paperGrad(feBasis.size());
adolcGradient = 0;
paperGrad = 0;
Dune::Timer time;
functionalAssembler.assemble(linAdolcAssembler,adolcGradient);
std::cout<<"ADOL-C functional assembler needed "<<time.stop()<<std::endl;
time.reset();
time.start();
functionalAssembler.assemble(linPaperAssembler,paperGrad);
std::cout<<"Paper&Pen functional assembler needed "<<time.stop()<<std::endl;
for (size_t i=0; i<adolcGradient.size();i++) {
FVector diff = adolcGradient[i];
diff -= paperGrad[i];
if (diff.two_norm()>1e-9)
DUNE_THROW(Dune::Exception,"Wrong local derivative, error is "<<diff.two_norm());
}
// ////////////////////////
// Test hessian assembler
// ////////////////////////
auto& hessPaperAssembler = material.secondDerivative(xGridFunc);
auto& hessAdolcAssembler = adolcMaterial.secondDerivative(xGridFunc);
OperatorAssembler<FEBasis,FEBasis> operatorAssembler(feBasis,feBasis);
MatrixType adolcHessian, paperHessian;
time.reset();
time.start();
operatorAssembler.assemble(hessAdolcAssembler, adolcHessian);
std::cout<<"ADOL-C operator assembler needed "<<time.stop()<<std::endl;
time.reset();
time.start();
operatorAssembler.assemble(hessPaperAssembler, paperHessian);
std::cout<<"Paper&Pen operator assembler needed "<<time.stop()<<std::endl;
for (size_t i=0; i<adolcHessian.N();i++) {
const auto& adolcRow = adolcHessian[i];
const auto& paperRow = paperHessian[i];
auto adolcIt = adolcRow.begin();
auto adolcEndIt = adolcRow.end();
auto paperIt = paperRow.begin();
auto paperEndIt = paperRow.end();
for (; adolcIt != adolcEndIt; ++adolcIt, ++paperIt) {
if (adolcIt.index() != paperIt.index())
DUNE_THROW(Dune::Exception,"Not the same sparsity pattern!"<<adolcIt.index()
<<"!="<<paperIt.index());
FMatrix diff = *adolcIt;
diff -= *paperIt;
if (diff.frobenius_norm()>1e-8)
DUNE_THROW(Dune::Exception,"Wrong local hessian, error is "<<diff.frobenius_norm());
}
assert(paperIt==paperEndIt);
}
// //////////////////////////////
} catch (Exception e) {
std::cout << e << std::endl;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment