diff -Naur -x '*.svn*' -x 'Makefile*' stokes/dgstokes.cc stokesnew/dgstokes.cc --- stokes/dgstokes.cc 2008-07-16 11:31:34.899597117 +0200 +++ stokesnew/dgstokes.cc 2008-07-16 11:29:30.951069151 +0200 @@ -19,7 +19,7 @@ Dune::GeometryType gt = ent.type(); //specify the quadrature order ? // #warning fixed quadrature order - int qord=12; + int qord=6; for (int nqp=0;nqp::rule(gt,qord).size();++nqp) { //local position of quad points @@ -148,12 +148,12 @@ //get parameter // DGStokesParameters parameter; //get the geometry type of the face - Dune::GeometryType gtface = isit.intersectionSelfLocal().type(); + Dune::GeometryType gtface = isit->intersectionSelfLocal().type(); //std::cout<<"----gtface: "<intersectionNeighborLocal().type(); //specify the quadrature order ? // #warning now fixed quadrature order - int qord=12; + int qord=6; // Grid grid; @@ -162,20 +162,20 @@ //quadrature position on the edge/face in local=facelocal //note that this is dim entity const Dune::FieldVector& local = Dune::QuadratureRules::rule(gtface,qord)[qedg].position(); - Dune:: FieldVector face_self_local = isit.intersectionSelfLocal().global(local); - Dune:: FieldVector face_neighbor_local = isit.intersectionNeighborLocal().global(local); + Dune:: FieldVector face_self_local = isit->intersectionSelfLocal().global(local); + Dune:: FieldVector face_neighbor_local = isit->intersectionNeighborLocal().global(local); - Dune::FieldVector global = isit.intersectionGlobal().global(local); + Dune::FieldVector global = isit->intersectionGlobal().global(local); //std::cout<<"local: "<::rule(gtface,qord)[qedg].weight(); - ctype detjacface = isit.intersectionGlobal().integrationElement(local); + ctype detjacface = isit->intersectionGlobal().integrationElement(local); // get the face normal: unit normal. - Dune::FieldVector normal = isit.unitOuterNormal(local); - ctype norm_e= isit.intersectionGlobal().integrationElement(local); + Dune::FieldVector normal = isit->unitOuterNormal(local); + ctype norm_e= isit->intersectionGlobal().integrationElement(local); @@ -435,23 +435,23 @@ //get the geometry type of the face - Dune::GeometryType gtboundary = isit.intersectionSelfLocal().type(); + Dune::GeometryType gtboundary = isit->intersectionSelfLocal().type(); //specify the quadrature order ? - int qord=12; + int qord=6; for(int bq=0;bq::rule(gtboundary,qord).size();++bq) { const Dune::FieldVector& boundlocal = Dune::QuadratureRules::rule(gtboundary,qord)[bq].position(); - Dune:: FieldVector blocal = isit.intersectionSelfLocal().global(boundlocal); - const Dune::FieldVector bglobal = isit.intersectionGlobal().global(boundlocal); - ctype norm_eb=isit.intersectionGlobal().integrationElement(boundlocal); + Dune:: FieldVector blocal = isit->intersectionSelfLocal().global(boundlocal); + const Dune::FieldVector bglobal = isit->intersectionGlobal().global(boundlocal); + ctype norm_eb=isit->intersectionGlobal().integrationElement(boundlocal); // calculating the inverse jacobian InverseJacobianMatrix inv_jac= ent.geometry().jacobianInverseTransposed(blocal); // get quadrature weight ctype quad_wt_bound = Dune::QuadratureRules::rule(gtboundary,qord)[bq].weight(); - ctype detjacbound = isit.intersectionGlobal().integrationElement(boundlocal); + ctype detjacbound = isit->intersectionGlobal().integrationElement(boundlocal); // get the boundary normal - Dune::FieldVector boundnormal = isit.unitOuterNormal(boundlocal); + Dune::FieldVector boundnormal = isit->unitOuterNormal(boundlocal); // velocity boundary condition // dirichlet boundary for(int i=0;ineighbor()) { - //mit.insert(grid.levelIndexSet(level).index(*iit.outside())); - mit.insert(grid.leafIndexSet().index(*iit.outside())); + //mit.insert(grid.levelIndexSet(level).index(*iit->outside())); + mit.insert(grid.leafIndexSet().index(*iit->outside())); } } ++mit; @@ -713,18 +713,18 @@ for(; is != endis; ++is) { - if(is.neighbor()) + if(is->neighbor()) { - // int eid = grid.levelIndexSet(level).index(*is.inside()); -// int fid = grid.levelIndexSet(level).index(*is.outside()); + // int eid = grid.levelIndexSet(level).index(*is->inside()); +// int fid = grid.levelIndexSet(level).index(*is->outside()); - int eid = grid.leafIndexSet().index(*is.inside()); - int fid = grid.leafIndexSet().index(*is.outside()); + int eid = grid.leafIndexSet().index(*is->inside()); + int fid = grid.leafIndexSet().index(*is->outside()); dgfem.assembleFaceTerm(*it,is,A[eid][eid],A[eid][fid],A[fid][eid],b[eid]); } - if (is.boundary()) + if (is->boundary()) { dgfem.assembleBoundaryTerm(*it,is,A[eid][eid],b[eid]); } @@ -828,15 +828,15 @@ //cgsolver.apply(solution,b,r); //printvector(std::cout,solution,"Solution",""); b=solution;//for l2 error calculation - for(typename LVector::iterator i=solution.begin();i!=solution.end();++i) - { - for(int n=0;n +void Dune::DGStokes::convertToCellData(int variable, BlockVector >& cellData) +{ + ElementLevelIterator it = grid.template lbegin<0>(level); + ElementLevelIterator itend = grid.template lend<0>(level); + for (; it != itend; ++it) + { + GeometryType gt = it->geometry().type(); + const FieldVector& local = ReferenceElements::general(gt).position(0, 0); + + int eid = grid.levelIndexSet(level).index(*it); + cellData[eid] = dgfem.evaluateSolution(variable, *it, local, b[eid]); + } + + return; +} + +template +void Dune::DGStokes::vtkout (const char* name, int k) +{ + VTKWriter::LevelIndexSet> + vtkwriter(grid, grid.levelIndexSet(level)); + typedef Dune::BlockVector > BlockVector; + BlockVector pressureCellData(grid.size(0)); + convertToCellData(dim, pressureCellData); + vtkwriter.addCellData(pressureCellData, "pressure"); + BlockVector velXCellData(grid.size(0)); + convertToCellData(0, velXCellData); + vtkwriter.addCellData(velXCellData, "velocity x-comp"); + BlockVector velYCellData(grid.size(0)); + convertToCellData(1, velYCellData); + vtkwriter.addCellData(velYCellData, "velocity y-comp"); + if (dim > 2) { + BlockVector velZCellData(grid.size(0)); + convertToCellData(2, velZCellData); + vtkwriter.addCellData(velZCellData, "velocity z-comp"); + } + char fname[128]; + sprintf(fname, "%s-%05d", name, k); + vtkwriter.write(fname, Dune::VTKOptions::ascii); +} diff -Naur -x '*.svn*' -x 'Makefile*' stokes/dgstokes.hh stokesnew/dgstokes.hh --- stokes/dgstokes.hh 2008-07-16 11:31:34.871596997 +0200 +++ stokesnew/dgstokes.hh 2008-07-16 11:29:30.951069151 +0200 @@ -15,8 +15,6 @@ #include #include -//#include -//#include #include"stokesparameters.hh" #include"boundaryconditions.hh" @@ -50,10 +48,8 @@ enum{p_ordr=p_order}; //local vector and matrix blocks - // local block size is sum of velocity dof and pressure dof - //block size = dim*vdof.size() + pdof.size() static const int BlockSize =dim*Dune::MonomialShapeFunctionSetSize::maxSize+Dune::MonomialShapeFunctionSetSize::maxSize; -static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize::maxSize)); + static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize::maxSize)); static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize::maxSize)); typedef Dune::FieldVector LocalVectorBlock; @@ -100,12 +96,9 @@ Grid& grid; Dune::MonomialShapeFunctionSetContainer vspace; Dune::MonomialShapeFunctionSetContainer pspace; -// ShapeFunctionSet vsfs(v_order); //for velocity -// ShapeFunctionSet psfs(p_order); // for pressure DGStokesParameters parameter; DirichletBoundary dirichletvalue; RightHandSide rhsvalue; - }; @@ -135,7 +128,7 @@ typedef typename Grid::template Codim<1>::EntityPointer InterSectionPointer; static const int BlockSize =((dim*Dune::MonomialShapeFunctionSetSize::maxSize)+(Dune::MonomialShapeFunctionSetSize::maxSize)); -static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize::maxSize)); + static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize::maxSize)); static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize::maxSize)); typedef typename DGFiniteElementMethod::Gradient Gradient; @@ -143,17 +136,11 @@ typedef typename DGFiniteElementMethod::LocalMatrixBlock LocalMatrixBlock; typedef Dune::BlockVector LVector; typedef Dune::BCRSMatrix LMatrix; - typedef Dune::LevelDGStokesVelocityFunction DGSVFunction; - typedef Dune::LevelDGStokesPressureFunction DGSPFunction; - //typedef Dune::LevelDGFunction DGFunction; public: - //enum {order = ordr}; - //inline constructor with initializer list - DGStokes(Grid &g, ExactSolution& ex, DGStokesParameters& par, DirichletBoundary& db, - RightHandSide& rh, int l ) : grid(g),exact(ex), dgfem(g,par,db,rh),level(l),xv(grid,level),xp(grid,level) + RightHandSide& rh, int l ) : grid(g),exact(ex), dgfem(g,par,db,rh),level(l) { if (par.sigma==1 & par.epsilon==1) std::cout<<"You are using NIPG scheme"< >& cellData); + void vtkout (const char* name, int k); - private: typedef typename DGFiniteElementMethod::ShapeFunctionSet ShapeFunctionSet; inline const ShapeFunctionSet & getVelocityShapeFunctionSet(const EntityPointer &ep) const; @@ -193,21 +173,11 @@ Grid & grid; int level; ExactSolution& exact; - //for velocity - // ShapeFunctionSet& vsfs(ordr); -// // for pressure -// ShapeFunctionSet& psfs(ordr-1); private: DGFiniteElementMethod dgfem; LMatrix A; LVector b; LVector solution; - DGSVFunction xv; - DGSPFunction xp; - //DGFunction x; - - - }; diff -Naur -x '*.svn*' -x 'Makefile*' stokes/l2error.hh stokesnew/l2error.hh --- stokes/l2error.hh 2008-07-16 11:31:34.879597033 +0200 +++ stokesnew/l2error.hh 2008-07-16 11:29:30.951069151 +0200 @@ -15,15 +15,12 @@ Dune::GeometryType gt = element.type(); // #warning fixed quadrature order int qord=12; -//G grid; int eid = grid.levelIndexSet(grid.maxLevel()).index(element); - //std::cout<<"EID:--> "<::rule(gt,qord).size();++qp) { qp_loc = Dune::QuadratureRules::rule(gt,qord)[qp].position(); - qp_glob =element.geometry().global(qp_loc); - //std::cout<<"qp_loc: "<::rule(gt,qord)[qp].weight(); double detjac = element.geometry().integrationElement(qp_loc); if (variable(level); ElementLevelIterator itend = grid.template lend<0>(level); - //ElementLeafIterator it = grid.template leafbegin<0>(); - //ElementLeafIterator itend = grid.template leafend<0>(); - //std::cout<<"level:::"< { typedef Dune::FieldVector< ct, dim > Point; + typedef Dune::FieldVector Gradient; public: Example1(){} @@ -102,6 +103,12 @@ if (variable==1) return -glob[1]*cos(glob[0])+glob[0];// -y*cos(x)+x if (variable==2) return 0.0; } + Gradient velocityGradient(int comp,const Point &glob)const + { + // //DUNE_THROW(NotImplemented, "velocityGradient not implemented yet"); + //return 0; + } + virtual ~Example1(){} }; @@ -110,6 +117,7 @@ class Example2 : public ExactSolution { typedef Dune::FieldVector< ct, dim > Point; + typedef Dune::FieldVector Gradient; public: Example2(){} @@ -131,6 +139,12 @@ if (variable==1) return 0.0; if (variable==2) return 0.0; } + Gradient velocityGradient(int comp,const Point &glob)const + { + // //DUNE_THROW(NotImplemented, "velocityGradient not implemented yet"); + //return 0; + } + virtual ~Example2(){} };