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(){}
 };
 
