dune-grid  2.2.0
yaspgrid.hh
Go to the documentation of this file.
00001 #ifndef DUNE_YASPGRID_HH
00002 #define DUNE_YASPGRID_HH
00003 
00004 #include<iostream>
00005 #include<vector>
00006 #include<algorithm>
00007 #include<stack>
00008 
00009 // either include stdint.h or provide fallback for uint8_t
00010 #if HAVE_STDINT_H
00011 #include<stdint.h>
00012 #else
00013 typedef unsigned char uint8_t;
00014 #endif
00015 
00016 #include <dune/grid/common/grid.hh>     // the grid base classes
00017 #include <dune/grid/yaspgrid/grids.hh>  // the yaspgrid base classes
00018 #include <dune/grid/common/capabilities.hh> // the capabilities
00019 #include <dune/common/misc.hh>
00020 #include <dune/common/bigunsignedint.hh>
00021 #include <dune/common/typetraits.hh>
00022 #include <dune/common/collectivecommunication.hh>
00023 #include <dune/common/mpihelper.hh>
00024 #include <dune/geometry/genericgeometry/topologytypes.hh>
00025 #include <dune/grid/common/indexidset.hh>
00026 #include <dune/grid/common/datahandleif.hh>
00027 
00028 
00029 #if HAVE_MPI
00030 #include <dune/common/mpicollectivecommunication.hh>
00031 #endif
00032 
00040 namespace Dune {
00041 
00042 //************************************************************************
00046 typedef double yaspgrid_ctype;
00047 static const yaspgrid_ctype yasptolerance=1E-13; // tolerance in coordinate computations
00048 
00049 /* some sizes for building global ids
00050  */
00051 const int yaspgrid_dim_bits = 24;   // bits for encoding each dimension
00052 const int yaspgrid_level_bits = 6;  // bits for encoding level number
00053 const int yaspgrid_codim_bits = 4;  // bits for encoding codimension
00054 
00055 
00056 //************************************************************************
00057 // forward declaration of templates
00058 
00059 template<int dim>                             class YaspGrid;
00060 template<int mydim, int cdim, class GridImp>  class YaspGeometry;
00061 template<int codim, int dim, class GridImp>   class YaspEntity;
00062 template<int codim, class GridImp>            class YaspEntityPointer;
00063 template<int codim, class GridImp>            class YaspEntitySeed;
00064 template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
00065 template<class GridImp>            class YaspIntersectionIterator;
00066 template<class GridImp>            class YaspIntersection;
00067 template<class GridImp>            class YaspHierarchicIterator;
00068 template<class GridImp>            class YaspLevelIndexSet;
00069 template<class GridImp>            class YaspLeafIndexSet;
00070 template<class GridImp>            class YaspGlobalIdSet;
00071 
00072 namespace FacadeOptions
00073 {
00074   
00075   template<int dim, int mydim, int cdim>
00076   struct StoreGeometryReference<mydim, cdim, YaspGrid<dim>, YaspGeometry>
00077   {
00078     static const bool v = false;
00079   };
00080   
00081   template<int dim, int mydim, int cdim>
00082   struct StoreGeometryReference<mydim, cdim, const YaspGrid<dim>, YaspGeometry>
00083   {
00084     static const bool v = false;
00085   };
00086   
00087 }
00088 
00089 //========================================================================
00090 // The transformation describing the refinement rule
00091 
00092 template<int dim, class GridImp>
00093 struct YaspFatherRelativeLocalElement {
00094   static const array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > _geo;
00095   static array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > initSons()
00096   {
00097     array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > geo;
00098     FieldVector<yaspgrid_ctype,dim> midpoint(0.25);
00099     FieldVector<yaspgrid_ctype,dim> extension(0.5);
00100     for (int i=0; i<(1<<dim); i++)
00101     {
00102       midpoint = 0.25;
00103       for (int k=0; k<dim; k++)
00104       {
00105         if (i&(1<<k))
00106           midpoint[k] = 0.75;
00107       }
00108       geo[i] = YaspGeometry<dim,dim,GridImp>(midpoint, extension);
00109     }
00110     return geo;
00111   }
00112 };
00113 
00114 template<int dim, class GridImp>
00115 const array<YaspGeometry<dim,dim,GridImp>, (1<<dim)>
00116 YaspFatherRelativeLocalElement<dim, GridImp>::_geo =
00117     YaspFatherRelativeLocalElement<dim, GridImp>::initSons();
00118 
00119 //========================================================================
00127 //========================================================================
00128 
00130 template<int mydim,int cdim, class GridImp>
00131 class YaspGeometry : public GeometryDefaultImplementation<mydim,cdim,GridImp,YaspGeometry>
00132 {
00133 public:
00135   typedef typename GridImp::ctype ctype;
00136   
00138   GeometryType type () const
00139   {
00140     return GeometryType(GeometryType::cube,mydim);
00141   }
00142 
00144   bool affine() const { return true; }
00145 
00147   int corners () const
00148   {
00149     return 1<<mydim;
00150   }
00151 
00153   FieldVector< ctype, cdim > corner ( const int i ) const
00154   {
00155     assert( i >= 0 && i < (int) coord_.N() ); 
00156     FieldVector<ctype, cdim>& c = coord_[i];
00157     int bit=0;
00158     for (int k=0; k<cdim; k++) // run over all directions in world
00159     {
00160       if (k==missing)
00161       {
00162         c[k] = midpoint[k];
00163         continue;
00164       }
00165       //k is not the missing direction
00166       if (i&(1<<bit))  // check whether bit is set or not
00167         c[k] = midpoint[k]+0.5*extension[k]; // bit is 1 in i
00168       else
00169         c[k] = midpoint[k]-0.5*extension[k]; // bit is 0 in i
00170       bit++; // we have processed a direction
00171     }
00172       
00173     return c;
00174   }
00175 
00177   FieldVector< ctype, cdim > center ( ) const
00178   {
00179     return midpoint;
00180   }
00181 
00183   FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const
00184   {
00185     FieldVector<ctype, cdim> g;
00186     int bit=0;
00187     for (int k=0; k<cdim; k++)
00188       if (k==missing)
00189         g[k] = midpoint[k];
00190       else
00191       {
00192         g[k] = midpoint[k] + (local[bit]-0.5)*extension[k];
00193         bit++;
00194       }
00195     return g;
00196   }
00197 
00199   FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const
00200   {
00201     FieldVector<ctype, mydim> l; // result
00202     int bit=0;
00203     for (int k=0; k<cdim; k++)
00204       if (k!=missing)
00205       {
00206         l[bit] = (global[k]-midpoint[k])/extension[k] + 0.5;
00207         bit++;
00208       }
00209     return l;
00210   }
00211 
00213   ctype volume () const
00214   {
00215     ctype volume=1.0;
00216     for (int k=0; k<cdim; k++)
00217       if (k!=missing) volume *= extension[k];
00218     return volume;
00219   }
00220 
00223   ctype integrationElement (const FieldVector<ctype, mydim>& local) const
00224   {
00225     return volume();
00226   }
00227 
00229   FieldMatrix<ctype,mydim,cdim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
00230   {
00231     JT = 0.0;
00232     int k=0;
00233     for (int i=0; i<cdim; ++i)
00234     {
00235       if (i != missing)
00236       {
00237         JT[k][i] = extension[i]; // set diagonal element
00238         k++;
00239       }
00240     }
00241     return JT;
00242   }
00244   FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
00245   {
00246     Jinv = 0.0;
00247     int k=0;
00248     for (int i=0; i<cdim; ++i)
00249     {
00250       if (i != missing)
00251       {
00252         Jinv[i][k] = 1.0/extension[i]; // set diagonal element
00253         k++;
00254       }
00255     }
00256     return Jinv;
00257   }
00258 
00260   YaspGeometry () {}
00261 
00263   YaspGeometry (const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, uint8_t& m)
00264     : midpoint(p), extension(h), missing(m)
00265   {
00266     if (cdim!=mydim+1)
00267       DUNE_THROW(GridError, "general YaspGeometry assumes cdim=mydim+1");
00268   }
00269 
00271   YaspGeometry (const YaspGeometry& other) 
00272     : midpoint(other.midpoint), 
00273       extension(other.extension), 
00274       missing(other.missing)
00275   {
00276   }
00277 
00279   void print (std::ostream& s) const
00280   {
00281     s << "YaspGeometry<"<<mydim<<","<<cdim<< "> ";
00282     s << "midpoint";
00283     for (int i=0; i<cdim; i++)
00284       s << " " << midpoint[i];
00285     s << " extension";
00286     for (int i=0; i<cdim; i++)
00287       s << " " << extension[i];
00288     s << " missing is " << missing;
00289   }
00290 
00291   // const YaspGeometry<mydim,cdim,GridImp>&
00292   // operator = (const YaspGeometry<mydim,cdim,GridImp>& g);
00293 
00294 private:
00295   // the element is fully defined by its midpoint the extension
00296   // in each direction and the missing direction.
00297   // Note cdim == mydim+1
00298 
00299   FieldVector<ctype, cdim> midpoint;    // the midpoint
00300   FieldVector<ctype, cdim> extension;   // the extension
00301   uint8_t missing;                      // the missing, i.e. constant direction
00302 
00303   // In addition we need memory in order to return references.
00304   // Possibly we should change this in the interface ...
00305   mutable FieldMatrix<ctype, mydim, cdim> JT;     // the transposed of the jacobian 
00306   mutable FieldMatrix<ctype, cdim, mydim> Jinv;   // the transposed of the jacobian inverse 
00307   mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, cdim> coord_;   // the coordinates 
00308 
00309 };
00310 
00311 
00312 
00314 template<int mydim, class GridImp>
00315 class YaspGeometry<mydim,mydim,GridImp> : public GeometryDefaultImplementation<mydim,mydim,GridImp,YaspGeometry>
00316 {
00317 public:
00318   typedef typename GridImp::ctype ctype;
00319   
00321   GeometryType type () const
00322   {
00323     return GeometryType(GeometryType::cube,mydim);
00324   }
00325 
00327   bool affine() const { return true; }
00328 
00330   int corners () const
00331   {
00332     return 1<<mydim;
00333   }
00334 
00336   const FieldVector<ctype, mydim>& operator[] (int i) const
00337   {
00338     return corner(i);
00339   }
00340 
00342   FieldVector< ctype, mydim > corner ( const int i ) const
00343   {
00344     assert( i >= 0 && i < (int) coord_.N() ); 
00345     FieldVector<ctype, mydim>& c = coord_[i];
00346     for (int k=0; k<mydim; k++)
00347       if (i&(1<<k))
00348         c[k] = midpoint[k]+0.5*extension[k]; // kth bit is 1 in i
00349       else
00350         c[k] = midpoint[k]-0.5*extension[k]; // kth bit is 0 in i
00351     return c;
00352   }
00353 
00355   FieldVector< ctype, mydim > center ( ) const
00356   {
00357     return midpoint;
00358   }
00359 
00361   FieldVector<ctype, mydim> global (const FieldVector<ctype, mydim>& local) const
00362   {
00363     FieldVector<ctype,mydim> g;
00364     for (int k=0; k<mydim; k++)
00365       g[k] = midpoint[k] + (local[k]-0.5)*extension[k];
00366     return g;
00367   }
00368 
00370   FieldVector<ctype, mydim> local (const FieldVector<ctype,mydim>& global) const
00371   {
00372     FieldVector<ctype, mydim> l; // result
00373     for (int k=0; k<mydim; k++)
00374       l[k] = (global[k]-midpoint[k])/extension[k] + 0.5;
00375     return l;
00376   }
00377 
00380   ctype integrationElement (const FieldVector<ctype, mydim>& local) const
00381   {
00382     return volume();
00383   }
00384 
00386   ctype volume () const
00387   {
00388     ctype vol=1.0;
00389     for (int k=0; k<mydim; k++) vol *= extension[k];
00390     return vol;
00391   }
00392 
00394   FieldMatrix<ctype,mydim,mydim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
00395   {
00396     for (int i=0; i<mydim; ++i)
00397     {
00398       JT[i] = 0.0;                // set column to zero
00399       JT[i][i] = extension[i]; // set diagonal element
00400     }
00401     return JT;
00402   }
00404   FieldMatrix<ctype,mydim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
00405   {
00406     for (int i=0; i<mydim; ++i)
00407     {
00408       Jinv[i] = 0.0;                // set column to zero
00409       Jinv[i][i] = 1.0/extension[i]; // set diagonal element
00410     }
00411     return Jinv;
00412   }
00413 
00415   YaspGeometry () {}
00416 
00418   YaspGeometry (const FieldVector<ctype, mydim>& p, const FieldVector<ctype, mydim>& h)
00419     : midpoint(p), extension(h)
00420   {}
00421 
00423   YaspGeometry (const YaspGeometry& other) 
00424     : midpoint(other.midpoint), 
00425       extension(other.extension) 
00426   {
00427   }
00428 
00430   void print (std::ostream& s) const
00431   {
00432     s << "YaspGeometry<"<<mydim<<","<<mydim<< "> ";
00433     s << "midpoint";
00434     for (int i=0; i<mydim; i++)
00435       s << " " << midpoint[i];
00436     s << " extension";
00437     for (int i=0; i<mydim; i++)
00438       s << " " << extension[i];
00439   }
00440 
00441   // const YaspGeometry<mydim,mydim,GridImp>&
00442   // operator = (const YaspGeometry<mydim,mydim,GridImp>& g);
00443 
00444 private:
00445   // the element is fully defined by midpoint and the extension
00446   // in each direction. References are used because this information
00447   // is known outside the element in many cases.
00448   // Note mydim==cdim
00449 
00450   FieldVector<ctype, mydim> midpoint;   // the midpoint
00451   FieldVector<ctype, mydim> extension;  // the extension
00452 
00453   // In addition we need memory in order to return references.
00454   // Possibly we should change this in the interface ...
00455   mutable FieldMatrix<ctype, mydim, mydim> Jinv,JT;   // the transpose of the jacobian and its inverse inverse
00456   mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, mydim> coord_;   // the coordinates 
00457 };
00458 
00460 template<int cdim, class GridImp>
00461 class YaspGeometry<0,cdim,GridImp> : public GeometryDefaultImplementation<0,cdim,GridImp,YaspGeometry>
00462 {
00463 public:
00464   typedef typename GridImp::ctype ctype;
00465   
00467   GeometryType type () const
00468   {
00469     return GeometryType(GeometryType::cube,0);
00470   }
00471 
00473   bool affine() const { return true; }
00474 
00476   int corners () const
00477   {
00478     return 1;
00479   }
00480 
00482   const FieldVector<ctype, cdim>& operator[] (int i) const
00483   {
00484     return position;
00485   }
00486 
00488   FieldVector< ctype, cdim > corner ( const int i ) const
00489   {
00490     return position;
00491   }
00492 
00494   FieldVector< ctype, cdim > center ( ) const
00495   {
00496     return position;
00497   }
00498 
00501   ctype integrationElement (const FieldVector<ctype, 0>& local) const
00502   {
00503     return 1.0;
00504   }
00505 
00507   FieldMatrix<ctype,0,cdim>& jacobianTransposed (const FieldVector<ctype, 0>& local) const
00508   {
00509     static FieldMatrix<ctype,0,cdim> JT(0.0);
00510     return JT;
00511   }
00513   FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
00514   {
00515     static FieldMatrix<ctype,cdim,0> Jinv(0.0);
00516     return Jinv;
00517   }
00518 
00520   YaspGeometry ()
00521   {}
00522   
00524   explicit YaspGeometry ( const FieldVector< ctype, cdim > &p )
00525     : position( p )
00526   {}
00527 
00528   YaspGeometry ( const FieldVector< ctype, cdim > &p, const FieldVector< ctype, cdim > &, uint8_t &)
00529     : position( p )
00530   {}
00531 
00533   void print (std::ostream& s) const
00534   {
00535     s << "YaspGeometry<"<<0<<","<<cdim<< "> ";
00536     s << "position " << position;
00537   }
00538 
00539   // const YaspGeometry<0,cdim,GridImp>&
00540   // operator = (const YaspGeometry<0,cdim,GridImp>& g);
00541 
00542 private:
00543   FieldVector<ctype, cdim> position; 
00544 };
00545 
00546 // operator<< for all YaspGeometrys
00547 template <int mydim, int cdim, class GridImp>
00548 inline
00549 std::ostream& operator<< (std::ostream& s, YaspGeometry<mydim,cdim,GridImp>& e)
00550 {
00551   e.print(s);
00552   return s;
00553 }
00554 
00555 //========================================================================
00563 //========================================================================
00564 
00565 template<int codim, int dim, class GridImp>
00566 class YaspSpecialEntity :
00567   public GridImp::template Codim<codim>::Entity
00568 {
00569 public:
00570   typedef typename GridImp::ctype ctype;
00571   
00572   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
00573   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
00574 
00575   YaspSpecialEntity(const GridImp* yg, const YGLI& g, const TSI& it) :
00576     GridImp::template Codim<codim>::Entity (YaspEntity<codim, dim, GridImp>(yg,g,it))
00577     {};
00578   YaspSpecialEntity(const YaspEntity<codim, dim, GridImp>& e) :
00579     GridImp::template Codim<codim>::Entity (e)
00580     {};
00581   const TSI& transformingsubiterator () const
00582   {
00583         return this->realEntity.transformingsubiterator();
00584   }
00585   const YGLI& gridlevel () const
00586   {
00587         return this->realEntity.gridlevel();
00588   }
00589   const GridImp * yaspgrid() const
00590   {
00591       return this->realEntity.yaspgrid();
00592   }
00593 };
00594 
00595 template<int codim, int dim, class GridImp>
00596 class YaspEntity 
00597 :  public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
00598 {
00599 public:
00600   typedef typename GridImp::ctype ctype;
00601   
00602   typedef typename GridImp::template Codim<codim>::Geometry Geometry;
00604   int level () const
00605   {
00606         DUNE_THROW(GridError, "YaspEntity not implemented");
00607   }
00608 
00610   int index () const
00611   {
00612         DUNE_THROW(GridError, "YaspEntity not implemented");
00613   }
00614 
00616   Geometry geometry () const
00617   {
00618         DUNE_THROW(GridError, "YaspEntity not implemented");
00619   }
00620 
00622   PartitionType partitionType () const
00623   {
00624         DUNE_THROW(GridError, "YaspEntity not implemented");
00625   }
00626 
00627   const GridImp * yaspgrid() const
00628   {
00629         DUNE_THROW(GridError, "YaspEntity not implemented");
00630   }
00631   
00632   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
00633   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
00634   YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
00635   {
00636         DUNE_THROW(GridError, "YaspEntity not implemented");
00637   }
00638 
00639   // IndexSets needs access to the private index methods
00640   friend class Dune::YaspLevelIndexSet<GridImp>;
00641   friend class Dune::YaspLeafIndexSet<GridImp>;
00642   friend class Dune::YaspGlobalIdSet<GridImp>;
00643   typedef typename GridImp::PersistentIndexType PersistentIndexType;
00644 
00646   PersistentIndexType persistentIndex () const 
00647   { 
00648     DUNE_THROW(GridError, "YaspEntity not implemented");
00649   }
00650 
00652   int compressedIndex () const 
00653   {
00654     DUNE_THROW(GridError, "YaspEntity not implemented");
00655   }
00656 
00658   int compressedLeafIndex () const 
00659   {
00660     DUNE_THROW(GridError, "YaspEntity not implemented");
00661   }
00662 };
00663 
00664 
00665 // specialization for codim=0
00666 template<int dim, class GridImp>
00667 class YaspEntity<0,dim,GridImp> 
00668 : public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
00669 {
00670   enum { dimworld = GridImp::dimensionworld };
00671 
00672   typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
00673 
00674 public:
00675   typedef typename GridImp::ctype ctype;
00676   
00677   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
00678   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
00679 
00680   typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
00681   typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
00682    
00683   template <int cd>
00684   struct Codim
00685   {
00686     typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
00687   };
00688 
00689   typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
00690   typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
00691   typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
00692   typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
00693   typedef typename GridImp::LeafIntersectionIterator  LeafIntersectionIterator;
00694   typedef typename GridImp::HierarchicIterator HierarchicIterator;
00695 
00697   typedef typename GridImp::PersistentIndexType PersistentIndexType;
00698 
00700   typedef typename YGrid<dim,ctype>::iTupel iTupel;
00701 
00702   // constructor
00703   YaspEntity (const GridImp * yg, const YGLI& g, const TSI& it)
00704     : _yg(yg), _it(it), _g(g)
00705   {
00706   }
00707 
00709   int level () const { return _g.level(); }
00710 
00712   int index () const { return _it.superindex(); } // superindex works also for iteration over subgrids
00713 
00715   int globalIndex () const {
00716     return _g.cell_global().index(_it.coord());
00717   }
00718 
00722   EntitySeed seed () const {
00723     return EntitySeed(_g.level(), _it.coord());
00724   }
00725 
00727   PartitionType partitionType () const
00728   {
00729         if (_g.cell_interior().inside(_it.coord())) return InteriorEntity;
00730         if (_g.cell_overlap().inside(_it.coord())) return OverlapEntity;
00731         DUNE_THROW(GridError, "Impossible GhostEntity " << _it.coord() << "\t"
00732                   << _g.cell_interior().origin() << "/" << _g.cell_interior().size());
00733         return GhostEntity;
00734   }
00735 
00737   Geometry geometry () const {
00738     // the element geometry
00739     GeometryImpl _geometry(_it.position(),_it.meshsize());
00740     return Geometry( _geometry );
00741   }
00742 
00745   template<int cc> int count () const
00746   {
00747     if (cc==dim) return 1<<dim;
00748     if (cc==1) return 2*dim;
00749     if (cc==dim-1) return dim*(1<<(dim-1));
00750     if (cc==0) return 1;
00751     DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
00752   }
00753 
00756   template<int cc>
00757   typename Codim<cc>::EntityPointer subEntity (int i) const
00758   {
00759       dune_static_assert( cc == dim || cc == 0 ,
00760           "YaspGrid only supports Entities with codim=dim and codim=0");
00761       // coordinates of the cell == coordinates of lower left corner
00762       if (cc==dim)
00763           {
00764                 iTupel coord = _it.coord();
00765 
00766                 // get corner from there
00767                 for (int k=0; k<dim; k++)
00768                   if (i&(1<<k)) (coord[k])++;
00769 
00770                 return YaspEntityPointer<cc,GridImp>(_yg,_g,_g.vertex_overlapfront().tsubbegin(coord));
00771           }
00772       if (cc==0)
00773           {
00774                 return YaspEntityPointer<cc,GridImp>(_yg,_g,_it);
00775           }
00776       DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
00777   }
00778 
00780   EntityPointer father () const
00781   {
00782         // check if coarse level exists
00783         if (_g.level()<=0)
00784           DUNE_THROW(GridError, "tried to call father on level 0");
00785 
00786         // yes, get iterator to it
00787         YGLI cg = _g.coarser();
00788 
00789         // coordinates of the cell
00790         iTupel coord = _it.coord();
00791 
00792         // get coordinates on next coarser level
00793         for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
00794 
00795         return YaspEntityPointer<0,GridImp>(_yg,cg,cg.cell_overlap().tsubbegin(coord));
00796   }
00797   
00799   bool hasFather () const
00800   {
00801       return (_g.level()>0);
00802   }
00803 
00813   LocalGeometry geometryInFather () const
00814   {
00815     // determine which son we are
00816     int son = 0;
00817     for (int k=0; k<dim; k++)
00818       if (_it.coord(k)%2)
00819         son += (1<<k);
00820     
00821     // configure one of the 2^dim transformations
00822     return LocalGeometry( YaspFatherRelativeLocalElement<dim,GridImp>::_geo[son] );
00823   }
00824 
00825   const TSI& transformingsubiterator () const
00826   {
00827         return _it;
00828   }
00829 
00830   const YGLI& gridlevel () const
00831   {
00832         return _g;
00833   }
00834 
00835   const GridImp* yaspgrid () const
00836   {
00837         return _yg;
00838   }
00839 
00840   bool isLeaf() const
00841   {
00842     return (_g.level() == _g.mg()->maxlevel());
00843   }
00844   
00847   bool isNew () const { return _yg->adaptRefCount > 0 && _g.mg()->maxlevel() < _g.level() + _yg->adaptRefCount; }
00848   
00851   bool mightVanish () const { return false; }
00852   // { return _yg->adaptRefCount < 0 && _g.mg()->maxlevel() < _g.level() - _yg->adaptRefCount; }
00853 
00855   IntersectionIterator ibegin () const
00856   {
00857     return YaspIntersectionIterator<GridImp>(*this,false);
00858   }
00859 
00861   LeafIntersectionIterator ileafbegin () const
00862   {
00863     // only if entity is leaf this iterator delivers intersections
00864     return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
00865   }
00866 
00868   LevelIntersectionIterator ilevelbegin () const
00869   {
00870     return ibegin(); 
00871   }
00872 
00874   IntersectionIterator iend () const
00875   {
00876     return YaspIntersectionIterator<GridImp>(*this,true);
00877   }
00878 
00880   LeafIntersectionIterator ileafend () const
00881   {
00882     return iend();
00883   }
00884 
00886   LevelIntersectionIterator ilevelend () const
00887   {
00888     return iend();
00889   }
00890 
00895   HierarchicIterator hbegin (int maxlevel) const
00896   {
00897     return YaspHierarchicIterator<GridImp>(_yg,_g,_it,maxlevel);
00898   }
00899 
00901   HierarchicIterator hend (int maxlevel) const
00902   {
00903     return YaspHierarchicIterator<GridImp>(_yg,_g,_it,_g.level());
00904   }
00905 
00906 private:
00907   // IndexSets needs access to the private index methods
00908   friend class Dune::YaspLevelIndexSet<GridImp>;
00909   friend class Dune::YaspLeafIndexSet<GridImp>;
00910   friend class Dune::YaspGlobalIdSet<GridImp>;
00911 
00913   PersistentIndexType persistentIndex () const 
00914   { 
00915     // get size of global grid
00916     const iTupel& size =  _g.cell_global().size();
00917 
00918     // get coordinate correction for periodic boundaries
00919     int coord[dim];
00920     for (int i=0; i<dim; i++)
00921       {
00922         coord[i] = _it.coord(i);
00923         if (coord[i]<0) coord[i] += size[i];
00924         if (coord[i]>=size[i]) coord[i] -= size[i];
00925       }
00926 
00927     // encode codim
00928     PersistentIndexType id(0);
00929 
00930     // encode level
00931     id = id << yaspgrid_level_bits;
00932     id = id+PersistentIndexType(_g.level());
00933     
00934 
00935     // encode coordinates
00936     for (int i=dim-1; i>=0; i--)
00937       {
00938         id = id << yaspgrid_dim_bits;
00939         id = id+PersistentIndexType(coord[i]);
00940       }
00941     
00942     return id;
00943   }
00944 
00946   int compressedIndex () const 
00947   {
00948     return _it.superindex();
00949   }
00950 
00952   int compressedLeafIndex () const 
00953   {
00954     return _it.superindex();
00955   }
00956 
00958   PersistentIndexType subPersistentIndex (int i, int cc) const
00959   {
00960     if (cc==0)
00961       return persistentIndex();
00962 
00963     // get position of cell, note that global origin is zero
00964     // adjust for periodic boundaries
00965     int coord[dim];
00966     for (int k=0; k<dim; k++)
00967       {
00968         coord[k] = _it.coord(k);
00969         if (coord[k]<0) coord[k] += _g.cell_global().size(k);
00970         if (coord[k]>=_g.cell_global().size(k)) coord[k] -= _g.cell_global().size(k);
00971       }
00972     
00973     if (cc==dim)
00974       {
00975         // transform to vertex coordinates
00976         for (int k=0; k<dim; k++)
00977           if (i&(1<<k)) (coord[k])++;
00978 
00979         // determine min number of trailing zeroes
00980         int trailing = 1000;
00981         for (int i=0; i<dim; i++)
00982           {
00983             // count trailing zeros
00984             int zeros = 0;
00985             for (int j=0; j<_g.level(); j++)
00986               if (coord[i]&(1<<j))
00987                 break;
00988               else
00989                 zeros++;
00990             trailing = std::min(trailing,zeros);
00991           } 
00992 
00993         // determine the level of this vertex
00994         int level = _g.level()-trailing;
00995 
00996         // encode codim
00997         PersistentIndexType id(dim);
00998         
00999         // encode level
01000         id = id << yaspgrid_level_bits;
01001         id = id+PersistentIndexType(level);
01002         
01003         // encode coordinates
01004         for (int i=dim-1; i>=0; i--)
01005           {
01006             id = id << yaspgrid_dim_bits;
01007             id = id+PersistentIndexType(coord[i]>>trailing);
01008           }
01009         
01010         return id;
01011       }
01012 
01013     if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
01014       {
01015         // Idea: Use the doubled grid to assign coordinates to faces
01016 
01017         // ivar is the direction that varies
01018         int ivar=i/2; 
01019 
01020         // compute position from cell position
01021         for (int k=0; k<dim; k++)
01022           coord[k] = coord[k]*2 + 1; // the doubled grid
01023         if (i%2) 
01024           coord[ivar] += 1;
01025         else
01026           coord[ivar] -= 1;
01027 
01028         // encode codim
01029         PersistentIndexType id(1);
01030         
01031         // encode level
01032         id = id << yaspgrid_level_bits;
01033         id = id+PersistentIndexType(_g.level());
01034         
01035         // encode coordinates
01036         for (int i=dim-1; i>=0; i--)
01037           {
01038             id = id << yaspgrid_dim_bits;
01039             id = id+PersistentIndexType(coord[i]);
01040           }
01041         
01042         return id;
01043       }
01044 
01045     // map to old numbering
01046     static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
01047     i = edge[i];
01048 
01049     if (cc==dim-1) // edges, exist only for dim>2
01050       {
01051         // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
01052 
01053         // number of entities per direction
01054         int m=1<<(dim-1);
01055 
01056         // ifix is the direction that is fixed
01057         int ifix=(dim-1)-(i/m); 
01058 
01059         // compute position from cell position
01060         int bit=1;
01061         for (int k=0; k<dim; k++)
01062           {
01063             coord[k] = coord[k]*2+1; // cell position in doubled grid
01064             if (k==ifix) continue;
01065             if ((i%m)&bit) coord[k] += 1; else coord[k] -= 1;
01066             bit *= 2;
01067           } 
01068 
01069         // encode codim
01070         PersistentIndexType id(dim-1);
01071         
01072         // encode level
01073         id = id << yaspgrid_level_bits;
01074         id = id+PersistentIndexType(_g.level());
01075         
01076         // encode coordinates
01077         for (int i=dim-1; i>=0; i--)
01078           {
01079             id = id << yaspgrid_dim_bits;
01080             id = id+PersistentIndexType(coord[i]);
01081           }
01082         
01083         return id;
01084       }
01085 
01086     DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
01087   }
01088 
01090   int subCompressedIndex (int i, int cc) const
01091   {
01092     if (cc==0)
01093       return compressedIndex();
01094     
01095     // get cell position relative to origin of local cell grid
01096     iTupel coord;
01097     for (int k=0; k<dim; ++k) 
01098       coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
01099 
01100     if (cc==dim) // vertices
01101       {
01102         // transform cell coordinate to corner coordinate
01103         for (int k=0; k<dim; k++)
01104           if (i&(1<<k)) (coord[k])++;
01105 
01106         // do lexicographic numbering
01107         int index = coord[dim-1];
01108         for (int k=dim-2; k>=0; --k)
01109           index = (index*(_g.cell_overlap().size(k)+1))+coord[k];
01110         return index;
01111       }
01112 
01113     if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
01114       {
01115         // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
01116 
01117         // ivar is the direction that varies
01118         int ivar=i/2; 
01119 
01120         // compute position from cell position
01121         if (i%2) coord[ivar] += 1; 
01122 
01123         // do lexicographic numbering
01124         int index = coord[dim-1];
01125         for (int k=dim-2; k>=0; --k)
01126           if (k==ivar)
01127             index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
01128           else
01129             index = (index*(_g.cell_overlap().size(k)))+coord[k];
01130 
01131         // add size of all subsets for smaller directions
01132         for (int j=0; j<ivar; j++)
01133           {
01134             int n=_g.cell_overlap().size(j)+1;
01135             for (int l=0; l<dim; l++)
01136               if (l!=j) n *= _g.cell_overlap().size(l);
01137             index += n;
01138           }
01139 
01140         return index;
01141       }
01142 
01143     // map to old numbering
01144     static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
01145     i = edge[i];
01146 
01147     if (cc==dim-1) // edges, exist only for dim>2
01148       {
01149         // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
01150 
01151         // number of entities per direction
01152         int m=1<<(dim-1);
01153 
01154         // ifix is the direction that is fixed
01155         int ifix=(dim-1)-(i/m); 
01156 
01157         // compute position from cell position
01158         int bit=1;
01159         for (int k=0; k<dim; k++)
01160           {
01161             if (k==ifix) continue;
01162             if ((i%m)&bit) coord[k] += 1;
01163             bit *= 2;
01164           } 
01165 
01166         // do lexicographic numbering
01167         int index = coord[dim-1];
01168         for (int k=dim-2; k>=0; --k)
01169           if (k!=ifix)
01170             index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
01171           else
01172             index = (index*(_g.cell_overlap().size(k)))+coord[k];
01173 
01174         // add size of all subsets for smaller directions
01175         for (int j=dim-1; j>ifix; j--)
01176           {
01177             int n=_g.cell_overlap().size(j);
01178             for (int l=0; l<dim; l++)
01179               if (l!=j) n *= _g.cell_overlap().size(l)+1;
01180             index += n;
01181           }
01182 
01183         return index;
01184       }
01185 
01186     DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
01187   }
01188 
01190   int subCompressedLeafIndex (int i, int cc) const
01191   {
01192     if (cc==0)
01193       return compressedIndex();
01194     
01195     // get cell position relative to origin of local cell grid
01196     iTupel coord;
01197     for (int k=0; k<dim; ++k) 
01198       coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
01199 
01200     if (cc==dim) // vertices
01201       {
01202         // transform cell coordinate to corner coordinate
01203         for (int k=0; k<dim; k++)
01204           if (i&(1<<k)) (coord[k])++;
01205 
01206         // move coordinates up to maxlevel
01207         for (int k=0; k<dim; k++)
01208           coord[k] = coord[k]<<(_g.mg()->maxlevel()-_g.level());
01209 
01210         // do lexicographic numbering
01211         int index = coord[dim-1];
01212         for (int k=dim-2; k>=0; --k)
01213           index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
01214         return index;
01215       }
01216 
01217     if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
01218       {
01219         // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
01220 
01221         // ivar is the direction that varies
01222         int ivar=i/2; 
01223 
01224         // compute position from cell position
01225         if (i%2) coord[ivar] += 1; 
01226 
01227         // do lexicographic numbering
01228         int index = coord[dim-1];
01229         for (int k=dim-2; k>=0; --k)
01230           if (k==ivar)
01231             index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
01232           else
01233             index = (index*(_g.cell_overlap().size(k)))+coord[k];
01234 
01235         // add size of all subsets for smaller directions
01236         for (int j=0; j<ivar; j++)
01237           {
01238             int n=_g.cell_overlap().size(j)+1;
01239             for (int l=0; l<dim; l++)
01240               if (l!=j) n *= _g.cell_overlap().size(l);
01241             index += n;
01242           }
01243 
01244         return index;
01245       }
01246 
01247     // map to old numbering
01248     static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
01249     i = edge[i];
01250 
01251     if (cc==dim-1) // edges, exist only for dim>2
01252       {
01253         // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
01254 
01255         // number of entities per direction
01256         int m=1<<(dim-1);
01257 
01258         // ifix is the direction that is fixed
01259         int ifix=(dim-1)-(i/m); 
01260 
01261         // compute position from cell position
01262         int bit=1;
01263         for (int k=0; k<dim; k++)
01264           {
01265             if (k==ifix) continue;
01266             if ((i%m)&bit) coord[k] += 1;
01267             bit *= 2;
01268           } 
01269 
01270         // do lexicographic numbering
01271         int index = coord[dim-1];
01272         for (int k=dim-2; k>=0; --k)
01273           if (k!=ifix)
01274             index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
01275           else
01276             index = (index*(_g.cell_overlap().size(k)))+coord[k];
01277 
01278         // add size of all subsets for smaller directions
01279         for (int j=dim-1; j>ifix; j--)
01280           {
01281             int n=_g.cell_overlap().size(j);
01282             for (int l=0; l<dim; l++)
01283               if (l!=j) n *= _g.cell_overlap().size(l)+1;
01284             index += n;
01285           }
01286 
01287         return index;
01288       }
01289 
01290     DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
01291   }
01292 
01293   const GridImp * _yg;      // access to YaspGrid
01294   const TSI& _it;           // position in the grid level
01295   const YGLI& _g;           // access to grid level
01296 };
01297 
01298 
01299 // specialization for codim=dim (vertex)
01300 template<int dim, class GridImp>
01301 class YaspEntity<dim,dim,GridImp> 
01302 : public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
01303 {
01304   enum { dimworld = GridImp::dimensionworld };
01305 
01306   typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
01307 
01308 public:
01309   typedef typename GridImp::ctype ctype;
01310   
01311   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
01312   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
01313 
01314   typedef typename GridImp::template Codim<dim>::Geometry Geometry;
01315 
01316   template <int cd>
01317   struct Codim
01318   {
01319     typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
01320   };
01321 
01322   typedef typename GridImp::template Codim<dim>::EntityPointer EntityPointer;
01323   typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
01324 
01326   typedef typename GridImp::PersistentIndexType PersistentIndexType;
01327 
01329   typedef typename YGrid<dim,ctype>::iTupel iTupel;
01330 
01331   // constructor
01332   YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
01333     : _yg(yg), _it(it), _g(g)
01334   {}
01335 
01337   int level () const {return _g.level();}
01338 
01340   int index () const {return _it.superindex();}
01341 
01343   int globalIndex () const { return _g.cell_global().index(_it.coord()); }
01344 
01348   EntitySeed seed () const {
01349     return EntitySeed(_g.level(), _it.coord());
01350   }
01351 
01353   Geometry geometry () const {
01354     GeometryImpl _geometry(_it.position());
01355     return Geometry( _geometry );
01356   }
01357 
01359   PartitionType partitionType () const
01360   {
01361     if (_g.vertex_interior().inside(_it.coord())) return InteriorEntity;
01362     if (_g.vertex_interiorborder().inside(_it.coord())) return BorderEntity;
01363     if (_g.vertex_overlap().inside(_it.coord())) return OverlapEntity;
01364     if (_g.vertex_overlapfront().inside(_it.coord())) return FrontEntity;
01365     return GhostEntity;
01366   }
01367 
01368 private:
01369   // IndexSets needs access to the private index methods
01370   friend class Dune::YaspLevelIndexSet<GridImp>;
01371   friend class Dune::YaspLeafIndexSet<GridImp>;
01372   friend class Dune::YaspGlobalIdSet<GridImp>;
01373 
01375   PersistentIndexType persistentIndex () const 
01376   { 
01377     // get coordinate and size of global grid
01378     const iTupel& size =  _g.vertex_global().size();
01379     int coord[dim];
01380 
01381     // correction for periodic boundaries
01382     for (int i=0; i<dim; i++)
01383       {
01384         coord[i] = _it.coord(i);
01385         if (coord[i]<0) coord[i] += size[i];
01386         if (coord[i]>=size[i]) coord[i] -= size[i];
01387       }
01388 
01389     // determine min number of trailing zeroes
01390     int trailing = 1000;
01391     for (int i=0; i<dim; i++)
01392       {
01393         // count trailing zeros
01394         int zeros = 0;
01395         for (int j=0; j<_g.level(); j++)
01396           if (coord[i]&(1<<j))
01397             break;
01398           else
01399             zeros++;
01400         trailing = std::min(trailing,zeros);
01401       } 
01402 
01403     // determine the level of this vertex
01404     int level = _g.level()-trailing;
01405 
01406     // encode codim
01407     PersistentIndexType id(dim);
01408 
01409     // encode level
01410     id = id << yaspgrid_level_bits;
01411     id = id+PersistentIndexType(level);
01412     
01413     // encode coordinates
01414     for (int i=dim-1; i>=0; i--)
01415       {
01416         id = id << yaspgrid_dim_bits;
01417         id = id+PersistentIndexType(coord[i]>>trailing);
01418       }
01419     
01420     return id;
01421   }
01422 
01424   int compressedIndex () const { return _it.superindex();}
01425 
01427   int compressedLeafIndex () const 
01428   { 
01429     if (_g.level()==_g.mg()->maxlevel())
01430       return _it.superindex();
01431 
01432     // not on leaf level, interpolate to finest grid
01433     int coord[dim];
01434     for (int i=0; i<dim; i++) coord[i] = _it.coord(i)-(_g).vertex_overlap().origin(i);
01435 
01436     // move coordinates up to maxlevel (multiply by 2 for each level
01437     for (int k=0; k<dim; k++)
01438       coord[k] = coord[k]*(1<<(_g.mg()->maxlevel()-_g.level()));
01439 
01440     // do lexicographic numbering
01441     int index = coord[dim-1];
01442     for (int k=dim-2; k>=0; --k)
01443       index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
01444     return index;
01445   }
01446 
01447 public:  
01448   const TSI& transformingsubiterator() const { return _it; }
01449   const YGLI& gridlevel() const { return _g; }
01450   const GridImp * yaspgrid() const { return _yg; }
01451 protected:
01452   const GridImp * _yg;            // access to YaspGrid
01453   const TSI& _it;                 // position in the grid level
01454   const YGLI& _g;                 // access to grid level
01455   // temporary object
01456   mutable FieldVector<ctype, dim> loc;   // always computed before being returned
01457 };
01458 
01459 //========================================================================
01464 //========================================================================
01465 
01466 template<class GridImp>
01467 class YaspIntersection
01468 {
01469     enum { dim=GridImp::dimension };
01470     enum { dimworld=GridImp::dimensionworld };  
01471     typedef typename GridImp::ctype ctype;
01472     YaspIntersection();
01473     YaspIntersection& operator = (const YaspIntersection&);
01474 
01475     typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
01476     typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
01477 
01478 public:
01479     // types used from grids
01480     typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
01481     typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
01482     typedef typename GridImp::template Codim<0>::Entity Entity;
01483     typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
01484     typedef typename GridImp::template Codim<1>::Geometry Geometry;
01485     typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
01486     typedef YaspSpecialEntity<0,dim,GridImp> SpecialEntity;
01487     typedef Dune::Intersection<const GridImp, Dune::YaspIntersectionIterator> Intersection;
01488   
01489     // void update() const {
01490     //     const_cast<YaspIntersection*>(this)->update();
01491     // }
01492     void update() const {
01493         if (_count == 2*_dir + _face || _count >= 2*dim)
01494             return;
01495 
01496         // cleanup old stuff
01497         _outside.transformingsubiterator().move(_dir,1-2*_face); // move home
01498         _pos_world[_dir] = _inside.transformingsubiterator().position(_dir);
01499                 
01500         // update face info
01501         _dir = _count / 2;
01502         _face = _count % 2;
01503                 
01504         // move transforming iterator
01505         _outside.transformingsubiterator().move(_dir,-1+2*_face);
01506                 
01507         // make up faces
01508         _pos_world[_dir] += (-0.5+_face)*_inside.transformingsubiterator().meshsize(_dir);
01509     }
01510 
01512     void increment() {
01513         _count += (_count < 2*dim);
01514     }
01515 
01517     bool equals (const YaspIntersection& other) const
01518     {
01519         return _inside.equals(other._inside) && _count == other._count;
01520     }
01521 
01526     bool boundary () const
01527     {
01528 #if 1
01529         return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 < _inside.gridlevel().cell_global().min(_count/2)
01530             ||
01531             _inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 > _inside.gridlevel().cell_global().max(_count/2));
01532 #else
01533         update();
01534         // The transforming iterator can be safely moved beyond the boundary.
01535         // So we only have to compare against the cell_global grid
01536         return (_outside.transformingsubiterator().coord(_dir) < _inside.gridlevel().cell_global().min(_dir)
01537             ||
01538             _outside.transformingsubiterator().coord(_dir) > _inside.gridlevel().cell_global().max(_dir));
01539 #endif
01540     }
01541 
01543     bool neighbor () const
01544     {
01545 #if 1
01546         return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 >= _inside.gridlevel().cell_overlap().min(_count/2)
01547             &&
01548             _inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 <= _inside.gridlevel().cell_overlap().max(_count/2));
01549 #else
01550         update();
01551         return (_outside.transformingsubiterator().coord(_dir) >= _inside.gridlevel().cell_overlap().min(_dir)
01552             &&
01553             _outside.transformingsubiterator().coord(_dir) <= _inside.gridlevel().cell_overlap().max(_dir));
01554 #endif
01555     }
01556 
01558     bool conforming () const
01559     {
01560         return true;
01561     }
01562 
01565     EntityPointer inside() const
01566     {
01567         return _inside;
01568     }
01569   
01572     EntityPointer outside() const
01573     {
01574         update();
01575         return _outside;
01576     }
01577 
01580     int boundaryId() const
01581     {
01582         if(boundary()) return indexInInside()+1;
01583         return 0;
01584     }
01585   
01588     int boundarySegmentIndex() const
01589     {
01590         if(! boundary())
01591             DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
01592         update();
01593         // size of local macro grid
01594         const FieldVector<int, dim> & size = _inside.gridlevel().mg()->begin().cell_overlap().size();
01595         const FieldVector<int, dim> & origin = _inside.gridlevel().mg()->begin().cell_overlap().origin();
01596         FieldVector<int, dim> sides;
01597         {
01598             for (int i=0; i<dim; i++)
01599             {
01600                 sides[i] =
01601                     ((_inside.gridlevel().mg()->begin().cell_overlap().origin(i)
01602                         == _inside.gridlevel().mg()->begin().cell_global().origin(i))+
01603                         (_inside.gridlevel().mg()->begin().cell_overlap().origin(i) +
01604                             _inside.gridlevel().mg()->begin().cell_overlap().size(i)
01605                             == _inside.gridlevel().mg()->begin().cell_global().origin(i) +
01606                             _inside.gridlevel().mg()->begin().cell_global().size(i)));
01607             }
01608         }
01609         // gobal position of the cell on macro grid
01610         FieldVector<int, dim> pos = _inside.transformingsubiterator().coord();
01611         pos /= (1<<_inside.level());
01612         pos -= origin;
01613         // compute unit-cube-face-sizes
01614         FieldVector<int, dim> fsize;
01615         {
01616             int vol = 1;
01617             for (int k=0; k<dim; k++)
01618                 vol *= size[k];
01619             for (int k=0; k<dim; k++)
01620                 fsize[k] = vol/size[k];
01621         }
01622         // compute index in the unit-cube-face
01623         int index = 0;
01624         {
01625             int localoffset = 1;
01626             for (int k=dim-1; k>=0; k--)
01627             {
01628                 if (k == _dir) continue;
01629                 index += (pos[k]) * localoffset;
01630                 localoffset *= size[k];
01631             }
01632         }
01633         // add unit-cube-face-offsets
01634         {
01635             for (int k=0; k<_dir; k++)
01636                 index += sides[k] * fsize[k];
01637             // add fsize if we are on the right face and there is a left-face-boundary on this processor
01638             index += _face * (sides[_dir]>1) * fsize[_dir];
01639         }
01640 
01641         // int rank = 0;
01642         // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
01643         // std::cout << rank << "... size: " << size << " sides: " << sides
01644         //           << " fsize: " << fsize
01645         //           << " pos: " << pos << " face: " << int(_dir) << "/" << int(_face)
01646         //           << " index: " << index << std::endl;
01647         
01648         return index;
01649     }
01650   
01652     FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dim-1>& local) const
01653     {
01654         return _faceInfo[_count].normal;
01655     }
01656 
01658     FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dim-1>& local) const
01659     {
01660         return _faceInfo[_count].normal;
01661     }
01662 
01664     FieldVector<ctype, dimworld> centerUnitOuterNormal () const
01665     {
01666         return _faceInfo[_count].normal;
01667     }
01668 
01672     FieldVector<ctype, dimworld> integrationOuterNormal (const FieldVector<ctype, dim-1>& local) const
01673     {
01674         FieldVector<ctype, dimworld> n = _faceInfo[_count].normal;
01675         n *= geometry().volume();
01676         return n;
01677     }
01678  
01682     LocalGeometry geometryInInside () const
01683     {
01684         return LocalGeometry( _faceInfo[_count].geom_inside );
01685     }
01686 
01690     LocalGeometry geometryInOutside () const
01691     {
01692         return LocalGeometry( _faceInfo[_count].geom_outside );
01693     }
01694 
01697     Geometry geometry () const
01698     {
01699         update();
01700         GeometryImpl
01701           _is_global(_pos_world,_inside.transformingsubiterator().meshsize(),_dir);
01702         return Geometry( _is_global );
01703     }
01704 
01706     GeometryType type () const
01707     {
01708         static const GeometryType cube(GeometryType::cube, dim-1);
01709         return cube;
01710     }
01711 
01713     int indexInInside () const
01714     {
01715         return _count;
01716     }
01717 
01719     int indexInOutside () const
01720     {
01721         // flip the last bit
01722         return _count^1;
01723     }
01724 
01726     YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
01727         _inside(myself.yaspgrid(), myself.gridlevel(),
01728             myself.transformingsubiterator()),
01729         _outside(myself.yaspgrid(), myself.gridlevel(),
01730             myself.transformingsubiterator()),
01731         // initialize to first neighbor
01732         _count(0),
01733         _dir(0),
01734         _face(0),
01735         _pos_world(myself.transformingsubiterator().position())
01736     {
01737         if (toend)
01738         {
01739             // initialize end iterator
01740             _count = 2*dim;
01741             return;
01742         }
01743         _count = 0;
01744         
01745         // move transforming iterator
01746         _outside.transformingsubiterator().move(_dir,-1);
01747 
01748         // make up faces
01749         _pos_world[0] -= 0.5*_inside.transformingsubiterator().meshsize(0);
01750     }
01751 
01753     YaspIntersection (const YaspIntersection& it) :
01754         _inside(it._inside),
01755         _outside(it._outside),
01756         _count(it._count),
01757         _dir(it._dir),
01758         _face(it._face),
01759         _pos_world(it._pos_world)
01760     {}
01761 
01763     void assign (const YaspIntersection& it)
01764     {
01765         _inside = it._inside;
01766         _outside = it._outside;
01767         _count = it._count;
01768         _dir = it._dir;
01769         _face = it._face;
01770         _pos_world = it._pos_world;
01771     }
01772 
01773 private:
01774     /* EntityPointers (get automatically updated) */
01775     mutable YaspEntityPointer<0,GridImp> _inside;  
01776     mutable YaspEntityPointer<0,GridImp> _outside; 
01777     /* current position */
01778     uint8_t _count;                                
01779     mutable uint8_t _dir;                          
01780     mutable uint8_t _face;                         
01781     /* current position */
01782     mutable FieldVector<ctype, dimworld> _pos_world;       
01783 
01784     /* static data */
01785     struct faceInfo
01786     {
01787         FieldVector<ctype, dimworld> normal;
01788         LocalGeometryImpl       geom_inside;   
01789         LocalGeometryImpl       geom_outside;  
01790     };
01791 
01792     /* static face info */
01793     static const array<faceInfo, 2*GridImp::dimension> _faceInfo;
01794     
01795     static array<faceInfo, 2*dim> initFaceInfo()
01796     {
01797         const FieldVector<typename GridImp::ctype, GridImp::dimension> ext_local(1.0);
01798         array<faceInfo, 2*dim> I;
01799         for (uint8_t i=0; i<dim; i++)
01800         {
01801             // center of face
01802             FieldVector<ctype, dim> a(0.5); a[i] = 0.0;
01803             FieldVector<ctype, dim> b(0.5); b[i] = 1.0;
01804             // normal vectors
01805             I[2*i].normal = 0.0;
01806             I[2*i+1].normal = 0.0;
01807             I[2*i].normal[i] = -1.0;
01808             I[2*i+1].normal[i] = +1.0;
01809             // geometries
01810             I[2*i].geom_inside =
01811               LocalGeometryImpl(a, ext_local, i);
01812             I[2*i].geom_outside =
01813               LocalGeometryImpl(b, ext_local, i);
01814             I[2*i+1].geom_inside =
01815               LocalGeometryImpl(b, ext_local, i);
01816             I[2*i+1].geom_outside =
01817               LocalGeometryImpl(a, ext_local, i);
01818         }
01819         return I;
01820     }
01821 };
01822 
01823 template<class GridImp>
01824 const array<typename YaspIntersection<GridImp>::faceInfo, 2*GridImp::dimension>
01825 YaspIntersection<GridImp>::_faceInfo =
01826     YaspIntersection<GridImp>::initFaceInfo();
01827 
01828 //========================================================================
01833 //========================================================================
01834 
01835 template<class GridImp>
01836 class YaspIntersectionIterator : public MakeableInterfaceObject< Dune::Intersection<const GridImp, Dune::YaspIntersection > >
01837 {
01838     enum { dim=GridImp::dimension };
01839     enum { dimworld=GridImp::dimensionworld };
01840     typedef typename GridImp::ctype ctype;
01841     YaspIntersectionIterator();
01842 public:
01843     // types used from grids
01844     typedef Dune::Intersection<const GridImp, Dune::YaspIntersection> Intersection;
01845     typedef MakeableInterfaceObject<Intersection> MakeableIntersection;
01846   
01848     void increment()
01849     {
01850         GridImp::getRealImplementation(*this).increment();
01851     }
01852 
01854     bool equals (const YaspIntersectionIterator& other) const
01855     {
01856         return GridImp::getRealImplementation(*this).equals(
01857             GridImp::getRealImplementation(other));
01858     }
01859 
01861     const Intersection & dereference() const
01862     {
01863         return *this;
01864     }
01865   
01867     YaspIntersectionIterator (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
01868         MakeableIntersection(YaspIntersection<GridImp>(myself, toend))
01869     {}
01870 
01872     YaspIntersectionIterator (const YaspIntersectionIterator& it) :
01873         MakeableIntersection(it)
01874     {}
01875 
01877     YaspIntersectionIterator & operator = (const YaspIntersectionIterator& it)
01878     {
01879         GridImp::getRealImplementation(*this).assign(
01880             GridImp::getRealImplementation(it));
01881         return *this;
01882     }
01883 };
01884 
01885 //========================================================================
01889 //========================================================================
01890 
01891 template<class GridImp>
01892 class YaspHierarchicIterator :
01893   public YaspEntityPointer<0,GridImp>
01894 {
01895   enum { dim=GridImp::dimension };
01896   enum { dimworld=GridImp::dimensionworld };  
01897   typedef typename GridImp::ctype ctype;
01898 public:
01899   // types used from grids
01900   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
01901   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
01902   typedef typename GridImp::template Codim<0>::Entity Entity;
01903   typedef YaspSpecialEntity<0,dim,GridImp> SpecialEntity;
01904 
01906   typedef typename YGrid<dim,ctype>::iTupel iTupel;
01907 
01909   YaspHierarchicIterator (const GridImp* yg, const YGLI& g, const TSI& it, int maxlevel) :
01910     YaspEntityPointer<0,GridImp>(yg,g,it)
01911   {
01912       // now iterator points to current cell
01913       StackElem se(this->_g);
01914       se.coord = this->_it.coord();
01915       stack.push(se);
01916 
01917       // determine maximum level
01918       _maxlevel = std::min(maxlevel,this->_g.mg()->maxlevel());
01919 
01920       // if maxlevel not reached then push yourself and sons
01921       if (this->_g.level()<_maxlevel)
01922       {
01923         push_sons();
01924       }
01925 
01926       // and make iterator point to first son if stack is not empty
01927       if (!stack.empty())
01928         pop_tos();
01929   }
01930 
01932   YaspHierarchicIterator (const YaspHierarchicIterator& it) :
01933     YaspEntityPointer<0,GridImp>(it),
01934     _maxlevel(it._maxlevel), stack(it.stack)
01935   {}
01936 
01938   void increment ()
01939   {
01940         // sanity check: do nothing when stack is empty
01941         if (stack.empty()) return;
01942 
01943         // if maxlevel not reached then push sons
01944         if (this->_g.level()<_maxlevel)
01945           push_sons();
01946 
01947         // in any case pop one element
01948         pop_tos();
01949   }
01950 
01951   void print (std::ostream& s) const
01952   {
01953         s << "HIER: " << "level=" << this->_g.level()
01954           << " position=" << this->_it.coord()
01955           << " superindex=" << this->_it.superindex()
01956           << " maxlevel=" << this->_maxlevel
01957           << " stacksize=" << stack.size()
01958           << std::endl;
01959   }
01960 
01961 private:
01962   int _maxlevel;         
01963 
01964   struct StackElem {
01965         YGLI g;       // grid level of the element
01966         iTupel coord; // and the coordinates
01967         StackElem(YGLI gg) : g(gg) {}
01968   };
01969   std::stack<StackElem> stack;      
01970 
01971   // push sons of current element on the stack
01972   void push_sons ()
01973   {
01974         // yes, process all 1<<dim sons
01975         StackElem se(this->_g.finer());
01976         for (int i=0; i<(1<<dim); i++)
01977           {
01978                 for (int k=0; k<dim; k++)
01979                   if (i&(1<<k))
01980                         se.coord[k] = this->_it.coord(k)*2+1;
01981                   else
01982                         se.coord[k] = this->_it.coord(k)*2;
01983                 stack.push(se);
01984           }
01985   }
01986 
01987   // make TOS the current element
01988   void pop_tos ()
01989   {
01990         StackElem se = stack.top();
01991         stack.pop();
01992         this->_g = se.g;
01993         this->_it.reinit(this->_g.cell_overlap(),se.coord);
01994   }
01995 };
01996 
01997 //========================================================================
02001 //========================================================================
02002 template<int codim, class GridImp>
02003 class YaspEntitySeed
02004 {
02006   enum { dim=GridImp::dimension };
02007 
02008 public:
02010   enum { codimension = codim };
02011   
02013   YaspEntitySeed (int level, FieldVector<int, dim> coord)
02014     : _l(level), _c(coord)
02015   {}
02016 
02018   YaspEntitySeed (const YaspEntitySeed& rhs)
02019     : _l(rhs._l), _c(rhs._c)
02020   {}
02021 
02022   int level () const { return _l; }
02023   const FieldVector<int, dim> & coord() const { return _c; }
02024 
02025 protected:
02026   int _l;                    // grid level
02027   FieldVector<int, dim> _c;  // coord in the global grid
02028 };
02029 
02030 //========================================================================
02040 //========================================================================
02041 template<int codim, class GridImp>
02042 class YaspEntityPointer
02043 {
02045   enum { dim=GridImp::dimension };
02047   enum { dimworld=GridImp::dimensionworld };
02048   typedef typename GridImp::ctype ctype;
02049 public:
02050   typedef typename GridImp::template Codim<codim>::Entity Entity;
02051   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
02052   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
02053   typedef YaspSpecialEntity<codim,dim,GridImp> SpecialEntity;
02054   typedef YaspEntityPointer<codim,GridImp> EntityPointerImp;
02055 protected:  
02056   typedef YaspEntity<codim, dim, GridImp> YaspEntityImp; 
02057 
02058 public:
02060   enum { codimension = codim };
02061   
02063   YaspEntityPointer (const GridImp * yg, const YGLI & g, const TSI & it)
02064     : _g(g), _it(it), _entity(yg, _g,_it)
02065   {
02066         if (codim>0 && codim<dim)
02067           {
02068                 DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
02069           }
02070   }
02071 
02073   YaspEntityPointer (const YaspEntityImp& entity) 
02074     : _g(entity.gridlevel()), 
02075       _it(entity.transformingsubiterator()), 
02076       _entity(entity.yaspgrid(),_g,_it)
02077   {
02078         if (codim>0 && codim<dim)
02079           {
02080                 DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
02081           }
02082   }
02083 
02085   YaspEntityPointer (const YaspEntityPointer& rhs)
02086     : _g(rhs._g), _it(rhs._it), _entity(rhs._entity.yaspgrid(),_g,_it)
02087   {
02088         if (codim>0 && codim<dim)
02089           {
02090                 DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
02091           }
02092   }
02093 
02095   bool equals (const YaspEntityPointer& rhs) const
02096   {
02097         return (_it==rhs._it && _g == rhs._g);
02098   }
02099 
02101   Entity& dereference() const
02102   {
02103         return _entity;
02104   }
02105 
02107   int level () const {return _g.level();}
02108 
02109   const YaspEntityPointer&
02110   operator = (const YaspEntityPointer& rhs)
02111     {
02112       _g = rhs._g;
02113       _it = rhs._it;
02114       /* _entity = i._entity
02115        * is done implicitely, as the entity is completely
02116        * defined via the interator it belongs to
02117        */
02118       return *this;
02119     }  
02120   
02121   const TSI& transformingsubiterator () const
02122   {
02123         return _it;
02124   }
02125 
02126   const YGLI& gridlevel () const
02127   {
02128         return _g;
02129   }
02130 
02131   TSI& transformingsubiterator ()
02132   {
02133         return _it;
02134   }
02135 
02136   YGLI& gridlevel ()
02137   {
02138         return _g;
02139   }
02140 
02141 protected:
02142   YGLI _g;               // access to grid level
02143   TSI _it;               // position in the grid level
02144   mutable SpecialEntity _entity; 
02145 };
02146 
02147 //========================================================================
02155 //========================================================================
02156 
02157 template<int codim, PartitionIteratorType pitype, class GridImp>
02158 class YaspLevelIterator :
02159   public YaspEntityPointer<codim,GridImp>
02160 {
02162   enum { dim=GridImp::dimension };
02164   enum { dimworld=GridImp::dimensionworld };
02165   typedef typename GridImp::ctype ctype;
02166 public:
02167   typedef typename GridImp::template Codim<codim>::Entity Entity;
02168   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
02169   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
02170   typedef YaspSpecialEntity<codim,dim,GridImp> SpecialEntity;
02171 
02173   YaspLevelIterator (const GridImp * yg, const YGLI & g, const TSI & it) :
02174     YaspEntityPointer<codim,GridImp>(yg,g,it) {}
02175 
02177   YaspLevelIterator (const YaspLevelIterator& i) :
02178     YaspEntityPointer<codim,GridImp>(i) {}
02179 
02181   void increment()
02182   {
02183         ++(this->_it);
02184   }
02185 };
02186 
02187 
02188 //========================================================================
02193 //========================================================================
02194 
02195 template<class GridImp>
02196 class YaspLevelIndexSet
02197   : public IndexSet< GridImp, YaspLevelIndexSet< GridImp >, unsigned int >
02198 {
02199   typedef YaspLevelIndexSet< GridImp > This;
02200   typedef IndexSet< GridImp, This, unsigned int > Base;
02201 
02202 public:
02203   typedef typename Base::IndexType IndexType;
02204 
02205   using Base::subIndex;
02206   
02208   YaspLevelIndexSet ( const GridImp &g, int l )
02209   : grid( g ),
02210     level( l )
02211   {
02212     // contains a single element type;
02213     for (int codim=0; codim<=GridImp::dimension; codim++)
02214       mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
02215   }
02216 
02218   template<int cc>
02219   IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const 
02220   {
02221     assert( cc == 0 || cc == GridImp::dimension );
02222     return grid.getRealImplementation(e).compressedIndex(); 
02223   }
02224 
02226   template< int cc >
02227   IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
02228                        int i, unsigned int codim ) const
02229   {
02230     assert( cc == 0 || cc == GridImp::dimension );
02231     if( cc == GridImp::dimension )
02232       return grid.getRealImplementation(e).compressedIndex(); 
02233     else
02234       return grid.getRealImplementation(e).subCompressedIndex(i,codim);
02235   }
02236 
02238   int size (GeometryType type) const
02239   {
02240     return grid.size( level, type );
02241   }
02242 
02244   int size (int codim) const
02245   {
02246     return grid.size( level, codim );
02247   }
02248   
02250   template<class EntityType>
02251   bool contains (const EntityType& e) const
02252   {
02253     return e.level() == level;
02254   }
02255   
02257   const std::vector<GeometryType>& geomTypes (int codim) const
02258   {
02259     return mytypes[codim];
02260   }
02261 
02262 private:
02263   const GridImp& grid;
02264   int level;
02265   std::vector<GeometryType> mytypes[GridImp::dimension+1];
02266 };
02267 
02268 
02269 // Leaf Index Set
02270 
02271 template<class GridImp>
02272 class YaspLeafIndexSet
02273   : public IndexSet< GridImp, YaspLeafIndexSet< GridImp >, unsigned int >
02274 {
02275   typedef YaspLeafIndexSet< GridImp > This;
02276   typedef IndexSet< GridImp, This, unsigned int > Base;
02277 
02278 public:
02279   typedef typename Base::IndexType IndexType;
02280 
02281   using Base::subIndex;
02282   
02284   explicit YaspLeafIndexSet ( const GridImp &g )
02285   : grid( g )
02286   {
02287     // contains a single element type;
02288     for (int codim=0; codim<=GridImp::dimension; codim++)
02289       mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
02290   }
02291 
02293   template<int cc>
02294   IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const 
02295   {
02296     assert( cc == 0 || cc == GridImp::dimension );
02297     return grid.getRealImplementation(e).compressedIndex(); 
02298   }
02299 
02301   template< int cc >
02302   IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
02303                        int i, unsigned int codim ) const
02304   {
02305     assert( cc == 0 || cc == GridImp::dimension );
02306     if( cc == GridImp::dimension )
02307       return grid.getRealImplementation(e).compressedIndex(); 
02308     else
02309       return grid.getRealImplementation(e).subCompressedIndex(i,codim);
02310   }
02311 
02313   int size (GeometryType type) const
02314   {
02315     return grid.size( grid.maxLevel(), type );
02316   }
02317 
02319   int size (int codim) const
02320   {
02321     return grid.size( grid.maxLevel(), codim );
02322   }
02323 
02325   template<class EntityType>
02326   bool contains (const EntityType& e) const
02327   {
02328     //return e.isLeaf();
02329     return (e.level() == grid.maxLevel());
02330   }
02331   
02333   const std::vector<GeometryType>& geomTypes (int codim) const
02334   {
02335     return mytypes[codim];
02336   }
02337 
02338 private:
02339   const GridImp& grid;
02340   enum { ncodim = remove_const<GridImp>::type::dimension+1 };
02341   std::vector<GeometryType> mytypes[ncodim];
02342 };
02343 
02344 
02345 
02346 
02347 //========================================================================
02352 //========================================================================
02353 
02354 template<class GridImp>
02355 class YaspGlobalIdSet : public IdSet<GridImp,YaspGlobalIdSet<GridImp>,
02356                                      typename remove_const<GridImp>::type::PersistentIndexType >
02357 /*
02358   We used the remove_const to extract the Type from the mutable class,
02359   because the const class is not instantiated yet.
02360 */
02361 {
02362   typedef YaspGlobalIdSet< GridImp > This;
02363 
02364 public:
02366   typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;
02367 
02368   using IdSet<GridImp, This, IdType>::subId;
02369 
02371   explicit YaspGlobalIdSet ( const GridImp &g )
02372   : grid( g )
02373   {}
02374 
02376   /*
02377     We use the remove_const to extract the Type from the mutable class,
02378     because the const class is not instantiated yet.
02379   */
02380   template<int cd>
02381   IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 
02382   {
02383     return grid.getRealImplementation(e).persistentIndex();
02384   }
02385 
02387   /*
02388     We use the remove_const to extract the Type from the mutable class,
02389     because the const class is not instantiated yet.
02390   */
02391   IdType subId (const typename remove_const<GridImp>::type::Traits::template Codim< 0 >::Entity &e,
02392                  int i, unsigned int codim ) const
02393   {
02394     return grid.getRealImplementation(e).subPersistentIndex(i,codim);
02395   }
02396 
02397 private:
02398   const GridImp& grid;
02399 };
02400 
02401 
02402 template<int dim, int dimworld>
02403 struct YaspGridFamily
02404 {
02405 #if HAVE_MPI
02406   typedef CollectiveCommunication<MPI_Comm> CCType;
02407 #else
02408   typedef CollectiveCommunication<Dune::YaspGrid<dim> > CCType;
02409 #endif
02410 
02411   typedef GridTraits<dim,dimworld,Dune::YaspGrid<dim>,
02412                      YaspGeometry,YaspEntity,
02413                      YaspEntityPointer,YaspLevelIterator,
02414                      YaspIntersection, // leaf  intersection 
02415                      YaspIntersection, // level intersection 
02416                      YaspIntersectionIterator, // leaf  intersection iter 
02417                      YaspIntersectionIterator, // level intersection iter 
02418                      YaspHierarchicIterator,
02419                      YaspLevelIterator,
02420                      YaspLevelIndexSet< const YaspGrid< dim > >,
02421                      YaspLeafIndexSet< const YaspGrid< dim > >,
02422                      YaspGlobalIdSet<const YaspGrid<dim> >, 
02423                      bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
02424                      YaspGlobalIdSet<const YaspGrid<dim> >, 
02425                      bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
02426                      CCType,
02427                      DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
02428                      YaspEntitySeed>
02429   Traits;  
02430 };
02431 
02432 template<int dim, int codim>
02433 struct YaspCommunicateMeta {
02434   template<class G, class DataHandle>
02435   static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
02436   {
02437     if (data.contains(dim,codim))
02438       {
02439         DUNE_THROW(GridError, "interface communication not implemented");
02440       }
02441     YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
02442   }
02443 };
02444 
02445 template<int dim>
02446 struct YaspCommunicateMeta<dim,dim> {
02447   template<class G, class DataHandle>
02448   static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
02449   {
02450     if (data.contains(dim,dim))
02451       g.template communicateCodim<DataHandle,dim>(data,iftype,dir,level);
02452     YaspCommunicateMeta<dim,dim-1>::comm(g,data,iftype,dir,level);
02453   }
02454 };
02455 
02456 template<int dim>
02457 struct YaspCommunicateMeta<dim,0> {
02458   template<class G, class DataHandle>
02459   static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
02460   {
02461     if (data.contains(dim,0))
02462       g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
02463   }
02464 };
02465 
02466 
02467 //************************************************************************
02483 template<int dim>
02484 class YaspGrid :
02485   public GridDefaultImplementation<dim,dim,yaspgrid_ctype,YaspGridFamily<dim,dim> >,
02486   public MultiYGrid<dim,yaspgrid_ctype>
02487 {
02488   typedef const YaspGrid<dim> GridImp;
02489 
02490   void init()
02491   {
02492     setsizes();
02493     indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,0) );
02494     theleafindexset.push_back( new YaspLeafIndexSet<const YaspGrid<dim> >(*this) );
02495     theglobalidset.push_back( new YaspGlobalIdSet<const YaspGrid<dim> >(*this) );
02496     boundarysegmentssize();
02497   }
02498 
02499   void boundarysegmentssize()
02500   {
02501       // sizes of local macro grid
02502       const FieldVector<int, dim> & size = MultiYGrid<dim,ctype>::begin().cell_overlap().size();
02503       FieldVector<int, dim> sides;
02504       {
02505           for (int i=0; i<dim; i++)
02506           {
02507               sides[i] =
02508                   ((MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i)
02509                       == MultiYGrid<dim,ctype>::begin().cell_global().origin(i))+
02510                       (MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i) +
02511                           MultiYGrid<dim,ctype>::begin().cell_overlap().size(i)
02512                           == MultiYGrid<dim,ctype>::begin().cell_global().origin(i) +
02513                           MultiYGrid<dim,ctype>::begin().cell_global().size(i)));
02514           }
02515       }
02516       nBSegments = 0;
02517       for (int k=0; k<dim; k++)
02518       {
02519           int offset = 1;
02520           for (int l=0; l<dim; l++)
02521           {
02522               if (l==k) continue;
02523               offset *= size[l];
02524           }
02525           nBSegments += sides[k]*offset;
02526       }
02527   }
02528 
02529 public:
02530 
02531   using MultiYGrid<dim,yaspgrid_ctype>::defaultLoadbalancer;
02532 
02534   typedef yaspgrid_ctype ctype;
02535 
02536   // define the persistent index type
02537   typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits> PersistentIndexType;
02538 
02540   typedef YaspGridFamily<dim,dim> GridFamily;
02541   // the Traits 
02542   typedef typename YaspGridFamily<dim,dim>::Traits Traits;
02543 
02544   // need for friend declarations in entity
02545   typedef YaspLevelIndexSet<YaspGrid<dim> > LevelIndexSetType;
02546   typedef YaspLeafIndexSet<YaspGrid<dim> > LeafIndexSetType;
02547   typedef YaspGlobalIdSet<YaspGrid<dim> > GlobalIdSetType;
02548 
02550   enum { MAXL=64 };
02551 
02553   typedef MultiYGrid<dim,ctype> YMG;
02554   typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
02555   typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
02556   typedef typename MultiYGrid<dim,ctype>::Intersection IS;
02557   typedef typename std::deque<IS>::const_iterator ISIT;
02558 
02567   YaspGrid (Dune::MPIHelper::MPICommunicator comm,
02568             Dune::FieldVector<ctype, dim> L, 
02569             Dune::FieldVector<int, dim> s, 
02570             Dune::FieldVector<bool, dim> periodic, int overlap,
02571             const YLoadBalance<dim>* lb = defaultLoadbalancer())
02572 #if HAVE_MPI
02573       : YMG(comm,L,s,periodic,overlap,lb), ccobj(comm),
02574         keep_ovlp(true), adaptRefCount(0), adaptActive(false)
02575 #else
02576       : YMG(L,s,periodic,overlap,lb),
02577         keep_ovlp(true), adaptRefCount(0), adaptActive(false)
02578 #endif
02579   {
02580     init();
02581   }
02582 
02583   
02596   YaspGrid (Dune::FieldVector<ctype, dim> L, 
02597             Dune::FieldVector<int, dim> s, 
02598             Dune::FieldVector<bool, dim> periodic, int overlap,
02599             const YLoadBalance<dim>* lb = YMG::defaultLoadbalancer())
02600 #if HAVE_MPI
02601       : YMG(MPI_COMM_SELF,L,s,periodic,overlap,lb), ccobj(MPI_COMM_SELF),
02602         keep_ovlp(true), adaptRefCount(0), adaptActive(false)
02603 #else
02604       : YMG(L,s,periodic,overlap,lb), 
02605         keep_ovlp(true), adaptRefCount(0), adaptActive(false)
02606 #endif
02607   {
02608     init();
02609   }
02610 
02611   ~YaspGrid()
02612   {
02613     deallocatePointers(indexsets);
02614     deallocatePointers(theleafindexset);
02615     deallocatePointers(theglobalidset);
02616   }
02617 
02618   private:
02619   // do not copy this class 
02620   YaspGrid(const YaspGrid&);
02621 
02622   public:
02623   
02627   int maxLevel() const {return MultiYGrid<dim,ctype>::maxlevel();} // delegate
02628 
02630   void globalRefine (int refCount)
02631   {
02632     if (refCount < -maxLevel())
02633       DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
02634                  "Coarsening " << -refCount << " levels requested!");
02635     for (int k=refCount; k<0; k++)
02636       {
02637         MultiYGrid<dim,ctype>::coarsen();
02638         setsizes();
02639         indexsets.pop_back();
02640       }
02641     for (int k=0; k<refCount; k++)
02642       {
02643         MultiYGrid<dim,ctype>::refine(keep_ovlp);
02644         setsizes();
02645         indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,maxLevel()) );
02646       }
02647   }
02648 
02653   void refineOptions (bool keepPhysicalOverlap)
02654   {
02655     keep_ovlp = keepPhysicalOverlap;
02656   }
02657 
02669   bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
02670   {
02671     assert(adaptActive == false);
02672     if (e.level() != maxLevel()) return false;
02673     adaptRefCount = std::max(adaptRefCount, refCount);
02674     return true;
02675   }
02676 
02683   int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
02684   {
02685     return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
02686   }
02687   
02689   bool adapt ()   
02690   {
02691     globalRefine(adaptRefCount);
02692     return (adaptRefCount > 0);
02693   }
02694   
02696   bool preAdapt ()
02697   {
02698     adaptActive = true;
02699     adaptRefCount = comm().max(adaptRefCount);
02700     return (adaptRefCount < 0);
02701   }
02702 
02704   void postAdapt()
02705   {
02706     adaptActive = false;
02707     adaptRefCount = 0;
02708   }
02709 
02711   template<int cd, PartitionIteratorType pitype>
02712   typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
02713     {
02714       return levelbegin<cd,pitype>(level);
02715     }
02716 
02718   template<int cd, PartitionIteratorType pitype>
02719   typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
02720     {
02721       return levelend<cd,pitype>(level);
02722     }
02723 
02725   template<int cd>
02726   typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
02727     {
02728       return levelbegin<cd,All_Partition>(level);
02729     }
02730 
02732   template<int cd>
02733   typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
02734     {
02735       return levelend<cd,All_Partition>(level);
02736     }
02737 
02739   template<int cd, PartitionIteratorType pitype>
02740   typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const
02741     {
02742       return levelbegin<cd,pitype>(maxLevel());
02743     }
02744 
02746   template<int cd, PartitionIteratorType pitype>
02747   typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const
02748     {
02749       return levelend<cd,pitype>(maxLevel());
02750     }
02751 
02753   template<int cd>
02754   typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
02755     {
02756       return levelbegin<cd,All_Partition>(maxLevel());
02757     }
02758 
02760   template<int cd>
02761   typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
02762     {
02763       return levelend<cd,All_Partition>(maxLevel());
02764     }
02765 
02766   // \brief obtain EntityPointer from EntitySeed. */
02767   template <typename Seed>
02768   typename Traits::template Codim<Seed::codimension>::EntityPointer
02769   entityPointer(const Seed& seed) const
02770   {
02771     static const int codim = Seed::codimension;
02772     YGLI g = MultiYGrid<dim,ctype>::begin(seed.level());
02773     switch (codim)
02774     {
02775       case 0:
02776         return YaspEntityPointer<codim,GridImp>(this,g,
02777           TSI(g.cell_overlap(), seed.coord()));
02778       case dim:
02779         return YaspEntityPointer<codim,GridImp>(this,g,
02780           TSI(g.vertex_overlap(), seed.coord()));
02781       default:
02782         DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
02783     }
02784   }
02785 
02787   int overlapSize (int level, int codim) const
02788   {
02789         YGLI g = MultiYGrid<dim,ctype>::begin(level);
02790         return g.overlap();
02791   }
02792 
02794   int overlapSize (int codim) const
02795   {
02796         YGLI g = MultiYGrid<dim,ctype>::begin(maxLevel());
02797         return g.overlap();
02798   }
02799 
02801   int ghostSize (int level, int codim) const
02802   {
02803         return 0;
02804   }
02805 
02807   int ghostSize (int codim) const
02808   {
02809       return 0;
02810   }
02811 
02813   int size (int level, int codim) const
02814   {
02815       return sizes[level][codim];
02816   }
02817 
02819   int size (int codim) const
02820   {
02821       return sizes[maxLevel()][codim];
02822   }
02823 
02825   int size (int level, GeometryType type) const
02826   {
02827       return (type.isCube()) ? sizes[level][dim-type.dim()] : 0;
02828   }
02829 
02831   int size (GeometryType type) const
02832   {
02833       return size(maxLevel(),type);
02834   }
02835 
02837   size_t numBoundarySegments () const
02838   {
02839     return nBSegments;
02840   }
02841 
02846   template<class DataHandleImp, class DataType>
02847   void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir, int level) const
02848   {
02849     YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
02850   }
02851 
02856   template<class DataHandleImp, class DataType>
02857   void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir) const
02858   {
02859     YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
02860   }
02861 
02866   template<class DataHandle, int codim>
02867   void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
02868   {
02869     // check input
02870     if (!data.contains(dim,codim)) return; // should have been checked outside
02871 
02872     // data types
02873     typedef typename DataHandle::DataType DataType;
02874 
02875     // access to grid level
02876     YGLI g = MultiYGrid<dim,ctype>::begin(level);
02877 
02878     // find send/recv lists or throw error
02879     const std::deque<IS>* sendlist=0;
02880     const std::deque<IS>* recvlist=0;
02881     if (codim==0) // the elements
02882       {
02883         if (iftype==InteriorBorder_InteriorBorder_Interface)
02884           return; // there is nothing to do in this case
02885         if (iftype==InteriorBorder_All_Interface)
02886           {
02887             sendlist = &g.send_cell_interior_overlap();
02888             recvlist = &g.recv_cell_overlap_interior();
02889           }
02890         if (iftype==Overlap_OverlapFront_Interface || iftype==Overlap_All_Interface || iftype==All_All_Interface)
02891           {
02892             sendlist = &g.send_cell_overlap_overlap();
02893             recvlist = &g.recv_cell_overlap_overlap();
02894           }
02895       }
02896     if (codim==dim) // the vertices
02897       {
02898         if (iftype==InteriorBorder_InteriorBorder_Interface)
02899           {
02900             sendlist = &g.send_vertex_interiorborder_interiorborder();
02901             recvlist = &g.recv_vertex_interiorborder_interiorborder();
02902           }
02903 
02904         if (iftype==InteriorBorder_All_Interface)
02905           {
02906             sendlist = &g.send_vertex_interiorborder_overlapfront();
02907             recvlist = &g.recv_vertex_overlapfront_interiorborder();
02908           }
02909         if (iftype==Overlap_OverlapFront_Interface || iftype==Overlap_All_Interface)
02910           {
02911             sendlist = &g.send_vertex_overlap_overlapfront();
02912             recvlist = &g.recv_vertex_overlapfront_overlap();
02913           }
02914         if (iftype==All_All_Interface)
02915           {
02916             sendlist = &g.send_vertex_overlapfront_overlapfront();
02917             recvlist = &g.recv_vertex_overlapfront_overlapfront();
02918           }
02919       }
02920 
02921     // change communication direction?
02922     if (dir==BackwardCommunication)
02923       std::swap(sendlist,recvlist);
02924 
02925     int cnt;
02926 
02927     // Size computation (requires communication if variable size)
02928     std::vector<int> send_size(sendlist->size(),-1);      // map rank to total number of objects (of type DataType) to be sent
02929     std::vector<int> recv_size(recvlist->size(),-1);      // map rank to total number of objects (of type DataType) to be recvd
02930     std::vector<size_t*> send_sizes(sendlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be sent
02931     std::vector<size_t*> recv_sizes(recvlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be recvd
02932     if (data.fixedsize(dim,codim))
02933       {
02934         // fixed size: just take a dummy entity, size can be computed without communication
02935         cnt=0;
02936         for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
02937           {
02938             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
02939               it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
02940             send_size[cnt] = is->grid.totalsize() * data.size(*it);
02941             cnt++;
02942           }
02943         cnt=0;
02944         for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
02945           {
02946             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
02947               it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
02948             recv_size[cnt] = is->grid.totalsize() * data.size(*it);
02949             cnt++;
02950           }
02951       }
02952     else
02953       {
02954         // variable size case: sender side determines the size
02955         cnt=0;
02956         for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
02957           {
02958             // allocate send buffer for sizes per entitiy
02959             size_t *buf = new size_t[is->grid.totalsize()];
02960             send_sizes[cnt] = buf;
02961 
02962             // loop over entities and ask for size
02963             int i=0; size_t n=0;
02964             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
02965               it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
02966             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
02967               tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
02968             for ( ; it!=tsubend; ++it)
02969               {
02970                 buf[i] = data.size(*it);
02971                 n += buf[i];
02972                 i++;
02973               }
02974 
02975             // now we know the size for this rank
02976             send_size[cnt] = n;
02977 
02978             // hand over send request to torus class
02979             MultiYGrid<dim,ctype>::torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
02980             cnt++;
02981           }
02982 
02983         // allocate recv buffers for sizes and store receive request
02984         cnt=0;
02985         for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
02986           {
02987             // allocate recv buffer
02988             size_t *buf = new size_t[is->grid.totalsize()];
02989             recv_sizes[cnt] = buf;
02990             
02991             // hand over recv request to torus class
02992             MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
02993             cnt++;
02994           }
02995 
02996         // exchange all size buffers now
02997         MultiYGrid<dim,ctype>::torus().exchange();
02998 
02999         // release send size buffers
03000         cnt=0;
03001         for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
03002           {
03003             delete[] send_sizes[cnt];
03004             send_sizes[cnt] = 0;
03005             cnt++;
03006           }
03007 
03008         // process receive size buffers
03009         cnt=0;
03010         for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
03011           {
03012             // get recv buffer
03013             size_t *buf = recv_sizes[cnt];
03014 
03015             // compute total size
03016             size_t n=0;
03017             for (int i=0; i<is->grid.totalsize(); ++i)
03018               n += buf[i];
03019 
03020             // ... and store it
03021             recv_size[cnt] = n;
03022             ++cnt;
03023           }
03024       }
03025 
03026 
03027     // allocate & fill the send buffers & store send request
03028     std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
03029     cnt=0;
03030     for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
03031       {
03032 //      std::cout << "[" << this->comm().rank() << "] "
03033 //                << " send " << " dest=" << is->rank
03034 //                << " size=" << send_size[cnt] 
03035 //                << std::endl;
03036 
03037         // allocate send buffer
03038         DataType *buf = new DataType[send_size[cnt]];
03039 
03040         // remember send buffer
03041         sends[cnt] = buf;
03042 
03043         // make a message buffer
03044         MessageBuffer<DataType> mb(buf);
03045 
03046         // fill send buffer; iterate over cells in intersection
03047         typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03048           it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
03049         typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03050           tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
03051         for ( ; it!=tsubend; ++it)
03052           data.gather(mb,*it);
03053 
03054         // hand over send request to torus class
03055         MultiYGrid<dim,ctype>::torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
03056         cnt++;
03057       }
03058 
03059     // allocate recv buffers and store receive request
03060     std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
03061     cnt=0;
03062     for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
03063       {
03064 //      std::cout << "[" << this->comm().rank() << "] "
03065 //                << " recv " << "  src=" << is->rank
03066 //                << " size=" << recv_size[cnt] 
03067 //                << std::endl;
03068 
03069         // allocate recv buffer
03070         DataType *buf = new DataType[recv_size[cnt]];
03071 
03072         // remember recv buffer
03073         recvs[cnt] = buf;
03074 
03075         // hand over recv request to torus class
03076         MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
03077         cnt++;
03078       }
03079 
03080     // exchange all buffers now
03081     MultiYGrid<dim,ctype>::torus().exchange();
03082 
03083     // release send buffers
03084     cnt=0;
03085     for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
03086     {
03087       delete[] sends[cnt];
03088       sends[cnt] = 0;
03089       cnt++;
03090     }
03091 
03092     // process receive buffers and delete them
03093     cnt=0;
03094     for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
03095       {
03096         // get recv buffer
03097         DataType *buf = recvs[cnt];
03098 
03099         // make a message buffer
03100         MessageBuffer<DataType> mb(buf);
03101 
03102         // copy data from receive buffer; iterate over cells in intersection
03103         if (data.fixedsize(dim,codim))
03104           {
03105             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03106               it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
03107             size_t n=data.size(*it);
03108             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03109               tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
03110             for ( ; it!=tsubend; ++it)
03111               data.scatter(mb,*it,n);
03112           }
03113         else
03114           {
03115             int i=0;
03116             size_t *sbuf = recv_sizes[cnt];
03117             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03118               it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
03119             typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
03120               tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
03121             for ( ; it!=tsubend; ++it)
03122               data.scatter(mb,*it,sbuf[i++]);
03123             delete[] sbuf;
03124           }
03125 
03126         // delete buffer
03127         delete[] buf; // hier krachts !
03128         cnt++;
03129       }
03130   }
03131 
03132   // The new index sets from DDM 11.07.2005
03133   const typename Traits::GlobalIdSet& globalIdSet() const
03134   {
03135     return *(theglobalidset[0]);
03136   }
03137   
03138   const typename Traits::LocalIdSet& localIdSet() const
03139   {
03140     return *(theglobalidset[0]);
03141   }
03142 
03143   const typename Traits::LevelIndexSet& levelIndexSet(int level) const
03144   {
03145     if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
03146     return *(indexsets[level]);
03147   }
03148 
03149   const typename Traits::LeafIndexSet& leafIndexSet() const
03150   {
03151     return *(theleafindexset[0]);
03152   }
03153 
03154 #if HAVE_MPI
03155 
03157   const CollectiveCommunication<MPI_Comm>& comm () const
03158   {
03159     return ccobj;
03160   }
03161 #else
03162 
03164   const CollectiveCommunication<YaspGrid>& comm () const
03165   {
03166     return ccobj;
03167   }
03168 #endif
03169 
03170   YaspIntersectionIterator<const YaspGrid<dim> >&
03171   getRealIntersectionIterator(typename Traits::LevelIntersectionIterator& it)
03172   {
03173     return this->getRealImplementation(it);
03174   }
03175 
03176   const YaspIntersectionIterator<const YaspGrid<dim> >&
03177   getRealIntersectionIterator(const typename Traits::LevelIntersectionIterator& it) const
03178   {
03179     return this->getRealImplementation(it);
03180   }
03181 
03182 private:
03183 
03184 #if HAVE_MPI
03185   CollectiveCommunication<MPI_Comm> ccobj; 
03186 #else
03187   CollectiveCommunication<YaspGrid> ccobj; 
03188 #endif
03189 
03190   std::vector<YaspLevelIndexSet<const YaspGrid<dim> >*> indexsets;
03191   std::vector<YaspLeafIndexSet<const YaspGrid<dim> >*> theleafindexset;
03192   std::vector<YaspGlobalIdSet<const YaspGrid<dim> >*> theglobalidset;
03193   int nBSegments;
03194 
03195   // Index classes need access to the real entity
03196   friend class Dune::YaspLevelIndexSet<const Dune::YaspGrid<dim> >;
03197   friend class Dune::YaspLeafIndexSet<const Dune::YaspGrid<dim> >;
03198   friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim> >;
03199 
03200   friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim> >;
03201   friend class Dune::YaspIntersection<const Dune::YaspGrid<dim> >;
03202   friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim> >;
03203 
03204   template<class T>
03205   void deallocatePointers(T& container)
03206   {
03207     typedef typename T::iterator Iterator;
03208     
03209     for(Iterator entry=container.begin(); entry != container.end(); ++entry)
03210       delete (*entry);
03211   }
03212   
03213   template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
03214   friend class Entity;
03215 
03216   template<class DT>
03217   class MessageBuffer {
03218   public:
03219     // Constructor
03220     MessageBuffer (DT *p)
03221     {
03222       a=p;
03223       i=0;
03224       j=0;
03225     }
03226 
03227     // write data to message buffer, acts like a stream !
03228     template<class Y>
03229     void write (const Y& data)
03230     {
03231       dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
03232       a[i++] = data;
03233     }
03234     
03235     // read data from message buffer, acts like a stream !
03236     template<class Y>
03237     void read (Y& data) const 
03238     {
03239       dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
03240       data = a[j++];
03241     }
03242 
03243   private:
03244     DT *a;
03245     int i;
03246     mutable int j;
03247   };
03248 
03249   void setsizes ()
03250   {
03251     for (YGLI g=MultiYGrid<dim,ctype>::begin(); g!=MultiYGrid<dim,ctype>::end(); ++g)
03252       {
03253         // codim 0 (elements)
03254         sizes[g.level()][0] = 1;
03255         for (int i=0; i<dim; ++i) 
03256           sizes[g.level()][0] *= g.cell_overlap().size(i);
03257 
03258         // codim 1 (faces)
03259         if (dim>1)
03260           {
03261             sizes[g.level()][1] = 0;
03262             for (int i=0; i<dim; ++i) 
03263               {
03264                 int s=g.cell_overlap().size(i)+1;
03265                 for (int j=0; j<dim; ++j)
03266                   if (j!=i)
03267                     s *= g.cell_overlap().size(j);
03268                 sizes[g.level()][1] += s;
03269               }
03270           }
03271 
03272         // codim dim-1 (edges)
03273         if (dim>2)
03274           {
03275             sizes[g.level()][dim-1] = 0;
03276             for (int i=0; i<dim; ++i) 
03277               {
03278                 int s=g.cell_overlap().size(i);
03279                 for (int j=0; j<dim; ++j)
03280                   if (j!=i)
03281                     s *= g.cell_overlap().size(j)+1;
03282                 sizes[g.level()][dim-1] += s;
03283               }
03284           }
03285 
03286         // codim dim (vertices)
03287         sizes[g.level()][dim] = 1;
03288         for (int i=0; i<dim; ++i) 
03289           sizes[g.level()][dim] *= g.vertex_overlapfront().size(i);
03290       }
03291   }
03292 
03294   template<int cd, PartitionIteratorType pitype>
03295   YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
03296   {
03297         dune_static_assert( cd == dim || cd == 0 ,
03298           "YaspGrid only supports Entities with codim=dim and codim=0");
03299         YGLI g = MultiYGrid<dim,ctype>::begin(level);
03300         if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
03301         if (pitype==Ghost_Partition)            
03302             return levelend <cd, pitype> (level);
03303         if (cd==0) // the elements
03304           {
03305                 if (pitype<=InteriorBorder_Partition) 
03306                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubbegin());
03307                 if (pitype<=All_Partition)            
03308                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubbegin());
03309           }
03310         if (cd==dim) // the vertices
03311           {
03312                 if (pitype==Interior_Partition)       
03313                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubbegin());
03314                 if (pitype==InteriorBorder_Partition) 
03315                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubbegin());
03316                 if (pitype==Overlap_Partition)        
03317                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubbegin());
03318                 if (pitype<=All_Partition)            
03319                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubbegin());
03320           }
03321         DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
03322   }
03323 
03325   template<int cd, PartitionIteratorType pitype>
03326   YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
03327   {
03328         dune_static_assert( cd == dim || cd == 0 ,
03329           "YaspGrid only supports Entities with codim=dim and codim=0");
03330         YGLI g = MultiYGrid<dim,ctype>::begin(level);
03331         if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
03332         if (cd==0) // the elements
03333           {
03334                 if (pitype<=InteriorBorder_Partition) 
03335                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubend());
03336                 if (pitype<=All_Partition || pitype == Ghost_Partition)
03337                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubend());
03338           }
03339         if (cd==dim) // the vertices
03340           {
03341                 if (pitype==Interior_Partition)       
03342                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubend());
03343                 if (pitype==InteriorBorder_Partition) 
03344                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubend());
03345                 if (pitype==Overlap_Partition)        
03346                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubend());
03347                 if (pitype<=All_Partition || pitype == Ghost_Partition)
03348                   return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubend());
03349           }
03350         DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
03351   }
03352 
03353   int sizes[MAXL][dim+1]; // total number of entities per level and codim
03354   bool keep_ovlp;
03355   int adaptRefCount;
03356   bool adaptActive;
03357 };
03358 
03359 namespace Capabilities
03360 {
03361 
03373   template<int dim>
03374   struct hasSingleGeometryType< YaspGrid<dim> >
03375   {
03376     static const bool v = true;
03377     static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
03378   };
03379 
03383   template<int dim>
03384   struct isCartesian< YaspGrid<dim> >
03385   {
03386     static const bool v = true;
03387   };
03388 
03392   template<int dim>
03393   struct hasEntity< YaspGrid<dim>, 0 >
03394   {
03395     static const bool v = true;
03396   };
03397   
03401   template<int dim>
03402   struct hasEntity< YaspGrid<dim>, dim >
03403   {
03404     static const bool v = true;
03405   };
03406 
03407   template< int dim >
03408   struct canCommunicate< YaspGrid< dim >, 0 >
03409   {
03410     static const bool v = true;
03411   };
03412 
03413   template< int dim >
03414   struct canCommunicate< YaspGrid< dim >, dim >
03415   {
03416     static const bool v = true;
03417   };
03418 
03422   template<int dim>
03423   struct isParallel< YaspGrid<dim> >
03424   {
03425     static const bool v = true;
03426   };
03427 
03431   template<int dim>
03432   struct isLevelwiseConforming< YaspGrid<dim> >
03433   {
03434     static const bool v = true;
03435   };
03436 
03440   template<int dim>
03441   struct isLeafwiseConforming< YaspGrid<dim> >
03442   {
03443     static const bool v = true;
03444   };
03445 
03446 }
03447 
03448 } // end namespace
03449 
03450 
03451 #endif
03452