Dune Core Modules (2.4.1)

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
13namespace 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<0>::EntityPointer EntityPointer;
36 typedef typename GridImp::template Codim<1>::Geometry Geometry;
37 typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
38
39 void update() {
40
41 // vector with per-direction movements
42 std::array<int,dim> dist{{0}};
43
44 // first move: back to center
45 dist[_dir] = 1 - 2*_face;
46
47 // update face info
48 _dir = _count / 2;
49 _face = _count % 2;
50
51 // second move: to new neighbor
52 dist[_dir] += -1 + 2*_face;
53
54 // move transforming iterator
55 _outside.transformingsubiterator().move(dist);
56 }
57
61 bool boundary () const
62 {
63 // Coordinate of intersection in its direction
64 int coord = _inside.transformingsubiterator().coord(_dir) + _face;
65 if (_inside.gridlevel()->mg->isPeriodic(_dir))
66 return false;
67 else
68 return coord == 0
69 ||
70 coord == _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),_dir);
71 }
72
74 bool neighbor () const
75 {
76 // Coordinate of intersection in its direction
77 int coord = _inside.transformingsubiterator().coord(_dir) + _face;
78 return coord > _inside.gridlevel()->overlap[0].dataBegin()->min(_dir)
79 &&
80 coord <= _inside.gridlevel()->overlap[0].dataBegin()->max(_dir);
81 }
82
84 bool conforming () const
85 {
86 return true;
87 }
88
91 Entity inside() const
92 {
93 return Entity(_inside);
94 }
95
97 Entity outside() const
98 {
99 return Entity(_outside);
100 }
101
102#if DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS
105 int boundaryId() const
106 {
107 if(boundary()) return indexInInside()+1;
108 return 0;
109 }
110#endif
111
115 {
116 if(! boundary())
117 DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
118 // size of local macro grid
119 const Dune::array<int, dim> & size = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size();
120 const Dune::array<int, dim> & origin = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin();
121 Dune::array<int, dim> sides;
122 {
123 for (int i=0; i<dim; i++)
124 {
125 sides[i] =
126 ((_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i)
127 == 0)+
128 (_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i) +
129 _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size(i)
130 ==
131 _inside.gridlevel()->mg->levelSize(0,i)));
132
133 }
134 }
135 // global position of the cell on macro grid
136 Dune::array<int, dim> pos = _inside.transformingsubiterator().coord();
137 for(int i=0; i<dim; i++)
138 {
139 pos[i] = pos[i] / (1<<_inside.level());
140 pos[i] = pos[i] - origin[i];
141 }
142 // compute unit-cube-face-sizes
143 Dune::array<int, dim> fsize;
144 {
145 int vol = 1;
146 for (int k=0; k<dim; k++)
147 vol *= size[k];
148 for (int k=0; k<dim; k++)
149 fsize[k] = vol/size[k];
150 }
151 // compute index in the unit-cube-face
152 int index = 0;
153 {
154 int localoffset = 1;
155 for (int k=dim-1; k>=0; k--)
156 {
157 if (k == _dir) continue;
158 index += (pos[k]) * localoffset;
159 localoffset *= size[k];
160 }
161 }
162 // add unit-cube-face-offsets
163 {
164 for (int k=0; k<_dir; k++)
165 index += sides[k] * fsize[k];
166 // add fsize if we are on the right face and there is a left-face-boundary on this processor
167 index += _face * (sides[_dir]>1) * fsize[_dir];
168 }
169
170 return index;
171 }
172
175 {
176 return _faceInfo[_count].normal;
177 }
178
181 {
182 return _faceInfo[_count].normal;
183 }
184
187 {
188 return _faceInfo[_count].normal;
189 }
190
195 {
196 FieldVector<ctype, dimworld> n = _faceInfo[_count].normal;
197 n *= geometry().volume();
198 return n;
199 }
200
204 LocalGeometry geometryInInside () const
205 {
206 return LocalGeometry( _faceInfo[_count].geom_inside );
207 }
208
212 LocalGeometry geometryInOutside () const
213 {
214 return LocalGeometry( _faceInfo[_count].geom_outside );
215 }
216
219 Geometry geometry () const
220 {
221
222 std::bitset<dim> shift;
223 shift.set();
224 shift[_dir] = false;
225
227 for (int i=0; i<dimworld; i++)
228 {
229 int coord = _inside.transformingsubiterator().coord(i);
230
231 if ((i == _dir) and (_face))
232 coord++;
233
234 ll[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
235 if (i != _dir)
236 coord++;
237 ur[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
238
239 // If on periodic overlap, transform coordinates by domain size
240 if (_inside.gridlevel()->mg->isPeriodic(i)) {
241 int coord = _inside.transformingsubiterator().coord(i);
242 if (coord < 0) {
243 auto size = _inside.gridlevel()->mg->domainSize()[i];
244 ll[i] += size;
245 ur[i] += size;
246 } else if (coord + 1 > _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),i)) {
247 auto size = _inside.gridlevel()->mg->domainSize()[i];
248 ll[i] -= size;
249 ur[i] -= size;
250 }
251 }
252 }
253
254 GeometryImpl _is_global(ll,ur,shift);
255 return Geometry( _is_global );
256 }
257
260 {
261 return GeometryType(GeometryType::cube, dim-1);
262 }
263
265 int indexInInside () const
266 {
267 return _count;
268 }
269
271 int indexInOutside () const
272 {
273 // flip the last bit
274 return _count^1;
275 }
276
278 : _count(~uint8_t(0)) // Use as marker for invalid intersection
279 , _dir(0)
280 , _face(0)
281 {}
282
284 YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
285 _inside(myself.gridlevel(),
286 myself.transformingsubiterator()),
287 _outside(myself.gridlevel(),
288 myself.transformingsubiterator()),
289 // initialize to first neighbor
290 _count(0),
291 _dir(0),
292 _face(0)
293 {
294 if (toend)
295 {
296 // initialize end iterator
297 _count = 2*dim;
298 return;
299 }
300 _count = 0;
301
302 // move transforming iterator
303 _outside.transformingsubiterator().move(_dir,-1);
304 }
305
307
309 void assign (const YaspIntersection& it)
310 {
311 *this = it;
312 }
313
314 bool equals(const YaspIntersection& other) const
315 {
316 // compare counts first -- that's cheaper if the test fails
317 return _count == other._count && _inside.equals(other._inside);
318 }
319
320 private:
321 /* EntityPointers (get automatically updated) */
322 YaspEntity<0,GridImp::dimension,GridImp> _inside;
323 YaspEntity<0,GridImp::dimension,GridImp> _outside;
324 /* current position */
325 uint8_t _count;
326 uint8_t _dir;
327 uint8_t _face;
328
329 /* static data */
330 struct faceInfo
331 {
332 FieldVector<ctype, dimworld> normal;
333 LocalGeometryImpl geom_inside;
334 LocalGeometryImpl geom_outside;
335 };
336
337 /* static face info */
338 static const array<faceInfo, 2*GridImp::dimension> _faceInfo;
339
340 static array<faceInfo, 2*dim> initFaceInfo()
341 {
342 array<faceInfo, 2*dim> I;
343 for (uint8_t i=0; i<dim; i++)
344 {
345 // compute normals
346 I[2*i].normal = 0.0;
347 I[2*i+1].normal = 0.0;
348 I[2*i].normal[i] = -1.0;
349 I[2*i+1].normal[i] = +1.0;
350
351 // determine the shift vector for these intersection
352 std::bitset<dim> s;
353 s.set();
354 s[i] = false;
355
356 // store intersection geometries
359 ur[i] = 0.0;
360
361 I[2*i].geom_inside = LocalGeometryImpl(ll,ur,s);
362 I[2*i+1].geom_outside = LocalGeometryImpl(ll,ur,s);
363
364 ll[i] = 1.0;
365 ur[i] = 1.0;
366
367 I[2*i].geom_outside = LocalGeometryImpl(ll,ur,s);
368 I[2*i+1].geom_inside = LocalGeometryImpl(ll,ur,s);
369 }
370
371 return I;
372 }
373 };
374
375 template<class GridImp>
376 const array<typename YaspIntersection<GridImp>::faceInfo, 2*GridImp::dimension>
377 YaspIntersection<GridImp>::_faceInfo =
378 YaspIntersection<GridImp>::initFaceInfo();
379
380} // namespace Dune
381
382#endif // DUNE_GRID_YASPGRIDINTERSECTION_HH
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:31
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
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 > centerUnitOuterNormal() const
return unit outer normal at center of intersection geometry
Definition: yaspgridintersection.hh:186
Entity inside() const
Definition: yaspgridintersection.hh:91
FieldVector< ctype, dimworld > unitOuterNormal(const FieldVector< ctype, dim-1 > &local) const
return unit outer normal, this should be dependent on local coordinates for higher order boundary
Definition: yaspgridintersection.hh:180
Geometry geometry() const
Definition: yaspgridintersection.hh:219
int indexInOutside() const
local index of codim 1 entity in neighbor where intersection is contained in
Definition: yaspgridintersection.hh:271
LocalGeometry geometryInInside() const
Definition: yaspgridintersection.hh:204
int boundarySegmentIndex() const
Definition: yaspgridintersection.hh:114
bool conforming() const
Yasp is always conform.
Definition: yaspgridintersection.hh:84
bool neighbor() const
return true if neighbor across intersection exists in this processor
Definition: yaspgridintersection.hh:74
YaspIntersection(const YaspEntity< 0, dim, GridImp > &myself, bool toend)
make intersection iterator from entity, initialize to first neighbor
Definition: yaspgridintersection.hh:284
GeometryType type() const
obtain the type of reference element for this intersection
Definition: yaspgridintersection.hh:259
void assign(const YaspIntersection &it)
copy constructor – use default
Definition: yaspgridintersection.hh:309
int indexInInside() const
local index of codim 1 entity in self where intersection is contained in
Definition: yaspgridintersection.hh:265
FieldVector< ctype, dimworld > integrationOuterNormal(const FieldVector< ctype, dim-1 > &local) const
Definition: yaspgridintersection.hh:194
LocalGeometry geometryInOutside() const
Definition: yaspgridintersection.hh:212
Entity outside() const
return EntityPointer to the Entity on the outside of this intersection
Definition: yaspgridintersection.hh:97
FieldVector< ctype, dimworld > outerNormal(const FieldVector< ctype, dim-1 > &local) const
return unit outer normal, this should be dependent on local coordinates for higher order boundary
Definition: yaspgridintersection.hh:174
bool boundary() const
Definition: yaspgridintersection.hh:61
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune namespace.
Definition: alignment.hh:10
Static tag representing a codimension.
Definition: dimension.hh:22
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)