[Dune] DG for Stokes
Sreejith Pulloor Kuttanikkad
sreejith at hal.iwr.uni-heidelberg.de
Wed Jul 16 12:55:06 CEST 2008
Hi Bernd,
thanks for the work.
I referred the following for The Implementation:
1. V. Girault, B. Riviere and M.F. Wheeler. ”A Discontinuous
Galerkin Method with Non-Overlapping
Domain Decomposition for the Stokes and Navier-Stokes Problems”,
Mathematics of Computation,
74, p. 53-84, 2005
2. B. Riviere and V. Girault. ”Discontinuous finite element
methods for incompressible flows
on subdomains with non-matching interfaces”, Computer Methods
in Applied Mechanics and
Engineering, 195 p.3274-3292, 2006
3. V. Girault, B. Riviere and M. Wheeler. ”A splitting
method using discontinuous Galerkin for the
transient incompressible Navier-Stokes equations”,
Mathematical Modelling and Numerical Analysis
(M2AN), 39 no 6, p. 1115-1148, 2005, also TR-MATH 04-08
regards,
Sreejith
On Wed, Jul 16, 2008 at 11:59:51AM +0200, Bernd Flemisch wrote:
> Hey Sreejith,
>
> thank you for your hints. I tested some of the 2D examples, it seems to
> work fine, at least the convergence rates look like they should. I also
> added a vtk-output. If you consider keeping it, please feel free to
> apply the attached patch via
> patch -p1 <stokes-patch
> in dune-disc/disc/stokes and commit it. I also attach my test file.
>
> One more question: Can you give me a reference based on which the
> implementation has been done?
>
> Best regards
> Bernd
>
> >As I said, the code is not updated for sometimes and some headers
> >r missing. It definitely has to be removed.
> >
> >>thank you for your quick reply. I just tried to include dgstokes.hh and
> >>it fails because it does not know LevelDGStokesVelocityFunction (a term
> >>which is not found anywhere else but in line 146 of dgstokes.hh). That's
> >>ok, I will try for myself then.
> >>
> >just looking at the code, if that is the only error u may comment out line
> >146,147 and 178,182 becz they are there for vtk visualization. But i cant
> >really assure.
> >
> >sreejith
> >
> >
> >
> >
> >
>
>
> --
> _____________________________________________________________________
>
> Bernd Flemisch phone: +49 711 685 69162
> IWS, Universit?t Stuttgart fax: +49 711 685 60430
> Pfaffenwaldring 61 email: bernd at iws.uni-stuttgart.de
> D-70569 Stuttgart url: www.iws.uni-stuttgart.de
> _____________________________________________________________________
>
> 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<Dune::QuadratureRules<ctype,dim>::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: "<<gtface<<std::endl;
> - Dune::GeometryType nbgtface = isit.intersectionNeighborLocal().type();
> + Dune::GeometryType nbgtface = isit->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<ctype,dim-1>& local = Dune::QuadratureRules<ctype,dim-1>::rule(gtface,qord)[qedg].position();
> - Dune:: FieldVector<ctype,dim> face_self_local = isit.intersectionSelfLocal().global(local);
> - Dune:: FieldVector<ctype,dim> face_neighbor_local = isit.intersectionNeighborLocal().global(local);
> + Dune:: FieldVector<ctype,dim> face_self_local = isit->intersectionSelfLocal().global(local);
> + Dune:: FieldVector<ctype,dim> face_neighbor_local = isit->intersectionNeighborLocal().global(local);
>
> - Dune::FieldVector<ctype,dim> global = isit.intersectionGlobal().global(local);
> + Dune::FieldVector<ctype,dim> global = isit->intersectionGlobal().global(local);
> //std::cout<<"local: "<<local<<" face_self_local: "<<face_self_local
> //<<" face_neighbor_local: "<<face_neighbor_local<<" glob: "<<global<<std::endl;
> // calculating the inverse jacobian
> InverseJacobianMatrix inv_jac= ent.geometry().jacobianInverseTransposed(face_self_local);
> // get quadrature weight
> ctype quad_wt_face = Dune::QuadratureRules<ctype,dim-1>::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<ctype,dim> normal = isit.unitOuterNormal(local);
> - ctype norm_e= isit.intersectionGlobal().integrationElement(local);
> + Dune::FieldVector<ctype,dim> 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<Dune::QuadratureRules<ctype,dim-1>::rule(gtboundary,qord).size();++bq)
> {
> const Dune::FieldVector<ctype,dim-1>& boundlocal = Dune::QuadratureRules<ctype,dim-1>::rule(gtboundary,qord)[bq].position();
> - Dune:: FieldVector<ctype,dim> blocal = isit.intersectionSelfLocal().global(boundlocal);
> - const Dune::FieldVector<ctype,dim> bglobal = isit.intersectionGlobal().global(boundlocal);
> - ctype norm_eb=isit.intersectionGlobal().integrationElement(boundlocal);
> + Dune:: FieldVector<ctype,dim> blocal = isit->intersectionSelfLocal().global(boundlocal);
> + const Dune::FieldVector<ctype,dim> 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<ctype,dim-1>::rule(gtboundary,qord)[bq].weight();
> - ctype detjacbound = isit.intersectionGlobal().integrationElement(boundlocal);
> + ctype detjacbound = isit->intersectionGlobal().integrationElement(boundlocal);
> // get the boundary normal
> - Dune::FieldVector<ctype,dim> boundnormal = isit.unitOuterNormal(boundlocal);
> + Dune::FieldVector<ctype,dim> boundnormal = isit->unitOuterNormal(boundlocal);
> // velocity boundary condition
> // dirichlet boundary
> for(int i=0;i<dim;++i)
> @@ -667,10 +667,10 @@
> for(; iit != endit; ++iit)
> {
>
> - if (iit.neighbor())
> + if (iit->neighbor())
> {
> - //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<VBlockSize;++n)
> - (*xv)[i.index()][n]=solution[i.index()][n];
> - for(int n=0;n<PBlockSize;++n)
> - (*xp)[i.index()][n]=solution[i.index()][VBlockSize+n];
> -
> - }
> - // printvector(std::cout,solution,"Solution","");
> +// for(typename LVector::iterator i=solution.begin();i!=solution.end();++i)
> +// {
> +// for(int n=0;n<VBlockSize;++n)
> +// (*xv)[i.index()][n]=solution[i.index()][n];
> +// for(int n=0;n<PBlockSize;++n)
> +// (*xp)[i.index()][n]=solution[i.index()][VBlockSize+n];
> +//
> +// }
> + //printvector(std::cout,solution,"Solution","");
> //printvector(std::cout,(*xv),"Velocity Coeff","");
> //printvector(std::cout,(*xp),"Pressure Coeff","");
> //------------------ISTL solver--------------------------//
> @@ -950,17 +950,51 @@
> grad[variable][sd] += grad_phi_ei[variable][sd] * xe[ii];
>
> }
> - // else if(variable==dim)
> -// for (int i=0; i<npsfs; ++i)
> -// {
> -// int ii=(dim*nvsfs)+i;
> -// //std::cout<<"val: "<<value[variable];
> -// value[variable] += xe[ii] * psfs[i].evaluateFunction(0, coord);
> -
> -// }
> +
> return grad[variable];
> }
>
> +template<class G,int v_order, int p_order>
> +void Dune::DGStokes<G,v_order,p_order>::convertToCellData(int variable, BlockVector<FieldVector<double, 1> >& 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<ctype,dim>& local = ReferenceElements<ctype,dim>::general(gt).position(0, 0);
> +
> + int eid = grid.levelIndexSet(level).index(*it);
> + cellData[eid] = dgfem.evaluateSolution(variable, *it, local, b[eid]);
> + }
> +
> + return;
> +}
> +
> +template<class G,int v_order, int p_order>
> +void Dune::DGStokes<G,v_order,p_order>::vtkout (const char* name, int k)
> +{
> + VTKWriter<G, typename G::template Codim<0>::LevelIndexSet>
> + vtkwriter(grid, grid.levelIndexSet(level));
> + typedef Dune::BlockVector<Dune::FieldVector<double, 1> > 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<dune/grid/common/quadraturerules.hh>
>
> #include<dune/disc/shapefunctions/dgspace/monomialshapefunctions.hh>
> -//#include <dune/disc/functions/dgstokesfunction.hh>
> -//#include <dune/disc/functions/dgfunction.hh>
>
> #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<dim,v_order>::maxSize+Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize;
> -static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
> + static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
> static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
>
> typedef Dune::FieldVector<double,BlockSize> LocalVectorBlock;
> @@ -100,12 +96,9 @@
> Grid& grid;
> Dune::MonomialShapeFunctionSetContainer<ctype,double,dim,v_order> vspace;
> Dune::MonomialShapeFunctionSetContainer<ctype,double,dim,(p_order)> pspace;
> -// ShapeFunctionSet vsfs(v_order); //for velocity
> -// ShapeFunctionSet psfs(p_order); // for pressure
> DGStokesParameters parameter;
> DirichletBoundary<G> dirichletvalue;
> RightHandSide<G> rhsvalue;
> -
> };
>
>
> @@ -135,7 +128,7 @@
> typedef typename Grid::template Codim<1>::EntityPointer InterSectionPointer;
> static const int BlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize)+(Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
>
> -static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
> + static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
> static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
>
> typedef typename DGFiniteElementMethod<G,v_order,p_order>::Gradient Gradient;
> @@ -143,17 +136,11 @@
> typedef typename DGFiniteElementMethod<G,v_order,p_order>::LocalMatrixBlock LocalMatrixBlock;
> typedef Dune::BlockVector<LocalVectorBlock> LVector;
> typedef Dune::BCRSMatrix<LocalMatrixBlock> LMatrix;
> - typedef Dune::LevelDGStokesVelocityFunction<Grid, double, v_order> DGSVFunction;
> - typedef Dune::LevelDGStokesPressureFunction<Grid, double, p_order> DGSPFunction;
>
> - //typedef Dune::LevelDGFunction<Grid, double, v_order> DGFunction;
> public:
> - //enum {order = ordr};
> - //inline constructor with initializer list
> -
> DGStokes(Grid &g, ExactSolution<ctype,dim>& ex,
> DGStokesParameters& par, DirichletBoundary<Grid>& db,
> - RightHandSide<Grid>& rh, int l ) : grid(g),exact(ex), dgfem(g,par,db,rh),level(l),xv(grid,level),xp(grid,level)
> + RightHandSide<Grid>& 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"<<std::endl;
> @@ -175,16 +162,9 @@
> // stokes system has dim+1 variables (dim velocity comps and 1 pressure)
> double l2errorStokesSystem(int variable) const;
> double h1errorStokesSystem(int variable) const;
> - const DGSVFunction & velocity() const
> - {
> - return xv;
> - }
> - const DGSPFunction & pressure() const
> - {
> - return xp;
> - }
> + void convertToCellData(int variable, BlockVector<FieldVector<double, 1> >& cellData);
> + void vtkout (const char* name, int k);
>
> -
> private:
> typedef typename DGFiniteElementMethod<G,v_order,p_order>::ShapeFunctionSet ShapeFunctionSet;
> inline const ShapeFunctionSet & getVelocityShapeFunctionSet(const EntityPointer &ep) const;
> @@ -193,21 +173,11 @@
> Grid & grid;
> int level;
> ExactSolution<ctype,dim>& exact;
> - //for velocity
> - // ShapeFunctionSet& vsfs(ordr);
> -// // for pressure
> -// ShapeFunctionSet& psfs(ordr-1);
> private:
> DGFiniteElementMethod<G,v_order,p_order> 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:--> "<<eid<<std::endl;
> +
> for (int qp=0;qp<Dune::QuadratureRules<ctype,dim>::rule(gt,qord).size();++qp)
> {
> qp_loc = Dune::QuadratureRules<ctype,dim>::rule(gt,qord)[qp].position();
> -
> qp_glob =element.geometry().global(qp_loc);
> - //std::cout<<"qp_loc: "<<qp_loc<<" qp_glob: "<<qp_glob<<std::endl;
> double weight = Dune::QuadratureRules<ctype,dim>::rule(gt,qord)[qp].weight();
> double detjac = element.geometry().integrationElement(qp_loc);
> if (variable<dim)
> @@ -31,15 +28,9 @@
> error[variable]+=weight*detjac
> *(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe))
> *(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe));
> -
> -
> - // if(variable==0)
> -// std::cout<<"u:"<<exact.velocity(0,qp_glob)<<" comp: "<<evaluateSolution(0,element,qp_loc,xe)<<std::endl;
> }
> if(variable==dim)
> {
> -
> - std::cout<<"qp: "<<qp_glob<<" p:"<<exact.pressure(qp_glob)<<" comp:"<<evaluateSolution(variable,element,qp_loc,xe)<<std::endl;
> error[variable]+=weight*detjac
> *(exact.pressure(qp_glob)-evaluateSolution(variable,element,qp_loc,xe))
> *(exact.pressure(qp_glob)-evaluateSolution(variable,element,qp_loc,xe));
> @@ -59,13 +50,9 @@
> error[variable]= 0.0;
> ElementLevelIterator it = grid.template lbegin<0>(level);
> ElementLevelIterator itend = grid.template lend<0>(level);
> - //ElementLeafIterator it = grid.template leafbegin<0>();
> - //ElementLeafIterator itend = grid.template leafend<0>();
> - //std::cout<<"level:::"<<level;
> for (; it != itend; ++it)
> {
> int eid = grid.levelIndexSet(level).index(*it);
> - //std::cout<<" eid: "<<eid<<std::endl;
> error[variable]+=dgfem.evaluateL2error(variable,exact,*it,b[eid]);
> }
> return sqrt(error[variable]);
> diff -Naur -x '*.svn*' -x 'Makefile*' stokes/testfunctions.cc stokesnew/testfunctions.cc
> --- stokes/testfunctions.cc 2008-07-16 11:31:34.939597289 +0200
> +++ stokesnew/testfunctions.cc 2008-07-16 11:29:30.951069151 +0200
> @@ -81,6 +81,7 @@
> class Example1 : public ExactSolution<ct, dim>
> {
> typedef Dune::FieldVector< ct, dim > Point;
> + typedef Dune::FieldVector<ct,dim> 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<ct, dim>
> {
> typedef Dune::FieldVector< ct, dim > Point;
> + typedef Dune::FieldVector<ct,dim> 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(){}
> };
>
> #include "config.h"
> #include <iostream>
> #include <iomanip>
> #include <dune/grid/utility/gridtype.hh>
> #include <dune/grid/common/gridinfo.hh>
> #include <dune/grid/sgrid.hh>
> #include <dune/grid/io/file/dgfparser/dgfparser.hh>
> #include <dune/grid/io/file/dgfparser/dgfs.hh>
> #include <dune/grid/io/file/vtk/vtkwriter.hh>
> #include <dune/istl/io.hh>
> #include <dune/common/timer.hh>
> #include <dune/disc/stokes/dgstokes.hh>
> #include <dune/disc/stokes/l2error.hh>
>
> int main(int argc, char** argv)
> {
> try{
> // define the problem dimensions
> const int dim=2;
> const int vOrder = 2;
> const int pOrder = 1;
>
> // create a grid object
> typedef double NumberType;
> typedef Dune::SGrid<dim,dim> GridType;
>
> if (argc != 2 && argc != 3) {
> std::cout << "Usage: test_stokes dgffilename [refinementsteps]" << std::endl;
> return (1);
> }
> int refinementSteps = 0;
> if (argc == 3) {
> std::string arg2(argv[2]);
> std::istringstream is2(arg2);
> is2 >> refinementSteps;
> }
>
> Dune::GridPtr<GridType> gridPtr( argv[1] );
> GridType& grid = *gridPtr;
>
> if (refinementSteps)
> grid.globalRefine(refinementSteps);
>
> DGStokesParameters parameters;
> Example<dim, NumberType> exactSolution;
> DirichletBoundary<GridType> dirichletBoundary(exactSolution);
> RightHandSide<GridType> rightHandSide(exactSolution);
>
> typedef Dune::DGStokes<GridType, vOrder, pOrder> DGStokes;
> DGStokes dGStokes(grid, exactSolution, parameters, dirichletBoundary, rightHandSide, refinementSteps);
> dGStokes.assembleStokesSystem();
> dGStokes.solveStokesSystem();
>
> std::cout << "L2Error: ";
> for (int i = 0; i < dim+1; i++)
> std::cout << dGStokes.l2errorStokesSystem(i) << ", ";
> std::cout << std::endl;
>
> dGStokes.vtkout("test_stokes", 0);
>
> return 0;
> }
> catch (Dune::Exception &e){
> std::cerr << "Dune reported error: " << e << std::endl;
> }
> catch (...){
> std::cerr << "Unknown exception thrown!" << std::endl;
> }
> }
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune
--
Sreejith Pulloor Kuttanikkad
Interdisciplinary Centre for Scientific Computing
University of Heidelberg
Room:336, IUP, Im Neuenheimer Feld-229
69120 Heidelberg,Germany.
Ph :+49-6221-54-6512 / +49-711-7816-4966
:+49-(0)17624228904(Mob)
http://hal.iwr.uni-heidelberg.de/~sreejith/
-------------------------------------------
More information about the Dune
mailing list