diff --git a/test/materialtest.cc b/test/materialtest.cc index 9c6eb16792928a5728a837b61bfb4b2fd8cccd2d..2b8a66ca3ce4acaa4d02095bd686f5175ce51b3d 100644 --- a/test/materialtest.cc +++ b/test/materialtest.cc @@ -9,10 +9,10 @@ #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/makering.hh> #include <dune/fufem/functions/basisgridfunction.hh> #include <dune/fufem/functions/deformationfunction.hh> #include <dune/fufem/functionspacebases/p1nodalbasis.hh> @@ -23,15 +23,17 @@ #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) { + const int dim = Basis::GridView::dimension; + + typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorType; + typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim> > MatrixType; + VectorType conf(basis.getGridView().size(dim)); conf = 0; typedef BasisGridFunction<Basis,VectorType> GridFunc; @@ -44,7 +46,7 @@ bool checkHessian(const Basis& basis, Material& material,double eps) 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); @@ -120,6 +122,10 @@ bool checkHessian(const Basis& basis, Material& material,double eps) template<class Basis, class Material> bool checkLinearization(const Basis& basis, Material& material,double eps) { + const int dim = Basis::GridView::dimension; + + typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorType; + typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim> > MatrixType; VectorType conf(basis.getGridView().size(dim)); conf = 0; @@ -188,12 +194,12 @@ int main (int argc, char *argv[]) try eps = atof(argv[1]); std::cout<<"eps "<<eps<<std::endl; - typedef UGGrid<dim> GridType; + typedef UGGrid<3> GridType; GridType* grid = new GridType; // make sphere grid - FieldVector<double,dim> center(0); + FieldVector<double,3> center(0); double radius = 1; makeSphere(*grid,center,radius); @@ -239,6 +245,31 @@ int main (int argc, char *argv[]) try checkLinearization(p1Basis,material,eps); checkHessian(p1Basis,material,eps); + std::cout<<"Checking 2D case!\n"; + + typedef UGGrid<2> GridType2; + GridType2* grid2 = new GridType2; + + FieldVector<double,2> cent2(0); + + makeRingSegment2D(cent2, 0.2, 1, 0, 1.5, *grid2 , 10, false); + + grid2->setRefinementType(GridType2::COPY); + for (int i=0; i<1; i++) + grid2->globalRefine(1); + + // Create the materials + typedef P1NodalBasis<GridType2::LeafGridView> P1Basis2; + P1Basis2 p1Basis2(grid2->leafView()); + + // make a material + typedef NeoHookeMaterial<P1Basis2> MaterialType2; + + MaterialType2 material2; + material2.setup(1e5,0.3,p1Basis2); + + checkLinearization(p1Basis2,material2,eps); + checkHessian(p1Basis2,material2,eps); } catch(Exception e) {