• Home
• About DUNE
  • What is DUNE
  • News
  • Gallery
  • Publications
• Download
  • Core Modules (Stable)
    • Release Notes (v1.1)
    • Release Notes (v1.0)
  • Recent Changes
  • External Modules
  • Dune Logo
• Documentation
  • Installation
  • Using DUNE
  • Class Documentation
• Community
  • People
  • Mailinglists
  • FAQ
• Development
  • Developing DUNE
  • Unstable Modules
    • Download via SVN
    • Class Documentation
  • Bug Tracking System
  • Automated Test Environment

massmatrix.hh

Go to the documentation of this file.

#ifndef DUNE_P1_MASSMATRIX_ASSEMBLER_HH
#define DUNE_P1_MASSMATRIX_ASSEMBLER_HH

#include<iostream>
#include<iomanip>
#include<fstream>
#include<sstream>

#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/common/exceptions.hh>
#include<dune/grid/common/grid.hh>
#include<dune/grid/common/referenceelements.hh>
#include<dune/grid/common/quadraturerules.hh>
#include<dune/istl/operators.hh>

#include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
#include<dune/disc/operators/p1operator.hh>
#include<dune/disc/operators/boundaryconditions.hh>
#include<dune/disc/functions/p1function.hh>

namespace Dune
{
    template<class GridView, class RT, int comp>
    class MassMatrixLocalStiffness 
      : public LinearLocalStiffness<GridView, RT, comp>
  {
        // grid types
      typedef typename GridView::Grid::ctype DT;
        typedef typename GridView::template Codim<0>::Entity Entity;
        typedef typename GridView::template Codim<0>::EntityPointer EEntityPointer;

        // some other sizes
        enum {dim=GridView::dimension};

  public:

      // Number of components for the global assembler to collect
      enum {m=comp};

      // types for matrics, vectors and boundary conditions
      typedef FieldMatrix<RT,comp,comp> MBlockType; // one entry in the stiffness matrix
      typedef FieldVector<RT,comp> VBlockType;   // one entry in the global vectors
      typedef array<BoundaryConditions::Flags,comp> BCBlockType;     // componentwise boundary conditions
      
      MassMatrixLocalStiffness (bool procBoundaryAsDirichlet_=true) {

          procBoundaryAsDirichlet = procBoundaryAsDirichlet_;

        }



      void assemble (const Entity& e, int k=1)
      {
          // extract some important parameters
          GeometryType gt = e.type();
          const typename LagrangeShapeFunctionSetContainer<DT,RT,dim>::value_type& sfs
              = LagrangeShapeFunctions<DT,RT,dim>::general(gt,k);
          
          // clear assemble data
          setcurrentsize(sfs.size());

          for (int i=0; i<sfs.size(); i++) {
              this->b[i] = 0;
              this->bctype[i][0] = BoundaryConditions::neumann;
              for (int j=0; j<sfs.size(); j++) 
                  this->A[i][j] = 0;
          }
          
          // Loop over all quadrature points and assemble matrix and right hand side
          
          int p=dim;
          if (gt.isPrism() || gt.isPyramid())
              p = 2;
          
          p *= k;
          
          for (size_t g=0; g<Dune::QuadratureRules<DT,dim>::rule(gt,p).size(); ++g) {
              
              // pos of integration point
              const Dune::FieldVector<DT,dim>& quadPos = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].position();
              
              // weight of quadrature point
              double weight = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].weight();
              
              // determinant of jacobian
              DT detjac = e.geometry().integrationElement(quadPos);

              RT factor = weight*detjac;

              
              RT v[sfs.size()];
              for(int i=0; i<sfs.size(); i++) 
                  v[i] = sfs[i].evaluateFunction(0,quadPos); 
                
              for(int i=0; i<sfs.size(); i++) 
                  for (int j=0; j<=i; j++ ) 
                      for (int k=0; k<comp; k++)
                          this->A[i][j][k][k] += ( v[i] * v[j] ) * factor;
              
          }

          for (int row=0; row<sfs.size(); row++)
              for (int col=0; col<=row; col++) {
                  
                  // Complete the symmetric matrix by copying the lower left triangular matrix
                  // to the upper right one
                  if (row!=col) {
                      for (int rcomp=0; rcomp<dim; rcomp++)
                          for (int ccomp=0; ccomp<dim; ccomp++)
                              this->A[col][row][ccomp][rcomp] = this->A[row][col][rcomp][ccomp];
                  }
              }
              
#if 0
          // evaluate boundary conditions via intersection iterator
          typedef typename IntersectionIteratorGetter<G,TypeTag>::IntersectionIterator
            IntersectionIterator;
          
