disc » stokes
l2error.hh

00001 #include<dune/disc/stokes/dgstokes.hh>
00002 
00003 
00004 template <class G,int v_order, int p_order>
00005 double
00006 Dune::DGFiniteElementMethod<G,v_order,p_order>::evaluateL2error(int variable,const ExactSolution<ctype, dim> & exact,
00007                                                                                            const Entity& element,const LocalVectorBlock& xe)const
00008 
00009 {
00010   // stokes system has dim+1 variables (dim velocity comps and 1 pressure)
00011   double error[dim+1];
00012   error[variable]=0.0;
00013   Dune::FieldVector<ctype, dim> qp_loc(0.0);
00014   Dune::FieldVector<ctype, dim> qp_glob(0.0);
00015   Dune::GeometryType gt = element.type();
00016   // #warning fixed quadrature order 
00017   int qord=12;
00018 //G grid;
00019   int eid = grid.levelIndexSet(grid.maxLevel()).index(element);
00020   //std::cout<<"EID:--> "<<eid<<std::endl;
00021   for (int qp=0;qp<Dune::QuadratureRules<ctype,dim>::rule(gt,qord).size();++qp) 
00022         {
00023           qp_loc = Dune::QuadratureRules<ctype,dim>::rule(gt,qord)[qp].position();
00024           
00025           qp_glob =element.geometry().global(qp_loc);
00026           //std::cout<<"qp_loc: "<<qp_loc<<" qp_glob: "<<qp_glob<<std::endl;
00027           double weight = Dune::QuadratureRules<ctype,dim>::rule(gt,qord)[qp].weight();
00028           double detjac = element.geometry().integrationElement(qp_loc);
00029           if (variable<dim)
00030                 {
00031                   error[variable]+=weight*detjac
00032                         *(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe))
00033                         *(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe));
00034                   
00035                   
00036                   //  if(variable==0)
00037 //                      std::cout<<"u:"<<exact.velocity(0,qp_glob)<<" comp: "<<evaluateSolution(0,element,qp_loc,xe)<<std::endl;
00038                 }
00039           if(variable==dim)
00040                 {
00041                            
00042                   std::cout<<"qp: "<<qp_glob<<"  p:"<<exact.pressure(qp_glob)<<"  comp:"<<evaluateSolution(variable,element,qp_loc,xe)<<std::endl;
00043                   error[variable]+=weight*detjac
00044                         *(exact.pressure(qp_glob)-evaluateSolution(variable,element,qp_loc,xe))
00045                         *(exact.pressure(qp_glob)-evaluateSolution(variable,element,qp_loc,xe));
00046                 }
00047            
00048         }
00049   return error[variable];
00050  
00051 }
00052 
00053 
00054 template<class G,int v_order, int p_order>
00055 double Dune::DGStokes<G,v_order,p_order>::l2errorStokesSystem(int variable) const
00056 {
00057   // stokes system has dim+1 variables (dim velocity comps and 1 pressure)
00058   double error[dim+1];
00059   error[variable]= 0.0;
00060    ElementLevelIterator it = grid.template lbegin<0>(level);
00061    ElementLevelIterator itend = grid.template lend<0>(level);
00062   //ElementLeafIterator it = grid.template leafbegin<0>();
00063   //ElementLeafIterator itend = grid.template leafend<0>();
00064   //std::cout<<"level:::"<<level;
00065   for (; it != itend; ++it)
00066         {
00067           int eid = grid.levelIndexSet(level).index(*it);
00068           //std::cout<<" eid: "<<eid<<std::endl;
00069           error[variable]+=dgfem.evaluateL2error(variable,exact,*it,b[eid]);
00070         }
00071   return sqrt(error[variable]);
00072 }
00073 
00074