• 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

p1mgtransfer.hh

// $Id: p1mgtransfer.hh 519 2008-09-30 18:52:38Z mblatt $

#ifndef DUNE_P1MGTRANSFER_HH
#define DUNE_P1MGTRANSFER_HH

#include<set>
#include<map>

#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/common/geometrytype.hh>

#include<dune/grid/common/grid.hh>
#include<dune/grid/common/mcmgmapper.hh>
#include<dune/grid/utility/intersectiongetter.hh>

#include<dune/istl/bvector.hh>
#include<dune/istl/bcrsmatrix.hh>

#include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
#include<dune/disc/operators/localstiffness.hh>

namespace Dune
{
  template<class G, class RT, int m=1>
  class P1MGTransfer
  {
  public:
        // export type used to store the matrix
        typedef FieldMatrix<RT,m,m> BlockType; 
        typedef BCRSMatrix<BlockType> RepresentationType;
        typedef typename RepresentationType::RowIterator rowiterator;
        typedef typename RepresentationType::ColIterator coliterator;

        // mapper: one data element per vertex
        template<int dim>
        struct P1Layout
        {
          bool contains (Dune::GeometryType gt)
          {
              return gt.dim() == 0;
          }
        }; 

        typedef typename G::ctype DT;
        typedef typename G::LevelGridView GridView;
        enum {n=G::dimension};
        typedef typename G::template Codim<0>::Entity Entity;
        typedef typename GridView::template Codim<0>::Iterator ElementLevelIterator;
        typedef typename GridView::template Codim<n>::Iterator VertexLevelIterator;
        typedef typename GridView::template Codim<0>::EntityPointer EntityPointer;
        typedef typename GridView::IndexSet IS;
        typedef MultipleCodimMultipleGeomTypeMapper<G,IS,P1Layout> VM;

        typedef std::set<int> IntSet;
        typedef std::map<int,IntSet> Graph;

        P1MGTransfer (const G& grid_, int level_)
          : grid(grid_), coarseview(grid.levelView(level_-1)),
            fineview(grid.levelView(level_)), level(level_)
        {
          // check that coarse grid exists
          if (level==0)
                DUNE_THROW(Exception,"P1MGTransfer: level greater 0 required");

          // allocate vertex mappers for the fine and coarse grid
          VM finemapper(grid,grid.levelIndexSet(level));
          VM coarsemapper(grid,grid.levelIndexSet(level-1));

          // allocate a graph structure for the sparsity pattern (yes we are lazy here)
          Graph graph;

          // allocate a flag vector to handle each vertex exactly once
          std::vector<bool> treated(finemapper.size(),false);


          // build the graph structure from the mesh
          for (ElementLevelIterator it = fineview.template begin<0>(); 
                   it!=fineview.template end<0>(); ++it)
                {
                  // get geometry type of entity 
                  GeometryType gt = it->type();

                  // get father entity
                  const EntityPointer father = it->father();

                  // get geometry type of father
                  GeometryType gtf = father->type();

                  // connect every vertex of the fine grid element with 
                  // with every vertex of the coarse grid element
                  for (int i=0; i<it->template count<n>(); i++)
                        {
                          // get index of fine grid vertex
                          int indexi = finemapper.template map<n>(*it,i);

                          // skip vertices that have already been treated
                          if (treated[indexi]) continue;

                          // get position of fine grid vertex in local coordinate system of father
                          const FieldVector<DT,n>& cpos=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].position();
                          FieldVector<DT,n> pos = it->geometryInFather().global(cpos);

                          // determine interpolation weights
                          for (int j=0; j<father-> template count<n>(); j++)
                                {
                                  // get value of j th basis function in father element 
                                  RT phi = Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1)[j].evaluateFunction(0,pos);

                                  // get global index of that basis function
                                  int indexj = coarsemapper.template map<n>(*father,j);

                                  // insert only if non-zero
                                  if (std::abs(phi)>1E-6) graph[indexi].insert(indexj);
                                }
                          treated[indexi] = true;
                        }
                }

          // compute number of nonzeroes
          int nnz=0;
          for (Graph::iterator i=graph.begin(); i!=graph.end(); ++i)
                nnz += i->second.size();

          // allocate the matrix
          std::cout << grid.comm().rank() << ": making " << finemapper.size() << "x"
                                << coarsemapper.size() << " transfer matrix with " << nnz 
                                << " nonzeroes" << std::endl;
          A = new RepresentationType(finemapper.size(),coarsemapper.size(),nnz,RepresentationType::random);

          // set row sizes
          for (Graph::iterator i=graph.begin(); i!=graph.end(); ++i)
                A->setrowsize(i->first,i->second.size());
          A->endrowsizes();

          // set the structure
          for (Graph::iterator i=graph.begin(); i!=graph.end(); ++i)
                for (IntSet::iterator j=i->second.begin(); j!=i->second.end(); ++j)
                  A->addindex(i->first,*j);
          A->endindices();

          // deallocate the map
          graph.clear();

