disc » stokes
dgstokes.hh

00001 
00002 #ifndef DUNE_DGSTOKES_HH
00003 #define DUNE_DGSTOKES_HH
00004 
00005 #include<dune/common/fvector.hh>
00006 #include<dune/common/fmatrix.hh>
00007 #include<dune/common/exceptions.hh>
00008 #include<dune/grid/common/grid.hh>
00009 #include<dune/grid/common/referenceelements.hh>
00010 #include<dune/istl/operators.hh>
00011 #include<dune/istl/bvector.hh>
00012 #include<dune/istl/bcrsmatrix.hh>
00013 #include<dune/istl/io.hh>
00014 
00015 #include<dune/grid/common/quadraturerules.hh>
00016 
00017 #include<dune/disc/shapefunctions/dgspace/monomialshapefunctions.hh>
00018 //#include <dune/disc/functions/dgstokesfunction.hh>
00019 //#include <dune/disc/functions/dgfunction.hh>
00020 
00021 #include"stokesparameters.hh"
00022 #include"boundaryconditions.hh"
00023 #include"testfunctions.hh"
00024 #include"rhs.hh"
00025 
00026 #include <dune/istl/solvers.hh>
00027 
00028 namespace Dune
00029 {
00030 template<class G,int v_order, int p_order>
00031   class DGFiniteElementMethod
00032   {
00033         //dimension of grid
00034         enum {dim=G::dimension};
00035         enum { dimw=G::dimensionworld };
00036   public:
00037         //Grid
00038         typedef G Grid;
00039         //coordinate type
00040     typedef typename Grid::ctype ctype;
00041         typedef Dune::FieldVector< double , dim> Gradient;
00042         typedef Dune::FieldMatrix< double , dim, dim> InverseJacobianMatrix;
00043         // "order" is order of velocity shapefn
00044         // order of pressure shapefn  is (order-1) i.e, assumed one order less than that of velocity
00045         enum{v_ordr = v_order};
00046         enum{p_ordr=p_order};
00047     
00048         //local vector and matrix blocks
00049         // local block size is sum of velocity dof and pressure dof
00050         //block size = dim*vdof.size() + pdof.size()
00051         static const int BlockSize =dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize+Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize;
00052 static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
00053         static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
00054         
00055         typedef Dune::FieldVector<double,BlockSize> LocalVectorBlock;
00056         typedef Dune::FieldMatrix<double,BlockSize,BlockSize> LocalMatrixBlock;
00057         //shapefn
00058         typedef Dune::MonomialShapeFunctionSet<ctype,double,dim> ShapeFunctionSet;
00059         inline const ShapeFunctionSet & getVelocityShapeFunctionSet(Dune::GeometryType gt) const;
00060         inline const ShapeFunctionSet & getPressureShapeFunctionSet(Dune::GeometryType gt) const;
00061         
00062         typedef typename Grid::template Codim<0>::LeafIterator ElementLeafIterator;
00063         typedef typename Grid::template Codim<0>::LevelIterator ElementLevelIterator;
00064         typedef typename Grid::template Codim<0>::EntityPointer EntityPointer;
00065         typedef typename Grid::template Codim<0>::Entity Entity;
00066         typedef typename Grid::template Codim<dim>::LevelIterator VertexIterator;
00067         typedef typename Grid::template Codim<dim>::Entity Vertex;
00068         typedef typename Grid::template Codim<0>::LevelIntersectionIterator IntersectionLevelIterator;
00069         typedef typename Grid::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
00070         typedef typename Grid::template Codim<1>::EntityPointer InterSectionPointer;
00071 
00072 
00073 
00074         
00075         DGFiniteElementMethod (Grid &g,const DGStokesParameters& par, DirichletBoundary<G>& db, RightHandSide<G>& rh): grid(g),parameter(par),dirichletvalue(db), rhsvalue(rh) {};
00076         //local assembly
00077         void assembleVolumeTerm(Entity& ep, LocalMatrixBlock& Aee,LocalVectorBlock& Be) const;
00078         void assembleFaceTerm(Entity& ep,IntersectionIterator& isp, LocalMatrixBlock& Aee,
00079                                                   LocalMatrixBlock& Aef,LocalMatrixBlock& Afe, LocalVectorBlock& Be) const;
00080         
00081         void assembleBoundaryTerm(Entity& ep, IntersectionIterator& isp, LocalMatrixBlock& Aee,LocalVectorBlock& Be)const ;
00082         // stokes system has dim+1 variables (dim velocity comps and 1 pressure)
00083         double evaluateSolution(int variable,const Entity& element,const Dune::FieldVector<ctype,dim>& local,
00084                                                         const LocalVectorBlock& xe) const;
00085         Gradient evaluateGradient(int variable,const Entity& element,const Dune::FieldVector<ctype,dim>& local,
00086                                                         const LocalVectorBlock& xe) const;
00087         double evaluateL2error(int variable,const ExactSolution<ctype, dim> & exact,const Entity& element,
00088                                                    const LocalVectorBlock& xe)const;
00089 
00090         double evaluateH1error(int variable,const ExactSolution<ctype, dim> & exact,const Entity& element,
00091                                                    const LocalVectorBlock& xe)const;
00092         
00093         
00094                 
00095   private:
00096         Grid& grid;
00097         Dune::MonomialShapeFunctionSetContainer<ctype,double,dim,v_order> vspace;
00098         Dune::MonomialShapeFunctionSetContainer<ctype,double,dim,(p_order)> pspace;
00099 //      ShapeFunctionSet vsfs(v_order); //for  velocity
00100 //      ShapeFunctionSet psfs(p_order); // for pressure
00101         DGStokesParameters parameter;
00102         DirichletBoundary<G> dirichletvalue;
00103         RightHandSide<G> rhsvalue;
00104 
00105   };
00106 
00107 
00108 
00109 template<class G,int v_order,int p_order>
00110   class DGStokes
00111   {
00112         
00113   public:
00114         //dimension of grid
00115         enum {dimension=G::dimension};
00116         enum { dimensionworld=G::dimensionworld };
00117         typedef G Grid;
00118   private:
00119         enum {dim=G::dimension};
00120         enum { dimw=G::dimensionworld };
00121     typedef typename Grid::ctype ctype;
00122         // Iterators
00123         typedef typename Grid::template Codim<0>::LeafIterator ElementLeafIterator;
00124         typedef typename Grid::template Codim<0>::LevelIterator ElementLevelIterator;
00125         typedef typename Grid::template Codim<0>::EntityPointer EntityPointer;
00126         typedef typename Grid::template Codim<0>::Entity Entity;
00127         typedef typename Grid::template Codim<dim>::Entity Vertex;
00128         typedef typename Grid::template Codim<dim>::LevelIterator VertexIterator;
00129         typedef typename Grid::template Codim<0>::LevelIntersectionIterator IntersectionLevelIterator;
00130         typedef typename Grid::template Codim<0>::LeafIntersectionIterator IntersectionIterator;
00131         typedef typename Grid::template Codim<1>::EntityPointer InterSectionPointer;    
00132         static const int BlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize)+(Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
00133 
00134 static const int VBlockSize =((dim*Dune::MonomialShapeFunctionSetSize<dim,v_order>::maxSize));
00135         static const int PBlockSize =((Dune::MonomialShapeFunctionSetSize<dim,p_order>::maxSize));
00136         
00137         typedef typename DGFiniteElementMethod<G,v_order,p_order>::Gradient Gradient;
00138         typedef typename DGFiniteElementMethod<G,v_order,p_order>::LocalVectorBlock LocalVectorBlock;
00139         typedef typename DGFiniteElementMethod<G,v_order,p_order>::LocalMatrixBlock LocalMatrixBlock;
00140         typedef Dune::BlockVector<LocalVectorBlock> LVector;
00141         typedef Dune::BCRSMatrix<LocalMatrixBlock> LMatrix;
00142         typedef Dune::LevelDGStokesVelocityFunction<Grid, double, v_order> DGSVFunction;
00143         typedef Dune::LevelDGStokesPressureFunction<Grid, double, p_order> DGSPFunction;
00144         
00145         //typedef Dune::LevelDGFunction<Grid, double, v_order> DGFunction;
00146   public:
00147         //enum {order = ordr};
00148         //inline constructor with initializer list
00149         
00150         DGStokes(Grid &g, ExactSolution<ctype,dim>& ex,
00151                          DGStokesParameters& par, DirichletBoundary<Grid>& db,
00152                          RightHandSide<Grid>& rh, int l ) : grid(g),exact(ex), dgfem(g,par,db,rh),level(l),xv(grid,level),xp(grid,level)
00153         {
00154           if (par.sigma==1 & par.epsilon==1)
00155         std::cout<<"You are using NIPG scheme"<<std::endl;
00156   else if(par.sigma==1 & par.epsilon==-1)
00157         std::cout<<"You are using SIPG scheme"<<std::endl;
00158   else if(par.sigma==0 & par.epsilon==1)
00159         std::cout<<"You are using OBB scheme"<<std::endl;
00160   else
00161         std::cout<<"check DG parameters epsilon and sigma"<<std::endl;
00162         };
00163 
00164         // global assembly and solving
00165         void assembleStokesSystem() ;
00166         void solveStokesSystem();
00167         double evaluateSolution(const EntityPointer & e,
00168                             const Dune::FieldVector<ctype, dim> & local) const;
00169         
00170         //l2error computation
00171         // stokes system has dim+1 variables (dim velocity comps and 1 pressure)
00172         double l2errorStokesSystem(int variable) const;
00173         double h1errorStokesSystem(int variable) const;
00174         const DGSVFunction & velocity() const
00175         {
00176           return xv;
00177         }
00178                 const DGSPFunction & pressure() const
00179         {
00180           return xp;
00181         }
00182         
00183 
00184   private:
00185         typedef typename DGFiniteElementMethod<G,v_order,p_order>::ShapeFunctionSet ShapeFunctionSet;
00186         inline const ShapeFunctionSet & getVelocityShapeFunctionSet(const EntityPointer &ep) const;
00187         inline const ShapeFunctionSet & getPressureShapeFunctionSet(const EntityPointer &ep) const;
00188   public:
00189         Grid & grid;
00190         int level;
00191         ExactSolution<ctype,dim>& exact;
00192         //for  velocity
00193   // ShapeFunctionSet& vsfs(ordr);
00194 //   // for pressure
00195 //   ShapeFunctionSet& psfs(ordr-1); 
00196   private:
00197         DGFiniteElementMethod<G,v_order,p_order> dgfem;
00198         LMatrix A;
00199     LVector b;
00200         LVector solution;
00201         DGSVFunction xv;
00202         DGSPFunction xp;
00203         //DGFunction x;
00204 
00205         
00206         
00207   };
00208 
00209    
00210 
00211 #include "dgstokes.cc"
00212 
00213 } // end of namespace
00214 
00215 #endif