← Back to branch summary

~ubuntu-branches/ubuntu/vivid/dune-grid/vivid

« back to all changes in this revision

Viewing changes to dune/grid/test/test-ug.cc

  • 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.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dune/grid/test/test-ug.cc
1
 
// $Id: test-ug.cc 7990 2012-04-10 15:34:00Z christi $
 
1
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 
2
// vi: set et ts=4 sw=2 sts=2:
 
3
// $Id$
2
4
 
3
5
#include <config.h>
4
6
 
6
8
 
7
9
/*
8
10
 
9
 
  Instantiate UG-Grid and feed it to the generic gridcheck()
 
11
   Instantiate UG-Grid and feed it to the generic gridcheck()
10
12
 
11
 
*/
 
13
 */
12
14
 
13
15
#include <dune/grid/uggrid.hh>
14
16
#include <doc/grids/gridfactory/hybridtestgrids.hh>
18
20
#include "checkgeometryinfather.cc"
19
21
#include "checkintersectionit.cc"
20
22
 
21
 
#include <dune/common/mpihelper.hh>
 
23
#include <dune/common/parallel/mpihelper.hh>
22
24
 
23
25
using namespace Dune;
24
26
 
25
27
class ArcOfCircle : public Dune::BoundarySegment<2>
26
28
{
27
 
    public:
 
29
public:
28
30
 
29
 
    ArcOfCircle(const Dune::FieldVector<double,2>& center, double radius,
30
 
                double fromAngle, double toAngle)
 
31
  ArcOfCircle(const Dune::FieldVector<double,2>& center, double radius,
 
32
              double fromAngle, double toAngle)
31
33
    : center_(center), radius_(radius), fromAngle_(fromAngle), toAngle_(toAngle)
32
 
    {}
33
 
 
34
 
    Dune::FieldVector<double,2> operator()(const Dune::FieldVector<double,1>& local) const {
35
 
 
36
 
        double angle = fromAngle_ + local[0]*(toAngle_ - fromAngle_);
37
 
 
38
 
        Dune::FieldVector<double,2> result = center_;
39
 
        result[0] += radius_ * std::cos(angle);
40
 
        result[1] += radius_ * std::sin(angle);
41
 
 
42
 
        return result;
43
 
    }
44
 
 
45
 
    Dune::FieldVector<double,2> center_;
46
 
 
47
 
    double radius_;
48
 
 
49
 
    double fromAngle_;
50
 
 
51
 
    double toAngle_;
 
34
  {}
 
35
 
 
36
  Dune::FieldVector<double,2> operator()(const Dune::FieldVector<double,1>& local) const {
 
37
 
 
38
    double angle = fromAngle_ + local[0]*(toAngle_ - fromAngle_);
 
39
 
 
40
    Dune::FieldVector<double,2> result = center_;
 
41
    result[0] += radius_ * std::cos(angle);
 
42
    result[1] += radius_ * std::sin(angle);
 
43
 
 
44
    return result;
 
45
  }
 
46
 
 
47
  Dune::FieldVector<double,2> center_;
 
48
 
 
49
  double radius_;
 
50
 
 
51
  double fromAngle_;
 
52
 
 
53
  double toAngle_;
52
54
};
53
55
 
54
56
 
55
57
void makeHalfCircleQuad(Dune::UGGrid<2>& grid, bool boundarySegments, bool parametrization)
56
58
{
57
 
    Dune::GridFactory<Dune::UGGrid<2> > factory(&grid);
58
 
 
59
 
    // /////////////////////////////
60
 
    //   Create boundary segments
61
 
    // /////////////////////////////
62
 
    if (boundarySegments) {
63
 
 
64
 
        FieldVector<double,2> center(0);
65
 
        center[1] = 15;
66
 
        
67
 
        std::vector<unsigned int> vertices(2);
68
 
            
69
 
        if (parametrization) {
70
 
            
71
 
            vertices[0] = 1;  vertices[1] = 2;
72
 
            factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI, M_PI*4/3)));
