laplace.hh

00001 #ifndef DUNE_P1_LAPLACE_ASSEMBLER_HH
00002 #define DUNE_P1_LAPLACE_ASSEMBLER_HH
00003 
00004 #include<iostream>
00005 #include<iomanip>
00006 #include<fstream>
00007 #include<sstream>
00008 
00009 #include<dune/common/fvector.hh>
00010 #include<dune/common/fmatrix.hh>
00011 #include<dune/common/exceptions.hh>
00012 #include<dune/grid/common/grid.hh>
00013 #include<dune/grid/common/referenceelements.hh>
00014 
00015 #include<dune/grid/common/quadraturerules.hh>
00016 
00017 #include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
00018 #include<dune/disc/operators/localstiffness.hh>
00019 
00027 namespace Dune
00028 {
00044     template<class GridView, class RT>
00045     class LaplaceLocalStiffness : public LinearLocalStiffness<GridView, RT, 1>
00046   {
00047         // grid types
00048       typedef typename GridView::Grid::ctype DT;
00049         typedef typename GridView::template Codim<0>::Entity Entity;
00050         typedef typename GridView::template Codim<0>::EntityPointer EEntityPointer;
00051 
00052         // some other sizes
00053         enum {dim=GridView::dimension};
00054 
00055   public:
00056 
00058       enum {m=1};
00059 
00060       // types for matrics, vectors and boundary conditions
00061 
00063       typedef FieldMatrix<RT,m,m> MBlockType;
00064 
00066       typedef FieldVector<RT,m> VBlockType;
00067 
00069       typedef array<BoundaryConditions::Flags,m> BCBlockType;
00070       
00072       LaplaceLocalStiffness (bool procBoundaryAsDirichlet_=true) {
00073 
00074           procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
00075 
00076         }
00077 
00078 
00080 
00089       void assemble (const Entity& e, int k=1)
00090       {
00091           // extract some important parameters
00092           GeometryType gt = e.type();
00093           const typename LagrangeShapeFunctionSetContainer<DT,RT,dim>::value_type& sfs
00094               = LagrangeShapeFunctions<DT,RT,dim>::general(gt,k);
00095           
00096           // clear assemble data
00097           setcurrentsize(sfs.size());
00098 
00099           for (int i=0; i<sfs.size(); i++) {
00100               this->b[i] = 0;
00101               this->bctype[i][0] = BoundaryConditions::neumann;
00102               for (int j=0; j<sfs.size(); j++) 
00103                   this->A[i][j] = 0;
00104           }
00105 
00106           // Loop over all quadrature points and assemble matrix and right hand side
00107           
00109           int p=dim*k;
00110           
00111           for (size_t g=0; g<Dune::QuadratureRules<DT,dim>::rule(gt,p).size(); ++g) {
00112               
00113               // pos of integration point
00114               const Dune::FieldVector<DT,dim>& quadPos = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].position();
00115               
00116               // weight of quadrature point
00117               double weight = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].weight();
00118               
00119               const FieldMatrix<DT,dim,dim>& inv = e.geometry().jacobianInverseTransposed(quadPos);
00120 
00121               // determinant of jacobian
00122               DT detjac = e.geometry().integrationElement(quadPos);
00123 
00124               RT factor = weight*detjac;
00125 
00126               // evaluate gradients at Gauss points
00127               Dune::FieldVector<DT,dim> grad[sfs.size()], temp;
00128               for (int i=0; i<sfs.size(); i++)
00129                   {
00130                       for (int l=0; l<dim; l++)
00131                           temp[l] = sfs[i].evaluateDerivative(0,l,quadPos);
00132                       grad[i] = 0;
00133                       inv.umv(temp,grad[i]); // transform gradient to global ooordinates
00134                   }
00135               
00136               // loop over test function number
00137               for (int i=0; i<sfs.size(); i++) {
00138 
00139                   // matrix
00140                   this->A[i][i] += (grad[i]*grad[i])*factor;
00141                   for (int j=0; j<i; j++)
00142                       {
00143                           RT t = (grad[j]*grad[i])*factor;
00144                           this->A[i][j] += t;
00145                           this->A[j][i] += t;
00146                       }
00147               }
00148           }
00149 
00150       }
00151 
00153         void assembleBoundaryCondition (const Entity& e, int k=1)
00154         {
00155         }
00156 
00157   private:
00158         // parameters given in constructor
00159         bool procBoundaryAsDirichlet;
00160 
00161   };
00162 
00164 }
00165 #endif