adolcmaterialtest.cc

#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/adolcneohookeanmaterial.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
    LocalNeoHookeanEnergy<GridType,Lfe> localNeoHookeEnergy(E,nu);
    typedef AdolcMaterial<FEBasis> AdolcMaterialType;
    AdolcMaterialType adolcMaterial(feBasis,localNeoHookeEnergy);

    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;

}