• 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

p1function.hh

Go to the documentation of this file.

// $Id: p1function.hh 523 2008-11-17 17:52:49Z christi $

#ifndef DUNE_P1FUNCTION_HH
#define DUNE_P1FUNCTION_HH

//C++ includes
#include<new>
#include<iostream>
#include<vector>
#include<list>
#include<map>
#include<set>

// Dune includes
#include<dune/common/fvector.hh>
#include<dune/common/exceptions.hh>
#include<dune/common/tuples.hh>
#include<dune/common/stdstreams.hh>
#include<dune/grid/common/grid.hh>
#include<dune/grid/common/mcmgmapper.hh>
#include<dune/grid/common/universalmapper.hh>
#include<dune/grid/io/file/vtk/vtkwriter.hh>
#include<dune/istl/bvector.hh>
#include<dune/istl/operators.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/owneroverlapcopy.hh>
#include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>

// same directory includes
#include"functions.hh"
#include"p0function.hh"

namespace Dune
{
  template<class G>
  class LeafCommunicate
  {
  public:
        LeafCommunicate (const G& g)
          : grid(g)
        {}

        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
        {
          grid.template communicate<DataHandle>(data,iftype,dir);
        }

  private:
        const G& grid;
  };

  template<class G>
  class LevelCommunicate
  {
  public:
        LevelCommunicate (const G& g, int l)
          : grid(g), level(l)
        {}

        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
        {
          grid.template communicate<DataHandle>(data,iftype,dir,level);
        }

  private:
        const G& grid;
        int level;
  };

  template<class G, class GV, class VM, class LC>
  class P1ExtendOverlap {

        // types
        typedef typename G::ctype DT;
        enum {n=G::dimension};
      typedef typename GV::IndexSet IS;
        typedef typename G::template Codim<0>::Entity Entity;
        typedef typename GV::template Codim<0>::Iterator Iterator;
        typedef typename GV::template Codim<n>::Iterator VIterator;
        typedef typename G::template Codim<0>::EntityPointer EEntityPointer;
        typedef typename G::Traits::GlobalIdSet IDS;
        typedef typename IDS::IdType IdType;
        typedef std::set<IdType> GIDSet;
        typedef std::pair<IdType,int> Pair;
        typedef std::set<int> ProcSet;

        // A DataHandle class to exchange border rows
        class IdExchange 
          : public CommDataHandleIF<IdExchange,Pair> {
        public:
          typedef Pair DataType;

          bool contains (int dim, int codim) const
          {
                return (codim==dim);
          }

          bool fixedsize (int dim, int codim) const
          {
                return false;
          }

          template<class EntityType>
          size_t size (EntityType& e) const
          {
                return myids[vertexmapper.map(e)].size();
          }

          template<class MessageBuffer, class EntityType>
          void gather (MessageBuffer& buff, const EntityType& e) const
          {
                int alpha=vertexmapper.map(e);
                GIDSet& thisset = myids[alpha];
                for (typename GIDSet::iterator i=thisset.begin(); i!=thisset.end(); ++i)
                  {
                        buff.write(Pair(*i,grid.comm().rank())); // I have these global ids
                        owner[*i] = grid.comm().rank();
                  }
                myprocs[alpha].insert(grid.comm().rank());
          }

          template<class MessageBuffer, class EntityType>
          void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
          {
                int alpha=vertexmapper.map(e);
                GIDSet& thisset = myids[alpha];
                int source;
                for (size_t i=0; i<n; i++)
                  {
                        Pair x;
                        buff.read(x);
                        thisset.insert(x.first);
                        source=x.second;
                        if (owner.find(x.first)==owner.end())
                          owner[x.first] = source;
                        else
                          owner[x.first] = std::min(owner[x.first],source);
                  }
                myprocs[alpha].insert(source);          
          }

          IdExchange (const G& g, const VM& vm, std::map<int,GIDSet>& ids, std::map<int,ProcSet>& procs, 
                                  std::map<IdType,int>& o) 
                : grid(g), vertexmapper(vm), myids(ids), myprocs(procs), owner(o)
          {}
 
        private:
          const G& grid;
          const VM& vertexmapper;
          std::map<int,GIDSet>& myids;
          std::map<int,ProcSet>& myprocs;
          std::map<IdType,int>& owner;
        };

