Test for material classes that compares the linearization and hessian of the...

Test for material classes that compares the linearization and hessian of the material energy with finite difference approximations

[[Imported from SVN: r11299]]

test/materialtest.cc

+#include<config.h>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.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/grid/io/file/amirameshreader.hh>
+
+
+#include <dune/fufem/makesphere.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/elasticity/materials/neohookeanmaterial.hh>
+#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
+#include <dune/elasticity/assemblers/neohookeassembler.hh>
+
+const int dim = 3;
+
+typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorType;
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim> > MatrixType;
+
+template<class Basis, class Material>
+bool checkHessian(const Basis& basis, Material& material,double eps)
+{
+
+    VectorType conf(basis.getGridView().size(dim));
+    conf = 0;
+    typedef BasisGridFunction<Basis,VectorType> GridFunc;
+    std::shared_ptr<GridFunc> confG(new GridFunc(basis, conf));
+
+    // assemble Hessian using the material assemblers
+    typename Material::LocalHessian& locHessian = material.secondDerivative();
+    locHessian.setConfiguration(confG);
+
+    MatrixType exactHessian;
+    OperatorAssembler<Basis,Basis> opAss(basis,basis);
+    opAss.assemble(locHessian, exactHessian);
+    /*
+    Dune::NeoHookeAssembler<Basis,Basis,6> neoAss(basis,basis);
+    neoAss.setEandNu(1e5,0.3);
+
+    MatrixType A;
+    VectorType b = conf;
+    neoAss.assembleProblem(A,conf,b);
+    */
+    // compute finite difference approximation
+ 
+    // assemble Hessian using the material assemblers
+    typename Material::LocalLinearization& locLin = material.firstDerivative();
+    locLin.setConfiguration(confG);
+
+    // 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);
+
+    // now the values
+    MatrixType 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++)  {
+
+
+            VectorType forward  = conf;
+            VectorType backward = conf;
+
+            forward[i][j]  += eps;
+            backward[i][j] -= eps;
+            
+            std::shared_ptr<GridFunc> fdConf(new GridFunc(basis,forward));
+            locLin.setConfiguration(fdConf);
+    
+            VectorType forw;
+            funcAss.assemble(locLin,forw);           
+    
+            std::shared_ptr<GridFunc> fdConfB(new GridFunc(basis,backward));
+            locLin.setConfiguration(fdConfB);
+
+            VectorType backw;
+            funcAss.assemble(locLin,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);
+
+            fdConfB.reset();
+            fdConf.reset();
+        }
+    
+    MatrixType diff = fdMat;
+    diff -= exactHessian;
+    std::cout<<"Maximal difference "<<std::setprecision(12)<<diff.infinity_norm()<<"\n";
+/*
+    MatrixType diff2 = fdMat;
+    diff2 -= A;
+    std::cout<<"Maximal difference 2"<<std::setprecision(12)<<diff2.infinity_norm()<<"\n";
+*/
+    if (diff.infinity_norm()>1e-8)
+        return false;
+
+    return true;
+
+}
+
+template<class Basis, class Material>
+bool checkLinearization(const Basis& basis, Material& material,double eps)
+{
+
+    VectorType conf(basis.getGridView().size(dim));
+    conf = 0;
+    typedef BasisGridFunction<Basis,VectorType> GridFunc;
+    std::shared_ptr<GridFunc> confG(new GridFunc(basis, conf));
+
+    // assemble Hessian using the material assemblers
+    typename Material::LocalLinearization& locLinear = material.firstDerivative();
+    locLinear.setConfiguration(confG);
+
+    VectorType exactLinear;
+    FunctionalAssembler<Basis> funcAss(basis);
+    funcAss.assemble(locLinear, exactLinear);
+    /*
+    Dune::NeoHookeAssembler<Basis,Basis,6> neoAss(basis,basis);
+    neoAss.setEandNu(1e5,0.3);
+
+    MatrixType A;
+    VectorType b = conf;
+    neoAss.assembleProblem(A,conf,b);
+    */
+    // compute finite difference approximation
+
+    VectorType fdLin = conf;
+
+    for (size_t i=0; i<conf.size(); i++)
+        for (int j=0; j<dim; j++)  {
+
+
+            VectorType forward  = conf;
+            VectorType backward = conf;
+
+            forward[i][j]  += eps;
+            backward[i][j] -= eps;
+
+            double forE = material.energy(forward);
+            double backE = material.energy(backward);
+
+            fdLin[i][j] = (forE-backE)/(2*eps);
+        }
+    
+    VectorType diff = fdLin;
+    diff -= exactLinear;
+    std::cout<<"Two norm "<<std::setprecision(12)<<diff.two_norm()<<"\n";
+    /*
+    VectorType diff2 = fdLin;
+    diff2 += b;
+    std::cout<<"Two norm 2 "<<std::setprecision(12)<<diff2.two_norm()<<"\n";
+    */
+    if (diff.two_norm()>1e-8)
+        return false;
+
+    return true;
+
+}
+
+
+using namespace Dune;
+
+int main (int argc, char *argv[]) try
+{
+
+    // parse eps
+    double eps = 1e-6;
+    if (argc==2)
+        eps = atof(argv[1]);
+    std::cout<<"eps "<<eps<<std::endl;
+
+    typedef UGGrid<dim> GridType;    
+    GridType* grid = new GridType;
+
+    
+    // make sphere grid
+    FieldVector<double,dim> center(0);
+    double radius = 1;
+            
+    makeSphere(*grid,center,radius);
+    
+    grid->setRefinementType(GridType::COPY);
+    for (int i=0; i<1; i++)
+        grid->globalRefine(1);
+
+    /*
+    typedef DeformationFunction<GridType::LeafGridView,VectorType> DeformationFunction;
+    typedef GeometryGrid<GridType,DeformationFunction> DeformedGridType;
+
+    grid = AmiraMeshReader<GridType>::read("/home/cocktail/akbib/data/tire/tire/cube.grid");
+
+      // initialize the deformed grids with zero displacements
+    VectorType displace(grid->size(dim));
+    displace=0;
+    DeformationFunction* deformation = new DeformationFunction(grid->leafView(),displace);
+    DeformedGridType* deformedGrid = new DeformedGridType(grid,deformation);
+
+    grid->setRefinementType(GridType::COPY);
+    //for (int i=0; i<2; i++)
+     //   grid->globalRefine(1);
+
+    // adapt deformed grids
+    displace.resize(grid->size(dim));
+    displace=0;
+    deformation->setDeformation(displace);
+    deformedGrid->update(); 
+    */
+
+    // Create the materials
+    typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
+    P1Basis p1Basis(grid->leafView());
+
+    // make a material
+    typedef NeoHookeMaterial<P1Basis> MaterialType;
+    //typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
+
+    MaterialType material;
+    material.setup(1e5,0.3,p1Basis);
+    
+    checkLinearization(p1Basis,material,eps);
+    checkHessian(p1Basis,material,eps);
+
+
+
+} catch(Exception e) {
+
+std::cout<<e<<std::endl;
+
+}