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;

}