          IntersectionIterator endit = IntersectionIteratorGetter<G,TypeTag>::end(e);
          for (IntersectionIterator it = IntersectionIteratorGetter<G,TypeTag>::begin(e); it!=endit; ++it)
                {
                  // if we have a neighbor then we assume there is no boundary
                  // (it might still be an interior boundary ... )
                  if (it.neighbor()) continue;

                  // determine boundary condition type for this face, initialize with processor boundary
                  typename BoundaryConditions::Flags bctypeface = BoundaryConditions::process;

                  // handle face on exterior boundary
                  if (it.boundary())
                        {
                          Dune::GeometryType gtface = it.intersectionSelfLocal().type();
                          for (size_t g=0; g<Dune::QuadratureRules<DT,n-1>::rule(gtface,p).size(); ++g)
                                {
                                  const Dune::FieldVector<DT,n-1>& facelocal = Dune::QuadratureRules<DT,n-1>::rule(gtface,p)[g].position();
                                  FieldVector<DT,dim> local = it.intersectionSelfLocal().global(facelocal);
                                  FieldVector<DT,dim> global = it.intersectionGlobal().global(facelocal);
                                  bctypeface = problem.bctype(global,e,local); // eval bctype

//                                std::cout << "=== Boundary found"
//                                                      << " facenumber=" << it.numberInSelf()
//                                                      << " quadpoint=" << g
//                                                      << " facelocal=" << facelocal
//                                                      << " local=" << local
//                                                      << " global=" << global
//                                                      << std::endl;

                                  if (bctypeface!=BoundaryConditions::neumann) break;

                                  RT J = problem.J(global,e,local);
                                  double weightface = Dune::QuadratureRules<DT,n-1>::rule(gtface,p)[g].weight();
                                  DT detjacface = it.intersectionGlobal().integrationElement(facelocal);
                                  for (int i=0; i<sfs.size(); i++) // loop over test function number
                                      if (bctype[i][0]==BoundaryConditions::neumann)
                                          {
                                                b[i] -= J*sfs[i].evaluateFunction(0,local)*weightface*detjacface;
//                                              std::cout << "Neumann BC found, accumulating" 
//                                                                << -J*sfs[i].evaluateFunction(0,local)*weightface*detjacface 
//                                                                << std::endl;
//                                              std::cout << "J=" << J << " shapef=" << sfs[i].evaluateFunction(0,local)
//                                                                << " weight=" << weightface
//                                                                << " detjac=" << detjacface << std::endl;
                                          }
                                }
                          if (bctypeface==BoundaryConditions::neumann) continue; // was a neumann face, go to next face
                        }

                  // If we are here, then its either an exterior boundary face with Dirichlet condition
                  // or a processor boundary (i.e. neither boundary() nor neighbor() was true)
                  // How processor boundaries are handled depends on the processor boundary mode
                  if (bctypeface==BoundaryConditions::process && procBoundaryAsDirichlet==false) 
                        continue; // then it acts like homogeneous Neumann

                  // now handle exterior or interior Dirichlet boundary
                  for (int i=0; i<sfs.size(); i++) // loop over test function number
                        {
                          if (sfs[i].codim()==0) continue; // skip interior dof
                          if (sfs[i].codim()==1) // handle face dofs
                                {
                                  if (sfs[i].entity()==it.numberInSelf())
                                        {
                                          if (bctype[i][0]<bctypeface)
                                                {
                                                  bctype[i][0] = bctypeface;
                                                  if (bctypeface==BoundaryConditions::process)
                                                        b[i] = 0;
                                                  if (bctypeface==BoundaryConditions::dirichlet)
                                                        {
                                                          Dune::FieldVector<DT,dim> global = e.geometry().global(sfs[i].position());
                                                          b[i] = problem.g(global,e,sfs[i].position());
                                                        }
                                                }
                                        }
                                  continue;
                                }
                          // handle subentities of this face
                          for (int j=0; j<ReferenceElements<DT,dim>::general(gt).size(it.numberInSelf(),1,sfs[i].codim()); j++)
                                if (sfs[i].entity()==ReferenceElements<DT,dim>::general(gt).subEntity(it.numberInSelf(),1,j,sfs[i].codim()))
                                  {
                                        if (bctype[i][0]<bctypeface)
                                          {
                                                bctype[i][0] = bctypeface;
                                                if (bctypeface==BoundaryConditions::process)
                                                  b[i] = 0;
                                                if (bctypeface==BoundaryConditions::dirichlet)
                                                  {
                                                        Dune::FieldVector<DT,dim> global = e.geometry().global(sfs[i].position());
                                                        b[i] = problem.g(global,e,sfs[i].position());
                                                  }
                                          }
                                  }
                        }
                }
#endif
      }

      void assembleBoundaryCondition (const Entity& e, int k=1)
      {
          DUNE_THROW(NotImplemented, "!");
      }

  private:
        // parameters given in constructor
        bool procBoundaryAsDirichlet;

  };

}
#endif

Generated on 6 Jan 2009 with Doxygen (ver 1.5.1) [logfile].