        // A DataHandle class to exchange border rows
        class BorderLinksExchange 
          : public CommDataHandleIF<BorderLinksExchange,IdType>{
        public:
          typedef IdType DataType;

          bool contains (int dim, int codim) const
          {
                return (codim==dim);
          }

          bool fixedsize (int dim, int codim) const
          {
                return false;
          }

          template<class EntityType>
          size_t size (EntityType& e) const
          {
                return borderlinks[vertexmapper.map(e)].size();
          }

          template<class MessageBuffer, class EntityType>
          void gather (MessageBuffer& buff, const EntityType& e) const
          {
                GIDSet& myset = borderlinks[vertexmapper.map(e)];
                for (typename GIDSet::iterator i=myset.begin(); i!=myset.end(); ++i)
                  buff.write(*i); 
          }

          template<class MessageBuffer, class EntityType>
          void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
          {
                GIDSet& myset = borderlinks[vertexmapper.map(e)];
                for (size_t i=0; i<n; i++)
                  {
                        DataType x;
                        buff.read(x);
                        myset.insert(x);
                  }
          }

          BorderLinksExchange (const G& g, std::map<int,GIDSet>& bl, const VM& vm) 
                : grid(g), borderlinks(bl), vertexmapper(vm) 
          {}
 
        private:
          const G& grid;        
          std::map<int,GIDSet>& borderlinks;
          const VM& vertexmapper;
        };

  public:

        enum Attributes {slave=OwnerOverlapCopyAttributeSet::copy, 
                                         master=OwnerOverlapCopyAttributeSet::owner, 
                                         overlap=OwnerOverlapCopyAttributeSet::overlap};

        typedef IndexInfoFromGrid<IdType,int> P1IndexInfoFromGrid;

