3 #ifndef DUNE_GRID_YASPGRID_HH
4 #define DUNE_GRID_YASPGRID_HH
10 #include <type_traits>
22 #include <dune/common/hybridutilities.hh>
23 #include <dune/common/power.hh>
24 #include <dune/common/bigunsignedint.hh>
25 #include <dune/common/typetraits.hh>
26 #include <dune/common/reservedvector.hh>
27 #include <dune/common/parallel/communication.hh>
28 #include <dune/common/parallel/mpihelper.hh>
29 #include <dune/common/deprecated.hh>
30 #include <dune/geometry/axisalignedcubegeometry.hh>
31 #include <dune/geometry/type.hh>
37 #include <dune/common/parallel/mpicommunication.hh>
58 template<
int dim,
class Coordinates>
class YaspGrid;
60 template<
int codim,
int dim,
class Gr
idImp>
class YaspEntity;
88 template<
int dim,
class Coordinates>
92 typedef CollectiveCommunication<MPI_Comm>
CCType;
94 typedef CollectiveCommunication<No_Comm>
CCType;
111 bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
113 bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+dim>,
121 template<
int dim,
int codim>
122 struct YaspCommunicateMeta {
123 template<
class G,
class DataHandle>
126 if (data.contains(dim,codim))
128 g.template communicateCodim<DataHandle,codim>(data,iftype,dir,level);
130 YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
135 struct YaspCommunicateMeta<dim,0> {
136 template<
class G,
class DataHandle>
139 if (data.contains(dim,0))
140 g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
162 template<
int dim,
class Coordinates = Equ
idistantCoordinates<
double, dim> >
164 :
public GridDefaultImplementation<dim,dim,typename Coordinates::ctype,YaspGridFamily<dim, Coordinates> >
167 template<
int, PartitionIteratorType,
typename>
179 typedef typename Coordinates::ctype
ctype;
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;
209 std::array<YGridComponent<Coordinates>, StaticPower<2,dim>::power> interior_data;
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;
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;
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;
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;
241 typedef std::array<int, dim> iTupel;
242 typedef FieldVector<ctype, dim> fTupel;
269 return _coarseSize[i] * (1 << l);
276 for (
int i=0; i<dim; ++i)
305 DUNE_THROW(
GridError,
"level not existing");
330 void makelevel (
const Coordinates& coords, std::bitset<dim> periodic, iTupel o_interior,
int overlap)
332 YGridLevel& g = _levels.back();
337 g.keepOverlap = keep_ovlp;
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();
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();
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();
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();
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();
366 std::array<int,dim> n;
367 std::fill(n.begin(), n.end(), 0);
370 std::bitset<dim> ovlp_low(0ULL);
371 std::bitset<dim> ovlp_up(0ULL);
377 for (
int i=0; i<dim; i++)
381 s_overlap[i] = g.coords.size(i);
386 o_overlap[i] = o_interior[i]-
overlap;
393 if (o_interior[i] -
overlap < 0)
397 o_overlap[i] = o_interior[i] -
overlap;
402 if (o_overlap[i] + g.coords.size(i) <
globalSize(i))
407 for (
unsigned int codim = 0; codim < dim + 1; codim++)
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);
424 for (
unsigned int index = 0; index < (1<<dim); index++)
427 std::bitset<dim> r(index);
428 if (r.count() != dim-codim)
432 std::array<int,dim> origin(o_overlap);
433 std::array<int,dim>
size(s_overlap);
437 for (
int i=0; i<dim; i++)
443 for (
int i=0; i<dim; i++)
459 for (
int i=0; i<dim; i++)
481 for (
int i=0; i<dim; i++)
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);
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;
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]);
541 struct mpifriendly_ygrid {
544 std::fill(origin.begin(), origin.end(), 0);
545 std::fill(
size.begin(),
size.end(), 0);
547 mpifriendly_ygrid (
const YGridComponent<Coordinates>& grid)
548 : origin(grid.origin()),
size(grid.
size())
564 std::deque<Intersection>& sendlist, std::deque<Intersection>& recvlist)
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());
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());
586 iTupel coord = _torus.
coord();
587 iTupel delta = i.delta();
589 for (
int k=0; k<dim; k++) nb[k] += delta[k];
591 std::fill(v.begin(), v.end(), 0);
593 for (
int k=0; k<dim; k++)
602 if (nb[k]>=_torus.
dims(k))
615 send_sendgrid[i.index()] = sendgrid.
move(v);
616 send_recvgrid[i.index()] = recvgrid.
move(v);
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));
634 _torus.
recv(i.rank(), &mpifriendly_recv_sendgrid[i.index()],
sizeof(mpifriendly_ygrid));
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));
648 _torus.
recv(i.rank(), &mpifriendly_recv_recvgrid[i.index()],
sizeof(mpifriendly_ygrid));
658 mpifriendly_ygrid yg = mpifriendly_recv_recvgrid[i.index()];
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);
666 yg = mpifriendly_recv_sendgrid[i.index()];
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);
688 std::array<int, dim> sides;
690 for (
int i=0; i<dim; i++)
693 ((
begin()->overlap[0].dataBegin()->origin(i) == 0)+
694 (
begin()->overlap[0].dataBegin()->origin(i) +
begin()->overlap[0].dataBegin()->size(i)
699 for (
int k=0; k<dim; k++)
702 for (
int l=0; l<dim; l++)
705 offset *=
begin()->overlap[0].dataBegin()->size(l);
707 nBSegments += sides[k]*offset;
734 std::bitset<dim> periodic = std::bitset<dim>(0ULL),
739 , leafIndexSet_(*this)
740 , _periodic(periodic)
749 for (std::size_t i=0; i<dim; i++)
750 _coarseSize[i] = coordinates.size(i);
753 _torus = decltype(_torus)(
comm,tag,_coarseSize,lb);
756 std::fill(o.begin(), o.end(), 0);
757 iTupel o_interior(o);
758 iTupel s_interior(_coarseSize);
760 _torus.
partition(_torus.
rank(),o,_coarseSize,o_interior,s_interior);
766 for (std::size_t i=0; i<dim; i++)
767 _L[i] = coordinates.size(i) * coordinates.meshsize(i,0);
772 for (
int i=0; i<dim; i++)
773 _L[i] = coordinates.coordinate(i,_coarseSize[i]) - coordinates.coordinate(i,0);
778 int mysteryFactor = (std::is_same<Coordinates,EquidistantCoordinates<ctype,dim> >::value) ? 1 : 2;
781 for (
int i=0; i<dim; i++)
784 int toosmall = (s_interior[i] <= mysteryFactor *
overlap) &&
785 (periodic[i] || (s_interior[i] != _coarseSize[i]));
788 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR,
comm);
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.");
794 #endif // #if HAVE_MPI
799 iTupel s_overlap(s_interior);
800 for (
int i=0; i<dim; i++)
802 if ((o_interior[i] -
overlap > 0) || (periodic[i]))
804 if ((o_interior[i] + s_interior[i] +
overlap <= _coarseSize[i]) || (periodic[i]))
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];
823 Dune::FieldVector<ctype,dim> lowerleft;
824 for (
int i=0; i<dim; i++)
825 lowerleft[i] =
id(coordinates).
origin(i);
837 std::array<std::vector<ctype>,dim> newCoords;
838 std::array<int, dim> offset(o_interior);
841 for (
int i=0; i<dim; ++i)
844 std::size_t
begin = o_interior[i];
845 std::size_t
end =
begin + s_interior[i] + 1;
849 if (o_interior[i] -
overlap > 0)
854 if (o_interior[i] + s_interior[i] +
overlap < _coarseSize[i])
859 auto newCoordsIt = newCoords[i].begin();
860 for (std::size_t j=
begin; j<
end; j++)
862 *newCoordsIt = coordinates.coordinate(i, j);
868 if ((periodic[i]) && (o_interior[i] + s_interior[i] +
overlap >= _coarseSize[i]))
872 newCoords[i].push_back(newCoords[i].back() - coordinates.coordinate(i,j) + coordinates.coordinate(i,j+1));
875 if ((periodic[i]) && (o_interior[i] -
overlap <= 0))
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));
905 std::array<int, dim> s,
906 std::bitset<dim> periodic = std::bitset<dim>(0ULL),
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)
916 "YaspGrid coordinate container template parameter and given constructor values do not match!");
921 std::fill(o.begin(), o.end(), 0);
922 iTupel o_interior(o);
923 iTupel s_interior(s);
929 for (
int i=0; i<dim; i++)
932 int toosmall = (s_interior[i] / 2 <=
overlap) &&
933 (periodic[i] || (s_interior[i] != s[i]));
936 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR,
comm);
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.");
942 #endif // #if HAVE_MPI
944 iTupel s_overlap(s_interior);
945 for (
int i=0; i<dim; i++)
947 if ((o_interior[i] -
overlap > 0) || (periodic[i]))
949 if ((o_interior[i] + s_interior[i] +
overlap <= _coarseSize[i]) || (periodic[i]))
953 FieldVector<ctype,dim> upperRightWithOverlap;
955 for (
int i=0; i<dim; i++)
956 upperRightWithOverlap[i] = (L[i] / s[i]) * s_overlap[i];
977 Dune::FieldVector<ctype, dim> upperright,
978 std::array<int, dim> s,
979 std::bitset<dim> periodic = std::bitset<dim>(0ULL),
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)
990 "YaspGrid coordinate container template parameter and given constructor values do not match!");
995 std::fill(o.begin(), o.end(), 0);
996 iTupel o_interior(o);
997 iTupel s_interior(s);
1003 for (
int i=0; i<dim; i++)
1006 int toosmall = (s_interior[i] / 2 <=
overlap) &&
1007 (periodic[i] || (s_interior[i] != s[i]));
1010 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR,
comm);
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.");
1016 #endif // #if HAVE_MPI
1018 iTupel s_overlap(s_interior);
1019 for (
int i=0; i<dim; i++)
1021 if ((o_interior[i] -
overlap > 0) || (periodic[i]))
1023 if ((o_interior[i] + s_interior[i] +
overlap <= _coarseSize[i]) || (periodic[i]))
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];
1048 std::bitset<dim> periodic = std::bitset<dim>(0ULL),
1053 leafIndexSet_(*this), _periodic(periodic), _overlap(
overlap),
1054 keep_ovlp(true), adaptRefCount(0), adaptActive(false)
1057 DUNE_THROW(
Dune::GridError,
"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1061 "YaspGrid coordinate container template parameter and given constructor values do not match!");
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];
1072 std::fill(o.begin(), o.end(), 0);
1073 iTupel o_interior(o);
1074 iTupel s_interior(_coarseSize);
1076 _torus.
partition(_torus.
rank(),o,_coarseSize,o_interior,s_interior);
1080 for (
int i=0; i<dim; i++)
1083 int toosmall = (s_interior[i] / 2 <=
overlap) &&
1084 (periodic[i] || (s_interior[i] != _coarseSize[i]));
1087 MPI_Allreduce(&toosmall, &global, 1, MPI_INT, MPI_LOR,
comm);
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.");
1093 #endif // #if HAVE_MPI
1096 std::array<std::vector<ctype>,dim> newcoords;
1097 std::array<int, dim> offset(o_interior);
1100 for (
int i=0; i<dim; ++i)
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;
1108 if (o_interior[i] -
overlap > 0)
1113 if (o_interior[i] + s_interior[i] +
overlap < _coarseSize[i])
1122 if ((periodic[i]) && (o_interior[i] + s_interior[i] +
overlap >= _coarseSize[i]))
1125 typename std::vector<ctype>::iterator it = coords[i].begin();
1127 newcoords[i].push_back(newcoords[i].back() - *it + *(++it));
1130 if ((periodic[i]) && (o_interior[i] -
overlap <= 0))
1135 typename std::vector<ctype>::iterator it = coords[i].end() - 1;
1137 newcoords[i].insert(newcoords[i].
begin(), newcoords[i].
front() - *it + *(--it));
1164 YaspGrid (std::array<std::vector<ctype>, dim> coords,
1165 std::bitset<dim> periodic,
1168 std::array<int,dim> coarseSize,
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)
1176 "YaspGrid coordinate container template parameter and given constructor values do not match!");
1179 DUNE_THROW(
Dune::GridError,
"Setup of a tensorproduct grid requires monotonous sequences of coordinates.");
1181 for (
int i=0; i<dim; i++)
1182 _L[i] = coords[i][coords[i].
size() - 1] - coords[i][0];
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);
1191 _torus.
partition(_torus.
rank(),o,coarseSize,o_interior,s_interior);
1194 std::array<int,dim> offset(o_interior);
1195 for (
int i=0; i<dim; i++)
1196 if ((periodic[i]) || (o_interior[i] > 0))
1220 return _levels.size()-1;
1228 "Coarsening " << -refCount <<
" levels requested!");
1231 for (
int k=refCount; k<0; k++)
1235 _levels.back() = empty;
1239 indexsets.pop_back();
1243 for (
int k=0; k<refCount; k++)
1246 YGridLevel& cg = _levels[
maxLevel()];
1248 std::bitset<dim> ovlp_low(0ULL), ovlp_up(0ULL);
1249 for (
int i=0; i<dim; i++)
1251 if (cg.overlap[0].dataBegin()->origin(i) > 0 || _periodic[i])
1253 if (cg.overlap[0].dataBegin()->max(i) + 1 <
globalSize(i) || _periodic[i])
1257 Coordinates newcont(cg.coords.refine(ovlp_low, ovlp_up, cg.overlapSize, keep_ovlp));
1259 int overlap = (keep_ovlp) ? 2*cg.overlapSize : cg.overlapSize;
1263 for (
int i=0; i<dim; i++)
1264 o_interior[i] = 2*cg.interior[0].dataBegin()->origin(i);
1267 _levels.resize(_levels.size() + 1);
1280 keep_ovlp = keepPhysicalOverlap;
1294 bool mark(
int refCount,
const typename Traits::template Codim<0>::Entity & e )
1296 assert(adaptActive ==
false);
1297 if (e.level() !=
maxLevel())
return false;
1298 adaptRefCount = std::max(adaptRefCount, refCount);
1308 int getMark (
const typename Traits::template Codim<0>::Entity &e )
const
1310 return ( e.level() ==
maxLevel() ) ? adaptRefCount : 0;
1317 return (adaptRefCount > 0);
1324 adaptRefCount =
comm().max(adaptRefCount);
1325 return (adaptRefCount < 0);
1331 adaptActive =
false;
1336 template<
int cd, PartitionIteratorType pitype>
1337 typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
lbegin (
int level)
const
1339 return levelbegin<cd,pitype>(level);
1343 template<
int cd, PartitionIteratorType pitype>
1344 typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
lend (
int level)
const
1346 return levelend<cd,pitype>(level);
1351 typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator
lbegin (
int level)
const
1353 return levelbegin<cd,All_Partition>(level);
1358 typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator
lend (
int level)
const
1360 return levelend<cd,All_Partition>(level);
1364 template<
int cd, PartitionIteratorType pitype>
1365 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator
leafbegin ()
const
1367 return levelbegin<cd,pitype>(
maxLevel());
1371 template<
int cd, PartitionIteratorType pitype>
1372 typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator
leafend ()
const
1374 return levelend<cd,pitype>(
maxLevel());
1379 typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator
leafbegin ()
const
1381 return levelbegin<cd,All_Partition>(
maxLevel());
1386 typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator
leafend ()
const
1388 return levelend<cd,All_Partition>(
maxLevel());
1392 template <
typename Seed>
1393 typename Traits::template Codim<Seed::codimension>::Entity
1396 const int codim = Seed::codimension;
1399 typedef typename Traits::template Codim<Seed::codimension>::Entity
Entity;
1403 return Entity(EntityImp(g,YIterator(g->overlapfront[codim],seed.impl().coord(),seed.impl().offset())));
1410 return g->overlapSize;
1417 return g->overlapSize;
1433 int size (
int level,
int codim)
const
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();
1455 return (type.isCube()) ?
size(level,dim-type.dim()) : 0;
1479 template<
class DataHandleImp,
class DataType>
1482 YaspCommunicateMeta<dim,dim>::comm(*
this,data,iftype,dir,level);
1489 template<
class DataHandleImp,
class DataType>
1492 YaspCommunicateMeta<dim,dim>::comm(*
this,data,iftype,dir,this->
maxLevel());
1499 template<
class DataHandle,
int codim>
1503 if (!data.contains(dim,codim))
return;
1506 typedef typename DataHandle::DataType DataType;
1517 sendlist = &g->send_interiorborder_interiorborder[codim];
1518 recvlist = &g->recv_interiorborder_interiorborder[codim];
1522 sendlist = &g->send_interiorborder_overlapfront[codim];
1523 recvlist = &g->recv_overlapfront_interiorborder[codim];
1527 sendlist = &g->send_overlap_overlapfront[codim];
1528 recvlist = &g->recv_overlapfront_overlap[codim];
1532 sendlist = &g->send_overlapfront_overlapfront[codim];
1533 recvlist = &g->recv_overlapfront_overlapfront[codim];
1538 std::swap(sendlist,recvlist);
1543 std::vector<int> send_size(sendlist->
size(),-1);
1544 std::vector<int> recv_size(recvlist->
size(),-1);
1545 std::vector<size_t*> send_sizes(sendlist->
size(),
static_cast<size_t*
>(0));
1546 std::vector<size_t*> recv_sizes(recvlist->
size(),
static_cast<size_t*
>(0));
1551 if (data.fixedSize(dim,codim))
1555 for (ListIt is=sendlist->
begin(); is!=sendlist->
end(); ++is)
1557 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1559 send_size[cnt] = is->grid.totalsize() * data.size(*it);
1563 for (ListIt is=recvlist->
begin(); is!=recvlist->
end(); ++is)
1565 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1567 recv_size[cnt] = is->grid.totalsize() * data.size(*it);
1575 for (ListIt is=sendlist->
begin(); is!=sendlist->
end(); ++is)
1578 size_t *buf =
new size_t[is->grid.totalsize()];
1579 send_sizes[cnt] = buf;
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
1587 for ( ; it!=itend; ++it)
1589 buf[i] = data.size(*it);
1598 torus().
send(is->rank,buf,is->grid.totalsize()*
sizeof(
size_t));
1604 for (ListIt is=recvlist->
begin(); is!=recvlist->
end(); ++is)
1607 size_t *buf =
new size_t[is->grid.totalsize()];
1608 recv_sizes[cnt] = buf;
1611 torus().
recv(is->rank,buf,is->grid.totalsize()*
sizeof(
size_t));
1620 for (ListIt is=sendlist->
begin(); is!=sendlist->
end(); ++is)
1622 delete[] send_sizes[cnt];
1623 send_sizes[cnt] = 0;
1629 for (ListIt is=recvlist->
begin(); is!=recvlist->
end(); ++is)
1632 size_t *buf = recv_sizes[cnt];
1636 for (
int i=0; i<is->grid.totalsize(); ++i)
1647 std::vector<DataType*> sends(sendlist->
size(),
static_cast<DataType*
>(0));
1649 for (ListIt is=sendlist->
begin(); is!=sendlist->
end(); ++is)
1652 DataType *buf =
new DataType[send_size[cnt]];
1658 MessageBuffer<DataType> mb(buf);
1661 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1663 typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1665 for ( ; it!=itend; ++it)
1666 data.gather(mb,*it);
1669 torus().
send(is->rank,buf,send_size[cnt]*
sizeof(DataType));
1674 std::vector<DataType*> recvs(recvlist->
size(),
static_cast<DataType*
>(0));
1676 for (ListIt is=recvlist->
begin(); is!=recvlist->
end(); ++is)
1679 DataType *buf =
new DataType[recv_size[cnt]];
1685 torus().
recv(is->rank,buf,recv_size[cnt]*
sizeof(DataType));
1694 for (ListIt is=sendlist->
begin(); is!=sendlist->
end(); ++is)
1696 delete[] sends[cnt];
1703 for (ListIt is=recvlist->
begin(); is!=recvlist->
end(); ++is)
1706 DataType *buf = recvs[cnt];
1709 MessageBuffer<DataType> mb(buf);
1712 if (data.fixedSize(dim,codim))
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
1719 for ( ; it!=itend; ++it)
1720 data.scatter(mb,*it,n);
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
1730 for ( ; it!=itend; ++it)
1731 data.scatter(mb,*it,sbuf[i++]);
1744 return theglobalidset;
1749 return theglobalidset;
1754 if (level<0 || level>
maxLevel()) DUNE_THROW(RangeError,
"level out of range");
1755 return *(indexsets[level]);
1760 return leafIndexSet_;
1785 template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1789 class MessageBuffer {
1792 MessageBuffer (DT *p)
1801 void write (
const Y& data)
1803 static_assert(( std::is_same<DT,Y>::value ),
"DataType mismatch");
1809 void read (Y& data)
const
1811 static_assert(( std::is_same<DT,Y>::value ),
"DataType mismatch");
1822 template<
int cd, PartitionIteratorType pitype>
1826 if (level<0 || level>
maxLevel()) DUNE_THROW(RangeError,
"level out of range");
1837 return levelend <cd, pitype> (level);
1839 DUNE_THROW(
GridError,
"YaspLevelIterator with this codim or partition type not implemented");
1843 template<
int cd, PartitionIteratorType pitype>
1844 YaspLevelIterator<cd,pitype,GridImp> levelend (
int level)
const
1847 if (level<0 || level>
maxLevel()) DUNE_THROW(RangeError,
"level out of range");
1850 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interior[cd].end());
1852 return YaspLevelIterator<cd,pitype,GridImp>(g,g->interiorborder[cd].end());
1854 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlap[cd].end());
1856 return YaspLevelIterator<cd,pitype,GridImp>(g,g->overlapfront[cd].end());
1858 DUNE_THROW(GridError,
"YaspLevelIterator with this codim or partition type not implemented");
1863 Torus<CollectiveCommunicationType,dim> _torus;
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;
1869 Dune::FieldVector<ctype, dim> _L;
1871 std::bitset<dim> _periodic;
1873 ReservedVector<YGridLevel,32> _levels;
1882 template <
int d,
class CC>
1885 int rank = grid.torus().rank();
1887 s <<
"[" << rank <<
"]:" <<
" YaspGrid maxlevel=" << grid.maxLevel() << std::endl;
1889 s <<
"Printing the torus: " <<std::endl;
1890 s << grid.torus() << std::endl;
1894 s <<
"[" << rank <<
"]: " << std::endl;
1895 s <<
"[" << rank <<
"]: " <<
"==========================================" << std::endl;
1896 s <<
"[" << rank <<
"]: " <<
"level=" << g->level() << std::endl;
1898 for (
int codim = 0; codim < d + 1; ++codim)
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;
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;
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;
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;
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;
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;
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;
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;
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;
1953 namespace Capabilities
1963 template<
int dim,
class Coordinates>
1966 static const bool v =
true;
1972 template<
int dim,
class Coordinates>
1975 static const bool v =
true;
1976 static const unsigned int topologyId = GeometryTypes::cube(dim).id();
1982 template<
int dim,
class Coordinates>
1985 static const bool v =
true;
1991 template<
int dim,
class Coordinates,
int codim>
1994 static const bool v =
true;
2001 template<
int dim,
class Coordinates,
int codim>
2004 static const bool v =
true;
2010 template<
int dim,
int codim,
class Coordinates>
2013 static const bool v =
true;
2019 template<
int dim,
class Coordinates>
2022 static const bool v =
true;
2028 template<
int dim,
class Coordinates>
2031 static const bool v =
true;