diff --git a/test/materialtest.cc b/test/materialtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9c6eb16792928a5728a837b61bfb4b2fd8cccd2d
--- /dev/null
+++ b/test/materialtest.cc
@@ -0,0 +1,248 @@
+#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;
+
+}