dune-grid  2.8.0
yaspgridintersection.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_YASPGRIDINTERSECTION_HH
4 #define DUNE_GRID_YASPGRIDINTERSECTION_HH
5 
13 namespace Dune {
14 
18  template<class GridImp>
20  {
21  enum { dim=GridImp::dimension };
22  enum { dimworld=GridImp::dimensionworld };
23  typedef typename GridImp::ctype ctype;
24 
25  typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
26  typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
27 
28  friend class YaspIntersectionIterator<GridImp>;
29 
30  public:
31  // types used from grids
32  typedef typename GridImp::YGridLevelIterator YGLI;
33  typedef typename GridImp::YGrid::Iterator I;
34  typedef typename GridImp::template Codim<0>::Entity Entity;
35  typedef typename GridImp::template Codim<1>::Geometry Geometry;
36  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
37 
38  void update() {
39 
40  // vector with per-direction movements
41  std::array<int,dim> dist{{0}};
42 
43  // first move: back to center
44  dist[_dir] = 1 - 2*_face;
45 
46  // update face info
47  _dir = _count / 2;
48  _face = _count % 2;
49 
50  // second move: to new neighbor
51  dist[_dir] += -1 + 2*_face;
52 
53  // move transforming iterator
54  _outside.transformingsubiterator().move(dist);
55  }
56 
60  bool boundary () const
61  {
62  // Coordinate of intersection in its direction
63  int coord = _inside.transformingsubiterator().coord(_dir) + _face;
64  if (_inside.gridlevel()->mg->isPeriodic(_dir))
65  return false;
66  else
67  return coord == 0
68  ||
69  coord == _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),_dir);
70  }
71 
73  bool neighbor () const
74  {
75  // Coordinate of intersection in its direction
76  int coord = _inside.transformingsubiterator().coord(_dir) + _face;
77  return coord > _inside.gridlevel()->overlap[0].dataBegin()->min(_dir)
78  &&
79  coord <= _inside.gridlevel()->overlap[0].dataBegin()->max(_dir);
80  }
81 
83  bool conforming () const
84  {
85  return true;
86  }
87 
90  Entity inside() const
91  {
92  return Entity(_inside);
93  }
94 
96  Entity outside() const
97  {
98  return Entity(_outside);
99  }
100 
104  {
105  if(! boundary())
106  DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
107  // size of local macro grid
108  const std::array<int, dim> & size = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size();
109  const std::array<int, dim> & origin = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin();
110  std::array<int, dim> sides;
111  {
112  for (int i=0; i<dim; i++)
113  {
114  sides[i] =
115  ((_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i)
116  == 0)+
117  (_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i) +
118  _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size(i)
119  ==
120  _inside.gridlevel()->mg->levelSize(0,i)));
121 
122  }
123  }
124  // global position of the cell on macro grid
125  std::array<int, dim> pos = _inside.transformingsubiterator().coord();
126  for(int i=0; i<dim; i++)
127  {
128  pos[i] = pos[i] / (1<<_inside.level());
129  pos[i] = pos[i] - origin[i];
130  }
131  // compute unit-cube-face-sizes
132  std::array<int, dim> fsize;
133  {
134  int vol = 1;
135  for (int k=0; k<dim; k++)
136  vol *= size[k];
137  for (int k=0; k<dim; k++)
138  fsize[k] = vol/size[k];
139  }
140  // compute index in the unit-cube-face
141  int index = 0;
142  {
143  int localoffset = 1;
144  for (int k=dim-1; k>=0; k--)
145  {
146  if (k == _dir) continue;
147  index += (pos[k]) * localoffset;
148  localoffset *= size[k];
149  }
150  }
151  // add unit-cube-face-offsets
152  {
153  for (int k=0; k<_dir; k++)
154  index += sides[k] * fsize[k];
155  // add fsize if we are on the right face and there is a left-face-boundary on this processor
156  index += _face * (sides[_dir]>1) * fsize[_dir];
157  }
158 
159  return index;
160  }
161 
163  FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dim-1>& /* local */) const
164  {
165  return centerUnitOuterNormal();
166  }
167 
169  FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dim-1>& /* local */) const
170  {
171  return centerUnitOuterNormal();
172  }
173 
175  FieldVector<ctype, dimworld> centerUnitOuterNormal () const
176  {
177  FieldVector<ctype, dimworld> normal(0);
178  normal[_dir] = (_face==0) ? -1.0 : 1.0;
179  return normal;
180  }
181 
185  FieldVector<ctype, dimworld> integrationOuterNormal ([[maybe_unused]] const FieldVector<ctype, dim-1>& local) const
186  {
187  return geometry().volume() * centerUnitOuterNormal();
188  }
189 
194  {
195  // set of dimensions that span the intersection
196  std::bitset<dim> s;
197  s.set();
198  s[_dir] = false;
199 
200  // lower-left and upper-right corners
201  Dune::FieldVector<ctype, dim> ll(0.0);
202  Dune::FieldVector<ctype, dim> ur(1.0);
203 
204  ll[_dir] = ur[_dir] = (_face==0) ? 0.0 : 1.0;
205 
206  return LocalGeometry(LocalGeometryImpl(ll,ur,s));
207  }
208 
213  {
214  // set of dimensions that span the intersection
215  std::bitset<dim> s;
216  s.set();
217  s[_dir] = false;
218 
219  // lower-left and upper-right corners
220  Dune::FieldVector<ctype, dim> ll(0.0);
221  Dune::FieldVector<ctype, dim> ur(1.0);
222 
223  ll[_dir] = ur[_dir] = (_face==1) ? 0.0 : 1.0;
224 
225  return LocalGeometry(LocalGeometryImpl(ll,ur,s));
226  }
227 
231  {
232 
233  std::bitset<dim> shift;
234  shift.set();
235  shift[_dir] = false;
236 
237  Dune::FieldVector<ctype,dimworld> ll, ur;
238  for (int i=0; i<dimworld; i++)
239  {
240  int coord = _inside.transformingsubiterator().coord(i);
241 
242  if ((i == _dir) and (_face))
243  coord++;
244 
245  ll[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
246  if (i != _dir)
247  coord++;
248  ur[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
249 
250  // If on periodic overlap, transform coordinates by domain size
251  if (_inside.gridlevel()->mg->isPeriodic(i)) {
252  int coordPeriodic = _inside.transformingsubiterator().coord(i);
253  if (coordPeriodic < 0) {
254  auto size = _inside.gridlevel()->mg->domainSize()[i];
255  ll[i] += size;
256  ur[i] += size;
257  } else if (coordPeriodic + 1 > _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),i)) {
258  auto size = _inside.gridlevel()->mg->domainSize()[i];
259  ll[i] -= size;
260  ur[i] -= size;
261  }
262  }
263  }
264 
265  GeometryImpl _is_global(ll,ur,shift);
266  return Geometry( _is_global );
267  }
268 
271  {
272  return GeometryTypes::cube(dim-1);
273  }
274 
276  int indexInInside () const
277  {
278  return _count;
279  }
280 
282  int indexInOutside () const
283  {
284  // flip the last bit
285  return _count^1;
286  }
287 
289  : _count(~std::uint8_t(0)) // Use as marker for invalid intersection
290  , _dir(0)
291  , _face(0)
292  {}
293 
295  YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
296  _inside(myself.gridlevel(),
297  myself.transformingsubiterator()),
298  _outside(myself.gridlevel(),
299  myself.transformingsubiterator()),
300  // initialize to first neighbor
301  _count(0),
302  _dir(0),
303  _face(0)
304  {
305  if (toend)
306  {
307  // initialize end iterator
308  _count = 2*dim;
309  return;
310  }
311  _count = 0;
312 
313  // move transforming iterator
314  _outside.transformingsubiterator().move(_dir,-1);
315  }
316 
318 
320  void assign (const YaspIntersection& it)
321  {
322  *this = it;
323  }
324 
325  bool equals(const YaspIntersection& other) const
326  {
327  // compare counts first -- that's cheaper if the test fails
328  return _count == other._count && _inside.equals(other._inside);
329  }
330 
331  private:
332  /* The two entities that make up the intersection */
335  /* current position */
336  std::uint8_t _count;
337  std::uint8_t _dir;
338  std::uint8_t _face;
339  };
340 } // namespace Dune
341 
342 #endif // DUNE_GRID_YASPGRIDINTERSECTION_HH
Include standard header files.
Definition: agrid.hh:58
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:130
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
int level() const
level of this element
Definition: yaspgridentity.hh:276
bool equals(const YaspEntity &e) const
Return true when two iterators over the same grid are equal (!).
Definition: yaspgridentity.hh:333
const I & transformingsubiterator() const
Definition: yaspgridentity.hh:405
const YGLI & gridlevel() const
Definition: yaspgridentity.hh:406
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
FieldVector< ctype, dimworld > unitOuterNormal(const FieldVector< ctype, dim-1 > &) const
return unit outer normal, this should be dependent on local coordinates for higher order boundary
Definition: yaspgridintersection.hh:169
bool equals(const YaspIntersection &other) const
Definition: yaspgridintersection.hh:325
Entity inside() const
Definition: yaspgridintersection.hh:90
Geometry geometry() const
Definition: yaspgridintersection.hh:230
GridImp::template Codim< 1 >::LocalGeometry LocalGeometry
Definition: yaspgridintersection.hh:36
YaspIntersection()
Definition: yaspgridintersection.hh:288
int indexInOutside() const
local index of codim 1 entity in neighbor where intersection is contained in
Definition: yaspgridintersection.hh:282
LocalGeometry geometryInInside() const
Definition: yaspgridintersection.hh:193
int boundarySegmentIndex() const
Definition: yaspgridintersection.hh:103
GridImp::template Codim< 1 >::Geometry Geometry
Definition: yaspgridintersection.hh:35
bool conforming() const
Yasp is always conform.
Definition: yaspgridintersection.hh:83
GridImp::template Codim< 0 >::Entity Entity
Definition: yaspgridintersection.hh:34
bool neighbor() const
return true if neighbor across intersection exists in this processor
Definition: yaspgridintersection.hh:73
YaspIntersection(const YaspEntity< 0, dim, GridImp > &myself, bool toend)
make intersection iterator from entity, initialize to first neighbor
Definition: yaspgridintersection.hh:295
GeometryType type() const
obtain the type of reference element for this intersection
Definition: yaspgridintersection.hh:270
void assign(const YaspIntersection &it)
copy constructor – use default
Definition: yaspgridintersection.hh:320
int indexInInside() const
local index of codim 1 entity in self where intersection is contained in
Definition: yaspgridintersection.hh:276
GridImp::YGrid::Iterator I
Definition: yaspgridintersection.hh:33
FieldVector< ctype, dimworld > integrationOuterNormal([[maybe_unused]] const FieldVector< ctype, dim-1 > &local) const
Definition: yaspgridintersection.hh:185
LocalGeometry geometryInOutside() const
Definition: yaspgridintersection.hh:212
GridImp::YGridLevelIterator YGLI
Definition: yaspgridintersection.hh:32
Entity outside() const
return Entity on the outside of this intersection
Definition: yaspgridintersection.hh:96
FieldVector< ctype, dimworld > outerNormal(const FieldVector< ctype, dim-1 > &) const
return unit outer normal, this should be dependent on local coordinates for higher order boundary
Definition: yaspgridintersection.hh:163
bool boundary() const
Definition: yaspgridintersection.hh:60
FieldVector< ctype, dimworld > centerUnitOuterNormal() const
return unit outer normal at center of intersection geometry
Definition: yaspgridintersection.hh:175
void update()
Definition: yaspgridintersection.hh:38
Definition: yaspgridentity.hh:420