          // allocate skip flag vectors
          for (int compi=0; compi<m; compi++)
                skipflag[compi] = new std::vector<bool>(finemapper.size(),false);
          for (int compi=0; compi<m; compi++)
                coarseskipflag[compi] = new std::vector<bool>(coarsemapper.size(),false);
        }


        void assemble (LocalStiffness<GridView,RT,m>& loc)
        {
          // allocate vertex mappers for the fine and coarse grid
          VM finemapper(grid,grid.levelIndexSet(level));
          VM coarsemapper(grid,grid.levelIndexSet(level-1));

          // clear data
          (*A) = 0;
          for (int compi=0; compi<m; compi++)
                for (int i=0; i<finemapper.size(); i++)
                  (*skipflag[compi])[i] = false;
          for (int compi=0; compi<m; compi++)
                for (int i=0; i<coarsemapper.size(); i++)
                  (*coarseskipflag[compi])[i] = false;

          // determine coarse grid skip flags !
          for (ElementLevelIterator it = coarseview.template begin<0>(); 
                   it!=coarseview.template end<0>(); ++it)
                {
                  // get geometry type of entity 
                  GeometryType gt = it->type();

                  // assemble boundary conditions for the given element
                  loc.assembleBoundaryCondition(*it);

                  // connect every vertex of the fine grid element with 
                  // with every vertex of the coarse grid element
                  for (int i=0; i<it->template count<n>(); i++)
                        {
                          // get index of fine grid vertex
                          int indexi = coarsemapper.template map<n>(*it,i);

                          // set skipflag if essential boundary condition is encountered
                          for (int compi=0; compi<m; compi++)
                                if (loc.bc(i)[compi]==BoundaryConditions::process || 
                                        loc.bc(i)[compi]==BoundaryConditions::dirichlet)
                                  (*coarseskipflag[compi])[indexi] = true;
                        }
                }

          // allocate a flag vector to handle each vertex exactly once
          std::vector<bool> treated(finemapper.size(),false);

          // assemble the entries; needs to consider Dirichlet boundary conditions later
          for (ElementLevelIterator it = fineview.template begin<0>(); 
                   it!=fineview.template end<0>(); ++it)
                {
                  // get geometry type of entity 
                  GeometryType gt = it->type();

                  // get father entity
                  const EntityPointer father = it->father();

                  // get geometry type of father
                  GeometryType gtf = father->type();

                  // assemble boundary conditions for the given element
                  loc.assembleBoundaryCondition(*it);

                  // connect every vertex of the fine grid element with 
                  // with every vertex of the coarse grid element
                  for (int i=0; i<it->template count<n>(); i++)
                        {
                          // get index of fine grid vertex
                          int indexi = finemapper.template map<n>(*it,i);

                          // set skipflag if essential boundary condition is encountered
                          for (int compi=0; compi<m; compi++)
                                if (loc.bc(i)[compi]==BoundaryConditions::process || 
                                        loc.bc(i)[compi]==BoundaryConditions::dirichlet)
                                  (*skipflag[compi])[indexi] = true;

                          // skip vertices that have already been treated
                          if (treated[indexi]) continue;

                          // get position of fine grid vertex in local coordinate system of father
                          const FieldVector<DT,n>& cpos=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].position();
                          FieldVector<DT,n> pos = it->geometryInFather().global(cpos);

                          // compute and set matrix values
                          for (int j=0; j<father->template count<n>(); ++j)
                                {
                                  // get value of j th basis function in father element 
                                  RT phi = Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1)[j].evaluateFunction(0,pos);

                                  // insert entry if large enough
                                  if (std::abs(phi)>1E-6)
                                        {
                                          // get global index of that basis function
                                          int indexj = coarsemapper.template map<n>(*father,j);

                                          // fill diagonal matrix block; consider coarse skip flags
                                          BlockType D;
                                          for (int compi=0; compi<m; compi++)
                                                {
                                                  double scale=1.0;
                                                  if ((*coarseskipflag[compi])[indexj]) scale=0;
                                                  for (int compj=0; compj<m; compj++)
                                                        if (compi==compj) 
                                                          D[compi][compj] = phi*scale; 
                                                        else 
                                                          D[compi][compj] = 0.0;
                                                }

                                          // set entry
                                          (*A)[indexi][indexj] = D;
                                        }
                                }

                          // mark this vertex as done
                          treated[indexi] = true;
                        }
                }
        }

        // return the skip flag
        bool skipFlag (int compi, int indexi) const
        {
          return (*skipflag[compi])[indexi];
        }

        // return the skip flag
        bool coarseSkipFlag (int compi, int indexi) const
        {
          return (*coarseskipflag[compi])[indexi];
        }

        const RepresentationType& operator* () const
        {
          return *A;
        }

        RepresentationType& operator* ()
        {
          return *A;
        }


        // destructor
        ~P1MGTransfer ()
        {
          delete A;
          for (int compi=0; compi<m; compi++)
                delete skipflag[compi];
          for (int compi=0; compi<m; compi++)
                delete coarseskipflag[compi];
        }

  private:

        // make copy constructor and assignment private
        P1MGTransfer (const P1MGTransfer&) {}
        P1MGTransfer& operator= (const P1MGTransfer&) {}

        // data members
        const G& grid;
        const GridView& fineview;
        const GridView& coarseview;
        int level;
        RepresentationType* A;
        std::vector<bool>* skipflag[m];
        std::vector<bool>* coarseskipflag[m];
  };



}
#endif

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