Dune Core Modules (2.5.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
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<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
101#if DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS
104 int boundaryId() const
105 {
106 if(boundary()) return indexInInside()+1;
107 return 0;
108 }
109#endif
110
114 {
115 if(! boundary())
116 DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
117 // size of local macro grid
118 const std::array<int, dim> & size = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size();
119 const std::array<int, dim> & origin = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin();
120 std::array<int, dim> sides;
121 {
122 for (int i=0; i<dim; i++)
123 {
124 sides[i] =
125 ((_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i)
126 == 0)+
127 (_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i) +
128 _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size(i)
129 ==
130 _inside.gridlevel()->mg->levelSize(0,i)));
131
132 }
133 }
134 // global position of the cell on macro grid
135 std::array<int, dim> pos = _inside.transformingsubiterator().coord();
136 for(int i=0; i<dim; i++)
137 {
138 pos[i] = pos[i] / (1<<_inside.level());
139 pos[i] = pos[i] - origin[i];
140 }
141 // compute unit-cube-face-sizes
142 std::array<int, dim> fsize;
143 {
144 int vol = 1;
145 for (int k=0; k<dim; k++)
146 vol *= size[k];
147 for (int k=0; k<dim; k++)
148 fsize[k] = vol/size[k];
149 }
150 // compute index in the unit-cube-face
151 int index = 0;
152 {
153 int localoffset = 1;
154 for (int k=dim-1; k>=0; k--)
155 {
156 if (k == _dir) continue;
157 index += (pos[k]) * localoffset;
158 localoffset *= size[k];
159 }
160 }
161 // add unit-cube-face-offsets
162 {
163 for (int k=0; k<_dir; k++)
164 index += sides[k] * fsize[k];
165 // add fsize if we are on the right face and there is a left-face-boundary on this processor
166 index += _face * (sides[_dir]>1) * fsize[_dir];
167 }
168
169 return index;
170 }
171
174 {
175 return _faceInfo[_count].normal;
176 }
177
180 {
181 return _faceInfo[_count].normal;
182 }
183
186 {
187 return _faceInfo[_count].normal;
188 }
189
194 {
195 FieldVector<ctype, dimworld> n = _faceInfo[_count].normal;
196 n *= geometry().volume();
197 return n;
198 }
199
203 LocalGeometry geometryInInside () const
204 {
205 return LocalGeometry( _faceInfo[_count].geom_inside );
206 }
207
211 LocalGeometry geometryInOutside () const
212 {
213 return LocalGeometry( _faceInfo[_count].geom_outside );
214 }
215
218 Geometry geometry () const
219 {
220
221 std::bitset<dim> shift;
222 shift.set();
223 shift[_dir] = false;
224
226 for (int i=0; i<dimworld; i++)
227 {
228 int coord = _inside.transformingsubiterator().coord(i);
229
230 if ((i == _dir) and (_face))
231 coord++;
232
233 ll[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
234 if (i != _dir)
235 coord++;
236 ur[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
237
238 // If on periodic overlap, transform coordinates by domain size
239 if (_inside.gridlevel()->mg->isPeriodic(i)) {
240 int coordPeriodic = _inside.transformingsubiterator().coord(i);
241 if (coordPeriodic < 0) {
242 auto size = _inside.gridlevel()->mg->domainSize()[i];
243 ll[i] += size;
244 ur[i] += size;
245 } else if (coordPeriodic + 1 > _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),i)) {
246 auto size = _inside.gridlevel()->mg->domainSize()[i];
247 ll[i] -= size;
248 ur[i] -= size;
249 }
250 }
251 }
252
253 GeometryImpl _is_global(ll,ur,shift);
254 return Geometry( _is_global );
255 }
256
259 {
260 return GeometryType(GeometryType::cube, dim-1);
261 }
262
264 int indexInInside () const
265 {
266 return _count;
267 }
268
270 int indexInOutside () const
271 {
272 // flip the last bit
273 return _count^1;
274 }
275
277 : _count(~uint8_t(0)) // Use as marker for invalid intersection
278 , _dir(0)
279 , _face(0)
280 {}
281
283 YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
284 _inside(myself.gridlevel(),
285 myself.transformingsubiterator()),
286 _outside(myself.gridlevel(),
287 myself.transformingsubiterator()),
288 // initialize to first neighbor
289 _count(0),
290 _dir(0),
291 _face(0)
292 {
293 if (toend)
294 {
295 // initialize end iterator
296 _count = 2*dim;
297 return;
298 }
299 _count = 0;
300
301 // move transforming iterator
302 _outside.transformingsubiterator().move(_dir,-1);
303 }
304
306
308 void assign (const YaspIntersection& it)
309 {
310 *this = it;
311 }
312
313 bool equals(const YaspIntersection& other) const
314 {
315 // compare counts first -- that's cheaper if the test fails
316 return _count == other._count && _inside.equals(other._inside);
317 }
318
319 private:
320 /* EntityPointers (get automatically updated) */
321 YaspEntity<0,GridImp::dimension,GridImp> _inside;
322 YaspEntity<0,GridImp::dimension,GridImp> _outside;
323 /* current position */
324 uint8_t _count;
325 uint8_t _dir;
326 uint8_t _face;
327
328 /* static data */
329 struct faceInfo
330 {
331 FieldVector<ctype, dimworld> normal;
332 LocalGeometryImpl geom_inside;
333 LocalGeometryImpl geom_outside;
334 };
335
336 /* static face info */
337 static const std::array<faceInfo, 2*GridImp::dimension> _faceInfo;
338
339 static std::array<faceInfo, 2*dim> initFaceInfo()
340 {
341 std::array<faceInfo, 2*dim> I;
342 for (uint8_t i=0; i<dim; i++)
343 {
344 // compute normals
345 I[2*i].normal = 0.0;
346 I[2*i+1].normal = 0.0;
347 I[2*i].normal[i] = -1.0;
348 I[2*i+1].normal[i] = +1.0;
349
350 // determine the shift vector for these intersection
351 std::bitset<dim> s;
352 s.set();
353 s[i] = false;
354
355 // store intersection geometries
358 ur[i] = 0.0;
359
360 I[2*i].geom_inside = LocalGeometryImpl(ll,ur,s);
361 I[2*i+1].geom_outside = LocalGeometryImpl(ll,ur,s);
362
363 ll[i] = 1.0;
364 ur[i] = 1.0;
365
366 I[2*i].geom_outside = LocalGeometryImpl(ll,ur,s);
367 I[2*i+1].geom_inside = LocalGeometryImpl(ll,ur,s);
368 }
369
370 return I;
371 }
372 };
373
374 template<class GridImp>
375 const std::array<typename YaspIntersection<GridImp>::faceInfo, 2*GridImp::dimension>
376 YaspIntersection<GridImp>::_faceInfo =
377 YaspIntersection<GridImp>::initFaceInfo();
378
379} // namespace Dune
380
381#endif // DUNE_GRID_YASPGRIDINTERSECTION_HH
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:268
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:274
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:185
Entity inside() const
Definition: yaspgridintersection.hh:90
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:179
Geometry geometry() const
Definition: yaspgridintersection.hh:218
int indexInOutside() const
local index of codim 1 entity in neighbor where intersection is contained in
Definition: yaspgridintersection.hh:270
LocalGeometry geometryInInside() const
Definition: yaspgridintersection.hh:203
int boundarySegmentIndex() const
Definition: yaspgridintersection.hh:113
bool conforming() const
Yasp is always conform.
Definition: yaspgridintersection.hh:83
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:283
GeometryType type() const
obtain the type of reference element for this intersection
Definition: yaspgridintersection.hh:258
void assign(const YaspIntersection &it)
copy constructor – use default
Definition: yaspgridintersection.hh:308
int indexInInside() const
local index of codim 1 entity in self where intersection is contained in
Definition: yaspgridintersection.hh:264
FieldVector< ctype, dimworld > integrationOuterNormal(const FieldVector< ctype, dim-1 > &local) const
Definition: yaspgridintersection.hh:193
LocalGeometry geometryInOutside() const
Definition: yaspgridintersection.hh:211
Entity outside() const
return Entity on the outside of this intersection
Definition: yaspgridintersection.hh:96
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:173
bool boundary() const
Definition: yaspgridintersection.hh:60
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:441
Dune namespace.
Definition: alignment.hh:11
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 23, 23:29, 2024)