        void fillIndexInfoFromGrid (const G& grid, const GV& gridview, const VM& vertexmapper, P1IndexInfoFromGrid& info)
        {
          // build a map of sets where each local index is assigned 
          // a set of global ids which are neighbors of this vertex
          // At the same time assign to each local index to a set of processors
          // and at the same time determine the owner of the gid
          std::map<int,GIDSet> myids;
          std::map<int,ProcSet> myprocs;
          std::map<IdType,int> owner;
          Iterator eendit = gridview.template end<0>();
          for (Iterator it = gridview.template begin<0>(); it!=eendit; ++it)
                {
                  Dune::GeometryType gt = it->type();
                  const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
                        refelem = ReferenceElements<DT,n>::general(gt);

                  if (it->partitionType()==InteriorEntity)
                        for (int i=0; i<refelem.size(n); i++)
                          {
                                if (it->template entity<n>(i)->partitionType()==BorderEntity)
                                  {
                                        int alpha = vertexmapper.template map<n>(*it,i);
                                        GIDSet& thisset = myids[alpha];
                                        for (int j=0; j<refelem.size(n); j++)
                                          {
                                                IdType beta = grid.globalIdSet().template subId<n>(*it,j);
                                                thisset.insert(beta);
                                          }
                                  }
                          }
                }
          IdExchange datahandle(grid,vertexmapper,myids,myprocs,owner);
          lc.template communicate<IdExchange>(datahandle,InteriorBorder_InteriorBorder_Interface,ForwardCommunication);

          // build map from global id to local index
          std::map<IdType,int> gid2index;
          for (typename std::map<int,GIDSet>::iterator i=myids.begin(); i!=myids.end(); ++i)
                for (typename GIDSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
                  gid2index[*j] = -1; // indicates "not assigned yet"
          VIterator vendit = gridview.template end<n>();
          for (VIterator it = gridview.template begin<n>(); it!=vendit; ++it)
                {
                  IdType beta = grid.globalIdSet().id(*it);
                  if (gid2index.find(beta)!=gid2index.end())
                        {
                          int alpha = vertexmapper.map(*it);
                          gid2index[beta] = alpha; // assign existing local index
                        }
                }
          int extraDOFs = 0;
          for (typename std::map<IdType,int>::iterator i=gid2index.begin(); i!=gid2index.end(); ++i)
                if (i->second==-1)
                  {
                        i->second = vertexmapper.size()+extraDOFs; // assign new local index
                        extraDOFs++;
                  }

          // build a set of all neighboring processors
          ProcSet neighbors;
          for (typename std::map<int,ProcSet>::iterator i=myprocs.begin(); i!=myprocs.end(); ++i)
                for (typename ProcSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
                  if (*j!=grid.comm().rank())
                        neighbors.insert(*j);

          // now all the necessary information is in place

          // application: for all neighbors build a list of global ids
          for (typename ProcSet::iterator p=neighbors.begin(); p!=neighbors.end(); ++p)
                {
                  GIDSet remote;
                  for (typename std::map<int,ProcSet>::iterator i=myprocs.begin(); i!=myprocs.end(); ++i)
                        if ((i->second).find(*p)!=(i->second).end())
                          {
                                GIDSet& thisset = myids[i->first];
                                for (typename GIDSet::iterator j=thisset.begin(); j!=thisset.end(); ++j)
                                  remote.insert(*j);
                          }
                }


          // fill the info object
          std::set< Tuple<IdType,int,int> > ownindices;
          for (typename std::map<int,GIDSet>::iterator i=myids.begin(); i!=myids.end(); ++i)
                for (typename GIDSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
                  {
                        int a=slave;
                        if (owner[*j]==grid.comm().rank()) a=master;
                        info.addLocalIndex(Tuple<IdType,int,int>(*j,gid2index[*j],a));
                  }
          std::set< Tuple<int,IdType,int> > remoteindices;
          for (typename std::map<int,ProcSet>::iterator i=myprocs.begin(); i!=myprocs.end(); ++i)
                {
                  GIDSet& thisset = myids[i->first];
                  for (typename GIDSet::iterator j=thisset.begin(); j!=thisset.end(); ++j)
                        for (typename ProcSet::iterator p=(i->second).begin(); p!=(i->second).end(); ++p)
                          {
                                int a=slave;
                                if (owner[*j]==(*p)) a=master;
                                if (*p!=grid.comm().rank()) info.addRemoteIndex(Tuple<int,IdType,int>(*p,*j,a));
                          }
                }

          // clear what is not needed anymore to save memory
          myids.clear();
          gid2index.clear();
          myprocs.clear();
          owner.clear();
          neighbors.clear();

          return;
        }


        void extend (const G& grid, const GV& gridView, const VM& vertexmapper,
                                 std::map<int,GIDSet>& borderlinks, int& extraDOFs, std::map<IdType,int>& gid2index)
        {
          // initialize output parameters
          borderlinks.clear();
          extraDOFs = 0;
          gid2index.clear();

          // build local borderlinks from mesh
          Iterator eendit = gridView.template end<0>();
          for (Iterator it = gridView.template begin<0>(); it!=eendit; ++it)
                {
                  Dune::GeometryType gt = it->type();
                  const typename Dune::ReferenceElementContainer<DT,n>::value_type& 
                        refelem = ReferenceElements<DT,n>::general(gt);

                  // generate set of neighbors in global ids for border vertices
                  if (it->partitionType()==InteriorEntity)
                        for (int i=0; i<refelem.size(n); i++)
                          if (it->template entity<n>(i)->partitionType()==BorderEntity)
                                {
                                  int alpha = vertexmapper.template map<n>(*it,i);
                                  GIDSet& myset = borderlinks[alpha];
                                  for (int j=0; j<refelem.size(n); j++)
                                        if (i!=j)
                                          {
                                                IdType beta = grid.globalIdSet().template subId<n>(*it,j);
                                                myset.insert(beta);
                                                //                                                std::cout << g.comm().rank() << ": "
                                                //                                                                      << "borderlink " << alpha
                                                //                                                                      << " " << vertexmapper.template map<n>(*it,j)
                                                //                                                                      << " " << beta
                                                //                                                                      << std::endl;
                                          }
                                }
                }

          // exchange neighbor info for border vertices
          BorderLinksExchange datahandle(grid,borderlinks,vertexmapper);
          lc.template communicate<BorderLinksExchange>(datahandle,
                                                                                                         InteriorBorder_InteriorBorder_Interface,
                                                                                                         ForwardCommunication);

          // initialize inverse map with ids we have
          for (typename std::map<int,GIDSet>::iterator i=borderlinks.begin(); i!=borderlinks.end(); ++i)
                for (typename GIDSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
                  gid2index[*j] = -1;

          // check with ids we already have in the grid to find out extra vertices
          VIterator vendit = gridView.template end<n>();
          for (VIterator it = gridView.template begin<n>(); it!=gridView.template end<n>(); ++it)
                {
                  IdType beta = grid.globalIdSet().id(*it);
                  if (gid2index.find(beta)!=gid2index.end())
                        {
                          int alpha = vertexmapper.map(*it);
                          gid2index[beta] = alpha;
                        }
                }

          // assign index to extra DOFs
          extraDOFs = 0;
          for (typename std::map<IdType,int>::iterator i=gid2index.begin(); i!=gid2index.end(); ++i)
                if (i->second==-1)
                  {
                        i->second = vertexmapper.size()+extraDOFs;
                        extraDOFs++;
                  }

          //              for (typename std::map<int,GIDSet>::iterator i=borderlinks.begin(); i!=borderlinks.end(); ++i)
          //                    for (typename GIDSet::iterator j=(i->second).begin(); j!=(i->second).end(); ++j)
          //                      std::cout << grid.comm().rank() << ": " << "comm borderlink " << i->first
          //                                            << " " << gid2index[*j] << " " << *j << std::endl;
        }

        P1ExtendOverlap (LC lcomm)
          : lc(lcomm)
        {}

  private:
        LC lc;
  };

  // forward declaration
  template<class G, class RT> class P1FunctionManager;



  template<class GV, class RT, class LC, int m=1>
  class P1Function
    : virtual public GridFunctionGlobalEvalDefault<GV,RT,m>
    , virtual public FunctionDefault<typename GV::Grid::ctype,RT,GV::Grid::dimension,m>
    , virtual public ElementwiseCInfinityFunction<typename GV::Grid,RT,m>
    , virtual public H1Function<typename GV::Grid::ctype,RT,GV::Grid::dimension,m>
    , virtual public C0GridFunction<typename GV::Grid,RT,m>
  {
        typedef typename GV::Grid G;

        typedef typename G::ctype DT;

        enum {n=G::dimension};

        typedef typename G::template Codim<0>::Entity Entity;

        template<int dim>
        struct P1Layout
        {
          bool contains (Dune::GeometryType gt)
          {
              return gt.dim() == 0;
          }
        }; 

        P1Function (const P1Function&);

        // types
        typedef typename G::Traits::GlobalIdSet IDS;
        typedef typename IDS::IdType IdType;
        typedef std::set<IdType> GIDSet;
    typedef typename GV::IndexSet IS;

  public:
    typedef GV GridView;
    typedef IS IndexSet;
    typedef G Grid;
        typedef FieldVector<RT,m> BlockType;
        typedef BlockVector<BlockType> RepresentationType;
    typedef MultipleCodimMultipleGeomTypeMapper<G,IS,P1Layout> VM;
        typedef typename P1ExtendOverlap<G,GV,VM,LC>::P1IndexInfoFromGrid P1IndexInfoFromGrid;

        P1Function (const GV& gridView, LC lcomm, bool extendoverlap=false) 
          : GridFunctionGlobalEvalDefault<GV,RT,m>(gridView)
          , gv(gridView), mapper_(gridView.grid(),gridView.indexSet()), lc(lcomm), oldcoeff(0)
        {
          // check if overlap extension is possible
          if (extendoverlap && gv.grid().overlapSize(0)>0)
                DUNE_THROW(GridError,"P1Function: extending overlap requires nonoverlapping grid");

          // no extra DOFs so far
          extraDOFs = 0;
          extendOverlap = extendoverlap;

          // overlap extension
          if (extendoverlap)
                {
                  // set of neighbors in global ids for border vertices
                  std::map<int,GIDSet> borderlinks;
                  std::map<IdType,int> gid2index;

                  // compute extension
                  P1ExtendOverlap<G,GV,VM,LC> extender(lc);
                  extender.extend(gv.grid(),gv,mapper_,borderlinks,extraDOFs,gid2index);
                }

          // allocate the vector
          oldcoeff = 0;
          try {
                coeff = new RepresentationType(mapper_.size()+extraDOFs);
          }
          catch (std::bad_alloc) {
                std::cerr << "not enough memory in P1Function" << std::endl;
                throw; // rethrow exception
          }
          dverb << "making FE function with " << mapper_.size()+extraDOFs << " components"
                        << "(" << extraDOFs << " extra degrees of freedom)" << std::endl;
        }

        ~P1Function ()
        {
          delete coeff;
          if (oldcoeff!=0) delete oldcoeff;
        }


        virtual RT derivative (int comp, const Dune::FieldVector<int,n>& d, const Dune::FieldVector<DT,n>& x) const
        {
          DUNE_THROW(NotImplemented, "global derivative not implemented yet");
        }


        virtual int order () const
        {
          return 1; // up to now only one derivative is implemented
        }


        virtual RT evallocal (int comp, const Entity& e, const Dune::FieldVector<DT,n>& xi) const
        {
          RT value=0;
          Dune::GeometryType gt = e.type(); // extract type of element
          for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1).size(); ++i)
                value += Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].evaluateFunction(0,xi)*(*coeff)[mapper_.template map<n>(e,i)][comp];
          return value;
        }


        virtual void evalalllocal (const Entity& e, const Dune::FieldVector<DT,G::dimension>& xi, 
                                                           Dune::FieldVector<RT,m>& y) const
        {
          Dune::GeometryType gt = e.type(); // extract type of element
          y = 0;
          for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1).size(); ++i)
                {
                  RT basefuncvalue=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].evaluateFunction(0,xi);
                  int index = mapper_.template map<n>(e,i);
                  for (int c=0; c<m; c++)
                        y[c] += basefuncvalue * (*coeff)[index][c];
                }
        }


        virtual RT derivativelocal (int comp, const Dune::FieldVector<int,n>& d, 
                                                                const Entity& e, const Dune::FieldVector<DT,n>& xi) const
        {
          int dir=-1; 
          int order=0;
          for (int i=0; i<n; i++)
                {
                  order += d[i];
                  if (d[i]>0) dir=i;
                }
      assert(dir != -1);
          if (order!=1) DUNE_THROW(GridError,"can only evaluate one derivative");

          RT value=0;
          Dune::GeometryType gt = e.type(); // extract type of element
          const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,n>::value_type& 
                sfs=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1);
          Dune::FieldMatrix<DT,n,n> jac = e.geometry().jacobianInverseTransposed(xi);
          for (int i=0; i<sfs.size(); ++i)
                {
                  Dune::FieldVector<DT,n> grad(0),temp;
                  for (int l=0; l<n; l++) 
                        temp[l] = sfs[i].evaluateDerivative(0,l,xi);
                  jac.umv(temp,grad); // transform gradient to global ooordinates
                  value += grad[dir] * (*coeff)[mapper_.template map<n>(e,i)][comp];
                }
          return value;
        }



        void interpolate (const C0GridFunction<G,RT,m>& u)
        {
      typedef typename GV::template Codim<0>::template Partition<All_Partition>::Iterator Iterator;
          std::vector<bool> visited(mapper_.size());
          for (int i=0; i<mapper_.size(); i++) visited[i] = false;

          Iterator eendit = gv.template end<0,All_Partition>();
          for (Iterator it = gv.template begin<0,All_Partition>(); it!=eendit; ++it)
                {
                  Dune::GeometryType gt = it->type();
                  for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1).size(); ++i)
                        if (!visited[mapper_.template map<n>(*it,i)])
                          {
                                for (int c=0; c<m; c++)
                                  (*coeff)[mapper_.template map<n>(*it,i)][c] = 
                                        u.evallocal(c,*it,Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].position());
                                visited[mapper_.template map<n>(*it,i)] = true;
                          }
                }
        }

        void interpolate (const P0Function<GV,RT,m>& u)
        {
          typedef typename GV::template Codim<0>::Iterator Iterator;
          std::vector<char> counter(mapper_.size());
          for (int i=0; i<mapper_.size(); i++) counter[i] = 0;

          for (int i=0; i<(*coeff).size(); i++)
                (*coeff)[i] = 0;
          Iterator eendit = gv.template end<0>();
          for (Iterator it = gv.template begin<0>(); it!=eendit; ++it)
                {
                  Dune::GeometryType gt = it->type();
                  for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1).size(); ++i)
                        {
                          for (int c=0; c<m; c++)
                                (*coeff)[mapper_.template map<n>(*it,i)][c] += 
                                  u.evallocal(c,*it,Dune::LagrangeShapeFunctions<DT,RT,n>::general(gt,1)[i].position());
                          counter[mapper_.template map<n>(*it,i)] += 1;
                        }
                }
          for (int i=0; i<counter.size(); i++)
                for (int c=0; c<m; c++)
                  (*coeff)[i][c] /= counter[i];
        }


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


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


        void fillIndexInfoFromGrid (P1IndexInfoFromGrid& info)
        {
          P1ExtendOverlap<G,GV,VM,LC> extender(lc);
          extender.fillIndexInfoFromGrid(gv.grid(),gv,mapper_,info);
        }


        void preAdapt ()
        {
        }

        void postAdapt (P1FunctionManager<G,RT>& manager)
        {
          typedef typename G::template Codim<n>::LeafIterator VLeafIterator;
          typedef typename G::template Codim<0>::LeafIterator ELeafIterator;
          typedef typename G::template Codim<n>::LevelIterator VLevelIterator;
          typedef typename G::template Codim<0>::LevelIterator ELevelIterator;
          typedef typename G::template Codim<0>::EntityPointer EEntityPointer;

          // \todo check that function is only called for data with respect to leafs
          // save the current representation
          oldcoeff = coeff;

          // allow mapper to recompute its internal sizes
          mapper_.update();

          // overlap extension, recompute extra DOFs
          if (extendOverlap)
                {
                  // set of neighbors in global ids for border vertices
                  std::map<int,GIDSet> borderlinks;
                  std::map<IdType,int> gid2index;
                  extraDOFs = 0;

                  // compute extension
                  P1ExtendOverlap<G,GV,VM,LC> extender(lc);
                  extender.extend(gv.grid(),gv,mapper_,borderlinks,extraDOFs,gid2index);
                }

          // allocate data with new size (while keeping the old data ...)
          try {
                coeff = new RepresentationType(mapper_.size()+extraDOFs); // allocate new representation
          }
          catch (std::bad_alloc) {
                std::cerr << "not enough memory in P1Function update" << std::endl;
                throw; // rethrow exception
          }
          std::cout << "P1  function enlarged to " << mapper_.size() << " components" << std::endl;

          // vector of flags to store which vertex has been handled already
          std::vector<bool> visited(mapper_.size());
          for (int i=0; i<mapper_.size(); i++) visited[i] = false;

          // now loop over the NEW mesh to copy the data that was already in the OLD mesh
          VLeafIterator veendit = gv.grid().template leafend<n>();
          for (VLeafIterator it = gv.grid().template leafbegin<n>(); it!=veendit; ++it)
                {
                  // lookup in mapper
                  int i;
                  if (manager.savedMap().contains(*it,i))
                        {
//                        std::cout << " found vertex=" << it->geometry()[0] 
//                                              << " at i=" << i
//                                              << " oldindex=" << manager.oldIndex()[i]
//                                              << " newindex=" << mapper_.map(*it)
//                                              << std::endl;
                          // the vertex existed already in the old mesh, copy data
                          for (int c=0; c<m; c++)
                                (*coeff)[mapper_.map(*it)][c] = (*oldcoeff)[manager.oldIndex()[i]][c];
                          // ... and mark as visited
                          visited[mapper_.map(*it)] = true;
                        }
                }

          // now loop the second time to interpolate the new coefficients
          // new implementation using interpolation on codim 0
          for (int level=1; level<=gv.grid().maxLevel(); ++level)
                {
                  ELevelIterator elendit = gv.grid().template lend<0>(level);
                  for (ELevelIterator it = gv.grid().template lbegin<0>(level); it!=elendit; ++it)
                        {
                          GeometryType gte = it->type();
                          for (int i=0; i<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gte,1).size(); ++i)
                                {
                                  int index = mapper_.template map<n>(*it,i);
                                  if (!visited[index])
                                        {
                                          // OK, this is a new vertex
                                          EEntityPointer father=it->father(); // the father element
                                          GeometryType gtf = father->type(); // fathers type
                                          const FieldVector<DT,n>& cpos=Dune::LagrangeShapeFunctions<DT,RT,n>::general(gte,1)[i].position();
                                          FieldVector<DT,n> pos = it->geometryInFather().global(cpos); // map corner to father element
                                          for (int c=0; c<m; c++)
                                                (*coeff)[index][c] = 0;
                                          for (int j=0; j<Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1).size(); ++j)
                                                {
                                                  RT basefuncvalue = Dune::LagrangeShapeFunctions<DT,RT,n>::general(gtf,1)[j].evaluateFunction(0,pos);
                                                  for (int c=0; c<m; c++)
                                                        (*coeff)[index][c] += basefuncvalue * (*coeff)[mapper_.template map<n>(*father,j)][c];
                                                  //                                      std::cout << "  corner=" << i 
                                                  //                                                            << " cpos=" << father->geometry()[i]
                                                  //                                                            << " u=" << (*coeff)[mapper_.template map<n>(*father,i)]
                                                  //                                                            << std::endl;
                                                }
                                          //                              std::cout << "index=" << mapper_.map(*it) << " value=" << value << std::endl;
                                          visited[index] = true;                                  
                                        }
                                }
                        }
                }

          // now really delete old representation
          if (oldcoeff!=0) delete oldcoeff;
          oldcoeff = 0;
        }

        const VM& mapper () const
        {
          return mapper_;
        }
 
        const GV& gridview () const
        {
          return gv;
        }
 
        void vtkout (VTKWriter<GV>& vtkwriter, std::string s) const
        {
          typename VTKWriter<GV>::VTKFunction *p = new VTKGridFunctionWrapper<GV,RT,m>(*this,s);
          vtkwriter.addVertexData(p);
        }

  private:
        // The GridView.  Don't use a reference here since Grid::levelView
        // returns a temporary
        const GV gv;

        // we need a mapper
        VM mapper_;

        // level or leafwise communication object
        LC lc;

        // extra DOFs from extending nonoverlapping to overlapping grid
        int extraDOFs;
        bool extendOverlap;

        // and a dynamically allocated vector
        RepresentationType* coeff;

        // saved pointer in update phase
        RepresentationType* oldcoeff;
  };


  template<class G, class RT, int m=1>
  class LeafP1Function : public P1Function<typename G::LeafGridView,RT,LeafCommunicate<G>,m>
  {
  public:
        LeafP1Function (const G& grid, bool extendoverlap=false) 
          : GridFunctionGlobalEvalDefault<typename G::LeafGridView,RT,m>(grid.leafView())
          , P1Function<typename G::LeafGridView,RT,LeafCommunicate<G>,m>(grid.leafView(),LeafCommunicate<G>(grid),extendoverlap)
        {}
  };


  template<class G, class RT, int m=1>
  class LevelP1Function : public P1Function<typename G::LevelGridView,RT,LevelCommunicate<G>,m>
  {
  public:
        LevelP1Function (const G& grid, int level, bool extendoverlap=false) 
          : GridFunctionGlobalEvalDefault<typename G::LevelGridView,RT,m>(grid.levelView(level))
          , P1Function<typename G::LevelGridView,RT,LevelCommunicate<G>,m>(grid.levelView(level),LevelCommunicate<G>(grid,level),extendoverlap)
        {}
  };


  template<class G, class RT>
  class P1FunctionManager {
        enum {dim=G::dimension};
        typedef typename G::ctype DT;
        typedef typename G::template Codim<dim>::LeafIterator VLeafIterator;
        typedef typename G::template Codim<0>::EntityPointer EEntityPointer;

        template<int dim>
        struct P1Layout
        {
          bool contains (Dune::GeometryType gt)
          {
              return gt.dim() == 0;
          }
        }; 

  public:

        P1FunctionManager (const G& g) : mapper(g,g.leafIndexSet()), grid(g), savedmap(g)
        {       
          // allocate index array to correct size (this possible for vertex data)
          oldindex.resize(mapper.size());

          // and allocate the universal mapper to acces the old indices
          savedmap.clear(); //should be empty already

          // now loop over all vertices and copy the index provided by the mapper
          VLeafIterator veendit = grid.template leafend<dim>();
          for (VLeafIterator it = grid.template leafbegin<dim>(); it!=veendit; ++it)
                {
                  oldindex[savedmap.map(*it)] = mapper.map(*it);
                }
        }

        const GlobalUniversalMapper<G>& savedMap ()
        {
          return savedmap;
        }

        const std::vector<int>& oldIndex ()
        {
          return oldindex;
        }

  private:
        // we need a mapper
        MultipleCodimMultipleGeomTypeMapper<G,typename G::template Codim<0>::LeafIndexSet,P1Layout> mapper;

        // store a reference to the grid that is managed
        const G& grid;

        // We need a persistent consecutive enumeration 
        GlobalUniversalMapper<G> savedmap;

        // The old leaf indices are stored in a dynamically allocated vector
        std::vector<int> oldindex;
  };


}
#endif

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