dune-grid  2.2.0
yaspgrid.hh
Go to the documentation of this file.
1 #ifndef DUNE_YASPGRID_HH
2 #define DUNE_YASPGRID_HH
3 
4 #include<iostream>
5 #include<vector>
6 #include<algorithm>
7 #include<stack>
8 
9 // either include stdint.h or provide fallback for uint8_t
10 #if HAVE_STDINT_H
11 #include<stdint.h>
12 #else
13 typedef unsigned char uint8_t;
14 #endif
15 
16 #include <dune/grid/common/grid.hh> // the grid base classes
17 #include <dune/grid/yaspgrid/grids.hh> // the yaspgrid base classes
18 #include <dune/grid/common/capabilities.hh> // the capabilities
19 #include <dune/common/misc.hh>
20 #include <dune/common/bigunsignedint.hh>
21 #include <dune/common/typetraits.hh>
22 #include <dune/common/collectivecommunication.hh>
23 #include <dune/common/mpihelper.hh>
24 #include <dune/geometry/genericgeometry/topologytypes.hh>
27 
28 
29 #if HAVE_MPI
30 #include <dune/common/mpicollectivecommunication.hh>
31 #endif
32 
40 namespace Dune {
41 
42 //************************************************************************
46 typedef double yaspgrid_ctype;
47 static const yaspgrid_ctype yasptolerance=1E-13; // tolerance in coordinate computations
48 
49 /* some sizes for building global ids
50  */
51 const int yaspgrid_dim_bits = 24; // bits for encoding each dimension
52 const int yaspgrid_level_bits = 6; // bits for encoding level number
53 const int yaspgrid_codim_bits = 4; // bits for encoding codimension
54 
55 
56 //************************************************************************
57 // forward declaration of templates
58 
59 template<int dim> class YaspGrid;
60 template<int mydim, int cdim, class GridImp> class YaspGeometry;
61 template<int codim, int dim, class GridImp> class YaspEntity;
62 template<int codim, class GridImp> class YaspEntityPointer;
63 template<int codim, class GridImp> class YaspEntitySeed;
64 template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
65 template<class GridImp> class YaspIntersectionIterator;
66 template<class GridImp> class YaspIntersection;
67 template<class GridImp> class YaspHierarchicIterator;
68 template<class GridImp> class YaspLevelIndexSet;
69 template<class GridImp> class YaspLeafIndexSet;
70 template<class GridImp> class YaspGlobalIdSet;
71 
72 namespace FacadeOptions
73 {
74 
75  template<int dim, int mydim, int cdim>
76  struct StoreGeometryReference<mydim, cdim, YaspGrid<dim>, YaspGeometry>
77  {
78  static const bool v = false;
79  };
80 
81  template<int dim, int mydim, int cdim>
82  struct StoreGeometryReference<mydim, cdim, const YaspGrid<dim>, YaspGeometry>
83  {
84  static const bool v = false;
85  };
86 
87 }
88 
89 //========================================================================
90 // The transformation describing the refinement rule
91 
92 template<int dim, class GridImp>
94  static const array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > _geo;
95  static array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > initSons()
96  {
97  array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > geo;
98  FieldVector<yaspgrid_ctype,dim> midpoint(0.25);
99  FieldVector<yaspgrid_ctype,dim> extension(0.5);
100  for (int i=0; i<(1<<dim); i++)
101  {
102  midpoint = 0.25;
103  for (int k=0; k<dim; k++)
104  {
105  if (i&(1<<k))
106  midpoint[k] = 0.75;
107  }
108  geo[i] = YaspGeometry<dim,dim,GridImp>(midpoint, extension);
109  }
110  return geo;
111  }
112 };
113 
114 template<int dim, class GridImp>
115 const array<YaspGeometry<dim,dim,GridImp>, (1<<dim)>
116 YaspFatherRelativeLocalElement<dim, GridImp>::_geo =
117  YaspFatherRelativeLocalElement<dim, GridImp>::initSons();
118 
119 //========================================================================
127 //========================================================================
128 
130 template<int mydim,int cdim, class GridImp>
131 class YaspGeometry : public GeometryDefaultImplementation<mydim,cdim,GridImp,YaspGeometry>
132 {
133 public:
135  typedef typename GridImp::ctype ctype;
136 
139  {
140  return GeometryType(GeometryType::cube,mydim);
141  }
142 
144  bool affine() const { return true; }
145 
147  int corners () const
148  {
149  return 1<<mydim;
150  }
151 
153  FieldVector< ctype, cdim > corner ( const int i ) const
154  {
155  assert( i >= 0 && i < (int) coord_.N() );
156  FieldVector<ctype, cdim>& c = coord_[i];
157  int bit=0;
158  for (int k=0; k<cdim; k++) // run over all directions in world
159  {
160  if (k==missing)
161  {
162  c[k] = midpoint[k];
163  continue;
164  }
165  //k is not the missing direction
166  if (i&(1<<bit)) // check whether bit is set or not
167  c[k] = midpoint[k]+0.5*extension[k]; // bit is 1 in i
168  else
169  c[k] = midpoint[k]-0.5*extension[k]; // bit is 0 in i
170  bit++; // we have processed a direction
171  }
172 
173  return c;
174  }
175 
177  FieldVector< ctype, cdim > center ( ) const
178  {
179  return midpoint;
180  }
181 
183  FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const
184  {
185  FieldVector<ctype, cdim> g;
186  int bit=0;
187  for (int k=0; k<cdim; k++)
188  if (k==missing)
189  g[k] = midpoint[k];
190  else
191  {
192  g[k] = midpoint[k] + (local[bit]-0.5)*extension[k];
193  bit++;
194  }
195  return g;
196  }
197 
199  FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const
200  {
201  FieldVector<ctype, mydim> l; // result
202  int bit=0;
203  for (int k=0; k<cdim; k++)
204  if (k!=missing)
205  {
206  l[bit] = (global[k]-midpoint[k])/extension[k] + 0.5;
207  bit++;
208  }
209  return l;
210  }
211 
213  ctype volume () const
214  {
215  ctype volume=1.0;
216  for (int k=0; k<cdim; k++)
217  if (k!=missing) volume *= extension[k];
218  return volume;
219  }
220 
223  ctype integrationElement (const FieldVector<ctype, mydim>& local) const
224  {
225  return volume();
226  }
227 
229  FieldMatrix<ctype,mydim,cdim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
230  {
231  JT = 0.0;
232  int k=0;
233  for (int i=0; i<cdim; ++i)
234  {
235  if (i != missing)
236  {
237  JT[k][i] = extension[i]; // set diagonal element
238  k++;
239  }
240  }
241  return JT;
242  }
244  FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
245  {
246  Jinv = 0.0;
247  int k=0;
248  for (int i=0; i<cdim; ++i)
249  {
250  if (i != missing)
251  {
252  Jinv[i][k] = 1.0/extension[i]; // set diagonal element
253  k++;
254  }
255  }
256  return Jinv;
257  }
258 
261 
263  YaspGeometry (const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, uint8_t& m)
264  : midpoint(p), extension(h), missing(m)
265  {
266  if (cdim!=mydim+1)
267  DUNE_THROW(GridError, "general YaspGeometry assumes cdim=mydim+1");
268  }
269 
271  YaspGeometry (const YaspGeometry& other)
272  : midpoint(other.midpoint),
273  extension(other.extension),
274  missing(other.missing)
275  {
276  }
277 
279  void print (std::ostream& s) const
280  {
281  s << "YaspGeometry<"<<mydim<<","<<cdim<< "> ";
282  s << "midpoint";
283  for (int i=0; i<cdim; i++)
284  s << " " << midpoint[i];
285  s << " extension";
286  for (int i=0; i<cdim; i++)
287  s << " " << extension[i];
288  s << " missing is " << missing;
289  }
290 
291  // const YaspGeometry<mydim,cdim,GridImp>&
292  // operator = (const YaspGeometry<mydim,cdim,GridImp>& g);
293 
294 private:
295  // the element is fully defined by its midpoint the extension
296  // in each direction and the missing direction.
297  // Note cdim == mydim+1
298 
299  FieldVector<ctype, cdim> midpoint; // the midpoint
300  FieldVector<ctype, cdim> extension; // the extension
301  uint8_t missing; // the missing, i.e. constant direction
302 
303  // In addition we need memory in order to return references.
304  // Possibly we should change this in the interface ...
305  mutable FieldMatrix<ctype, mydim, cdim> JT; // the transposed of the jacobian
306  mutable FieldMatrix<ctype, cdim, mydim> Jinv; // the transposed of the jacobian inverse
307  mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, cdim> coord_; // the coordinates
308 
309 };
310 
311 
312 
314 template<int mydim, class GridImp>
315 class YaspGeometry<mydim,mydim,GridImp> : public GeometryDefaultImplementation<mydim,mydim,GridImp,YaspGeometry>
316 {
317 public:
318  typedef typename GridImp::ctype ctype;
319 
322  {
323  return GeometryType(GeometryType::cube,mydim);
324  }
325 
327  bool affine() const { return true; }
328 
330  int corners () const
331  {
332  return 1<<mydim;
333  }
334 
336  const FieldVector<ctype, mydim>& operator[] (int i) const
337  {
338  return corner(i);
339  }
340 
342  FieldVector< ctype, mydim > corner ( const int i ) const
343  {
344  assert( i >= 0 && i < (int) coord_.N() );
345  FieldVector<ctype, mydim>& c = coord_[i];
346  for (int k=0; k<mydim; k++)
347  if (i&(1<<k))
348  c[k] = midpoint[k]+0.5*extension[k]; // kth bit is 1 in i
349  else
350  c[k] = midpoint[k]-0.5*extension[k]; // kth bit is 0 in i
351  return c;
352  }
353 
355  FieldVector< ctype, mydim > center ( ) const
356  {
357  return midpoint;
358  }
359 
361  FieldVector<ctype, mydim> global (const FieldVector<ctype, mydim>& local) const
362  {
363  FieldVector<ctype,mydim> g;
364  for (int k=0; k<mydim; k++)
365  g[k] = midpoint[k] + (local[k]-0.5)*extension[k];
366  return g;
367  }
368 
370  FieldVector<ctype, mydim> local (const FieldVector<ctype,mydim>& global) const
371  {
372  FieldVector<ctype, mydim> l; // result
373  for (int k=0; k<mydim; k++)
374  l[k] = (global[k]-midpoint[k])/extension[k] + 0.5;
375  return l;
376  }
377 
380  ctype integrationElement (const FieldVector<ctype, mydim>& local) const
381  {
382  return volume();
383  }
384 
386  ctype volume () const
387  {
388  ctype vol=1.0;
389  for (int k=0; k<mydim; k++) vol *= extension[k];
390  return vol;
391  }
392 
394  FieldMatrix<ctype,mydim,mydim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
395  {
396  for (int i=0; i<mydim; ++i)
397  {
398  JT[i] = 0.0; // set column to zero
399  JT[i][i] = extension[i]; // set diagonal element
400  }
401  return JT;
402  }
404  FieldMatrix<ctype,mydim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
405  {
406  for (int i=0; i<mydim; ++i)
407  {
408  Jinv[i] = 0.0; // set column to zero
409  Jinv[i][i] = 1.0/extension[i]; // set diagonal element
410  }
411  return Jinv;
412  }
413 
416 
418  YaspGeometry (const FieldVector<ctype, mydim>& p, const FieldVector<ctype, mydim>& h)
419  : midpoint(p), extension(h)
420  {}
421 
423  YaspGeometry (const YaspGeometry& other)
424  : midpoint(other.midpoint),
425  extension(other.extension)
426  {
427  }
428 
430  void print (std::ostream& s) const
431  {
432  s << "YaspGeometry<"<<mydim<<","<<mydim<< "> ";
433  s << "midpoint";
434  for (int i=0; i<mydim; i++)
435  s << " " << midpoint[i];
436  s << " extension";
437  for (int i=0; i<mydim; i++)
438  s << " " << extension[i];
439  }
440 
441  // const YaspGeometry<mydim,mydim,GridImp>&
442  // operator = (const YaspGeometry<mydim,mydim,GridImp>& g);
443 
444 private:
445  // the element is fully defined by midpoint and the extension
446  // in each direction. References are used because this information
447  // is known outside the element in many cases.
448  // Note mydim==cdim
449 
450  FieldVector<ctype, mydim> midpoint; // the midpoint
451  FieldVector<ctype, mydim> extension; // the extension
452 
453  // In addition we need memory in order to return references.
454  // Possibly we should change this in the interface ...
455  mutable FieldMatrix<ctype, mydim, mydim> Jinv,JT; // the transpose of the jacobian and its inverse inverse
456  mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, mydim> coord_; // the coordinates
457 };
458 
460 template<int cdim, class GridImp>
461 class YaspGeometry<0,cdim,GridImp> : public GeometryDefaultImplementation<0,cdim,GridImp,YaspGeometry>
462 {
463 public:
464  typedef typename GridImp::ctype ctype;
465 
468  {
470  }
471 
473  bool affine() const { return true; }
474 
476  int corners () const
477  {
478  return 1;
479  }
480 
482  const FieldVector<ctype, cdim>& operator[] (int i) const
483  {
484  return position;
485  }
486 
488  FieldVector< ctype, cdim > corner ( const int i ) const
489  {
490  return position;
491  }
492 
494  FieldVector< ctype, cdim > center ( ) const
495  {
496  return position;
497  }
498 
501  ctype integrationElement (const FieldVector<ctype, 0>& local) const
502  {
503  return 1.0;
504  }
505 
507  FieldMatrix<ctype,0,cdim>& jacobianTransposed (const FieldVector<ctype, 0>& local) const
508  {
509  static FieldMatrix<ctype,0,cdim> JT(0.0);
510  return JT;
511  }
513  FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
514  {
515  static FieldMatrix<ctype,cdim,0> Jinv(0.0);
516  return Jinv;
517  }
518 
521  {}
522 
524  explicit YaspGeometry ( const FieldVector< ctype, cdim > &p )
525  : position( p )
526  {}
527 
528  YaspGeometry ( const FieldVector< ctype, cdim > &p, const FieldVector< ctype, cdim > &, uint8_t &)
529  : position( p )
530  {}
531 
533  void print (std::ostream& s) const
534  {
535  s << "YaspGeometry<"<<0<<","<<cdim<< "> ";
536  s << "position " << position;
537  }
538 
539  // const YaspGeometry<0,cdim,GridImp>&
540  // operator = (const YaspGeometry<0,cdim,GridImp>& g);
541 
542 private:
543  FieldVector<ctype, cdim> position;
544 };
545 
546 // operator<< for all YaspGeometrys
547 template <int mydim, int cdim, class GridImp>
548 inline
549 std::ostream& operator<< (std::ostream& s, YaspGeometry<mydim,cdim,GridImp>& e)
550 {
551  e.print(s);
552  return s;
553 }
554 
555 //========================================================================
563 //========================================================================
564 
565 template<int codim, int dim, class GridImp>
567  public GridImp::template Codim<codim>::Entity
568 {
569 public:
570  typedef typename GridImp::ctype ctype;
571 
572  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
573  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
574 
575  YaspSpecialEntity(const GridImp* yg, const YGLI& g, const TSI& it) :
576  GridImp::template Codim<codim>::Entity (YaspEntity<codim, dim, GridImp>(yg,g,it))
577  {};
579  GridImp::template Codim<codim>::Entity (e)
580  {};
581  const TSI& transformingsubiterator () const
582  {
583  return this->realEntity.transformingsubiterator();
584  }
585  const YGLI& gridlevel () const
586  {
587  return this->realEntity.gridlevel();
588  }
589  const GridImp * yaspgrid() const
590  {
591  return this->realEntity.yaspgrid();
592  }
593 };
594 
595 template<int codim, int dim, class GridImp>
597 : public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
598 {
599 public:
600  typedef typename GridImp::ctype ctype;
601 
602  typedef typename GridImp::template Codim<codim>::Geometry Geometry;
604  int level () const
605  {
606  DUNE_THROW(GridError, "YaspEntity not implemented");
607  }
608 
610  int index () const
611  {
612  DUNE_THROW(GridError, "YaspEntity not implemented");
613  }
614 
617  {
618  DUNE_THROW(GridError, "YaspEntity not implemented");
619  }
620 
623  {
624  DUNE_THROW(GridError, "YaspEntity not implemented");
625  }
626 
627  const GridImp * yaspgrid() const
628  {
629  DUNE_THROW(GridError, "YaspEntity not implemented");
630  }
631 
632  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
633  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
634  YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
635  {
636  DUNE_THROW(GridError, "YaspEntity not implemented");
637  }
638 
639  // IndexSets needs access to the private index methods
640  friend class Dune::YaspLevelIndexSet<GridImp>;
641  friend class Dune::YaspLeafIndexSet<GridImp>;
642  friend class Dune::YaspGlobalIdSet<GridImp>;
643  typedef typename GridImp::PersistentIndexType PersistentIndexType;
644 
647  {
648  DUNE_THROW(GridError, "YaspEntity not implemented");
649  }
650 
652  int compressedIndex () const
653  {
654  DUNE_THROW(GridError, "YaspEntity not implemented");
655  }
656 
658  int compressedLeafIndex () const
659  {
660  DUNE_THROW(GridError, "YaspEntity not implemented");
661  }
662 };
663 
664 
665 // specialization for codim=0
666 template<int dim, class GridImp>
667 class YaspEntity<0,dim,GridImp>
668 : public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
669 {
670  enum { dimworld = GridImp::dimensionworld };
671 
672  typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
673 
674 public:
675  typedef typename GridImp::ctype ctype;
676 
677  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
678  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
679 
680  typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
681  typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
682 
683  template <int cd>
684  struct Codim
685  {
686  typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
687  };
688 
689  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
690  typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
691  typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
692  typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
693  typedef typename GridImp::LeafIntersectionIterator LeafIntersectionIterator;
694  typedef typename GridImp::HierarchicIterator HierarchicIterator;
695 
697  typedef typename GridImp::PersistentIndexType PersistentIndexType;
698 
700  typedef typename YGrid<dim,ctype>::iTupel iTupel;
701 
702  // constructor
703  YaspEntity (const GridImp * yg, const YGLI& g, const TSI& it)
704  : _yg(yg), _it(it), _g(g)
705  {
706  }
707 
709  int level () const { return _g.level(); }
710 
712  int index () const { return _it.superindex(); } // superindex works also for iteration over subgrids
713 
715  int globalIndex () const {
716  return _g.cell_global().index(_it.coord());
717  }
718 
722  EntitySeed seed () const {
723  return EntitySeed(_g.level(), _it.coord());
724  }
725 
728  {
729  if (_g.cell_interior().inside(_it.coord())) return InteriorEntity;
730  if (_g.cell_overlap().inside(_it.coord())) return OverlapEntity;
731  DUNE_THROW(GridError, "Impossible GhostEntity " << _it.coord() << "\t"
732  << _g.cell_interior().origin() << "/" << _g.cell_interior().size());
733  return GhostEntity;
734  }
735 
737  Geometry geometry () const {
738  // the element geometry
739  GeometryImpl _geometry(_it.position(),_it.meshsize());
740  return Geometry( _geometry );
741  }
742 
745  template<int cc> int count () const
746  {
747  if (cc==dim) return 1<<dim;
748  if (cc==1) return 2*dim;
749  if (cc==dim-1) return dim*(1<<(dim-1));
750  if (cc==0) return 1;
751  DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
752  }
753 
756  template<int cc>
757  typename Codim<cc>::EntityPointer subEntity (int i) const
758  {
759  dune_static_assert( cc == dim || cc == 0 ,
760  "YaspGrid only supports Entities with codim=dim and codim=0");
761  // coordinates of the cell == coordinates of lower left corner
762  if (cc==dim)
763  {
764  iTupel coord = _it.coord();
765 
766  // get corner from there
767  for (int k=0; k<dim; k++)
768  if (i&(1<<k)) (coord[k])++;
769 
770  return YaspEntityPointer<cc,GridImp>(_yg,_g,_g.vertex_overlapfront().tsubbegin(coord));
771  }
772  if (cc==0)
773  {
774  return YaspEntityPointer<cc,GridImp>(_yg,_g,_it);
775  }
776  DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
777  }
778 
780  EntityPointer father () const
781  {
782  // check if coarse level exists
783  if (_g.level()<=0)
784  DUNE_THROW(GridError, "tried to call father on level 0");
785 
786  // yes, get iterator to it
787  YGLI cg = _g.coarser();
788 
789  // coordinates of the cell
790  iTupel coord = _it.coord();
791 
792  // get coordinates on next coarser level
793  for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
794 
795  return YaspEntityPointer<0,GridImp>(_yg,cg,cg.cell_overlap().tsubbegin(coord));
796  }
797 
799  bool hasFather () const
800  {
801  return (_g.level()>0);
802  }
803 
813  LocalGeometry geometryInFather () const
814  {
815  // determine which son we are
816  int son = 0;
817  for (int k=0; k<dim; k++)
818  if (_it.coord(k)%2)
819  son += (1<<k);
820 
821  // configure one of the 2^dim transformations
823  }
824 
825  const TSI& transformingsubiterator () const
826  {
827  return _it;
828  }
829 
830  const YGLI& gridlevel () const
831  {
832  return _g;
833  }
834 
835  const GridImp* yaspgrid () const
836  {
837  return _yg;
838  }
839 
840  bool isLeaf() const
841  {
842  return (_g.level() == _g.mg()->maxlevel());
843  }
844 
847  bool isNew () const { return _yg->adaptRefCount > 0 && _g.mg()->maxlevel() < _g.level() + _yg->adaptRefCount; }
848 
851  bool mightVanish () const { return false; }
852  // { return _yg->adaptRefCount < 0 && _g.mg()->maxlevel() < _g.level() - _yg->adaptRefCount; }
853 
855  IntersectionIterator ibegin () const
856  {
857  return YaspIntersectionIterator<GridImp>(*this,false);
858  }
859 
861  LeafIntersectionIterator ileafbegin () const
862  {
863  // only if entity is leaf this iterator delivers intersections
864  return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
865  }
866 
868  LevelIntersectionIterator ilevelbegin () const
869  {
870  return ibegin();
871  }
872 
874  IntersectionIterator iend () const
875  {
876  return YaspIntersectionIterator<GridImp>(*this,true);
877  }
878 
880  LeafIntersectionIterator ileafend () const
881  {
882  return iend();
883  }
884 
886  LevelIntersectionIterator ilevelend () const
887  {
888  return iend();
889  }
890 
895  HierarchicIterator hbegin (int maxlevel) const
896  {
897  return YaspHierarchicIterator<GridImp>(_yg,_g,_it,maxlevel);
898  }
899 
901  HierarchicIterator hend (int maxlevel) const
902  {
903  return YaspHierarchicIterator<GridImp>(_yg,_g,_it,_g.level());
904  }
905 
906 private:
907  // IndexSets needs access to the private index methods
908  friend class Dune::YaspLevelIndexSet<GridImp>;
909  friend class Dune::YaspLeafIndexSet<GridImp>;
910  friend class Dune::YaspGlobalIdSet<GridImp>;
911 
914  {
915  // get size of global grid
916  const iTupel& size = _g.cell_global().size();
917 
918  // get coordinate correction for periodic boundaries
919  int coord[dim];
920  for (int i=0; i<dim; i++)
921  {
922  coord[i] = _it.coord(i);
923  if (coord[i]<0) coord[i] += size[i];
924  if (coord[i]>=size[i]) coord[i] -= size[i];
925  }
926 
927  // encode codim
929 
930  // encode level
931  id = id << yaspgrid_level_bits;
932  id = id+PersistentIndexType(_g.level());
933 
934 
935  // encode coordinates
936  for (int i=dim-1; i>=0; i--)
937  {
938  id = id << yaspgrid_dim_bits;
939  id = id+PersistentIndexType(coord[i]);
940  }
941 
942  return id;
943  }
944 
946  int compressedIndex () const
947  {
948  return _it.superindex();
949  }
950 
952  int compressedLeafIndex () const
953  {
954  return _it.superindex();
955  }
956 
958  PersistentIndexType subPersistentIndex (int i, int cc) const
959  {
960  if (cc==0)
961  return persistentIndex();
962 
963  // get position of cell, note that global origin is zero
964  // adjust for periodic boundaries
965  int coord[dim];
966  for (int k=0; k<dim; k++)
967  {
968  coord[k] = _it.coord(k);
969  if (coord[k]<0) coord[k] += _g.cell_global().size(k);
970  if (coord[k]>=_g.cell_global().size(k)) coord[k] -= _g.cell_global().size(k);
971  }
972 
973  if (cc==dim)
974  {
975  // transform to vertex coordinates
976  for (int k=0; k<dim; k++)
977  if (i&(1<<k)) (coord[k])++;
978 
979  // determine min number of trailing zeroes
980  int trailing = 1000;
981  for (int i=0; i<dim; i++)
982  {
983  // count trailing zeros
984  int zeros = 0;
985  for (int j=0; j<_g.level(); j++)
986  if (coord[i]&(1<<j))
987  break;
988  else
989  zeros++;
990  trailing = std::min(trailing,zeros);
991  }
992 
993  // determine the level of this vertex
994  int level = _g.level()-trailing;
995 
996  // encode codim
997  PersistentIndexType id(dim);
998 
999  // encode level
1000  id = id << yaspgrid_level_bits;
1001  id = id+PersistentIndexType(level);
1002 
1003  // encode coordinates
1004  for (int i=dim-1; i>=0; i--)
1005  {
1006  id = id << yaspgrid_dim_bits;
1007  id = id+PersistentIndexType(coord[i]>>trailing);
1008  }
1009 
1010  return id;
1011  }
1012 
1013  if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1014  {
1015  // Idea: Use the doubled grid to assign coordinates to faces
1016 
1017  // ivar is the direction that varies
1018  int ivar=i/2;
1019 
1020  // compute position from cell position
1021  for (int k=0; k<dim; k++)
1022  coord[k] = coord[k]*2 + 1; // the doubled grid
1023  if (i%2)
1024  coord[ivar] += 1;
1025  else
1026  coord[ivar] -= 1;
1027 
1028  // encode codim
1030 
1031  // encode level
1032  id = id << yaspgrid_level_bits;
1033  id = id+PersistentIndexType(_g.level());
1034 
1035  // encode coordinates
1036  for (int i=dim-1; i>=0; i--)
1037  {
1038  id = id << yaspgrid_dim_bits;
1039  id = id+PersistentIndexType(coord[i]);
1040  }
1041 
1042  return id;
1043  }
1044 
1045  // map to old numbering
1046  static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1047  i = edge[i];
1048 
1049  if (cc==dim-1) // edges, exist only for dim>2
1050  {
1051  // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1052 
1053  // number of entities per direction
1054  int m=1<<(dim-1);
1055 
1056  // ifix is the direction that is fixed
1057  int ifix=(dim-1)-(i/m);
1058 
1059  // compute position from cell position
1060  int bit=1;
1061  for (int k=0; k<dim; k++)
1062  {
1063  coord[k] = coord[k]*2+1; // cell position in doubled grid
1064  if (k==ifix) continue;
1065  if ((i%m)&bit) coord[k] += 1; else coord[k] -= 1;
1066  bit *= 2;
1067  }
1068 
1069  // encode codim
1070  PersistentIndexType id(dim-1);
1071 
1072  // encode level
1073  id = id << yaspgrid_level_bits;
1074  id = id+PersistentIndexType(_g.level());
1075 
1076  // encode coordinates
1077  for (int i=dim-1; i>=0; i--)
1078  {
1079  id = id << yaspgrid_dim_bits;
1080  id = id+PersistentIndexType(coord[i]);
1081  }
1082 
1083  return id;
1084  }
1085 
1086  DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1087  }
1088 
1090  int subCompressedIndex (int i, int cc) const
1091  {
1092  if (cc==0)
1093  return compressedIndex();
1094 
1095  // get cell position relative to origin of local cell grid
1096  iTupel coord;
1097  for (int k=0; k<dim; ++k)
1098  coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
1099 
1100  if (cc==dim) // vertices
1101  {
1102  // transform cell coordinate to corner coordinate
1103  for (int k=0; k<dim; k++)
1104  if (i&(1<<k)) (coord[k])++;
1105 
1106  // do lexicographic numbering
1107  int index = coord[dim-1];
1108  for (int k=dim-2; k>=0; --k)
1109  index = (index*(_g.cell_overlap().size(k)+1))+coord[k];
1110  return index;
1111  }
1112 
1113  if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1114  {
1115  // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
1116 
1117  // ivar is the direction that varies
1118  int ivar=i/2;
1119 
1120  // compute position from cell position
1121  if (i%2) coord[ivar] += 1;
1122 
1123  // do lexicographic numbering
1124  int index = coord[dim-1];
1125  for (int k=dim-2; k>=0; --k)
1126  if (k==ivar)
1127  index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1128  else
1129  index = (index*(_g.cell_overlap().size(k)))+coord[k];
1130 
1131  // add size of all subsets for smaller directions
1132  for (int j=0; j<ivar; j++)
1133  {
1134  int n=_g.cell_overlap().size(j)+1;
1135  for (int l=0; l<dim; l++)
1136  if (l!=j) n *= _g.cell_overlap().size(l);
1137  index += n;
1138  }
1139 
1140  return index;
1141  }
1142 
1143  // map to old numbering
1144  static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1145  i = edge[i];
1146 
1147  if (cc==dim-1) // edges, exist only for dim>2
1148  {
1149  // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1150 
1151  // number of entities per direction
1152  int m=1<<(dim-1);
1153 
1154  // ifix is the direction that is fixed
1155  int ifix=(dim-1)-(i/m);
1156 
1157  // compute position from cell position
1158  int bit=1;
1159  for (int k=0; k<dim; k++)
1160  {
1161  if (k==ifix) continue;
1162  if ((i%m)&bit) coord[k] += 1;
1163  bit *= 2;
1164  }
1165 
1166  // do lexicographic numbering
1167  int index = coord[dim-1];
1168  for (int k=dim-2; k>=0; --k)
1169  if (k!=ifix)
1170  index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1171  else
1172  index = (index*(_g.cell_overlap().size(k)))+coord[k];
1173 
1174  // add size of all subsets for smaller directions
1175  for (int j=dim-1; j>ifix; j--)
1176  {
1177  int n=_g.cell_overlap().size(j);
1178  for (int l=0; l<dim; l++)
1179  if (l!=j) n *= _g.cell_overlap().size(l)+1;
1180  index += n;
1181  }
1182 
1183  return index;
1184  }
1185 
1186  DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1187  }
1188 
1190  int subCompressedLeafIndex (int i, int cc) const
1191  {
1192  if (cc==0)
1193  return compressedIndex();
1194 
1195  // get cell position relative to origin of local cell grid
1196  iTupel coord;
1197  for (int k=0; k<dim; ++k)
1198  coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
1199 
1200  if (cc==dim) // vertices
1201  {
1202  // transform cell coordinate to corner coordinate
1203  for (int k=0; k<dim; k++)
1204  if (i&(1<<k)) (coord[k])++;
1205 
1206  // move coordinates up to maxlevel
1207  for (int k=0; k<dim; k++)
1208  coord[k] = coord[k]<<(_g.mg()->maxlevel()-_g.level());
1209 
1210  // do lexicographic numbering
1211  int index = coord[dim-1];
1212  for (int k=dim-2; k>=0; --k)
1213  index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
1214  return index;
1215  }
1216 
1217  if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1218  {
1219  // Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
1220 
1221  // ivar is the direction that varies
1222  int ivar=i/2;
1223 
1224  // compute position from cell position
1225  if (i%2) coord[ivar] += 1;
1226 
1227  // do lexicographic numbering
1228  int index = coord[dim-1];
1229  for (int k=dim-2; k>=0; --k)
1230  if (k==ivar)
1231  index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1232  else
1233  index = (index*(_g.cell_overlap().size(k)))+coord[k];
1234 
1235  // add size of all subsets for smaller directions
1236  for (int j=0; j<ivar; j++)
1237  {
1238  int n=_g.cell_overlap().size(j)+1;
1239  for (int l=0; l<dim; l++)
1240  if (l!=j) n *= _g.cell_overlap().size(l);
1241  index += n;
1242  }
1243 
1244  return index;
1245  }
1246 
1247  // map to old numbering
1248  static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1249  i = edge[i];
1250 
1251  if (cc==dim-1) // edges, exist only for dim>2
1252  {
1253  // Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1254 
1255  // number of entities per direction
1256  int m=1<<(dim-1);
1257 
1258  // ifix is the direction that is fixed
1259  int ifix=(dim-1)-(i/m);
1260 
1261  // compute position from cell position
1262  int bit=1;
1263  for (int k=0; k<dim; k++)
1264  {
1265  if (k==ifix) continue;
1266  if ((i%m)&bit) coord[k] += 1;
1267  bit *= 2;
1268  }
1269 
1270  // do lexicographic numbering
1271  int index = coord[dim-1];
1272  for (int k=dim-2; k>=0; --k)
1273  if (k!=ifix)
1274  index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1275  else
1276  index = (index*(_g.cell_overlap().size(k)))+coord[k];
1277 
1278  // add size of all subsets for smaller directions
1279  for (int j=dim-1; j>ifix; j--)
1280  {
1281  int n=_g.cell_overlap().size(j);
1282  for (int l=0; l<dim; l++)
1283  if (l!=j) n *= _g.cell_overlap().size(l)+1;
1284  index += n;
1285  }
1286 
1287  return index;
1288  }
1289 
1290  DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1291  }
1292 
1293  const GridImp * _yg; // access to YaspGrid
1294  const TSI& _it; // position in the grid level
1295  const YGLI& _g; // access to grid level
1296 };
1297 
1298 
1299 // specialization for codim=dim (vertex)
1300 template<int dim, class GridImp>
1301 class YaspEntity<dim,dim,GridImp>
1302 : public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
1303 {
1304  enum { dimworld = GridImp::dimensionworld };
1305 
1306  typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
1307 
1308 public:
1309  typedef typename GridImp::ctype ctype;
1310 
1311  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1312  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1313 
1314  typedef typename GridImp::template Codim<dim>::Geometry Geometry;
1315 
1316  template <int cd>
1317  struct Codim
1318  {
1319  typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
1320  };
1321 
1322  typedef typename GridImp::template Codim<dim>::EntityPointer EntityPointer;
1323  typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
1324 
1326  typedef typename GridImp::PersistentIndexType PersistentIndexType;
1327 
1329  typedef typename YGrid<dim,ctype>::iTupel iTupel;
1330 
1331  // constructor
1332  YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
1333  : _yg(yg), _it(it), _g(g)
1334  {}
1335 
1337  int level () const {return _g.level();}
1338 
1340  int index () const {return _it.superindex();}
1341 
1343  int globalIndex () const { return _g.cell_global().index(_it.coord()); }
1344 
1348  EntitySeed seed () const {
1349  return EntitySeed(_g.level(), _it.coord());
1350  }
1351 
1353  Geometry geometry () const {
1354  GeometryImpl _geometry(_it.position());
1355  return Geometry( _geometry );
1356  }
1357 
1360  {
1361  if (_g.vertex_interior().inside(_it.coord())) return InteriorEntity;
1362  if (_g.vertex_interiorborder().inside(_it.coord())) return BorderEntity;
1363  if (_g.vertex_overlap().inside(_it.coord())) return OverlapEntity;
1364  if (_g.vertex_overlapfront().inside(_it.coord())) return FrontEntity;
1365  return GhostEntity;
1366  }
1367 
1368 private:
1369  // IndexSets needs access to the private index methods
1370  friend class Dune::YaspLevelIndexSet<GridImp>;
1371  friend class Dune::YaspLeafIndexSet<GridImp>;
1372  friend class Dune::YaspGlobalIdSet<GridImp>;
1373 
1376  {
1377  // get coordinate and size of global grid
1378  const iTupel& size = _g.vertex_global().size();
1379  int coord[dim];
1380 
1381  // correction for periodic boundaries
1382  for (int i=0; i<dim; i++)
1383  {
1384  coord[i] = _it.coord(i);
1385  if (coord[i]<0) coord[i] += size[i];
1386  if (coord[i]>=size[i]) coord[i] -= size[i];
1387  }
1388 
1389  // determine min number of trailing zeroes
1390  int trailing = 1000;
1391  for (int i=0; i<dim; i++)
1392  {
1393  // count trailing zeros
1394  int zeros = 0;
1395  for (int j=0; j<_g.level(); j++)
1396  if (coord[i]&(1<<j))
1397  break;
1398  else
1399  zeros++;
1400  trailing = std::min(trailing,zeros);
1401  }
1402 
1403  // determine the level of this vertex
1404  int level = _g.level()-trailing;
1405 
1406  // encode codim
1407  PersistentIndexType id(dim);
1408 
1409  // encode level
1410  id = id << yaspgrid_level_bits;
1411  id = id+PersistentIndexType(level);
1412 
1413  // encode coordinates
1414  for (int i=dim-1; i>=0; i--)
1415  {
1416  id = id << yaspgrid_dim_bits;
1417  id = id+PersistentIndexType(coord[i]>>trailing);
1418  }
1419 
1420  return id;
1421  }
1422 
1424  int compressedIndex () const { return _it.superindex();}
1425 
1427  int compressedLeafIndex () const
1428  {
1429  if (_g.level()==_g.mg()->maxlevel())
1430  return _it.superindex();
1431 
1432  // not on leaf level, interpolate to finest grid
1433  int coord[dim];
1434  for (int i=0; i<dim; i++) coord[i] = _it.coord(i)-(_g).vertex_overlap().origin(i);
1435 
1436  // move coordinates up to maxlevel (multiply by 2 for each level
1437  for (int k=0; k<dim; k++)
1438  coord[k] = coord[k]*(1<<(_g.mg()->maxlevel()-_g.level()));
1439 
1440  // do lexicographic numbering
1441  int index = coord[dim-1];
1442  for (int k=dim-2; k>=0; --k)
1443  index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
1444  return index;
1445  }
1446 
1447 public:
1448  const TSI& transformingsubiterator() const { return _it; }
1449  const YGLI& gridlevel() const { return _g; }
1450  const GridImp * yaspgrid() const { return _yg; }
1451 protected:
1452  const GridImp * _yg; // access to YaspGrid
1453  const TSI& _it; // position in the grid level
1454  const YGLI& _g; // access to grid level
1455  // temporary object
1456  mutable FieldVector<ctype, dim> loc; // always computed before being returned
1457 };
1458 
1459 //========================================================================
1464 //========================================================================
1465 
1466 template<class GridImp>
1468 {
1469  enum { dim=GridImp::dimension };
1470  enum { dimworld=GridImp::dimensionworld };
1471  typedef typename GridImp::ctype ctype;
1472  YaspIntersection();
1473  YaspIntersection& operator = (const YaspIntersection&);
1474 
1475  typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
1476  typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
1477 
1478 public:
1479  // types used from grids
1480  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1481  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1482  typedef typename GridImp::template Codim<0>::Entity Entity;
1483  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
1484  typedef typename GridImp::template Codim<1>::Geometry Geometry;
1485  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
1488 
1489  // void update() const {
1490  // const_cast<YaspIntersection*>(this)->update();
1491  // }
1492  void update() const {
1493  if (_count == 2*_dir + _face || _count >= 2*dim)
1494  return;
1495 
1496  // cleanup old stuff
1497  _outside.transformingsubiterator().move(_dir,1-2*_face); // move home
1498  _pos_world[_dir] = _inside.transformingsubiterator().position(_dir);
1499 
1500  // update face info
1501  _dir = _count / 2;
1502  _face = _count % 2;
1503 
1504  // move transforming iterator
1505  _outside.transformingsubiterator().move(_dir,-1+2*_face);
1506 
1507  // make up faces
1508  _pos_world[_dir] += (-0.5+_face)*_inside.transformingsubiterator().meshsize(_dir);
1509  }
1510 
1512  void increment() {
1513  _count += (_count < 2*dim);
1514  }
1515 
1517  bool equals (const YaspIntersection& other) const
1518  {
1519  return _inside.equals(other._inside) && _count == other._count;
1520  }
1521 
1526  bool boundary () const
1527  {
1528 #if 1
1529  return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 < _inside.gridlevel().cell_global().min(_count/2)
1530  ||
1531  _inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 > _inside.gridlevel().cell_global().max(_count/2));
1532 #else
1533  update();
1534  // The transforming iterator can be safely moved beyond the boundary.
1535  // So we only have to compare against the cell_global grid
1536  return (_outside.transformingsubiterator().coord(_dir) < _inside.gridlevel().cell_global().min(_dir)
1537  ||
1538  _outside.transformingsubiterator().coord(_dir) > _inside.gridlevel().cell_global().max(_dir));
1539 #endif
1540  }
1541 
1543  bool neighbor () const
1544  {
1545 #if 1
1546  return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 >= _inside.gridlevel().cell_overlap().min(_count/2)
1547  &&
1548  _inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 <= _inside.gridlevel().cell_overlap().max(_count/2));
1549 #else
1550  update();
1551  return (_outside.transformingsubiterator().coord(_dir) >= _inside.gridlevel().cell_overlap().min(_dir)
1552  &&
1553  _outside.transformingsubiterator().coord(_dir) <= _inside.gridlevel().cell_overlap().max(_dir));
1554 #endif
1555  }
1556 
1558  bool conforming () const
1559  {
1560  return true;
1561  }
1562 
1566  {
1567  return _inside;
1568  }
1569 
1573  {
1574  update();
1575  return _outside;
1576  }
1577 
1580  int boundaryId() const
1581  {
1582  if(boundary()) return indexInInside()+1;
1583  return 0;
1584  }
1585 
1589  {
1590  if(! boundary())
1591  DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
1592  update();
1593  // size of local macro grid
1594  const FieldVector<int, dim> & size = _inside.gridlevel().mg()->begin().cell_overlap().size();
1595  const FieldVector<int, dim> & origin = _inside.gridlevel().mg()->begin().cell_overlap().origin();
1596  FieldVector<int, dim> sides;
1597  {
1598  for (int i=0; i<dim; i++)
1599  {
1600  sides[i] =
1601  ((_inside.gridlevel().mg()->begin().cell_overlap().origin(i)
1602  == _inside.gridlevel().mg()->begin().cell_global().origin(i))+
1603  (_inside.gridlevel().mg()->begin().cell_overlap().origin(i) +
1604  _inside.gridlevel().mg()->begin().cell_overlap().size(i)
1605  == _inside.gridlevel().mg()->begin().cell_global().origin(i) +
1606  _inside.gridlevel().mg()->begin().cell_global().size(i)));
1607  }
1608  }
1609  // gobal position of the cell on macro grid
1610  FieldVector<int, dim> pos = _inside.transformingsubiterator().coord();
1611  pos /= (1<<_inside.level());
1612  pos -= origin;
1613  // compute unit-cube-face-sizes
1614  FieldVector<int, dim> fsize;
1615  {
1616  int vol = 1;
1617  for (int k=0; k<dim; k++)
1618  vol *= size[k];
1619  for (int k=0; k<dim; k++)
1620  fsize[k] = vol/size[k];
1621  }
1622  // compute index in the unit-cube-face
1623  int index = 0;
1624  {
1625  int localoffset = 1;
1626  for (int k=dim-1; k>=0; k--)
1627  {
1628  if (k == _dir) continue;
1629  index += (pos[k]) * localoffset;
1630  localoffset *= size[k];
1631  }
1632  }
1633  // add unit-cube-face-offsets
1634  {
1635  for (int k=0; k<_dir; k++)
1636  index += sides[k] * fsize[k];
1637  // add fsize if we are on the right face and there is a left-face-boundary on this processor
1638  index += _face * (sides[_dir]>1) * fsize[_dir];
1639  }
1640 
1641  // int rank = 0;
1642  // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1643  // std::cout << rank << "... size: " << size << " sides: " << sides
1644  // << " fsize: " << fsize
1645  // << " pos: " << pos << " face: " << int(_dir) << "/" << int(_face)
1646  // << " index: " << index << std::endl;
1647 
1648  return index;
1649  }
1650 
1652  FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dim-1>& local) const
1653  {
1654  return _faceInfo[_count].normal;
1655  }
1656 
1658  FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dim-1>& local) const
1659  {
1660  return _faceInfo[_count].normal;
1661  }
1662 
1664  FieldVector<ctype, dimworld> centerUnitOuterNormal () const
1665  {
1666  return _faceInfo[_count].normal;
1667  }
1668 
1672  FieldVector<ctype, dimworld> integrationOuterNormal (const FieldVector<ctype, dim-1>& local) const
1673  {
1674  FieldVector<ctype, dimworld> n = _faceInfo[_count].normal;
1675  n *= geometry().volume();
1676  return n;
1677  }
1678 
1683  {
1684  return LocalGeometry( _faceInfo[_count].geom_inside );
1685  }
1686 
1691  {
1692  return LocalGeometry( _faceInfo[_count].geom_outside );
1693  }
1694 
1698  {
1699  update();
1700  GeometryImpl
1701  _is_global(_pos_world,_inside.transformingsubiterator().meshsize(),_dir);
1702  return Geometry( _is_global );
1703  }
1704 
1707  {
1708  static const GeometryType cube(GeometryType::cube, dim-1);
1709  return cube;
1710  }
1711 
1713  int indexInInside () const
1714  {
1715  return _count;
1716  }
1717 
1719  int indexInOutside () const
1720  {
1721  // flip the last bit
1722  return _count^1;
1723  }
1724 
1726  YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
1727  _inside(myself.yaspgrid(), myself.gridlevel(),
1728  myself.transformingsubiterator()),
1729  _outside(myself.yaspgrid(), myself.gridlevel(),
1730  myself.transformingsubiterator()),
1731  // initialize to first neighbor
1732  _count(0),
1733  _dir(0),
1734  _face(0),
1735  _pos_world(myself.transformingsubiterator().position())
1736  {
1737  if (toend)
1738  {
1739  // initialize end iterator
1740  _count = 2*dim;
1741  return;
1742  }
1743  _count = 0;
1744 
1745  // move transforming iterator
1746  _outside.transformingsubiterator().move(_dir,-1);
1747 
1748  // make up faces
1749  _pos_world[0] -= 0.5*_inside.transformingsubiterator().meshsize(0);
1750  }
1751 
1754  _inside(it._inside),
1755  _outside(it._outside),
1756  _count(it._count),
1757  _dir(it._dir),
1758  _face(it._face),
1759  _pos_world(it._pos_world)
1760  {}
1761 
1763  void assign (const YaspIntersection& it)
1764  {
1765  _inside = it._inside;
1766  _outside = it._outside;
1767  _count = it._count;
1768  _dir = it._dir;
1769  _face = it._face;
1770  _pos_world = it._pos_world;
1771  }
1772 
1773 private:
1774  /* EntityPointers (get automatically updated) */
1775  mutable YaspEntityPointer<0,GridImp> _inside;
1776  mutable YaspEntityPointer<0,GridImp> _outside;
1777  /* current position */
1778  uint8_t _count;
1779  mutable uint8_t _dir;
1780  mutable uint8_t _face;
1781  /* current position */
1782  mutable FieldVector<ctype, dimworld> _pos_world;
1783 
1784  /* static data */
1785  struct faceInfo
1786  {
1787  FieldVector<ctype, dimworld> normal;
1788  LocalGeometryImpl geom_inside;
1789  LocalGeometryImpl geom_outside;
1790  };
1791 
1792  /* static face info */
1793  static const array<faceInfo, 2*GridImp::dimension> _faceInfo;
1794 
1795  static array<faceInfo, 2*dim> initFaceInfo()
1796  {
1797  const FieldVector<typename GridImp::ctype, GridImp::dimension> ext_local(1.0);
1798  array<faceInfo, 2*dim> I;
1799  for (uint8_t i=0; i<dim; i++)
1800  {
1801  // center of face
1802  FieldVector<ctype, dim> a(0.5); a[i] = 0.0;
1803  FieldVector<ctype, dim> b(0.5); b[i] = 1.0;
1804  // normal vectors
1805  I[2*i].normal = 0.0;
1806  I[2*i+1].normal = 0.0;
1807  I[2*i].normal[i] = -1.0;
1808  I[2*i+1].normal[i] = +1.0;
1809  // geometries
1810  I[2*i].geom_inside =
1811  LocalGeometryImpl(a, ext_local, i);
1812  I[2*i].geom_outside =
1813  LocalGeometryImpl(b, ext_local, i);
1814  I[2*i+1].geom_inside =
1815  LocalGeometryImpl(b, ext_local, i);
1816  I[2*i+1].geom_outside =
1817  LocalGeometryImpl(a, ext_local, i);
1818  }
1819  return I;
1820  }
1821 };
1822 
1823 template<class GridImp>
1824 const array<typename YaspIntersection<GridImp>::faceInfo, 2*GridImp::dimension>
1825 YaspIntersection<GridImp>::_faceInfo =
1826  YaspIntersection<GridImp>::initFaceInfo();
1827 
1828 //========================================================================
1833 //========================================================================
1834 
1835 template<class GridImp>
1836 class YaspIntersectionIterator : public MakeableInterfaceObject< Dune::Intersection<const GridImp, Dune::YaspIntersection > >
1837 {
1838  enum { dim=GridImp::dimension };
1839  enum { dimworld=GridImp::dimensionworld };
1840  typedef typename GridImp::ctype ctype;
1842 public:
1843  // types used from grids
1846 
1848  void increment()
1849  {
1850  GridImp::getRealImplementation(*this).increment();
1851  }
1852 
1854  bool equals (const YaspIntersectionIterator& other) const
1855  {
1856  return GridImp::getRealImplementation(*this).equals(
1857  GridImp::getRealImplementation(other));
1858  }
1859 
1861  const Intersection & dereference() const
1862  {
1863  return *this;
1864  }
1865 
1868  MakeableIntersection(YaspIntersection<GridImp>(myself, toend))
1869  {}
1870 
1874  {}
1875 
1878  {
1879  GridImp::getRealImplementation(*this).assign(
1880  GridImp::getRealImplementation(it));
1881  return *this;
1882  }
1883 };
1884 
1885 //========================================================================
1889 //========================================================================
1890 
1891 template<class GridImp>
1893  public YaspEntityPointer<0,GridImp>
1894 {
1895  enum { dim=GridImp::dimension };
1896  enum { dimworld=GridImp::dimensionworld };
1897  typedef typename GridImp::ctype ctype;
1898 public:
1899  // types used from grids
1900  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1901  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1902  typedef typename GridImp::template Codim<0>::Entity Entity;
1904 
1906  typedef typename YGrid<dim,ctype>::iTupel iTupel;
1907 
1909  YaspHierarchicIterator (const GridImp* yg, const YGLI& g, const TSI& it, int maxlevel) :
1910  YaspEntityPointer<0,GridImp>(yg,g,it)
1911  {
1912  // now iterator points to current cell
1913  StackElem se(this->_g);
1914  se.coord = this->_it.coord();
1915  stack.push(se);
1916 
1917  // determine maximum level
1918  _maxlevel = std::min(maxlevel,this->_g.mg()->maxlevel());
1919 
1920  // if maxlevel not reached then push yourself and sons
1921  if (this->_g.level()<_maxlevel)
1922  {
1923  push_sons();
1924  }
1925 
1926  // and make iterator point to first son if stack is not empty
1927  if (!stack.empty())
1928  pop_tos();
1929  }
1930 
1933  YaspEntityPointer<0,GridImp>(it),
1934  _maxlevel(it._maxlevel), stack(it.stack)
1935  {}
1936 
1938  void increment ()
1939  {
1940  // sanity check: do nothing when stack is empty
1941  if (stack.empty()) return;
1942 
1943  // if maxlevel not reached then push sons
1944  if (this->_g.level()<_maxlevel)
1945  push_sons();
1946 
1947  // in any case pop one element
1948  pop_tos();
1949  }
1950 
1951  void print (std::ostream& s) const
1952  {
1953  s << "HIER: " << "level=" << this->_g.level()
1954  << " position=" << this->_it.coord()
1955  << " superindex=" << this->_it.superindex()
1956  << " maxlevel=" << this->_maxlevel
1957  << " stacksize=" << stack.size()
1958  << std::endl;
1959  }
1960 
1961 private:
1962  int _maxlevel;
1963 
1964  struct StackElem {
1965  YGLI g; // grid level of the element
1966  iTupel coord; // and the coordinates
1967  StackElem(YGLI gg) : g(gg) {}
1968  };
1969  std::stack<StackElem> stack;
1970 
1971  // push sons of current element on the stack
1972  void push_sons ()
1973  {
1974  // yes, process all 1<<dim sons
1975  StackElem se(this->_g.finer());
1976  for (int i=0; i<(1<<dim); i++)
1977  {
1978  for (int k=0; k<dim; k++)
1979  if (i&(1<<k))
1980  se.coord[k] = this->_it.coord(k)*2+1;
1981  else
1982  se.coord[k] = this->_it.coord(k)*2;
1983  stack.push(se);
1984  }
1985  }
1986 
1987  // make TOS the current element
1988  void pop_tos ()
1989  {
1990  StackElem se = stack.top();
1991  stack.pop();
1992  this->_g = se.g;
1993  this->_it.reinit(this->_g.cell_overlap(),se.coord);
1994  }
1995 };
1996 
1997 //========================================================================
2001 //========================================================================
2002 template<int codim, class GridImp>
2004 {
2006  enum { dim=GridImp::dimension };
2007 
2008 public:
2010  enum { codimension = codim };
2011 
2013  YaspEntitySeed (int level, FieldVector<int, dim> coord)
2014  : _l(level), _c(coord)
2015  {}
2016 
2019  : _l(rhs._l), _c(rhs._c)
2020  {}
2021 
2022  int level () const { return _l; }
2023  const FieldVector<int, dim> & coord() const { return _c; }
2024 
2025 protected:
2026  int _l; // grid level
2027  FieldVector<int, dim> _c; // coord in the global grid
2028 };
2029 
2030 //========================================================================
2040 //========================================================================
2041 template<int codim, class GridImp>
2043 {
2045  enum { dim=GridImp::dimension };
2047  enum { dimworld=GridImp::dimensionworld };
2048  typedef typename GridImp::ctype ctype;
2049 public:
2050  typedef typename GridImp::template Codim<codim>::Entity Entity;
2051  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2052  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2055 protected:
2057 
2058 public:
2060  enum { codimension = codim };
2061 
2063  YaspEntityPointer (const GridImp * yg, const YGLI & g, const TSI & it)
2064  : _g(g), _it(it), _entity(yg, _g,_it)
2065  {
2066  if (codim>0 && codim<dim)
2067  {
2068  DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2069  }
2070  }
2071 
2074  : _g(entity.gridlevel()),
2075  _it(entity.transformingsubiterator()),
2076  _entity(entity.yaspgrid(),_g,_it)
2077  {
2078  if (codim>0 && codim<dim)
2079  {
2080  DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2081  }
2082  }
2083 
2086  : _g(rhs._g), _it(rhs._it), _entity(rhs._entity.yaspgrid(),_g,_it)
2087  {
2088  if (codim>0 && codim<dim)
2089  {
2090  DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2091  }
2092  }
2093 
2095  bool equals (const YaspEntityPointer& rhs) const
2096  {
2097  return (_it==rhs._it && _g == rhs._g);
2098  }
2099 
2102  {
2103  return _entity;
2104  }
2105 
2107  int level () const {return _g.level();}
2108 
2109  const YaspEntityPointer&
2111  {
2112  _g = rhs._g;
2113  _it = rhs._it;
2114  /* _entity = i._entity
2115  * is done implicitely, as the entity is completely
2116  * defined via the interator it belongs to
2117  */
2118  return *this;
2119  }
2120 
2122  {
2123  return _it;
2124  }
2125 
2126  const YGLI& gridlevel () const
2127  {
2128  return _g;
2129  }
2130 
2132  {
2133  return _it;
2134  }
2135 
2137  {
2138  return _g;
2139  }
2140 
2141 protected:
2142  YGLI _g; // access to grid level
2143  TSI _it; // position in the grid level
2145 };
2146 
2147 //========================================================================
2155 //========================================================================
2156 
2157 template<int codim, PartitionIteratorType pitype, class GridImp>
2159  public YaspEntityPointer<codim,GridImp>
2160 {
2162  enum { dim=GridImp::dimension };
2164  enum { dimworld=GridImp::dimensionworld };
2165  typedef typename GridImp::ctype ctype;
2166 public:
2167  typedef typename GridImp::template Codim<codim>::Entity Entity;
2168  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2169  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2171 
2173  YaspLevelIterator (const GridImp * yg, const YGLI & g, const TSI & it) :
2174  YaspEntityPointer<codim,GridImp>(yg,g,it) {}
2175 
2178  YaspEntityPointer<codim,GridImp>(i) {}
2179 
2181  void increment()
2182  {
2183  ++(this->_it);
2184  }
2185 };
2186 
2187 
2188 //========================================================================
2193 //========================================================================
2194 
2195 template<class GridImp>
2197  : public IndexSet< GridImp, YaspLevelIndexSet< GridImp >, unsigned int >
2198 {
2201 
2202 public:
2203  typedef typename Base::IndexType IndexType;
2204 
2205  using Base::subIndex;
2206 
2208  YaspLevelIndexSet ( const GridImp &g, int l )
2209  : grid( g ),
2210  level( l )
2211  {
2212  // contains a single element type;
2213  for (int codim=0; codim<=GridImp::dimension; codim++)
2214  mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
2215  }
2216 
2218  template<int cc>
2219  IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const
2220  {
2221  assert( cc == 0 || cc == GridImp::dimension );
2222  return grid.getRealImplementation(e).compressedIndex();
2223  }
2224 
2226  template< int cc >
2227  IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
2228  int i, unsigned int codim ) const
2229  {
2230  assert( cc == 0 || cc == GridImp::dimension );
2231  if( cc == GridImp::dimension )
2232  return grid.getRealImplementation(e).compressedIndex();
2233  else
2234  return grid.getRealImplementation(e).subCompressedIndex(i,codim);
2235  }
2236 
2238  int size (GeometryType type) const
2239  {
2240  return grid.size( level, type );
2241  }
2242 
2244  int size (int codim) const
2245  {
2246  return grid.size( level, codim );
2247  }
2248 
2250  template<class EntityType>
2251  bool contains (const EntityType& e) const
2252  {
2253  return e.level() == level;
2254  }
2255 
2257  const std::vector<GeometryType>& geomTypes (int codim) const
2258  {
2259  return mytypes[codim];
2260  }
2261 
2262 private:
2263  const GridImp& grid;
2264  int level;
2265  std::vector<GeometryType> mytypes[GridImp::dimension+1];
2266 };
2267 
2268 
2269 // Leaf Index Set
2270 
2271 template<class GridImp>
2273  : public IndexSet< GridImp, YaspLeafIndexSet< GridImp >, unsigned int >
2274 {
2277 
2278 public:
2279  typedef typename Base::IndexType IndexType;
2280 
2281  using Base::subIndex;
2282 
2284  explicit YaspLeafIndexSet ( const GridImp &g )
2285  : grid( g )
2286  {
2287  // contains a single element type;
2288  for (int codim=0; codim<=GridImp::dimension; codim++)
2289  mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
2290  }
2291 
2293  template<int cc>
2294  IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const
2295  {
2296  assert( cc == 0 || cc == GridImp::dimension );
2297  return grid.getRealImplementation(e).compressedIndex();
2298  }
2299 
2301  template< int cc >
2302  IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
2303  int i, unsigned int codim ) const
2304  {
2305  assert( cc == 0 || cc == GridImp::dimension );
2306  if( cc == GridImp::dimension )
2307  return grid.getRealImplementation(e).compressedIndex();
2308  else
2309  return grid.getRealImplementation(e).subCompressedIndex(i,codim);
2310  }
2311 
2313  int size (GeometryType type) const
2314  {
2315  return grid.size( grid.maxLevel(), type );
2316  }
2317 
2319  int size (int codim) const
2320  {
2321  return grid.size( grid.maxLevel(), codim );
2322  }
2323 
2325  template<class EntityType>
2326  bool contains (const EntityType& e) const
2327  {
2328  //return e.isLeaf();
2329  return (e.level() == grid.maxLevel());
2330  }
2331 
2333  const std::vector<GeometryType>& geomTypes (int codim) const
2334  {
2335  return mytypes[codim];
2336  }
2337 
2338 private:
2339  const GridImp& grid;
2340  enum { ncodim = remove_const<GridImp>::type::dimension+1 };
2341  std::vector<GeometryType> mytypes[ncodim];
2342 };
2343 
2344 
2345 
2346 
2347 //========================================================================
2352 //========================================================================
2353 
2354 template<class GridImp>
2355 class YaspGlobalIdSet : public IdSet<GridImp,YaspGlobalIdSet<GridImp>,
2356  typename remove_const<GridImp>::type::PersistentIndexType >
2357 /*
2358  We used the remove_const to extract the Type from the mutable class,
2359  because the const class is not instantiated yet.
2360 */
2361 {
2363 
2364 public:
2366  typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;
2367 
2369 
2371  explicit YaspGlobalIdSet ( const GridImp &g )
2372  : grid( g )
2373  {}
2374 
2376  /*
2377  We use the remove_const to extract the Type from the mutable class,
2378  because the const class is not instantiated yet.
2379  */
2380  template<int cd>
2381  IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
2382  {
2383  return grid.getRealImplementation(e).persistentIndex();
2384  }
2385 
2387  /*
2388  We use the remove_const to extract the Type from the mutable class,
2389  because the const class is not instantiated yet.
2390  */
2391  IdType subId (const typename remove_const<GridImp>::type::Traits::template Codim< 0 >::Entity &e,
2392  int i, unsigned int codim ) const
2393  {
2394  return grid.getRealImplementation(e).subPersistentIndex(i,codim);
2395  }
2396 
2397 private:
2398  const GridImp& grid;
2399 };
2400 
2401 
2402 template<int dim, int dimworld>
2404 {
2405 #if HAVE_MPI
2406  typedef CollectiveCommunication<MPI_Comm> CCType;
2407 #else
2408  typedef CollectiveCommunication<Dune::YaspGrid<dim> > CCType;
2409 #endif
2410 
2414  YaspIntersection, // leaf intersection
2415  YaspIntersection, // level intersection
2416  YaspIntersectionIterator, // leaf intersection iter
2417  YaspIntersectionIterator, // level intersection iter
2423  bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
2425  bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
2426  CCType,
2430 };
2431 
2432 template<int dim, int codim>
2434  template<class G, class DataHandle>
2435  static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2436  {
2437  if (data.contains(dim,codim))
2438  {
2439  DUNE_THROW(GridError, "interface communication not implemented");
2440  }
2441  YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
2442  }
2443 };
2444 
2445 template<int dim>
2446 struct YaspCommunicateMeta<dim,dim> {
2447  template<class G, class DataHandle>
2448  static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2449  {
2450  if (data.contains(dim,dim))
2451  g.template communicateCodim<DataHandle,dim>(data,iftype,dir,level);
2452  YaspCommunicateMeta<dim,dim-1>::comm(g,data,iftype,dir,level);
2453  }
2454 };
2455 
2456 template<int dim>
2457 struct YaspCommunicateMeta<dim,0> {
2458  template<class G, class DataHandle>
2459  static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2460  {
2461  if (data.contains(dim,0))
2462  g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
2463  }
2464 };
2465 
2466 
2467 //************************************************************************
2483 template<int dim>
2484 class YaspGrid :
2485  public GridDefaultImplementation<dim,dim,yaspgrid_ctype,YaspGridFamily<dim,dim> >,
2486  public MultiYGrid<dim,yaspgrid_ctype>
2487 {
2488  typedef const YaspGrid<dim> GridImp;
2489 
2490  void init()
2491  {
2492  setsizes();
2493  indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,0) );
2494  theleafindexset.push_back( new YaspLeafIndexSet<const YaspGrid<dim> >(*this) );
2495  theglobalidset.push_back( new YaspGlobalIdSet<const YaspGrid<dim> >(*this) );
2496  boundarysegmentssize();
2497  }
2498 
2499  void boundarysegmentssize()
2500  {
2501  // sizes of local macro grid
2502  const FieldVector<int, dim> & size = MultiYGrid<dim,ctype>::begin().cell_overlap().size();
2503  FieldVector<int, dim> sides;
2504  {
2505  for (int i=0; i<dim; i++)
2506  {
2507  sides[i] =
2508  ((MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i)
2509  == MultiYGrid<dim,ctype>::begin().cell_global().origin(i))+
2510  (MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i) +
2511  MultiYGrid<dim,ctype>::begin().cell_overlap().size(i)
2512  == MultiYGrid<dim,ctype>::begin().cell_global().origin(i) +
2513  MultiYGrid<dim,ctype>::begin().cell_global().size(i)));
2514  }
2515  }
2516  nBSegments = 0;
2517  for (int k=0; k<dim; k++)
2518  {
2519  int offset = 1;
2520  for (int l=0; l<dim; l++)
2521  {
2522  if (l==k) continue;
2523  offset *= size[l];
2524  }
2525  nBSegments += sides[k]*offset;
2526  }
2527  }
2528 
2529 public:
2530 
2531  using MultiYGrid<dim,yaspgrid_ctype>::defaultLoadbalancer;
2532 
2535 
2536  // define the persistent index type
2537  typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits> PersistentIndexType;
2538 
2541  // the Traits
2543 
2544  // need for friend declarations in entity
2548 
2550  enum { MAXL=64 };
2551 
2553  typedef MultiYGrid<dim,ctype> YMG;
2554  typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2555  typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2556  typedef typename MultiYGrid<dim,ctype>::Intersection IS;
2557  typedef typename std::deque<IS>::const_iterator ISIT;
2558 
2567  YaspGrid (Dune::MPIHelper::MPICommunicator comm,
2568  Dune::FieldVector<ctype, dim> L,
2569  Dune::FieldVector<int, dim> s,
2570  Dune::FieldVector<bool, dim> periodic, int overlap,
2571  const YLoadBalance<dim>* lb = defaultLoadbalancer())
2572 #if HAVE_MPI
2573  : YMG(comm,L,s,periodic,overlap,lb), ccobj(comm),
2574  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2575 #else
2576  : YMG(L,s,periodic,overlap,lb),
2577  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2578 #endif
2579  {
2580  init();
2581  }
2582 
2583 
2596  YaspGrid (Dune::FieldVector<ctype, dim> L,
2597  Dune::FieldVector<int, dim> s,
2598  Dune::FieldVector<bool, dim> periodic, int overlap,
2599  const YLoadBalance<dim>* lb = YMG::defaultLoadbalancer())
2600 #if HAVE_MPI
2601  : YMG(MPI_COMM_SELF,L,s,periodic,overlap,lb), ccobj(MPI_COMM_SELF),
2602  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2603 #else
2604  : YMG(L,s,periodic,overlap,lb),
2605  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2606 #endif
2607  {
2608  init();
2609  }
2610 
2612  {
2613  deallocatePointers(indexsets);
2614  deallocatePointers(theleafindexset);
2615  deallocatePointers(theglobalidset);
2616  }
2617 
2618  private:
2619  // do not copy this class
2620  YaspGrid(const YaspGrid&);
2621 
2622  public:
2623 
2627  int maxLevel() const {return MultiYGrid<dim,ctype>::maxlevel();} // delegate
2628 
2630  void globalRefine (int refCount)
2631  {
2632  if (refCount < -maxLevel())
2633  DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
2634  "Coarsening " << -refCount << " levels requested!");
2635  for (int k=refCount; k<0; k++)
2636  {
2637  MultiYGrid<dim,ctype>::coarsen();
2638  setsizes();
2639  indexsets.pop_back();
2640  }
2641  for (int k=0; k<refCount; k++)
2642  {
2643  MultiYGrid<dim,ctype>::refine(keep_ovlp);
2644  setsizes();
2645  indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,maxLevel()) );
2646  }
2647  }
2648 
2653  void refineOptions (bool keepPhysicalOverlap)
2654  {
2655  keep_ovlp = keepPhysicalOverlap;
2656  }
2657 
2669  bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
2670  {
2671  assert(adaptActive == false);
2672  if (e.level() != maxLevel()) return false;
2673  adaptRefCount = std::max(adaptRefCount, refCount);
2674  return true;
2675  }
2676 
2683  int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
2684  {
2685  return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
2686  }
2687 
2689  bool adapt ()
2690  {
2691  globalRefine(adaptRefCount);
2692  return (adaptRefCount > 0);
2693  }
2694 
2696  bool preAdapt ()
2697  {
2698  adaptActive = true;
2699  adaptRefCount = comm().max(adaptRefCount);
2700  return (adaptRefCount < 0);
2701  }
2702 
2704  void postAdapt()
2705  {
2706  adaptActive = false;
2707  adaptRefCount = 0;
2708  }
2709 
2711  template<int cd, PartitionIteratorType pitype>
2712  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
2713  {
2714  return levelbegin<cd,pitype>(level);
2715  }
2716 
2718  template<int cd, PartitionIteratorType pitype>
2719  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
2720  {
2721  return levelend<cd,pitype>(level);
2722  }
2723 
2725  template<int cd>
2726  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
2727  {
2728  return levelbegin<cd,All_Partition>(level);
2729  }
2730 
2732  template<int cd>
2733  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
2734  {
2735  return levelend<cd,All_Partition>(level);
2736  }
2737 
2739  template<int cd, PartitionIteratorType pitype>
2741  {
2742  return levelbegin<cd,pitype>(maxLevel());
2743  }
2744 
2746  template<int cd, PartitionIteratorType pitype>
2748  {
2749  return levelend<cd,pitype>(maxLevel());
2750  }
2751 
2753  template<int cd>
2755  {
2756  return levelbegin<cd,All_Partition>(maxLevel());
2757  }
2758 
2760  template<int cd>
2762  {
2763  return levelend<cd,All_Partition>(maxLevel());
2764  }
2765 
2766  // \brief obtain EntityPointer from EntitySeed. */
2767  template <typename Seed>
2768  typename Traits::template Codim<Seed::codimension>::EntityPointer
2769  entityPointer(const Seed& seed) const
2770  {
2771  static const int codim = Seed::codimension;
2772  YGLI g = MultiYGrid<dim,ctype>::begin(seed.level());
2773  switch (codim)
2774  {
2775  case 0:
2776  return YaspEntityPointer<codim,GridImp>(this,g,
2777  TSI(g.cell_overlap(), seed.coord()));
2778  case dim:
2779  return YaspEntityPointer<codim,GridImp>(this,g,
2780  TSI(g.vertex_overlap(), seed.coord()));
2781  default:
2782  DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2783  }
2784  }
2785 
2787  int overlapSize (int level, int codim) const
2788  {
2789  YGLI g = MultiYGrid<dim,ctype>::begin(level);
2790  return g.overlap();
2791  }
2792 
2794  int overlapSize (int codim) const
2795  {
2796  YGLI g = MultiYGrid<dim,ctype>::begin(maxLevel());
2797  return g.overlap();
2798  }
2799 
2801  int ghostSize (int level, int codim) const
2802  {
2803  return 0;
2804  }
2805 
2807  int ghostSize (int codim) const
2808  {
2809  return 0;
2810  }
2811 
2813  int size (int level, int codim) const
2814  {
2815  return sizes[level][codim];
2816  }
2817 
2819  int size (int codim) const
2820  {
2821  return sizes[maxLevel()][codim];
2822  }
2823 
2825  int size (int level, GeometryType type) const
2826  {
2827  return (type.isCube()) ? sizes[level][dim-type.dim()] : 0;
2828  }
2829 
2831  int size (GeometryType type) const
2832  {
2833  return size(maxLevel(),type);
2834  }
2835 
2837  size_t numBoundarySegments () const
2838  {
2839  return nBSegments;
2840  }
2841 
2846  template<class DataHandleImp, class DataType>
2848  {
2849  YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
2850  }
2851 
2856  template<class DataHandleImp, class DataType>
2858  {
2859  YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
2860  }
2861 
2866  template<class DataHandle, int codim>
2867  void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
2868  {
2869  // check input
2870  if (!data.contains(dim,codim)) return; // should have been checked outside
2871 
2872  // data types
2873  typedef typename DataHandle::DataType DataType;
2874 
2875  // access to grid level
2876  YGLI g = MultiYGrid<dim,ctype>::begin(level);
2877 
2878  // find send/recv lists or throw error
2879  const std::deque<IS>* sendlist=0;
2880  const std::deque<IS>* recvlist=0;
2881  if (codim==0) // the elements
2882  {
2884  return; // there is nothing to do in this case
2885  if (iftype==InteriorBorder_All_Interface)
2886  {
2887  sendlist = &g.send_cell_interior_overlap();
2888  recvlist = &g.recv_cell_overlap_interior();
2889  }
2891  {
2892  sendlist = &g.send_cell_overlap_overlap();
2893  recvlist = &g.recv_cell_overlap_overlap();
2894  }
2895  }
2896  if (codim==dim) // the vertices
2897  {
2899  {
2900  sendlist = &g.send_vertex_interiorborder_interiorborder();
2901  recvlist = &g.recv_vertex_interiorborder_interiorborder();
2902  }
2903 
2904  if (iftype==InteriorBorder_All_Interface)
2905  {
2906  sendlist = &g.send_vertex_interiorborder_overlapfront();
2907  recvlist = &g.recv_vertex_overlapfront_interiorborder();
2908  }
2910  {
2911  sendlist = &g.send_vertex_overlap_overlapfront();
2912  recvlist = &g.recv_vertex_overlapfront_overlap();
2913  }
2914  if (iftype==All_All_Interface)
2915  {
2916  sendlist = &g.send_vertex_overlapfront_overlapfront();
2917  recvlist = &g.recv_vertex_overlapfront_overlapfront();
2918  }
2919  }
2920 
2921  // change communication direction?
2922  if (dir==BackwardCommunication)
2923  std::swap(sendlist,recvlist);
2924 
2925  int cnt;
2926 
2927  // Size computation (requires communication if variable size)
2928  std::vector<int> send_size(sendlist->size(),-1); // map rank to total number of objects (of type DataType) to be sent
2929  std::vector<int> recv_size(recvlist->size(),-1); // map rank to total number of objects (of type DataType) to be recvd
2930  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
2931  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
2932  if (data.fixedsize(dim,codim))
2933  {
2934  // fixed size: just take a dummy entity, size can be computed without communication
2935  cnt=0;
2936  for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
2937  {
2939  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2940  send_size[cnt] = is->grid.totalsize() * data.size(*it);
2941  cnt++;
2942  }
2943  cnt=0;
2944  for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
2945  {
2947  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2948  recv_size[cnt] = is->grid.totalsize() * data.size(*it);
2949  cnt++;
2950  }
2951  }
2952  else
2953  {
2954  // variable size case: sender side determines the size
2955  cnt=0;
2956  for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
2957  {
2958  // allocate send buffer for sizes per entitiy
2959  size_t *buf = new size_t[is->grid.totalsize()];
2960  send_sizes[cnt] = buf;
2961 
2962  // loop over entities and ask for size
2963  int i=0; size_t n=0;
2965  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2967  tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
2968  for ( ; it!=tsubend; ++it)
2969  {
2970  buf[i] = data.size(*it);
2971  n += buf[i];
2972  i++;
2973  }
2974 
2975  // now we know the size for this rank
2976  send_size[cnt] = n;
2977 
2978  // hand over send request to torus class
2979  MultiYGrid<dim,ctype>::torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
2980  cnt++;
2981  }
2982 
2983  // allocate recv buffers for sizes and store receive request
2984  cnt=0;
2985  for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
2986  {
2987  // allocate recv buffer
2988  size_t *buf = new size_t[is->grid.totalsize()];
2989  recv_sizes[cnt] = buf;
2990 
2991  // hand over recv request to torus class
2992  MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
2993  cnt++;
2994  }
2995 
2996  // exchange all size buffers now
2997  MultiYGrid<dim,ctype>::torus().exchange();
2998 
2999  // release send size buffers
3000  cnt=0;
3001  for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3002  {
3003  delete[] send_sizes[cnt];
3004  send_sizes[cnt] = 0;
3005  cnt++;
3006  }
3007 
3008  // process receive size buffers
3009  cnt=0;
3010  for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3011  {
3012  // get recv buffer
3013  size_t *buf = recv_sizes[cnt];
3014 
3015  // compute total size
3016  size_t n=0;
3017  for (int i=0; i<is->grid.totalsize(); ++i)
3018  n += buf[i];
3019 
3020  // ... and store it
3021  recv_size[cnt] = n;
3022  ++cnt;
3023  }
3024  }
3025 
3026 
3027  // allocate & fill the send buffers & store send request
3028  std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
3029  cnt=0;
3030  for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3031  {
3032 // std::cout << "[" << this->comm().rank() << "] "
3033 // << " send " << " dest=" << is->rank
3034 // << " size=" << send_size[cnt]
3035 // << std::endl;
3036 
3037  // allocate send buffer
3038  DataType *buf = new DataType[send_size[cnt]];
3039 
3040  // remember send buffer
3041  sends[cnt] = buf;
3042 
3043  // make a message buffer
3044  MessageBuffer<DataType> mb(buf);
3045 
3046  // fill send buffer; iterate over cells in intersection
3048  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3050  tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3051  for ( ; it!=tsubend; ++it)
3052  data.gather(mb,*it);
3053 
3054  // hand over send request to torus class
3055  MultiYGrid<dim,ctype>::torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
3056  cnt++;
3057  }
3058 
3059  // allocate recv buffers and store receive request
3060  std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
3061  cnt=0;
3062  for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3063  {
3064 // std::cout << "[" << this->comm().rank() << "] "
3065 // << " recv " << " src=" << is->rank
3066 // << " size=" << recv_size[cnt]
3067 // << std::endl;
3068 
3069  // allocate recv buffer
3070  DataType *buf = new DataType[recv_size[cnt]];
3071 
3072  // remember recv buffer
3073  recvs[cnt] = buf;
3074 
3075  // hand over recv request to torus class
3076  MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
3077  cnt++;
3078  }
3079 
3080  // exchange all buffers now
3081  MultiYGrid<dim,ctype>::torus().exchange();
3082 
3083  // release send buffers
3084  cnt=0;
3085  for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3086  {
3087  delete[] sends[cnt];
3088  sends[cnt] = 0;
3089  cnt++;
3090  }
3091 
3092  // process receive buffers and delete them
3093  cnt=0;
3094  for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3095  {
3096  // get recv buffer
3097  DataType *buf = recvs[cnt];
3098 
3099  // make a message buffer
3100  MessageBuffer<DataType> mb(buf);
3101 
3102  // copy data from receive buffer; iterate over cells in intersection
3103  if (data.fixedsize(dim,codim))
3104  {
3106  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3107  size_t n=data.size(*it);
3109  tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3110  for ( ; it!=tsubend; ++it)
3111  data.scatter(mb,*it,n);
3112  }
3113  else
3114  {
3115  int i=0;
3116  size_t *sbuf = recv_sizes[cnt];
3118  it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3120  tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3121  for ( ; it!=tsubend; ++it)
3122  data.scatter(mb,*it,sbuf[i++]);
3123  delete[] sbuf;
3124  }
3125 
3126  // delete buffer
3127  delete[] buf; // hier krachts !
3128  cnt++;
3129  }
3130  }
3131 
3132  // The new index sets from DDM 11.07.2005
3133  const typename Traits::GlobalIdSet& globalIdSet() const
3134  {
3135  return *(theglobalidset[0]);
3136  }
3137 
3138  const typename Traits::LocalIdSet& localIdSet() const
3139  {
3140  return *(theglobalidset[0]);
3141  }
3142 
3143  const typename Traits::LevelIndexSet& levelIndexSet(int level) const
3144  {
3145  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3146  return *(indexsets[level]);
3147  }
3148 
3149  const typename Traits::LeafIndexSet& leafIndexSet() const
3150  {
3151  return *(theleafindexset[0]);
3152  }
3153 
3154 #if HAVE_MPI
3155 
3157  const CollectiveCommunication<MPI_Comm>& comm () const
3158  {
3159  return ccobj;
3160  }
3161 #else
3162 
3164  const CollectiveCommunication<YaspGrid>& comm () const
3165  {
3166  return ccobj;
3167  }
3168 #endif
3169 
3170  YaspIntersectionIterator<const YaspGrid<dim> >&
3172  {
3173  return this->getRealImplementation(it);
3174  }
3175 
3178  {
3179  return this->getRealImplementation(it);
3180  }
3181 
3182 private:
3183 
3184 #if HAVE_MPI
3185  CollectiveCommunication<MPI_Comm> ccobj;
3186 #else
3187  CollectiveCommunication<YaspGrid> ccobj;
3188 #endif
3189 
3190  std::vector<YaspLevelIndexSet<const YaspGrid<dim> >*> indexsets;
3191  std::vector<YaspLeafIndexSet<const YaspGrid<dim> >*> theleafindexset;
3192  std::vector<YaspGlobalIdSet<const YaspGrid<dim> >*> theglobalidset;
3193  int nBSegments;
3194 
3195  // Index classes need access to the real entity
3196  friend class Dune::YaspLevelIndexSet<const Dune::YaspGrid<dim> >;
3197  friend class Dune::YaspLeafIndexSet<const Dune::YaspGrid<dim> >;
3198  friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim> >;
3199 
3200  friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim> >;
3201  friend class Dune::YaspIntersection<const Dune::YaspGrid<dim> >;
3202  friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim> >;
3203 
3204  template<class T>
3205  void deallocatePointers(T& container)
3206  {
3207  typedef typename T::iterator Iterator;
3208 
3209  for(Iterator entry=container.begin(); entry != container.end(); ++entry)
3210  delete (*entry);
3211  }
3212 
3213  template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
3214  friend class Entity;
3215 
3216  template<class DT>
3217  class MessageBuffer {
3218  public:
3219  // Constructor
3220  MessageBuffer (DT *p)
3221  {
3222  a=p;
3223  i=0;
3224  j=0;
3225  }
3226 
3227  // write data to message buffer, acts like a stream !
3228  template<class Y>
3229  void write (const Y& data)
3230  {
3231  dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
3232  a[i++] = data;
3233  }
3234 
3235  // read data from message buffer, acts like a stream !
3236  template<class Y>
3237  void read (Y& data) const
3238  {
3239  dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
3240  data = a[j++];
3241  }
3242 
3243  private:
3244  DT *a;
3245  int i;
3246  mutable int j;
3247  };
3248 
3249  void setsizes ()
3250  {
3251  for (YGLI g=MultiYGrid<dim,ctype>::begin(); g!=MultiYGrid<dim,ctype>::end(); ++g)
3252  {
3253  // codim 0 (elements)
3254  sizes[g.level()][0] = 1;
3255  for (int i=0; i<dim; ++i)
3256  sizes[g.level()][0] *= g.cell_overlap().size(i);
3257 
3258  // codim 1 (faces)
3259  if (dim>1)
3260  {
3261  sizes[g.level()][1] = 0;
3262  for (int i=0; i<dim; ++i)
3263  {
3264  int s=g.cell_overlap().size(i)+1;
3265  for (int j=0; j<dim; ++j)
3266  if (j!=i)
3267  s *= g.cell_overlap().size(j);
3268  sizes[g.level()][1] += s;
3269  }
3270  }
3271 
3272  // codim dim-1 (edges)
3273  if (dim>2)
3274  {
3275  sizes[g.level()][dim-1] = 0;
3276  for (int i=0; i<dim; ++i)
3277  {
3278  int s=g.cell_overlap().size(i);
3279  for (int j=0; j<dim; ++j)
3280  if (j!=i)
3281  s *= g.cell_overlap().size(j)+1;
3282  sizes[g.level()][dim-1] += s;
3283  }
3284  }
3285 
3286  // codim dim (vertices)
3287  sizes[g.level()][dim] = 1;
3288  for (int i=0; i<dim; ++i)
3289  sizes[g.level()][dim] *= g.vertex_overlapfront().size(i);
3290  }
3291  }
3292 
3294  template<int cd, PartitionIteratorType pitype>
3295  YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
3296  {
3297  dune_static_assert( cd == dim || cd == 0 ,
3298  "YaspGrid only supports Entities with codim=dim and codim=0");
3299  YGLI g = MultiYGrid<dim,ctype>::begin(level);
3300  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3301  if (pitype==Ghost_Partition)
3302  return levelend <cd, pitype> (level);
3303  if (cd==0) // the elements
3304  {
3305  if (pitype<=InteriorBorder_Partition)
3306  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubbegin());
3307  if (pitype<=All_Partition)
3308  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubbegin());
3309  }
3310  if (cd==dim) // the vertices
3311  {
3312  if (pitype==Interior_Partition)
3313  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubbegin());
3314  if (pitype==InteriorBorder_Partition)
3315  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubbegin());
3316  if (pitype==Overlap_Partition)
3317  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubbegin());
3318  if (pitype<=All_Partition)
3319  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubbegin());
3320  }
3321  DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
3322  }
3323 
3325  template<int cd, PartitionIteratorType pitype>
3326  YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
3327  {
3328  dune_static_assert( cd == dim || cd == 0 ,
3329  "YaspGrid only supports Entities with codim=dim and codim=0");
3330  YGLI g = MultiYGrid<dim,ctype>::begin(level);
3331  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3332  if (cd==0) // the elements
3333  {
3334  if (pitype<=InteriorBorder_Partition)
3335  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubend());
3336  if (pitype<=All_Partition || pitype == Ghost_Partition)
3337  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubend());
3338  }
3339  if (cd==dim) // the vertices
3340  {
3341  if (pitype==Interior_Partition)
3342  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubend());
3343  if (pitype==InteriorBorder_Partition)
3344  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubend());
3345  if (pitype==Overlap_Partition)
3346  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubend());
3347  if (pitype<=All_Partition || pitype == Ghost_Partition)
3348  return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubend());
3349  }
3350  DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
3351  }
3352 
3353  int sizes[MAXL][dim+1]; // total number of entities per level and codim
3354  bool keep_ovlp;
3355  int adaptRefCount;
3356  bool adaptActive;
3357 };
3358 
3359 namespace Capabilities
3360 {
3361 
3373  template<int dim>
3375  {
3376  static const bool v = true;
3377  static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
3378  };
3379 
3383  template<int dim>
3384  struct isCartesian< YaspGrid<dim> >
3385  {
3386  static const bool v = true;
3387  };
3388 
3392  template<int dim>
3393  struct hasEntity< YaspGrid<dim>, 0 >
3394  {
3395  static const bool v = true;
3396  };
3397 
3401  template<int dim>
3402  struct hasEntity< YaspGrid<dim>, dim >
3403  {
3404  static const bool v = true;
3405  };
3406 
3407  template< int dim >
3408  struct canCommunicate< YaspGrid< dim >, 0 >
3409  {
3410  static const bool v = true;
3411  };
3412 
3413  template< int dim >
3414  struct canCommunicate< YaspGrid< dim >, dim >
3415  {
3416  static const bool v = true;
3417  };
3418 
3422  template<int dim>
3423  struct isParallel< YaspGrid<dim> >
3424  {
3425  static const bool v = true;
3426  };
3427 
3431  template<int dim>
3433  {
3434  static const bool v = true;
3435  };
3436 
3440  template<int dim>
3442  {
3443  static const bool v = true;
3444  };
3445 
3446 }
3447 
3448 } // end namespace
3449 
3450 
3451 #endif
3452