p1groundwater.hh

00001 // $Id: p1groundwater.hh 517 2008-09-30 14:46:52Z mblatt $
00002 
00003 #ifndef DUNE_P1GROUNDWATER_HH
00004 #define DUNE_P1GROUNDWATER_HH
00005 
00006 #include<map>
00007 #include<iostream>
00008 #include<iomanip>
00009 #include<fstream>
00010 #include<vector>
00011 #include<sstream>
00012 
00013 #include<dune/common/exceptions.hh>
00014 #include<dune/common/geometrytype.hh>
00015 
00016 #include<dune/grid/common/grid.hh>
00017 #include<dune/grid/common/referenceelements.hh>
00018 #include<dune/grid/common/quadraturerules.hh>
00019 #include<dune/grid/utility/intersectiongetter.hh>
00020 
00021 #include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
00022 #include<dune/disc/operators/boundaryconditions.hh>
00023 #include<dune/disc/operators/localstiffness.hh>
00024 
00025 #include"groundwater.hh"
00026 
00034 namespace Dune
00035 {
00046 
00047 
00064   template<class GV, class RT>
00065   class GroundwaterEquationLocalStiffness 
00066     : public LinearLocalStiffness<GV,RT,1>
00067   {
00069     typedef GV GridView;
00071     typedef typename GV::Grid Grid;
00073     typedef typename Grid::ctype DT;
00075     typedef typename GV::template Codim<0>::Entity Entity;
00076 
00077   public:
00078         // define the number of components of your system, this is used outside
00079         // to allocate the correct size of (dense) blocks with a FieldMatrix
00080         enum {n=Grid::dimension};
00081         enum {m=1};
00082         enum {SIZE=LagrangeShapeFunctionSetContainer<DT,RT,n>::maxsize};
00083 
00085         GroundwaterEquationLocalStiffness (const GV& gv,
00086                                            const GroundwaterEquationParameters<Grid,RT>& params,
00087                                                                            bool levelBoundaryAsDirichlet_,
00088                                                                            bool procBoundaryAsDirichlet_=true)
00089           : gridview(&gv), problem(params),levelBoundaryAsDirichlet(levelBoundaryAsDirichlet_),
00090                 procBoundaryAsDirichlet(procBoundaryAsDirichlet_)
00091         {}
00092 
00093     void setGridView(const GV& gv)
00094     {
00095       gridview=&gv;
00096     }
00097     
00099 
00108         void assemble (const Entity& e, int k=1)
00109         {
00110           // extract some important parameters
00111           Dune::GeometryType gt = e.type();
00112           const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
00113                 sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
00114           setcurrentsize(sfs.size());
00115 
00116           // clear assemble data
00117           for (int i=0; i<sfs.size(); i++)
00118                 {
00119                   this->b[i] = 0;
00120                   this->bctype[i][0] = BoundaryConditions::neumann;
00121                   for (int j=0; j<sfs.size(); j++) 
00122                         this->A[i][j] = 0;
00123                 }
00124 
00125           assembleV(e,k);
00126           assembleBC(e,k);
00127         }
00128 
00130 
00137         void assembleBoundaryCondition (const Entity& e, int k=1)
00138         {
00139           // extract some important parameters
00140           Dune::GeometryType gt = e.type();
00141           const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
00142                 sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
00143           setcurrentsize(sfs.size());
00144 
00145           // clear assemble data
00146           for (int i=0; i<sfs.size(); i++)
00147                 {
00148                   this->b[i] = 0;
00149                   this->bctype[i][0] = BoundaryConditions::neumann;
00150                 }
00151 
00152           this->template assembleBC(e,k);
00153         }
00154 
00155   private:
00156 
00157         void assembleV (const Entity& e, int k=1)
00158         {
00159           // extract some important parameters
00160           Dune::GeometryType gt = e.type();
00161           const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
00162                 sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
00163           setcurrentsize(sfs.size());
00164 
00165           // Loop over all quadrature points and assemble matrix and right hand side
00166           int p=2;
00167           if (gt.isSimplex()) p=1;
00168           if (k>1) p=2*(k-1);
00169           for (size_t g=0; g<Dune::QuadratureRules<DT,n>::rule(gt,p).size(); ++g) // run through all quadrature points
00170                 {
00171                   const Dune::FieldVector<DT,n>& 
00172                         local = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].position(); // pos of integration point
00173                   Dune::FieldVector<DT,n> global = e.geometry().global(local);     // ip in global coordinates
00174                   const Dune::FieldMatrix<DT,n,n> 
00175                         jac = e.geometry().jacobianInverseTransposed(local);           // eval jacobian inverse
00176                   const Dune::FieldMatrix<DT,n,n> K = problem.K(global,e,local);   // eval diffusion tensor
00177                   double weight = Dune::QuadratureRules<DT,n>::rule(gt,p)[g].weight();// weight of quadrature point
00178                   DT detjac = e.geometry().integrationElement(local);              // determinant of jacobian
00179                   RT q = problem.q(global,e,local);                                // source term
00180                   RT factor = weight*detjac;
00181 
00182                   // evaluate gradients at Gauss points
00183                   Dune::FieldVector<DT,n> grad[SIZE], temp, gv;
00184                   for (int i=0; i<sfs.size(); i++)
00185                         {
00186                           for (int l=0; l<n; l++) 
00187                                 temp[l] = sfs[i].evaluateDerivative(0,l,local);
00188                           grad[i] = 0;
00189                           jac.umv(temp,grad[i]); // transform gradient to global ooordinates
00190                         }
00191 
00192                   for (int i=0; i<sfs.size(); i++) // loop over test function number
00193                   {
00194                         // rhs
00195                         this->b[i] += q*sfs[i].evaluateFunction(0,local)*factor;
00196 
00197                         // matrix
00198                         gv = 0; K.umv(grad[i],gv); // multiply with diffusion tensor
00199                         this->A[i][i] += (grad[i]*gv)*factor;
00200                         for (int j=0; j<i; j++)
00201                           {
00202                                 RT t = (grad[j]*gv)*factor;
00203                                 this->A[i][j] += t;
00204                                 this->A[j][i] += t;
00205                           }
00206                   }
00207                 }
00208         }
00209 
00210         void assembleBC (const Entity& e, int k=1)
00211         {
00212           // extract some important parameters
00213           Dune::GeometryType gt = e.type();
00214           const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
00215                 sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,k);
00216           setcurrentsize(sfs.size());
00217 
00218           // determine quadrature order
00219           int p=2;
00220           if (gt.isSimplex()) p=1;
00221           if (k>1) p=2*(k-1);
00222 
00223           // evaluate boundary conditions via intersection iterator
00224           
00225           typename GridView::IntersectionIterator endit = gridview->iend(e);
00226           for (typename GridView::IntersectionIterator it = gridview->ibegin(e); 
00227                it!=endit; ++it)
00228                 {
00229                   // if we have a neighbor then we assume there is no boundary (forget interior boundaries)
00230                   // in level assemble treat non-level neighbors as boundary
00231                   if (it.neighbor())
00232                         {
00233                           if (levelBoundaryAsDirichlet && it.outside()->level()==e.level()) 
00234                                 continue;
00235                           if (!levelBoundaryAsDirichlet)
00236                                 continue;
00237                         }
00238 
00239                   // determine boundary condition type for this face, initialize with processor boundary
00240                   typename BoundaryConditions::Flags bctypeface = BoundaryConditions::process;
00241 
00242                   // handle face on exterior boundary, this assumes there are no interior boundaries
00243                   if (it.boundary())
00244                         {
00245                           Dune::GeometryType gtface = it.intersectionSelfLocal().type();
00246                           for (size_t g=0; g<Dune::QuadratureRules<DT,n-1>::rule(gtface,p).size(); ++g)
00247                                 {
00248                                   const Dune::FieldVector<DT,n-1>& facelocal = Dune::QuadratureRules<DT,n-1>::rule(gtface,p)[g].position();
00249                                   FieldVector<DT,n> local = it.intersectionSelfLocal().global(facelocal);
00250                                   FieldVector<DT,n> global = it.intersectionGlobal().global(facelocal);
00251                                   bctypeface = problem.bctype(global,e,local); // eval bctype
00252 
00253 
00254                                   if (bctypeface!=BoundaryConditions::neumann) break;
00255 
00256                                   RT J = problem.J(global,e,local);
00257                                   double weightface = Dune::QuadratureRules<DT,n-1>::rule(gtface,p)[g].weight();
00258                                   DT detjacface = it->intersectionGlobal().integrationElement(facelocal);
00259                                   for (int i=0; i<sfs.size(); i++) // loop over test function number
00260                                         if (this->bctype[i][0]==BoundaryConditions::neumann)
00261                                           {
00262                                                 this->b[i] -= J*sfs[i].evaluateFunction(0,local)*weightface*detjacface;
00263                                           }
00264                                 }
00265                           if (bctypeface==BoundaryConditions::neumann) continue; // was a neumann face, go to next face
00266                         }
00267 
00268                   // If we are here, then it is 
00269                   // (i)   an exterior boundary face with Dirichlet condition, or
00270                   // (ii)  a processor boundary (i.e. neither boundary() nor neighbor() was true), or
00271                   // (iii) a level boundary in case of level-wise assemble
00272                   // How processor boundaries are handled depends on the processor boundary mode
00273                   if (bctypeface==BoundaryConditions::process && procBoundaryAsDirichlet==false 
00274                           && levelBoundaryAsDirichlet==false) 
00275                         continue; // then it acts like homogeneous Neumann
00276 
00277                   // now handle exterior or interior Dirichlet boundary
00278                   for (int i=0; i<sfs.size(); i++) // loop over test function number
00279                         {
00280                           if (sfs[i].codim()==0) continue; // skip interior dof
00281                           if (sfs[i].codim()==1) // handle face dofs
00282                                 {
00283                                   if (sfs[i].entity()==it->numberInSelf())
00284                                         {
00285                                           if (this->bctype[i][0]<bctypeface)
00286                                                 {
00287                                                   this->bctype[i][0] = bctypeface;
00288                                                   if (bctypeface==BoundaryConditions::process)
00289                                                         this->b[i] = 0;
00290                                                   if (bctypeface==BoundaryConditions::dirichlet)
00291                                                         {
00292                                                           Dune::FieldVector<DT,n> global = e.geometry().global(sfs[i].position());
00293                                                           this->b[i] = problem.g(global,e,sfs[i].position());
00294                                                         }
00295                                                 }
00296                                         }
00297                                   continue;
00298                                 }
00299                           // handle subentities of this face
00300                           for (int j=0; j<ReferenceElements<DT,n>::general(gt).size(it->numberInSelf(),1,sfs[i].codim()); j++)
00301                                 if (sfs[i].entity()==ReferenceElements<DT,n>::general(gt).subEntity(it->numberInSelf(),1,j,sfs[i].codim()))
00302                                   {
00303                                         if (this->bctype[i][0]<bctypeface)
00304                                           {
00305                                                 this->bctype[i][0] = bctypeface;
00306                                                 if (bctypeface==BoundaryConditions::process)
00307                                                   this->b[i] = 0;
00308                                                 if (bctypeface==BoundaryConditions::dirichlet)
00309                                                   {
00310                                                         Dune::FieldVector<DT,n> global = e.geometry().global(sfs[i].position());
00311                                                         this->b[i] = problem.g(global,e,sfs[i].position());
00312                                                   }
00313                                           }
00314                                   }
00315                         }
00316                 }
00317         }
00318 
00319         // parameters given in constructor
00320        const GV* gridview;
00321         const GroundwaterEquationParameters<Grid,RT>& problem;
00322         bool levelBoundaryAsDirichlet;
00323         bool procBoundaryAsDirichlet;
00324   };
00325 
00327 }
00328 #endif