Dune Core Modules (2.8.0)

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 <cstdint>
7#include <iostream>
8#include <vector>
9#include <algorithm>
10#include <stack>
11#include <type_traits>
12
13#include <dune/grid/common/backuprestore.hh>
14#include <dune/grid/common/grid.hh> // the grid base classes
15#include <dune/grid/common/capabilities.hh> // the capabilities
16#include <dune/common/hybridutilities.hh>
17#include <dune/common/power.hh>
24#include <dune/geometry/type.hh>
27
28
29#if HAVE_MPI
31#endif
32
40namespace Dune {
41
42 /* some sizes for building global ids
43 */
44 const int yaspgrid_dim_bits = 24; // bits for encoding each dimension
45 const int yaspgrid_level_bits = 5; // bits for encoding level number
46
47
48 //************************************************************************
49 // forward declaration of templates
50
51 template<int dim, class Coordinates> class YaspGrid;
52 template<int mydim, int cdim, class GridImp> class YaspGeometry;
53 template<int codim, int dim, class GridImp> class YaspEntity;
54 template<int codim, class GridImp> class YaspEntitySeed;
55 template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
56 template<class GridImp> class YaspIntersectionIterator;
57 template<class GridImp> class YaspIntersection;
58 template<class GridImp> class YaspHierarchicIterator;
59 template<class GridImp, bool isLeafIndexSet> class YaspIndexSet;
60 template<class GridImp> class YaspGlobalIdSet;
61 template<class GridImp> class YaspPersistentContainerIndex;
62
63} // namespace Dune
64
76#include <dune/grid/yaspgrid/yaspgrididset.hh>
78
79namespace Dune {
80
81#if HAVE_MPI
82 using YaspCollectiveCommunication = CollectiveCommunication<MPI_Comm>;
83#else
84 using YaspCollectiveCommunication = CollectiveCommunication<No_Comm>;
85#endif
86
87 template<int dim, class Coordinates>
88 struct YaspGridFamily
89 {
90 typedef YaspCollectiveCommunication CCType;
91
92 typedef GridTraits<dim, // dimension of the grid
93 dim, // dimension of the world space
95 YaspGeometry,YaspEntity,
96 YaspLevelIterator, // type used for the level iterator
97 YaspIntersection, // leaf intersection
98 YaspIntersection, // level intersection
99 YaspIntersectionIterator, // leaf intersection iter
100 YaspIntersectionIterator, // level intersection iter
101 YaspHierarchicIterator,
102 YaspLevelIterator, // type used for the leaf(!) iterator
103 YaspIndexSet< const YaspGrid< dim, Coordinates >, false >, // level index set
104 YaspIndexSet< const YaspGrid< dim, Coordinates >, true >, // leaf index set
105 YaspGlobalIdSet<const YaspGrid<dim, Coordinates> >,
106 bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
107 YaspGlobalIdSet<const YaspGrid<dim, Coordinates> >,
108 bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
109 CCType,
110 DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
111 YaspEntitySeed>
112 Traits;
113 };
114
115#ifndef DOXYGEN
116 template<int dim, int codim>
117 struct YaspCommunicateMeta {
118 template<class G, class DataHandle>
119 static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
120 {
121 if (data.contains(dim,codim))
122 {
123 g.template communicateCodim<DataHandle,codim>(data,iftype,dir,level);
124 }
125 YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
126 }
127 };
128
129 template<int dim>
130 struct YaspCommunicateMeta<dim,0> {
131 template<class G, class DataHandle>
132 static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
133 {
134 if (data.contains(dim,0))
135 g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
136 }
137 };
138#endif
139
140 //************************************************************************
157 template<int dim, class Coordinates = EquidistantCoordinates<double, dim> >
159 : public GridDefaultImplementation<dim,dim,typename Coordinates::ctype,YaspGridFamily<dim, Coordinates> >
160 {
161
162 template<int, PartitionIteratorType, typename>
163 friend class YaspLevelIterator;
164
165 template<typename>
166 friend class YaspHierarchicIterator;
167
168 protected:
169
170 public:
172 typedef typename Coordinates::ctype ctype;
174
175#ifndef DOXYGEN
176 typedef typename Dune::YGrid<Coordinates> YGrid;
178
181 struct YGridLevel {
182
184 int level() const
185 {
186 return level_;
187 }
188
189 Coordinates coords;
190
191 std::array<YGrid, dim+1> overlapfront;
192 std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> overlapfront_data;
193 std::array<YGrid, dim+1> overlap;
194 std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> overlap_data;
195 std::array<YGrid, dim+1> interiorborder;
196 std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> interiorborder_data;
197 std::array<YGrid, dim+1> interior;
198 std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> interior_data;
199
200 std::array<YGridList<Coordinates>,dim+1> send_overlapfront_overlapfront;
201 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_overlapfront_overlapfront_data;
202 std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_overlapfront;
203 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_overlapfront_data;
204
205 std::array<YGridList<Coordinates>,dim+1> send_overlap_overlapfront;
206 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_overlap_overlapfront_data;
207 std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_overlap;
208 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_overlap_data;
209
210 std::array<YGridList<Coordinates>,dim+1> send_interiorborder_interiorborder;
211 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_interiorborder_interiorborder_data;
212 std::array<YGridList<Coordinates>,dim+1> recv_interiorborder_interiorborder;
213 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_interiorborder_interiorborder_data;
214
215 std::array<YGridList<Coordinates>,dim+1> send_interiorborder_overlapfront;
216 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> send_interiorborder_overlapfront_data;
217 std::array<YGridList<Coordinates>,dim+1> recv_overlapfront_interiorborder;
218 std::array<std::deque<Intersection>, StaticPower<2,dim>::power> recv_overlapfront_interiorborder_data;
219
220 // general
221 YaspGrid<dim,Coordinates>* mg; // each grid level knows its multigrid
222 int overlapSize; // in mesh cells on this level
223 bool keepOverlap;
224
226 int level_;
227 };
228
230 typedef std::array<int, dim> iTupel;
231 typedef FieldVector<ctype, dim> fTupel;
232
233 // communication tag used by multigrid
234 enum { tag = 17 };
235#endif
236
239 {
240 return _torus;
241 }
242
244 int globalSize(int i) const
245 {
246 return levelSize(maxLevel(),i);
247 }
248
250 iTupel globalSize() const
251 {
252 return levelSize(maxLevel());
253 }
254
256 int levelSize(int l, int i) const
257 {
258 return _coarseSize[i] * (1 << l);
259 }
260
262 iTupel levelSize(int l) const
263 {
264 iTupel s;
265 for (int i=0; i<dim; ++i)
266 s[i] = levelSize(l,i);
267 return s;
268 }
269
271 bool isPeriodic(int i) const
272 {
273 return _periodic[i];
274 }
275
276 bool getRefineOption() const
277 {
278 return keep_ovlp;
279 }
280
283
286 {
287 return YGridLevelIterator(_levels,0);
288 }
289
292 {
293 if (i<0 || i>maxLevel())
294 DUNE_THROW(GridError, "level not existing");
295 return YGridLevelIterator(_levels,i);
296 }
297
300 {
301 return YGridLevelIterator(_levels,maxLevel()+1);
302 }
303
304 // static method to create the default load balance strategy
305 static const YLoadBalanceDefault<dim>* defaultLoadbalancer()
306 {
307 static YLoadBalanceDefault<dim> lb;
308 return & lb;
309 }
310
311 protected:
319 void makelevel (const Coordinates& coords, std::bitset<dim> periodic, iTupel o_interior, int overlap)
320 {
321 YGridLevel& g = _levels.back();
322 g.overlapSize = overlap;
323 g.mg = this;
324 g.level_ = maxLevel();
325 g.coords = coords;
326 g.keepOverlap = keep_ovlp;
327
328 // set the inserting positions in the corresponding arrays of YGridLevelStructure
329 typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator overlapfront_it = g.overlapfront_data.begin();
330 typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator overlap_it = g.overlap_data.begin();
331 typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator interiorborder_it = g.interiorborder_data.begin();
332 typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator interior_it = g.interior_data.begin();
333
334 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
335 send_overlapfront_overlapfront_it = g.send_overlapfront_overlapfront_data.begin();
336 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
337 recv_overlapfront_overlapfront_it = g.recv_overlapfront_overlapfront_data.begin();
338
339 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
340 send_overlap_overlapfront_it = g.send_overlap_overlapfront_data.begin();
341 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
342 recv_overlapfront_overlap_it = g.recv_overlapfront_overlap_data.begin();
343
344 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
345 send_interiorborder_interiorborder_it = g.send_interiorborder_interiorborder_data.begin();
346 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
347 recv_interiorborder_interiorborder_it = g.recv_interiorborder_interiorborder_data.begin();
348
349 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
350 send_interiorborder_overlapfront_it = g.send_interiorborder_overlapfront_data.begin();
351 typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator
352 recv_overlapfront_interiorborder_it = g.recv_overlapfront_interiorborder_data.begin();
353
354 // have a null array for constructor calls around
355 std::array<int,dim> n;
356 std::fill(n.begin(), n.end(), 0);
357
358 // determine origin of the grid with overlap and store whether an overlap area exists in direction i.
359 std::bitset<dim> ovlp_low(0ULL);
360 std::bitset<dim> ovlp_up(0ULL);
361
362 iTupel o_overlap;
363 iTupel s_overlap;
364
365 // determine at where we have overlap and how big the size of the overlap partition is
366 for (int i=0; i<dim; i++)
367 {
368 // the coordinate container has been contructed to hold the entire grid on
369 // this processor, including overlap. this is the element size.
370 s_overlap[i] = g.coords.size(i);
371
372 //in the periodic case there is always overlap
373 if (periodic[i])
374 {
375 o_overlap[i] = o_interior[i]-overlap;
376 ovlp_low[i] = true;
377 ovlp_up[i] = true;
378 }
379 else
380 {
381 //check lower boundary
382 if (o_interior[i] - overlap < 0)
383 o_overlap[i] = 0;
384 else
385 {
386 o_overlap[i] = o_interior[i] - overlap;
387 ovlp_low[i] = true;
388 }
389
390 //check upper boundary
391 if (o_overlap[i] + g.coords.size(i) < globalSize(i))
392 ovlp_up[i] = true;
393 }
394 }
395
396 for (unsigned int codim = 0; codim < dim + 1; codim++)
397 {
398 // set the begin iterator for the corresponding ygrids
399 g.overlapfront[codim].setBegin(overlapfront_it);
400 g.overlap[codim].setBegin(overlap_it);
401 g.interiorborder[codim].setBegin(interiorborder_it);
402 g.interior[codim].setBegin(interior_it);
403 g.send_overlapfront_overlapfront[codim].setBegin(send_overlapfront_overlapfront_it);
404 g.recv_overlapfront_overlapfront[codim].setBegin(recv_overlapfront_overlapfront_it);
405 g.send_overlap_overlapfront[codim].setBegin(send_overlap_overlapfront_it);
406 g.recv_overlapfront_overlap[codim].setBegin(recv_overlapfront_overlap_it);
407 g.send_interiorborder_interiorborder[codim].setBegin(send_interiorborder_interiorborder_it);
408 g.recv_interiorborder_interiorborder[codim].setBegin(recv_interiorborder_interiorborder_it);
409 g.send_interiorborder_overlapfront[codim].setBegin(send_interiorborder_overlapfront_it);
410 g.recv_overlapfront_interiorborder[codim].setBegin(recv_overlapfront_interiorborder_it);
411
412 // find all combinations of unit vectors that span entities of the given codimension
413 for (unsigned int index = 0; index < (1<<dim); index++)
414 {
415 // check whether the given shift is of our codimension
416 std::bitset<dim> r(index);
417 if (r.count() != dim-codim)
418 continue;
419
420 // get an origin and a size array for subsequent modification
421 std::array<int,dim> origin(o_overlap);
422 std::array<int,dim> size(s_overlap);
423
424 // build overlapfront
425 // we have to extend the element size by one in all directions without shift.
426 for (int i=0; i<dim; i++)
427 if (!r[i])
428 size[i]++;
429 *overlapfront_it = YGridComponent<Coordinates>(origin, r, &g.coords, size, n, size);
430
431 // build overlap
432 for (int i=0; i<dim; i++)
433 {
434 if (!r[i])
435 {
436 if (ovlp_low[i])
437 {
438 origin[i]++;
439 size[i]--;
440 }
441 if (ovlp_up[i])
442 size[i]--;
443 }
444 }
445 *overlap_it = YGridComponent<Coordinates>(origin,size,*overlapfront_it);
446
447 // build interiorborder
448 for (int i=0; i<dim; i++)
449 {
450 if (ovlp_low[i])
451 {
452 origin[i] += overlap;
453 size[i] -= overlap;
454 if (!r[i])
455 {
456 origin[i]--;
457 size[i]++;
458 }
459 }
460 if (ovlp_up[i])
461 {
462 size[i] -= overlap;
463 if (!r[i])
464 size[i]++;
465 }
466 }
467 *interiorborder_it = YGridComponent<Coordinates>(origin,size,*overlapfront_it);
468
469 // build interior
470 for (int i=0; i<dim; i++)
471 {
472 if (!r[i])
473 {
474 if (ovlp_low[i])
475 {
476 origin[i]++;
477 size[i]--;
478 }
479 if (ovlp_up[i])
480 size[i]--;
481 }
482 }
483 *interior_it = YGridComponent<Coordinates>(origin, size, *overlapfront_it);
484
485 intersections(*overlapfront_it,*overlapfront_it,*send_overlapfront_overlapfront_it, *recv_overlapfront_overlapfront_it);
486 intersections(*overlap_it,*overlapfront_it,*send_overlap_overlapfront_it, *recv_overlapfront_overlap_it);
487 intersections(*interiorborder_it,*interiorborder_it,*send_interiorborder_interiorborder_it,*recv_interiorborder_interiorborder_it);
488 intersections(*interiorborder_it,*overlapfront_it,*send_interiorborder_overlapfront_it,*recv_overlapfront_interiorborder_it);
489
490 // advance all iterators pointing to the next insertion point
491 ++overlapfront_it;
492 ++overlap_it;
493 ++interiorborder_it;
494 ++interior_it;
495 ++send_overlapfront_overlapfront_it;
496 ++recv_overlapfront_overlapfront_it;
497 ++send_overlap_overlapfront_it;
498 ++recv_overlapfront_overlap_it;
499 ++send_interiorborder_interiorborder_it;
500 ++recv_interiorborder_interiorborder_it;
501 ++send_interiorborder_overlapfront_it;
502 ++recv_overlapfront_interiorborder_it;
503 }
504
505 // set end iterators in the corresonding ygrids
506 g.overlapfront[codim].finalize(overlapfront_it);
507 g.overlap[codim].finalize(overlap_it);
508 g.interiorborder[codim].finalize(interiorborder_it);
509 g.interior[codim].finalize(interior_it);
510 g.send_overlapfront_overlapfront[codim].finalize(send_overlapfront_overlapfront_it,g.overlapfront[codim]);
511 g.recv_overlapfront_overlapfront[codim].finalize(recv_overlapfront_overlapfront_it,g.overlapfront[codim]);
512 g.send_overlap_overlapfront[codim].finalize(send_overlap_overlapfront_it,g.overlapfront[codim]);
513 g.recv_overlapfront_overlap[codim].finalize(recv_overlapfront_overlap_it,g.overlapfront[codim]);
514 g.send_interiorborder_interiorborder[codim].finalize(send_interiorborder_interiorborder_it,g.overlapfront[codim]);
515 g.recv_interiorborder_interiorborder[codim].finalize(recv_interiorborder_interiorborder_it,g.overlapfront[codim]);
516 g.send_interiorborder_overlapfront[codim].finalize(send_interiorborder_overlapfront_it,g.overlapfront[codim]);
517 g.recv_overlapfront_interiorborder[codim].finalize(recv_overlapfront_interiorborder_it,g.overlapfront[codim]);
518 }
519 }
520
521#ifndef DOXYGEN
530 struct mpifriendly_ygrid {
531 mpifriendly_ygrid ()
532 {
533 std::fill(origin.begin(), origin.end(), 0);
534 std::fill(size.begin(), size.end(), 0);
535 }
536 mpifriendly_ygrid (const YGridComponent<Coordinates>& grid)
537 : origin(grid.origin()), size(grid.size())
538 {}
539 iTupel origin;
540 iTupel size;
541 };
542#endif
543
553 std::deque<Intersection>& sendlist, std::deque<Intersection>& recvlist)
554 {
555 iTupel size = globalSize();
556
557 // the exchange buffers
558 std::vector<YGridComponent<Coordinates> > send_recvgrid(_torus.neighbors());
559 std::vector<YGridComponent<Coordinates> > recv_recvgrid(_torus.neighbors());
560 std::vector<YGridComponent<Coordinates> > send_sendgrid(_torus.neighbors());
561 std::vector<YGridComponent<Coordinates> > recv_sendgrid(_torus.neighbors());
562
563 // new exchange buffers to send simple struct without virtual functions
564 std::vector<mpifriendly_ygrid> mpifriendly_send_recvgrid(_torus.neighbors());
565 std::vector<mpifriendly_ygrid> mpifriendly_recv_recvgrid(_torus.neighbors());
566 std::vector<mpifriendly_ygrid> mpifriendly_send_sendgrid(_torus.neighbors());
567 std::vector<mpifriendly_ygrid> mpifriendly_recv_sendgrid(_torus.neighbors());
568
569 // fill send buffers; iterate over all neighboring processes
570 // non-periodic case is handled automatically because intersection will be zero
571 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
572 {
573 // determine if we communicate with this neighbor (and what)
574 bool skip = false;
575 iTupel coord = _torus.coord(); // my coordinates
576 iTupel delta = i.delta(); // delta to neighbor
577 iTupel nb = coord; // the neighbor
578 for (int k=0; k<dim; k++) nb[k] += delta[k];
579 iTupel v; // grid movement
580 std::fill(v.begin(), v.end(), 0);
581
582 for (int k=0; k<dim; k++)
583 {
584 if (nb[k]<0)
585 {
586 if (_periodic[k])
587 v[k] += size[k];
588 else
589 skip = true;
590 }
591 if (nb[k]>=_torus.dims(k))
592 {
593 if (_periodic[k])
594 v[k] -= size[k];
595 else
596 skip = true;
597 }
598 // neither might be true, then v=0
599 }
600
601 // store moved grids in send buffers
602 if (!skip)
603 {
604 send_sendgrid[i.index()] = sendgrid.move(v);
605 send_recvgrid[i.index()] = recvgrid.move(v);
606 }
607 else
608 {
609 send_sendgrid[i.index()] = YGridComponent<Coordinates>();
610 send_recvgrid[i.index()] = YGridComponent<Coordinates>();
611 }
612 }
613
614 // issue send requests for sendgrid being sent to all neighbors
615 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
616 {
617 mpifriendly_send_sendgrid[i.index()] = mpifriendly_ygrid(send_sendgrid[i.index()]);
618 _torus.send(i.rank(), &mpifriendly_send_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
619 }
620
621 // issue recv requests for sendgrids of neighbors
622 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
623 _torus.recv(i.rank(), &mpifriendly_recv_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
624
625 // exchange the sendgrids
626 _torus.exchange();
627
628 // issue send requests for recvgrid being sent to all neighbors
629 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
630 {
631 mpifriendly_send_recvgrid[i.index()] = mpifriendly_ygrid(send_recvgrid[i.index()]);
632 _torus.send(i.rank(), &mpifriendly_send_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
633 }
634
635 // issue recv requests for recvgrid of neighbors
636 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
637 _torus.recv(i.rank(), &mpifriendly_recv_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
638
639 // exchange the recvgrid
640 _torus.exchange();
641
642 // process receive buffers and compute intersections
643 for (typename Torus<CollectiveCommunicationType,dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
644 {
645 // what must be sent to this neighbor
646 Intersection send_intersection;
647 mpifriendly_ygrid yg = mpifriendly_recv_recvgrid[i.index()];
648 recv_recvgrid[i.index()] = YGridComponent<Coordinates>(yg.origin,yg.size);
649 send_intersection.grid = sendgrid.intersection(recv_recvgrid[i.index()]);
650 send_intersection.rank = i.rank();
651 send_intersection.distance = i.distance();
652 if (!send_intersection.grid.empty()) sendlist.push_front(send_intersection);
653
654 Intersection recv_intersection;
655 yg = mpifriendly_recv_sendgrid[i.index()];
656 recv_sendgrid[i.index()] = YGridComponent<Coordinates>(yg.origin,yg.size);
657 recv_intersection.grid = recvgrid.intersection(recv_sendgrid[i.index()]);
658 recv_intersection.rank = i.rank();
659 recv_intersection.distance = i.distance();
660 if(!recv_intersection.grid.empty()) recvlist.push_back(recv_intersection);
661 }
662 }
663
664 protected:
665
666 typedef const YaspGrid<dim,Coordinates> GridImp;
667
668 void init()
669 {
670 indexsets.push_back( std::make_shared< YaspIndexSet<const YaspGrid<dim, Coordinates>, false > >(*this,0) );
671 boundarysegmentssize();
672 }
673
674 void boundarysegmentssize()
675 {
676 // sizes of local macro grid
677 std::array<int, dim> sides;
678 {
679 for (int i=0; i<dim; i++)
680 {
681 sides[i] =
682 ((begin()->overlap[0].dataBegin()->origin(i) == 0)+
683 (begin()->overlap[0].dataBegin()->origin(i) + begin()->overlap[0].dataBegin()->size(i)
684 == levelSize(0,i)));
685 }
686 }
687 nBSegments = 0;
688 for (int k=0; k<dim; k++)
689 {
690 int offset = 1;
691 for (int l=0; l<dim; l++)
692 {
693 if (l==k) continue;
694 offset *= begin()->overlap[0].dataBegin()->size(l);
695 }
696 nBSegments += sides[k]*offset;
697 }
698 }
699
700 public:
701
702 // define the persistent index type
703 typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim> PersistentIndexType;
704
706 typedef YaspGridFamily<dim, Coordinates> GridFamily;
707 // the Traits
709
710 // need for friend declarations in entity
714
722 YaspGrid (const Coordinates& coordinates,
723 std::bitset<dim> periodic = std::bitset<dim>(0ULL),
724 int overlap = 1,
726 const YLoadBalance<dim>* lb = defaultLoadbalancer())
727 : ccobj(comm)
728 , leafIndexSet_(*this)
729 , _periodic(periodic)
730 , _overlap(overlap)
731 , keep_ovlp(true)
732 , adaptRefCount(0)
733 , adaptActive(false)
734 {
735 _levels.resize(1);
736
737 // Number of elements per coordinate direction on the coarsest level
738 for (std::size_t i=0; i<dim; i++)
739 _coarseSize[i] = coordinates.size(i);
740
741 // Construct the communication torus
742 _torus = decltype(_torus)(comm,tag,_coarseSize,lb);
743
744 iTupel o;
745 std::fill(o.begin(), o.end(), 0);
746 iTupel o_interior(o);
747 iTupel s_interior(_coarseSize);
748
749 _torus.partition(_torus.rank(),o,_coarseSize,o_interior,s_interior);
750
751 // Set domain size
752 if (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value
753 || std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >::value)
754 {
755 for (std::size_t i=0; i<dim; i++)
756 _L[i] = coordinates.size(i) * coordinates.meshsize(i,0);
757 }
758 if (std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >::value)
759 {
760 //determine sizes of vector to correctly construct torus structure and store for later size requests
761 for (int i=0; i<dim; i++)
762 _L[i] = coordinates.coordinate(i,_coarseSize[i]) - coordinates.coordinate(i,0);
763 }
764
765#if HAVE_MPI
766 // TODO: Settle on a single value for all coordinate types
767 int mysteryFactor = (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value) ? 1 : 2;
768
769 // check whether the grid is large enough to be overlapping
770 for (int i=0; i<dim; i++)
771 {
772 // find out whether the grid is too small to
773 int toosmall = (s_interior[i] <= mysteryFactor * overlap) && // interior is very small
774 (periodic[i] || (s_interior[i] != _coarseSize[i])); // there is an overlap in that direction
775 // communicate the result to all those processes to have all processors error out if one process failed.
776 int global = 0;
777 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
778 if (global)
779 DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
780 " Note that this also holds for DOFs on subdomain boundaries."
781 " Increase grid elements or decrease overlap accordingly.");
782 }
783#endif // #if HAVE_MPI
784
785 if (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value
786 || std::is_same<Coordinates,EquidistantOffsetCoordinates<ctype,dim> >::value)
787 {
788 iTupel s_overlap(s_interior);
789 for (int i=0; i<dim; i++)
790 {
791 if ((o_interior[i] - overlap > 0) || (periodic[i]))
792 s_overlap[i] += overlap;
793 if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
794 s_overlap[i] += overlap;
795 }
796
797 FieldVector<ctype,dim> upperRightWithOverlap;
798 for (int i=0; i<dim; i++)
799 upperRightWithOverlap[i] = coordinates.coordinate(i,0) + coordinates.meshsize(i,0) * s_overlap[i];
800
801 if constexpr (std::is_same_v<Coordinates,EquidistantCoordinates<ctype,dim>>)
802 {
803 // New coordinate object that additionally contains the overlap elements
804 EquidistantCoordinates<ctype,dim> coordinatesWithOverlap(upperRightWithOverlap,s_overlap);
805
806 // add level (the this-> is needed to make g++-6 happy)
807 this->makelevel(coordinatesWithOverlap,periodic,o_interior,overlap);
808 }
809
810 if constexpr (std::is_same_v<Coordinates,EquidistantOffsetCoordinates<ctype,dim>>)
811 {
813 for (int i=0; i<dim; i++)
814 lowerleft[i] = coordinates.origin(i);
815
816 // New coordinate object that additionally contains the overlap elements
817 EquidistantOffsetCoordinates<ctype,dim> coordinatesWithOverlap(lowerleft,upperRightWithOverlap,s_overlap);
818
819 // add level (the this-> is needed to make g++-6 happy)
820 this->makelevel(coordinatesWithOverlap,periodic,o_interior,overlap);
821 }
822 }
823
824 if constexpr (std::is_same_v<Coordinates,TensorProductCoordinates<ctype,dim>>)
825 {
826 std::array<std::vector<ctype>,dim> newCoords;
827 std::array<int, dim> offset(o_interior);
828
829 // find the relevant part of the coords vector for this processor and copy it to newCoords
830 for (int i=0; i<dim; ++i)
831 {
832 //define the coordinate range to be used
833 std::size_t begin = o_interior[i];
834 std::size_t end = begin + s_interior[i] + 1;
835
836 // check whether we are not at the physical boundary. In that case overlap is a simple
837 // extension of the coordinate range to be used
838 if (o_interior[i] - overlap > 0)
839 {
840 begin = begin - overlap;
841 offset[i] -= overlap;
842 }
843 if (o_interior[i] + s_interior[i] + overlap < _coarseSize[i])
844 end = end + overlap;
845
846 //copy the selected part in the new coord vector
847 newCoords[i].resize(end-begin);
848 auto newCoordsIt = newCoords[i].begin();
849 for (std::size_t j=begin; j<end; j++)
850 {
851 *newCoordsIt = coordinates.coordinate(i, j);
852 newCoordsIt++;
853 }
854
855 // Check whether we are at the physical boundary and have a periodic grid.
856 // In this case the coordinate vector has to be tweaked manually.
857 if ((periodic[i]) && (o_interior[i] + s_interior[i] + overlap >= _coarseSize[i]))
858 {
859 // we need to add the first <overlap> cells to the end of newcoords
860 for (int j=0; j<overlap; ++j)
861 newCoords[i].push_back(newCoords[i].back() - coordinates.coordinate(i,j) + coordinates.coordinate(i,j+1));
862 }
863
864 if ((periodic[i]) && (o_interior[i] - overlap <= 0))
865 {
866 offset[i] -= overlap;
867
868 // we need to add the last <overlap> cells to the begin of newcoords
869 std::size_t reverseCounter = coordinates.size(i);
870 for (int j=0; j<overlap; ++j, --reverseCounter)
871 newCoords[i].insert(newCoords[i].begin(), newCoords[i].front()
872 - coordinates.coordinate(i,reverseCounter) + coordinates.coordinate(i,reverseCounter-1));
873 }
874 }
875
876 TensorProductCoordinates<ctype,dim> coordinatesWithOverlap(newCoords, offset);
877
878 // add level (the this-> is needed to make g++-6 happy)
879 this->makelevel(coordinatesWithOverlap,periodic,o_interior,overlap);
880 }
881
882 init();
883 }
884
893 template<class C = Coordinates,
894 typename std::enable_if_t< std::is_same_v<C, EquidistantCoordinates<ctype,dim> >, int> = 0>
896 std::array<int, std::size_t{dim}> s,
897 std::bitset<std::size_t{dim}> periodic = std::bitset<std::size_t{dim}>{0ULL},
898 int overlap = 1,
899 CollectiveCommunicationType comm = CollectiveCommunicationType(),
900 const YLoadBalance<dim>* lb = defaultLoadbalancer())
901 : ccobj(comm), _torus(comm,tag,s,lb), leafIndexSet_(*this),
902 _L(L), _periodic(periodic), _coarseSize(s), _overlap(overlap),
903 keep_ovlp(true), adaptRefCount(0), adaptActive(false)
904 {
905 _levels.resize(1);
906
907 iTupel o;
908 std::fill(o.begin(), o.end(), 0);
909 iTupel o_interior(o);
910 iTupel s_interior(s);
911
912 _torus.partition(_torus.rank(),o,s,o_interior,s_interior);
913
914#if HAVE_MPI
915 // check whether the grid is large enough to be overlapping
916 for (int i=0; i<dim; i++)
917 {
918 // find out whether the grid is too small to
919 int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
920 (periodic[i] || (s_interior[i] != s[i])); // there is an overlap in that direction
921 // communicate the result to all those processes to have all processors error out if one process failed.
922 int global = 0;
923 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
924 if (global)
925 DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
926 " Note that this also holds for DOFs on subdomain boundaries."
927 " Increase grid elements or decrease overlap accordingly.");
928 }
929#endif // #if HAVE_MPI
930
931 iTupel s_overlap(s_interior);
932 for (int i=0; i<dim; i++)
933 {
934 if ((o_interior[i] - overlap > 0) || (periodic[i]))
935 s_overlap[i] += overlap;
936 if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
937 s_overlap[i] += overlap;
938 }
939
940 FieldVector<ctype,dim> upperRightWithOverlap;
941
942 for (int i=0; i<dim; i++)
943 upperRightWithOverlap[i] = (L[i] / s[i]) * s_overlap[i];
944
945 // New coordinate object that additionally contains the overlap elements
946 EquidistantCoordinates<ctype,dim> cc(upperRightWithOverlap,s_overlap);
947
948 // add level
949 makelevel(cc,periodic,o_interior,overlap);
950
951 init();
952 }
953
963 template<class C = Coordinates,
964 typename std::enable_if_t< std::is_same_v<C, EquidistantOffsetCoordinates<ctype,dim> >, int> = 0>
967 std::array<int, std::size_t{dim}> s,
968 std::bitset<std::size_t{dim}> periodic = std::bitset<std::size_t{dim}>(0ULL),
969 int overlap = 1,
970 CollectiveCommunicationType comm = CollectiveCommunicationType(),
971 const YLoadBalance<dim>* lb = defaultLoadbalancer())
972 : ccobj(comm), _torus(comm,tag,s,lb), leafIndexSet_(*this),
973 _L(upperright - lowerleft),
974 _periodic(periodic), _coarseSize(s), _overlap(overlap),
975 keep_ovlp(true), adaptRefCount(0), adaptActive(false)
976 {
977 _levels.resize(1);
978
979 iTupel o;
980 std::fill(o.begin(), o.end(), 0);
981 iTupel o_interior(o);
982 iTupel s_interior(s);
983
984 _torus.partition(_torus.rank(),o,s,o_interior,s_interior);
985
986#if HAVE_MPI
987 // check whether the grid is large enough to be overlapping
988 for (int i=0; i<dim; i++)
989 {
990 // find out whether the grid is too small to
991 int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
992 (periodic[i] || (s_interior[i] != s[i])); // there is an overlap in that direction
993 // communicate the result to all those processes to have all processors error out if one process failed.
994 int global = 0;
995 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
996 if (global)
997 DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
998 " Note that this also holds for DOFs on subdomain boundaries."
999 " Increase grid elements or decrease overlap accordingly.");
1000 }
1001#endif // #if HAVE_MPI
1002
1003 iTupel s_overlap(s_interior);
1004 for (int i=0; i<dim; i++)
1005 {
1006 if ((o_interior[i] - overlap > 0) || (periodic[i]))
1007 s_overlap[i] += overlap;
1008 if ((o_interior[i] + s_interior[i] + overlap <= _coarseSize[i]) || (periodic[i]))
1009 s_overlap[i] += overlap;
1010 }
1011
1012 FieldVector<ctype,dim> upperRightWithOverlap;
1013 for (int i=0; i<dim; i++)
1014 upperRightWithOverlap[i] = lowerleft[i]
1015 + s_overlap[i] * (upperright[i]-lowerleft[i]) / s[i];
1016
1017 EquidistantOffsetCoordinates<ctype,dim> cc(lowerleft,upperRightWithOverlap,s_overlap);
1018
1019 // add level
1020 makelevel(cc,periodic,o_interior,overlap);
1021
1022 init();
1023 }
1024
1032 template<class C = Coordinates,
1033 typename std::enable_if_t< std::is_same_v<C, TensorProductCoordinates<ctype,dim> >, int> = 0>
1034 YaspGrid (std::array<std::vector<ctype>, std::size_t{dim}> coords,
1035 std::bitset<std::size_t{dim}> periodic = std::bitset<std::size_t{dim}>(0ULL),
1036 int overlap = 1,
1037 CollectiveCommunicationType comm = CollectiveCommunicationType(),
1038 const YLoadBalance<dim>* lb = defaultLoadbalancer())
1039 : ccobj(comm), _torus(comm,tag,Dune::Yasp::sizeArray<dim>(coords),lb),
1040 leafIndexSet_(*this), _periodic(periodic), _overlap(overlap),
1041 keep_ovlp(true), adaptRefCount(0), adaptActive(false)
1042 {
1043 if (!Dune::Yasp::checkIfMonotonous(coords))
1044 DUNE_THROW(Dune::GridError,"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1045
1046 _levels.resize(1);
1047
1048 //determine sizes of vector to correctly construct torus structure and store for later size requests
1049 for (int i=0; i<dim; i++) {
1050 _coarseSize[i] = coords[i].size() - 1;
1051 _L[i] = coords[i][_coarseSize[i]] - coords[i][0];
1052 }
1053
1054 iTupel o;
1055 std::fill(o.begin(), o.end(), 0);
1056 iTupel o_interior(o);
1057 iTupel s_interior(_coarseSize);
1058
1059 _torus.partition(_torus.rank(),o,_coarseSize,o_interior,s_interior);
1060
1061#if HAVE_MPI
1062 // check whether the grid is large enough to be overlapping
1063 for (int i=0; i<dim; i++)
1064 {
1065 // find out whether the grid is too small to
1066 int toosmall = (s_interior[i] / 2 <= overlap) && // interior is very small
1067 (periodic[i] || (s_interior[i] != _coarseSize[i])); // there is an overlap in that direction
1068 // communicate the result to all those processes to have all processors error out if one process failed.
1069 int global = 0;
1070 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR, comm);
1071 if (global)
1072 DUNE_THROW(Dune::GridError,"YaspGrid does not support degrees of freedom shared by more than immediately neighboring subdomains."
1073 " Note that this also holds for DOFs on subdomain boundaries."
1074 " Increase grid elements or decrease overlap accordingly.");
1075 }
1076#endif // #if HAVE_MPI
1077
1078
1079 std::array<std::vector<ctype>,dim> newcoords;
1080 std::array<int, dim> offset(o_interior);
1081
1082 // find the relevant part of the coords vector for this processor and copy it to newcoords
1083 for (int i=0; i<dim; ++i)
1084 {
1085 //define iterators on coords that specify the coordinate range to be used
1086 typename std::vector<ctype>::iterator begin = coords[i].begin() + o_interior[i];
1087 typename std::vector<ctype>::iterator end = begin + s_interior[i] + 1;
1088
1089 // check whether we are not at the physical boundary. In that case overlap is a simple
1090 // extension of the coordinate range to be used
1091 if (o_interior[i] - overlap > 0)
1092 {
1093 begin = begin - overlap;
1094 offset[i] -= overlap;
1095 }
1096 if (o_interior[i] + s_interior[i] + overlap < _coarseSize[i])
1097 end = end + overlap;
1098
1099 //copy the selected part in the new coord vector
1100 newcoords[i].resize(end-begin);
1101 std::copy(begin, end, newcoords[i].begin());
1102
1103 // check whether we are at the physical boundary and a have a periodic grid.
1104 // In this case the coordinate vector has to be tweaked manually.
1105 if ((periodic[i]) && (o_interior[i] + s_interior[i] + overlap >= _coarseSize[i]))
1106 {
1107 // we need to add the first <overlap> cells to the end of newcoords
1108 typename std::vector<ctype>::iterator it = coords[i].begin();
1109 for (int j=0; j<overlap; ++j)
1110 newcoords[i].push_back(newcoords[i].back() - *it + *(++it));
1111 }
1112
1113 if ((periodic[i]) && (o_interior[i] - overlap <= 0))
1114 {
1115 offset[i] -= overlap;
1116
1117 // we need to add the last <overlap> cells to the begin of newcoords
1118 typename std::vector<ctype>::iterator it = coords[i].end() - 1;
1119 for (int j=0; j<overlap; ++j)
1120 newcoords[i].insert(newcoords[i].begin(), newcoords[i].front() - *it + *(--it));
1121 }
1122 }
1123
1124 TensorProductCoordinates<ctype,dim> cc(newcoords, offset);
1125
1126 // add level
1127 makelevel(cc,periodic,o_interior,overlap);
1128 init();
1129 }
1130
1131 private:
1132
1147 YaspGrid (std::array<std::vector<ctype>, std::size_t{dim}> coords,
1148 std::bitset<std::size_t{dim}> periodic,
1149 int overlap,
1150 CollectiveCommunicationType comm,
1151 std::array<int,dim> coarseSize,
1152 const YLoadBalance<dim>* lb = defaultLoadbalancer())
1153 : ccobj(comm), _torus(comm,tag,coarseSize,lb), leafIndexSet_(*this),
1154 _periodic(periodic), _coarseSize(coarseSize), _overlap(overlap),
1155 keep_ovlp(true), adaptRefCount(0), adaptActive(false)
1156 {
1157 // check whether YaspGrid has been given the correct template parameter
1158 static_assert(std::is_same<Coordinates,TensorProductCoordinates<ctype,dim> >::value,
1159 "YaspGrid coordinate container template parameter and given constructor values do not match!");
1160
1161 if (!Dune::Yasp::checkIfMonotonous(coords))
1162 DUNE_THROW(Dune::GridError,"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1163
1164 for (int i=0; i<dim; i++)
1165 _L[i] = coords[i][coords[i].size() - 1] - coords[i][0];
1166
1167 _levels.resize(1);
1168
1169 std::array<int,dim> o;
1170 std::fill(o.begin(), o.end(), 0);
1171 std::array<int,dim> o_interior(o);
1172 std::array<int,dim> s_interior(coarseSize);
1173
1174 _torus.partition(_torus.rank(),o,coarseSize,o_interior,s_interior);
1175
1176 // get offset by modifying o_interior according to overlap
1177 std::array<int,dim> offset(o_interior);
1178 for (int i=0; i<dim; i++)
1179 if ((periodic[i]) || (o_interior[i] > 0))
1180 offset[i] -= overlap;
1181
1182 TensorProductCoordinates<ctype,dim> cc(coords, offset);
1183
1184 // add level
1185 makelevel(cc,periodic,o_interior,overlap);
1186
1187 init();
1188 }
1189
1190 // the backup restore facility needs to be able to use above constructor
1191 friend struct BackupRestoreFacility<YaspGrid<dim,Coordinates> >;
1192
1193 // do not copy this class
1194 YaspGrid(const YaspGrid&);
1195
1196 public:
1197
1201 int maxLevel() const
1202 {
1203 return _levels.size()-1;
1204 }
1205
1207 void globalRefine (int refCount)
1208 {
1209 if (refCount < -maxLevel())
1210 DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
1211 "Coarsening " << -refCount << " levels requested!");
1212
1213 // If refCount is negative then coarsen the grid
1214 for (int k=refCount; k<0; k++)
1215 {
1216 // create an empty grid level
1217 YGridLevel empty;
1218 _levels.back() = empty;
1219 // reduce maxlevel
1220 _levels.pop_back();
1221
1222 indexsets.pop_back();
1223 }
1224
1225 // If refCount is positive refine the grid
1226 for (int k=0; k<refCount; k++)
1227 {
1228 // access to coarser grid level
1229 YGridLevel& cg = _levels[maxLevel()];
1230
1231 std::bitset<dim> ovlp_low(0ULL), ovlp_up(0ULL);
1232 for (int i=0; i<dim; i++)
1233 {
1234 if (cg.overlap[0].dataBegin()->origin(i) > 0 || _periodic[i])
1235 ovlp_low[i] = true;
1236 if (cg.overlap[0].dataBegin()->max(i) + 1 < globalSize(i) || _periodic[i])
1237 ovlp_up[i] = true;
1238 }
1239
1240 Coordinates newcont(cg.coords.refine(ovlp_low, ovlp_up, cg.overlapSize, keep_ovlp));
1241
1242 int overlap = (keep_ovlp) ? 2*cg.overlapSize : cg.overlapSize;
1243
1244 //determine new origin
1245 iTupel o_interior;
1246 for (int i=0; i<dim; i++)
1247 o_interior[i] = 2*cg.interior[0].dataBegin()->origin(i);
1248
1249 // add level
1250 _levels.resize(_levels.size() + 1);
1251 makelevel(newcont,_periodic,o_interior,overlap);
1252
1253 indexsets.push_back( std::make_shared<YaspIndexSet<const YaspGrid<dim,Coordinates>, false > >(*this,maxLevel()) );
1254 }
1255 }
1256
1261 void refineOptions (bool keepPhysicalOverlap)
1262 {
1263 keep_ovlp = keepPhysicalOverlap;
1264 }
1265
1277 bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
1278 {
1279 assert(adaptActive == false);
1280 if (e.level() != maxLevel()) return false;
1281 adaptRefCount = std::max(adaptRefCount, refCount);
1282 return true;
1283 }
1284
1291 int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
1292 {
1293 return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
1294 }
1295
1297 bool adapt ()
1298 {
1299 globalRefine(adaptRefCount);
1300 return (adaptRefCount > 0);
1301 }
1302
1304 bool preAdapt ()
1305 {
1306 adaptActive = true;
1307 adaptRefCount = comm().max(adaptRefCount);
1308 return (adaptRefCount < 0);
1309 }
1310
1313 {
1314 adaptActive = false;
1315 adaptRefCount = 0;
1316 }
1317
1319 template<int cd, PartitionIteratorType pitype>
1320 typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
1321 {
1322 return levelbegin<cd,pitype>(level);
1323 }
1324
1326 template<int cd, PartitionIteratorType pitype>
1327 typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
1328 {
1329 return levelend<cd,pitype>(level);
1330 }
1331
1333 template<int cd>
1334 typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
1335 {
1336 return levelbegin<cd,All_Partition>(level);
1337 }
1338
1340 template<int cd>
1341 typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
1342 {
1343 return levelend<cd,All_Partition>(level);
1344 }
1345
1347 template<int cd, PartitionIteratorType pitype>
1348 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const
1349 {
1350 return levelbegin<cd,pitype>(maxLevel());
1351 }
1352
1354 template<int cd, PartitionIteratorType pitype>
1355 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const
1356 {
1357 return levelend<cd,pitype>(maxLevel());
1358 }
1359
1361 template<int cd>
1362 typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
1363 {
1364 return levelbegin<cd,All_Partition>(maxLevel());
1365 }
1366
1368 template<int cd>
1369 typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
1370 {
1371 return levelend<cd,All_Partition>(maxLevel());
1372 }
1373
1374 // \brief obtain Entity from EntitySeed. */
1375 template <typename Seed>
1376 typename Traits::template Codim<Seed::codimension>::Entity
1377 entity(const Seed& seed) const
1378 {
1379 const int codim = Seed::codimension;
1380 YGridLevelIterator g = begin(seed.impl().level());
1381
1382 typedef typename Traits::template Codim<Seed::codimension>::Entity Entity;
1383 typedef YaspEntity<codim,dim,const YaspGrid> EntityImp;
1384 typedef typename YGrid::Iterator YIterator;
1385
1386 return Entity(EntityImp(g,YIterator(g->overlapfront[codim],seed.impl().coord(),seed.impl().offset())));
1387 }
1388
1390 int overlapSize (int level, [[maybe_unused]] int codim) const
1391 {
1392 YGridLevelIterator g = begin(level);
1393 return g->overlapSize;
1394 }
1395
1397 int overlapSize ([[maybe_unused]] int odim) const
1398 {
1400 return g->overlapSize;
1401 }
1402
1404 int ghostSize ([[maybe_unused]] int level, [[maybe_unused]] int codim) const
1405 {
1406 return 0;
1407 }
1408
1410 int ghostSize ([[maybe_unused]] int codim) const
1411 {
1412 return 0;
1413 }
1414
1416 int size (int level, int codim) const
1417 {
1418 YGridLevelIterator g = begin(level);
1419
1420 // sum over all components of the codimension
1421 int count = 0;
1422 typedef typename std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power>::iterator DAI;
1423 for (DAI it = g->overlapfront[codim].dataBegin(); it != g->overlapfront[codim].dataEnd(); ++it)
1424 count += it->totalsize();
1425
1426 return count;
1427 }
1428
1430 int size (int codim) const
1431 {
1432 return size(maxLevel(),codim);
1433 }
1434
1436 int size (int level, GeometryType type) const
1437 {
1438 return (type.isCube()) ? size(level,dim-type.dim()) : 0;
1439 }
1440
1442 int size (GeometryType type) const
1443 {
1444 return size(maxLevel(),type);
1445 }
1446
1448 size_t numBoundarySegments () const
1449 {
1450 return nBSegments;
1451 }
1452
1455 return _L;
1456 }
1457
1462 template<class DataHandleImp, class DataType>
1464 {
1465 YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
1466 }
1467
1472 template<class DataHandleImp, class DataType>
1474 {
1475 YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
1476 }
1477
1482 template<class DataHandle, int codim>
1483 void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
1484 {
1485 // check input
1486 if (!data.contains(dim,codim)) return; // should have been checked outside
1487
1488 // data types
1489 typedef typename DataHandle::DataType DataType;
1490
1491 // access to grid level
1492 YGridLevelIterator g = begin(level);
1493
1494 // find send/recv lists or throw error
1495 const YGridList<Coordinates>* sendlist = 0;
1496 const YGridList<Coordinates>* recvlist = 0;
1497
1499 {
1500 sendlist = &g->send_interiorborder_interiorborder[codim];
1501 recvlist = &g->recv_interiorborder_interiorborder[codim];
1502 }
1503 if (iftype==InteriorBorder_All_Interface)
1504 {
1505 sendlist = &g->send_interiorborder_overlapfront[codim];
1506 recvlist = &g->recv_overlapfront_interiorborder[codim];
1507 }
1509 {
1510 sendlist = &g->send_overlap_overlapfront[codim];
1511 recvlist = &g->recv_overlapfront_overlap[codim];
1512 }
1513 if (iftype==All_All_Interface)
1514 {
1515 sendlist = &g->send_overlapfront_overlapfront[codim];
1516 recvlist = &g->recv_overlapfront_overlapfront[codim];
1517 }
1518
1519 // change communication direction?
1520 if (dir==BackwardCommunication)
1521 std::swap(sendlist,recvlist);
1522
1523 int cnt;
1524
1525 // Size computation (requires communication if variable size)
1526 std::vector<int> send_size(sendlist->size(),-1); // map rank to total number of objects (of type DataType) to be sent
1527 std::vector<int> recv_size(recvlist->size(),-1); // map rank to total number of objects (of type DataType) to be recvd
1528 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
1529 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
1530
1531 // define type to iterate over send and recv lists
1532 typedef typename YGridList<Coordinates>::Iterator ListIt;
1533
1534 if (data.fixedSize(dim,codim))
1535 {
1536 // fixed size: just take a dummy entity, size can be computed without communication
1537 cnt=0;
1538 for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1539 {
1540 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1542 send_size[cnt] = is->grid.totalsize() * data.size(*it);
1543 cnt++;
1544 }
1545 cnt=0;
1546 for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1547 {
1548 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1550 recv_size[cnt] = is->grid.totalsize() * data.size(*it);
1551 cnt++;
1552 }
1553 }
1554 else
1555 {
1556 // variable size case: sender side determines the size
1557 cnt=0;
1558 for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1559 {
1560 // allocate send buffer for sizes per entitiy
1561 size_t *buf = new size_t[is->grid.totalsize()];
1562 send_sizes[cnt] = buf;
1563
1564 // loop over entities and ask for size
1565 int i=0; size_t n=0;
1566 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1568 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1569 itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1570 for ( ; it!=itend; ++it)
1571 {
1572 buf[i] = data.size(*it);
1573 n += buf[i];
1574 i++;
1575 }
1576
1577 // now we know the size for this rank
1578 send_size[cnt] = n;
1579
1580 // hand over send request to torus class
1581 torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1582 cnt++;
1583 }
1584
1585 // allocate recv buffers for sizes and store receive request
1586 cnt=0;
1587 for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1588 {
1589 // allocate recv buffer
1590 size_t *buf = new size_t[is->grid.totalsize()];
1591 recv_sizes[cnt] = buf;
1592
1593 // hand over recv request to torus class
1594 torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1595 cnt++;
1596 }
1597
1598 // exchange all size buffers now
1599 torus().exchange();
1600
1601 // release send size buffers
1602 cnt=0;
1603 for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1604 {
1605 delete[] send_sizes[cnt];
1606 send_sizes[cnt] = 0;
1607 cnt++;
1608 }
1609
1610 // process receive size buffers
1611 cnt=0;
1612 for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1613 {
1614 // get recv buffer
1615 size_t *buf = recv_sizes[cnt];
1616
1617 // compute total size
1618 size_t n=0;
1619 for (int i=0; i<is->grid.totalsize(); ++i)
1620 n += buf[i];
1621
1622 // ... and store it
1623 recv_size[cnt] = n;
1624 ++cnt;
1625 }
1626 }
1627
1628
1629 // allocate & fill the send buffers & store send request
1630 std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
1631 cnt=0;
1632 for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1633 {
1634 // allocate send buffer
1635 DataType *buf = new DataType[send_size[cnt]];
1636
1637 // remember send buffer
1638 sends[cnt] = buf;
1639
1640 // make a message buffer
1641 MessageBuffer<DataType> mb(buf);
1642
1643 // fill send buffer; iterate over cells in intersection
1644 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1646 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1647 itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1648 for ( ; it!=itend; ++it)
1649 data.gather(mb,*it);
1650
1651 // hand over send request to torus class
1652 torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
1653 cnt++;
1654 }
1655
1656 // allocate recv buffers and store receive request
1657 std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
1658 cnt=0;
1659 for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1660 {
1661 // allocate recv buffer
1662 DataType *buf = new DataType[recv_size[cnt]];
1663
1664 // remember recv buffer
1665 recvs[cnt] = buf;
1666
1667 // hand over recv request to torus class
1668 torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
1669 cnt++;
1670 }
1671
1672 // exchange all buffers now
1673 torus().exchange();
1674
1675 // release send buffers
1676 cnt=0;
1677 for (ListIt is=sendlist->begin(); is!=sendlist->end(); ++is)
1678 {
1679 delete[] sends[cnt];
1680 sends[cnt] = 0;
1681 cnt++;
1682 }
1683
1684 // process receive buffers and delete them
1685 cnt=0;
1686 for (ListIt is=recvlist->begin(); is!=recvlist->end(); ++is)
1687 {
1688 // get recv buffer
1689 DataType *buf = recvs[cnt];
1690
1691 // make a message buffer
1692 MessageBuffer<DataType> mb(buf);
1693
1694 // copy data from receive buffer; iterate over cells in intersection
1695 if (data.fixedSize(dim,codim))
1696 {
1697 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1699 size_t n=data.size(*it);
1700 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1701 itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1702 for ( ; it!=itend; ++it)
1703 data.scatter(mb,*it,n);
1704 }
1705 else
1706 {
1707 int i=0;
1708 size_t *sbuf = recv_sizes[cnt];
1709 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1711 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1712 itend(YaspLevelIterator<codim,All_Partition,GridImp>(g, typename YGrid::Iterator(is->yg,true)));
1713 for ( ; it!=itend; ++it)
1714 data.scatter(mb,*it,sbuf[i++]);
1715 delete[] sbuf;
1716 }
1717
1718 // delete buffer
1719 delete[] buf; // hier krachts !
1720 cnt++;
1721 }
1722 }
1723
1724 // The new index sets from DDM 11.07.2005
1725 const typename Traits::GlobalIdSet& globalIdSet() const
1726 {
1727 return theglobalidset;
1728 }
1729
1730 const typename Traits::LocalIdSet& localIdSet() const
1731 {
1732 return theglobalidset;
1733 }
1734
1735 const typename Traits::LevelIndexSet& levelIndexSet(int level) const
1736 {
1737 if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1738 return *(indexsets[level]);
1739 }
1740
1741 const typename Traits::LeafIndexSet& leafIndexSet() const
1742 {
1743 return leafIndexSet_;
1744 }
1745
1749 {
1750 return ccobj;
1751 }
1752
1753 private:
1754
1755 // number of boundary segments of the level 0 grid
1756 int nBSegments;
1757
1758 // Index classes need access to the real entity
1759 friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim, Coordinates>, true >;
1760 friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim, Coordinates>, false >;
1761 friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim, Coordinates> >;
1762 friend class Dune::YaspPersistentContainerIndex<const Dune::YaspGrid<dim, Coordinates> >;
1763
1764 friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim, Coordinates> >;
1765 friend class Dune::YaspIntersection<const Dune::YaspGrid<dim, Coordinates> >;
1766 friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim, Coordinates> >;
1767
1768 template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1769 friend class Entity;
1770
1771 template<class DT>
1772 class MessageBuffer {
1773 public:
1774 // Constructor
1775 MessageBuffer (DT *p)
1776 {
1777 a=p;
1778 i=0;
1779 j=0;
1780 }
1781
1782 // write data to message buffer, acts like a stream !
1783 template<class Y>
1784 void write (const Y& data)
1785 {
1786 static_assert(( std::is_same<DT,Y>::value ), "DataType mismatch");
1787 a[i++] = data;
1788 }
1789
1790 // read data from message buffer, acts like a stream !
1791 template<class Y>
1792 void read (Y& data) const
1793 {
1794 static_assert(( std::is_same<DT,Y>::value ), "DataType mismatch");
1795 data = a[j++];
1796 }
1797
1798 private:
1799 DT *a;
1800 int i;
1801 mutable int j;
1802 };
1803
1805 template<int cd, PartitionIteratorType pitype>
1806 YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
1807 {
1808 YGridLevelIterator g = begin(level);
1809 if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1810
1811 if (pitype==Interior_Partition)
1812 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interior[cd].begin());
1813 if (pitype==InteriorBorder_Partition)
1814 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interiorborder[cd].begin());
1815 if (pitype==Overlap_Partition)
1816 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlap[cd].begin());
1817 if (pitype<=All_Partition)
1818 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlapfront[cd].begin());
1819 if (pitype==Ghost_Partition)
1820 return levelend <cd, pitype> (level);
1821
1822 DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1823 }
1824
1826 template<int cd, PartitionIteratorType pitype>
1827 YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
1828 {
1829 YGridLevelIterator g = begin(level);
1830 if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1831
1832 if (pitype==Interior_Partition)
1833 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interior[cd].end());
1834 if (pitype==InteriorBorder_Partition)
1835 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interiorborder[cd].end());
1836 if (pitype==Overlap_Partition)
1837 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlap[cd].end());
1838 if (pitype<=All_Partition || pitype == Ghost_Partition)
1839 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlapfront[cd].end());
1840
1841 DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1842 }
1843
1844 CollectiveCommunicationType ccobj;
1845
1846 Torus<CollectiveCommunicationType,dim> _torus;
1847
1848 std::vector< std::shared_ptr< YaspIndexSet<const YaspGrid<dim,Coordinates>, false > > > indexsets;
1849 YaspIndexSet<const YaspGrid<dim,Coordinates>, true> leafIndexSet_;
1850 YaspGlobalIdSet<const YaspGrid<dim,Coordinates> > theglobalidset;
1851
1853 iTupel _s;
1854 std::bitset<dim> _periodic;
1855 iTupel _coarseSize;
1856 ReservedVector<YGridLevel,32> _levels;
1857 int _overlap;
1858 bool keep_ovlp;
1859 int adaptRefCount;
1860 bool adaptActive;
1861 };
1862
1863#if __cpp_deduction_guides >= 201611
1864 // Class template deduction guides
1865 template<typename ctype, int dim>
1866 YaspGrid(FieldVector<ctype, dim>,
1867 std::array<int, std::size_t{dim}>,
1868 std::bitset<std::size_t{dim}> = std::bitset<std::size_t{dim}>{0ULL},
1869 int = 1,
1870 YaspCollectiveCommunication = YaspCollectiveCommunication(),
1871 const YLoadBalance<dim>* = YaspGrid< dim, EquidistantCoordinates<ctype, dim> >::defaultLoadbalancer())
1872 -> YaspGrid< dim, EquidistantCoordinates<ctype, dim> >;
1873
1874 template<typename ctype, int dim>
1875 YaspGrid(FieldVector<ctype, dim>,
1876 FieldVector<ctype, dim>,
1877 std::array<int, std::size_t{dim}>,
1878 std::bitset<std::size_t{dim}> = std::bitset<std::size_t{dim}>{0ULL},
1879 int = 1,
1880 YaspCollectiveCommunication = YaspCollectiveCommunication(),
1881 const YLoadBalance<dim>* = YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >::defaultLoadbalancer())
1882 -> YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >;
1883
1884 template<typename ctype, std::size_t dim>
1885 YaspGrid(std::array<std::vector<ctype>, dim>,
1886 std::bitset<dim> = std::bitset<dim>{0ULL},
1887 int = 1,
1888 YaspCollectiveCommunication = YaspCollectiveCommunication(),
1889 const YLoadBalance<int{dim}>* = YaspGrid< int{dim}, TensorProductCoordinates<ctype, int{dim}> >::defaultLoadbalancer())
1890 -> YaspGrid< int{dim}, TensorProductCoordinates<ctype, int{dim}> >;
1891#endif
1892
1894 template <int d, class CC>
1895 std::ostream& operator<< (std::ostream& s, const YaspGrid<d,CC>& grid)
1896 {
1897 int rank = grid.torus().rank();
1898
1899 s << "[" << rank << "]:" << " YaspGrid maxlevel=" << grid.maxLevel() << std::endl;
1900
1901 s << "Printing the torus: " <<std::endl;
1902 s << grid.torus() << std::endl;
1903
1904 for (typename YaspGrid<d,CC>::YGridLevelIterator g=grid.begin(); g!=grid.end(); ++g)
1905 {
1906 s << "[" << rank << "]: " << std::endl;
1907 s << "[" << rank << "]: " << "==========================================" << std::endl;
1908 s << "[" << rank << "]: " << "level=" << g->level() << std::endl;
1909
1910 for (int codim = 0; codim < d + 1; ++codim)
1911 {
1912 s << "[" << rank << "]: " << "overlapfront[" << codim << "]: " << g->overlapfront[codim] << std::endl;
1913 s << "[" << rank << "]: " << "overlap[" << codim << "]: " << g->overlap[codim] << std::endl;
1914 s << "[" << rank << "]: " << "interiorborder[" << codim << "]: " << g->interiorborder[codim] << std::endl;
1915 s << "[" << rank << "]: " << "interior[" << codim << "]: " << g->interior[codim] << std::endl;
1916
1917 typedef typename YGridList<CC>::Iterator I;
1918 for (I i=g->send_overlapfront_overlapfront[codim].begin();
1919 i!=g->send_overlapfront_overlapfront[codim].end(); ++i)
1920 s << "[" << rank << "]: " << " s_of_of[" << codim << "] to rank "
1921 << i->rank << " " << i->grid << std::endl;
1922
1923 for (I i=g->recv_overlapfront_overlapfront[codim].begin();
1924 i!=g->recv_overlapfront_overlapfront[codim].end(); ++i)
1925 s << "[" << rank << "]: " << " r_of_of[" << codim << "] to rank "
1926 << i->rank << " " << i->grid << std::endl;
1927
1928 for (I i=g->send_overlap_overlapfront[codim].begin();
1929 i!=g->send_overlap_overlapfront[codim].end(); ++i)
1930 s << "[" << rank << "]: " << " s_o_of[" << codim << "] to rank "
1931 << i->rank << " " << i->grid << std::endl;
1932
1933 for (I i=g->recv_overlapfront_overlap[codim].begin();
1934 i!=g->recv_overlapfront_overlap[codim].end(); ++i)
1935 s << "[" << rank << "]: " << " r_of_o[" << codim << "] to rank "
1936 << i->rank << " " << i->grid << std::endl;
1937
1938 for (I i=g->send_interiorborder_interiorborder[codim].begin();
1939 i!=g->send_interiorborder_interiorborder[codim].end(); ++i)
1940 s << "[" << rank << "]: " << " s_ib_ib[" << codim << "] to rank "
1941 << i->rank << " " << i->grid << std::endl;
1942
1943 for (I i=g->recv_interiorborder_interiorborder[codim].begin();
1944 i!=g->recv_interiorborder_interiorborder[codim].end(); ++i)
1945 s << "[" << rank << "]: " << " r_ib_ib[" << codim << "] to rank "
1946 << i->rank << " " << i->grid << std::endl;
1947
1948 for (I i=g->send_interiorborder_overlapfront[codim].begin();
1949 i!=g->send_interiorborder_overlapfront[codim].end(); ++i)
1950 s << "[" << rank << "]: " << " s_ib_of[" << codim << "] to rank "
1951 << i->rank << " " << i->grid << std::endl;
1952
1953 for (I i=g->recv_overlapfront_interiorborder[codim].begin();
1954 i!=g->recv_overlapfront_interiorborder[codim].end(); ++i)
1955 s << "[" << rank << "]: " << " r_of_ib[" << codim << "] to rank "
1956 << i->rank << " " << i->grid << std::endl;
1957 }
1958 }
1959
1960 s << std::endl;
1961
1962 return s;
1963 }
1964
1965 namespace Capabilities
1966 {
1967
1975 template<int dim, class Coordinates>
1976 struct hasBackupRestoreFacilities< YaspGrid<dim, Coordinates> >
1977 {
1978 static const bool v = true;
1979 };
1980
1984 template<int dim, class Coordinates>
1985 struct hasSingleGeometryType< YaspGrid<dim, Coordinates> >
1986 {
1987 static const bool v = true;
1988 static const unsigned int topologyId = GeometryTypes::cube(dim).id();
1989 };
1990
1994 template<int dim, class Coordinates>
1995 struct isCartesian< YaspGrid<dim, Coordinates> >
1996 {
1997 static const bool v = true;
1998 };
1999
2003 template<int dim, class Coordinates, int codim>
2004 struct hasEntity< YaspGrid<dim, Coordinates>, codim>
2005 {
2006 static const bool v = true;
2007 };
2008
2013 template<int dim, class Coordinates, int codim>
2014 struct hasEntityIterator<YaspGrid<dim, Coordinates>, codim>
2015 {
2016 static const bool v = true;
2017 };
2018
2022 template<int dim, int codim, class Coordinates>
2023 struct canCommunicate< YaspGrid< dim, Coordinates>, codim >
2024 {
2025 static const bool v = true;
2026 };
2027
2031 template<int dim, class Coordinates>
2032 struct isLevelwiseConforming< YaspGrid<dim, Coordinates> >
2033 {
2034 static const bool v = true;
2035 };
2036
2040 template<int dim, class Coordinates>
2041 struct isLeafwiseConforming< YaspGrid<dim, Coordinates> >
2042 {
2043 static const bool v = true;
2044 };
2045
2046 }
2047
2048} // end namespace
2049
2050// Include the specialization of the StructuredGridFactory class for YaspGrid
2052// Include the specialization of the BackupRestoreFacility class for YaspGrid
2053#include <dune/grid/yaspgrid/backuprestore.hh>
2054
2055#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
T max(const T &in) const
Compute the maximum of the argument over all processes and return the result in every process....
Definition: communication.hh:246
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:123
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:369
constexpr bool isCube() const
Return true if entity is a cube of any dimension.
Definition: type.hh:333
constexpr unsigned int id() const
Return the topology id of the type.
Definition: type.hh:374
Definition: grid.hh:851
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
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:160
YaspGridFamily< dim, Coordinates > GridFamily
the GridFamily of this grid
Definition: yaspgrid.hh:706
int ghostSize(int level, int codim) const
return size (= distance in graph) of ghost region
Definition: yaspgrid.hh:1404
YaspGrid(Dune::FieldVector< ctype, dim > L, std::array< int, std::size_t{dim}> s, std::bitset< std::size_t{dim}> periodic=std::bitset< std::size_t{dim}>{0ULL}, int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Definition: yaspgrid.hh:895
Traits::template Codim< cd >::template Partition< All_Partition >::LevelIterator lbegin(int level) const
version without second template parameter for convenience
Definition: yaspgrid.hh:1334
int size(int level, GeometryType type) const
number of entities per level and geometry type in this process
Definition: yaspgrid.hh:1436
Traits::template Codim< cd >::template Partition< All_Partition >::LeafIterator leafbegin() const
return LeafIterator which points to the first entity in maxLevel
Definition: yaspgrid.hh:1362
void globalRefine(int refCount)
refine the grid refCount times.
Definition: yaspgrid.hh:1207
const Torus< CollectiveCommunicationType, dim > & torus() const
return reference to torus
Definition: yaspgrid.hh:238
int getMark(const typename Traits::template Codim< 0 >::Entity &e) const
returns adaptation mark for given entity
Definition: yaspgrid.hh:1291
int ghostSize(int codim) const
return size (= distance in graph) of ghost region
Definition: yaspgrid.hh:1410
int globalSize(int i) const
return number of cells on finest level in given direction on all processors
Definition: yaspgrid.hh:244
Traits::template Codim< cd >::template Partition< pitype >::LeafIterator leafend() const
return LeafIterator which points behind the last entity in maxLevel
Definition: yaspgrid.hh:1355
void postAdapt()
clean up some markers
Definition: yaspgrid.hh:1312
const Dune::FieldVector< ctype, dim > & domainSize() const
returns the size of the physical domain
Definition: yaspgrid.hh:1454
YGridLevelIterator end() const
return iterator pointing to one past the finest level
Definition: yaspgrid.hh:299
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition: yaspgrid.hh:1442
int maxLevel() const
Definition: yaspgrid.hh:1201
Traits::template Codim< cd >::template Partition< All_Partition >::LeafIterator leafend() const
return LeafIterator which points behind the last entity in maxLevel
Definition: yaspgrid.hh:1369
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir, int level) const
Definition: yaspgrid.hh:1463
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:552
Traits::template Codim< cd >::template Partition< pitype >::LeafIterator leafbegin() const
return LeafIterator which points to the first entity in maxLevel
Definition: yaspgrid.hh:1348
bool preAdapt()
returns true, if the grid will be coarsened
Definition: yaspgrid.hh:1304
iTupel levelSize(int l) const
return size vector of the grid (in cells) on level l
Definition: yaspgrid.hh:262
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:1277
int overlapSize(int level, int codim) const
return size (= distance in graph) of overlap region
Definition: yaspgrid.hh:1390
const CollectiveCommunicationType & comm() const
return a collective communication object
Definition: yaspgrid.hh:1748
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir) const
Definition: yaspgrid.hh:1473
int size(int codim) const
number of leaf entities per codim in this process
Definition: yaspgrid.hh:1430
bool isPeriodic(int i) const
return whether the grid is periodic in direction i
Definition: yaspgrid.hh:271
void refineOptions(bool keepPhysicalOverlap)
set options for refinement
Definition: yaspgrid.hh:1261
int size(int level, int codim) const
number of entities per level and codim in this process
Definition: yaspgrid.hh:1416
bool adapt()
map adapt to global refine
Definition: yaspgrid.hh:1297
Traits::template Codim< cd >::template Partition< All_Partition >::LevelIterator lend(int level) const
version without second template parameter for convenience
Definition: yaspgrid.hh:1341
ReservedVector< YGridLevel, 32 >::const_iterator YGridLevelIterator
Iterator over the grid levels.
Definition: yaspgrid.hh:282
size_t numBoundarySegments() const
returns the number of boundary segments within the macro grid
Definition: yaspgrid.hh:1448
void communicateCodim(DataHandle &data, InterfaceType iftype, CommunicationDirection dir, int level) const
Definition: yaspgrid.hh:1483
void makelevel(const Coordinates &coords, std::bitset< dim > periodic, iTupel o_interior, int overlap)
Make a new YGridLevel structure.
Definition: yaspgrid.hh:319
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:1327
Traits::template Codim< cd >::template Partition< pitype >::LevelIterator lbegin(int level) const
one past the end on this level
Definition: yaspgrid.hh:1320
YaspGrid(std::array< std::vector< ctype >, std::size_t{dim}> coords, std::bitset< std::size_t{dim}> periodic=std::bitset< std::size_t{dim}>(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Standard constructor for a tensorproduct YaspGrid.
Definition: yaspgrid.hh:1034
YaspGrid(Dune::FieldVector< ctype, dim > lowerleft, Dune::FieldVector< ctype, dim > upperright, std::array< int, std::size_t{dim}> s, std::bitset< std::size_t{dim}> periodic=std::bitset< std::size_t{dim}>(0ULL), int overlap=1, CollectiveCommunicationType comm=CollectiveCommunicationType(), const YLoadBalance< dim > *lb=defaultLoadbalancer())
Definition: yaspgrid.hh:965
iTupel globalSize() const
return number of cells on finest level on all processors
Definition: yaspgrid.hh:250
int overlapSize(int odim) const
return size (= distance in graph) of overlap region
Definition: yaspgrid.hh:1397
YGridLevelIterator begin(int i) const
return iterator pointing to given level
Definition: yaspgrid.hh:291
Coordinates::ctype ctype
Type used for coordinates.
Definition: yaspgrid.hh:172
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:722
int levelSize(int l, int i) const
return size of the grid (in cells) on level l in direction i
Definition: yaspgrid.hh:256
YGridLevelIterator begin() const
return iterator pointing to coarsest level
Definition: yaspgrid.hh:285
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.
Various implementations of the power function for run-time and static arguments.
Traits for type conversions and type information.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto periodic(RawPreBasisIndicator &&rawPreBasisIndicator, PIS &&periodicIndexSet)
Create a pre-basis factory that can create a periodic pre-basis.
Definition: periodicbasis.hh:195
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:470
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
constexpr HybridTreePath< T..., std::size_t > push_back(const HybridTreePath< T... > &tp, std::size_t i)
Appends a run time index to a HybridTreePath.
Definition: treepath.hh:278
constexpr auto back(const HybridTreePath< T... > &tp) -> decltype(treePathEntry< sizeof...(T) -1 >(tp))
Returns a copy of the last element of the HybridTreePath.
Definition: treepath.hh:254
constexpr auto front(const HybridTreePath< T... > &tp) -> decltype(treePathEntry< 0 >(tp))
Returns a copy of the first element of the HybridTreePath.
Definition: treepath.hh:267
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:11
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=hasEntity<codim>:...
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:984
IdSet< const GridImp, LocalIdSetImp, LIDType > LocalIdSet
The type of the local id set.
Definition: grid.hh:1057
IdSet< const GridImp, GlobalIdSetImp, GIDType > GlobalIdSet
The type of the global id set.
Definition: grid.hh:1055
IndexSet< const GridImp, LeafIndexSetImp > LeafIndexSet
The type of the leaf index set.
Definition: grid.hh:1053
IndexSet< const GridImp, LevelIndexSetImp > LevelIndexSet
The type of the level index set.
Definition: grid.hh:1051
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.
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.111.3 (Dec 22, 23:30, 2024)