« back to all changes in this revision
Viewing changes to dune/grid/yaspgrid.hh
- browse files at revision 8
- compare with another revision
- download diff
- download tarball
- view history from revision 8
- Committer: Package Import Robot
- Author(s): Ansgar Burchardt
- Date: 2014-02-14 10:49:16 UTC
- mfrom: (1.3.1) (5.1.4 experimental)
- Revision ID: package-import@ubuntu.com-20140214104916-2clchnny3eu7ks4w
Tags: 2.3.0-2
Install /usr/share/dune-grid.
- files added:
- .gitignore
- cmake
- cmake/modules
- files removed:
- .pc/install-grids-as-doc.patch
- .pc/install-grids-as-doc.patch/doc
- .pc/install-grids-as-doc.patch/doc/grids
- .pc/install-grids-as-doc.patch/doc/grids/amc
- .pc/install-grids-as-doc.patch/doc/grids/amiramesh
- .pc/install-grids-as-doc.patch/doc/grids/dgf
- .pc/install-grids-as-doc.patch/doc/grids/gmsh
- .pc/install-grids-as-doc.patch/doc/grids/gridfactory
- .pc/install-grids-as-doc.patch/doc/grids/starcd
- doc/refinement
- dune/grid/common/quadraturerules
- dune/grid/common/refinement
- dune/grid/genericgeometry
- dune/grid/genericgeometry/test
- files modified:
- .pc/applied-patches
added
removed
dune/grid/yaspgrid.hh
1
#ifndef DUNE_YASPGRID_HH
2
#define DUNE_YASPGRID_HH
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_YASPGRID_HH
4
#define DUNE_GRID_YASPGRID_HH
3
5
4
#include<iostream>
5
#include<vector>
6
#include<algorithm>
7
#include<stack>
6
#include <iostream>
7
#include <vector>
8
#include <algorithm>
9
#include <stack>
8
10
9
11
// either include stdint.h or provide fallback for uint8_t
10
12
#if HAVE_STDINT_H
11
#include<stdint.h>
13
#include <stdint.h>
12
14
#else
13
15
typedef unsigned char uint8_t;
14
16
#endif
16
18
#include <dune/grid/common/grid.hh> // the grid base classes
17
19
#include <dune/grid/yaspgrid/grids.hh> // the yaspgrid base classes
18
20
#include <dune/grid/common/capabilities.hh> // the capabilities
19
#include <dune/common/misc.hh>
21
#include <dune/common/shared_ptr.hh>
20
22
#include <dune/common/bigunsignedint.hh>
21
23
#include <dune/common/typetraits.hh>
22
#include <dune/common/collectivecommunication.hh>
23
#include <dune/common/mpihelper.hh>
24
#include <dune/common/reservedvector.hh>
25
#include <dune/common/parallel/collectivecommunication.hh>
26
#include <dune/common/parallel/mpihelper.hh>
24
27
#include <dune/geometry/genericgeometry/topologytypes.hh>
28
#include <dune/geometry/axisalignedcubegeometry.hh>
25
29
#include <dune/grid/common/indexidset.hh>
26
30
#include <dune/grid/common/datahandleif.hh>
27
31
28
32
29
33
#if HAVE_MPI
30
#include <dune/common/mpicollectivecommunication.hh>
34
#include <dune/common/parallel/mpicollectivecommunication.hh>
31
35
#endif
32
36
33
37
/*! \file yaspgrid.hh
34
YaspGrid stands for yet another structured parallel grid.
35
It will implement the dune grid interface for structured grids with codim 0
36
and dim, with arbitrary overlap, parallel features with two overlap
37
models, periodic boundaries and fast a implementation allowing on-the-fly computations.
38
*/
39
40
namespace Dune {
41
42
//************************************************************************
43
/*! define name for floating point type used for coordinates in yaspgrid.
44
You can change the type for coordinates by changing this single typedef.
45
*/
46
typedef double yaspgrid_ctype;
47
static const yaspgrid_ctype yasptolerance=1E-13; // tolerance in coordinate computations
48
49
/* some sizes for building global ids
50
*/
51
const int yaspgrid_dim_bits = 24; // bits for encoding each dimension
52
const int yaspgrid_level_bits = 6; // bits for encoding level number
53
const int yaspgrid_codim_bits = 4; // bits for encoding codimension
54
55
56
//************************************************************************
57
// forward declaration of templates
58
59
template<int dim> class YaspGrid;
60
template<int mydim, int cdim, class GridImp> class YaspGeometry;
61
template<int codim, int dim, class GridImp> class YaspEntity;
62
template<int codim, class GridImp> class YaspEntityPointer;
63
template<int codim, class GridImp> class YaspEntitySeed;
64
template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
65
template<class GridImp> class YaspIntersectionIterator;
66
template<class GridImp> class YaspIntersection;
67
template<class GridImp> class YaspHierarchicIterator;
68
template<class GridImp> class YaspLevelIndexSet;
69
template<class GridImp> class YaspLeafIndexSet;
70
template<class GridImp> class YaspGlobalIdSet;
71
72
namespace FacadeOptions
73
{
74
75
template<int dim, int mydim, int cdim>
76
struct StoreGeometryReference<mydim, cdim, YaspGrid<dim>, YaspGeometry>
77
{
78
static const bool v = false;
79
};
80
81
template<int dim, int mydim, int cdim>
82
struct StoreGeometryReference<mydim, cdim, const YaspGrid<dim>, YaspGeometry>
83
{
84
static const bool v = false;
85
};
86
87
}
88
89
//========================================================================
90
// The transformation describing the refinement rule
91
92
template<int dim, class GridImp>
93
struct YaspFatherRelativeLocalElement {
94
static const array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > _geo;
95
static array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > initSons()
96
{
97
array<YaspGeometry<dim,dim,GridImp>, (1<<dim) > geo;
98
FieldVector<yaspgrid_ctype,dim> midpoint(0.25);
99
FieldVector<yaspgrid_ctype,dim> extension(0.5);
100
for (int i=0; i<(1<<dim); i++)
101
{
102
midpoint = 0.25;
103
for (int k=0; k<dim; k++)
104
{
105
if (i&(1<<k))
106
midpoint[k] = 0.75;
107
}
108
geo[i] = YaspGeometry<dim,dim,GridImp>(midpoint, extension);
109
}
110
return geo;
111
}
112
};
113
114
template<int dim, class GridImp>
115
const array<YaspGeometry<dim,dim,GridImp>, (1<<dim)>
116
YaspFatherRelativeLocalElement<dim, GridImp>::_geo =
117
YaspFatherRelativeLocalElement<dim, GridImp>::initSons();
118
119
//========================================================================
120
/*!
121
YaspGeometry realizes the concept of the geometric part of a mesh entity.
122
123
We have specializations for dim == dimworld (elements) and dim == 0
124
(vertices). The general version implements dim == dimworld-1 (faces)
125
and otherwise throws a GridError.
126
*/
127
//========================================================================
128
129
//! The general version implements dim==dimworld-1. If this is not the case an error is thrown
130
template<int mydim,int cdim, class GridImp>
131
class YaspGeometry : public GeometryDefaultImplementation<mydim,cdim,GridImp,YaspGeometry>
132
{
133
public:
134
//! define type used for coordinates in grid module
135
typedef typename GridImp::ctype ctype;
136
137
//! return the element type identifier
138
GeometryType type () const
139
{
140
return GeometryType(GeometryType::cube,mydim);
141
}
142
143
//! here we have always an affine geometry
144
bool affine() const { return true; }
145
146
//! return the number of corners of this element. Corners are numbered 0...n-1
147
int corners () const
148
{
149
return 1<<mydim;
150
}
151
152
//! access to coordinates of corners. Index is the number of the corner
153
FieldVector< ctype, cdim > corner ( const int i ) const
154
{
155
assert( i >= 0 && i < (int) coord_.N() );
156
FieldVector<ctype, cdim>& c = coord_[i];
157
int bit=0;
158
for (int k=0; k<cdim; k++) // run over all directions in world
159
{
160
if (k==missing)
161
{
162
c[k] = midpoint[k];
163
continue;
164
}
165
//k is not the missing direction
166
if (i&(1<<bit)) // check whether bit is set or not
167
c[k] = midpoint[k]+0.5*extension[k]; // bit is 1 in i
168
else
169
c[k] = midpoint[k]-0.5*extension[k]; // bit is 0 in i
170
bit++; // we have processed a direction
171
}
172
173
return c;
174
}
175
176
//! access to the center/centroid
177
FieldVector< ctype, cdim > center ( ) const
178
{
179
return midpoint;
180
}
181
182
//! maps a local coordinate within reference element to global coordinate in element
183
FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const
184
{
185
FieldVector<ctype, cdim> g;
186
int bit=0;
187
for (int k=0; k<cdim; k++)
188
if (k==missing)
189
g[k] = midpoint[k];
190
else
191
{
192
g[k] = midpoint[k] + (local[bit]-0.5)*extension[k];
193
bit++;
194
}
195
return g;
196
}
197
198
//! maps a global coordinate within the element to a local coordinate in its reference element
199
FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const
200
{
201
FieldVector<ctype, mydim> l; // result
202
int bit=0;
203
for (int k=0; k<cdim; k++)
204
if (k!=missing)
205
{
206
l[bit] = (global[k]-midpoint[k])/extension[k] + 0.5;
207
bit++;
208
}
209
return l;
210
}
211
212
//! return volume of geometry
213
ctype volume () const
214
{
215
ctype volume=1.0;
216
for (int k=0; k<cdim; k++)
217
if (k!=missing) volume *= extension[k];
218
return volume;
219
}
220
221
/*! determinant of the jacobian of the mapping
222
*/
223
ctype integrationElement (const FieldVector<ctype, mydim>& local) const
224
{
225
return volume();
226
}
227
228
//! Compute the transposed of the jacobi matrix
229
FieldMatrix<ctype,mydim,cdim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
230
{
231
JT = 0.0;
232
int k=0;
233
for (int i=0; i<cdim; ++i)
234
{
235
if (i != missing)
236
{
237
JT[k][i] = extension[i]; // set diagonal element
238
k++;
239
}
240
}
241
return JT;
242
}
243
//! Compute the transposed of the inverse jacobi matrix
244
FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
245
{
246
Jinv = 0.0;
247
int k=0;
248
for (int i=0; i<cdim; ++i)
249
{
250
if (i != missing)
251
{
252
Jinv[i][k] = 1.0/extension[i]; // set diagonal element
253
k++;
254
}
255
}
256
return Jinv;
257
}
258
259
//! default constructor
260
YaspGeometry () {}
261
262
//! constructor from midpoint and extension and missing direction number
263
YaspGeometry (const FieldVector<ctype, cdim>& p, const FieldVector<ctype, cdim>& h, uint8_t& m)
264
: midpoint(p), extension(h), missing(m)
265
{
266
if (cdim!=mydim+1)
267
DUNE_THROW(GridError, "general YaspGeometry assumes cdim=mydim+1");
268
}
269
270
//! copy constructor (skipping temporary variables)
271
YaspGeometry (const YaspGeometry& other)
272
: midpoint(other.midpoint),
273
extension(other.extension),
274
missing(other.missing)
275
{
276
}
277
278
//! print function
279
void print (std::ostream& s) const
280
{
281
s << "YaspGeometry<"<<mydim<<","<<cdim<< "> ";
282
s << "midpoint";
283
for (int i=0; i<cdim; i++)
284
s << " " << midpoint[i];
285
s << " extension";
286
for (int i=0; i<cdim; i++)
287
s << " " << extension[i];
288
s << " missing is " << missing;
289
}
290
291
// const YaspGeometry<mydim,cdim,GridImp>&
292
// operator = (const YaspGeometry<mydim,cdim,GridImp>& g);
293
294
private:
295
// the element is fully defined by its midpoint the extension
296
// in each direction and the missing direction.
297
// Note cdim == mydim+1
298
299
FieldVector<ctype, cdim> midpoint; // the midpoint
300
FieldVector<ctype, cdim> extension; // the extension
301
uint8_t missing; // the missing, i.e. constant direction
302
303
// In addition we need memory in order to return references.
304
// Possibly we should change this in the interface ...
305
mutable FieldMatrix<ctype, mydim, cdim> JT; // the transposed of the jacobian
306
mutable FieldMatrix<ctype, cdim, mydim> Jinv; // the transposed of the jacobian inverse
307
mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, cdim> coord_; // the coordinates
308
309
};
310
311
312
313
//! specialize for dim=dimworld, i.e. a volume element
314
template<int mydim, class GridImp>
315
class YaspGeometry<mydim,mydim,GridImp> : public GeometryDefaultImplementation<mydim,mydim,GridImp,YaspGeometry>
316
{
317
public:
318
typedef typename GridImp::ctype ctype;
319
320
//! return the element type identifier
321
GeometryType type () const
322
{
323
return GeometryType(GeometryType::cube,mydim);
324
}
325
326
//! here we have always an affine geometry
327
bool affine() const { return true; }
328
329
//! return the number of corners of this element. Corners are numbered 0...n-1
330
int corners () const
331
{
332
return 1<<mydim;
333
}
334
335
//! access to coordinates of corners. Index is the number of the corner
336
const FieldVector<ctype, mydim>& operator[] (int i) const
337
{
338
return corner(i);
339
}
340
341
//! access to coordinates of corners. Index is the number of the corner
342
FieldVector< ctype, mydim > corner ( const int i ) const
343
{
344
assert( i >= 0 && i < (int) coord_.N() );
345
FieldVector<ctype, mydim>& c = coord_[i];
346
for (int k=0; k<mydim; k++)
347
if (i&(1<<k))
348
c[k] = midpoint[k]+0.5*extension[k]; // kth bit is 1 in i
349
else
350
c[k] = midpoint[k]-0.5*extension[k]; // kth bit is 0 in i
351
return c;
352
}
353
354
//! access to the center/centroid
355
FieldVector< ctype, mydim > center ( ) const
356
{
357
return midpoint;
358
}
359
360
//! maps a local coordinate within reference element to global coordinate in element
361
FieldVector<ctype, mydim> global (const FieldVector<ctype, mydim>& local) const
362
{
363
FieldVector<ctype,mydim> g;
364
for (int k=0; k<mydim; k++)
365
g[k] = midpoint[k] + (local[k]-0.5)*extension[k];
366
return g;
367
}
368
369
//! maps a global coordinate within the element to a local coordinate in its reference element
370
FieldVector<ctype, mydim> local (const FieldVector<ctype,mydim>& global) const
371
{
372
FieldVector<ctype, mydim> l; // result
373
for (int k=0; k<mydim; k++)
374
l[k] = (global[k]-midpoint[k])/extension[k] + 0.5;
375
return l;
376
}
377
378
/*! determinant of the jacobian of the mapping
379
*/
380
ctype integrationElement (const FieldVector<ctype, mydim>& local) const
381
{
382
return volume();
383
}
384
385
//! return volume of geometry
386
ctype volume () const
387
{
388
ctype vol=1.0;
389
for (int k=0; k<mydim; k++) vol *= extension[k];
390
return vol;
391
}
392
393
//! Compute the transposed of the jacobi matrix
394
FieldMatrix<ctype,mydim,mydim>& jacobianTransposed (const FieldVector<ctype, mydim>& local) const
395
{
396
for (int i=0; i<mydim; ++i)
397
{
398
JT[i] = 0.0; // set column to zero
399
JT[i][i] = extension[i]; // set diagonal element
400
}
401
return JT;
402
}
403
//! Compute the transposed of the inverse jacobi matrix
404
FieldMatrix<ctype,mydim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
405
{
406
for (int i=0; i<mydim; ++i)
407
{
408
Jinv[i] = 0.0; // set column to zero
409
Jinv[i][i] = 1.0/extension[i]; // set diagonal element
410
}
411
return Jinv;
412
}
413
414
//! default constructor
415
YaspGeometry () {}
416
417
//! constructor from midpoint and extension
418
YaspGeometry (const FieldVector<ctype, mydim>& p, const FieldVector<ctype, mydim>& h)
419
: midpoint(p), extension(h)
420
{}
421
422
//! copy constructor (skipping temporary variables)
423
YaspGeometry (const YaspGeometry& other)
424
: midpoint(other.midpoint),
425
extension(other.extension)
426
{
427
}
428
429
//! print function
430
void print (std::ostream& s) const
431
{
432
s << "YaspGeometry<"<<mydim<<","<<mydim<< "> ";
433
s << "midpoint";
434
for (int i=0; i<mydim; i++)
435
s << " " << midpoint[i];
436
s << " extension";
437
for (int i=0; i<mydim; i++)
438
s << " " << extension[i];
439
}
440
441
// const YaspGeometry<mydim,mydim,GridImp>&
442
// operator = (const YaspGeometry<mydim,mydim,GridImp>& g);
443
444
private:
445
// the element is fully defined by midpoint and the extension
446
// in each direction. References are used because this information
447
// is known outside the element in many cases.
448
// Note mydim==cdim
449
450
FieldVector<ctype, mydim> midpoint; // the midpoint
451
FieldVector<ctype, mydim> extension; // the extension
452
453
// In addition we need memory in order to return references.
454
// Possibly we should change this in the interface ...
455
mutable FieldMatrix<ctype, mydim, mydim> Jinv,JT; // the transpose of the jacobian and its inverse inverse
456
mutable FieldMatrix<ctype, Power_m_p<2,mydim>::power, mydim> coord_; // the coordinates
457
};
458
459
//! specialization for dim=0, this is a vertex
460
template<int cdim, class GridImp>
461
class YaspGeometry<0,cdim,GridImp> : public GeometryDefaultImplementation<0,cdim,GridImp,YaspGeometry>
462
{
463
public:
464
typedef typename GridImp::ctype ctype;
465
466
//! return the element type identifier
467
GeometryType type () const
468
{
469
return GeometryType(GeometryType::cube,0);
470
}
471
472
//! here we have always an affine geometry
473
bool affine() const { return true; }
474
475
//! return the number of corners of this element. Corners are numbered 0...n-1
476
int corners () const
477
{
478
return 1;
479
}
480
481
//! access to coordinates of corners. Index is the number of the corner
482
const FieldVector<ctype, cdim>& operator[] (int i) const
483
{
484
return position;
485
}
486
487
//! access to coordinates of corners. Index is the number of the corner
488
FieldVector< ctype, cdim > corner ( const int i ) const
489
{
490
return position;
491
}
492
493
//! access to the center/centroid
494
FieldVector< ctype, cdim > center ( ) const
495
{
496
return position;
497
}
498
499
/*! determinant of the jacobian of the mapping
500
*/
501
ctype integrationElement (const FieldVector<ctype, 0>& local) const
502
{
503
return 1.0;
504
}
505
506
//! Compute the transposed of the jacobi matrix
507
FieldMatrix<ctype,0,cdim>& jacobianTransposed (const FieldVector<ctype, 0>& local) const
508
{
509
static FieldMatrix<ctype,0,cdim> JT(0.0);
510
return JT;
511
}
512
//! Compute the transposed of the inverse jacobi matrix
513
FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
514
{
515
static FieldMatrix<ctype,cdim,0> Jinv(0.0);
516
return Jinv;
517
}
518
519
//! default constructor
520
YaspGeometry ()
521
{}
522
523
//! constructor
524
explicit YaspGeometry ( const FieldVector< ctype, cdim > &p )
525
: position( p )
526
{}
527
528
YaspGeometry ( const FieldVector< ctype, cdim > &p, const FieldVector< ctype, cdim > &, uint8_t &)
529
: position( p )
530
{}
531
532
//! print function
533
void print (std::ostream& s) const
534
{
535
s << "YaspGeometry<"<<0<<","<<cdim<< "> ";
536
s << "position " << position;
537
}
538
539
// const YaspGeometry<0,cdim,GridImp>&
540
// operator = (const YaspGeometry<0,cdim,GridImp>& g);
541
542
private:
543
FieldVector<ctype, cdim> position; //!< position of the vertex
544
};
545
546
// operator<< for all YaspGeometrys
547
template <int mydim, int cdim, class GridImp>
548
inline
549
std::ostream& operator<< (std::ostream& s, YaspGeometry<mydim,cdim,GridImp>& e)
550
{
551
e.print(s);
552
return s;
553
}
554
555
//========================================================================
556
/*!
557
YaspEntity realizes the concept a mesh entity.
558
559
We have specializations for codim==0 (elements) and
560
codim=dim (vertices).
561
The general version throws a GridError.
562
*/
563
//========================================================================
564
565
template<int codim, int dim, class GridImp>
566
class YaspSpecialEntity :
567
public GridImp::template Codim<codim>::Entity
568
{
569
public:
570
typedef typename GridImp::ctype ctype;
571
572
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
573
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
574
575
YaspSpecialEntity(const GridImp* yg, const YGLI& g, const TSI& it) :
576
GridImp::template Codim<codim>::Entity (YaspEntity<codim, dim, GridImp>(yg,g,it))
577
{}
578
YaspSpecialEntity(const YaspEntity<codim, dim, GridImp>& e) :
579
GridImp::template Codim<codim>::Entity (e)
580
{}
581
const TSI& transformingsubiterator () const
582
{
583
return this->realEntity.transformingsubiterator();
584
}
585
const YGLI& gridlevel () const
586
{
587
return this->realEntity.gridlevel();
588
}
589
const GridImp * yaspgrid() const
590
{
591
return this->realEntity.yaspgrid();
592
}
593
};
594
595
template<int codim, int dim, class GridImp>
596
class YaspEntity
597
: public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
598
{
599
public:
600
typedef typename GridImp::ctype ctype;
601
602
typedef typename GridImp::template Codim<codim>::Geometry Geometry;
603
//! level of this element
604
int level () const
605
{
606
DUNE_THROW(GridError, "YaspEntity not implemented");
607
}
608
609
//! index is unique and consecutive per level and codim used for access to degrees of freedom
610
int index () const
611
{
612
DUNE_THROW(GridError, "YaspEntity not implemented");
613
}
614
615
//! geometry of this entity
616
Geometry geometry () const
617
{
618
DUNE_THROW(GridError, "YaspEntity not implemented");
619
}
620
621
//! return partition type attribute
622
PartitionType partitionType () const
623
{
624
DUNE_THROW(GridError, "YaspEntity not implemented");
625
}
626
627
const GridImp * yaspgrid() const
628
{
629
DUNE_THROW(GridError, "YaspEntity not implemented");
630
}
631
632
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
633
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
634
YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
635
{
636
DUNE_THROW(GridError, "YaspEntity not implemented");
637
}
638
639
// IndexSets needs access to the private index methods
640
friend class Dune::YaspLevelIndexSet<GridImp>;
641
friend class Dune::YaspLeafIndexSet<GridImp>;
642
friend class Dune::YaspGlobalIdSet<GridImp>;
643
typedef typename GridImp::PersistentIndexType PersistentIndexType;
644
645
//! globally unique, persistent index
646
PersistentIndexType persistentIndex () const
647
{
648
DUNE_THROW(GridError, "YaspEntity not implemented");
649
}
650
651
//! consecutive, codim-wise, level-wise index
652
int compressedIndex () const
653
{
654
DUNE_THROW(GridError, "YaspEntity not implemented");
655
}
656
657
//! consecutive, codim-wise, level-wise index
658
int compressedLeafIndex () const
659
{
660
DUNE_THROW(GridError, "YaspEntity not implemented");
661
}
662
};
663
664
665
// specialization for codim=0
666
template<int dim, class GridImp>
667
class YaspEntity<0,dim,GridImp>
668
: public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
669
{
670
enum { dimworld = GridImp::dimensionworld };
671
672
typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
673
674
public:
675
typedef typename GridImp::ctype ctype;
676
677
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
678
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
679
680
typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
681
typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
682
683
template <int cd>
684
struct Codim
685
{
686
typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
687
};
688
689
typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
690
typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
691
typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
692
typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
693
typedef typename GridImp::LeafIntersectionIterator LeafIntersectionIterator;
694
typedef typename GridImp::HierarchicIterator HierarchicIterator;
695
696
//! define the type used for persisitent indices
697
typedef typename GridImp::PersistentIndexType PersistentIndexType;
698
699
//! define type used for coordinates in grid module
700
typedef typename YGrid<dim,ctype>::iTupel iTupel;
701
702
// constructor
703
YaspEntity (const GridImp * yg, const YGLI& g, const TSI& it)
704
: _yg(yg), _it(it), _g(g)
705
{
706
}
707
708
//! level of this element
709
int level () const { return _g.level(); }
710
711
//! index is unique and consecutive per level
712
int index () const { return _it.superindex(); } // superindex works also for iteration over subgrids
713
714
//! globalIndex is unique and consecutive per global level
715
int globalIndex () const {
716
return _g.cell_global().index(_it.coord());
717
}
718
719
/** \brief Return the entity seed which contains sufficient information
720
* to generate the entity again and uses as less memory as possible
721
*/
722
EntitySeed seed () const {
723
return EntitySeed(_g.level(), _it.coord());
724
}
725
726
//! return partition type attribute
727
PartitionType partitionType () const
728
{
729
if (_g.cell_interior().inside(_it.coord())) return InteriorEntity;
730
if (_g.cell_overlap().inside(_it.coord())) return OverlapEntity;
731
DUNE_THROW(GridError, "Impossible GhostEntity " << _it.coord() << "\t"
732
<< _g.cell_interior().origin() << "/" << _g.cell_interior().size());
733
return GhostEntity;
734
}
735
736
//! geometry of this entity
737
Geometry geometry () const {
738
// the element geometry
739
GeometryImpl _geometry(_it.position(),_it.meshsize());
740
return Geometry( _geometry );
741
}
742
743
/*! Return number of subentities with codimension cc.
744
*/
745
template<int cc> int count () const
746
{
747
if (cc==dim) return 1<<dim;
748
if (cc==1) return 2*dim;
749
if (cc==dim-1) return dim*(1<<(dim-1));
750
if (cc==0) return 1;
751
DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
752
}
753
754
/*! Intra-element access to subentities of codimension cc > codim.
755
*/
756
template<int cc>
757
typename Codim<cc>::EntityPointer subEntity (int i) const
758
{
759
dune_static_assert( cc == dim || cc == 0 ,
760
"YaspGrid only supports Entities with codim=dim and codim=0");
761
// coordinates of the cell == coordinates of lower left corner
762
if (cc==dim)
763
{
764
iTupel coord = _it.coord();
765
766
// get corner from there
767
for (int k=0; k<dim; k++)
768
if (i&(1<<k)) (coord[k])++;
769
770
return YaspEntityPointer<cc,GridImp>(_yg,_g,_g.vertex_overlapfront().tsubbegin(coord));
771
}
772
if (cc==0)
773
{
774
return YaspEntityPointer<cc,GridImp>(_yg,_g,_it);
775
}
776
DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
777
}
778
779
//! Inter-level access to father element on coarser grid. Assumes that meshes are nested.
780
EntityPointer father () const
781
{
782
// check if coarse level exists
783
if (_g.level()<=0)
784
DUNE_THROW(GridError, "tried to call father on level 0");
785
786
// yes, get iterator to it
787
YGLI cg = _g.coarser();
788
789
// coordinates of the cell
790
iTupel coord = _it.coord();
791
792
// get coordinates on next coarser level
793
for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
794
795
return YaspEntityPointer<0,GridImp>(_yg,cg,cg.cell_overlap().tsubbegin(coord));
796
}
797
798
//! returns true if father entity exists
799
bool hasFather () const
800
{
801
return (_g.level()>0);
802
}
803
804
/*! Location of this element relative to the reference element element of the father.
805
This is sufficient to interpolate all dofs in conforming case.
806
Nonconforming case may require access to neighbors of father and
807
computations with local coordinates.
808
On the fly case is somewhat inefficient since dofs are visited several times.
809
If we store interpolation matrices, this is tolerable. We assume that on-the-fly
810
implementation of numerical algorithms is only done for simple discretizations.
811
Assumes that meshes are nested.
812
*/
813
LocalGeometry geometryInFather () const
814
{
815
// determine which son we are
816
int son = 0;
817
for (int k=0; k<dim; k++)
818
if (_it.coord(k)%2)
819
son += (1<<k);
820
821
// configure one of the 2^dim transformations
822
return LocalGeometry( YaspFatherRelativeLocalElement<dim,GridImp>::_geo[son] );
823
}
824
825
const TSI& transformingsubiterator () const
826
{
827
return _it;
828
}
829
830
const YGLI& gridlevel () const
831
{
832
return _g;
833
}
834
835
const GridImp* yaspgrid () const
836
{
837
return _yg;
838
}
839
840
bool isLeaf() const
841
{
842
return (_g.level() == _g.mg()->maxlevel());
843
}
844
845
/**\brief Returns true, if the entity has been created during the last call to adapt()
846
*/
847
bool isNew () const { return _yg->adaptRefCount > 0 && _g.mg()->maxlevel() < _g.level() + _yg->adaptRefCount; }
848
849
/**\brief Returns true, if entity might disappear during the next call to adapt()
850
*/
851
bool mightVanish () const { return false; }
852
// { return _yg->adaptRefCount < 0 && _g.mg()->maxlevel() < _g.level() - _yg->adaptRefCount; }
853
854
//! returns intersection iterator for first intersection
855
IntersectionIterator ibegin () const
856
{
857
return YaspIntersectionIterator<GridImp>(*this,false);
858
}
859
860
//! returns intersection iterator for first intersection
861
LeafIntersectionIterator ileafbegin () const
862
{
863
// only if entity is leaf this iterator delivers intersections
864
return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
865
}
866
867
//! returns intersection iterator for first intersection
868
LevelIntersectionIterator ilevelbegin () const
869
{
870
return ibegin();
871
}
872
873
//! Reference to one past the last neighbor
874
IntersectionIterator iend () const
875
{
876
return YaspIntersectionIterator<GridImp>(*this,true);
877
}
878
879
//! Reference to one past the last neighbor
880
LeafIntersectionIterator ileafend () const
881
{
882
return iend();
883
}
884
885
//! Reference to one past the last neighbor
886
LevelIntersectionIterator ilevelend () const
887
{
888
return iend();
889
}
890
891
/*! Inter-level access to son elements on higher levels<=maxlevel.
892
This is provided for sparsely stored nested unstructured meshes.
893
Returns iterator to first son.
894
*/
895
HierarchicIterator hbegin (int maxlevel) const
896
{
897
return YaspHierarchicIterator<GridImp>(_yg,_g,_it,maxlevel);
898
}
899
900
//! Returns iterator to one past the last son
901
HierarchicIterator hend (int maxlevel) const
902
{
903
return YaspHierarchicIterator<GridImp>(_yg,_g,_it,_g.level());
904
}
905
906
private:
907
// IndexSets needs access to the private index methods
908
friend class Dune::YaspLevelIndexSet<GridImp>;
909
friend class Dune::YaspLeafIndexSet<GridImp>;
910
friend class Dune::YaspGlobalIdSet<GridImp>;
911
912
//! globally unique, persistent index
913
PersistentIndexType persistentIndex () const
914
{
915
// get size of global grid
916
const iTupel& size = _g.cell_global().size();
917
918
// get coordinate correction for periodic boundaries
919
int coord[dim];
920
for (int i=0; i<dim; i++)
921
{
922
coord[i] = _it.coord(i);
923
if (coord[i]<0) coord[i] += size[i];
924
if (coord[i]>=size[i]) coord[i] -= size[i];
925
}
926
927
// encode codim
928
PersistentIndexType id(0);
929
930
// encode level
931
id = id << yaspgrid_level_bits;
932
id = id+PersistentIndexType(_g.level());
933
934
935
// encode coordinates
936
for (int i=dim-1; i>=0; i--)
937
{
938
id = id << yaspgrid_dim_bits;
939
id = id+PersistentIndexType(coord[i]);
940
}
941
942
return id;
943
}
944
945
//! consecutive, codim-wise, level-wise index
946
int compressedIndex () const
947
{
948
return _it.superindex();
949
}
950
951
//! consecutive, codim-wise, level-wise index
952
int compressedLeafIndex () const
953
{
954
return _it.superindex();
955
}
956
957
//! subentity persistent index
958
PersistentIndexType subPersistentIndex (int i, int cc) const
959
{
960
if (cc==0)
961
return persistentIndex();
962
963
// get position of cell, note that global origin is zero
964
// adjust for periodic boundaries
965
int coord[dim];
966
for (int k=0; k<dim; k++)
967
{
968
coord[k] = _it.coord(k);
969
if (coord[k]<0) coord[k] += _g.cell_global().size(k);
970
if (coord[k]>=_g.cell_global().size(k)) coord[k] -= _g.cell_global().size(k);
971
}
972
973
if (cc==dim)
974
{
975
// transform to vertex coordinates
976
for (int k=0; k<dim; k++)
977
if (i&(1<<k)) (coord[k])++;
978
979
// determine min number of trailing zeroes
980
int trailing = 1000;
981
for (int i=0; i<dim; i++)
982
{
983
// count trailing zeros
984
int zeros = 0;
985
for (int j=0; j<_g.level(); j++)
986
if (coord[i]&(1<<j))
987
break;
988
else
989
zeros++;
990
trailing = std::min(trailing,zeros);
991
}
992
993
// determine the level of this vertex
994
int level = _g.level()-trailing;
995
996
// encode codim
997
PersistentIndexType id(dim);
998
999
// encode level
1000
id = id << yaspgrid_level_bits;
1001
id = id+PersistentIndexType(level);
1002
1003
// encode coordinates
1004
for (int i=dim-1; i>=0; i--)
1005
{
1006
id = id << yaspgrid_dim_bits;
1007
id = id+PersistentIndexType(coord[i]>>trailing);
1008
}
1009
1010
return id;
1011
}
1012
1013
if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1014
{
1015
// Idea: Use the doubled grid to assign coordinates to faces
1016
1017
// ivar is the direction that varies
1018
int ivar=i/2;
1019
1020
// compute position from cell position
1021
for (int k=0; k<dim; k++)
1022
coord[k] = coord[k]*2 + 1; // the doubled grid
1023
if (i%2)
1024
coord[ivar] += 1;
1025
else
1026
coord[ivar] -= 1;
1027
1028
// encode codim
1029
PersistentIndexType id(1);
1030
1031
// encode level
1032
id = id << yaspgrid_level_bits;
1033
id = id+PersistentIndexType(_g.level());
1034
1035
// encode coordinates
1036
for (int i=dim-1; i>=0; i--)
1037
{
1038
id = id << yaspgrid_dim_bits;
1039
id = id+PersistentIndexType(coord[i]);
1040
}
1041
1042
return id;
1043
}
1044
1045
// map to old numbering
1046
static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1047
i = edge[i];
1048
1049
if (cc==dim-1) // edges, exist only for dim>2
1050
{
1051
// Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1052
1053
// number of entities per direction
1054
int m=1<<(dim-1);
1055
1056
// ifix is the direction that is fixed
1057
int ifix=(dim-1)-(i/m);
1058
1059
// compute position from cell position
1060
int bit=1;
1061
for (int k=0; k<dim; k++)
1062
{
1063
coord[k] = coord[k]*2+1; // cell position in doubled grid
1064
if (k==ifix) continue;
1065
if ((i%m)&bit) coord[k] += 1; else coord[k] -= 1;
1066
bit *= 2;
1067
}
1068
1069
// encode codim
1070
PersistentIndexType id(dim-1);
1071
1072
// encode level
1073
id = id << yaspgrid_level_bits;
1074
id = id+PersistentIndexType(_g.level());
1075
1076
// encode coordinates
1077
for (int i=dim-1; i>=0; i--)
1078
{
1079
id = id << yaspgrid_dim_bits;
1080
id = id+PersistentIndexType(coord[i]);
1081
}
1082
1083
return id;
1084
}
1085
1086
DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1087
}
1088
1089
//! subentity compressed index
1090
int subCompressedIndex (int i, int cc) const
1091
{
1092
if (cc==0)
1093
return compressedIndex();
1094
1095
// get cell position relative to origin of local cell grid
1096
iTupel coord;
1097
for (int k=0; k<dim; ++k)
1098
coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
1099
1100
if (cc==dim) // vertices
1101
{
1102
// transform cell coordinate to corner coordinate
1103
for (int k=0; k<dim; k++)
1104
if (i&(1<<k)) (coord[k])++;
1105
1106
// do lexicographic numbering
1107
int index = coord[dim-1];
1108
for (int k=dim-2; k>=0; --k)
1109
index = (index*(_g.cell_overlap().size(k)+1))+coord[k];
1110
return index;
1111
}
1112
1113
if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1114
{
1115
// Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
1116
1117
// ivar is the direction that varies
1118
int ivar=i/2;
1119
1120
// compute position from cell position
1121
if (i%2) coord[ivar] += 1;
1122
1123
// do lexicographic numbering
1124
int index = coord[dim-1];
1125
for (int k=dim-2; k>=0; --k)
1126
if (k==ivar)
1127
index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1128
else
1129
index = (index*(_g.cell_overlap().size(k)))+coord[k];
1130
1131
// add size of all subsets for smaller directions
1132
for (int j=0; j<ivar; j++)
1133
{
1134
int n=_g.cell_overlap().size(j)+1;
1135
for (int l=0; l<dim; l++)
1136
if (l!=j) n *= _g.cell_overlap().size(l);
1137
index += n;
1138
}
1139
1140
return index;
1141
}
1142
1143
// map to old numbering
1144
static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1145
i = edge[i];
1146
1147
if (cc==dim-1) // edges, exist only for dim>2
1148
{
1149
// Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1150
1151
// number of entities per direction
1152
int m=1<<(dim-1);
1153
1154
// ifix is the direction that is fixed
1155
int ifix=(dim-1)-(i/m);
1156
1157
// compute position from cell position
1158
int bit=1;
1159
for (int k=0; k<dim; k++)
1160
{
1161
if (k==ifix) continue;
1162
if ((i%m)&bit) coord[k] += 1;
1163
bit *= 2;
1164
}
1165
1166
// do lexicographic numbering
1167
int index = coord[dim-1];
1168
for (int k=dim-2; k>=0; --k)
1169
if (k!=ifix)
1170
index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1171
else
1172
index = (index*(_g.cell_overlap().size(k)))+coord[k];
1173
1174
// add size of all subsets for smaller directions
1175
for (int j=dim-1; j>ifix; j--)
1176
{
1177
int n=_g.cell_overlap().size(j);
1178
for (int l=0; l<dim; l++)
1179
if (l!=j) n *= _g.cell_overlap().size(l)+1;
1180
index += n;
1181
}
1182
1183
return index;
1184
}
1185
1186
DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1187
}
1188
1189
//! subentity compressed index
1190
int subCompressedLeafIndex (int i, int cc) const
1191
{
1192
if (cc==0)
1193
return compressedIndex();
1194
1195
// get cell position relative to origin of local cell grid
1196
iTupel coord;
1197
for (int k=0; k<dim; ++k)
1198
coord[k] = _it.coord(k)-_g.cell_overlap().origin(k);
1199
1200
if (cc==dim) // vertices
1201
{
1202
// transform cell coordinate to corner coordinate
1203
for (int k=0; k<dim; k++)
1204
if (i&(1<<k)) (coord[k])++;
1205
1206
// move coordinates up to maxlevel
1207
for (int k=0; k<dim; k++)
1208
coord[k] = coord[k]<<(_g.mg()->maxlevel()-_g.level());
1209
1210
// do lexicographic numbering
1211
int index = coord[dim-1];
1212
for (int k=dim-2; k>=0; --k)
1213
index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
1214
return index;
1215
}
1216
1217
if (cc==1) // faces, i.e. for dim=2 codim=1 is treated as a face
1218
{
1219
// Idea: direction ivar varies, all others are fixed, i.e. 2 possibilities per direction
1220
1221
// ivar is the direction that varies
1222
int ivar=i/2;
1223
1224
// compute position from cell position
1225
if (i%2) coord[ivar] += 1;
1226
1227
// do lexicographic numbering
1228
int index = coord[dim-1];
1229
for (int k=dim-2; k>=0; --k)
1230
if (k==ivar)
1231
index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1232
else
1233
index = (index*(_g.cell_overlap().size(k)))+coord[k];
1234
1235
// add size of all subsets for smaller directions
1236
for (int j=0; j<ivar; j++)
1237
{
1238
int n=_g.cell_overlap().size(j)+1;
1239
for (int l=0; l<dim; l++)
1240
if (l!=j) n *= _g.cell_overlap().size(l);
1241
index += n;
1242
}
1243
1244
return index;
1245
}
1246
1247
// map to old numbering
1248
static unsigned int edge[ 12 ] = { 0, 1, 2, 3, 4, 5, 8, 9, 6, 7, 10, 11 };
1249
i = edge[i];
1250
1251
if (cc==dim-1) // edges, exist only for dim>2
1252
{
1253
// Idea: direction i is fixed, all others are vary, i.e. 2^(dim-1) possibilities per direction
1254
1255
// number of entities per direction
1256
int m=1<<(dim-1);
1257
1258
// ifix is the direction that is fixed
1259
int ifix=(dim-1)-(i/m);
1260
1261
// compute position from cell position
1262
int bit=1;
1263
for (int k=0; k<dim; k++)
1264
{
1265
if (k==ifix) continue;
1266
if ((i%m)&bit) coord[k] += 1;
1267
bit *= 2;
1268
}
1269
1270
// do lexicographic numbering
1271
int index = coord[dim-1];
1272
for (int k=dim-2; k>=0; --k)
1273
if (k!=ifix)
1274
index = (index*(_g.cell_overlap().size(k)+1))+coord[k]; // one more
1275
else
1276
index = (index*(_g.cell_overlap().size(k)))+coord[k];
1277
1278
// add size of all subsets for smaller directions
1279
for (int j=dim-1; j>ifix; j--)
1280
{
1281
int n=_g.cell_overlap().size(j);
1282
for (int l=0; l<dim; l++)
1283
if (l!=j) n *= _g.cell_overlap().size(l)+1;
1284
index += n;
1285
}
1286
1287
return index;
1288
}
1289
1290
DUNE_THROW(GridError, "codim " << cc << " (dim=" << dim << ") not (yet) implemented");
1291
}
1292
1293
const GridImp * _yg; // access to YaspGrid
1294
const TSI& _it; // position in the grid level
1295
const YGLI& _g; // access to grid level
1296
};
1297
1298
1299
// specialization for codim=dim (vertex)
1300
template<int dim, class GridImp>
1301
class YaspEntity<dim,dim,GridImp>
1302
: public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
1303
{
1304
enum { dimworld = GridImp::dimensionworld };
1305
1306
typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
1307
1308
public:
1309
typedef typename GridImp::ctype ctype;
1310
1311
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1312
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1313
1314
typedef typename GridImp::template Codim<dim>::Geometry Geometry;
1315
1316
template <int cd>
1317
struct Codim
1318
{
1319
typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
1320
};
1321
1322
typedef typename GridImp::template Codim<dim>::EntityPointer EntityPointer;
1323
typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
1324
1325
//! define the type used for persisitent indices
1326
typedef typename GridImp::PersistentIndexType PersistentIndexType;
1327
1328
//! define type used for coordinates in grid module
1329
typedef typename YGrid<dim,ctype>::iTupel iTupel;
1330
1331
// constructor
1332
YaspEntity (const GridImp* yg, const YGLI& g, const TSI& it)
1333
: _yg(yg), _it(it), _g(g)
1334
{}
1335
1336
//! level of this element
1337
int level () const {return _g.level();}
1338
1339
//! index is unique and consecutive per level
1340
int index () const {return _it.superindex();}
1341
1342
//! globally unique, persistent index
1343
int globalIndex () const { return _g.cell_global().index(_it.coord()); }
1344
1345
/** \brief Return the entity seed which contains sufficient information
1346
* to generate the entity again and uses as less memory as possible
1347
*/
1348
EntitySeed seed () const {
1349
return EntitySeed(_g.level(), _it.coord());
1350
}
1351
1352
//! geometry of this entity
1353
Geometry geometry () const {
1354
GeometryImpl _geometry(_it.position());
1355
return Geometry( _geometry );
1356
}
1357
1358
//! return partition type attribute
1359
PartitionType partitionType () const
1360
{
1361
if (_g.vertex_interior().inside(_it.coord())) return InteriorEntity;
1362
if (_g.vertex_interiorborder().inside(_it.coord())) return BorderEntity;
1363
if (_g.vertex_overlap().inside(_it.coord())) return OverlapEntity;
1364
if (_g.vertex_overlapfront().inside(_it.coord())) return FrontEntity;
1365
return GhostEntity;
1366
}
1367
1368
private:
1369
// IndexSets needs access to the private index methods
1370
friend class Dune::YaspLevelIndexSet<GridImp>;
1371
friend class Dune::YaspLeafIndexSet<GridImp>;
1372
friend class Dune::YaspGlobalIdSet<GridImp>;
1373
1374
//! globally unique, persistent index
1375
PersistentIndexType persistentIndex () const
1376
{
1377
// get coordinate and size of global grid
1378
const iTupel& size = _g.vertex_global().size();
1379
int coord[dim];
1380
1381
// correction for periodic boundaries
1382
for (int i=0; i<dim; i++)
1383
{
1384
coord[i] = _it.coord(i);
1385
if (coord[i]<0) coord[i] += size[i];
1386
if (coord[i]>=size[i]) coord[i] -= size[i];
1387
}
1388
1389
// determine min number of trailing zeroes
1390
int trailing = 1000;
1391
for (int i=0; i<dim; i++)
1392
{
1393
// count trailing zeros
1394
int zeros = 0;
1395
for (int j=0; j<_g.level(); j++)
1396
if (coord[i]&(1<<j))
1397
break;
1398
else
1399
zeros++;
1400
trailing = std::min(trailing,zeros);
1401
}
1402
1403
// determine the level of this vertex
1404
int level = _g.level()-trailing;
1405
1406
// encode codim
1407
PersistentIndexType id(dim);
1408
1409
// encode level
1410
id = id << yaspgrid_level_bits;
1411
id = id+PersistentIndexType(level);
1412
1413
// encode coordinates
1414
for (int i=dim-1; i>=0; i--)
1415
{
1416
id = id << yaspgrid_dim_bits;
1417
id = id+PersistentIndexType(coord[i]>>trailing);
1418
}
1419
1420
return id;
1421
}
1422
1423
//! consecutive, codim-wise, level-wise index
1424
int compressedIndex () const { return _it.superindex();}
1425
1426
//! consecutive, codim-wise, level-wise index
1427
int compressedLeafIndex () const
1428
{
1429
if (_g.level()==_g.mg()->maxlevel())
1430
return _it.superindex();
1431
1432
// not on leaf level, interpolate to finest grid
1433
int coord[dim];
1434
for (int i=0; i<dim; i++) coord[i] = _it.coord(i)-(_g).vertex_overlap().origin(i);
1435
1436
// move coordinates up to maxlevel (multiply by 2 for each level
1437
for (int k=0; k<dim; k++)
1438
coord[k] = coord[k]*(1<<(_g.mg()->maxlevel()-_g.level()));
1439
1440
// do lexicographic numbering
1441
int index = coord[dim-1];
1442
for (int k=dim-2; k>=0; --k)
1443
index = (index*(_g.mg()->rbegin().cell_overlap().size(k)+1))+coord[k];
1444
return index;
1445
}
1446
1447
public:
1448
const TSI& transformingsubiterator() const { return _it; }
1449
const YGLI& gridlevel() const { return _g; }
1450
const GridImp * yaspgrid() const { return _yg; }
1451
protected:
1452
const GridImp * _yg; // access to YaspGrid
1453
const TSI& _it; // position in the grid level
1454
const YGLI& _g; // access to grid level
1455
// temporary object
1456
mutable FieldVector<ctype, dim> loc; // always computed before being returned
1457
};
1458
1459
//========================================================================
1460
/*!
1461
YaspIntersection provides data about intersection with
1462
neighboring codim 0 entities.
1463
*/
1464
//========================================================================
1465
1466
template<class GridImp>
1467
class YaspIntersection
1468
{
1469
enum { dim=GridImp::dimension };
1470
enum { dimworld=GridImp::dimensionworld };
1471
typedef typename GridImp::ctype ctype;
1472
YaspIntersection();
1473
YaspIntersection& operator = (const YaspIntersection&);
1474
1475
typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
1476
typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
1477
1478
public:
1479
// types used from grids
1480
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1481
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1482
typedef typename GridImp::template Codim<0>::Entity Entity;
1483
typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
1484
typedef typename GridImp::template Codim<1>::Geometry Geometry;
1485
typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
1486
typedef YaspSpecialEntity<0,dim,GridImp> SpecialEntity;
1487
typedef Dune::Intersection<const GridImp, Dune::YaspIntersectionIterator> Intersection;
1488
1489
// void update() const {
1490
// const_cast<YaspIntersection*>(this)->update();
1491
// }
1492
void update() const {
1493
if (_count == 2*_dir + _face || _count >= 2*dim)
1494
return;
1495
1496
// cleanup old stuff
1497
_outside.transformingsubiterator().move(_dir,1-2*_face); // move home
1498
_pos_world[_dir] = _inside.transformingsubiterator().position(_dir);
1499
1500
// update face info
1501
_dir = _count / 2;
1502
_face = _count % 2;
1503
1504
// move transforming iterator
1505
_outside.transformingsubiterator().move(_dir,-1+2*_face);
1506
1507
// make up faces
1508
_pos_world[_dir] += (-0.5+_face)*_inside.transformingsubiterator().meshsize(_dir);
1509
}
1510
1511
//! increment
1512
void increment() {
1513
_count += (_count < 2*dim);
1514
}
1515
1516
//! equality
1517
bool equals (const YaspIntersection& other) const
1518
{
1519
return _inside.equals(other._inside) && _count == other._count;
1520
}
1521
1522
/*! return true if neighbor ist outside the domain. Still the neighbor might
1523
exist in case of periodic boundary conditions, i.e. true is returned
1524
if the neighbor is outside the periodic unit cell
1525
*/
1526
bool boundary () const
1527
{
1528
#if 1
1529
return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 < _inside.gridlevel().cell_global().min(_count/2)
1530
||
1531
_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 > _inside.gridlevel().cell_global().max(_count/2));
1532
#else
1533
update();
1534
// The transforming iterator can be safely moved beyond the boundary.
1535
// So we only have to compare against the cell_global grid
1536
return (_outside.transformingsubiterator().coord(_dir) < _inside.gridlevel().cell_global().min(_dir)
1537
||
1538
_outside.transformingsubiterator().coord(_dir) > _inside.gridlevel().cell_global().max(_dir));
1539
#endif
1540
}
1541
1542
//! return true if neighbor across intersection exists in this processor
1543
bool neighbor () const
1544
{
1545
#if 1
1546
return (_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 >= _inside.gridlevel().cell_overlap().min(_count/2)
1547
&&
1548
_inside.transformingsubiterator().coord(_count/2) + 2*(_count%2) - 1 <= _inside.gridlevel().cell_overlap().max(_count/2));
1549
#else
1550
update();
1551
return (_outside.transformingsubiterator().coord(_dir) >= _inside.gridlevel().cell_overlap().min(_dir)
1552
&&
1553
_outside.transformingsubiterator().coord(_dir) <= _inside.gridlevel().cell_overlap().max(_dir));
1554
#endif
1555
}
1556
1557
//! Yasp is always conform
1558
bool conforming () const
1559
{
1560
return true;
1561
}
1562
1563
//! return EntityPointer to the Entity on the inside of this intersection
1564
//! (that is the Entity where we started this Iterator)
1565
EntityPointer inside() const
1566
{
1567
return _inside;
1568
}
1569
1570
//! return EntityPointer to the Entity on the outside of this intersection
1571
//! (that is the neighboring Entity)
1572
EntityPointer outside() const
1573
{
1574
update();
1575
return _outside;
1576
}
1577
1578
//! identifier for boundary segment from macro grid
1579
//! (attach your boundary condition as needed)
1580
int boundaryId() const
1581
{
1582
if(boundary()) return indexInInside()+1;
1583
return 0;
1584
}
1585
1586
//! identifier for boundary segment from macro grid
1587
//! (attach your boundary condition as needed)
1588
int boundarySegmentIndex() const
1589
{
1590
if(! boundary())
1591
DUNE_THROW(GridError, "called boundarySegmentIndex while boundary() == false");
1592
update();
1593
// size of local macro grid
1594
const FieldVector<int, dim> & size = _inside.gridlevel().mg()->begin().cell_overlap().size();
1595
const FieldVector<int, dim> & origin = _inside.gridlevel().mg()->begin().cell_overlap().origin();
1596
FieldVector<int, dim> sides;
1597
{
1598
for (int i=0; i<dim; i++)
1599
{
1600
sides[i] =
1601
((_inside.gridlevel().mg()->begin().cell_overlap().origin(i)
1602
== _inside.gridlevel().mg()->begin().cell_global().origin(i))+
1603
(_inside.gridlevel().mg()->begin().cell_overlap().origin(i) +
1604
_inside.gridlevel().mg()->begin().cell_overlap().size(i)
1605
== _inside.gridlevel().mg()->begin().cell_global().origin(i) +
1606
_inside.gridlevel().mg()->begin().cell_global().size(i)));
1607
}
1608
}
1609
// gobal position of the cell on macro grid
1610
FieldVector<int, dim> pos = _inside.transformingsubiterator().coord();
1611
pos /= (1<<_inside.level());
1612
pos -= origin;
1613
// compute unit-cube-face-sizes
1614
FieldVector<int, dim> fsize;
1615
{
1616
int vol = 1;
1617
for (int k=0; k<dim; k++)
1618
vol *= size[k];
1619
for (int k=0; k<dim; k++)
1620
fsize[k] = vol/size[k];
1621
}
1622
// compute index in the unit-cube-face
1623
int index = 0;
1624
{
1625
int localoffset = 1;
1626
for (int k=dim-1; k>=0; k--)
1627
{
1628
if (k == _dir) continue;
1629
index += (pos[k]) * localoffset;
1630
localoffset *= size[k];
1631
}
1632
}
1633
// add unit-cube-face-offsets
1634
{
1635
for (int k=0; k<_dir; k++)
1636
index += sides[k] * fsize[k];
1637
// add fsize if we are on the right face and there is a left-face-boundary on this processor
1638
index += _face * (sides[_dir]>1) * fsize[_dir];
1639
}
1640
1641
// int rank = 0;
1642
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1643
// std::cout << rank << "... size: " << size << " sides: " << sides
1644
// << " fsize: " << fsize
1645
// << " pos: " << pos << " face: " << int(_dir) << "/" << int(_face)
1646
// << " index: " << index << std::endl;
1647
1648
return index;
1649
}
1650
1651
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
1652
FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dim-1>& local) const
1653
{
1654
return _faceInfo[_count].normal;
1655
}
1656
1657
//! return unit outer normal, this should be dependent on local coordinates for higher order boundary
1658
FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dim-1>& local) const
1659
{
1660
return _faceInfo[_count].normal;
1661
}
1662
1663
//! return unit outer normal at center of intersection geometry
1664
FieldVector<ctype, dimworld> centerUnitOuterNormal () const
1665
{
1666
return _faceInfo[_count].normal;
1667
}
1668
1669
//! return unit outer normal, this should be dependent on
1670
//! local coordinates for higher order boundary
1671
//! the normal is scaled with the integration element of the intersection.
1672
FieldVector<ctype, dimworld> integrationOuterNormal (const FieldVector<ctype, dim-1>& local) const
1673
{
1674
FieldVector<ctype, dimworld> n = _faceInfo[_count].normal;
1675
n *= geometry().volume();
1676
return n;
1677
}
1678
1679
/*! intersection of codimension 1 of this neighbor with element where iteration started.
1680
Here returned element is in LOCAL coordinates of the element where iteration started.
1681
*/
1682
LocalGeometry geometryInInside () const
1683
{
1684
return LocalGeometry( _faceInfo[_count].geom_inside );
1685
}
1686
1687
/*! intersection of codimension 1 of this neighbor with element where iteration started.
1688
Here returned element is in LOCAL coordinates of neighbor
1689
*/
1690
LocalGeometry geometryInOutside () const
1691
{
1692
return LocalGeometry( _faceInfo[_count].geom_outside );
1693
}
1694
1695
/*! intersection of codimension 1 of this neighbor with element where iteration started.
1696
*/
1697
Geometry geometry () const
1698
{
1699
update();
1700
GeometryImpl
1701
_is_global(_pos_world,_inside.transformingsubiterator().meshsize(),_dir);
1702
return Geometry( _is_global );
1703
}
1704
1705
/** \brief obtain the type of reference element for this intersection */
1706
GeometryType type () const
1707
{
1708
static const GeometryType cube(GeometryType::cube, dim-1);
1709
return cube;
1710
}
1711
1712
//! local index of codim 1 entity in self where intersection is contained in
1713
int indexInInside () const
1714
{
1715
return _count;
1716
}
1717
1718
//! local index of codim 1 entity in neighbor where intersection is contained in
1719
int indexInOutside () const
1720
{
1721
// flip the last bit
1722
return _count^1;
1723
}
1724
1725
//! make intersection iterator from entity, initialize to first neighbor
1726
YaspIntersection (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
1727
_inside(myself.yaspgrid(), myself.gridlevel(),
1728
myself.transformingsubiterator()),
1729
_outside(myself.yaspgrid(), myself.gridlevel(),
1730
myself.transformingsubiterator()),
1731
// initialize to first neighbor
1732
_count(0),
1733
_dir(0),
1734
_face(0),
1735
_pos_world(myself.transformingsubiterator().position())
1736
{
1737
if (toend)
1738
{
1739
// initialize end iterator
1740
_count = 2*dim;
1741
return;
1742
}
1743
_count = 0;
1744
1745
// move transforming iterator
1746
_outside.transformingsubiterator().move(_dir,-1);
1747
1748
// make up faces
1749
_pos_world[0] -= 0.5*_inside.transformingsubiterator().meshsize(0);
1750
}
1751
1752
//! copy constructor
1753
YaspIntersection (const YaspIntersection& it) :
1754
_inside(it._inside),
1755
_outside(it._outside),
1756
_count(it._count),
1757
_dir(it._dir),
1758
_face(it._face),
1759
_pos_world(it._pos_world)
1760
{}
1761
1762
//! copy operator
1763
void assign (const YaspIntersection& it)
1764
{
1765
_inside = it._inside;
1766
_outside = it._outside;
1767
_count = it._count;
1768
_dir = it._dir;
1769
_face = it._face;
1770
_pos_world = it._pos_world;
1771
}
1772
1773
private:
1774
/* EntityPointers (get automatically updated) */
1775
mutable YaspEntityPointer<0,GridImp> _inside; //!< entitypointer to myself
1776
mutable YaspEntityPointer<0,GridImp> _outside; //!< outside entitypointer
1777
/* current position */
1778
uint8_t _count; //!< valid neighbor count in 0 .. 2*dim-1
1779
mutable uint8_t _dir; //!< count/2
1780
mutable uint8_t _face; //!< count%2
1781
/* current position */
1782
mutable FieldVector<ctype, dimworld> _pos_world; //!< center of face in world coordinates
1783
1784
/* static data */
1785
struct faceInfo
1786
{
1787
FieldVector<ctype, dimworld> normal;
1788
LocalGeometryImpl geom_inside; //!< intersection in own local coordinates
1789
LocalGeometryImpl geom_outside; //!< intersection in neighbors local coordinates
1790
};
1791
1792
/* static face info */
1793
static const array<faceInfo, 2*GridImp::dimension> _faceInfo;
1794
1795
static array<faceInfo, 2*dim> initFaceInfo()
1796
{
1797
const FieldVector<typename GridImp::ctype, GridImp::dimension> ext_local(1.0);
1798
array<faceInfo, 2*dim> I;
1799
for (uint8_t i=0; i<dim; i++)
1800
{
1801
// center of face
1802
FieldVector<ctype, dim> a(0.5); a[i] = 0.0;
1803
FieldVector<ctype, dim> b(0.5); b[i] = 1.0;
1804
// normal vectors
1805
I[2*i].normal = 0.0;
1806
I[2*i+1].normal = 0.0;
1807
I[2*i].normal[i] = -1.0;
1808
I[2*i+1].normal[i] = +1.0;
1809
// geometries
1810
I[2*i].geom_inside =
1811
LocalGeometryImpl(a, ext_local, i);
1812
I[2*i].geom_outside =
1813
LocalGeometryImpl(b, ext_local, i);
1814
I[2*i+1].geom_inside =
1815
LocalGeometryImpl(b, ext_local, i);
1816
I[2*i+1].geom_outside =
1817
LocalGeometryImpl(a, ext_local, i);
1818
}
1819
return I;
1820
}
1821
};
1822
1823
template<class GridImp>
1824
const array<typename YaspIntersection<GridImp>::faceInfo, 2*GridImp::dimension>
1825
YaspIntersection<GridImp>::_faceInfo =
1826
YaspIntersection<GridImp>::initFaceInfo();
1827
1828
//========================================================================
1829
/*!
1830
YaspIntersectionIterator enables iteration over intersection with
1831
neighboring codim 0 entities.
1832
*/
1833
//========================================================================
1834
1835
template<class GridImp>
1836
class YaspIntersectionIterator : public MakeableInterfaceObject< Dune::Intersection<const GridImp, Dune::YaspIntersection > >
1837
{
1838
enum { dim=GridImp::dimension };
1839
enum { dimworld=GridImp::dimensionworld };
1840
typedef typename GridImp::ctype ctype;
1841
YaspIntersectionIterator();
1842
public:
1843
// types used from grids
1844
typedef Dune::Intersection<const GridImp, Dune::YaspIntersection> Intersection;
1845
typedef MakeableInterfaceObject<Intersection> MakeableIntersection;
1846
1847
//! increment
1848
void increment()
1849
{
1850
GridImp::getRealImplementation(*this).increment();
1851
}
1852
1853
//! equality
1854
bool equals (const YaspIntersectionIterator& other) const
1855
{
1856
return GridImp::getRealImplementation(*this).equals(
1857
GridImp::getRealImplementation(other));
1858
}
1859
1860
//! \brief dereferencing
1861
const Intersection & dereference() const
1862
{
1863
return *this;
1864
}
1865
1866
//! make intersection iterator from entity
1867
YaspIntersectionIterator (const YaspEntity<0,dim,GridImp>& myself, bool toend) :
1868
MakeableIntersection(YaspIntersection<GridImp>(myself, toend))
1869
{}
1870
1871
//! copy constructor
1872
YaspIntersectionIterator (const YaspIntersectionIterator& it) :
1873
MakeableIntersection(it)
1874
{}
1875
1876
//! assignment
1877
YaspIntersectionIterator & operator = (const YaspIntersectionIterator& it)
1878
{
1879
GridImp::getRealImplementation(*this).assign(
1880
GridImp::getRealImplementation(it));
1881
return *this;
1882
}
1883
};
1884
1885
//========================================================================
1886
/*!
1887
YaspHierarchicIterator enables iteration over son entities of codim 0
1888
*/
1889
//========================================================================
1890
1891
template<class GridImp>
1892
class YaspHierarchicIterator :
1893
public YaspEntityPointer<0,GridImp>
1894
{
1895
enum { dim=GridImp::dimension };
1896
enum { dimworld=GridImp::dimensionworld };
1897
typedef typename GridImp::ctype ctype;
1898
public:
1899
// types used from grids
1900
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
1901
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
1902
typedef typename GridImp::template Codim<0>::Entity Entity;
1903
typedef YaspSpecialEntity<0,dim,GridImp> SpecialEntity;
1904
1905
//! define type used for coordinates in grid module
1906
typedef typename YGrid<dim,ctype>::iTupel iTupel;
1907
1908
//! constructor
1909
YaspHierarchicIterator (const GridImp* yg, const YGLI& g, const TSI& it, int maxlevel) :
1910
YaspEntityPointer<0,GridImp>(yg,g,it)
1911
{
1912
// now iterator points to current cell
1913
StackElem se(this->_g);
1914
se.coord = this->_it.coord();
1915
stack.push(se);
1916
1917
// determine maximum level
1918
_maxlevel = std::min(maxlevel,this->_g.mg()->maxlevel());
1919
1920
// if maxlevel not reached then push yourself and sons
1921
if (this->_g.level()<_maxlevel)
1922
{
1923
push_sons();
1924
}
1925
1926
// and make iterator point to first son if stack is not empty
1927
if (!stack.empty())
1928
pop_tos();
1929
}
1930
1931
//! constructor
1932
YaspHierarchicIterator (const YaspHierarchicIterator& it) :
1933
YaspEntityPointer<0,GridImp>(it),
1934
_maxlevel(it._maxlevel), stack(it.stack)
1935
{}
1936
1937
//! increment
1938
void increment ()
1939
{
1940
// sanity check: do nothing when stack is empty
1941
if (stack.empty()) return;
1942
1943
// if maxlevel not reached then push sons
1944
if (this->_g.level()<_maxlevel)
1945
push_sons();
1946
1947
// in any case pop one element
1948
pop_tos();
1949
}
1950
1951
void print (std::ostream& s) const
1952
{
1953
s << "HIER: " << "level=" << this->_g.level()
1954
<< " position=" << this->_it.coord()
1955
<< " superindex=" << this->_it.superindex()
1956
<< " maxlevel=" << this->_maxlevel
1957
<< " stacksize=" << stack.size()
1958
<< std::endl;
1959
}
1960
1961
private:
1962
int _maxlevel; //!< maximum level of elements to be processed
1963
1964
struct StackElem {
1965
YGLI g; // grid level of the element
1966
iTupel coord; // and the coordinates
1967
StackElem(YGLI gg) : g(gg) {}
1968
};
1969
std::stack<StackElem> stack; //!< stack holding elements to be processed
1970
1971
// push sons of current element on the stack
1972
void push_sons ()
1973
{
1974
// yes, process all 1<<dim sons
1975
StackElem se(this->_g.finer());
1976
for (int i=0; i<(1<<dim); i++)
1977
{
1978
for (int k=0; k<dim; k++)
1979
if (i&(1<<k))
1980
se.coord[k] = this->_it.coord(k)*2+1;
1981
else
1982
se.coord[k] = this->_it.coord(k)*2;
1983
stack.push(se);
1984
}
1985
}
1986
1987
// make TOS the current element
1988
void pop_tos ()
1989
{
1990
StackElem se = stack.top();
1991
stack.pop();
1992
this->_g = se.g;
1993
this->_it.reinit(this->_g.cell_overlap(),se.coord);
1994
}
1995
};
1996
1997
//========================================================================
1998
/*!
1999
YaspEntitySeed describes the minimal information necessary to create a fully functional YaspEntity
2000
*/
2001
//========================================================================
2002
template<int codim, class GridImp>
2003
class YaspEntitySeed
2004
{
2005
//! know your own dimension
2006
enum { dim=GridImp::dimension };
2007
2008
public:
2009
//! codimension of entity pointer
2010
enum { codimension = codim };
2011
2012
//! constructor
2013
YaspEntitySeed (int level, FieldVector<int, dim> coord)
2014
: _l(level), _c(coord)
2015
{}
2016
2017
//! copy constructor
2018
YaspEntitySeed (const YaspEntitySeed& rhs)
2019
: _l(rhs._l), _c(rhs._c)
2020
{}
2021
2022
int level () const { return _l; }
2023
const FieldVector<int, dim> & coord() const { return _c; }
2024
2025
protected:
2026
int _l; // grid level
2027
FieldVector<int, dim> _c; // coord in the global grid
2028
};
2029
2030
//========================================================================
2031
/*!
2032
YaspEntityPointer serves as a Reference or Pointer to a YaspGrid::Entity.
2033
It can also be initialized from Yasp::LevelIterator, Yasp::LeafIterator,
2034
Yasp::HierarchicIterator and Yasp::IntersectionIterator.
2035
2036
We have specializations for codim==0 (elements) and
2037
codim=dim (vertices).
2038
The general version throws a GridError.
2039
*/
2040
//========================================================================
2041
template<int codim, class GridImp>
2042
class YaspEntityPointer
2043
{
2044
//! know your own dimension
2045
enum { dim=GridImp::dimension };
2046
//! know your own dimension of world
2047
enum { dimworld=GridImp::dimensionworld };
2048
typedef typename GridImp::ctype ctype;
2049
public:
2050
typedef typename GridImp::template Codim<codim>::Entity Entity;
2051
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2052
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2053
typedef YaspSpecialEntity<codim,dim,GridImp> SpecialEntity;
2054
typedef YaspEntityPointer<codim,GridImp> EntityPointerImp;
2055
protected:
2056
typedef YaspEntity<codim, dim, GridImp> YaspEntityImp;
2057
2058
public:
2059
//! codimension of entity pointer
2060
enum { codimension = codim };
2061
2062
//! constructor
2063
YaspEntityPointer (const GridImp * yg, const YGLI & g, const TSI & it)
2064
: _g(g), _it(it), _entity(yg, _g,_it)
2065
{
2066
if (codim>0 && codim<dim)
2067
{
2068
DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2069
}
2070
}
2071
2072
//! copy constructor
2073
YaspEntityPointer (const YaspEntityImp& entity)
2074
: _g(entity.gridlevel()),
2075
_it(entity.transformingsubiterator()),
2076
_entity(entity.yaspgrid(),_g,_it)
2077
{
2078
if (codim>0 && codim<dim)
2079
{
2080
DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2081
}
2082
}
2083
2084
//! copy constructor
2085
YaspEntityPointer (const YaspEntityPointer& rhs)
2086
: _g(rhs._g), _it(rhs._it), _entity(rhs._entity.yaspgrid(),_g,_it)
2087
{
2088
if (codim>0 && codim<dim)
2089
{
2090
DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2091
}
2092
}
2093
2094
//! equality
2095
bool equals (const YaspEntityPointer& rhs) const
2096
{
2097
return (_it==rhs._it && _g == rhs._g);
2098
}
2099
2100
//! dereferencing
2101
Entity& dereference() const
2102
{
2103
return _entity;
2104
}
2105
2106
//! ask for level of entity
2107
int level () const {return _g.level();}
2108
2109
const YaspEntityPointer&
2110
operator = (const YaspEntityPointer& rhs)
2111
{
2112
_g = rhs._g;
2113
_it = rhs._it;
2114
/* _entity = i._entity
2115
* is done implicitely, as the entity is completely
2116
* defined via the interator it belongs to
2117
*/
2118
return *this;
2119
}
2120
2121
const TSI& transformingsubiterator () const
2122
{
2123
return _it;
2124
}
2125
2126
const YGLI& gridlevel () const
2127
{
2128
return _g;
2129
}
2130
2131
TSI& transformingsubiterator ()
2132
{
2133
return _it;
2134
}
2135
2136
YGLI& gridlevel ()
2137
{
2138
return _g;
2139
}
2140
2141
protected:
2142
YGLI _g; // access to grid level
2143
TSI _it; // position in the grid level
2144
mutable SpecialEntity _entity; //!< virtual entity
2145
};
2146
2147
//========================================================================
2148
/*!
2149
YaspLevelIterator enables iteration over entities of one grid level
2150
2151
We have specializations for codim==0 (elements) and
2152
codim=dim (vertices).
2153
The general version throws a GridError.
2154
*/
2155
//========================================================================
2156
2157
template<int codim, PartitionIteratorType pitype, class GridImp>
2158
class YaspLevelIterator :
2159
public YaspEntityPointer<codim,GridImp>
2160
{
2161
//! know your own dimension
2162
enum { dim=GridImp::dimension };
2163
//! know your own dimension of world
2164
enum { dimworld=GridImp::dimensionworld };
2165
typedef typename GridImp::ctype ctype;
2166
public:
2167
typedef typename GridImp::template Codim<codim>::Entity Entity;
2168
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2169
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2170
typedef YaspSpecialEntity<codim,dim,GridImp> SpecialEntity;
2171
2172
//! constructor
2173
YaspLevelIterator (const GridImp * yg, const YGLI & g, const TSI & it) :
2174
YaspEntityPointer<codim,GridImp>(yg,g,it) {}
2175
2176
//! copy constructor
2177
YaspLevelIterator (const YaspLevelIterator& i) :
2178
YaspEntityPointer<codim,GridImp>(i) {}
2179
2180
//! increment
2181
void increment()
2182
{
2183
++(this->_it);
2184
}
2185
};
2186
2187
2188
//========================================================================
2189
/*!
2190
\brief level-wise, non-persistent, consecutive index
2191
2192
*/
2193
//========================================================================
2194
2195
template<class GridImp>
2196
class YaspLevelIndexSet
2197
: public IndexSet< GridImp, YaspLevelIndexSet< GridImp >, unsigned int >
2198
{
2199
typedef YaspLevelIndexSet< GridImp > This;
2200
typedef IndexSet< GridImp, This, unsigned int > Base;
2201
2202
public:
2203
typedef typename Base::IndexType IndexType;
2204
2205
using Base::subIndex;
2206
2207
//! constructor stores reference to a grid and level
2208
YaspLevelIndexSet ( const GridImp &g, int l )
2209
: grid( g ),
2210
level( l )
2211
{
2212
// contains a single element type;
2213
for (int codim=0; codim<=GridImp::dimension; codim++)
2214
mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
2215
}
2216
2217
//! get index of an entity
2218
template<int cc>
2219
IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const
2220
{
2221
assert( cc == 0 || cc == GridImp::dimension );
2222
return grid.getRealImplementation(e).compressedIndex();
2223
}
2224
2225
//! get index of subentity of an entity
2226
template< int cc >
2227
IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
2228
int i, unsigned int codim ) const
2229
{
2230
assert( cc == 0 || cc == GridImp::dimension );
2231
if( cc == GridImp::dimension )
2232
return grid.getRealImplementation(e).compressedIndex();
2233
else
2234
return grid.getRealImplementation(e).subCompressedIndex(i,codim);
2235
}
2236
2237
//! get number of entities of given type and level (the level is known to the object)
2238
int size (GeometryType type) const
2239
{
2240
return grid.size( level, type );
2241
}
2242
2243
//! return size of set for a given codim
2244
int size (int codim) const
2245
{
2246
return grid.size( level, codim );
2247
}
2248
2249
//! return true if the given entity is contained in \f$E\f$.
2250
template<class EntityType>
2251
bool contains (const EntityType& e) const
2252
{
2253
return e.level() == level;
2254
}
2255
2256
//! deliver all geometry types used in this grid
2257
const std::vector<GeometryType>& geomTypes (int codim) const
2258
{
2259
return mytypes[codim];
2260
}
2261
2262
private:
2263
const GridImp& grid;
2264
int level;
2265
std::vector<GeometryType> mytypes[GridImp::dimension+1];
2266
};
2267
2268
2269
// Leaf Index Set
2270
2271
template<class GridImp>
2272
class YaspLeafIndexSet
2273
: public IndexSet< GridImp, YaspLeafIndexSet< GridImp >, unsigned int >
2274
{
2275
typedef YaspLeafIndexSet< GridImp > This;
2276
typedef IndexSet< GridImp, This, unsigned int > Base;
2277
2278
public:
2279
typedef typename Base::IndexType IndexType;
2280
2281
using Base::subIndex;
2282
2283
//! constructor stores reference to a grid
2284
explicit YaspLeafIndexSet ( const GridImp &g )
2285
: grid( g )
2286
{
2287
// contains a single element type;
2288
for (int codim=0; codim<=GridImp::dimension; codim++)
2289
mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
2290
}
2291
2292
//! get index of an entity
2293
template<int cc>
2294
IndexType index (const typename GridImp::Traits::template Codim<cc>::Entity& e) const
2295
{
2296
assert( cc == 0 || cc == GridImp::dimension );
2297
return grid.getRealImplementation(e).compressedIndex();
2298
}
2299
2300
//! get index of subentity of an entity
2301
template< int cc >
2302
IndexType subIndex ( const typename remove_const< GridImp >::type::Traits::template Codim< cc >::Entity &e,
2303
int i, unsigned int codim ) const
2304
{
2305
assert( cc == 0 || cc == GridImp::dimension );
2306
if( cc == GridImp::dimension )
2307
return grid.getRealImplementation(e).compressedIndex();
2308
else
2309
return grid.getRealImplementation(e).subCompressedIndex(i,codim);
2310
}
2311
2312
//! get number of entities of given type
2313
int size (GeometryType type) const
2314
{
2315
return grid.size( grid.maxLevel(), type );
2316
}
2317
2318
//! return size of set for a given codim
2319
int size (int codim) const
2320
{
2321
return grid.size( grid.maxLevel(), codim );
2322
}
2323
2324
//! return true if the given entity is contained in \f$E\f$.
2325
template<class EntityType>
2326
bool contains (const EntityType& e) const
2327
{
2328
//return e.isLeaf();
2329
return (e.level() == grid.maxLevel());
2330
}
2331
2332
//! deliver all geometry types used in this grid
2333
const std::vector<GeometryType>& geomTypes (int codim) const
2334
{
2335
return mytypes[codim];
2336
}
2337
2338
private:
2339
const GridImp& grid;
2340
enum { ncodim = remove_const<GridImp>::type::dimension+1 };
2341
std::vector<GeometryType> mytypes[ncodim];
2342
};
2343
2344
2345
2346
2347
//========================================================================
2348
/*!
2349
\brief persistent, globally unique Ids
2350
2351
*/
2352
//========================================================================
2353
2354
template<class GridImp>
2355
class YaspGlobalIdSet : public IdSet<GridImp,YaspGlobalIdSet<GridImp>,
2356
typename remove_const<GridImp>::type::PersistentIndexType >
2357
/*
2358
We used the remove_const to extract the Type from the mutable class,
2359
because the const class is not instantiated yet.
2360
*/
2361
{
2362
typedef YaspGlobalIdSet< GridImp > This;
2363
2364
public:
2365
//! define the type used for persisitent indices
2366
typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;
2367
2368
using IdSet<GridImp, This, IdType>::subId;
2369
2370
//! constructor stores reference to a grid
2371
explicit YaspGlobalIdSet ( const GridImp &g )
2372
: grid( g )
2373
{}
2374
2375
//! get id of an entity
2376
/*
2377
We use the remove_const to extract the Type from the mutable class,
2378
because the const class is not instantiated yet.
2379
*/
2380
template<int cd>
2381
IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
2382
{
2383
return grid.getRealImplementation(e).persistentIndex();
2384
}
2385
2386
//! get id of subentity
2387
/*
2388
We use the remove_const to extract the Type from the mutable class,
2389
because the const class is not instantiated yet.
2390
*/
2391
IdType subId (const typename remove_const<GridImp>::type::Traits::template Codim< 0 >::Entity &e,
2392
int i, unsigned int codim ) const
2393
{
2394
return grid.getRealImplementation(e).subPersistentIndex(i,codim);
2395
}
2396
2397
private:
2398
const GridImp& grid;
2399
};
2400
2401
2402
template<int dim, int dimworld>
2403
struct YaspGridFamily
2404
{
38
YaspGrid stands for yet another structured parallel grid.
39
It will implement the dune grid interface for structured grids with codim 0
40
and dim, with arbitrary overlap, parallel features with two overlap
41
models, periodic boundaries and fast a implementation allowing on-the-fly computations.
42
*/
43
44
namespace Dune {
45
46
//************************************************************************
47
/*! define name for floating point type used for coordinates in yaspgrid.
48
You can change the type for coordinates by changing this single typedef.
49
*/
50
typedef double yaspgrid_ctype;
51
52
/* some sizes for building global ids
53
*/
54
const int yaspgrid_dim_bits = 24; // bits for encoding each dimension
55
const int yaspgrid_level_bits = 6; // bits for encoding level number
56
const int yaspgrid_codim_bits = 4; // bits for encoding codimension
57
58
59
//************************************************************************
60
// forward declaration of templates
61
62
template<int dim> class YaspGrid;
63
template<int mydim, int cdim, class GridImp> class YaspGeometry;
64
template<int codim, int dim, class GridImp> class YaspEntity;
65
template<int codim, class GridImp> class YaspEntityPointer;
66
template<int codim, class GridImp> class YaspEntitySeed;
67
template<int codim, PartitionIteratorType pitype, class GridImp> class YaspLevelIterator;
68
template<class GridImp> class YaspIntersectionIterator;
69
template<class GridImp> class YaspIntersection;
70
template<class GridImp> class YaspHierarchicIterator;
71
template<class GridImp, bool isLeafIndexSet> class YaspIndexSet;
72
template<class GridImp> class YaspGlobalIdSet;
73
74
namespace FacadeOptions
75
{
76
77
template<int dim, int mydim, int cdim>
78
struct StoreGeometryReference<mydim, cdim, YaspGrid<dim>, YaspGeometry>
79
{
80
static const bool v = false;
81
};
82
83
template<int dim, int mydim, int cdim>
84
struct StoreGeometryReference<mydim, cdim, const YaspGrid<dim>, YaspGeometry>
85
{
86
static const bool v = false;
87
};
88
89
}
90
91
} // namespace Dune
92
93
#include <dune/grid/yaspgrid/yaspgridgeometry.hh>
94
#include <dune/grid/yaspgrid/yaspgridentity.hh>
95
#include <dune/grid/yaspgrid/yaspgridintersection.hh>
96
#include <dune/grid/yaspgrid/yaspgridintersectioniterator.hh>
97
#include <dune/grid/yaspgrid/yaspgridhierarchiciterator.hh>
98
#include <dune/grid/yaspgrid/yaspgridentityseed.hh>
99
#include <dune/grid/yaspgrid/yaspgridentitypointer.hh>
100
#include <dune/grid/yaspgrid/yaspgridleveliterator.hh>
101
#include <dune/grid/yaspgrid/yaspgridindexsets.hh>
102
#include <dune/grid/yaspgrid/yaspgrididset.hh>
103
104
namespace Dune {
105
106
template<int dim>
107
struct YaspGridFamily
108
{
2405
109
#if HAVE_MPI
2406
typedef CollectiveCommunication<MPI_Comm> CCType;
110
typedef CollectiveCommunication<MPI_Comm> CCType;
2407
111
#else
2408
typedef CollectiveCommunication<Dune::YaspGrid<dim> > CCType;
112
typedef CollectiveCommunication<Dune::YaspGrid<dim> > CCType;
2409
113
#endif
2410
114
2411
typedef GridTraits<dim,dimworld,Dune::YaspGrid<dim>,
2412
YaspGeometry,YaspEntity,
2413
YaspEntityPointer,YaspLevelIterator,
2414
YaspIntersection, // leaf intersection
2415
YaspIntersection, // level intersection
2416
YaspIntersectionIterator, // leaf intersection iter
2417
YaspIntersectionIterator, // level intersection iter
2418
YaspHierarchicIterator,
2419
YaspLevelIterator,
2420
YaspLevelIndexSet< const YaspGrid< dim > >,
2421
YaspLeafIndexSet< const YaspGrid< dim > >,
2422
YaspGlobalIdSet<const YaspGrid<dim> >,
2423
bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
2424
YaspGlobalIdSet<const YaspGrid<dim> >,
2425
bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
2426
CCType,
2427
DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
2428
YaspEntitySeed>
2429
Traits;
2430
};
115
typedef GridTraits<dim, // dimension of the grid
116
dim, // dimension of the world space
117
Dune::YaspGrid<dim>,
118
YaspGeometry,YaspEntity,
119
YaspEntityPointer,
120
YaspLevelIterator, // type used for the level iterator
121
YaspIntersection, // leaf intersection
122
YaspIntersection, // level intersection
123
YaspIntersectionIterator, // leaf intersection iter
124
YaspIntersectionIterator, // level intersection iter
125
YaspHierarchicIterator,
126
YaspLevelIterator, // type used for the leaf(!) iterator
127
YaspIndexSet< const YaspGrid< dim >, false >, // level index set
128
YaspIndexSet< const YaspGrid< dim >, true >, // leaf index set
129
YaspGlobalIdSet<const YaspGrid<dim> >,
130
bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
131
YaspGlobalIdSet<const YaspGrid<dim> >,
132
bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits>,
133
CCType,
134
DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
135
YaspEntitySeed>
136
Traits;
137
};
2431
138
2432
template<int dim, int codim>
2433
struct YaspCommunicateMeta {
2434
template<class G, class DataHandle>
2435
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2436
{
2437
if (data.contains(dim,codim))
139
template<int dim, int codim>
140
struct YaspCommunicateMeta {
141
template<class G, class DataHandle>
142
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
143
{
144
if (data.contains(dim,codim))
2438
145
{
2439
146
DUNE_THROW(GridError, "interface communication not implemented");
2440
147
}
2441
YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
2442
}
2443
};
2444
2445
template<int dim>
2446
struct YaspCommunicateMeta<dim,dim> {
2447
template<class G, class DataHandle>
2448
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2449
{
2450
if (data.contains(dim,dim))
2451
g.template communicateCodim<DataHandle,dim>(data,iftype,dir,level);
2452
YaspCommunicateMeta<dim,dim-1>::comm(g,data,iftype,dir,level);
2453
}
2454
};
2455
2456
template<int dim>
2457
struct YaspCommunicateMeta<dim,0> {
2458
template<class G, class DataHandle>
2459
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
2460
{
2461
if (data.contains(dim,0))
2462
g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
2463
}
2464
};
2465
2466
2467
//************************************************************************
2468
/*!
2469
\brief [<em> provides \ref Dune::Grid </em>]
2470
\brief Provides a distributed structured cube mesh.
2471
\ingroup GridImplementations
2472
2473
YaspGrid stands for yet another structured parallel grid.
2474
It implements the dune grid interface for structured grids with codim 0
2475
and dim, with arbitrary overlap (including zero),
2476
periodic boundaries and fast implementation allowing on-the-fly computations.
2477
2478
\tparam dim The dimension of the grid and its surrounding world
2479
2480
\par History:
2481
\li started on July 31, 2004 by PB based on abstractions developed in summer 2003
2482
*/
2483
template<int dim>
2484
class YaspGrid :
2485
public GridDefaultImplementation<dim,dim,yaspgrid_ctype,YaspGridFamily<dim,dim> >,
2486
public MultiYGrid<dim,yaspgrid_ctype>
2487
{
2488
typedef const YaspGrid<dim> GridImp;
2489
2490
void init()
2491
{
2492
setsizes();
2493
indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,0) );
2494
theleafindexset.push_back( new YaspLeafIndexSet<const YaspGrid<dim> >(*this) );
2495
theglobalidset.push_back( new YaspGlobalIdSet<const YaspGrid<dim> >(*this) );
2496
boundarysegmentssize();
2497
}
2498
2499
void boundarysegmentssize()
2500
{
148
YaspCommunicateMeta<dim,codim-1>::comm(g,data,iftype,dir,level);
149
}
150
};
151
152
template<int dim>
153
struct YaspCommunicateMeta<dim,dim> {
154
template<class G, class DataHandle>
155
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
156
{
157
if (data.contains(dim,dim))
158
g.template communicateCodim<DataHandle,dim>(data,iftype,dir,level);
159
YaspCommunicateMeta<dim,dim-1>::comm(g,data,iftype,dir,level);
160
}
161
};
162
163
template<int dim>
164
struct YaspCommunicateMeta<dim,0> {
165
template<class G, class DataHandle>
166
static void comm (const G& g, DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level)
167
{
168
if (data.contains(dim,0))
169
g.template communicateCodim<DataHandle,0>(data,iftype,dir,level);
170
}
171
};
172
173
//************************************************************************
174
/*!
175
\brief [<em> provides \ref Dune::Grid </em>]
176
\brief Provides a distributed structured cube mesh.
177
\ingroup GridImplementations
178
\ingroup YaspGrid
179
180
YaspGrid stands for yet another structured parallel grid.
181
It implements the dune grid interface for structured grids with codim 0
182
and dim, with arbitrary overlap (including zero),
183
periodic boundaries and fast implementation allowing on-the-fly computations.
184
185
\tparam dim The dimension of the grid and its surrounding world
186
187
\par History:
188
\li started on July 31, 2004 by PB based on abstractions developed in summer 2003
189
*/
190
template<int dim>
191
class YaspGrid
192
: public GridDefaultImplementation<dim,dim,yaspgrid_ctype,YaspGridFamily<dim> >
193
{
194
public:
195
//! Type used for coordinates
196
typedef yaspgrid_ctype ctype;
197
198
struct Intersection {
199
/** \brief The intersection as a subgrid of the local grid */
200
SubYGrid<dim,ctype> grid;
201
202
/** \brief Rank of the process where the other grid is stored */
203
int rank;
204
205
/** \brief Manhattan distance to the other grid */
206
int distance;
207
};
208
209
/** \brief A single grid level within a YaspGrid
210
*/
211
struct YGridLevel {
212
213
/** \brief Level number of this level grid */
214
int level() const
215
{
216
return level_;
217
}
218
219
// cell (codim 0) data
220
YGrid<dim,ctype> cell_global; // the whole cell grid on that level
221
SubYGrid<dim,ctype> cell_overlap; // we have no ghost cells, so our part is overlap completely
222
SubYGrid<dim,ctype> cell_interior; // interior cells are a subgrid of all cells
223
224
std::deque<Intersection> send_cell_overlap_overlap; // each intersection is a subgrid of overlap
225
std::deque<Intersection> recv_cell_overlap_overlap; // each intersection is a subgrid of overlap
226
227
std::deque<Intersection> send_cell_interior_overlap; // each intersection is a subgrid of overlap
228
std::deque<Intersection> recv_cell_overlap_interior; // each intersection is a subgrid of overlap
229
230
// vertex (codim dim) data
231
YGrid<dim,ctype> vertex_global; // the whole vertex grid on that level
232
SubYGrid<dim,ctype> vertex_overlapfront; // all our vertices are overlap and front
233
SubYGrid<dim,ctype> vertex_overlap; // subgrid containing only overlap
234
SubYGrid<dim,ctype> vertex_interiorborder; // subgrid containing only interior and border
235
SubYGrid<dim,ctype> vertex_interior; // subgrid containing only interior
236
237
std::deque<Intersection> send_vertex_overlapfront_overlapfront; // each intersection is a subgrid of overlapfront
238
std::deque<Intersection> recv_vertex_overlapfront_overlapfront; // each intersection is a subgrid of overlapfront
239
240
std::deque<Intersection> send_vertex_overlap_overlapfront; // each intersection is a subgrid of overlapfront
241
std::deque<Intersection> recv_vertex_overlapfront_overlap; // each intersection is a subgrid of overlapfront
242
243
std::deque<Intersection> send_vertex_interiorborder_interiorborder; // each intersection is a subgrid of overlapfront
244
std::deque<Intersection> recv_vertex_interiorborder_interiorborder; // each intersection is a subgrid of overlapfront
245
246
std::deque<Intersection> send_vertex_interiorborder_overlapfront; // each intersection is a subgrid of overlapfront
247
std::deque<Intersection> recv_vertex_overlapfront_interiorborder; // each intersection is a subgrid of overlapfront
248
249
// general
250
YaspGrid<dim>* mg; // each grid level knows its multigrid
251
int overlap; // in mesh cells on this level
252
/** \brief The level number within the YaspGrid level hierarchy */
253
int level_;
254
};
255
256
//! define types used for arguments
257
typedef FieldVector<int, dim> iTupel;
258
typedef FieldVector<ctype, dim> fTupel;
259
260
// communication tag used by multigrid
261
enum { tag = 17 };
262
263
//! return reference to torus
264
const Torus<dim>& torus () const
265
{
266
return _torus;
267
}
268
269
//! Iterator over the grid levels
270
typedef typename ReservedVector<YGridLevel,32>::const_iterator YGridLevelIterator;
271
272
//! return iterator pointing to coarsest level
273
YGridLevelIterator begin () const
274
{
275
return YGridLevelIterator(_levels,0);
276
}
277
278
//! return iterator pointing to given level
279
YGridLevelIterator begin (int i) const
280
{
281
if (i<0 || i>maxLevel())
282
DUNE_THROW(GridError, "level not existing");
283
return YGridLevelIterator(_levels,i);
284
}
285
286
//! return iterator pointing to one past the finest level
287
YGridLevelIterator end () const
288
{
289
return YGridLevelIterator(_levels,maxLevel()+1);
290
}
291
292
// static method to create the default load balance strategy
293
static const YLoadBalance<dim>* defaultLoadbalancer()
294
{
295
static YLoadBalance<dim> lb;
296
return & lb;
297
}
298
299
protected:
300
/** \brief Make a new YGridLevel structure
301
*
302
* \param L size of the whole domain in each direction
303
* \param s number of cells in each direction
304
* \param periodic indicate periodicity for each direction
305
* \param o_interior origin of interior (non-overlapping) cell decomposition
306
* \param s_interior size of interior cell decomposition
307
* \param overlap to be used on this grid level
308
*/
309
YGridLevel makelevel (int level, fTupel L, iTupel s, std::bitset<dim> periodic, iTupel o_interior, iTupel s_interior, int overlap)
310
{
311
// first, lets allocate a new structure
312
YGridLevel g;
313
g.overlap = overlap;
314
g.mg = this;
315
g.level_ = level;
316
317
// the global cell grid
318
iTupel o = iTupel(0); // logical origin is always 0, that is not a restriction
319
fTupel h;
320
fTupel r;
321
for (int i=0; i<dim; i++) h[i] = L[i]/s[i]; // the mesh size in each direction
322
for (int i=0; i<dim; i++) r[i] = 0.5*h[i]; // the shift for cell centers
323
g.cell_global = YGrid<dim,ctype>(o,s,h,r); // this is the global cell grid
324
325
// extend the cell interior grid by overlap considering periodicity
326
iTupel o_overlap;
327
iTupel s_overlap;
328
for (int i=0; i<dim; i++)
329
{
330
if (periodic[i])
331
{
332
// easy case, extend by 2 overlaps in total
333
o_overlap[i] = o_interior[i]-overlap; // Note: origin might be negative now
334
s_overlap[i] = s_interior[i]+2*overlap; // Note: might be larger than global size
335
}
336
else
337
{
338
// nonperiodic case, intersect with global size
339
int min = std::max(0,o_interior[i]-overlap);
340
int max = std::min(s[i]-1,o_interior[i]+s_interior[i]-1+overlap);
341
o_overlap[i] = min;
342
s_overlap[i] = max-min+1;
343
}
344
}
345
g.cell_overlap = SubYGrid<dim,ctype>(YGrid<dim,ctype>(o_overlap,s_overlap,h,r));
346
347
// now make the interior grid a subgrid of the overlapping grid
348
iTupel offset;
349
for (int i=0; i<dim; i++) offset[i] = o_interior[i]-o_overlap[i];
350
g.cell_interior = SubYGrid<dim,ctype>(o_interior,s_interior,offset,s_overlap,h,r);
351
352
// compute cell intersections
353
intersections(g.cell_overlap,g.cell_overlap,g.cell_global.size(),g.send_cell_overlap_overlap,g.recv_cell_overlap_overlap);
354
intersections(g.cell_interior,g.cell_overlap,g.cell_global.size(),g.send_cell_interior_overlap,g.recv_cell_overlap_interior);
355
356
// now we can do the vertex grids. They are derived completely from the cell grids
357
iTupel o_vertex_global, s_vertex_global;
358
for (int i=0; i<dim; i++) r[i] = 0.0; // the shift for vertices is zero, and the mesh size is same as for cells
359
360
// first let's make the global grid
361
for (int i=0; i<dim; i++) o_vertex_global[i] = g.cell_global.origin(i);
362
for (int i=0; i<dim; i++) s_vertex_global[i] = g.cell_global.size(i)+1; // one more vertices than cells ...
363
g.vertex_global = YGrid<dim,ctype>(o_vertex_global,s_vertex_global,h,r);
364
365
// now the local grid stored in this processor. All other grids are subgrids of this
366
iTupel o_vertex_overlapfront;
367
iTupel s_vertex_overlapfront;
368
for (int i=0; i<dim; i++) o_vertex_overlapfront[i] = g.cell_overlap.origin(i);
369
for (int i=0; i<dim; i++) s_vertex_overlapfront[i] = g.cell_overlap.size(i)+1; // one more vertices than cells ...
370
g.vertex_overlapfront = SubYGrid<dim,ctype>(YGrid<dim,ctype>(o_vertex_overlapfront,s_vertex_overlapfront,h,r));
371
372
// now overlap only (i.e. without front), is subgrid of overlapfront
373
iTupel o_vertex_overlap;
374
iTupel s_vertex_overlap;
375
for (int i=0; i<dim; i++)
376
{
377
o_vertex_overlap[i] = g.cell_overlap.origin(i);
378
s_vertex_overlap[i] = g.cell_overlap.size(i)+1;
379
380
if (!periodic[i] && g.cell_overlap.origin(i)>g.cell_global.origin(i))
381
{
382
// not at the lower boundary
383
o_vertex_overlap[i] += 1;
384
s_vertex_overlap[i] -= 1;
385
}
386
387
if (!periodic[i] && g.cell_overlap.origin(i)+g.cell_overlap.size(i)<g.cell_global.origin(i)+g.cell_global.size(i))
388
{
389
// not at the upper boundary
390
s_vertex_overlap[i] -= 1;
391
}
392
393
394
offset[i] = o_vertex_overlap[i]-o_vertex_overlapfront[i];
395
}
396
g.vertex_overlap = SubYGrid<dim,ctype>(o_vertex_overlap,s_vertex_overlap,offset,s_vertex_overlapfront,h,r);
397
398
// now interior with border
399
iTupel o_vertex_interiorborder;
400
iTupel s_vertex_interiorborder;
401
for (int i=0; i<dim; i++) o_vertex_interiorborder[i] = g.cell_interior.origin(i);
402
for (int i=0; i<dim; i++) s_vertex_interiorborder[i] = g.cell_interior.size(i)+1;
403
for (int i=0; i<dim; i++) offset[i] = o_vertex_interiorborder[i]-o_vertex_overlapfront[i];
404
g.vertex_interiorborder = SubYGrid<dim,ctype>(o_vertex_interiorborder,s_vertex_interiorborder,offset,s_vertex_overlapfront,h,r);
405
406
// now only interior
407
iTupel o_vertex_interior;
408
iTupel s_vertex_interior;
409
for (int i=0; i<dim; i++)
410
{
411
o_vertex_interior[i] = g.cell_interior.origin(i);
412
s_vertex_interior[i] = g.cell_interior.size(i)+1;
413
414
if (!periodic[i] && g.cell_interior.origin(i)>g.cell_global.origin(i))
415
{
416
// not at the lower boundary
417
o_vertex_interior[i] += 1;
418
s_vertex_interior[i] -= 1;
419
}
420
421
if (!periodic[i] && g.cell_interior.origin(i)+g.cell_interior.size(i)<g.cell_global.origin(i)+g.cell_global.size(i))
422
{
423
// not at the upper boundary
424
s_vertex_interior[i] -= 1;
425
}
426
427
offset[i] = o_vertex_interior[i]-o_vertex_overlapfront[i];
428
}
429
g.vertex_interior = SubYGrid<dim,ctype>(o_vertex_interior,s_vertex_interior,offset,s_vertex_overlapfront,h,r);
430
431
// compute vertex intersections
432
intersections(g.vertex_overlapfront,g.vertex_overlapfront,g.cell_global.size(),
433
g.send_vertex_overlapfront_overlapfront,g.recv_vertex_overlapfront_overlapfront);
434
intersections(g.vertex_overlap,g.vertex_overlapfront,g.cell_global.size(),
435
g.send_vertex_overlap_overlapfront,g.recv_vertex_overlapfront_overlap);
436
intersections(g.vertex_interiorborder,g.vertex_interiorborder,g.cell_global.size(),
437
g.send_vertex_interiorborder_interiorborder,g.recv_vertex_interiorborder_interiorborder);
438
intersections(g.vertex_interiorborder,g.vertex_overlapfront,g.cell_global.size(),
439
g.send_vertex_interiorborder_overlapfront,g.recv_vertex_overlapfront_interiorborder);
440
441
// return the whole thing
442
return g;
443
}
444
445
446
struct mpifriendly_ygrid {
447
mpifriendly_ygrid ()
448
: origin(0), size(0), h(0.0), r(0.0)
449
{}
450
mpifriendly_ygrid (const YGrid<dim,ctype>& grid)
451
: origin(grid.origin()), size(grid.size()), h(grid.meshsize()), r(grid.shift())
452
{}
453
iTupel origin;
454
iTupel size;
455
fTupel h;
456
fTupel r;
457
};
458
459
/** \brief Construct list of intersections with neighboring processors
460
*
461
* \param recvgrid the grid stored in this processor
462
* \param sendgrid the subgrid to be sent to neighboring processors
463
* \param size needed to shift local grid in periodic case
464
* \returns two lists: Intersections to be sent and Intersections to be received
465
* \note sendgrid/recvgrid may be SubYGrids. Since intersection method is virtual it should work properly
466
*/
467
void intersections (const SubYGrid<dim,ctype>& sendgrid, const SubYGrid<dim,ctype>& recvgrid, const iTupel& size,
468
std::deque<Intersection>& sendlist, std::deque<Intersection>& recvlist)
469
{
470
// the exchange buffers
471
std::vector<YGrid<dim,ctype> > send_recvgrid(_torus.neighbors());
472
std::vector<YGrid<dim,ctype> > recv_recvgrid(_torus.neighbors());
473
std::vector<YGrid<dim,ctype> > send_sendgrid(_torus.neighbors());
474
std::vector<YGrid<dim,ctype> > recv_sendgrid(_torus.neighbors());
475
476
// new exchange buffers to send simple struct without virtual functions
477
std::vector<mpifriendly_ygrid> mpifriendly_send_recvgrid(_torus.neighbors());
478
std::vector<mpifriendly_ygrid> mpifriendly_recv_recvgrid(_torus.neighbors());
479
std::vector<mpifriendly_ygrid> mpifriendly_send_sendgrid(_torus.neighbors());
480
std::vector<mpifriendly_ygrid> mpifriendly_recv_sendgrid(_torus.neighbors());
481
482
// fill send buffers; iterate over all neighboring processes
483
// non-periodic case is handled automatically because intersection will be zero
484
for (typename Torus<dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
485
{
486
// determine if we communicate with this neighbor (and what)
487
bool skip = false;
488
iTupel coord = _torus.coord(); // my coordinates
489
iTupel delta = i.delta(); // delta to neighbor
490
iTupel nb = coord; // the neighbor
491
for (int k=0; k<dim; k++) nb[k] += delta[k];
492
iTupel v = iTupel(0); // grid movement
493
494
for (int k=0; k<dim; k++)
495
{
496
if (nb[k]<0)
497
{
498
if (_periodic[k])
499
v[k] += size[k];
500
else
501
skip = true;
502
}
503
if (nb[k]>=_torus.dims(k))
504
{
505
if (_periodic[k])
506
v[k] -= size[k];
507
else
508
skip = true;
509
}
510
// neither might be true, then v=0
511
}
512
513
// store moved grids in send buffers
514
if (!skip)
515
{
516
send_sendgrid[i.index()] = sendgrid.move(v);
517
send_recvgrid[i.index()] = recvgrid.move(v);
518
}
519
else
520
{
521
send_sendgrid[i.index()] = YGrid<dim,ctype>(iTupel(0),iTupel(0),fTupel(0.0),fTupel(0.0));
522
send_recvgrid[i.index()] = YGrid<dim,ctype>(iTupel(0),iTupel(0),fTupel(0.0),fTupel(0.0));
523
}
524
}
525
526
// issue send requests for sendgrid being sent to all neighbors
527
for (typename Torus<dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
528
{
529
mpifriendly_send_sendgrid[i.index()] = mpifriendly_ygrid(send_sendgrid[i.index()]);
530
_torus.send(i.rank(), &mpifriendly_send_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
531
}
532
533
// issue recv requests for sendgrids of neighbors
534
for (typename Torus<dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
535
_torus.recv(i.rank(), &mpifriendly_recv_sendgrid[i.index()], sizeof(mpifriendly_ygrid));
536
537
// exchange the sendgrids
538
_torus.exchange();
539
540
// issue send requests for recvgrid being sent to all neighbors
541
for (typename Torus<dim>::ProcListIterator i=_torus.sendbegin(); i!=_torus.sendend(); ++i)
542
{
543
mpifriendly_send_recvgrid[i.index()] = mpifriendly_ygrid(send_recvgrid[i.index()]);
544
_torus.send(i.rank(), &mpifriendly_send_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
545
}
546
547
// issue recv requests for recvgrid of neighbors
548
for (typename Torus<dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
549
_torus.recv(i.rank(), &mpifriendly_recv_recvgrid[i.index()], sizeof(mpifriendly_ygrid));
550
551
// exchange the recvgrid
552
_torus.exchange();
553
554
// process receive buffers and compute intersections
555
for (typename Torus<dim>::ProcListIterator i=_torus.recvbegin(); i!=_torus.recvend(); ++i)
556
{
557
// what must be sent to this neighbor
558
Intersection send_intersection;
559
mpifriendly_ygrid yg = mpifriendly_recv_recvgrid[i.index()];
560
recv_recvgrid[i.index()] = YGrid<dim,ctype>(yg.origin,yg.size,yg.h,yg.r);
561
send_intersection.grid = sendgrid.intersection(recv_recvgrid[i.index()]);
562
send_intersection.rank = i.rank();
563
send_intersection.distance = i.distance();
564
if (!send_intersection.grid.empty()) sendlist.push_front(send_intersection);
565
566
Intersection recv_intersection;
567
yg = mpifriendly_recv_sendgrid[i.index()];
568
recv_sendgrid[i.index()] = YGrid<dim,ctype>(yg.origin,yg.size,yg.h,yg.r);
569
recv_intersection.grid = recvgrid.intersection(recv_sendgrid[i.index()]);
570
recv_intersection.rank = i.rank();
571
recv_intersection.distance = i.distance();
572
if(!recv_intersection.grid.empty()) recvlist.push_back(recv_intersection);
573
}
574
}
575
576
protected:
577
578
typedef const YaspGrid<dim> GridImp;
579
580
void init()
581
{
582
setsizes();
583
indexsets.push_back( make_shared< YaspIndexSet<const YaspGrid<dim>, false > >(*this,0) );
584
boundarysegmentssize();
585
}
586
587
void boundarysegmentssize()
588
{
2501
589
// sizes of local macro grid
2502
const FieldVector<int, dim> & size = MultiYGrid<dim,ctype>::begin().cell_overlap().size();
2503
FieldVector<int, dim> sides;
590
const FieldVector<int, dim> & size = begin()->cell_overlap.size();
591
Dune::array<int, dim> sides;
2504
592
{
2505
for (int i=0; i<dim; i++)
2506
{
2507
sides[i] =
2508
((MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i)
2509
== MultiYGrid<dim,ctype>::begin().cell_global().origin(i))+
2510
(MultiYGrid<dim,ctype>::begin().cell_overlap().origin(i) +
2511
MultiYGrid<dim,ctype>::begin().cell_overlap().size(i)
2512
== MultiYGrid<dim,ctype>::begin().cell_global().origin(i) +
2513
MultiYGrid<dim,ctype>::begin().cell_global().size(i)));
2514
}
593
for (int i=0; i<dim; i++)
594
{
595
sides[i] =
596
((begin()->cell_overlap.origin(i) == begin()->cell_global.origin(i))+
597
(begin()->cell_overlap.origin(i) + begin()->cell_overlap.size(i)
598
== begin()->cell_global.origin(i) + begin()->cell_global.size(i)));
599
}
2515
600
}
2516
601
nBSegments = 0;
2517
602
for (int k=0; k<dim; k++)
2518
603
{
2519
int offset = 1;
2520
for (int l=0; l<dim; l++)
2521
{
2522
if (l==k) continue;
2523
offset *= size[l];
2524
}
2525
nBSegments += sides[k]*offset;
604
int offset = 1;
605
for (int l=0; l<dim; l++)
606
{
607
if (l==k) continue;
608
offset *= size[l];
609
}
610
nBSegments += sides[k]*offset;
2526
611
}
2527
}
2528
2529
public:
2530
2531
using MultiYGrid<dim,yaspgrid_ctype>::defaultLoadbalancer;
2532
2533
//! define type used for coordinates in grid module
2534
typedef yaspgrid_ctype ctype;
2535
2536
// define the persistent index type
2537
typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits> PersistentIndexType;
2538
2539
//! the GridFamily of this grid
2540
typedef YaspGridFamily<dim,dim> GridFamily;
2541
// the Traits
2542
typedef typename YaspGridFamily<dim,dim>::Traits Traits;
2543
2544
// need for friend declarations in entity
2545
typedef YaspLevelIndexSet<YaspGrid<dim> > LevelIndexSetType;
2546
typedef YaspLeafIndexSet<YaspGrid<dim> > LeafIndexSetType;
2547
typedef YaspGlobalIdSet<YaspGrid<dim> > GlobalIdSetType;
2548
2549
//! maximum number of levels allowed
2550
enum { MAXL=64 };
2551
2552
//! shorthand for base class data types
2553
typedef MultiYGrid<dim,ctype> YMG;
2554
typedef typename MultiYGrid<dim,ctype>::YGridLevelIterator YGLI;
2555
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
2556
typedef typename MultiYGrid<dim,ctype>::Intersection IS;
2557
typedef typename std::deque<IS>::const_iterator ISIT;
2558
2559
/*! Constructor for a YaspGrid, they are all forwarded to the base class
2560
@param comm MPI communicator where this mesh is distributed to
2561
@param L extension of the domain
2562
@param s number of cells on coarse mesh in each direction
2563
@param periodic tells if direction is periodic or not
2564
@param overlap size of overlap on coarsest grid (same in all directions)
2565
@param lb pointer to an overloaded YLoadBalance instance
2566
*/
2567
YaspGrid (Dune::MPIHelper::MPICommunicator comm,
2568
Dune::FieldVector<ctype, dim> L,
2569
Dune::FieldVector<int, dim> s,
2570
Dune::FieldVector<bool, dim> periodic, int overlap,
2571
const YLoadBalance<dim>* lb = defaultLoadbalancer())
2572
#if HAVE_MPI
2573
: YMG(comm,L,s,periodic,overlap,lb), ccobj(comm),
2574
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2575
#else
2576
: YMG(L,s,periodic,overlap,lb),
2577
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2578
#endif
2579
{
2580
init();
2581
}
2582
2583
2584
/*! Constructor for a sequential YaspGrid, they are all forwarded to the base class.
2585
2586
Sequential here means that the whole grid is living on one process even if your program is running
2587
in parallel.
2588
@see YaspGrid(Dune::MPIHelper::MPICommunicator, Dune::FieldVector<ctype, dim>, Dune::FieldVector<int, dim>, Dune::FieldVector<bool, dim>, int)
2589
for constructing one parallel grid decomposed between the processors.
2590
@param L extension of the domain
2591
@param s number of cells on coarse mesh in each direction
2592
@param periodic tells if direction is periodic or not
2593
@param overlap size of overlap on coarsest grid (same in all directions)
2594
@param lb pointer to an overloaded YLoadBalance instance
2595
*/
2596
YaspGrid (Dune::FieldVector<ctype, dim> L,
2597
Dune::FieldVector<int, dim> s,
2598
Dune::FieldVector<bool, dim> periodic, int overlap,
2599
const YLoadBalance<dim>* lb = YMG::defaultLoadbalancer())
2600
#if HAVE_MPI
2601
: YMG(MPI_COMM_SELF,L,s,periodic,overlap,lb), ccobj(MPI_COMM_SELF),
2602
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2603
#else
2604
: YMG(L,s,periodic,overlap,lb),
2605
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
2606
#endif
2607
{
2608
init();
2609
}
2610
2611
~YaspGrid()
2612
{
2613
deallocatePointers(indexsets);
2614
deallocatePointers(theleafindexset);
2615
deallocatePointers(theglobalidset);
2616
}
612
}
613
614
public:
615
616
// define the persistent index type
617
typedef bigunsignedint<dim*yaspgrid_dim_bits+yaspgrid_level_bits+yaspgrid_codim_bits> PersistentIndexType;
618
619
//! the GridFamily of this grid
620
typedef YaspGridFamily<dim> GridFamily;
621
// the Traits
622
typedef typename YaspGridFamily<dim>::Traits Traits;
623
624
// need for friend declarations in entity
625
typedef YaspIndexSet<YaspGrid<dim>, false > LevelIndexSetType;
626
typedef YaspIndexSet<YaspGrid<dim>, true > LeafIndexSetType;
627
typedef YaspGlobalIdSet<YaspGrid<dim> > GlobalIdSetType;
628
629
//! shorthand for some data types
630
typedef typename SubYGrid<dim,ctype>::TransformingSubIterator TSI;
631
typedef typename std::deque<Intersection>::const_iterator ISIT;
632
633
//! The constructor of the old MultiYGrid class
634
void MultiYGridSetup (
635
fTupel L, iTupel s, std::bitset<dim> periodic, int overlap, const YLoadBalance<dim>* lb = defaultLoadbalancer())
636
{
637
_LL = L;
638
_s = s;
639
_periodic = periodic;
640
_levels.resize(1);
641
_overlap = overlap;
642
643
// coarse cell interior grid obtained through partitioning of global grid
644
#if HAVE_MPI
645
iTupel o_interior;
646
iTupel s_interior;
647
iTupel o = iTupel(0);
648
array<int,dim> sArray;
649
std::copy(s.begin(), s.end(), sArray.begin());
650
double imbal = _torus.partition(_torus.rank(),o,sArray,o_interior,s_interior);
651
imbal = _torus.global_max(imbal);
652
#else
653
iTupel o = iTupel(0);
654
iTupel o_interior(o);
655
iTupel s_interior(s);
656
#endif
657
// add level
658
_levels[0] = makelevel(0,L,s,periodic,o_interior,s_interior,overlap);
659
}
660
661
//! The constructor of the old MultiYGrid class
662
void MultiYGridSetup (
663
fTupel L,
664
Dune::array<int,dim> s,
665
std::bitset<dim> periodic,
666
int overlap,
667
const YLoadBalance<dim>* lb = defaultLoadbalancer())
668
{
669
_LL = L;
670
_periodic = periodic;
671
_levels.resize(1);
672
_overlap = overlap;
673
674
std::copy(s.begin(), s.end(), this->_s.begin());
675
676
// coarse cell interior grid obtained through partitioning of global grid
677
iTupel o = iTupel(0);
678
iTupel o_interior(o);
679
iTupel s_interior;
680
std::copy(s.begin(), s.end(), s_interior.begin());
681
#if HAVE_MPI
682
double imbal = _torus.partition(_torus.rank(),o,s,o_interior,s_interior);
683
imbal = _torus.global_max(imbal);
684
#endif
685
686
// add level
687
_levels[0] = makelevel(0,L,_s,periodic,o_interior,s_interior,overlap);
688
}
689
690
/*! Constructor
691
@param comm MPI communicator where this mesh is distributed to
692
@param L extension of the domain
693
@param s number of cells on coarse mesh in each direction
694
@param periodic tells if direction is periodic or not
695
@param overlap size of overlap on coarsest grid (same in all directions)
696
@param lb pointer to an overloaded YLoadBalance instance
697
698
\deprecated Will be removed after dune-grid 2.3.
699
Use the corresponding constructor taking array<int> and std::bitset instead.
700
*/
701
YaspGrid (Dune::MPIHelper::MPICommunicator comm,
702
Dune::FieldVector<ctype, dim> L,
703
Dune::FieldVector<int, dim> s,
704
Dune::FieldVector<bool, dim> periodic, int overlap,
705
const YLoadBalance<dim>* lb = defaultLoadbalancer())
706
DUNE_DEPRECATED_MSG("Use the corresponding constructor taking array<int> and std::bitset")
707
#if HAVE_MPI
708
: ccobj(comm),
709
_torus(comm,tag,s,lb),
710
#else
711
: _torus(tag,s,lb),
712
#endif
713
leafIndexSet_(*this),
714
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
715
{
716
MultiYGridSetup(L,s,std::bitset<dim>(),overlap,lb);
717
718
// hack: copy input bitfield (in FieldVector<bool>) into std::bitset
719
for (size_t i=0; i<dim; i++)
720
this->_periodic[i] = periodic[i];
721
init();
722
}
723
724
725
/*! Constructor for a sequential YaspGrid
726
727
Sequential here means that the whole grid is living on one process even if your program is running
728
in parallel.
729
@see YaspGrid(Dune::MPIHelper::MPICommunicator, Dune::FieldVector<ctype, dim>, Dune::FieldVector<int, dim>, Dune::FieldVector<bool, dim>, int)
730
for constructing one parallel grid decomposed between the processors.
731
@param L extension of the domain
732
@param s number of cells on coarse mesh in each direction
733
@param periodic tells if direction is periodic or not
734
@param overlap size of overlap on coarsest grid (same in all directions)
735
@param lb pointer to an overloaded YLoadBalance instance
736
737
\deprecated Will be removed after dune-grid 2.3.
738
Use the corresponding constructor taking array<int> and std::bitset instead.
739
*/
740
YaspGrid (Dune::FieldVector<ctype, dim> L,
741
Dune::FieldVector<int, dim> s,
742
Dune::FieldVector<bool, dim> periodic, int overlap,
743
const YLoadBalance<dim>* lb = defaultLoadbalancer())
744
DUNE_DEPRECATED_MSG("Use the corresponding constructor taking array<int> and std::bitset")
745
#if HAVE_MPI
746
: ccobj(MPI_COMM_SELF),
747
_torus(MPI_COMM_SELF,tag,s,lb),
748
#else
749
: _torus(tag,s,lb),
750
#endif
751
leafIndexSet_(*this),
752
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
753
{
754
MultiYGridSetup(L,s,std::bitset<dim>(),overlap,lb);
755
756
// hack: copy input bitfield (in FieldVector<bool>) into std::bitset
757
for (size_t i=0; i<dim; i++)
758
this->_periodic[i] = periodic[i];
759
init();
760
}
761
762
/*! Constructor
763
@param comm MPI communicator where this mesh is distributed to
764
@param L extension of the domain
765
@param s number of cells on coarse mesh in each direction
766
@param periodic tells if direction is periodic or not
767
@param overlap size of overlap on coarsest grid (same in all directions)
768
@param lb pointer to an overloaded YLoadBalance instance
769
*/
770
YaspGrid (Dune::MPIHelper::MPICommunicator comm,
771
Dune::FieldVector<ctype, dim> L,
772
Dune::array<int, dim> s,
773
std::bitset<dim> periodic,
774
int overlap,
775
const YLoadBalance<dim>* lb = defaultLoadbalancer())
776
#if HAVE_MPI
777
: ccobj(comm),
778
_torus(comm,tag,s,lb),
779
#else
780
: _torus(tag,s,lb),
781
#endif
782
leafIndexSet_(*this),
783
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
784
{
785
MultiYGridSetup(L,s,periodic,overlap,lb);
786
787
init();
788
}
789
790
791
/*! Constructor for a sequential YaspGrid
792
793
Sequential here means that the whole grid is living on one process even if your program is running
794
in parallel.
795
@see YaspGrid(Dune::MPIHelper::MPICommunicator, Dune::FieldVector<ctype, dim>, Dune::FieldVector<int, dim>, Dune::FieldVector<bool, dim>, int)
796
for constructing one parallel grid decomposed between the processors.
797
@param L extension of the domain
798
@param s number of cells on coarse mesh in each direction
799
@param periodic tells if direction is periodic or not
800
@param overlap size of overlap on coarsest grid (same in all directions)
801
@param lb pointer to an overloaded YLoadBalance instance
802
*/
803
YaspGrid (Dune::FieldVector<ctype, dim> L,
804
Dune::array<int, dim> s,
805
std::bitset<dim> periodic,
806
int overlap,
807
const YLoadBalance<dim>* lb = defaultLoadbalancer())
808
#if HAVE_MPI
809
: ccobj(MPI_COMM_SELF),
810
_torus(MPI_COMM_SELF,tag,s,lb),
811
#else
812
: _torus(tag,s,lb),
813
#endif
814
leafIndexSet_(*this),
815
keep_ovlp(true), adaptRefCount(0), adaptActive(false)
816
{
817
MultiYGridSetup(L,s,periodic,overlap,lb);
818
819
init();
820
}
821
822
/*! Constructor for a sequential YaspGrid without periodicity
823
824
Sequential here means that the whole grid is living on one process even if your program is running
825
in parallel.
826
@see YaspGrid(Dune::MPIHelper::MPICommunicator, Dune::FieldVector<ctype, dim>, Dune::FieldVector<int, dim>, Dune::FieldVector<bool, dim>, int)
827
for constructing one parallel grid decomposed between the processors.
828
@param L extension of the domain (lower left is always (0,...,0)
829
@param elements number of cells on coarse mesh in each direction
830
*/
831
YaspGrid (Dune::FieldVector<ctype, dim> L,
832
Dune::array<int, dim> elements)
833
#if HAVE_MPI
834
: ccobj(MPI_COMM_SELF),
835
_torus(MPI_COMM_SELF,tag,elements,defaultLoadbalancer()),
836
#else
837
: _torus(tag,elements,defaultLoadbalancer()),
838
#endif
839
leafIndexSet_(*this),
840
_LL(L),
841
_overlap(0),
842
keep_ovlp(true),
843
adaptRefCount(0), adaptActive(false)
844
{
845
_levels.resize(1);
846
847
std::copy(elements.begin(), elements.end(), _s.begin());
848
849
// coarse cell interior grid obtained through partitioning of global grid
850
iTupel o = iTupel(0);
851
iTupel o_interior(o);
852
iTupel s_interior;
853
std::copy(elements.begin(), elements.end(), s_interior.begin());
854
#if HAVE_MPI
855
double imbal = _torus.partition(_torus.rank(),o,elements,o_interior,s_interior);
856
imbal = _torus.global_max(imbal);
857
#endif
858
859
// add level
860
_levels[0] = makelevel(0,L,_s,_periodic,o_interior,s_interior,0);
861
862
init();
863
}
2617
864
2618
865
private:
2619
// do not copy this class
2620
YaspGrid(const YaspGrid&);
866
// do not copy this class
867
YaspGrid(const YaspGrid&);
2621
868
2622
869
public:
2623
2624
/*! Return maximum level defined in this grid. Levels are numbered
2625
0 ... maxlevel with 0 the coarsest level.
2626
*/
2627
int maxLevel() const {return MultiYGrid<dim,ctype>::maxlevel();} // delegate
2628
2629
//! refine the grid refCount times. What about overlap?
2630
void globalRefine (int refCount)
2631
{
2632
if (refCount < -maxLevel())
2633
DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
2634
"Coarsening " << -refCount << " levels requested!");
2635
for (int k=refCount; k<0; k++)
870
871
/*! Return maximum level defined in this grid. Levels are numbered
872
0 ... maxlevel with 0 the coarsest level.
873
*/
874
int maxLevel() const
875
{
876
return _levels.size()-1;
877
}
878
879
//! refine the grid refCount times. What about overlap?
880
void globalRefine (int refCount)
881
{
882
if (refCount < -maxLevel())
883
DUNE_THROW(GridError, "Only " << maxLevel() << " levels left. " <<
884
"Coarsening " << -refCount << " levels requested!");
885
886
// If refCount is negative then coarsen the grid
887
for (int k=refCount; k<0; k++)
2636
888
{
2637
MultiYGrid<dim,ctype>::coarsen();
889
// create an empty grid level
890
YGridLevel empty;
891
_levels.back() = empty;
892
// reduce maxlevel
893
_levels.pop_back();
894
2638
895
setsizes();
2639
896
indexsets.pop_back();
2640
897
}
2641
for (int k=0; k<refCount; k++)
898
899
// If refCount is positive refine the grid
900
for (int k=0; k<refCount; k++)
2642
901
{
2643
MultiYGrid<dim,ctype>::refine(keep_ovlp);
902
// access to coarser grid level
903
YGridLevel& cg = _levels[maxLevel()];
904
905
// compute size of new global grid
906
iTupel s;
907
for (int i=0; i<dim; i++)
908
s[i] = 2*cg.cell_global.size(i);
909
910
// compute overlap
911
int overlap = (keep_ovlp) ? 2*cg.overlap : cg.overlap;
912
913
// the cell interior grid obtained from coarse cell interior grid
914
iTupel o_interior;
915
iTupel s_interior;
916
for (int i=0; i<dim; i++)
917
o_interior[i] = 2*cg.cell_interior.origin(i);
918
for (int i=0; i<dim; i++)
919
s_interior[i] = 2*cg.cell_interior.size(i);
920
921
// add level
922
_levels.push_back( makelevel(_levels.size(),_LL,s,_periodic,o_interior,s_interior,overlap) );
923
2644
924
setsizes();
2645
indexsets.push_back( new YaspLevelIndexSet<const YaspGrid<dim> >(*this,maxLevel()) );
925
indexsets.push_back( make_shared<YaspIndexSet<const YaspGrid<dim>, false > >(*this,maxLevel()) );
2646
926
}
2647
}
2648
2649
/**
2650
\brief set options for refinement
2651
@param keepPhysicalOverlap [true] keep the physical size of the overlap, [false] keep the number of cells in the overlap. Default is [true].
2652
*/
2653
void refineOptions (bool keepPhysicalOverlap)
2654
{
2655
keep_ovlp = keepPhysicalOverlap;
2656
}
2657
2658
/** \brief Marks an entity to be refined/coarsened in a subsequent adapt.
2659
2660
\param[in] refCount Number of subdivisions that should be applied. Negative value means coarsening.
2661
\param[in] e Entity to Entity that should be refined
2662
2663
\return true if Entity was marked, false otherwise.
2664
2665
\note
2666
- On yaspgrid marking one element will mark all other elements of the level aswell
2667
- If refCount is lower than refCount of a previous mark-call, nothing is changed
2668
*/
2669
bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
2670
{
2671
assert(adaptActive == false);
2672
if (e.level() != maxLevel()) return false;
2673
adaptRefCount = std::max(adaptRefCount, refCount);
2674
return true;
2675
}
2676
2677
/** \brief returns adaptation mark for given entity
2678
2679
\param[in] e Entity for which adaptation mark should be determined
2680
2681
\return int adaptation mark, here the default value 0 is returned
2682
*/
2683
int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
2684
{
2685
return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
2686
}
2687
2688
//! map adapt to global refine
2689
bool adapt ()
2690
{
2691
globalRefine(adaptRefCount);
2692
return (adaptRefCount > 0);
2693
}
2694
2695
//! returns true, if the grid will be coarsened
2696
bool preAdapt ()
2697
{
2698
adaptActive = true;
2699
adaptRefCount = comm().max(adaptRefCount);
2700
return (adaptRefCount < 0);
2701
}
2702
2703
//! clean up some markers
2704
void postAdapt()
2705
{
2706
adaptActive = false;
2707
adaptRefCount = 0;
2708
}
2709
2710
//! one past the end on this level
2711
template<int cd, PartitionIteratorType pitype>
2712
typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
927
}
928
929
/**
930
\brief set options for refinement
931
@param keepPhysicalOverlap [true] keep the physical size of the overlap, [false] keep the number of cells in the overlap. Default is [true].
932
*/
933
void refineOptions (bool keepPhysicalOverlap)
934
{
935
keep_ovlp = keepPhysicalOverlap;
936
}
937
938
/** \brief Marks an entity to be refined/coarsened in a subsequent adapt.
939
940
\param[in] refCount Number of subdivisions that should be applied. Negative value means coarsening.
941
\param[in] e Entity to Entity that should be refined
942
943
\return true if Entity was marked, false otherwise.
944
945
\note
946
- On yaspgrid marking one element will mark all other elements of the level as well
947
- If refCount is lower than refCount of a previous mark-call, nothing is changed
948
*/
949
bool mark( int refCount, const typename Traits::template Codim<0>::Entity & e )
950
{
951
assert(adaptActive == false);
952
if (e.level() != maxLevel()) return false;
953
adaptRefCount = std::max(adaptRefCount, refCount);
954
return true;
955
}
956
957
/** \brief returns adaptation mark for given entity
958
959
\param[in] e Entity for which adaptation mark should be determined
960
961
\return int adaptation mark, here the default value 0 is returned
962
*/
963
int getMark ( const typename Traits::template Codim<0>::Entity &e ) const
964
{
965
return ( e.level() == maxLevel() ) ? adaptRefCount : 0;
966
}
967
968
//! map adapt to global refine
969
bool adapt ()
970
{
971
globalRefine(adaptRefCount);
972
return (adaptRefCount > 0);
973
}
974
975
//! returns true, if the grid will be coarsened
976
bool preAdapt ()
977
{
978
adaptActive = true;
979
adaptRefCount = comm().max(adaptRefCount);
980
return (adaptRefCount < 0);
981
}
982
983
//! clean up some markers
984
void postAdapt()
985
{
986
adaptActive = false;
987
adaptRefCount = 0;
988
}
989
990
//! one past the end on this level
991
template<int cd, PartitionIteratorType pitype>
992
typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const
2713
993
{
2714
994
return levelbegin<cd,pitype>(level);
2715
995
}
2716
996
2717
//! Iterator to one past the last entity of given codim on level for partition type
2718
template<int cd, PartitionIteratorType pitype>
2719
typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
997
//! Iterator to one past the last entity of given codim on level for partition type
998
template<int cd, PartitionIteratorType pitype>
999
typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const
2720
1000
{
2721
1001
return levelend<cd,pitype>(level);
2722
1002
}
2723
1003
2724
//! version without second template parameter for convenience
2725
template<int cd>
2726
typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
1004
//! version without second template parameter for convenience
1005
template<int cd>
1006
typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
2727
1007
{
2728
1008
return levelbegin<cd,All_Partition>(level);
2729
1009
}
2730
1010
2731
//! version without second template parameter for convenience
2732
template<int cd>
2733
typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
1011
//! version without second template parameter for convenience
1012
template<int cd>
1013
typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
2734
1014
{
2735
1015
return levelend<cd,All_Partition>(level);
2736
1016
}
2737
1017
2738
//! return LeafIterator which points to the first entity in maxLevel
2739
template<int cd, PartitionIteratorType pitype>
2740
typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const
1018
//! return LeafIterator which points to the first entity in maxLevel
1019
template<int cd, PartitionIteratorType pitype>
1020
typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const
2741
1021
{
2742
1022
return levelbegin<cd,pitype>(maxLevel());
2743
1023
}
2744
1024
2745
//! return LeafIterator which points behind the last entity in maxLevel
2746
template<int cd, PartitionIteratorType pitype>
2747
typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const
1025
//! return LeafIterator which points behind the last entity in maxLevel
1026
template<int cd, PartitionIteratorType pitype>
1027
typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const
2748
1028
{
2749
1029
return levelend<cd,pitype>(maxLevel());
2750
1030
}
2751
1031
2752
//! return LeafIterator which points to the first entity in maxLevel
2753
template<int cd>
2754
typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
1032
//! return LeafIterator which points to the first entity in maxLevel
1033
template<int cd>
1034
typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
2755
1035
{
2756
1036
return levelbegin<cd,All_Partition>(maxLevel());
2757
1037
}
2758
1038
2759
//! return LeafIterator which points behind the last entity in maxLevel
2760
template<int cd>
2761
typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
1039
//! return LeafIterator which points behind the last entity in maxLevel
1040
template<int cd>
1041
typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
2762
1042
{
2763
1043
return levelend<cd,All_Partition>(maxLevel());
2764
1044
}
2765
1045
2766
// \brief obtain EntityPointer from EntitySeed. */
2767
template <typename Seed>
2768
typename Traits::template Codim<Seed::codimension>::EntityPointer
2769
entityPointer(const Seed& seed) const
2770
{
2771
const int codim = Seed::codimension;
2772
YGLI g = MultiYGrid<dim,ctype>::begin(seed.level());
2773
switch (codim)
1046
// \brief obtain EntityPointer from EntitySeed. */
1047
template <typename Seed>
1048
typename Traits::template Codim<Seed::codimension>::EntityPointer
1049
entityPointer(const Seed& seed) const
2774
1050
{
2775
case 0:
2776
return YaspEntityPointer<codim,GridImp>(this,g,
2777
TSI(g.cell_overlap(), seed.coord()));
2778
case dim:
2779
return YaspEntityPointer<codim,GridImp>(this,g,
2780
TSI(g.vertex_overlap(), seed.coord()));
2781
default:
1051
const int codim = Seed::codimension;
1052
YGridLevelIterator g = begin(this->getRealImplementation(seed).level());
1053
switch (codim)
1054
{
1055
case 0 :
1056
return YaspEntityPointer<codim,GridImp>(this,g,
1057
TSI(g->cell_overlap, this->getRealImplementation(seed).coord()));
1058
case dim :
1059
return YaspEntityPointer<codim,GridImp>(this,g,
1060
TSI(g->vertex_overlap, this->getRealImplementation(seed).coord()));
1061
default :
2782
1062
DUNE_THROW(GridError, "YaspEntityPointer: codim not implemented");
2783
}
2784
}
2785
2786
//! return size (= distance in graph) of overlap region
2787
int overlapSize (int level, int codim) const
2788
{
2789
YGLI g = MultiYGrid<dim,ctype>::begin(level);
2790
return g.overlap();
2791
}
2792
2793
//! return size (= distance in graph) of overlap region
2794
int overlapSize (int codim) const
2795
{
2796
YGLI g = MultiYGrid<dim,ctype>::begin(maxLevel());
2797
return g.overlap();
2798
}
2799
2800
//! return size (= distance in graph) of ghost region
2801
int ghostSize (int level, int codim) const
2802
{
2803
return 0;
2804
}
2805
2806
//! return size (= distance in graph) of ghost region
2807
int ghostSize (int codim) const
2808
{
2809
return 0;
2810
}
2811
2812
//! number of entities per level and codim in this process
2813
int size (int level, int codim) const
2814
{
1063
}
1064
}
1065
1066
//! return size (= distance in graph) of overlap region
1067
int overlapSize (int level, int codim) const
1068
{
1069
YGridLevelIterator g = begin(level);
1070
return g->overlap;
1071
}
1072
1073
//! return size (= distance in graph) of overlap region
1074
int overlapSize (int codim) const
1075
{
1076
YGridLevelIterator g = begin(maxLevel());
1077
return g->overlap;
1078
}
1079
1080
//! return size (= distance in graph) of ghost region
1081
int ghostSize (int level, int codim) const
1082
{
1083
return 0;
1084
}
1085
1086
//! return size (= distance in graph) of ghost region
1087
int ghostSize (int codim) const
1088
{
1089
return 0;
1090
}
1091
1092
//! number of entities per level and codim in this process
1093
int size (int level, int codim) const
1094
{
2815
1095
return sizes[level][codim];
2816
}
1096
}
2817
1097
2818
//! number of leaf entities per codim in this process
2819
int size (int codim) const
2820
{
1098
//! number of leaf entities per codim in this process
1099
int size (int codim) const
1100
{
2821
1101
return sizes[maxLevel()][codim];
2822
}
1102
}
2823
1103
2824
//! number of entities per level and geometry type in this process
2825
int size (int level, GeometryType type) const
2826
{
1104
//! number of entities per level and geometry type in this process
1105
int size (int level, GeometryType type) const
1106
{
2827
1107
return (type.isCube()) ? sizes[level][dim-type.dim()] : 0;
2828
}
1108
}
2829
1109
2830
//! number of leaf entities per geometry type in this process
2831
int size (GeometryType type) const
2832
{
1110
//! number of leaf entities per geometry type in this process
1111
int size (GeometryType type) const
1112
{
2833
1113
return size(maxLevel(),type);
2834
}
2835
2836
//! \brief returns the number of boundary segments within the macro grid
2837
size_t numBoundarySegments () const
2838
{
2839
return nBSegments;
2840
}
2841
2842
/*! The new communication interface
2843
2844
communicate objects for all codims on a given level
2845
*/
2846
template<class DataHandleImp, class DataType>
2847
void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir, int level) const
2848
{
2849
YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
2850
}
2851
2852
/*! The new communication interface
2853
2854
communicate objects for all codims on the leaf grid
2855
*/
2856
template<class DataHandleImp, class DataType>
2857
void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir) const
2858
{
2859
YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
2860
}
2861
2862
/*! The new communication interface
2863
2864
communicate objects for one codim
2865
*/
2866
template<class DataHandle, int codim>
2867
void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
2868
{
2869
// check input
2870
if (!data.contains(dim,codim)) return; // should have been checked outside
2871
2872
// data types
2873
typedef typename DataHandle::DataType DataType;
2874
2875
// access to grid level
2876
YGLI g = MultiYGrid<dim,ctype>::begin(level);
2877
2878
// find send/recv lists or throw error
2879
const std::deque<IS>* sendlist=0;
2880
const std::deque<IS>* recvlist=0;
2881
if (codim==0) // the elements
1114
}
1115
1116
//! \brief returns the number of boundary segments within the macro grid
1117
size_t numBoundarySegments () const
1118
{
1119
return nBSegments;
1120
}
1121
1122
/*! The new communication interface
1123
1124
communicate objects for all codims on a given level
1125
*/
1126
template<class DataHandleImp, class DataType>
1127
void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir, int level) const
1128
{
1129
YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,level);
1130
}
1131
1132
/*! The new communication interface
1133
1134
communicate objects for all codims on the leaf grid
1135
*/
1136
template<class DataHandleImp, class DataType>
1137
void communicate (CommDataHandleIF<DataHandleImp,DataType> & data, InterfaceType iftype, CommunicationDirection dir) const
1138
{
1139
YaspCommunicateMeta<dim,dim>::comm(*this,data,iftype,dir,this->maxLevel());
1140
}
1141
1142
/*! The new communication interface
1143
1144
communicate objects for one codim
1145
*/
1146
template<class DataHandle, int codim>
1147
void communicateCodim (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
1148
{
1149
// check input
1150
if (!data.contains(dim,codim)) return; // should have been checked outside
1151
1152
// data types
1153
typedef typename DataHandle::DataType DataType;
1154
1155
// access to grid level
1156
YGridLevelIterator g = begin(level);
1157
1158
// find send/recv lists or throw error
1159
const std::deque<Intersection>* sendlist=0;
1160
const std::deque<Intersection>* recvlist=0;
1161
if (codim==0) // the elements
2882
1162
{
2883
1163
if (iftype==InteriorBorder_InteriorBorder_Interface)
2884
1164
return; // there is nothing to do in this case
2885
1165
if (iftype==InteriorBorder_All_Interface)
2886
{
2887
sendlist = &g.send_cell_interior_overlap();
2888
recvlist = &g.recv_cell_overlap_interior();
2889
}
1166
{
1167
sendlist = &g->send_cell_interior_overlap;
1168
recvlist = &g->recv_cell_overlap_interior;
1169
}
2890
1170
if (iftype==Overlap_OverlapFront_Interface || iftype==Overlap_All_Interface || iftype==All_All_Interface)
2891
{
2892
sendlist = &g.send_cell_overlap_overlap();
2893
recvlist = &g.recv_cell_overlap_overlap();
2894
}
1171
{
1172
sendlist = &g->send_cell_overlap_overlap;
1173
recvlist = &g->recv_cell_overlap_overlap;
1174
}
2895
1175
}
2896
if (codim==dim) // the vertices
1176
if (codim==dim) // the vertices
2897
1177
{
2898
1178
if (iftype==InteriorBorder_InteriorBorder_Interface)
2899
{
2900
sendlist = &g.send_vertex_interiorborder_interiorborder();
2901
recvlist = &g.recv_vertex_interiorborder_interiorborder();
2902
}
1179
{
1180
sendlist = &g->send_vertex_interiorborder_interiorborder;
1181
recvlist = &g->recv_vertex_interiorborder_interiorborder;
1182
}
2903
1183
2904
1184
if (iftype==InteriorBorder_All_Interface)
2905
{
2906
sendlist = &g.send_vertex_interiorborder_overlapfront();
2907
recvlist = &g.recv_vertex_overlapfront_interiorborder();
2908
}
1185
{
1186
sendlist = &g->send_vertex_interiorborder_overlapfront;
1187
recvlist = &g->recv_vertex_overlapfront_interiorborder;
1188
}
2909
1189
if (iftype==Overlap_OverlapFront_Interface || iftype==Overlap_All_Interface)
2910
{
2911
sendlist = &g.send_vertex_overlap_overlapfront();
2912
recvlist = &g.recv_vertex_overlapfront_overlap();
2913
}
1190
{
1191
sendlist = &g->send_vertex_overlap_overlapfront;
1192
recvlist = &g->recv_vertex_overlapfront_overlap;
1193
}
2914
1194
if (iftype==All_All_Interface)
2915
{
2916
sendlist = &g.send_vertex_overlapfront_overlapfront();
2917
recvlist = &g.recv_vertex_overlapfront_overlapfront();
2918
}
1195
{
1196
sendlist = &g->send_vertex_overlapfront_overlapfront;
1197
recvlist = &g->recv_vertex_overlapfront_overlapfront;
1198
}
2919
1199
}
2920
1200
2921
// change communication direction?
2922
if (dir==BackwardCommunication)
2923
std::swap(sendlist,recvlist);
2924
2925
int cnt;
2926
2927
// Size computation (requires communication if variable size)
2928
std::vector<int> send_size(sendlist->size(),-1); // map rank to total number of objects (of type DataType) to be sent
2929
std::vector<int> recv_size(recvlist->size(),-1); // map rank to total number of objects (of type DataType) to be recvd
2930
std::vector<size_t*> send_sizes(sendlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be sent
2931
std::vector<size_t*> recv_sizes(recvlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be recvd
2932
if (data.fixedsize(dim,codim))
1201
// change communication direction?
1202
if (dir==BackwardCommunication)
1203
std::swap(sendlist,recvlist);
1204
1205
int cnt;
1206
1207
// Size computation (requires communication if variable size)
1208
std::vector<int> send_size(sendlist->size(),-1); // map rank to total number of objects (of type DataType) to be sent
1209
std::vector<int> recv_size(recvlist->size(),-1); // map rank to total number of objects (of type DataType) to be recvd
1210
std::vector<size_t*> send_sizes(sendlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be sent
1211
std::vector<size_t*> recv_sizes(recvlist->size(),static_cast<size_t*>(0)); // map rank to array giving number of objects per entity to be recvd
1212
if (data.fixedsize(dim,codim))
2933
1213
{
2934
1214
// fixed size: just take a dummy entity, size can be computed without communication
2935
1215
cnt=0;
2936
1216
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
2937
{
2938
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
2939
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2940
send_size[cnt] = is->grid.totalsize() * data.size(*it);
2941
cnt++;
2942
}
1217
{
1218
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1219
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1220
send_size[cnt] = is->grid.totalsize() * data.size(*it);
1221
cnt++;
1222
}
2943
1223
cnt=0;
2944
1224
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
2945
{
2946
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
2947
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2948
recv_size[cnt] = is->grid.totalsize() * data.size(*it);
2949
cnt++;
2950
}
1225
{
1226
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1227
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1228
recv_size[cnt] = is->grid.totalsize() * data.size(*it);
1229
cnt++;
1230
}
2951
1231
}
2952
else
1232
else
2953
1233
{
2954
1234
// variable size case: sender side determines the size
2955
1235
cnt=0;
2956
1236
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
1237
{
1238
// allocate send buffer for sizes per entitiy
1239
size_t *buf = new size_t[is->grid.totalsize()];
1240
send_sizes[cnt] = buf;
1241
1242
// loop over entities and ask for size
1243
int i=0; size_t n=0;
1244
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1245
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1246
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1247
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
1248
for ( ; it!=tsubend; ++it)
2957
1249
{
2958
// allocate send buffer for sizes per entitiy
2959
size_t *buf = new size_t[is->grid.totalsize()];
2960
send_sizes[cnt] = buf;
2961
2962
// loop over entities and ask for size
2963
int i=0; size_t n=0;
2964
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
2965
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
2966
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
2967
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
2968
for ( ; it!=tsubend; ++it)
2969
{
2970
buf[i] = data.size(*it);
2971
n += buf[i];
2972
i++;
2973
}
2974
2975
// now we know the size for this rank
2976
send_size[cnt] = n;
2977
2978
// hand over send request to torus class
2979
MultiYGrid<dim,ctype>::torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
2980
cnt++;
1250
buf[i] = data.size(*it);
1251
n += buf[i];
1252
i++;
2981
1253
}
2982
1254
1255
// now we know the size for this rank
1256
send_size[cnt] = n;
1257
1258
// hand over send request to torus class
1259
torus().send(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1260
cnt++;
1261
}
1262
2983
1263
// allocate recv buffers for sizes and store receive request
2984
1264
cnt=0;
2985
1265
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
2986
{
2987
// allocate recv buffer
2988
size_t *buf = new size_t[is->grid.totalsize()];
2989
recv_sizes[cnt] = buf;
2990
2991
// hand over recv request to torus class
2992
MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
2993
cnt++;
2994
}
1266
{
1267
// allocate recv buffer
1268
size_t *buf = new size_t[is->grid.totalsize()];
1269
recv_sizes[cnt] = buf;
1270
1271
// hand over recv request to torus class
1272
torus().recv(is->rank,buf,is->grid.totalsize()*sizeof(size_t));
1273
cnt++;
1274
}
2995
1275
2996
1276
// exchange all size buffers now
2997
MultiYGrid<dim,ctype>::torus().exchange();
1277
torus().exchange();
2998
1278
2999
1279
// release send size buffers
3000
1280
cnt=0;
3001
1281
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3002
{
3003
delete[] send_sizes[cnt];
3004
send_sizes[cnt] = 0;
3005
cnt++;
3006
}
1282
{
1283
delete[] send_sizes[cnt];
1284
send_sizes[cnt] = 0;
1285
cnt++;
1286
}
3007
1287
3008
1288
// process receive size buffers
3009
1289
cnt=0;
3010
1290
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3011
{
3012
// get recv buffer
3013
size_t *buf = recv_sizes[cnt];
3014
3015
// compute total size
3016
size_t n=0;
3017
for (int i=0; i<is->grid.totalsize(); ++i)
3018
n += buf[i];
3019
3020
// ... and store it
3021
recv_size[cnt] = n;
3022
++cnt;
3023
}
1291
{
1292
// get recv buffer
1293
size_t *buf = recv_sizes[cnt];
1294
1295
// compute total size
1296
size_t n=0;
1297
for (int i=0; i<is->grid.totalsize(); ++i)
1298
n += buf[i];
1299
1300
// ... and store it
1301
recv_size[cnt] = n;
1302
++cnt;
1303
}
3024
1304
}
3025
1305
3026
1306
3027
// allocate & fill the send buffers & store send request
3028
std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
3029
cnt=0;
3030
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
1307
// allocate & fill the send buffers & store send request
1308
std::vector<DataType*> sends(sendlist->size(), static_cast<DataType*>(0)); // store pointers to send buffers
1309
cnt=0;
1310
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3031
1311
{
3032
// std::cout << "[" << this->comm().rank() << "] "
3033
// << " send " << " dest=" << is->rank
3034
// << " size=" << send_size[cnt]
3035
// << std::endl;
3036
3037
1312
// allocate send buffer
3038
1313
DataType *buf = new DataType[send_size[cnt]];
3039
1314
3045
1320
3046
1321
// fill send buffer; iterate over cells in intersection
3047
1322
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3048
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1323
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3049
1324
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3050
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
1325
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3051
1326
for ( ; it!=tsubend; ++it)
3052
1327
data.gather(mb,*it);
3053
1328
3054
1329
// hand over send request to torus class
3055
MultiYGrid<dim,ctype>::torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
1330
torus().send(is->rank,buf,send_size[cnt]*sizeof(DataType));
3056
1331
cnt++;
3057
1332
}
3058
1333
3059
// allocate recv buffers and store receive request
3060
std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
3061
cnt=0;
3062
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
1334
// allocate recv buffers and store receive request
1335
std::vector<DataType*> recvs(recvlist->size(),static_cast<DataType*>(0)); // store pointers to send buffers
1336
cnt=0;
1337
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3063
1338
{
3064
// std::cout << "[" << this->comm().rank() << "] "
3065
// << " recv " << " src=" << is->rank
3066
// << " size=" << recv_size[cnt]
3067
// << std::endl;
3068
3069
1339
// allocate recv buffer
3070
1340
DataType *buf = new DataType[recv_size[cnt]];
3071
1341
3073
1343
recvs[cnt] = buf;
3074
1344
3075
1345
// hand over recv request to torus class
3076
MultiYGrid<dim,ctype>::torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
3077
cnt++;
3078
}
3079
3080
// exchange all buffers now
3081
MultiYGrid<dim,ctype>::torus().exchange();
3082
3083
// release send buffers
3084
cnt=0;
3085
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
3086
{
3087
delete[] sends[cnt];
3088
sends[cnt] = 0;
3089
cnt++;
3090
}
3091
3092
// process receive buffers and delete them
3093
cnt=0;
3094
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
1346
torus().recv(is->rank,buf,recv_size[cnt]*sizeof(DataType));
1347
cnt++;
1348
}
1349
1350
// exchange all buffers now
1351
torus().exchange();
1352
1353
// release send buffers
1354
cnt=0;
1355
for (ISIT is=sendlist->begin(); is!=sendlist->end(); ++is)
1356
{
1357
delete[] sends[cnt];
1358
sends[cnt] = 0;
1359
cnt++;
1360
}
1361
1362
// process receive buffers and delete them
1363
cnt=0;
1364
for (ISIT is=recvlist->begin(); is!=recvlist->end(); ++is)
3095
1365
{
3096
1366
// get recv buffer
3097
1367
DataType *buf = recvs[cnt];
3101
1371
3102
1372
// copy data from receive buffer; iterate over cells in intersection
3103
1373
if (data.fixedsize(dim,codim))
3104
{
3105
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3106
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3107
size_t n=data.size(*it);
3108
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3109
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3110
for ( ; it!=tsubend; ++it)
3111
data.scatter(mb,*it,n);
3112
}
1374
{
1375
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1376
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1377
size_t n=data.size(*it);
1378
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1379
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
1380
for ( ; it!=tsubend; ++it)
1381
data.scatter(mb,*it,n);
1382
}
3113
1383
else
3114
{
3115
int i=0;
3116
size_t *sbuf = recv_sizes[cnt];
3117
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3118
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
3119
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
3120
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
3121
for ( ; it!=tsubend; ++it)
3122
data.scatter(mb,*it,sbuf[i++]);
3123
delete[] sbuf;
3124
}
1384
{
1385
int i=0;
1386
size_t *sbuf = recv_sizes[cnt];
1387
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1388
it(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubbegin()));
1389
typename Traits::template Codim<codim>::template Partition<All_Partition>::LevelIterator
1390
tsubend(YaspLevelIterator<codim,All_Partition,GridImp>(this,g,is->grid.tsubend()));
1391
for ( ; it!=tsubend; ++it)
1392
data.scatter(mb,*it,sbuf[i++]);
1393
delete[] sbuf;
1394
}
3125
1395
3126
1396
// delete buffer
3127
1397
delete[] buf; // hier krachts !
3128
1398
cnt++;
3129
1399
}
3130
}
3131
3132
// The new index sets from DDM 11.07.2005
3133
const typename Traits::GlobalIdSet& globalIdSet() const
3134
{
3135
return *(theglobalidset[0]);
3136
}
3137
3138
const typename Traits::LocalIdSet& localIdSet() const
3139
{
3140
return *(theglobalidset[0]);
3141
}
3142
3143
const typename Traits::LevelIndexSet& levelIndexSet(int level) const
3144
{
3145
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3146
return *(indexsets[level]);
3147
}
3148
3149
const typename Traits::LeafIndexSet& leafIndexSet() const
3150
{
3151
return *(theleafindexset[0]);
3152
}
3153
3154
#if HAVE_MPI
3155
/*! @brief return a collective communication object
3156
*/
3157
const CollectiveCommunication<MPI_Comm>& comm () const
3158
{
3159
return ccobj;
3160
}
3161
#else
3162
/*! @brief return a collective communication object
3163
*/
3164
const CollectiveCommunication<YaspGrid>& comm () const
3165
{
3166
return ccobj;
3167
}
3168
#endif
3169
3170
YaspIntersectionIterator<const YaspGrid<dim> >&
3171
getRealIntersectionIterator(typename Traits::LevelIntersectionIterator& it)
3172
{
3173
return this->getRealImplementation(it);
3174
}
3175
3176
const YaspIntersectionIterator<const YaspGrid<dim> >&
3177
getRealIntersectionIterator(const typename Traits::LevelIntersectionIterator& it) const
3178
{
3179
return this->getRealImplementation(it);
3180
}
3181
3182
private:
3183
3184
#if HAVE_MPI
3185
CollectiveCommunication<MPI_Comm> ccobj;
3186
#else
3187
CollectiveCommunication<YaspGrid> ccobj;
3188
#endif
3189
3190
std::vector<YaspLevelIndexSet<const YaspGrid<dim> >*> indexsets;
3191
std::vector<YaspLeafIndexSet<const YaspGrid<dim> >*> theleafindexset;
3192
std::vector<YaspGlobalIdSet<const YaspGrid<dim> >*> theglobalidset;
3193
int nBSegments;
3194
3195
// Index classes need access to the real entity
3196
friend class Dune::YaspLevelIndexSet<const Dune::YaspGrid<dim> >;
3197
friend class Dune::YaspLeafIndexSet<const Dune::YaspGrid<dim> >;
3198
friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim> >;
3199
3200
friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim> >;
3201
friend class Dune::YaspIntersection<const Dune::YaspGrid<dim> >;
3202
friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim> >;
3203
3204
template<class T>
3205
void deallocatePointers(T& container)
3206
{
3207
typedef typename T::iterator Iterator;
3208
3209
for(Iterator entry=container.begin(); entry != container.end(); ++entry)
3210
delete (*entry);
3211
}
3212
3213
template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
3214
friend class Entity;
3215
3216
template<class DT>
3217
class MessageBuffer {
3218
public:
3219
// Constructor
3220
MessageBuffer (DT *p)
3221
{
3222
a=p;
3223
i=0;
3224
j=0;
3225
}
3226
3227
// write data to message buffer, acts like a stream !
3228
template<class Y>
3229
void write (const Y& data)
3230
{
3231
dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
3232
a[i++] = data;
3233
}
3234
3235
// read data from message buffer, acts like a stream !
3236
template<class Y>
3237
void read (Y& data) const
3238
{
3239
dune_static_assert(( is_same<DT,Y>::value ), "DataType missmatch");
3240
data = a[j++];
3241
}
1400
}
1401
1402
// The new index sets from DDM 11.07.2005
1403
const typename Traits::GlobalIdSet& globalIdSet() const
1404
{
1405
return theglobalidset;
1406
}
1407
1408
const typename Traits::LocalIdSet& localIdSet() const
1409
{
1410
return theglobalidset;
1411
}
1412
1413
const typename Traits::LevelIndexSet& levelIndexSet(int level) const
1414
{
1415
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1416
return *(indexsets[level]);
1417
}
1418
1419
const typename Traits::LeafIndexSet& leafIndexSet() const
1420
{
1421
return leafIndexSet_;
1422
}
1423
1424
#if HAVE_MPI
1425
/*! @brief return a collective communication object
1426
*/
1427
const CollectiveCommunication<MPI_Comm>& comm () const
1428
{
1429
return ccobj;
1430
}
1431
#else
1432
/*! @brief return a collective communication object
1433
*/
1434
const CollectiveCommunication<YaspGrid>& comm () const
1435
{
1436
return ccobj;
1437
}
1438
#endif
3242
1439
3243
1440
private:
3244
DT *a;
3245
int i;
3246
mutable int j;
3247
};
3248
3249
void setsizes ()
3250
{
3251
for (YGLI g=MultiYGrid<dim,ctype>::begin(); g!=MultiYGrid<dim,ctype>::end(); ++g)
1441
1442
// number of boundary segments of the level 0 grid
1443
int nBSegments;
1444
1445
// Index classes need access to the real entity
1446
friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim>, true >;
1447
friend class Dune::YaspIndexSet<const Dune::YaspGrid<dim>, false >;
1448
friend class Dune::YaspGlobalIdSet<const Dune::YaspGrid<dim> >;
1449
1450
friend class Dune::YaspIntersectionIterator<const Dune::YaspGrid<dim> >;
1451
friend class Dune::YaspIntersection<const Dune::YaspGrid<dim> >;
1452
friend class Dune::YaspEntity<0, dim, const Dune::YaspGrid<dim> >;
1453
1454
template <int codim_, class GridImp_>
1455
friend class Dune::YaspEntityPointer;
1456
1457
template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
1458
friend class Entity;
1459
1460
template<class DT>
1461
class MessageBuffer {
1462
public:
1463
// Constructor
1464
MessageBuffer (DT *p)
1465
{
1466
a=p;
1467
i=0;
1468
j=0;
1469
}
1470
1471
// write data to message buffer, acts like a stream !
1472
template<class Y>
1473
void write (const Y& data)
1474
{
1475
dune_static_assert(( is_same<DT,Y>::value ), "DataType mismatch");
1476
a[i++] = data;
1477
}
1478
1479
// read data from message buffer, acts like a stream !
1480
template<class Y>
1481
void read (Y& data) const
1482
{
1483
dune_static_assert(( is_same<DT,Y>::value ), "DataType mismatch");
1484
data = a[j++];
1485
}
1486
1487
private:
1488
DT *a;
1489
int i;
1490
mutable int j;
1491
};
1492
1493
void setsizes ()
1494
{
1495
for (YGridLevelIterator g=begin(); g!=end(); ++g)
3252
1496
{
3253
1497
// codim 0 (elements)
3254
sizes[g.level()][0] = 1;
3255
for (int i=0; i<dim; ++i)
3256
sizes[g.level()][0] *= g.cell_overlap().size(i);
1498
sizes[g->level()][0] = 1;
1499
for (int i=0; i<dim; ++i)
1500
sizes[g->level()][0] *= g->cell_overlap.size(i);
3257
1501
3258
1502
// codim 1 (faces)
3259
1503
if (dim>1)
1504
{
1505
sizes[g->level()][1] = 0;
1506
for (int i=0; i<dim; ++i)
3260
1507
{
3261
sizes[g.level()][1] = 0;
3262
for (int i=0; i<dim; ++i)
3263
{
3264
int s=g.cell_overlap().size(i)+1;
3265
for (int j=0; j<dim; ++j)
3266
if (j!=i)
3267
s *= g.cell_overlap().size(j);
3268
sizes[g.level()][1] += s;
3269
}
1508
int s=g->cell_overlap.size(i)+1;
1509
for (int j=0; j<dim; ++j)
1510
if (j!=i)
1511
s *= g->cell_overlap.size(j);
1512
sizes[g->level()][1] += s;
3270
1513
}
1514
}
3271
1515
3272
1516
// codim dim-1 (edges)
3273
1517
if (dim>2)
1518
{
1519
sizes[g->level()][dim-1] = 0;
1520
for (int i=0; i<dim; ++i)
3274
1521
{
3275
sizes[g.level()][dim-1] = 0;
3276
for (int i=0; i<dim; ++i)
3277
{
3278
int s=g.cell_overlap().size(i);
3279
for (int j=0; j<dim; ++j)
3280
if (j!=i)
3281
s *= g.cell_overlap().size(j)+1;
3282
sizes[g.level()][dim-1] += s;
3283
}
1522
int s=g->cell_overlap.size(i);
1523
for (int j=0; j<dim; ++j)
1524
if (j!=i)
1525
s *= g->cell_overlap.size(j)+1;
1526
sizes[g->level()][dim-1] += s;
3284
1527
}
1528
}
3285
1529
3286
1530
// codim dim (vertices)
3287
sizes[g.level()][dim] = 1;
3288
for (int i=0; i<dim; ++i)
3289
sizes[g.level()][dim] *= g.vertex_overlapfront().size(i);
3290
}
3291
}
3292
3293
//! one past the end on this level
3294
template<int cd, PartitionIteratorType pitype>
3295
YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
3296
{
3297
dune_static_assert( cd == dim || cd == 0 ,
3298
"YaspGrid only supports Entities with codim=dim and codim=0");
3299
YGLI g = MultiYGrid<dim,ctype>::begin(level);
3300
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3301
if (pitype==Ghost_Partition)
3302
return levelend <cd, pitype> (level);
3303
if (cd==0) // the elements
3304
{
3305
if (pitype<=InteriorBorder_Partition)
3306
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubbegin());
3307
if (pitype<=All_Partition)
3308
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubbegin());
3309
}
3310
if (cd==dim) // the vertices
3311
{
3312
if (pitype==Interior_Partition)
3313
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubbegin());
3314
if (pitype==InteriorBorder_Partition)
3315
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubbegin());
3316
if (pitype==Overlap_Partition)
3317
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubbegin());
3318
if (pitype<=All_Partition)
3319
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubbegin());
3320
}
3321
DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
3322
}
3323
3324
//! Iterator to one past the last entity of given codim on level for partition type
3325
template<int cd, PartitionIteratorType pitype>
3326
YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
3327
{
3328
dune_static_assert( cd == dim || cd == 0 ,
3329
"YaspGrid only supports Entities with codim=dim and codim=0");
3330
YGLI g = MultiYGrid<dim,ctype>::begin(level);
3331
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
3332
if (cd==0) // the elements
3333
{
3334
if (pitype<=InteriorBorder_Partition)
3335
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_interior().tsubend());
3336
if (pitype<=All_Partition || pitype == Ghost_Partition)
3337
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.cell_overlap().tsubend());
3338
}
3339
if (cd==dim) // the vertices
3340
{
3341
if (pitype==Interior_Partition)
3342
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interior().tsubend());
3343
if (pitype==InteriorBorder_Partition)
3344
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_interiorborder().tsubend());
3345
if (pitype==Overlap_Partition)
3346
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlap().tsubend());
3347
if (pitype<=All_Partition || pitype == Ghost_Partition)
3348
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g.vertex_overlapfront().tsubend());
3349
}
3350
DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
3351
}
3352
3353
int sizes[MAXL][dim+1]; // total number of entities per level and codim
3354
bool keep_ovlp;
3355
int adaptRefCount;
3356
bool adaptActive;
3357
};
3358
3359
namespace Capabilities
3360
{
3361
3362
/** \struct hasEntity
3363
\ingroup YaspGrid
3364
*/
3365
3366
/** \struct hasBackupRestoreFacilities
3367
\ingroup YaspGrid
3368
*/
3369
3370
/** \brief YaspGrid has only one geometry type for codim 0 entities
3371
\ingroup YaspGrid
3372
*/
3373
template<int dim>
3374
struct hasSingleGeometryType< YaspGrid<dim> >
3375
{
3376
static const bool v = true;
3377
static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
3378
};
3379
3380
/** \brief YaspGrid is a Cartesian grid
3381
\ingroup YaspGrid
3382
*/
3383
template<int dim>
3384
struct isCartesian< YaspGrid<dim> >
3385
{
3386
static const bool v = true;
3387
};
3388
3389
/** \brief YaspGrid has codim=0 entities (elements)
3390
\ingroup YaspGrid
3391
*/
3392
template<int dim>
3393
struct hasEntity< YaspGrid<dim>, 0 >
3394
{
3395
static const bool v = true;
3396
};
3397
3398
/** \brief YaspGrid has codim=dim entities (vertices)
3399
\ingroup YaspGrid
3400
*/
3401
template<int dim>
3402
struct hasEntity< YaspGrid<dim>, dim >
3403
{
3404
static const bool v = true;
3405
};
3406
3407
template< int dim >
3408
struct canCommunicate< YaspGrid< dim >, 0 >
3409
{
3410
static const bool v = true;
3411
};
3412
3413
template< int dim >
3414
struct canCommunicate< YaspGrid< dim >, dim >
3415
{
3416
static const bool v = true;
3417
};
3418
3419
/** \brief YaspGrid is parallel
3420
\ingroup YaspGrid
3421
*/
3422
template<int dim>
3423
struct isParallel< YaspGrid<dim> >
3424
{
3425
static const bool v = true;
3426
};
3427
3428
/** \brief YaspGrid is levelwise conforming
3429
\ingroup YaspGrid
3430
*/
3431
template<int dim>
3432
struct isLevelwiseConforming< YaspGrid<dim> >
3433
{
3434
static const bool v = true;
3435
};
3436
3437
/** \brief YaspGrid is leafwise conforming
3438
\ingroup YaspGrid
3439
*/
3440
template<int dim>
3441
struct isLeafwiseConforming< YaspGrid<dim> >
3442
{
3443
static const bool v = true;
3444
};
3445
3446
}
1531
sizes[g->level()][dim] = 1;
1532
for (int i=0; i<dim; ++i)
1533
sizes[g->level()][dim] *= g->vertex_overlapfront.size(i);
1534
}
1535
}
1536
1537
//! one past the end on this level
1538
template<int cd, PartitionIteratorType pitype>
1539
YaspLevelIterator<cd,pitype,GridImp> levelbegin (int level) const
1540
{
1541
dune_static_assert( cd == dim || cd == 0 ,
1542
"YaspGrid only supports Entities with codim=dim and codim=0");
1543
YGridLevelIterator g = begin(level);
1544
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1545
if (pitype==Ghost_Partition)
1546
return levelend <cd, pitype> (level);
1547
if (cd==0) // the elements
1548
{
1549
if (pitype<=InteriorBorder_Partition)
1550
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->cell_interior.tsubbegin());
1551
if (pitype<=All_Partition)
1552
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->cell_overlap.tsubbegin());
1553
}
1554
if (cd==dim) // the vertices
1555
{
1556
if (pitype==Interior_Partition)
1557
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_interior.tsubbegin());
1558
if (pitype==InteriorBorder_Partition)
1559
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_interiorborder.tsubbegin());
1560
if (pitype==Overlap_Partition)
1561
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_overlap.tsubbegin());
1562
if (pitype<=All_Partition)
1563
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_overlapfront.tsubbegin());
1564
}
1565
DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1566
}
1567
1568
//! Iterator to one past the last entity of given codim on level for partition type
1569
template<int cd, PartitionIteratorType pitype>
1570
YaspLevelIterator<cd,pitype,GridImp> levelend (int level) const
1571
{
1572
dune_static_assert( cd == dim || cd == 0 ,
1573
"YaspGrid only supports Entities with codim=dim and codim=0");
1574
YGridLevelIterator g = begin(level);
1575
if (level<0 || level>maxLevel()) DUNE_THROW(RangeError, "level out of range");
1576
if (cd==0) // the elements
1577
{
1578
if (pitype<=InteriorBorder_Partition)
1579
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->cell_interior.tsubend());
1580
if (pitype<=All_Partition || pitype == Ghost_Partition)
1581
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->cell_overlap.tsubend());
1582
}
1583
if (cd==dim) // the vertices
1584
{
1585
if (pitype==Interior_Partition)
1586
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_interior.tsubend());
1587
if (pitype==InteriorBorder_Partition)
1588
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_interiorborder.tsubend());
1589
if (pitype==Overlap_Partition)
1590
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_overlap.tsubend());
1591
if (pitype<=All_Partition || pitype == Ghost_Partition)
1592
return YaspLevelIterator<cd,pitype,GridImp>(this,g,g->vertex_overlapfront.tsubend());
1593
}
1594
DUNE_THROW(GridError, "YaspLevelIterator with this codim or partition type not implemented");
1595
}
1596
1597
#if HAVE_MPI
1598
CollectiveCommunication<MPI_Comm> ccobj;
1599
#else
1600
CollectiveCommunication<YaspGrid> ccobj;
1601
#endif
1602
1603
Torus<dim> _torus;
1604
1605
std::vector< shared_ptr< YaspIndexSet<const YaspGrid<dim>, false > > > indexsets;
1606
YaspIndexSet<const YaspGrid<dim>, true> leafIndexSet_;
1607
YaspGlobalIdSet<const YaspGrid<dim> > theglobalidset;
1608
1609
fTupel _LL;
1610
iTupel _s;
1611
std::bitset<dim> _periodic;
1612
ReservedVector<YGridLevel,32> _levels;
1613
int _overlap;
1614
int sizes[32][dim+1]; // total number of entities per level and codim
1615
bool keep_ovlp;
1616
int adaptRefCount;
1617
bool adaptActive;
1618
};
1619
1620
//! Output operator for multigrids
1621
1622
template <int d>
1623
inline std::ostream& operator<< (std::ostream& s, YaspGrid<d>& grid)
1624
{
1625
int rank = grid.torus().rank();
1626
1627
s << "[" << rank << "]:" << " YaspGrid maxlevel=" << grid.maxLevel() << std::endl;
1628
1629
for (typename YaspGrid<d>::YGridLevelIterator g=grid.begin(); g!=grid.end(); ++g)
1630
{
1631
s << "[" << rank << "]: " << std::endl;
1632
s << "[" << rank << "]: " << "==========================================" << std::endl;
1633
s << "[" << rank << "]: " << "level=" << g->level() << std::endl;
1634
s << "[" << rank << "]: " << "cell_global=" << g->cell_global << std::endl;
1635
s << "[" << rank << "]: " << "cell_overlap=" << g->cell_overlap << std::endl;
1636
s << "[" << rank << "]: " << "cell_interior=" << g->cell_interior << std::endl;
1637
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_cell_overlap_overlap.begin();
1638
i!=g->send_cell_overlap_overlap.end(); ++i)
1639
{
1640
s << "[" << rank << "]: " << " s_c_o_o "
1641
<< i->rank << " " << i->grid << std::endl;
1642
}
1643
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_cell_overlap_overlap.begin();
1644
i!=g->recv_cell_overlap_overlap.end(); ++i)
1645
{
1646
s << "[" << rank << "]: " << " r_c_o_o "
1647
<< i->rank << " " << i->grid << std::endl;
1648
}
1649
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_cell_interior_overlap.begin();
1650
i!=g->send_cell_interior_overlap.end(); ++i)
1651
{
1652
s << "[" << rank << "]: " << " s_c_i_o "
1653
<< i->rank << " " << i->grid << std::endl;
1654
}
1655
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_cell_overlap_interior.begin();
1656
i!=g->recv_cell_overlap_interior.end(); ++i)
1657
{
1658
s << "[" << rank << "]: " << " r_c_o_i "
1659
<< i->rank << " " << i->grid << std::endl;
1660
}
1661
1662
s << "[" << rank << "]: " << "-----------------------------------------------" << std::endl;
1663
s << "[" << rank << "]: " << "vertex_global=" << g->vertex_global << std::endl;
1664
s << "[" << rank << "]: " << "vertex_overlapfront=" << g->vertex_overlapfront << std::endl;
1665
s << "[" << rank << "]: " << "vertex_overlap=" << g->vertex_overlap << std::endl;
1666
s << "[" << rank << "]: " << "vertex_interiorborder=" << g->vertex_interiorborder << std::endl;
1667
s << "[" << rank << "]: " << "vertex_interior=" << g->vertex_interior << std::endl;
1668
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_vertex_overlapfront_overlapfront.begin();
1669
i!=g->send_vertex_overlapfront_overlapfront.end(); ++i)
1670
{
1671
s << "[" << rank << "]: " << " s_v_of_of "
1672
<< i->rank << " " << i->grid << std::endl;
1673
}
1674
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_vertex_overlapfront_overlapfront.begin();
1675
i!=g->recv_vertex_overlapfront_overlapfront.end(); ++i)
1676
{
1677
s << "[" << rank << "]: " << " r_v_of_of "
1678
<< i->rank << " " << i->grid << std::endl;
1679
}
1680
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_vertex_overlap_overlapfront.begin();
1681
i!=g->send_vertex_overlap_overlapfront.end(); ++i)
1682
{
1683
s << "[" << rank << "]: " << " s_v_o_of "
1684
<< i->rank << " " << i->grid << std::endl;
1685
}
1686
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_vertex_overlapfront_overlap.begin();
1687
i!=g->recv_vertex_overlapfront_overlap.end(); ++i)
1688
{
1689
s << "[" << rank << "]: " << " r_v_of_o "
1690
<< i->rank << " " << i->grid << std::endl;
1691
}
1692
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_vertex_interiorborder_interiorborder.begin();
1693
i!=g->send_vertex_interiorborder_interiorborder.end(); ++i)
1694
{
1695
s << "[" << rank << "]: " << " s_v_ib_ib "
1696
<< i->rank << " " << i->grid << std::endl;
1697
}
1698
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_vertex_interiorborder_interiorborder.begin();
1699
i!=g->recv_vertex_interiorborder_interiorborder.end(); ++i)
1700
{
1701
s << "[" << rank << "]: " << " r_v_ib_ib "
1702
<< i->rank << " " << i->grid << std::endl;
1703
}
1704
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->send_vertex_interiorborder_overlapfront.begin();
1705
i!=g->send_vertex_interiorborder_overlapfront.end(); ++i)
1706
{
1707
s << "[" << rank << "]: " << " s_v_ib_of "
1708
<< i->rank << " " << i->grid << std::endl;
1709
}
1710
for (typename std::deque<typename YaspGrid<d>::Intersection>::const_iterator i=g->recv_vertex_overlapfront_interiorborder.begin();
1711
i!=g->recv_vertex_overlapfront_interiorborder.end(); ++i)
1712
{
1713
s << "[" << rank << "]: " << " s_v_of_ib "
1714
<< i->rank << " " << i->grid << std::endl;
1715
}
1716
}
1717
1718
s << std::endl;
1719
1720
return s;
1721
}
1722
1723
namespace Capabilities
1724
{
1725
1726
/** \struct hasEntity
1727
\ingroup YaspGrid
1728
*/
1729
1730
/** \struct hasBackupRestoreFacilities
1731
\ingroup YaspGrid
1732
*/
1733
1734
/** \brief YaspGrid has only one geometry type for codim 0 entities
1735
\ingroup YaspGrid
1736
*/
1737
template<int dim>
1738
struct hasSingleGeometryType< YaspGrid<dim> >
1739
{
1740
static const bool v = true;
1741
static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
1742
};
1743
1744
/** \brief YaspGrid is a Cartesian grid
1745
\ingroup YaspGrid
1746
*/
1747
template<int dim>
1748
struct isCartesian< YaspGrid<dim> >
1749
{
1750
static const bool v = true;
1751
};
1752
1753
/** \brief YaspGrid has codim=0 entities (elements)
1754
\ingroup YaspGrid
1755
*/
1756
template<int dim>
1757
struct hasEntity< YaspGrid<dim>, 0 >
1758
{
1759
static const bool v = true;
1760
};
1761
1762
/** \brief YaspGrid has codim=dim entities (vertices)
1763
\ingroup YaspGrid
1764
*/
1765
template<int dim>
1766
struct hasEntity< YaspGrid<dim>, dim >
1767
{
1768
static const bool v = true;
1769
};
1770
1771
template< int dim >
1772
struct canCommunicate< YaspGrid< dim >, 0 >
1773
{
1774
static const bool v = true;
1775
};
1776
1777
template< int dim >
1778
struct canCommunicate< YaspGrid< dim >, dim >
1779
{
1780
static const bool v = true;
1781
};
1782
1783
/** \brief YaspGrid is parallel
1784
\ingroup YaspGrid
1785
*/
1786
template<int dim>
1787
struct isParallel< YaspGrid<dim> >
1788
{
1789
static const bool v = true;
1790
};
1791
1792
/** \brief YaspGrid is levelwise conforming
1793
\ingroup YaspGrid
1794
*/
1795
template<int dim>
1796
struct isLevelwiseConforming< YaspGrid<dim> >
1797
{
1798
static const bool v = true;
1799
};
1800
1801
/** \brief YaspGrid is leafwise conforming
1802
\ingroup YaspGrid
1803
*/
1804
template<int dim>
1805
struct isLeafwiseConforming< YaspGrid<dim> >
1806
{
1807
static const bool v = true;
1808
};
1809
1810
}
3447
1811
3448
1812
} // end namespace
3449
1813
3450
1814
3451
1815
#endif
3452
Loggerhead is a web-based interface for Breezy
Version: 2.0.1