disc » stokes
h1error.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>::evaluateH1error(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],error1[dim+1],error2[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                   error1[variable]=(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe))
00032                         *(exact.velocity(variable,qp_glob)-evaluateSolution(variable,element,qp_loc,xe));
00033 
00034                   error2[variable]=(exact.velocityGradient(variable,qp_glob)-evaluateGradient(variable,element,qp_loc,xe)).two_norm2();
00035                   error[variable]+=weight*detjac*(error1[variable]+error2[variable]);
00036                   
00037                   
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>::h1errorStokesSystem(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   //std::cout<<"level:::"<<level;
00063   for (; it != itend; ++it)
00064         {
00065           int eid = grid.levelIndexSet(level).index(*it);
00066           //std::cout<<" eid: "<<eid<<std::endl;
00067           error[variable]+=dgfem.evaluateH1error(variable,exact,*it,b[eid]);
00068         }
00069   return sqrt(error[variable]);
00070 }
00071 
00072