DUNE-FEM (unstable)

yaspgridintersection.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_GRID_YASPGRIDINTERSECTION_HH
6#define DUNE_GRID_YASPGRIDINTERSECTION_HH
7
15namespace Dune {
16
20 template<class GridImp>
22 {
23 constexpr static int dim = GridImp::dimension;
24 constexpr static int dimworld = GridImp::dimensionworld;
25 typedef typename GridImp::ctype ctype;
26
27 typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
28 typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
29
30 friend class YaspIntersectionIterator<GridImp>;
31
32 public:
33 // types used from grids
34 typedef typename GridImp::YGridLevelIterator YGLI;
35 typedef typename GridImp::YGrid::Iterator I;
36 typedef typename GridImp::template Codim<0>::Entity Entity;
37 typedef typename GridImp::template Codim<1>::Geometry Geometry;
38 typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
39
40 void update() {
41
42 // vector with per-direction movements
43 std::array<int,dim> dist{{0}};
44
45 // first move: back to center
46 dist[_dir] = 1 - 2*_face;
47
48 // update face info
49 _dir = _count / 2;
50 _face = _count % 2;
51
52 // second move: to new neighbor
53 dist[_dir] += -1 + 2*_face;
54
55 // move transforming iterator
56 _outside.transformingsubiterator().move(dist);
57 }
58
62 bool boundary () const
63 {
64 // Coordinate of intersection in its direction
65 int coord = _inside.transformingsubiterator().coord(_dir) + _face;
66 if (_inside.gridlevel()->mg->isPeriodic(_dir))
67 return false;
68 else
69 return coord == 0
70 ||
71 coord == _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),_dir);
72 }
73
75 bool neighbor () const
76 {
77 // Coordinate of intersection in its direction
78 int coord = _inside.transformingsubiterator().coord(_dir) + _face;
79 return coord > _inside.gridlevel()->overlap[0].dataBegin()->min(_dir)
80 &&
81 coord <= _inside.gridlevel()->overlap[0].dataBegin()->max(_dir);
82 }
83
85 bool conforming () const
86 {
87 return true;
88 }
89
92 Entity inside() const
93 {
94 return Entity(_inside);
95 }
96
98 Entity outside() const
99 {
100 return Entity(_outside);
101 }
102
106 {
107 if(! boundary())
108 DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
109 // size of local macro grid
110 const std::array<int, dim> & size = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size();
111 const std::array<int, dim> & origin = _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin();
112 std::array<int, dim> sides;
113 {
114 for (int i=0; i<dim; i++)
115 {
116 sides[i] =
117 ((_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i)
118 == 0)+
119 (_inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->origin(i) +
120 _inside.gridlevel()->mg->begin()->overlap[0].dataBegin()->size(i)
121 ==
122 _inside.gridlevel()->mg->levelSize(0,i)));
123
124 }
125 }
126 // global position of the cell on macro grid
127 std::array<int, dim> pos = _inside.transformingsubiterator().coord();
128 for(int i=0; i<dim; i++)
129 {
130 pos[i] = pos[i] / (1<<_inside.level());
131 pos[i] = pos[i] - origin[i];
132 }
133 // compute unit-cube-face-sizes
134 std::array<int, dim> fsize;
135 {
136 int vol = 1;
137 for (int k=0; k<dim; k++)
138 vol *= size[k];
139 for (int k=0; k<dim; k++)
140 fsize[k] = vol/size[k];
141 }
142 // compute index in the unit-cube-face
143 int index = 0;
144 {
145 int localoffset = 1;
146 for (int k=dim-1; k>=0; k--)
147 {
148 if (k == _dir) continue;
149 index += (pos[k]) * localoffset;
150 localoffset *= size[k];
151 }
152 }
153 // add unit-cube-face-offsets
154 {
155 for (int k=0; k<_dir; k++)
156 index += sides[k] * fsize[k];
157 // add fsize if we are on the right face and there is a left-face-boundary on this processor
158 index += _face * (sides[_dir]>1) * fsize[_dir];
159 }
160
161 return index;
162 }
163
166 {
167 return centerUnitOuterNormal();
168 }
169
172 {
173 return centerUnitOuterNormal();
174 }
175
178 {
180 normal[_dir] = (_face==0) ? -1.0 : 1.0;
181 return normal;
182 }
183
188 {
189 return geometry().volume() * centerUnitOuterNormal();
190 }
191
195 LocalGeometry geometryInInside () const
196 {
197 // set of dimensions that span the intersection
198 std::bitset<dim> s;
199 s.set();
200 s[_dir] = false;
201
202 // lower-left and upper-right corners
205
206 ll[_dir] = ur[_dir] = (_face==0) ? 0.0 : 1.0;
207
208 return LocalGeometry(LocalGeometryImpl(ll,ur,s));
209 }
210
214 LocalGeometry geometryInOutside () const
215 {
216 // set of dimensions that span the intersection
217 std::bitset<dim> s;
218 s.set();
219 s[_dir] = false;
220
221 // lower-left and upper-right corners
224
225 ll[_dir] = ur[_dir] = (_face==1) ? 0.0 : 1.0;
226
227 return LocalGeometry(LocalGeometryImpl(ll,ur,s));
228 }
229
232 Geometry geometry () const
233 {
234
235 std::bitset<dim> shift;
236 shift.set();
237 shift[_dir] = false;
238
240 for (int i=0; i<dimworld; i++)
241 {
242 int coord = _inside.transformingsubiterator().coord(i);
243
244 if ((i == _dir) and (_face))
245 coord++;
246
247 ll[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
248 if (i != _dir)
249 coord++;
250 ur[i] = _inside.transformingsubiterator().coordCont()->coordinate(i,coord);
251
252 // If on periodic overlap, transform coordinates by domain size
253 if (_inside.gridlevel()->mg->isPeriodic(i)) {
254 int coordPeriodic = _inside.transformingsubiterator().coord(i);
255 if (coordPeriodic < 0) {
256 auto size = _inside.gridlevel()->mg->domainSize()[i];
257 ll[i] += size;
258 ur[i] += size;
259 } else if (coordPeriodic + 1 > _inside.gridlevel()->mg->levelSize(_inside.gridlevel()->level(),i)) {
260 auto size = _inside.gridlevel()->mg->domainSize()[i];
261 ll[i] -= size;
262 ur[i] -= size;
263 }
264 }
265 }
266
267 GeometryImpl _is_global(ll,ur,shift);
268 return Geometry( _is_global );
269 }
270
273 {
274 return GeometryTypes::cube(dim-1);
275 }
276
278 int indexInInside () const
279 {
280 return _count;
281 }
282
284 int indexInOutside () const
285 {
286 // flip the last bit
287 return _count^1;
288 }
289
291 : _count(~std::uint8_t(0)) // Use as marker for invalid intersection
292 , _dir(0)
293 , _face(0)
294 {}
295
297 YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
298 _inside(myself.gridlevel(),
299 myself.transformingsubiterator()),
300 _outside(myself.gridlevel(),
301 myself.transformingsubiterator()),
302 // initialize to first neighbor
303 _count(0),
304 _dir(0),
305 _face(0)
306 {
307 if (toend)
308 {
309 // initialize end iterator
310 _count = 2*dim;
311 return;
312 }
313 _count = 0;
314
315 // move transforming iterator
316 _outside.transformingsubiterator().move(_dir,-1);
317 }
318
320
322 void assign (const YaspIntersection& it)
323 {
324 *this = it;
325 }
326
327 bool equals(const YaspIntersection& other) const
328 {
329 // compare counts first -- that's cheaper if the test fails
330 return _count == other._count && _inside.equals(other._inside);
331 }
332
333 private:
334 /* The two entities that make up the intersection */
335 YaspEntity<0,GridImp::dimension,GridImp> _inside;
336 YaspEntity<0,GridImp::dimension,GridImp> _outside;
337 /* current position */
338 std::uint8_t _count;
339 std::uint8_t _dir;
340 std::uint8_t _face;
341 };
342} // namespace Dune
343
344#endif // DUNE_GRID_YASPGRIDINTERSECTION_HH
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:20
YaspIntersectionIterator enables iteration over intersections with neighboring codim 0 entities.
Definition: yaspgridintersectioniterator.hh:22
YaspIntersection provides data about intersection with neighboring codim 0 entities.
Definition: yaspgridintersection.hh:22
FieldVector< ctype, dimworld > centerUnitOuterNormal() const
return unit outer normal at center of intersection geometry
Definition: yaspgridintersection.hh:177
Entity inside() const
Definition: yaspgridintersection.hh:92
Geometry geometry() const
Definition: yaspgridintersection.hh:232
int indexInOutside() const
local index of codim 1 entity in neighbor where intersection is contained in
Definition: yaspgridintersection.hh:284
LocalGeometry geometryInInside() const
Definition: yaspgridintersection.hh:195
int boundarySegmentIndex() const
Definition: yaspgridintersection.hh:105
bool conforming() const
Yasp is always conform.
Definition: yaspgridintersection.hh:85
bool neighbor() const
return true if neighbor across intersection exists in this processor
Definition: yaspgridintersection.hh:75
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:165
YaspIntersection(const YaspEntity< 0, dim, GridImp > &myself, bool toend)
make intersection iterator from entity, initialize to first neighbor
Definition: yaspgridintersection.hh:297
GeometryType type() const
obtain the type of reference element for this intersection
Definition: yaspgridintersection.hh:272
void assign(const YaspIntersection &it)
copy constructor – use default
Definition: yaspgridintersection.hh:322
int indexInInside() const
local index of codim 1 entity in self where intersection is contained in
Definition: yaspgridintersection.hh:278
FieldVector< ctype, dimworld > integrationOuterNormal(const FieldVector< ctype, dim-1 > &local) const
Definition: yaspgridintersection.hh:187
LocalGeometry geometryInOutside() const
Definition: yaspgridintersection.hh:214
Entity outside() const
return Entity on the outside of this intersection
Definition: yaspgridintersection.hh:98
bool boundary() const
Definition: yaspgridintersection.hh:62
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:171
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:587
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
STL namespace.
Static tag representing a codimension.
Definition: dimension.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)