Commit 6935f934 authored by akbib's avatar akbib Committed by akbib@FU-BERLIN.DE
Browse files

also test the 2D case

[[Imported from SVN: r11302]]
parent 5de5af51
......@@ -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) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment