Dune Core Modules (2.7.1)

yaspgrid.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GRID_YASPGRID_HH
4 #define DUNE_GRID_YASPGRID_HH
5 
6 #include <iostream>
7 #include <vector>
8 #include <algorithm>
9 #include <stack>
10 #include <type_traits>
11 
12 // either include stdint.h or provide fallback for uint8_t
13 #if HAVE_STDINT_H
14 #include <stdint.h>
15 #else
16 typedef unsigned char uint8_t;
17 #endif
18 
19 #include <dune/grid/common/backuprestore.hh>
20 #include <dune/grid/common/grid.hh> // the grid base classes
21 #include <dune/grid/common/capabilities.hh> // the capabilities
22 #include <dune/common/hybridutilities.hh>
23 #include <dune/common/power.hh>
31 #include <dune/geometry/type.hh>
34 
35 
36 #if HAVE_MPI
38 #endif
39 
47 namespace Dune {
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 = 5; // bits for encoding level number
53 
54 
55  //************************************************************************
56  // forward declaration of templates
57 
58  template<int dim, class Coordinates> class YaspGrid;
59  template<int mydim, int cdim, class GridImp> class YaspGeometry;
60  template<int codim, int dim, class GridImp> class YaspEntity;
61  template<int codim, class GridImp> class YaspEntitySeed;
62  template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
63  template<class GridImp> class YaspIntersectionIterator;
64  template<class GridImp> class YaspIntersection;
65  template<class GridImp> class YaspHierarchicIterator;
66  template<class GridImp, bool isLeafIndexSet> class YaspIndexSet;
67  template<class GridImp> class YaspGlobalIdSet;
68  template<class GridImp> class YaspPersistentContainerIndex;
69 
70 } // namespace Dune
71 
83 #include <dune/grid/yaspgrid/yaspgrididset.hh>
85 
86 namespace Dune {
87 
88  template<int dim, class Coordinates>
89  struct YaspGridFamily
90  {
91 #if HAVE_MPI
92  typedef CollectiveCommunication<MPI_Comm> CCType;
93 #else
94  typedef CollectiveCommunication<No_Comm> CCType;
95 #endif
96 
97  typedef GridTraits<dim, // dimension of the grid
98  dim, // dimension of the world space
100  YaspGeometry,YaspEntity,
101  YaspLevelIterator, // type used for the level iterator
102  YaspIntersection, // leaf intersection
103  YaspIntersection, // level intersection
104  YaspIntersectionIterator, // leaf intersection iter
105  YaspIntersectionIterator, // level intersection iter
106  YaspHierarchicIterator,
107  YaspLevelIterator, // type used for the leaf(!) iterator
108  YaspIndexSet< const YaspGrid< dim, Coordinates >, false >, // level index set
109  YaspIndexSet< const YaspGrid< dim, Coordinates >, true >, // leaf index set
110  YaspGlobalIdSet<const YaspGrid<dim, Coordinates> >,
111  bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
112  YaspGlobalIdSet<const YaspGrid<dim, Coordinates> >,
113  bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
114  CCType,
115  DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
116  YaspEntitySeed>
117  Traits;
118  };
119 
120 #ifndef DOXYGEN
121  template<int dim, int codim>
122  struct YaspCommunicateMeta {
123  template<class G, class DataHandle>
124  static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
125  {
126  if (data.contains(dim,codim))
127  {
128  g.template communicateCodim<DataHandle,codim>(data,iftype,dir,level);
129  }
130  YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
131  }
132  };
133 
134  template<int dim>
135  struct YaspCommunicateMeta<dim,0> {
136  template<class G, class DataHandle>
137  static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
138  {
139  if (data.contains(dim,0))
140  g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
141  }
142  };
143 #endif
144 
145  //************************************************************************
162  template<int dim, class Coordinates = EquidistantCoordinates<double, dim> >
163  class YaspGrid
164  : public GridDefaultImplementation<dim,dim,typename Coordinates::ctype,YaspGridFamily<dim, Coordinates> >
165  {
166 
167  template<int, PartitionIteratorType, typename>
168  friend class YaspLevelIterator;
169 
170  template<typename>
171  friend class YaspHierarchicIterator;
172 
173  protected:
174 
176 
177  public:
179  typedef typename Coordinates::ctype ctype;
180 #if HAVE_MPI
181  typedef CollectiveCommunication<MPI_Comm> CollectiveCommunicationType;
182 #else
183  typedef CollectiveCommunication<No_Comm> CollectiveCommunicationType;
184 #endif
185 
186 #ifndef DOXYGEN
187  typedef typename Dune::YGrid<Coordinates> YGrid;
189 
192  struct YGridLevel {
193 
195  int level() const
196  {
197  return level_;
198  }
199 
200  Coordinates coords;
201 
202  std::array<YGrid, dim+1> overlapfront;
203  std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> overlapfront_data;
204  std::array<YGrid, dim+1> overlap;
205  std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> overlap_data;
206  std::array<YGrid, dim+1> interiorborder;
207  std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> interiorborder_data;
208  std::array<YGrid, dim+1> interior;
209  std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> interior_data;
210 
211  std::array<YGridList<Coordinates>,dim+1> send_overlapfront_overlapfront;
212  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_overlapfront_overlapfront_data;
213  std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_overlapfront;
214  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_overlapfront_data;
215 
216  std::array<YGridList<Coordinates>,dim+1> send_overlap_overlapfront;
217  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_overlap_overlapfront_data;
218  std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_overlap;
219  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_overlap_data;
220 
221  std::array<YGridList<Coordinates>,dim+1> send_interiorborder_interiorborder;
222  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_interiorborder_interiorborder_data;
223  std::array<YGridList<Coordinates>,dim+1> recv_interiorborder_interiorborder;
224  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_interiorborder_interiorborder_data;
225 
226  std::array<YGridList<Coordinates>,dim+1> send_interiorborder_overlapfront;
227  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_interiorborder_overlapfront_data;
228  std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_interiorborder;
229  std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_interiorborder_data;
230 
231  // general
232  YaspGrid<dim,Coordinates>* mg; // each grid level knows its multigrid
233  int overlapSize; // in mesh cells on this level
234  bool keepOverlap;
235 
237  int level_;
238  };
239 
241  typedef std::array<int, dim> iTupel;
242  typedef FieldVector<ctype, dim> fTupel;
243 
244  // communication tag used by multigrid
245  enum { tag = 17 };
246 #endif
247 
250  {
251  return _torus;
252  }
253 
255  int globalSize(int i) const
256  {
257  return levelSize(maxLevel(),i);
258  }
259 
261  iTupel globalSize() const
262  {
263  return levelSize(maxLevel());
264  }
265 
267  int levelSize(int l, int i) const
268  {
269  return _coarseSize[i] * (1 << l);
270  }
271 
273  iTupel levelSize(int l) const
274  {
275  iTupel s;
276  for (int i=0; i<dim; ++i)
277  s[i] = levelSize(l,i);
278  return s;
279  }
280 
282  bool isPeriodic(int i) const
283  {
284  return _periodic[i];
285  }
286 
287  bool getRefineOption() const
288  {
289  return keep_ovlp;
290  }
291 
294 
297  {
298  return YGridLevelIterator(_levels,0);
299  }
300 
302  YGridLevelIterator begin (int i) const
303  {
304  if (i<0 || i>maxLevel())
305  DUNE_THROW(GridError, "level not existing");
306  return YGridLevelIterator(_levels,i);
307  }
308 
311  {
312  return YGridLevelIterator(_levels,maxLevel()+1);
313  }
314 
315  // static method to create the default load balance strategy
316  static const YLoadBalanceDefault<dim>* defaultLoadbalancer()
317  {
318  static YLoadBalanceDefault<dim> lb;
319  return & lb;
320  }
321 
322  protected:
330  void makelevel (const Coordinates& coords, std::bitset<dim> periodic, iTupel o_interior, int overlap)
331  {
332  YGridLevel& g = _levels.back();
333  g.overlapSize = overlap;
334  g.mg = this;
335  g.level_ = maxLevel();
336  g.coords = coords;
337  g.keepOverlap = keep_ovlp;
338 
339  // set the inserting positions in the corresponding arrays of YGridLevelStructure
340  typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator overlapfront_it = g.overlapfront_data.begin();
341  typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator overlap_it = g.overlap_data.begin();
342  typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator interiorborder_it = g.interiorborder_data.begin();
343  typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator interior_it = g.interior_data.begin();
344 
345  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
346  send_overlapfront_overlapfront_it = g.send_overlapfront_overlapfront_data.begin();
347  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
348  recv_overlapfront_overlapfront_it = g.recv_overlapfront_overlapfront_data.begin();
349 
350  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
351  send_overlap_overlapfront_it = g.send_overlap_overlapfront_data.begin();
352  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
353  recv_overlapfront_overlap_it = g.recv_overlapfront_overlap_data.begin();
354 
355  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
356  send_interiorborder_interiorborder_it = g.send_interiorborder_interiorborder_data.begin();
357  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
358  recv_interiorborder_interiorborder_it = g.recv_interiorborder_interiorborder_data.begin();
359 
360  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
361  send_interiorborder_overlapfront_it = g.send_interiorborder_overlapfront_data.begin();
362  typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
363  recv_overlapfront_interiorborder_it = g.recv_overlapfront_interiorborder_data.begin();
364 
365  // have a null array for constructor calls around
366  std::array<int,dim> n;
367  std::fill(n.begin(), n.end(), 0);
368 
369  // determine origin of the grid with overlap and store whether an overlap area exists in direction i.
370  std::bitset<dim> ovlp_low(0ULL);
371  std::bitset<dim> ovlp_up(0ULL);
372 
373  iTupel o_overlap;
374  iTupel s_overlap;
375 
376  // determine at where we have overlap and how big the size of the overlap partition is
377  for (int i=0; i<dim; i++)
378  {
379  // the coordinate container has been contructed to hold the entire grid on
380  // this processor, including overlap. this is the element size.
381  s_overlap[i] = g.coords.size(i);
382 
383  //in the periodic case there is always overlap
384  if (periodic[i])
385  {
386  o_overlap[i] = o_interior[i]-overlap;
387  ovlp_low[i] = true;
388  ovlp_up[i] = true;
389  }
390  else
391  {
392  //check lower boundary
393  if (o_interior[i] - overlap < 0)
394  o_overlap[i] = 0;
395  else
396  {
397  o_overlap[i] = o_interior[i] - overlap;
398  ovlp_low[i] = true;
399  }
400 
401  //check upper boundary
402  if (o_overlap[i] + g.coords.size(i) < globalSize(i))
403  ovlp_up[i] = true;
404  }
405  }
406 
407  for (unsigned int codim = 0; codim < dim + 1; codim++)
408  {
409  // set the begin iterator for the corresponding ygrids
410  g.overlapfront[codim].setBegin(overlapfront_it);
411  g.overlap[codim].setBegin(overlap_it);
412  g.interiorborder[codim].setBegin(interiorborder_it);
413  g.interior[codim].setBegin(interior_it);
414  g.send_overlapfront_overlapfront[codim].setBegin(send_overlapfront_overlapfront_it);
415  g.recv_overlapfront_overlapfront[codim].setBegin(recv_overlapfront_overlapfront_it);
416  g.send_overlap_overlapfront[codim].setBegin(send_overlap_overlapfront_it);
417  g.recv_overlapfront_overlap[codim].setBegin(recv_overlapfront_overlap_it);
418  g.send_interiorborder_interiorborder[codim].setBegin(send_interiorborder_interiorborder_it);
419  g.recv_interiorborder_interiorborder[codim].setBegin(recv_interiorborder_interiorborder_it);
420  g.send_interiorborder_overlapfront[codim].setBegin(send_interiorborder_overlapfront_it);
421  g.recv_overlapfront_interiorborder[codim].setBegin(recv_overlapfront_interiorborder_it);
422 
423  // find all combinations of unit vectors that span entities of the given codimension
424  for (unsigned int index = 0; index < (1<<dim); index++)
425  {
426  // check whether the given shift is of our codimension
427  std::bitset<dim> r(index);
428  if (r.count() != dim-codim)
429  continue;
430 
431  // get an origin and a size array for subsequent modification
432  std::array<int,dim> origin(o_overlap);
433  std::array<int,dim> size(s_overlap);
434 
435  // build overlapfront
436  // we have to extend the element size by one in all directions without shift.
437  for (int i=0; i<dim; i++)
438  if (!r[i])
439  size[i]++;
440  *overlapfront_it = YGridComponent<Coordinates>(origin, r, &g.coords, size, n, size);
441 
442  // build overlap
443  for (int i=0; i<dim; i++)
444  {
445  if (!r[i])
446  {
447  if (ovlp_low[i])
448  {
449  origin[i]++;
450  size[i]--;
451  }
452  if (ovlp_up[i])
453  size[i]--;
454  }
455  }
456  *overlap_it = YGridComponent<Coordinates>(origin,size,*overlapfront_it);
457 
458  // build interiorborder
459  for (int i=0; i<dim; i++)
460  {
461  if (ovlp_low[i])
462  {
463  origin[i] += overlap;
464  size[i] -= overlap;
465  if (!r[i])
466  {
467  origin[i]--;
468  size[i]++;
469  }
470  }
471  if (ovlp_up[i])
472  {
473  size[i] -= overlap;
474  if (!r[i])
475  size[i]++;
476  }
477  }
478  *interiorborder_it = YGridComponent<Coordinates>(origin,size,*overlapfront_it);
479 
480  // build interior
481  for (int i=0; i<dim; i++)
482  {
483  if (!r[i])
484  {
485  if (ovlp_low[i])
486  {
487  origin[i]++;
488  size[i]--;
489  }
490  if (ovlp_up[i])
491  size[i]--;
492  }
493  }
494  *interior_it = YGridComponent<Coordinates>(origin, size, *overlapfront_it);
495 
496  intersections(*overlapfront_it,*overlapfront_it,*send_overlapfront_overlapfront_it, *recv_overlapfront_overlapfront_it);
497  intersections(*overlap_it,*overlapfront_it,*send_overlap_overlapfront_it, *recv_overlapfront_overlap_it);
498  intersections(*interiorborder_it,*interiorborder_it,*send_interiorborder_interiorborder_it,*recv_interiorborder_interiorborder_it);
499  intersections(*interiorborder_it,*overlapfront_it,*send_interiorborder_overlapfront_it,*recv_overlapfront_interiorborder_it);
500 
501  // advance all iterators pointing to the next insertion point
502  ++overlapfront_it;
503  ++overlap_it;
504  ++interiorborder_it;
505  ++interior_it;
506  ++send_overlapfront_overlapfront_it;
507  ++recv_overlapfront_overlapfront_it;
508  ++send_overlap_overlapfront_it;
509  ++recv_overlapfront_overlap_it;
510  ++send_interiorborder_interiorborder_it;
511  ++recv_interiorborder_interiorborder_it;
512  ++send_interiorborder_overlapfront_it;
513  ++recv_overlapfront_interiorborder_it;
514  }
515 
516  // set end iterators in the corresonding ygrids
517  g.overlapfront[codim].finalize(overlapfront_it);
518  g.overlap[codim].finalize(overlap_it);
519  g.interiorborder[codim].finalize(interiorborder_it);
520  g.interior[codim].finalize(interior_it);
521  g.send_overlapfront_overlapfront[codim].finalize(send_overlapfront_overlapfront_it,g.overlapfront[codim]);
522  g.recv_overlapfront_overlapfront[codim].finalize(recv_overlapfront_overlapfront_it,g.overlapfront[codim]);
523  g.send_overlap_overlapfront[codim].finalize(send_overlap_overlapfront_it,g.overlapfront[codim]);
524  g.recv_overlapfront_overlap[codim].finalize(recv_overlapfront_overlap_it,g.overlapfront[codim]);
525  g.send_interiorborder_interiorborder[codim].finalize(send_interiorborder_interiorborder_it,g.overlapfront[codim]);
526  g.recv_interiorborder_interiorborder[codim].finalize(recv_interiorborder_interiorborder_it,g.overlapfront[codim]);
527  g.send_interiorborder_overlapfront[codim].finalize(send_interiorborder_overlapfront_it,g.overlapfront[codim]);
528  g.recv_overlapfront_interiorborder[codim].finalize(recv_overlapfront_interiorborder_it,g.overlapfront[codim]);
529  }
530  }
531 
532 #ifndef DOXYGEN
541  struct mpifriendly_ygrid {
542  mpifriendly_ygrid ()
543  {
544  std::fill(origin.begin(), origin.end(), 0);
545  std::fill(size.begin(), size.end(), 0);
546  }
547  mpifriendly_ygrid (const YGridComponent<Coordinates>& grid)
548  : origin(grid.origin()), size(grid.size())
549  {}
550  iTupel origin;
551  iTupel size;
552  };
553 #endif
554 
564  std::deque<Intersection>& sendlist, std::deque<Intersection>& recvlist)
565  {
566  iTupel size = globalSize();
567 
568  // the exchange buffers
569  std::vector<YGridComponent<Coordinates> > send_recvgrid(_torus.neighbors());
570  std::vector<YGridComponent<Coordinates> > recv_recvgrid(_torus.neighbors());
571  std::vector<YGridComponent<Coordinates> > send_sendgrid(_torus.neighbors());
572  std::vector<YGridComponent<Coordinates> > recv_sendgrid(_torus.neighbors());
573 
574  // new exchange buffers to send simple struct without virtual functions
575  std::vector<mpifriendly_ygrid> mpifriendly_send_recvgrid(_torus.neighbors());
576  std::vector<mpifriendly_ygrid> mpifriendly_recv_recvgrid(_torus.neighbors());
577  std::vector<mpifriendly_ygrid> mpifriendly_send_sendgrid(_torus.neighbors());
578  std::vector<mpifriendly_ygrid> mpifriendly_recv_sendgrid(_torus.neighbors());
579 
580  // fill send buffers; iterate over all neighboring processes
581  // non-periodic case is handled automatically because intersection will be zero
582  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
583  {
584  // determine if we communicate with this neighbor (and what)
585  bool skip = false;
586  iTupel coord = _torus.coord(); // my coordinates
587  iTupel delta = i.delta(); // delta to neighbor
588  iTupel nb = coord; // the neighbor
589  for (int k=0; k<dim; k++) nb[k] += delta[k];
590  iTupel v; // grid movement
591  std::fill(v.begin(), v.end(), 0);
592 
593  for (int k=0; k<dim; k++)
594  {
595  if (nb[k]<0)
596  {
597  if (_periodic[k])
598  v[k] += size[k];
599  else
600  skip = true;
601  }
602  if (nb[k]>=_torus.dims(k))
603  {
604  if (_periodic[k])
605  v[k] -= size[k];
606  else
607  skip = true;
608  }
609  // neither might be true, then v=0
610  }
611 
612  // store moved grids in send buffers
613  if (!skip)
614  {
615  send_sendgrid[i.index()] = sendgrid.move(v);
616  send_recvgrid[i.index()] = recvgrid.move(v);
617  }
618  else
619  {
620  send_sendgrid[i.index()] = YGridComponent<Coordinates>();
621  send_recvgrid[i.index()] = YGridComponent<Coordinates>();
622  }
623  }
624 
625  // issue send requests for sendgrid being sent to all neighbors
626  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
627  {
628  mpifriendly_send_sendgrid[i.index()] = mpifriendly_ygrid(send_sendgrid[i.index()]);
629  _torus.send(i.rank(), &mpifriendly_send_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
630  }
631 
632  // issue recv requests for sendgrids of neighbors
633  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
634  _torus.recv(i.rank(), &mpifriendly_recv_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
635 
636  // exchange the sendgrids
637  _torus.exchange();
638 
639  // issue send requests for recvgrid being sent to all neighbors
640  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
641  {
642  mpifriendly_send_recvgrid[i.index()] = mpifriendly_ygrid(send_recvgrid[i.index()]);
643  _torus.send(i.rank(), &mpifriendly_send_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
644  }
645 
646  // issue recv requests for recvgrid of neighbors
647  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
648  _torus.recv(i.rank(), &mpifriendly_recv_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
649 
650  // exchange the recvgrid
651  _torus.exchange();
652 
653  // process receive buffers and compute intersections
654  for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
655  {
656  // what must be sent to this neighbor
657  Intersection send_intersection;
658  mpifriendly_ygrid yg = mpifriendly_recv_recvgrid[i.index()];
659  recv_recvgrid[i.index()] = YGridComponent<Coordinates>(yg.origin,yg.size);
660  send_intersection.grid = sendgrid.intersection(recv_recvgrid[i.index()]);
661  send_intersection.rank = i.rank();
662  send_intersection.distance = i.distance();
663  if (!send_intersection.grid.empty()) sendlist.push_front(send_intersection);
664 
665  Intersection recv_intersection;
666  yg = mpifriendly_recv_sendgrid[i.index()];
667  recv_sendgrid[i.index()] = YGridComponent<Coordinates>(yg.origin,yg.size);
668  recv_intersection.grid = recvgrid.intersection(recv_sendgrid[i.index()]);
669  recv_intersection.rank = i.rank();
670  recv_intersection.distance = i.distance();
671  if(!recv_intersection.grid.empty()) recvlist.push_back(recv_intersection);
672  }
673  }
674 
675  protected:
676 
677  typedef const YaspGrid<dim,Coordinates> GridImp;
678 
679  void init()
680  {
681  indexsets.push_back( std::make_shared< YaspIndexSet<const YaspGrid<dim, Coordinates>, false > >(*this,0) );
682  boundarysegmentssize();
683  }
684 
685  void boundarysegmentssize()
686  {
687  // sizes of local macro grid
688  std::array<int, dim> sides;
689  {
690  for (int i=0; i<dim; i++)
691  {
692  sides[i] =
693  ((begin()->overlap[0].dataBegin()->origin(i) == 0)+
694  (begin()->overlap[0].dataBegin()->origin(i) + begin()->overlap[0].dataBegin()->size(i)
695  == levelSize(0,i)));
696  }
697  }
698  nBSegments = 0;
699  for (int k=0; k<dim; k++)
700  {
701  int offset = 1;
702  for (int l=0; l<dim; l++)
703  {
704  if (l==k) continue;
705  offset *= begin()->overlap[0].dataBegin()->size(l);
706  }
707  nBSegments += sides[k]*offset;
708  }
709  }
710 
711  public:
712 
713  // define the persistent index type
714  typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim> PersistentIndexType;
715 
717  typedef YaspGridFamily<dim, Coordinates> GridFamily;
718  // the Traits
720 
721  // need for friend declarations in entity
725 
733  YaspGrid (const Coordinates& coordinates,
734  std::bitset<dim> periodic = std::bitset<dim>(0ULL),
735  int overlap = 1,
736  CollectiveCommunicationType comm = CollectiveCommunicationType(),
737  const YLoadBalance<dim>* lb = defaultLoadbalancer())
738  : ccobj(comm)
739  , leafIndexSet_(*this)
740  , _periodic(periodic)
741  , _overlap(overlap)
742  , keep_ovlp(true)
743  , adaptRefCount(0)
744  , adaptActive(false)
745  {
746  _levels.resize(1);
747 
748  // Number of elements per coordinate direction on the coarsest level
749  for (std::size_t i=0; i<dim; i++)
750  _coarseSize[i] = coordinates.size(i);
751 
752  // Construct the communication torus
753  _torus = decltype(_torus)(comm,tag,_coarseSize,lb);
754 
755  iTupel o;
756  std::fill(o.begin(), o.end(), 0);
757  iTupel o_interior(o);
758  iTupel s_interior(_coarseSize);
759 
760  _torus.partition(_torus.rank(),o,_coarseSize,o_interior,s_interior);
761 
762  // Set domain size
763  if (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value
764  || std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >::value)
765  {
766  for (std::size_t i=0; i<dim; i++)
767  _L[i] = coordinates.size(i) * coordinates.meshsize(i,0);
768  }
769  if (std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >::value)
770  {
771  //determine sizes of vector to correctly construct torus structure and store for later size requests
772  for (int i=0; i<dim; i++)
773  _L[i] = coordinates.coordinate(i,_coarseSize[i]) - coordinates.coordinate(i,0);
774  }
775 
776 #if HAVE_MPI
777  // TODO: Settle on a single value for all coordinate types
778  int mysteryFactor = (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value) ? 1 : 2;
779 
780  // check whether the grid is large enough to be overlapping
781  for (int i=0; i<dim; i++)
782  {
783  // find out whether the grid is too small to
784  int toosmall = (s_interior[i] <= mysteryFactor * overlap) && // interior is very small
785  (periodic[i] || (s_interior[i] != _coarseSize[i])); // there is an overlap in that direction
786  // communicate the result to all those processes to have all processors error out if one process failed.
787  int global = 0;
788  MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
789  if (global)
790  DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
791  " Note that this also holds for DOFs on subdomain boundaries."
792  " Increase grid elements or decrease overlap accordingly.");
793  }
794 #endif // #if HAVE_MPI
795 
796  if (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value
797  || std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >::value)
798  {
799  iTupel s_overlap(s_interior);
800  for (int i=0; i<dim; i++)
801  {
802  if ((o_interior[i] - overlap > 0) || (periodic[i]))
803  s_overlap[i] += overlap;
804  if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
805  s_overlap[i] += overlap;
806  }
807 
808  FieldVector<ctype,dim> upperRightWithOverlap;
809  for (int i=0; i<dim; i++)
810  upperRightWithOverlap[i] = coordinates.coordinate(i,0) + coordinates.meshsize(i,0) * s_overlap[i];
811 
812  Hybrid::ifElse(std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >{}, [&](auto id)
813  {
814  // New coordinate object that additionally contains the overlap elements
815  EquidistantCoordinates<ctype,dim> coordinatesWithOverlap(upperRightWithOverlap,s_overlap);
816 
817  // add level (the this-> is needed to make g++-6 happy)
818  this->makelevel(id(coordinatesWithOverlap),periodic,o_interior,overlap);
819  });
820 
821  Hybrid::ifElse(std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >{}, [&](auto id)
822  {
824  for (int i=0; i<dim; i++)
825  lowerleft[i] = id(coordinates).origin(i);
826 
827  // New coordinate object that additionally contains the overlap elements
828  EquidistantOffsetCoordinates<ctype,dim> coordinatesWithOverlap(lowerleft,upperRightWithOverlap,s_overlap);
829 
830  // add level (the this-> is needed to make g++-6 happy)
831  this->makelevel(id(coordinatesWithOverlap),periodic,o_interior,overlap);
832  });
833  }
834 
835  Hybrid::ifElse(std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >{}, [&](auto id)
836  {
837  std::array<std::vector<ctype>,dim> newCoords;
838  std::array<int, dim> offset(o_interior);
839 
840  // find the relevant part of the coords vector for this processor and copy it to newCoords
841  for (int i=0; i<dim; ++i)
842  {
843  //define the coordinate range to be used
844  std::size_t begin = o_interior[i];
845  std::size_t end = begin + s_interior[i] + 1;
846 
847  // check whether we are not at the physical boundary. In that case overlap is a simple
848  // extension of the coordinate range to be used
849  if (o_interior[i] - overlap > 0)
850  {
851  begin = begin - overlap;
852  offset[i] -= overlap;
853  }
854  if (o_interior[i] + s_interior[i] + overlap < _coarseSize[i])
855  end = end + overlap;
856 
857  //copy the selected part in the new coord vector
858  newCoords[i].resize(end-begin);
859  auto newCoordsIt = newCoords[i].begin();
860  for (std::size_t j=begin; j<end; j++)
861  {
862  *newCoordsIt = coordinates.coordinate(i, j);
863  newCoordsIt++;
864  }
865 
866  // Check whether we are at the physical boundary and have a periodic grid.
867  // In this case the coordinate vector has to be tweaked manually.
868  if ((periodic[i]) && (o_interior[i] + s_interior[i] + overlap >= _coarseSize[i]))
869  {
870  // we need to add the first <overlap> cells to the end of newcoords
871  for (int j=0; j<overlap; ++j)
872  newCoords[i].push_back(newCoords[i].back() - coordinates.coordinate(i,j) + coordinates.coordinate(i,j+1));
873  }
874 
875  if ((periodic[i]) && (o_interior[i] - overlap <= 0))
876  {
877  offset[i] -= overlap;
878 
879  // we need to add the last <overlap> cells to the begin of newcoords
880  std::size_t reverseCounter = coordinates.size(i);
881  for (int j=0; j<overlap; ++j, --reverseCounter)
882  newCoords[i].insert(newCoords[i].begin(), newCoords[i].front()
883  - coordinates.coordinate(i,reverseCounter) + coordinates.coordinate(i,reverseCounter-1));
884  }
885  }
886 
887  TensorProductCoordinates<ctype,dim> coordinatesWithOverlap(newCoords, offset);
888 
889  // add level (the this-> is needed to make g++-6 happy)
890  this->makelevel(id(coordinatesWithOverlap),periodic,o_interior,overlap);
891  });
892 
893  init();
894  }
895 
905  std::array<int, dim> s,
906  std::bitset<dim> periodic = std::bitset<dim>(0ULL),
907  int overlap = 1,
908  CollectiveCommunicationType comm = CollectiveCommunicationType(),
909  const YLoadBalance<dim>* lb = defaultLoadbalancer())
910  : ccobj(comm), _torus(comm,tag,s,lb), leafIndexSet_(*this),
911  _L(L), _periodic(periodic), _coarseSize(s), _overlap(overlap),
912  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
913  {
914  // check whether YaspGrid has been given the correct template parameter
915  static_assert(std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value,
916  "YaspGrid coordinate container template parameter and given constructor values do not match!");
917 
918  _levels.resize(1);
919 
920  iTupel o;
921  std::fill(o.begin(), o.end(), 0);
922  iTupel o_interior(o);
923  iTupel s_interior(s);
924 
925  _torus.partition(_torus.rank(),o,s,o_interior,s_interior);
926 
927 #if HAVE_MPI
928  // check whether the grid is large enough to be overlapping
929  for (int i=0; i<dim; i++)
930  {
931  // find out whether the grid is too small to
932  int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
933  (periodic[i] || (s_interior[i] != s[i])); // there is an overlap in that direction
934  // communicate the result to all those processes to have all processors error out if one process failed.
935  int global = 0;
936  MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
937  if (global)
938  DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
939  " Note that this also holds for DOFs on subdomain boundaries."
940  " Increase grid elements or decrease overlap accordingly.");
941  }
942 #endif // #if HAVE_MPI
943 
944  iTupel s_overlap(s_interior);
945  for (int i=0; i<dim; i++)
946  {
947  if ((o_interior[i] - overlap > 0) || (periodic[i]))
948  s_overlap[i] += overlap;
949  if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
950  s_overlap[i] += overlap;
951  }
952 
953  FieldVector<ctype,dim> upperRightWithOverlap;
954 
955  for (int i=0; i<dim; i++)
956  upperRightWithOverlap[i] = (L[i] / s[i]) * s_overlap[i];
957 
958  // New coordinate object that additionally contains the overlap elements
959  EquidistantCoordinates<ctype,dim> cc(upperRightWithOverlap,s_overlap);
960 
961  // add level
962  makelevel(cc,periodic,o_interior,overlap);
963 
964  init();
965  }
966 
978  std::array<int, dim> s,
979  std::bitset<dim> periodic = std::bitset<dim>(0ULL),
980  int overlap = 1,
981  CollectiveCommunicationType comm = CollectiveCommunicationType(),
982  const YLoadBalance<dim>* lb = defaultLoadbalancer())
983  : ccobj(comm), _torus(comm,tag,s,lb), leafIndexSet_(*this),
984  _L(upperright - lowerleft),
985  _periodic(periodic), _coarseSize(s), _overlap(overlap),
986  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
987  {
988  // check whether YaspGrid has been given the correct template parameter
989  static_assert(std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >::value,
990  "YaspGrid coordinate container template parameter and given constructor values do not match!");
991 
992  _levels.resize(1);
993 
994  iTupel o;
995  std::fill(o.begin(), o.end(), 0);
996  iTupel o_interior(o);
997  iTupel s_interior(s);
998 
999  _torus.partition(_torus.rank(),o,s,o_interior,s_interior);
1000 
1001 #if HAVE_MPI
1002  // check whether the grid is large enough to be overlapping
1003  for (int i=0; i<dim; i++)
1004  {
1005  // find out whether the grid is too small to
1006  int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
1007  (periodic[i] || (s_interior[i] != s[i])); // there is an overlap in that direction
1008  // communicate the result to all those processes to have all processors error out if one process failed.
1009  int global = 0;
1010  MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
1011  if (global)
1012  DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
1013  " Note that this also holds for DOFs on subdomain boundaries."
1014  " Increase grid elements or decrease overlap accordingly.");
1015  }
1016 #endif // #if HAVE_MPI
1017 
1018  iTupel s_overlap(s_interior);
1019  for (int i=0; i<dim; i++)
1020  {
1021  if ((o_interior[i] - overlap > 0) || (periodic[i]))
1022  s_overlap[i] += overlap;
1023  if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
1024  s_overlap[i] += overlap;
1025  }
1026 
1027  FieldVector<ctype,dim> upperRightWithOverlap;
1028  for (int i=0; i<dim; i++)
1029  upperRightWithOverlap[i] = lowerleft[i]
1030  + s_overlap[i] * (upperright[i]-lowerleft[i]) / s[i];
1031 
1032  EquidistantOffsetCoordinates<ctype,dim> cc(lowerleft,upperRightWithOverlap,s_overlap);
1033 
1034  // add level
1035  makelevel(cc,periodic,o_interior,overlap);
1036 
1037  init();
1038  }
1039 
1047  YaspGrid (std::array<std::vector<ctype>, dim> coords,
1048  std::bitset<dim> periodic = std::bitset<dim>(0ULL),
1049  int overlap = 1,
1050  CollectiveCommunicationType comm = CollectiveCommunicationType(),
1051  const YLoadBalance<dim>* lb = defaultLoadbalancer())
1052  : ccobj(comm), _torus(comm,tag,Dune::Yasp::sizeArray<dim>(coords),lb),
1053  leafIndexSet_(*this), _periodic(periodic), _overlap(overlap),
1054  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
1055  {
1056  if (!Dune::Yasp::checkIfMonotonous(coords))
1057  DUNE_THROW(Dune::GridError,"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1058 
1059  // check whether YaspGrid has been given the correct template parameter
1060  static_assert(std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >::value,
1061  "YaspGrid coordinate container template parameter and given constructor values do not match!");
1062 
1063  _levels.resize(1);
1064 
1065  //determine sizes of vector to correctly construct torus structure and store for later size requests
1066  for (int i=0; i<dim; i++) {
1067  _coarseSize[i] = coords[i].size() - 1;
1068  _L[i] = coords[i][_coarseSize[i]] - coords[i][0];
1069  }
1070 
1071  iTupel o;
1072  std::fill(o.begin(), o.end(), 0);
1073  iTupel o_interior(o);
1074  iTupel s_interior(_coarseSize);
1075 
1076  _torus.partition(_torus.rank(),o,_coarseSize,o_interior,s_interior);
1077 
1078 #if HAVE_MPI
1079  // check whether the grid is large enough to be overlapping
1080  for (int i=0; i<dim; i++)
1081  {
1082  // find out whether the grid is too small to
1083  int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
1084  (periodic[i] || (s_interior[i] != _coarseSize[i])); // there is an overlap in that direction
1085  // communicate the result to all those processes to have all processors error out if one process failed.
1086  int global = 0;
1087  MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
1088  if (global)
1089  DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
1090  " Note that this also holds for DOFs on subdomain boundaries."
1091  " Increase grid elements or decrease overlap accordingly.");
1092  }
1093 #endif // #if HAVE_MPI
1094 
1095 
1096  std::array<std::vector<ctype>,dim> newcoords;
1097  std::array<int, dim> offset(o_interior);
1098 
1099  // find the relevant part of the coords vector for this processor and copy it to newcoords
1100  for (int i=0; i<dim; ++i)
1101  {
1102  //define iterators on coords that specify the coordinate range to be used
1103  typename std::vector<ctype>::iterator begin = coords[i].begin() + o_interior[i];
1104  typename std::vector<ctype>::iterator end = begin + s_interior[i] + 1;
1105 
1106  // check whether we are not at the physical boundary. In that case overlap is a simple
1107  // extension of the coordinate range to be used
1108  if (o_interior[i] - overlap > 0)
1109  {
1110  begin = begin - overlap;
1111  offset[i] -= overlap;
1112  }
1113  if (o_interior[i] + s_interior[i] + overlap < _coarseSize[i])
1114  end = end + overlap;
1115 
1116  //copy the selected part in the new coord vector
1117  newcoords[i].resize(end-begin);
1118  std::copy(begin, end, newcoords[i].begin());
1119 
1120  // check whether we are at the physical boundary and a have a periodic grid.
1121  // In this case the coordinate vector has to be tweaked manually.
1122  if ((periodic[i]) && (o_interior[i] + s_interior[i] + overlap >= _coarseSize[i]))
1123  {
1124  // we need to add the first <overlap> cells to the end of newcoords
1125  typename std::vector<ctype>::iterator it = coords[i].begin();
1126  for (int j=0; j<overlap; ++j)
1127  newcoords[i].push_back(newcoords[i].back() - *it + *(++it));
1128  }
1129 
1130  if ((periodic[i]) && (o_interior[i] - overlap <= 0))
1131  {
1132  offset[i] -= overlap;
1133 
1134  // we need to add the last <overlap> cells to the begin of newcoords
1135  typename std::vector<ctype>::iterator it = coords[i].end() - 1;
1136  for (int j=0; j<overlap; ++j)
1137  newcoords[i].insert(newcoords[i].begin(), newcoords[i].front() - *it + *(--it));
1138  }
1139  }
1140 
1141  TensorProductCoordinates<ctype,dim> cc(newcoords, offset);
1142 
1143  // add level
1144  makelevel(cc,periodic,o_interior,overlap);
1145  init();
1146  }
1147 
1148  private:
1149 
1164  YaspGrid (std::array<std::vector<ctype>, dim> coords,
1165  std::bitset<dim> periodic,
1166  int overlap,
1167  CollectiveCommunicationType comm,
1168  std::array<int,dim> coarseSize,
1169  const YLoadBalance<dim>* lb = defaultLoadbalancer())
1170  : ccobj(comm), _torus(comm,tag,coarseSize,lb), leafIndexSet_(*this),
1171  _periodic(periodic), _coarseSize(coarseSize), _overlap(overlap),
1172  keep_ovlp(true), adaptRefCount(0), adaptActive(false)
1173  {
1174  // check whether YaspGrid has been given the correct template parameter
1175  static_assert(std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >::value,
1176  "YaspGrid coordinate container template parameter and given constructor values do not match!");
1177 
1178  if (!Dune::Yasp::checkIfMonotonous(coords))
1179  DUNE_THROW(Dune::GridError,"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1180 
1181  for (int i=0; i<dim; i++)
1182  _L[i] = coords[i][coords[i].size() - 1] - coords[i][0];
1183 
1184  _levels.resize(1);
1185 
1186  std::array<int,dim> o;
1187  std::fill(o.begin(), o.end(), 0);
1188  std::array<int,dim> o_interior(o);
1189  std::array<int,dim> s_interior(coarseSize);
1190 
1191  _torus.partition(_torus.rank(),o,coarseSize,o_interior,s_interior);
1192 
1193  // get offset by modifying o_interior according to overlap
1194  std::array<int,dim> offset(o_interior);
1195  for (int i=0; i<dim; i++)
1196  if ((periodic[i]) || (o_interior[i] > 0))
1197  offset[i] -= overlap;
1198 
1199  TensorProductCoordinates<ctype,dim> cc(coords, offset);
1200 
1201  // add level
1202  makelevel(cc,periodic,o_interior,overlap);
1203 
1204  init();
1205  }
1206 
1207  // the backup restore facility needs to be able to use above constructor
1208  friend struct BackupRestoreFacility<YaspGrid<dim,Coordinates> >;
1209 
1210  // do not copy this class
1211  YaspGrid(const YaspGrid&);
1212 
1213  public:
1214 
1218  int maxLevel() const
1219  {
1220  return _levels.size()-1;
1221  }
1222 
1224  void globalRefine (int refCount)
1225  {
1226  if (refCount < -maxLevel())
1227  DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
1228  "Coarsening " << -refCount << " levels requested!");
1229 
1230  // If refCount is negative then coarsen the grid
1231  for (int k=refCount; k<0; k++)
1232  {
1233  // create an empty grid level
1234  YGridLevel empty;
1235  _levels.back() = empty;
1236  // reduce maxlevel
1237  _levels.pop_back();
1238 
1239  indexsets.pop_back();
1240  }
1241 
1242  // If refCount is positive refine the grid
1243  for (int k=0; k<refCount; k++)
1244  {
1245  // access to coarser grid level
1246  YGridLevel& cg = _levels[maxLevel()];
1247 
1248  std::bitset<dim> ovlp_low(0ULL), ovlp_up(0ULL);
1249  for (int i=0; i<dim; i++)
1250  {
1251  if (cg.overlap[0].dataBegin()->origin(i) > 0 || _periodic[i])
1252  ovlp_low[i] = true;
1253  if (cg.overlap[0].dataBegin()->max(i) + 1 < globalSize(i) || _periodic[i])
1254  ovlp_up[i] = true;
1255  }
1256 
1257  Coordinates newcont(cg.coords.refine(ovlp_low, ovlp_up, cg.overlapSize, keep_ovlp));
1258 
1259  int overlap = (keep_ovlp) ? 2*cg.overlapSize : cg.overlapSize;
1260 
1261  //determine new origin
1262  iTupel o_interior;
1263  for (int i=0; i<dim; i++)
1264  o_interior[i] = 2*cg.interior[0].dataBegin()->origin(i);
1265 
1266  // add level
1267  _levels.resize(_levels.size() + 1);
1268  makelevel(newcont,_periodic,o_interior,overlap);
1269 
1270  indexsets.push_back( std::make_shared<YaspIndexSet<const YaspGrid<dim,Coordinates>, false > >(*this,maxLevel()) );
1271  }
1272  }
1273 
1278  void refineOptions (bool keepPhysicalOverlap)
1279  {
1280  keep_ovlp = keepPhysicalOverlap;
1281  }
1282 
1294  bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
1295  {
1296  assert(adaptActive == false);
1297  if (e.level() != maxLevel()) return false;
1298  adaptRefCount = std::max(adaptRefCount, refCount);
1299  return true;
1300  }
1301 
1308  int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
1309  {
1310  return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
1311  }
1312 
1314  bool adapt ()
1315  {
1316  globalRefine(adaptRefCount);
1317  return (adaptRefCount > 0);
1318  }
1319 
1321  bool preAdapt ()
1322  {
1323  adaptActive = true;
1324  adaptRefCount = comm().max(adaptRefCount);
1325  return (adaptRefCount < 0);
1326  }
1327 
1329  void postAdapt()
1330  {
1331  adaptActive = false;
1332  adaptRefCount = 0;
1333  }
1334 
1336  template<int cd, PartitionIteratorType pitype>
1337  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
1338  {
1339  return levelbegin<cd,pitype>(level);
1340  }
1341 
1343  template<int cd, PartitionIteratorType pitype>
1344  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
1345  {
1346  return levelend<cd,pitype>(level);
1347  }
1348 
1350  template<int cd>
1351  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
1352  {
1353  return levelbegin<cd,All_Partition>(level);
1354  }
1355 
1357  template<int cd>
1358  typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
1359  {
1360  return levelend<cd,All_Partition>(level);
1361  }
1362 
1364  template<int cd, PartitionIteratorType pitype>
1365  typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const
1366  {
1367  return levelbegin<cd,pitype>(maxLevel());
1368  }
1369 
1371  template<int cd, PartitionIteratorType pitype>
1372  typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const
1373  {
1374  return levelend<cd,pitype>(maxLevel());
1375  }
1376 
1378  template<int cd>
1379  typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
1380  {
1381  return levelbegin<cd,All_Partition>(maxLevel());
1382  }
1383 
1385  template<int cd>
1386  typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
1387  {
1388  return levelend<cd,All_Partition>(maxLevel());
1389  }
1390 
1391  // \brief obtain Entity from EntitySeed. */
1392  template <typename Seed>
1393  typename Traits::template Codim<Seed::codimension>::Entity
1394  entity(const Seed& seed) const
1395  {
1396  const int codim = Seed::codimension;
1397  YGridLevelIterator g = begin(seed.impl().level());
1398 
1399  typedef typename Traits::template Codim<Seed::codimension>::Entity Entity;
1400  typedef YaspEntity<codim,dim,const YaspGrid> EntityImp;
1401  typedef typename YGrid::Iterator YIterator;
1402 
1403  return Entity(EntityImp(g,YIterator(g->overlapfront[codim],seed.impl().coord(),seed.impl().offset())));
1404  }
1405 
1407  int overlapSize (int level, int codim) const
1408  {
1409  YGridLevelIterator g = begin(level);
1410  return g->overlapSize;
1411  }
1412 
1414  int overlapSize (int codim) const
1415  {
1417  return g->overlapSize;
1418  }
1419 
1421  int ghostSize (int level, int codim) const
1422  {
1423  return 0;
1424  }
1425 
1427  int ghostSize (int codim) const
1428  {
1429  return 0;
1430  }
1431 
1433  int size (int level, int codim) const
1434  {
1435  YGridLevelIterator g = begin(level);
1436 
1437  // sum over all components of the codimension
1438  int count = 0;
1439  typedef typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator DAI;
1440  for (DAI it = g->overlapfront[codim].dataBegin(); it != g->overlapfront[codim].dataEnd(); ++it)
1441  count += it->totalsize();
1442 
1443  return count;
1444  }
1445 
1447  int size (int codim) const
1448  {
1449  return size(maxLevel(),codim);
1450  }
1451 
1453  int size (int level, GeometryType type) const
1454  {
1455  return (type.isCube()) ? size(level,dim-type.dim()) : 0;
1456  }
1457 
1459  int size (GeometryType type) const
1460  {
1461  return size(maxLevel(),type);
1462  }
1463 
1465  size_t numBoundarySegments () const
1466  {
1467  return nBSegments;
1468  }
1469 
1472  return _L;
1473  }
1474 
1479  template<class DataHandleImp, class DataType>
1481  {
1482  YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
1483  }
1484 
1489  template<class DataHandleImp, class DataType>
1491  {
1492  YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
1493  }
1494 
1499  template<class DataHandle, int codim>
1500  void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
1501  {
1502  // check input
1503  if (!data.contains(dim,codim)) return; // should have been checked outside
1504 
1505  // data types
1506  typedef typename DataHandle::DataType DataType;
1507 
1508  // access to grid level
1509  YGridLevelIterator g = begin(level);
1510 
1511  // find send/recv lists or throw error
1512  const YGridList<Coordinates>* sendlist = 0;
1513  const YGridList<Coordinates>* recvlist = 0;
1514 
1516  {
1517  sendlist = &g->send_interiorborder_interiorborder[codim];
1518  recvlist = &g->recv_interiorborder_interiorborder[codim];
1519  }
1520  if (iftype==InteriorBorder_All_Interface)
1521  {
1522  sendlist = &g->send_interiorborder_overlapfront[codim];
1523  recvlist = &g->recv_overlapfront_interiorborder[codim];
1524  }
1526  {
1527  sendlist = &g->send_overlap_overlapfront[codim];
1528  recvlist = &g->recv_overlapfront_overlap[codim];
1529  }
1530  if (iftype==All_All_Interface)
1531  {
1532  sendlist = &g->send_overlapfront_overlapfront[codim];
1533  recvlist = &g->recv_overlapfront_overlapfront[codim];
1534  }
1535 
1536  // change communication direction?
1537  if (dir==BackwardCommunication)
1538  std::swap(sendlist,recvlist);
1539 
1540  int cnt;
1541 
1542  // Size computation (requires communication if variable size)
1543  std::vector<int> send_size(sendlist->size(),-1); // map rank to total number of objects (of type DataType) to be sent
1544  std::vector<int> recv_size(recvlist->size(),-1); // map rank to total number of objects (of type DataType) to be recvd
1545  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
1546  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
1547 
1548  // define type to iterate over send and recv lists
1549  typedef typename YGridList<Coordinates>::Iterator ListIt;
1550 
1551  if (data.fixedSize(dim,codim))
1552  {
1553  // fixed size: just take a dummy entity, size can be computed without communication
1554  cnt=0;
1555  for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1556  {
1557  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1559  send_size[cnt] = is->grid.totalsize() * data.size(*it);
1560  cnt++;
1561  }
1562  cnt=0;
1563  for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1564  {
1565  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1567  recv_size[cnt] = is->grid.totalsize() * data.size(*it);
1568  cnt++;
1569  }
1570  }
1571  else
1572  {
1573  // variable size case: sender side determines the size
1574  cnt=0;
1575  for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1576  {
1577  // allocate send buffer for sizes per entitiy
1578  size_t *buf = new size_t[is->grid.totalsize()];
1579  send_sizes[cnt] = buf;
1580 
1581  // loop over entities and ask for size
1582  int i=0; size_t n=0;
1583  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1585  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1586  itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1587  for ( ; it!=itend; ++it)
1588  {
1589  buf[i] = data.size(*it);
1590  n += buf[i];
1591  i++;
1592  }
1593 
1594  // now we know the size for this rank
1595  send_size[cnt] = n;
1596 
1597  // hand over send request to torus class
1598  torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1599  cnt++;
1600  }
1601 
1602  // allocate recv buffers for sizes and store receive request
1603  cnt=0;
1604  for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1605  {
1606  // allocate recv buffer
1607  size_t *buf = new size_t[is->grid.totalsize()];
1608  recv_sizes[cnt] = buf;
1609 
1610  // hand over recv request to torus class
1611  torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1612  cnt++;
1613  }
1614 
1615  // exchange all size buffers now
1616  torus().exchange();
1617 
1618  // release send size buffers
1619  cnt=0;
1620  for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1621  {
1622  delete[] send_sizes[cnt];
1623  send_sizes[cnt] = 0;
1624  cnt++;
1625  }
1626 
1627  // process receive size buffers
1628  cnt=0;
1629  for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1630  {
1631  // get recv buffer
1632  size_t *buf = recv_sizes[cnt];
1633 
1634  // compute total size
1635  size_t n=0;
1636  for (int i=0; i<is->grid.totalsize(); ++i)
1637  n += buf[i];
1638 
1639  // ... and store it
1640  recv_size[cnt] = n;
1641  ++cnt;
1642  }
1643  }
1644 
1645 
1646  // allocate & fill the send buffers & store send request
1647  std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
1648  cnt=0;
1649  for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1650  {
1651  // allocate send buffer
1652  DataType *buf = new DataType[send_size[cnt]];
1653 
1654  // remember send buffer
1655  sends[cnt] = buf;
1656 
1657  // make a message buffer
1658  MessageBuffer<DataType> mb(buf);
1659 
1660  // fill send buffer; iterate over cells in intersection
1661  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1663  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1664  itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1665  for ( ; it!=itend; ++it)
1666  data.gather(mb,*it);
1667 
1668  // hand over send request to torus class
1669  torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
1670  cnt++;
1671  }
1672 
1673  // allocate recv buffers and store receive request
1674  std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
1675  cnt=0;
1676  for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1677  {
1678  // allocate recv buffer
1679  DataType *buf = new DataType[recv_size[cnt]];
1680 
1681  // remember recv buffer
1682  recvs[cnt] = buf;
1683 
1684  // hand over recv request to torus class
1685  torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
1686  cnt++;
1687  }
1688 
1689  // exchange all buffers now
1690  torus().exchange();
1691 
1692  // release send buffers
1693  cnt=0;
1694  for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1695  {
1696  delete[] sends[cnt];
1697  sends[cnt] = 0;
1698  cnt++;
1699  }
1700 
1701  // process receive buffers and delete them
1702  cnt=0;
1703  for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1704  {
1705  // get recv buffer
1706  DataType *buf = recvs[cnt];
1707 
1708  // make a message buffer
1709  MessageBuffer<DataType> mb(buf);
1710 
1711  // copy data from receive buffer; iterate over cells in intersection
1712  if (data.fixedSize(dim,codim))
1713  {
1714  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1716  size_t n=data.size(*it);
1717  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1718  itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1719  for ( ; it!=itend; ++it)
1720  data.scatter(mb,*it,n);
1721  }
1722  else
1723  {
1724  int i=0;
1725  size_t *sbuf = recv_sizes[cnt];
1726  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1728  typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1729  itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1730  for ( ; it!=itend; ++it)
1731  data.scatter(mb,*it,sbuf[i++]);
1732  delete[] sbuf;
1733  }
1734 
1735  // delete buffer
1736  delete[] buf; // hier krachts !
1737  cnt++;
1738  }
1739  }
1740 
1741  // The new index sets from DDM 11.07.2005
1742  const typename Traits::GlobalIdSet& globalIdSet() const
1743  {
1744  return theglobalidset;
1745  }
1746 
1747  const typename Traits::LocalIdSet& localIdSet() const
1748  {
1749  return theglobalidset;
1750  }
1751 
1752  const typename Traits::LevelIndexSet& levelIndexSet(int level) const
1753  {
1754  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1755  return *(indexsets[level]);
1756  }
1757 
1758  const typename Traits::LeafIndexSet& leafIndexSet() const
1759  {
1760  return leafIndexSet_;
1761  }
1762 
1765  const CollectiveCommunicationType& comm () const
1766  {
1767  return ccobj;
1768  }
1769 
1770  private:
1771 
1772  // number of boundary segments of the level 0 grid
1773  int nBSegments;
1774 
1775  // Index classes need access to the real entity
1776  friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim, Coordinates>, true >;
1777  friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim, Coordinates>, false >;
1778  friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim, Coordinates> >;
1779  friend class Dune::YaspPersistentContainerIndex<const Dune::YaspGrid<dim, Coordinates> >;
1780 
1781  friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim, Coordinates> >;
1782  friend class Dune::YaspIntersection<const Dune::YaspGrid<dim, Coordinates> >;
1783  friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim, Coordinates> >;
1784 
1785  template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1786  friend class Entity;
1787 
1788  template<class DT>
1789  class MessageBuffer {
1790  public:
1791  // Constructor
1792  MessageBuffer (DT *p)
1793  {
1794  a=p;
1795  i=0;
1796  j=0;
1797  }
1798 
1799  // write data to message buffer, acts like a stream !
1800  template<class Y>
1801  void write (const Y& data)
1802  {
1803  static_assert(( std::is_same<DT,Y>::value ), "DataType mismatch");
1804  a[i++] = data;
1805  }
1806 
1807  // read data from message buffer, acts like a stream !
1808  template<class Y>
1809  void read (Y& data) const
1810  {
1811  static_assert(( std::is_same<DT,Y>::value ), "DataType mismatch");
1812  data = a[j++];
1813  }
1814 
1815  private:
1816  DT *a;
1817  int i;
1818  mutable int j;
1819  };
1820 
1822  template<int cd, PartitionIteratorType pitype>
1823  YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
1824  {
1825  YGridLevelIterator g = begin(level);
1826  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1827 
1828  if (pitype==Interior_Partition)
1829  return YaspLevelIterator<cd,pitype,GridImp>(g,g->interior[cd].begin());
1830  if (pitype==InteriorBorder_Partition)
1831  return YaspLevelIterator<cd,pitype,GridImp>(g,g->interiorborder[cd].begin());
1832  if (pitype==Overlap_Partition)
1833  return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlap[cd].begin());
1834  if (pitype<=All_Partition)
1835  return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlapfront[cd].begin());
1836  if (pitype==Ghost_Partition)
1837  return levelend <cd, pitype> (level);
1838 
1839  DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1840  }
1841 
1843  template<int cd, PartitionIteratorType pitype>
1844  YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
1845  {
1846  YGridLevelIterator g = begin(level);
1847  if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1848 
1849  if (pitype==Interior_Partition)
1850  return YaspLevelIterator<cd,pitype,GridImp>(g,g->interior[cd].end());
1851  if (pitype==InteriorBorder_Partition)
1852  return YaspLevelIterator<cd,pitype,GridImp>(g,g->interiorborder[cd].end());
1853  if (pitype==Overlap_Partition)
1854  return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlap[cd].end());
1855  if (pitype<=All_Partition || pitype == Ghost_Partition)
1856  return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlapfront[cd].end());
1857 
1858  DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1859  }
1860 
1861  CollectiveCommunicationType ccobj;
1862 
1863  Torus<CollectiveCommunicationType,dim> _torus;
1864 
1865  std::vector< std::shared_ptr< YaspIndexSet<const YaspGrid<dim,Coordinates>, false > > > indexsets;
1866  YaspIndexSet<const YaspGrid<dim,Coordinates>, true> leafIndexSet_;
1867  YaspGlobalIdSet<const YaspGrid<dim,Coordinates> > theglobalidset;
1868 
1870  iTupel _s;
1871  std::bitset<dim> _periodic;
1872  iTupel _coarseSize;
1873  ReservedVector<YGridLevel,32> _levels;
1874  int _overlap;
1875  bool keep_ovlp;
1876  int adaptRefCount;
1877  bool adaptActive;
1878  };
1879 
1881 
1882  template <int d, class CC>
1883  std::ostream& operator<< (std::ostream& s, const YaspGrid<d,CC>& grid)
1884  {
1885  int rank = grid.torus().rank();
1886 
1887  s << "[" << rank << "]:" << " YaspGrid maxlevel=" << grid.maxLevel() << std::endl;
1888 
1889  s << "Printing the torus: " <<std::endl;
1890  s << grid.torus() << std::endl;
1891 
1892  for (typename YaspGrid<d,CC>::YGridLevelIterator g=grid.begin(); g!=grid.end(); ++g)
1893  {
1894  s << "[" << rank << "]: " << std::endl;
1895  s << "[" << rank << "]: " << "==========================================" << std::endl;
1896  s << "[" << rank << "]: " << "level=" << g->level() << std::endl;
1897 
1898  for (int codim = 0; codim < d + 1; ++codim)
1899  {
1900  s << "[" << rank << "]: " << "overlapfront[" << codim << "]: " << g->overlapfront[codim] << std::endl;
1901  s << "[" << rank << "]: " << "overlap[" << codim << "]: " << g->overlap[codim] << std::endl;
1902  s << "[" << rank << "]: " << "interiorborder[" << codim << "]: " << g->interiorborder[codim] << std::endl;
1903  s << "[" << rank << "]: " << "interior[" << codim << "]: " << g->interior[codim] << std::endl;
1904 
1905  typedef typename YGridList<CC>::Iterator I;
1906  for (I i=g->send_overlapfront_overlapfront[codim].begin();
1907  i!=g->send_overlapfront_overlapfront[codim].end(); ++i)
1908  s << "[" << rank << "]: " << " s_of_of[" << codim << "] to rank "
1909  << i->rank << " " << i->grid << std::endl;
1910 
1911  for (I i=g->recv_overlapfront_overlapfront[codim].begin();
1912  i!=g->recv_overlapfront_overlapfront[codim].end(); ++i)
1913  s << "[" << rank << "]: " << " r_of_of[" << codim << "] to rank "
1914  << i->rank << " " << i->grid << std::endl;
1915 
1916  for (I i=g->send_overlap_overlapfront[codim].begin();
1917  i!=g->send_overlap_overlapfront[codim].end(); ++i)
1918  s << "[" << rank << "]: " << " s_o_of[" << codim << "] to rank "
1919  << i->rank << " " << i->grid << std::endl;
1920 
1921  for (I i=g->recv_overlapfront_overlap[codim].begin();
1922  i!=g->recv_overlapfront_overlap[codim].end(); ++i)
1923  s << "[" << rank << "]: " << " r_of_o[" << codim << "] to rank "
1924  << i->rank << " " << i->grid << std::endl;
1925 
1926  for (I i=g->send_interiorborder_interiorborder[codim].begin();
1927  i!=g->send_interiorborder_interiorborder[codim].end(); ++i)
1928  s << "[" << rank << "]: " << " s_ib_ib[" << codim << "] to rank "
1929  << i->rank << " " << i->grid << std::endl;
1930 
1931  for (I i=g->recv_interiorborder_interiorborder[codim].begin();
1932  i!=g->recv_interiorborder_interiorborder[codim].end(); ++i)
1933  s << "[" << rank << "]: " << " r_ib_ib[" << codim << "] to rank "
1934  << i->rank << " " << i->grid << std::endl;
1935 
1936  for (I i=g->send_interiorborder_overlapfront[codim].begin();
1937  i!=g->send_interiorborder_overlapfront[codim].end(); ++i)
1938  s << "[" << rank << "]: " << " s_ib_of[" << codim << "] to rank "
1939  << i->rank << " " << i->grid << std::endl;
1940 
1941  for (I i=g->recv_overlapfront_interiorborder[codim].begin();
1942  i!=g->recv_overlapfront_interiorborder[codim].end(); ++i)
1943  s << "[" << rank << "]: " << " r_of_ib[" << codim << "] to rank "
1944  << i->rank << " " << i->grid << std::endl;
1945  }
1946  }
1947 
1948  s << std::endl;
1949 
1950  return s;
1951  }
1952 
1953  namespace Capabilities
1954  {
1955 
1963  template<int dim, class Coordinates>
1964  struct hasBackupRestoreFacilities< YaspGrid<dim, Coordinates> >
1965  {
1966  static const bool v = true;
1967  };
1968 
1972  template<int dim, class Coordinates>
1973  struct hasSingleGeometryType< YaspGrid<dim, Coordinates> >
1974  {
1975  static const bool v = true;
1976  static const unsigned int topologyId = GeometryTypes::cube(dim).id();
1977  };
1978 
1982  template<int dim, class Coordinates>
1983  struct isCartesian< YaspGrid<dim, Coordinates> >
1984  {
1985  static const bool v = true;
1986  };
1987 
1991  template<int dim, class Coordinates, int codim>
1992  struct hasEntity< YaspGrid<dim, Coordinates>, codim>
1993  {
1994  static const bool v = true;
1995  };
1996 
2001  template<int dim, class Coordinates, int codim>
2002  struct hasEntityIterator<YaspGrid<dim, Coordinates>, codim>
2003  {
2004  static const bool v = true;
2005  };
2006 
2010  template<int dim, int codim, class Coordinates>
2011  struct canCommunicate< YaspGrid< dim, Coordinates>, codim >
2012  {
2013  static const bool v = true;
2014  };
2015 
2019  template<int dim, class Coordinates>
2020  struct isLevelwiseConforming< YaspGrid<dim, Coordinates> >
2021  {
2022  static const bool v = true;
2023  };
2024 
2028  template<int dim, class Coordinates>
2029  struct isLeafwiseConforming< YaspGrid<dim, Coordinates> >
2030  {
2031  static const bool v = true;
2032  };
2033 
2034  }
2035 
2036 } // end namespace
2037 
2038 // Include the specialization of the StructuredGridFactory class for YaspGrid
2040 // Include the specialization of the BackupRestoreFacility class for YaspGrid
2041 #include <dune/grid/yaspgrid/backuprestore.hh>
2042 
2043 #endif
A geometry implementation for axis-aligned hypercubes.
Portable very large unsigned integers.
CommDataHandleIF describes the features of a data handle for communication in parallel runs using the...
Definition: datahandleif.hh:76
Wrapper class for entities.
Definition: entity.hh:64
Container for equidistant coordinates in a YaspGrid.
Definition: coordinates.hh:27
Container for equidistant coordinates in a YaspGrid with non-trivial origin.
Definition: coordinates.hh:129
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:279
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:644
constexpr bool isCube() const
Return true if entity is a cube of any dimension.
Definition: type.hh:634
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:649
Definition: grid.hh:855
static std::conditional< std::is_reference< InterfaceType >::value, typename std::add_lvalue_reference< typename ReturnImplementationType< typename std::remove_reference< InterfaceType >::type >::ImplementationType >::type, typename std::remove_const< typename ReturnImplementationType< typename std::remove_reference< InterfaceType >::type >::ImplementationType >::type >::type getRealImplementation(InterfaceType &&i)
return real implementation of interface class
Definition: grid.hh:1027
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
GridFamily::Traits::CollectiveCommunication CollectiveCommunication
A type that is a model of Dune::CollectiveCommunication. It provides a portable way for collective co...
Definition: grid.hh:519
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: intersection.hh:162
A Vector class with statically reserved memory.
Definition: reservedvector.hh:43
size_type size() const
Returns number of elements in the vector.
Definition: reservedvector.hh:184
void resize(size_t s)
Specifies a new size for the vector.
Definition: reservedvector.hh:98
reference back()
Returns reference to last element of vector.
Definition: reservedvector.hh:166
void pop_back()
Erases the last element of the vector, O(1) time.
Definition: reservedvector.hh:112
Coordinate container for a tensor product YaspGrid.
Definition: coordinates.hh:243
Definition: torus.hh:274
Definition: torus.hh:44
Definition: ygrid.hh:73
YGridComponent< Coordinates > move(iTupel v) const
return grid moved by the vector v
Definition: ygrid.hh:261
YGridComponent< Coordinates > intersection(const YGridComponent< Coordinates > &r) const
Return YGridComponent of supergrid of self which is the intersection of self and another YGridCompone...
Definition: ygrid.hh:269
implements a collection of multiple std::deque<Intersection> Intersections with neighboring processor...
Definition: ygrid.hh:821
int size() const
return the size of the container, this is the sum of the sizes of all deques
Definition: ygrid.hh:951
Iterator end() const
return iterator pointing to the end of the container
Definition: ygrid.hh:927
Iterator begin() const
return iterator pointing to the begin of the container
Definition: ygrid.hh:921
Iterator over a collection o YGrids A YGrid::Iterator is the heart of an entity in YaspGrid.
Definition: ygrid.hh:592
implements a collection of YGridComponents which form a codimension Entities of given codimension c n...
Definition: ygrid.hh:549
Implement the default load balance strategy of yaspgrid.
Definition: partitioning.hh:35
a base class for the yaspgrid partitioning strategy The name might be irritating. It will probably ch...
Definition: partitioning.hh:24
persistent, globally unique Ids
Definition: yaspgrididset.hh:23
[ provides Dune::Grid ]
Definition: yaspgrid.hh:165
YaspGridFamily< dim, Coordinates > GridFamily
the GridFamily of this grid
Definition: yaspgrid.hh:717
int ghostSize(int level, int codim) const
return size (= distance in graph) of ghost region
Definition: yaspgrid.hh:1421
YaspGrid(std::array< std::vector< ctype >, dim > coords, std::bitset< dim > periodic=std::bitset< dim >(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Standard constructor for a tensorproduct YaspGrid.
Definition: yaspgrid.hh:1047
int size(int level, GeometryType type) const
number of entities per level and geometry type in this process
Definition: yaspgrid.hh:1453
Traits::template Codim< cd >::template Partition< All_Partition >::LeafIterator leafend() const
return LeafIterator which points behind the last entity in maxLevel
Definition: yaspgrid.hh:1386
void globalRefine(int refCount)
refine the grid refCount times.
Definition: yaspgrid.hh:1224
int getMark(const typename Traits::template Codim< 0 >::Entity &e) const
returns adaptation mark for given entity
Definition: yaspgrid.hh:1308
Traits::template Codim< cd >::template Partition< All_Partition >::LevelIterator lbegin(int level) const
version without second template parameter for convenience
Definition: yaspgrid.hh:1351
int ghostSize(int codim) const
return size (= distance in graph) of ghost region
Definition: yaspgrid.hh:1427
int overlapSize(int codim) const
return size (= distance in graph) of overlap region
Definition: yaspgrid.hh:1414
int globalSize(int i) const
return number of cells on finest level in given direction on all processors
Definition: yaspgrid.hh:255
Traits::template Codim< cd >::template Partition< All_Partition >::LeafIterator leafbegin() const
return LeafIterator which points to the first entity in maxLevel
Definition: yaspgrid.hh:1379
const Dune::FieldVector< ctype, dim > & domainSize() const
returns the size of the physical domain
Definition: yaspgrid.hh:1471
void postAdapt()
clean up some markers
Definition: yaspgrid.hh:1329
YGridLevelIterator end() const
return iterator pointing to one past the finest level
Definition: yaspgrid.hh:310
const CollectiveCommunicationType & comm() const
return a collective communication object
Definition: yaspgrid.hh:1765
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition: yaspgrid.hh:1459
YaspGrid(Dune::FieldVector< ctype, dim > lowerleft, Dune::FieldVector< ctype, dim > upperright, std::array< int, dim > s, std::bitset< dim > periodic=std::bitset< dim >(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Definition: yaspgrid.hh:976
int maxLevel() const
Definition: yaspgrid.hh:1218
Traits::template Codim< cd >::template Partition< pitype >::LeafIterator leafend() const
return LeafIterator which points behind the last entity in maxLevel
Definition: yaspgrid.hh:1372
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir, int level) const
Definition: yaspgrid.hh:1480
void intersections(const YGridComponent< Coordinates > &sendgrid, const YGridComponent< Coordinates > &recvgrid, std::deque< Intersection > &sendlist, std::deque< Intersection > &recvlist)
Construct list of intersections with neighboring processors.
Definition: yaspgrid.hh:563
bool preAdapt()
returns true, if the grid will be coarsened
Definition: yaspgrid.hh:1321
iTupel levelSize(int l) const
return size vector of the grid (in cells) on level l
Definition: yaspgrid.hh:273
bool mark(int refCount, const typename Traits::template Codim< 0 >::Entity &e)
Marks an entity to be refined/coarsened in a subsequent adapt.
Definition: yaspgrid.hh:1294
int overlapSize(int level, int codim) const
return size (= distance in graph) of overlap region
Definition: yaspgrid.hh:1407
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir) const
Definition: yaspgrid.hh:1490
int size(int codim) const
number of leaf entities per codim in this process
Definition: yaspgrid.hh:1447
bool isPeriodic(int i) const
return whether the grid is periodic in direction i
Definition: yaspgrid.hh:282
void refineOptions(bool keepPhysicalOverlap)
set options for refinement
Definition: yaspgrid.hh:1278
YaspGrid(Dune::FieldVector< ctype, dim > L, std::array< int, dim > s, std::bitset< dim > periodic=std::bitset< dim >(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Definition: yaspgrid.hh:904
int size(int level, int codim) const
number of entities per level and codim in this process
Definition: yaspgrid.hh:1433
const Torus< CollectiveCommunicationType, dim > & torus() const
return reference to torus
Definition: yaspgrid.hh:249
Traits::template Codim< cd >::template Partition< All_Partition >::LevelIterator lend(int level) const
version without second template parameter for convenience
Definition: yaspgrid.hh:1358
bool adapt()
map adapt to global refine
Definition: yaspgrid.hh:1314
ReservedVector< YGridLevel, 32 >::const_iterator YGridLevelIterator
Iterator over the grid levels.
Definition: yaspgrid.hh:293
size_t numBoundarySegments() const
returns the number of boundary segments within the macro grid
Definition: yaspgrid.hh:1465
void communicateCodim(DataHandle &data, InterfaceType iftype, CommunicationDirection dir, int level) const
Definition: yaspgrid.hh:1500
Traits::template Codim< cd >::template Partition< pitype >::LevelIterator lend(int level) const
Iterator to one past the last entity of given codim on level for partition type.
Definition: yaspgrid.hh:1344
void makelevel(const Coordinates &coords, std::bitset< dim > periodic, iTupel o_interior, int overlap)
Make a new YGridLevel structure.
Definition: yaspgrid.hh:330
Traits::template Codim< cd >::template Partition< pitype >::LeafIterator leafbegin() const
return LeafIterator which points to the first entity in maxLevel
Definition: yaspgrid.hh:1365
Traits::template Codim< cd >::template Partition< pitype >::LevelIterator lbegin(int level) const
one past the end on this level
Definition: yaspgrid.hh:1337
iTupel globalSize() const
return number of cells on finest level on all processors
Definition: yaspgrid.hh:261
YGridLevelIterator begin(int i) const
return iterator pointing to given level
Definition: yaspgrid.hh:302
Coordinates::ctype ctype
Type used for coordinates.
Definition: yaspgrid.hh:179
YaspGrid(const Coordinates &coordinates, std::bitset< dim > periodic=std::bitset< dim >(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Definition: yaspgrid.hh:733
int levelSize(int l, int i) const
return size of the grid (in cells) on level l in direction i
Definition: yaspgrid.hh:267
YGridLevelIterator begin() const
return iterator pointing to coarsest level
Definition: yaspgrid.hh:296
YaspHierarchicIterator enables iteration over son entities of codim 0.
Definition: yaspgridhierarchiciterator.hh:18
Implementation of Level- and LeafIndexSets for YaspGrid.
Definition: yaspgridindexsets.hh:23
YaspIntersectionIterator enables iteration over intersections with neighboring codim 0 entities.
Definition: yaspgridintersectioniterator.hh:20
YaspIntersection provides data about intersection with neighboring codim 0 entities.
Definition: yaspgridintersection.hh:20
Iterates over entities of one grid level.
Definition: yaspgridleveliterator.hh:17
implement a consecutive index for all entities of given codim of a YaspGrid
Definition: yaspgridpersistentcontainer.hh:33
A set of traits classes to store static information about grid implementation.
Different resources needed by all grid implementations.
Implements an utility class that provides collective communication methods for sequential programs.
This provides container classes for the coordinates to be used in YaspGrid Upon implementation of the...
Describes the parallel communication interface class for MessageBuffers and DataHandles.
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
Various implementations of the power function for run-time and static arguments.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
CommunicationDirection
Define a type for communication direction parameter.
Definition: gridenums.hh:168
InterfaceType
Parameter to be used for the communication functions.
Definition: gridenums.hh:84
@ All_Partition
all entities
Definition: gridenums.hh:139
@ Interior_Partition
only interior entities
Definition: gridenums.hh:135
@ InteriorBorder_Partition
interior and border entities
Definition: gridenums.hh:136
@ Overlap_Partition
interior, border, and overlap entities
Definition: gridenums.hh:137
@ Ghost_Partition
only ghost entities
Definition: gridenums.hh:140
@ BackwardCommunication
reverse communication direction
Definition: gridenums.hh:170
@ InteriorBorder_All_Interface
send interior and border, receive all entities
Definition: gridenums.hh:86
@ All_All_Interface
send all and receive all entities
Definition: gridenums.hh:89
@ Overlap_All_Interface
send overlap, receive all entities
Definition: gridenums.hh:88
@ Overlap_OverlapFront_Interface
send overlap, receive overlap and front entities
Definition: gridenums.hh:87
@ InteriorBorder_InteriorBorder_Interface
send/receive interior and border entities
Definition: gridenums.hh:85
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:775
decltype(auto) ifElse(const Condition &condition, IfFunc &&ifFunc, ElseFunc &&elseFunc)
A conditional expression.
Definition: hybridutilities.hh:355
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Provides base classes for index and id sets.
Implements an utility class that provides MPI's collective communication methods.
Helpers for dealing with MPI.
Dune namespace.
Definition: alignedallocator.hh:14
An stl-compliant random-access container which stores everything on the stack.
specialize with 'true' for all codims that a grid can communicate data on (default=false)
Definition: capabilities.hh:95
specialize with 'true' for all codims that a grid provides an iterator for (default=false)
Definition: capabilities.hh:72
Specialize with 'true' for all codims that a grid implements entities for. (default=false)
Definition: capabilities.hh:56
Specialize with 'true' for if the codimension 0 entity of the grid has only one possible geometry typ...
Definition: capabilities.hh:25
Specialize with 'true' if the grid is a Cartesian grid. Cartesian grids satisfy the following propert...
Definition: capabilities.hh:46
Specialize with 'true' if implementation guarantees a conforming leaf grid. (default=false)
Definition: capabilities.hh:113
Specialize with 'true' if implementation guarantees conforming level grids. (default=false)
Definition: capabilities.hh:104
Static tag representing a codimension.
Definition: dimension.hh:22
A traits struct that collects all associated types of one grid model.
Definition: grid.hh:1065
IdSet< const GridImp, LocalIdSetImp, LIDType > LocalIdSet
The type of the local id set.
Definition: grid.hh:1138
IdSet< const GridImp, GlobalIdSetImp, GIDType > GlobalIdSet
The type of the global id set.
Definition: grid.hh:1136
IndexSet< const GridImp, LeafIndexSetImp > LeafIndexSet
The type of the leaf index set.
Definition: grid.hh:1134
IndexSet< const GridImp, LevelIndexSetImp > LevelIndexSet
The type of the level index set.
Definition: grid.hh:1132
Calculates m^p at compile time.
Definition: power.hh:24
type describing an intersection with a neighboring processor
Definition: ygrid.hh:827
Specialization of the StructuredGridFactory class for YaspGrid.
This file provides the infrastructure for toroidal communication in YaspGrid.
A unique label for each type of element that can occur in a grid.
Traits for type conversions and type information.
the YaspEntity class and its specializations
The YaspEntitySeed class.
The YaspGeometry class and its specializations.
level-wise, non-persistent, consecutive indices for YaspGrid
The YaspIntersection class.
The YaspIntersectionIterator class.
The YaspLevelIterator class.
Specialization of the PersistentContainer for YaspGrid.
This provides a YGrid, the elemental component of the yaspgrid implementation.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 13, 22:30, 2024)