dune-grid  2.2.0
sgrid.hh
Go to the documentation of this file.
1 #ifndef DUNE_SGRID_HH
2 #define DUNE_SGRID_HH
3 
4 #include <limits>
5 #include <vector>
6 #include <stack>
7 
8 #include <dune/common/fvector.hh>
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/bigunsignedint.hh>
11 #include <dune/common/collectivecommunication.hh>
12 #include <dune/common/reservedvector.hh>
13 #include <dune/geometry/genericgeometry/topologytypes.hh>
15 #include <dune/grid/common/grid.hh>
16 #include <dune/grid/sgrid/numbering.hh>
18 
24 namespace Dune {
25 
26 //************************************************************************
30  typedef double sgrid_ctype;
31 
32  // globally define the persistent index type
33  const int sgrid_dim_bits = 24; // bits for encoding each dimension
34  const int sgrid_level_bits = 6; // bits for encoding level number
35  const int sgrid_codim_bits = 4; // bits for encoding codimension
36 
37 //************************************************************************
38 // forward declaration of templates
39 
40 template<int dim, int dimworld, class GridImp> class SGeometry;
41 template<int codim, int dim, class GridImp> class SEntity;
42 template<int codim, class GridImp> class SEntityPointer;
43 template<int codim, class GridImp> class SEntitySeed;
44 template<int codim, PartitionIteratorType, class GridImp> class SLevelIterator;
45 template<int dim, int dimworld, class ctype> class SGrid;
46 template<class GridImp> class SIntersection;
47 template<class GridImp> class SIntersectionIterator;
48 template<class GridImp> class SHierarchicIterator;
49 
50 namespace FacadeOptions
51 {
52 
53  template<int dim, int dimworld, class ctype, int mydim, int cdim>
54  struct StoreGeometryReference<mydim, cdim,
55  SGrid<dim,dimworld,ctype>, SGeometry>
56  {
57  static const bool v = false;
58  };
59 
60  template<int dim, int dimworld, class ctype, int mydim, int cdim>
61  struct StoreGeometryReference<mydim, cdim,
62  const SGrid<dim,dimworld,ctype>, SGeometry>
63  {
64  static const bool v = false;
65  };
66 
67 }
68 
69 //************************************************************************
100 template<int mydim, int cdim, class GridImp>
102 : public GeometryDefaultImplementation<mydim,cdim,GridImp,SGeometry>
103 {
104 public:
106  typedef typename GridImp::ctype ctype;
107 
110  {
111  static const GeometryType cubeType(GeometryType::cube,mydim);
112  return cubeType;
113  }
114 
116  bool affine() const { return true ; }
117 
119  int corners () const
120  {
121  return 1<<mydim;
122  }
123 
125  FieldVector< ctype, cdim > corner ( const int i ) const
126  {
127  return c[i];
128  }
129 
131  FieldVector<ctype, cdim > center ( ) const
132  {
133  return centroid;
134  }
135 
137  FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const;
138 
140  FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const;
141 
161  ctype integrationElement (const FieldVector<ctype, mydim>& local) const
162  {
163  return volume();
164  }
165 
167  ctype volume() const;
168 
169  const FieldMatrix<ctype, mydim, cdim > &jacobianTransposed ( const FieldVector< ctype, mydim > &local ) const;
170  const FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const;
171 
173  void print (std::ostream& ss, int indent) const;
174 
179  void make (FieldMatrix<ctype,mydim+1,cdim>& __As);
180 
182  SGeometry () : builtinverse(false) {};
183 
184 private:
185  FieldVector<ctype, cdim> s;
186  FieldVector<ctype, cdim> centroid;
187  FieldMatrix<ctype,mydim,cdim> A;
188  array<FieldVector<ctype, cdim>, 1<<mydim> c;
189  mutable FieldMatrix<ctype,cdim,mydim> Jinv;
190  mutable bool builtinverse;
191 };
192 
194 template<int cdim, class GridImp>
195 class SGeometry<0,cdim,GridImp>
196 : public GeometryDefaultImplementation<0,cdim,GridImp,SGeometry>
197 {
198 public:
200  typedef typename GridImp::ctype ctype;
201 
204  {
205  static const GeometryType cubeType(GeometryType::cube,0);
206  return cubeType;
207  }
208 
210  bool affine() const { return true ; }
211 
213  int corners () const
214  {
215  return 1;
216  }
217 
219  FieldVector<ctype, cdim > corner ( const int i ) const
220  {
221  return s;
222  }
223 
225  FieldVector<ctype, cdim > center ( ) const
226  {
227  return s;
228  }
229 
231  void print (std::ostream& ss, int indent) const;
232 
234  void make (FieldMatrix<ctype,1,cdim>& __As);
235 
237  FieldVector<ctype, cdim> global (const FieldVector<ctype, 0>& local) const { return corner(0); }
238 
240  FieldVector<ctype, 0> local (const FieldVector<ctype, cdim>& global) const { return FieldVector<ctype,0> (0.0); }
241 
263 
264  SGeometry () {};
265 
271  ctype volume() const
272  {
273  return 1;
274  }
275 
281  ctype integrationElement(const FieldVector<ctype, 0>& local) const {
282  return volume();
283  }
284 
285  const FieldMatrix<ctype, 0, cdim > &jacobianTransposed ( const FieldVector<ctype, 0 > &local ) const
286  {
287  static const FieldMatrix<ctype, 0, cdim > dummy ( ctype( 0 ) );
288  return dummy;
289  }
290 
291  const FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
292  {
293  static const FieldMatrix<ctype,cdim,0> dummy( ctype(0) );
294  return dummy;
295  }
296 
297 protected:
298  FieldVector<ctype, cdim> s;
299 };
300 
301 template <int mydim, int cdim, class GridImp>
302 inline std::ostream& operator<< (std::ostream& s, SGeometry<mydim,cdim,GridImp>& e)
303 {
304  e.print(s,0);
305  return s;
306 }
307 
308 //************************************************************************
313 template<int codim, int dim, class GridImp, template<int,int,class> class EntityImp>
314 class SEntityBase :
315  public EntityDefaultImplementation<codim,dim,GridImp,EntityImp>
316 {
317  friend class SEntityPointer<codim,GridImp>;
318  friend class SIntersectionIterator<GridImp>;
319  enum { dimworld = GridImp::dimensionworld };
320 
321  typedef typename GridImp::Traits::template Codim< codim >::GeometryImpl GeometryImpl;
322 
323 public:
324  typedef typename GridImp::ctype ctype;
325  typedef typename GridImp::template Codim<codim>::Geometry Geometry;
326  typedef typename GridImp::PersistentIndexType PersistentIndexType;
327 
329  int level () const
330  {
331  return l;
332  }
333 
335  int globalIndex() const;
336 
342  }
343 
346  {
347  static const GeometryType cubeType(GeometryType::cube,dim-codim);
348  return cubeType;
349  }
350 
353  {
354  if (!builtgeometry) makegeometry();
355 
356  // return result
357  return Geometry( geo );
358  }
359 
361 
363  SEntityBase (GridImp* _grid, int _l, int _index) :
364  grid(_grid),
365  l(_l),
366  index(_index),
367  z(grid->z(l,index,codim)),
368  builtgeometry(false) {}
369 
372  builtgeometry(false) // mark geometry as not built
373  {}
374 
376  SEntityBase ( const SEntityBase& other ) :
377  grid(other.grid),
378  l(other.l),
379  index(other.index),
380  z(other.z),
381  geo(), // do not copy geometry
382  builtgeometry(false) // mark geometry as not built
383  {}
384 
386  void make (GridImp* _grid, int _l, int _id);
387 
389  void make (int _l, int _id);
390 
392  void makegeometry () const;
393 
396  {
397  return grid->persistentIndex(l, codim, z);
398  }
399 
401  int compressedIndex () const
402  {
403  return index;
404  }
405 
407  int compressedLeafIndex () const
408  {
409  // codim != dim -> there are no copies of entities
410  // maxlevel -> ids are fine
411  if (codim<dim || l==grid->maxLevel())
412  return compressedIndex();
413 
414  // this is a vertex which is not on the finest level
415  // move coordinates up to maxlevel (multiply by 2 for each level
416  array<int,dim> coord;
417  for (int k=0; k<dim; k++)
418  coord[k] = z[k]*(1<<(grid->maxLevel()-l));
419 
420  // compute number with respect to maxLevel
421  return grid->n(grid->maxLevel(),coord);
422  }
423 
424 protected:
425  // this is how we implement our elements
426  GridImp* grid;
427  int l;
428  int index;
429  array<int,dim> z;
430  mutable GeometryImpl geo;
431  mutable bool builtgeometry;
432 };
433 
434 
440 template<int codim, int dim, class GridImp>
441 class SEntity : public SEntityBase<codim,dim,GridImp,SEntity>
442 {
444  friend class SEntityPointer<codim,GridImp>;
445  friend class SIntersectionIterator<GridImp>;
446 public:
448  SEntity (GridImp* _grid, int _l, int _id) :
449  SEntityBase(_grid,_l,_id) {};
450 };
451 
478 template<int dim, class GridImp>
479 class SEntity<0,dim,GridImp> : public SEntityBase<0,dim,GridImp,SEntity>
480 {
481  enum { dimworld = GridImp::dimensionworld };
483  using SEntityBase::grid;
484  using SEntityBase::l;
485  using SEntityBase::index;
486  using SEntityBase::z;
487 
488  typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
489  typedef typename GridImp::Traits::template Codim< 0 >::LocalGeometryImpl LocalGeometryImpl;
490 
491  friend class SEntityPointer<0,GridImp>;
492  friend class SIntersectionIterator<GridImp>;
493 
494 public:
495  typedef typename GridImp::ctype ctype;
496  typedef typename GridImp::template Codim<0>::Geometry Geometry;
497  typedef typename GridImp::template Codim<0>::LocalGeometry LocalGeometry;
498  template <int cd>
499  struct Codim
500  {
501  typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
502  };
503  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
504  typedef typename GridImp::LeafIntersectionIterator IntersectionIterator;
505  typedef typename GridImp::HierarchicIterator HierarchicIterator;
506  typedef typename GridImp::PersistentIndexType PersistentIndexType;
507 
509  friend class SHierarchicIterator<GridImp>;
510 
515  template<int cc> int count () const;
516 
521  template<int cc> typename Codim<cc>::EntityPointer subEntity (int i) const;
522 
524  int subCompressedIndex (int codim, int i) const
525  {
526  if (codim==0) return this->compressedIndex();
527  // compute subIndex
528  return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
529  }
530 
534  int subCompressedLeafIndex (int codim, int i) const
535  {
536  if (codim==0) return this->compressedLeafIndex();
537 
538  assert(this->l == this->grid->maxLevel());
539  // compute subIndex
540  return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
541  }
542 
544  PersistentIndexType subPersistentIndex (int codim, int i) const
545  {
546  if (codim==0) return this->persistentIndex();
547  // compute subId
548  return this->grid->persistentIndex(this->l, codim, this->grid->subz(this->z,i,codim));
549  }
550 
558  IntersectionIterator ibegin () const;
559  IntersectionIterator ileafbegin () const;
560  IntersectionIterator ilevelbegin () const;
562  IntersectionIterator iend () const;
563  IntersectionIterator ileafend () const;
564  IntersectionIterator ilevelend () const;
565 
571  EntityPointer father () const;
572 
574  bool hasFather () const
575  {
576  return (this->level()>0);
577  }
578 
580  bool isLeaf () const
581  {
582  return ( this->grid->maxLevel() == this->level() );
583  }
584 
596  LocalGeometry geometryInFather () const;
597 
604  HierarchicIterator hbegin (int maxLevel) const;
605 
607  HierarchicIterator hend (int maxLevel) const;
608 
609  // members specific to SEntity
611  SEntity (GridImp* _grid, int _l, int _index) :
612  SEntityBase(_grid,_l,_index),
613  built_father(false)
614  {}
615 
616  SEntity (const SEntity& other ) :
617  SEntityBase(other.grid, other.l, other.index ),
618  built_father(false)
619  {}
620 
622  void make (GridImp* _grid, int _l, int _id)
623  {
624  SEntityBase::make(_grid,_l,_id);
625  built_father = false;
626  }
627 
629  void make (int _l, int _id)
630  {
631  SEntityBase::make(_l,_id);
632  built_father = false;
633  }
634 
635 private:
636 
637  SEntity();
638 
639  mutable bool built_father;
640  mutable int father_index;
641  mutable LocalGeometryImpl in_father_local;
642  void make_father() const;
643 };
644 
645 
646 //************************************************************************
655  int l;
656  int index;
657  SHierarchicStackElem () : l(-1), index(-1) {}
658  SHierarchicStackElem (int _l, int _index) {l=_l; index=_index;}
659  bool operator== (const SHierarchicStackElem& s) const {return !operator!=(s);}
660  bool operator!= (const SHierarchicStackElem& s) const {return l!=s.l || index!=s.index;}
661 };
662 
663 template<class GridImp>
665  public Dune::SEntityPointer <0,GridImp>
666 {
667  friend class SHierarchicIterator<const GridImp>;
668  enum { dim = GridImp::dimension };
669  enum { dimworld = GridImp::dimensionworld };
670  typedef Dune::SEntityPointer<0,GridImp> SEntityPointer;
672  using SEntityPointer::grid;
673  using SEntityPointer::l;
674  using SEntityPointer::index;
675 public:
676  typedef typename GridImp::template Codim<0>::Entity Entity;
677  typedef typename GridImp::ctype ctype;
678 
680  void increment();
681 
688  SHierarchicIterator (GridImp* _grid,
690  int _maxLevel, bool makeend) :
691  SEntityPointer(_grid,_e.level(),_e.compressedIndex())
692  {
693  // without sons, we are done
694  // (the end iterator is equal to the calling iterator)
695  if (makeend) return;
696 
697  // remember element where begin has been called
698  orig_l = this->entity().level();
699  orig_index = _grid->getRealImplementation(this->entity()).compressedIndex();
700 
701  // push original element on stack
702  SHierarchicStackElem originalElement(orig_l, orig_index);
703  stack.push(originalElement);
704 
705  // compute maxLevel
706  maxLevel = std::min(_maxLevel,this->grid->maxLevel());
707 
708  // ok, push all the sons as well
709  push_sons(orig_l,orig_index);
710 
711  // and pop the first son
712  increment();
713  }
714 
715 private:
716  int maxLevel;
717  int orig_l, orig_index;
718 
720  std::stack<SHierarchicStackElem, Dune::ReservedVector<SHierarchicStackElem,GridImp::MAXL> > stack;
721 
722  void push_sons (int level, int fatherid);
723 };
724 
725 //************************************************************************
732 template<class GridImp>
734 {
735  enum { dim=GridImp::dimension };
736  enum { dimworld=GridImp::dimensionworld };
737 
738  typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
739  typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
740 
741 public:
742  typedef typename GridImp::template Codim<0>::Entity Entity;
743  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
744  typedef typename GridImp::template Codim<1>::Geometry Geometry;
745  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
749  enum { dimension=dim };
751  enum { dimensionworld=dimworld };
753  typedef typename GridImp::ctype ctype;
754 
756  bool equals(const SIntersectionIterator<GridImp>& i) const;
758  void increment();
759 
761  const Intersection & dereference() const
762  {
763  return intersection;
764  }
765 
768  EntityPointer inside() const;
769 
772  EntityPointer outside() const;
773 
775  bool boundary () const;
776 
778  bool conforming () const;
779 
780  int boundaryId () const {
781  if (boundary()) return count + 1;
782  return 0;
783  };
784 
785  int boundarySegmentIndex () const {
786  if (boundary())
787  return grid->boundarySegmentIndex(self.level(), count, zred);
788  return -1;
789  };
790 
792  bool neighbor () const;
793 
795  FieldVector<ctype, GridImp::dimensionworld> outerNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
796  {
797  return unitOuterNormal(local);
798  }
800  FieldVector<ctype, GridImp::dimensionworld> unitOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
801  {
802  return centerUnitOuterNormal();
803  }
805  FieldVector<ctype, GridImp::dimensionworld> centerUnitOuterNormal () const
806  {
807  // while we are at it, compute normal direction
808  FieldVector<ctype, dimworld> normal(0.0);
809  if (count%2)
810  normal[count/2] = 1.0; // odd
811  else
812  normal[count/2] = -1.0; // even
813 
814  return normal;
815  }
817  FieldVector<ctype, GridImp::dimensionworld> integrationOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
818  {
819  FieldVector<ctype, dimworld> n = unitOuterNormal(local);
820  n *= geometry().integrationElement(local);
821  return n;
822  }
823 
835  Geometry geometry () const;
836 
839  {
840  static const GeometryType cubeType(GeometryType::cube,dim-1);
841  return cubeType;
842  }
843 
845  int indexInInside () const;
847  int indexInOutside () const;
848 
850  SIntersectionIterator (GridImp* _grid, const SEntity<0,dim,GridImp >* _self, int _count) :
851  self(*_self), ne(self), grid(_grid),
852  partition(_grid->partition(grid->getRealImplementation(ne).l,_self->z)),
853  zred(_grid->compress(grid->getRealImplementation(ne).l,_self->z)),
854  intersection(IntersectionImp(*this))
855  {
856  // make neighbor
857  make(_count);
858  }
859 
861  self(other.self), ne(other.ne), grid(other.grid),
862  partition(other.partition), zred(other.zred),
863  count(other.count), valid_count(other.valid_count),
864  valid_nb(other.valid_nb), is_on_boundary(other.is_on_boundary),
865  built_intersections(false),
866  intersection(IntersectionImp(*this))
867  {
868  }
869 
872  {
873  /* We can't assign the grid */
874  assert(grid == other.grid);
875 
876  /* Assign data from other */
877  self = other.self;
878  ne = other.ne;
879  partition = other.partition;
880  zred = other.zred;
881  count = other.count;
882  valid_count = other.valid_count;
883  valid_nb = other.valid_nb;
884  is_on_boundary = other.is_on_boundary;
885 
886  /* mark cached data as invalid */
887  built_intersections = false;
888 
889  return *this;
890  }
891 
892 private:
893  void make (int _count) const;
894  void makeintersections () const;
895  EntityPointer self;
896  mutable EntityPointer ne;
897  const GridImp * grid;
898  int partition;
899  array<int,dim> zred;
900  mutable int count;
901  mutable bool valid_count;
902  mutable bool valid_nb;
903  mutable bool is_on_boundary;
904  mutable bool built_intersections;
905  mutable LocalGeometryImpl is_self_local;
906  mutable GeometryImpl is_global;
907  mutable LocalGeometryImpl is_nb_local;
908  Intersection intersection;
909 };
910 
911 template<class GridImp>
913 {
914  enum { dim=GridImp::dimension };
915  enum { dimworld=GridImp::dimensionworld };
916 public:
917  typedef typename GridImp::template Codim<0>::Entity Entity;
918  typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
919  typedef typename GridImp::template Codim<1>::Geometry Geometry;
922  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
925  enum { dimension=dim };
927  enum { dimensionworld=dimworld };
929  typedef typename GridImp::ctype ctype;
930 
931  bool boundary () const
932  {
933  return is.boundary();
934  }
935 
937  int boundaryId () const
938  {
939  return is.boundaryId();
940  }
941 
943  size_t boundarySegmentIndex () const
944  {
945  return is.boundarySegmentIndex();
946  }
947 
949  bool neighbor () const
950  {
951  return is.neighbor();
952  }
953 
956  {
957  return is.inside();
958  }
959 
962  {
963  return is.outside();
964  }
965 
967  bool conforming () const
968  {
969  return is.conforming();
970  }
971 
974  {
975  return is.geometryInInside();
976  }
977 
980  {
981  return is.geometryInOutside();
982  }
983 
986  {
987  return is.geometry();
988  }
989 
992  {
993  return is.type();
994  }
995 
997  int indexInInside () const
998  {
999  return is.indexInInside();
1000  }
1001 
1003  int indexInOutside () const
1004  {
1005  return is.indexInOutside();
1006  }
1007 
1010  {
1011  return is.outerNormal(local);
1012  }
1013 
1016  {
1017  return is.integrationOuterNormal(local);
1018  }
1019 
1022  {
1023  return is.unitOuterNormal(local);
1024  }
1025 
1028  {
1029  return is.centerUnitOuterNormal();
1030  }
1031 
1034 
1035 private:
1036 #ifndef DOXYGEN // doxygen can't handle this recursive usage
1037  const SIntersectionIterator<GridImp> & is;
1038 #endif
1039 };
1040 
1041 //************************************************************************
1042 
1046 template <class T>
1047 class AutoPtrStack : public std::stack<T*>
1048 {
1049 public:
1051  {
1052  while(! this->empty())
1053  {
1054  T* e = this->top();
1055  delete e;
1056  this->pop();
1057  }
1058  }
1059 };
1060 
1063 template<int codim, class GridImp>
1065 {
1066  enum { dim = GridImp::dimension };
1067  friend class SIntersectionIterator<GridImp>;
1068 public:
1070  typedef typename GridImp::template Codim<codim>::Entity Entity;
1072  enum { codimension = codim };
1073 
1075  bool equals(const SEntityPointer<codim,GridImp>& i) const;
1077  Entity& dereference() const;
1079  int level () const;
1080 
1082  SEntityPointer (GridImp * _grid, int _l, int _index) :
1083  grid(_grid), l(_l), index(_index),
1084  e(0)
1085  {
1086  }
1087 
1090  grid(_e.grid), l(_e.l), index(_e.index),
1091  e(0)
1092  {
1093  }
1094 
1097  grid(other.grid), l(other.l), index(other.index),
1098  e( 0 )
1099  {
1100  }
1101 
1104  {
1105  if( e )
1106  enStack().push( e );
1107 #ifndef NDEBUG
1108  index = -1;
1109 #endif
1110  }
1111 
1114  {
1115  grid = other.grid;
1116  l = other.l;
1117  index = other.index;
1118 
1119  // free current entity
1120  if( e )
1121  enStack().push( e );
1122  e = 0;
1123 
1124  return *this;
1125  }
1126 
1127 protected:
1129  {
1130  return grid->getRealImplementation(entity());
1131  }
1132 
1133  inline Entity& entity() const
1134  {
1135  if( ! e )
1136  {
1137  e = getEntity( grid, l, index );
1138  }
1139  return *e;
1140  }
1141 
1143  static inline EntityStackType& enStack()
1144  {
1145  static EntityStackType eStack;
1146  return eStack;
1147  }
1148 
1149  inline Entity* getEntity(GridImp* _grid, int _l, int _id ) const
1150  {
1151  // get stack reference
1152  EntityStackType& enSt = enStack();
1153 
1154  if( enSt.empty() )
1155  {
1156  return (new Entity(SEntity<codim,dim,GridImp>(_grid, _l, _id)));
1157  }
1158  else
1159  {
1160  Entity* e = enSt.top();
1161  enSt.pop();
1162  grid->getRealImplementation(*e).make(_grid, _l,_id);
1163  return e;
1164  }
1165  }
1166 
1167  GridImp* grid;
1168  int l;
1169  mutable int index;
1170  mutable Entity* e;
1171 };
1172 
1175 template<int codim, class GridImp>
1177 {
1178  enum { dim = GridImp::dimension };
1179 public:
1180  enum { codimension = codim };
1181 
1183  SEntitySeed (int l, int index) :
1184  _l(l), _index(index)
1185  {}
1186 
1187  int level () const { return this->_l; }
1188  int index () const { return this->_index; }
1189 
1190 private:
1191  int _l;
1192  int _index;
1193 };
1194 
1195 //************************************************************************
1196 
1197 
1200 template<int codim, PartitionIteratorType pitype, class GridImp>
1202  public Dune::SEntityPointer <codim,GridImp>
1203 {
1204  friend class SLevelIterator<codim, pitype,const GridImp>;
1205  enum { dim = GridImp::dimension };
1208  using SEntityPointer::l;
1209  using SEntityPointer::index;
1210 public:
1211  typedef typename GridImp::template Codim<codim>::Entity Entity;
1212 
1214  void increment();
1215 
1217  SLevelIterator (GridImp * _grid, int _l, int _id) :
1218  SEntityPointer(_grid,_l,_id) {}
1219 };
1220 
1221 
1222 //========================================================================
1227 //========================================================================
1228 
1229 template<class GridImp>
1230 class SGridLevelIndexSet : public IndexSet<GridImp,SGridLevelIndexSet<GridImp> >
1231 {
1234 
1235  enum { dim = GridImp::dimension };
1236 
1237 public:
1238 
1240  SGridLevelIndexSet ( const GridImp &g, int l )
1241  : grid( g ),
1242  level( l )
1243  {
1244  // TODO move list of geometrytypes to grid, can be computed static (singleton)
1245  // contains a single element type;
1246  for (int codim=0; codim<=GridImp::dimension; codim++)
1247  mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
1248  }
1249 
1251  template<int cd>
1252  int index (const typename GridImp::Traits::template Codim<cd>::Entity& e) const
1253  {
1254  return grid.getRealImplementation(e).compressedIndex();
1255  }
1256 
1257  template< int cc >
1258  int subIndex ( const typename GridImp::Traits::template Codim< cc >::Entity &e,
1259  int i, unsigned int codim ) const
1260  {
1261  if( cc == 0 )
1262  return grid.getRealImplementation(e).subCompressedIndex(codim, i);
1263  else
1264  DUNE_THROW( NotImplemented, "subIndex for higher codimension entity not implemented for SGrid." );
1265  }
1266 
1267  // return true if the given entity is contained in \f$E\f$.
1268  template< class EntityType >
1269  bool contains ( const EntityType &e ) const
1270  {
1271  return (e.level() == level);
1272  }
1273 
1275  int size (GeometryType type) const
1276  {
1277  return grid.size( level, type );
1278  }
1279 
1281  int size (int codim) const
1282  {
1283  return grid.size( level, codim );
1284  }
1285 
1287  const std::vector<GeometryType>& geomTypes (int codim) const
1288  {
1289  return mytypes[codim];
1290  }
1291 
1292 private:
1293  const GridImp& grid;
1294  int level;
1295  std::vector<GeometryType> mytypes[GridImp::dimension+1];
1296 };
1297 
1298 // Leaf Index Set
1299 
1300 template<class GridImp>
1301 class SGridLeafIndexSet : public IndexSet<GridImp,SGridLeafIndexSet<GridImp> >
1302 {
1305 
1306  enum { dim = GridImp::dimension };
1307 
1308 public:
1309 
1311  explicit SGridLeafIndexSet ( const GridImp &g )
1312  : grid( g )
1313  {
1314  // TODO move list of geometrytypes to grid, can be computed static (singleton)
1315  // contains a single element type;
1316  for (int codim=0; codim<=dim; codim++)
1317  mytypes[codim].push_back(GeometryType(GeometryType::cube,dim-codim));
1318  }
1319 
1321  /*
1322  We use the remove_const to extract the Type from the mutable class,
1323  because the const class is not instantiated yet.
1324  */
1325  template<int cd>
1326  int index (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
1327  {
1328  return grid.getRealImplementation(e).compressedLeafIndex();
1329  }
1330 
1331  template< int cc >
1332  int subIndex ( const typename GridImp::Traits::template Codim< cc >::Entity &e,
1333  int i, unsigned int codim ) const
1334  {
1335  if( cc == 0 )
1336  return grid.getRealImplementation(e).subCompressedIndex(codim, i);
1337  else
1338  DUNE_THROW( NotImplemented, "subIndex for higher codimension entity not implemented for SGrid." );
1339  }
1340 
1342  int size (GeometryType type) const
1343  {
1344  return grid.size( grid.maxLevel(), type );
1345  }
1346 
1348  int size (int codim) const
1349  {
1350  return grid.size( grid.maxLevel(), codim );
1351  }
1352 
1353  // return true if the given entity is contained in \f$E\f$.
1354  template< class EntityType >
1355  bool contains ( const EntityType &e ) const
1356  {
1357  return (e.level() == grid.maxLevel());
1358  }
1359 
1361  const std::vector<GeometryType>& geomTypes (int codim) const
1362  {
1363  return mytypes[codim];
1364  }
1365 
1366 private:
1367  const GridImp& grid;
1368  std::vector<GeometryType> mytypes[dim+1];
1369 };
1370 
1371 
1372 //========================================================================
1377 //========================================================================
1378 
1379 template<class GridImp>
1381  public IdSet<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>
1382  /*
1383  We used the remove_const to extract the Type from the mutable class,
1384  because the const class is not instantiated yet.
1385  */
1386 {
1388 
1389 public:
1390 
1393  using IdSet<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>::subId;
1394 
1396  /*
1397  We use the remove_const to extract the Type from the mutable class,
1398  because the const class is not instantiated yet.
1399  */
1400  typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;
1401 
1403  explicit SGridGlobalIdSet ( const GridImp &g )
1404  : grid( g )
1405  {}
1406 
1408  /*
1409  We use the remove_const to extract the Type from the mutable class,
1410  because the const class is not instantiated yet.
1411  */
1412  template<int cd>
1413  IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
1414  {
1415  return grid.getRealImplementation(e).persistentIndex();
1416  }
1417 
1419  /*
1420  We use the remove_const to extract the Type from the mutable class,
1421  because the const class is not instantiated yet.
1422  */
1423  IdType subId ( const typename remove_const< GridImp >::type::Traits::template Codim< 0 >::Entity &e,
1424  int i, unsigned int codim ) const
1425  {
1426  return grid.getRealImplementation(e).subPersistentIndex(codim, i);
1427  }
1428 
1429 private:
1430  const GridImp& grid;
1431 };
1432 
1433 
1434 template<int dim, int dimworld, class ctype>
1436 {
1440  SIntersection, // leaf intersection
1441  SIntersection, // level intersection
1442  SIntersectionIterator, // leaf intersection iter
1443  SIntersectionIterator, // level intersection iter
1449  bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
1451  bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
1452  CollectiveCommunication<Dune::SGrid<dim,dimworld,ctype> >,
1454  SEntitySeed>
1456 };
1457 
1458 
1459 //************************************************************************
1513 template<int dim, int dimworld, typename _ctype = sgrid_ctype>
1514 class SGrid : public GridDefaultImplementation <dim,dimworld,_ctype,SGridFamily<dim,dimworld,_ctype> >
1515 {
1516 public:
1518  typedef bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits> PersistentIndexType;
1519 
1520  // need for friend declarations in entity
1524 
1526 
1528  enum { MAXL=32 };
1529 
1531  typedef _ctype ctype;
1532 
1533  // constructors
1534 
1542  SGrid (const int * const N_, const ctype * const H_);
1543 
1551  SGrid (const int * const N_, const ctype * const L_, const ctype * const H_);
1552 
1562  SGrid (FieldVector<int,dim> N_, FieldVector<ctype,dim> L_, FieldVector<ctype,dim> H_);
1563 
1565  SGrid ();
1566 
1568  ~SGrid ();
1569 
1572  int maxLevel() const;
1573 
1575  template<int cd, PartitionIteratorType pitype>
1576  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const;
1577 
1579  template<int cd, PartitionIteratorType pitype>
1580  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const;
1581 
1583  template<int cd>
1584  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
1585  {
1586  return lbegin<cd,All_Partition>(level);
1587  }
1588 
1590  template<int cd>
1591  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
1592  {
1593  return lend<cd,All_Partition>(level);
1594  }
1595 
1597  template<int cd, PartitionIteratorType pitype>
1598  typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const;
1599 
1601  template<int cd, PartitionIteratorType pitype>
1602  typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const;
1603 
1605  template<int cd>
1607  {
1608  return leafbegin<cd,All_Partition>();
1609  }
1610 
1612  template<int cd>
1614  {
1615  return leafend<cd,All_Partition>();
1616  }
1617 
1618  // \brief obtain EntityPointer from EntitySeed. */
1619  template <typename Seed>
1620  typename Traits::template Codim<Seed::codimension>::EntityPointer
1621  entityPointer(const Seed& seed) const
1622  {
1623  enum { codim = Seed::codimension };
1624  return SEntityPointer<codim,const SGrid<dim,dimworld> >(this,seed.level(),seed.index());
1625  }
1626 
1640  template<class T, template<class> class P, int codim>
1641  void communicate (T& t, InterfaceType iftype, CommunicationDirection dir, int level)
1642  {
1643  // SGrid is sequential and has no periodic boundaries, so do nothing ...
1644  return;
1645  }
1646 
1648  int size (int level, int codim) const;
1649 
1651  int size (int codim) const
1652  {
1653  return size(maxLevel(),codim);
1654  }
1655 
1657  int size (int level, GeometryType type) const
1658  {
1659  return (type.isCube()) ? size(level,dim-type.dim()) : 0;
1660  }
1661 
1663  int size (GeometryType type) const
1664  {
1665  return size(maxLevel(),type);
1666  }
1667 
1669  size_t numBoundarySegments () const
1670  {
1671  return boundarysize;
1672  }
1673 
1675  int global_size (int codim) const;
1676 
1678  int overlapSize (int level, int codim)
1679  {
1680  return 0;
1681  }
1682 
1684  int ghostSize (int level, int codim)
1685  {
1686  return 0;
1687  }
1688 
1689  // these are all members specific to sgrid
1690 
1692  void globalRefine (int refCount);
1693 
1695  const array<int, dim>& dims(int level) const {
1696  return N[level];
1697  }
1698 
1700  const FieldVector<ctype, dimworld>& lowerLeft() const {
1701  return low;
1702  }
1703 
1705  FieldVector<ctype, dimworld> upperRight() const {
1706  return H;
1707  }
1708 
1710  bool adapt ()
1711  {
1712  globalRefine(1);
1713  return true;
1714  }
1715 
1716  // The new index sets from DDM 11.07.2005
1717  const typename Traits::GlobalIdSet& globalIdSet() const
1718  {
1719  return *theglobalidset;
1720  }
1721 
1722  const typename Traits::LocalIdSet& localIdSet() const
1723  {
1724  return *theglobalidset;
1725  }
1726 
1727  const typename Traits::LevelIndexSet& levelIndexSet(int level) const
1728  {
1729  assert(level>=0 && level<=maxLevel());
1730  return *(indexsets[level]);
1731  }
1732 
1733  const typename Traits::LeafIndexSet& leafIndexSet() const
1734  {
1735  return *theleafindexset;
1736  }
1737 
1742  template<class DataHandle>
1743  void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
1744  {
1745  }
1746 
1747  template<class DataHandle>
1748  void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
1749  {
1750  }
1751 
1753  {
1754  return ccobj;
1755  }
1756 
1758  int overlapSize (int level, int codim) const
1759  {
1760  return 0;
1761  }
1762 
1764  int overlapSize (int codim) const
1765  {
1766  return 0;
1767  }
1768 
1770  int ghostSize (int level, int codim) const
1771  {
1772  return 0;
1773  }
1774 
1776  int ghostSize (int codim) const
1777  {
1778  return 0;
1779  }
1780 
1781  /*
1782  @}
1783  */
1784 
1785 private:
1786  /*
1787  Make associated classes friends to grant access to the real entity
1788  */
1789  friend class Dune::SGridLevelIndexSet<Dune::SGrid<dim,dimworld> >;
1790  friend class Dune::SGridLeafIndexSet<Dune::SGrid<dim,dimworld> >;
1791  friend class Dune::SGridGlobalIdSet<Dune::SGrid<dim,dimworld> >;
1792  friend class Dune::SIntersectionIterator<Dune::SGrid<dim,dimworld> >;
1793  friend class Dune::SHierarchicIterator<Dune::SGrid<dim,dimworld> >;
1794  friend class Dune::SEntity<0,dim,Dune::SGrid<dim,dimworld> >;
1795 
1796  friend class Dune::SGridLevelIndexSet<const Dune::SGrid<dim,dimworld> >;
1797  friend class Dune::SGridLeafIndexSet<const Dune::SGrid<dim,dimworld> >;
1798  friend class Dune::SGridGlobalIdSet<const Dune::SGrid<dim,dimworld> >;
1799  friend class Dune::SIntersectionIterator<const Dune::SGrid<dim,dimworld> >;
1800  friend class Dune::SHierarchicIterator<const Dune::SGrid<dim,dimworld> >;
1801  friend class Dune::SEntity<0,dim,const Dune::SGrid<dim,dimworld> >;
1802 
1803  template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1804  friend class Dune::SEntityBase;
1805 
1806  template<int codim_, class GridImp_>
1807  friend class Dune::SEntityPointer;
1808 
1809  template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1810  friend class Entity;
1811 
1813  FieldVector<ctype, dimworld> pos (int level, array<int,dim>& z) const;
1814 
1816  int calc_codim (int level, const array<int,dim>& z) const;
1817 
1819  int n (int level, const array<int,dim>& z) const;
1820 
1822  array<int,dim> z (int level, int i, int codim) const;
1823 
1825  array<int,dim> subz (const array<int,dim> & z, int i, int codim) const;
1826 
1828  array<int,dim> compress (int level, const array<int,dim>& z) const;
1829 
1831  array<int,dim> expand (int level, const array<int,dim>& r, int b) const;
1832 
1836  int partition (int level, const array<int,dim>& z) const;
1837 
1839  bool exists (int level, const array<int,dim>& zred) const;
1840 
1841  // compute boundary segment index for a given zentity and a face
1842  int boundarySegmentIndex (int l, int face, const array<int,dim> & zentity) const
1843  {
1844  array<int,dim-1> zface;
1845  int dir = face/2;
1846  int side = face%2;
1847  // compute z inside the global face
1848  for (int i=0; i<dir; i++) zface[i] = zentity[i]/(1<<l);
1849  for (int i=dir+1; i<dim; i++) zface[i-1] = zentity[i]/(1<<l);
1850  zface = boundarymapper[dir].expand(zface, 0);
1851  // compute index in the face
1852  int index = boundarymapper[dir].n(zface);
1853  // compute offset
1854  for (int i=0; i<dir; i++)
1855  index += 2*boundarymapper[i].elements(0);
1856  index += side*boundarymapper[dir].elements(0);
1857  return index;
1858  }
1859 
1860  // compute persistent index for a given zentity
1861  PersistentIndexType persistentIndex (int l, int codim, const array<int,dim> & zentity) const
1862  {
1863  if (codim!=dim)
1864  {
1865  // encode codim, this would actually not be necessary
1866  // because z is unique in codim
1867  PersistentIndexType id(codim);
1868 
1869  // encode level
1870  id = id << sgrid_level_bits;
1871  id = id+PersistentIndexType(l);
1872 
1873  // encode coordinates
1874  for (int i=dim-1; i>=0; i--)
1875  {
1876  id = id << sgrid_dim_bits;
1877  id = id+PersistentIndexType(zentity[i]);
1878  }
1879 
1880  return id;
1881  }
1882  else
1883  {
1884  // determine min number of trailing zeroes
1885  // consider that z is on the doubled grid !
1886  int trailing = 1000;
1887  for (int i=0; i<dim; i++)
1888  {
1889  // count trailing zeros
1890  int zeros = 0;
1891  for (int j=0; j<l; j++)
1892  if (zentity[i]&(1<<(j+1)))
1893  break;
1894  else
1895  zeros++;
1896  trailing = std::min(trailing,zeros);
1897  }
1898 
1899  // determine the level of this vertex
1900  int level = l-trailing;
1901 
1902  // encode codim
1903  PersistentIndexType id(dim);
1904 
1905  // encode level
1906  id = id << sgrid_level_bits;
1907  id = id+PersistentIndexType(level);
1908 
1909  // encode coordinates
1910  for (int i=dim-1; i>=0; i--)
1911  {
1912  id = id << sgrid_dim_bits;
1913  id = id+PersistentIndexType(zentity[i]>>trailing);
1914  }
1915 
1916  return id;
1917  }
1918  }
1919 
1920  // disable copy and assign
1921  SGrid(const SGrid &) {};
1922  SGrid & operator = (const SGrid &) { return *this; };
1923  // generate SGrid
1924  void makeSGrid (const array<int,dim>& N_, const FieldVector<ctype, dim>& L_, const FieldVector<ctype, dim>& H_);
1925 
1926  /*
1927  internal data
1928  */
1929  CollectiveCommunication<SGrid> ccobj;
1930 
1931  ReservedVector<SGridLevelIndexSet<const SGrid<dim,dimworld> >*, MAXL> indexsets;
1932  SGridLeafIndexSet<const SGrid<dim,dimworld> > *theleafindexset;
1933  SGridGlobalIdSet<const SGrid<dim,dimworld> > *theglobalidset;
1934 
1935  int L; // number of levels in hierarchic mesh 0<=level<L
1936  FieldVector<ctype, dim> low; // lower left corner of the grid
1937  FieldVector<ctype, dim> H; // length of cube per direction
1938  array<int,dim> *N; // number of elements per direction
1939  FieldVector<ctype, dim> *h; // mesh size per direction
1940  mutable CubeMapper<dim> *mapper;// a mapper for each level
1941 
1942  // boundary segement index set
1943  array<CubeMapper<dim-1>, dim> boundarymapper; // a mapper for each coarse grid face
1944  int boundarysize;
1945 
1946  // faster implementation of subIndex
1947  mutable array <int,dim> zrefStatic; // for subIndex of SEntity
1948  mutable array <int,dim> zentityStatic; // for subIndex of SEntity
1949 };
1950 
1951 namespace Capabilities
1952 {
1953 
1965  template<int dim, int dimw>
1966  struct hasSingleGeometryType< SGrid<dim,dimw> >
1967  {
1968  static const bool v = true;
1969  static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
1970  };
1971 
1975  template<int dim, int dimw>
1976  struct isCartesian< SGrid<dim,dimw> >
1977  {
1978  static const bool v = true;
1979  };
1980 
1984  template<int dim, int dimw, int cdim>
1985  struct hasEntity< SGrid<dim,dimw>, cdim>
1986  {
1987  static const bool v = true;
1988  };
1989 
1993  template<int dim, int dimw>
1994  struct isLevelwiseConforming< SGrid<dim,dimw> >
1995  {
1996  static const bool v = true;
1997  };
1998 
2002  template<int dim, int dimw>
2003  struct isLeafwiseConforming< SGrid<dim,dimw> >
2004  {
2005  static const bool v = true;
2006  };
2007 
2008 } // end namespace Capabilities
2009 
2010 } // end namespace Dune
2011 
2012 #include"sgrid/sgrid.cc"
2013 
2014 #endif