73
 
            
74
 
            vertices[0] = 2;  vertices[1] = 3;
75
 
            factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI*4/3, M_PI*5/3)));
76
 
            
77
 
            vertices[0] = 3;  vertices[1] = 0;
78
 
            factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI*5/3, M_PI*2)));
79
 
            
80
 
        } else {
81
 
 
82
 
            vertices[0] = 1;  vertices[1] = 2;
83
 
            factory.insertBoundarySegment(vertices);
84
 
            
85
 
            vertices[0] = 2;  vertices[1] = 3;
86
 
            factory.insertBoundarySegment(vertices);
87
 
            
88
 
            vertices[0] = 3;  vertices[1] = 0;
89
 
            factory.insertBoundarySegment(vertices);
90
 
            
91
 
        }
 
59
  Dune::GridFactory<Dune::UGGrid<2> > factory(&grid);
 
60
 
 
61
  // /////////////////////////////
 
62
  //   Create boundary segments
 
63
  // /////////////////////////////
 
64
  if (boundarySegments) {
 
65
 
 
66
    FieldVector<double,2> center(0);
 
67
    center[1] = 15;
 
68
 
 
69
    std::vector<unsigned int> vertices(2);
 
70
 
 
71
    if (parametrization) {
 
72
 
 
73
      vertices[0] = 1;  vertices[1] = 2;
 
74
      factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI, M_PI*4/3)));
 
75
 
 
76
      vertices[0] = 2;  vertices[1] = 3;
 
77
      factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI*4/3, M_PI*5/3)));
 
78
 
 
79
      vertices[0] = 3;  vertices[1] = 0;
 
80
      factory.insertBoundarySegment(vertices, shared_ptr<BoundarySegment<2,2> >(new ArcOfCircle(center, 15, M_PI*5/3, M_PI*2)));
 
81
 
 
82
    } else {
 
83
 
 
84
      vertices[0] = 1;  vertices[1] = 2;
 
85
      factory.insertBoundarySegment(vertices);
 
86
 
 
87
      vertices[0] = 2;  vertices[1] = 3;
 
88
      factory.insertBoundarySegment(vertices);
 
89
 
 
90
      vertices[0] = 3;  vertices[1] = 0;
 
91
      factory.insertBoundarySegment(vertices);
92
92
 
93
93
    }
94
94
 
95
 
    // ///////////////////////
96
 
    //   Insert vertices
97
 
    // ///////////////////////
98
 
    FieldVector<double,2> pos;
99
 
 
100
 
    pos[0] = 15;  pos[1] = 15;
101
 
    factory.insertVertex(pos);
102
 
 
103
 
    pos[0] = -15; pos[1] = 15;
104
 
    factory.insertVertex(pos);
105
 
 
106
 
    pos[0] = -7.5; pos[1] = 2.00962;
107
 
    factory.insertVertex(pos);
108
 
 
109
 
    pos[0] = 7.5; pos[1] = 2.00962;
110
 
    factory.insertVertex(pos);
111
 
 
112
 
    // /////////////////
113
 
    // Insert elements
114
 
    // /////////////////
115
 
 
116
 
    std::vector<unsigned int> cornerIDs(4);
117
 
    cornerIDs[0] = 0;
118
 
    cornerIDs[1] = 1;
119
 
    cornerIDs[2] = 3;
120
 
    cornerIDs[3] = 2;
121
 
 
122
 
    factory.insertElement(GeometryType(GeometryType::cube,2), cornerIDs);
123
 
 
124
 
    // //////////////////////////////////////
125
 
    //   Finish initialization
126
 
    // //////////////////////////////////////
127
 
    factory.createGrid();
 
95
  }
 
96
 
 
97
  // ///////////////////////
 
98
  //   Insert vertices
 
99
  // ///////////////////////
 
100
  FieldVector<double,2> pos;
 
101
 
 
102
  pos[0] = 15;  pos[1] = 15;
 
103
  factory.insertVertex(pos);
 
104
 
 
105
  pos[0] = -15; pos[1] = 15;
 
106
  factory.insertVertex(pos);
 
107
 
 
108
  pos[0] = -7.5; pos[1] = 2.00962;
 
109
  factory.insertVertex(pos);
 
110
 
 
111
  pos[0] = 7.5; pos[1] = 2.00962;
 
112
  factory.insertVertex(pos);
 
113
 
 
114
  // /////////////////
 
115
  // Insert elements
 
116
  // /////////////////
 
117
 
 
118
  std::vector<unsigned int> cornerIDs(4);
 
119
  cornerIDs[0] = 0;
 
120
  cornerIDs[1] = 1;
 
121
  cornerIDs[2] = 3;
 
122
  cornerIDs[3] = 2;
 
123
 
 
124
  factory.insertElement(GeometryType(GeometryType::cube,2), cornerIDs);
 
125
 
 
126
  // //////////////////////////////////////
 
127
  //   Finish initialization
 
128
  // //////////////////////////////////////
 
129
  factory.createGrid();
128
130
 
129
131
}
130
132
 
135
137
  typedef typename GridType::template Codim<0>::LeafIterator LeafIterator;
136
138
 
137
139
  int count = 0;
138
 
  
 
140
 
139
141
  LeafIterator endit = grid.template leafend  <0> ();
140
142
  for(LeafIterator it = grid.template leafbegin<0> (); it != endit ; ++it )
141
143
  {
150
152
 
151
153
void generalTests(bool greenClosure)
152
154
{
153
 
    // /////////////////////////////////////////////////////////////////
154
 
    //   Prelude: the UGGrid implementation relies on the face that std::array,
155
 
    //   Dune::FieldVector, and C-arrays have the same memory layout.  Here
156
 
    //   we make sure this is really so.
157
 
    // /////////////////////////////////////////////////////////////////
158
 
 
159
 
    double cArray[3] = {1,2,3};
160
 
 
161
 
    for (int i=0; i<3; i++) {
162
 
        assert(cArray[i] == (*((FieldVector<double,3>*)&cArray))[i]);
163
 
        assert(cArray[i] == (*((array<double,3>*)&cArray))[i]);
164
 
    }
165
 
 
166
 
    // //////////////////////////////////////////////////////////
167
 
    //   Make some grids for testing
168
 
    // //////////////////////////////////////////////////////////
169
 
 
170
 
    std::auto_ptr<Dune::UGGrid<2> > grid2d(make2DHybridTestGrid<Dune::UGGrid<2> >());
171
 
    std::auto_ptr<Dune::UGGrid<3> > grid3d(make3DHybridTestGrid<Dune::UGGrid<3> >());
172
 
    
173
 
    // Switch of the green closure, if requested
174
 
    if (!greenClosure) {
175
 
        grid2d->setClosureType(UGGrid<2>::NONE);
176
 
        grid3d->setClosureType(UGGrid<3>::NONE);
177
 
    }
178
 
 
179
 
    // check macro grid 
180
 
    gridcheck(*grid2d);
181
 
    gridcheck(*grid3d);
182
 
    
183
 
    // check communication interface 
184
 
    checkCommunication(*grid2d,-1,Dune::dvverb);
185
 
    checkCommunication(*grid3d,-1,Dune::dvverb);
186
 
    for(int l=0; l<=grid2d->maxLevel(); ++l)
187
 
        checkCommunication(*grid2d,l,Dune::dvverb);
188
 
    for(int l=0; l<=grid3d->maxLevel(); ++l)
189
 
        checkCommunication(*grid3d,l,Dune::dvverb);
190
 
 
191
 
    // create hybrid grid 
192
 
    markOne(*grid2d,0,1) ;
193
 
    markOne(*grid3d,0,1) ;
194
 
    
195
 
    gridcheck(*grid2d);
196
 
    gridcheck(*grid3d);
197
 
 
198
 
    // check communication interface 
199
 
    checkCommunication(*grid2d,-1,Dune::dvverb);
200
 
    checkCommunication(*grid3d,-1,Dune::dvverb);
201
 
    for(int l=0; l<=grid2d->maxLevel(); ++l)
202
 
        checkCommunication(*grid2d,l,Dune::dvverb);
203
 
    for(int l=0; l<=grid3d->maxLevel(); ++l)
204
 
        checkCommunication(*grid3d,l,Dune::dvverb);
205
 
    
206
 
    grid2d->globalRefine(1);
207
 
    grid3d->globalRefine(1);
208
 
 
209
 
    gridcheck(*grid2d);
210
 
    gridcheck(*grid3d);
211
 
    // check communication interface 
212
 
    checkCommunication(*grid2d,-1,Dune::dvverb);
213
 
    checkCommunication(*grid3d,-1,Dune::dvverb);
214
 
    for(int l=0; l<=grid2d->maxLevel(); ++l)
215
 
        checkCommunication(*grid2d,l,Dune::dvverb);
216
 
    for(int l=0; l<=grid3d->maxLevel(); ++l)
217
 
        checkCommunication(*grid3d,l,Dune::dvverb);
218
 
   
219
 
    // check geometry lifetime
220
 
    checkGeometryLifetime( grid2d->leafView() );
221
 
    checkGeometryLifetime( grid3d->leafView() );
222
 
  
223
 
    // check the method geometryInFather()
224
 
    checkGeometryInFather(*grid2d);
225
 
    checkGeometryInFather(*grid3d);
226
 
    
227
 
    // check the intersection iterator
228
 
    checkIntersectionIterator(*grid2d);
229
 
    checkIntersectionIterator(*grid3d);
230
 
    
 
155
  // /////////////////////////////////////////////////////////////////
 
156
  //   Prelude: the UGGrid implementation relies on the face that std::array,
 
157
  //   Dune::FieldVector, and C-arrays have the same memory layout.  Here
 
158
  //   we make sure this is really so.
 
159
  // /////////////////////////////////////////////////////////////////
 
160
 
 
161
  double cArray[3] = {1,2,3};
 
162
 
 
163
  for (int i=0; i<3; i++) {
 
164
    assert(cArray[i] == (*((FieldVector<double,3>*)&cArray))[i]);
 
165
    assert(cArray[i] == (*((array<double,3>*)&cArray))[i]);
 
166
  }
 
167
 
 
168
  // //////////////////////////////////////////////////////////
 
169
  //   Make some grids for testing
 
170
  // //////////////////////////////////////////////////////////
 
171
 
 
172
  std::auto_ptr<Dune::UGGrid<2> > grid2d(make2DHybridTestGrid<Dune::UGGrid<2> >());
 
173
  std::auto_ptr<Dune::UGGrid<3> > grid3d(make3DHybridTestGrid<Dune::UGGrid<3> >());
 
174
 
 
175
  // Switch of the green closure, if requested
 
176
  if (!greenClosure) {
 
177
    grid2d->setClosureType(UGGrid<2>::NONE);
 
178
    grid3d->setClosureType(UGGrid<3>::NONE);
 
179
  }
 
180
 
 
181
  // check macro grid
 
182
  gridcheck(*grid2d);
 
183
  gridcheck(*grid3d);
 
184
 
 
185
  // check communication interface
 
186
  checkCommunication(*grid2d,-1,Dune::dvverb);
 
187
  checkCommunication(*grid3d,-1,Dune::dvverb);
 
188
  for(int l=0; l<=grid2d->maxLevel(); ++l)
 
189
    checkCommunication(*grid2d,l,Dune::dvverb);
 
190
  for(int l=0; l<=grid3d->maxLevel(); ++l)
 
191
    checkCommunication(*grid3d,l,Dune::dvverb);
 
192
 
 
193
  // create hybrid grid
 
194
  markOne(*grid2d,0,1) ;
 
195
  markOne(*grid3d,0,1) ;
 
196
 
 
197
  gridcheck(*grid2d);
 
198
  gridcheck(*grid3d);
 
199
 
 
200
  // check communication interface
 
201
  checkCommunication(*grid2d,-1,Dune::dvverb);
 
202
  checkCommunication(*grid3d,-1,Dune::dvverb);
 
203
  for(int l=0; l<=grid2d->maxLevel(); ++l)
 
204
    checkCommunication(*grid2d,l,Dune::dvverb);
 
205
  for(int l=0; l<=grid3d->maxLevel(); ++l)
 
206
    checkCommunication(*grid3d,l,Dune::dvverb);
 
207
 
 
208
  grid2d->globalRefine(1);
 
209
  grid3d->globalRefine(1);
 
210
 
 
211
  gridcheck(*grid2d);
 
212
  gridcheck(*grid3d);
 
213
  // check communication interface
 
214
  checkCommunication(*grid2d,-1,Dune::dvverb);
 
215
  checkCommunication(*grid3d,-1,Dune::dvverb);
 
216
  for(int l=0; l<=grid2d->maxLevel(); ++l)
 
217
    checkCommunication(*grid2d,l,Dune::dvverb);
 
218
  for(int l=0; l<=grid3d->maxLevel(); ++l)
 
219
    checkCommunication(*grid3d,l,Dune::dvverb);
 
220
 
 
221
  // check geometry lifetime
 
222
  checkGeometryLifetime( grid2d->leafView() );
 
223
  checkGeometryLifetime( grid3d->leafView() );
 
224
 
 
225
  // check the method geometryInFather()
 
226
  checkGeometryInFather(*grid2d);
 
227
  checkGeometryInFather(*grid3d);
 
228
 
 
229
  // check the intersection iterator
 
230
  checkIntersectionIterator(*grid2d);
 
231
  checkIntersectionIterator(*grid3d);
 
232
 
231
233
}
232
234
 
233
235
int main (int argc , char **argv) try
234
236
{
235
 
    // use MPI helper to initialize MPI 
236
 
    MPIHelper :: instance( argc, argv );
237
 
 
238
 
    // ////////////////////////////////////////////////////////////////////////
239
 
    //  Do the standard grid test for a 2d and a 3d UGGrid
240
 
    // ////////////////////////////////////////////////////////////////////////
241
 
 
242
 
    // Do the general tests for red/green refinement
243
 
    std::cout << "Testing UGGrid<2> and UGGrid<3> with red/green refinement" << std::endl;
244
 
    generalTests(true);
245
 
 
246
 
    // Do the general tests for nonconforming refinement
247
 
    std::cout << "Testing UGGrid<2> and UGGrid<3> with nonconforming refinement" << std::endl;
248
 
    generalTests(false);
249
 
 
250
 
    // ////////////////////////////////////////////////////////////////////////////
251
 
    //   Test whether I can create a grid with explict boundary segment ordering,
252
 
    //   but not parametrization functions (only 2d, so far)
253
 
    // ////////////////////////////////////////////////////////////////////////////
254
 
 
255
 
    Dune::UGGrid<2> gridWithOrderedBoundarySegments;
256
 
    makeHalfCircleQuad(gridWithOrderedBoundarySegments, true, false);
257
 
 
258
 
    gridWithOrderedBoundarySegments.globalRefine(1);
259
 
    gridcheck(gridWithOrderedBoundarySegments);
260
 
 
261
 
    // ////////////////////////////////////////////////////////////////////////
262
 
    //   Check whether geometryInFather returns equal results with and
263
 
    //   without parametrized boundaries
264
 
    // ////////////////////////////////////////////////////////////////////////
265
 
 
266
 
    Dune::UGGrid<2> gridWithParametrization, gridWithoutParametrization;
 
237
  // use MPI helper to initialize MPI
 
238
  MPIHelper :: instance( argc, argv );
 
239
 
 
240
  // ////////////////////////////////////////////////////////////////////////
 
241
  //  Do the standard grid test for a 2d and a 3d UGGrid
 
242
  // ////////////////////////////////////////////////////////////////////////
 
243
 
 
244
  // Do the general tests for red/green refinement
 
245
  std::cout << "Testing UGGrid<2> and UGGrid<3> with red/green refinement" << std::endl;
 
246
  generalTests(true);
 
247
 
 
248
  // Do the general tests for nonconforming refinement
 
249
  std::cout << "Testing UGGrid<2> and UGGrid<3> with nonconforming refinement" << std::endl;
 
250
  generalTests(false);
 
251
 
 
252
  // ////////////////////////////////////////////////////////////////////////////
 
253
  //   Test whether I can create a grid with explict boundary segment ordering,
 
254
  //   but not parametrization functions (only 2d, so far)
 
255
  // ////////////////////////////////////////////////////////////////////////////
 
256
#ifdef ModelP
 
257
  {  // Force destruction of test grid right after the next test,
 
258
     // because for parallel UG there can be only one (2d-)grid at once.
 
259
#endif
 
260
  Dune::UGGrid<2> gridWithOrderedBoundarySegments;
 
261
  makeHalfCircleQuad(gridWithOrderedBoundarySegments, true, false);
 
262
 
 
263
  gridWithOrderedBoundarySegments.globalRefine(1);
 
264
  gridcheck(gridWithOrderedBoundarySegments);
 
265
#ifdef ModelP
 
266
  }
 
267
#endif
 
268
  // ////////////////////////////////////////////////////////////////////////
 
269
  //   Check whether geometryInFather returns equal results with and
 
270
  //   without parametrized boundaries
 
271
  // ////////////////////////////////////////////////////////////////////////
 
272
 
 
273
  // Only the sequential UG can provide more than one 2d- or 3d-grid at once.
 
274
  // Therefore we do not perform the following test for parallel UGGrid.
 
275
#if ! defined ModelP
 
276
  Dune::UGGrid<2> gridWithParametrization, gridWithoutParametrization;
 
277
 
 
278
  // make grids
 
279
  makeHalfCircleQuad(gridWithoutParametrization, false, false);
 
280
  makeHalfCircleQuad(gridWithParametrization, true, true);
 
281
 
 
282
  // make grids again just to check this is possible
 
283
  makeHalfCircleQuad(gridWithoutParametrization, false, false);
 
284
  makeHalfCircleQuad(gridWithParametrization, true, true);
 
285
 
 
286
  gridWithParametrization.globalRefine(1);
 
287
  gridWithoutParametrization.globalRefine(1);
 
288
 
 
289
  typedef Dune::UGGrid<2>::Codim<0>::LevelIterator ElementIterator;
 
290
  ElementIterator eIt    = gridWithParametrization.lbegin<0>(1);
 
291
  ElementIterator eWoIt  = gridWithoutParametrization.lbegin<0>(1);
 
292
  ElementIterator eEndIt = gridWithParametrization.lend<0>(1);
 
293
 
 
294
  for (; eIt!=eEndIt; ++eIt, ++eWoIt) {
 
295
 
 
296
    // The grids where constructed identically and they are traversed identically
 
297
    // Thus their respective output from geometryInFather should be the same
 
298
    for (int i=0; i<eIt->geometry().corners(); i++) {
 
299
 
 
300
      Dune::FieldVector<double,2> diff = eIt->geometryInFather().corner(i) - eWoIt->geometryInFather().corner(i);
 
301
 
 
302
      if ( diff.two_norm() > 1e-5 )
 
303
        DUNE_THROW(Dune::GridError, "output of geometryInFather() depends on boundary parametrization!");
 
304
 
 
305
    }
 
306
 
 
307
  }
 
308
#endif
 
309
  // ////////////////////////////////////////////////////////////////////////
 
310
  //   Check whether copies of elements have the same global ID
 
311
  // ////////////////////////////////////////////////////////////////////////
 
312
 
 
313
  {
 
314
    std::cout << "Testing if copies of elements have the same id." << std::endl;
 
315
    Dune::UGGrid<2> locallyRefinedGrid;
 
316
 
 
317
    locallyRefinedGrid.setRefinementType(Dune::UGGrid<2>::COPY);
 
318
 
 
319
    typedef Dune::UGGrid<2>::Codim<0>::LevelIterator ElementIterator;
 
320
    typedef Dune::UGGrid<2>::HierarchicIterator HierarchicIterator;
 
321
    typedef Dune::UGGrid<2>::Traits::GlobalIdSet GlobalIdSet;
 
322
    typedef Dune::UGGrid<2>::Traits::LocalIdSet LocalIdSet;
267
323
 
268
324
    // make grids
269
 
    makeHalfCircleQuad(gridWithoutParametrization, false, false);
270
 
    makeHalfCircleQuad(gridWithParametrization, true, true);
271
 
 
272
 
    // make grids again just to check this is possible
273
 
    makeHalfCircleQuad(gridWithoutParametrization, false, false);
274
 
    makeHalfCircleQuad(gridWithParametrization, true, true);
275
 
 
276
 
    gridWithParametrization.globalRefine(1);
277
 
    gridWithoutParametrization.globalRefine(1);
278
 
 
279
 
    typedef Dune::UGGrid<2>::Codim<0>::LevelIterator ElementIterator;
280
 
    ElementIterator eIt    = gridWithParametrization.lbegin<0>(1);
281
 
    ElementIterator eWoIt  = gridWithoutParametrization.lbegin<0>(1);
282
 
    ElementIterator eEndIt = gridWithParametrization.lend<0>(1);
283
 
 
284
 
    for (; eIt!=eEndIt; ++eIt, ++eWoIt) {
285
 
 
286
 
        // The grids where constructed identically and they are traversed identically
287
 
        // Thus their respective output from geometryInFather should be the same
288
 
        for (int i=0; i<eIt->geometry().corners(); i++) {
289
 
 
290
 
            Dune::FieldVector<double,2> diff = eIt->geometryInFather().corner(i) - eWoIt->geometryInFather().corner(i);
291
 
            
292
 
            if ( diff.two_norm() > 1e-5 )
293
 
                DUNE_THROW(Dune::GridError, "output of geometryInFather() depends on boundary parametrization!");
294
 
 
295
 
        }
296
 
 
297
 
    }
298
 
 
299
 
    // ////////////////////////////////////////////////////////////////////////
300
 
    //   Check whether copies of elements have the same global ID
301
 
    // ////////////////////////////////////////////////////////////////////////
302
 
 
 
325
    makeHalfCircleQuad(locallyRefinedGrid, false, false);
 
326
 
 
327
    markOne(locallyRefinedGrid,0,1);
 
328
    markOne(locallyRefinedGrid,0,1);
 
329
 
 
330
    const GlobalIdSet& globalIdSet = locallyRefinedGrid.globalIdSet();
 
331
    const LocalIdSet&  localIdSet  = locallyRefinedGrid.localIdSet();
 
332
 
 
333
    for (int level=0; level<locallyRefinedGrid.maxLevel(); ++level)
303
334
    {
304
 
        std::cout << "Testing if copies of elements have the same id." << std::endl;
305
 
        Dune::UGGrid<2> locallyRefinedGrid;
306
 
        
307
 
        locallyRefinedGrid.setRefinementType(Dune::UGGrid<2>::COPY);
308
 
        
309
 
        typedef Dune::UGGrid<2>::Codim<0>::LevelIterator ElementIterator;
310
 
        typedef Dune::UGGrid<2>::HierarchicIterator HierarchicIterator;
311
 
        typedef Dune::UGGrid<2>::Traits::GlobalIdSet GlobalIdSet;
312
 
        typedef Dune::UGGrid<2>::Traits::LocalIdSet LocalIdSet;
313
 
        
314
 
        // make grids
315
 
        makeHalfCircleQuad(locallyRefinedGrid, false, false);
316
 
        
317
 
        markOne(locallyRefinedGrid,0,1);
318
 
        markOne(locallyRefinedGrid,0,1);
319
 
        
320
 
        const GlobalIdSet& globalIdSet = locallyRefinedGrid.globalIdSet();
321
 
        const LocalIdSet&  localIdSet  = locallyRefinedGrid.localIdSet();
322
 
        
323
 
        for (int level=0; level<locallyRefinedGrid.maxLevel(); ++level)
 
335
      ElementIterator eIt = locallyRefinedGrid.lbegin<0>(level);
 
336
      ElementIterator eEnd = locallyRefinedGrid.lend<0>(level);
 
337
      for(; eIt!=eEnd; ++eIt)
 
338
      {
 
339
        int children = 0;
 
340
        GlobalIdSet::IdType globalChildId;
 
341
        LocalIdSet::IdType localChildId;
 
342
 
 
343
        HierarchicIterator hIt = eIt->hbegin(level+1);
 
344
        HierarchicIterator hEnd = eIt->hend(level+1);
 
345
 
 
346
 
 
347
        for( ; hIt!=hEnd; ++hIt)
324
348
        {
325
 
            ElementIterator eIt = locallyRefinedGrid.lbegin<0>(level);
326
 
            ElementIterator eEnd = locallyRefinedGrid.lend<0>(level);
327
 
            for(; eIt!=eEnd; ++eIt)
328
 
            {
329
 
                int children = 0;
330
 
                GlobalIdSet::IdType globalChildId;
331
 
                LocalIdSet::IdType  localChildId;
332
 
                
333
 
                HierarchicIterator hIt = eIt->hbegin(level+1);
334
 
                HierarchicIterator hEnd = eIt->hend(level+1);
335
 
                
336
 
                
337
 
                for( ;hIt!=hEnd; ++hIt)
338
 
                {
339
 
                    globalChildId = globalIdSet.id<0>(*hIt);
340
 
                    localChildId =  localIdSet.id<0>(*hIt);
341
 
                    ++children;
342
 
                }
343
 
                if (children != 1) 
344
 
                    continue;
345
 
 
346
 
                if (globalIdSet.id<0>(*eIt) != globalChildId)
347
 
                    DUNE_THROW(Dune::GridError, "Copy of element has different globalId!");
348
 
 
349
 
                if (localIdSet.id<0>(*eIt) != localChildId)
350
 
                    DUNE_THROW(Dune::GridError, "Copy of element has different localId!");
351
 
            }
 
349
          globalChildId = globalIdSet.id<0>(*hIt);
 
350
          localChildId =  localIdSet.id<0>(*hIt);
 
351
          ++children;
352
352
        }
 
353
        if (children != 1)
 
354
          continue;
 
355
 
 
356
        if (globalIdSet.id<0>(*eIt) != globalChildId)
 
357
          DUNE_THROW(Dune::GridError, "Copy of element has different globalId!");
 
358
 
 
359
        if (localIdSet.id<0>(*eIt) != localChildId)
 
360
          DUNE_THROW(Dune::GridError, "Copy of element has different localId!");
 
361
      }
353
362
    }
 
363
  }
354
364
 
355
365
  return 0;
356
 
} catch (Dune::Exception& e) {
357
 
    std::cerr << e << std::endl;
358
 
    return 1;
359
 
 } catch (...) {
360
 
    std::cerr << "Generic exception!" << std::endl;
361
 
    return 2;
362
 
 }
 
366
}
 
367
catch (Dune::Exception& e) {
 
368
  std::cerr << e << std::endl;
 
369
  return 1;
 
370
} catch (...) {
 
371
  std::cerr << "Generic exception!" << std::endl;
 
372
  return 2;
 
373
}