← Back to branch summary

~ubuntu-branches/ubuntu/trusty/dune-grid/trusty

« back to all changes in this revision

Viewing changes to dune/grid/albertagrid/misc.hh

  • Committer: Package Import Robot
  • Author(s): Ansgar Burchardt
  • Date: 2012-04-06 11:31:20 UTC
  • Revision ID: package-import@ubuntu.com-20120406113120-x0z4e0qqgfhmaj2a
Tags: upstream-2.2~svn7982
Import upstream version 2.2~svn7982

Show diffs side-by-side

added added

removed removed

Lines of Context:
dune/grid/albertagrid/misc.hh
 
1
#ifndef DUNE_ALBERTA_MISC_HH
 
2
#define DUNE_ALBERTA_MISC_HH
 
3
 
 
4
#include <cassert>
 
5
 
 
6
#include <dune/common/exceptions.hh>
 
7
#include <dune/common/typetraits.hh>
 
8
#include <dune/common/forloop.hh>
 
9
 
 
10
#include <dune/grid/albertagrid/albertaheader.hh>
 
11
 
 
12
#if HAVE_ALBERTA
 
13
 
 
14
// should the coordinates be cached in a vector (required for ALBERTA 2.0)?
 
15
#ifndef DUNE_ALBERTA_CACHE_COORDINATES
 
16
#define DUNE_ALBERTA_CACHE_COORDINATES 1
 
17
#endif
 
18
 
 
19
// set to 1 to use generic geometries in AlbertaGrid
 
20
#ifndef DUNE_ALBERTA_USE_GENERICGEOMETRY
 
21
#define DUNE_ALBERTA_USE_GENERICGEOMETRY 0
 
22
#endif
 
23
 
 
24
namespace Dune
 
25
{
 
26
 
 
27
  // Exceptions
 
28
  // ----------
 
29
 
 
30
  class AlbertaError
 
31
  : public Exception
 
32
  {};
 
33
 
 
34
  class AlbertaIOError
 
35
  : public IOError
 
36
  {};
 
37
 
 
38
 
 
39
 
 
40
  namespace Alberta
 
41
  {
 
42
 
 
43
    // Import Types
 
44
    // ------------
 
45
 
 
46
    static const int dimWorld = DIM_OF_WORLD;
 
47
 
 
48
    typedef ALBERTA REAL Real;
 
49
    typedef ALBERTA REAL_B LocalVector; // in barycentric coordinates
 
50
    typedef ALBERTA REAL_D GlobalVector;
 
51
    typedef ALBERTA REAL_DD GlobalMatrix;
 
52
 
 
53
#if DUNE_ALBERTA_VERSION >= 0x300
 
54
    typedef ALBERTA AFF_TRAFO AffineTransformation;
 
55
#else
 
56
    struct AffineTransformation
 
57
    {
 
58
      GlobalMatrix M;
 
59
      GlobalVector t;
 
60
    };
 
61
#endif
 
62
 
 
63
    typedef ALBERTA MESH Mesh;
 
64
    typedef ALBERTA EL Element;
 
65
 
 
66
    static const int meshRefined = MESH_REFINED;
 
67
    static const int meshCoarsened = MESH_COARSENED;
 
68
 
 
69
    static const int InteriorBoundary = INTERIOR;
 
70
    static const int DirichletBoundary = DIRICHLET;
 
71
#if DUNE_ALBERTA_VERSION >= 0x300
 
72
    typedef ALBERTA BNDRY_TYPE BoundaryId;
 
73
#else
 
74
    typedef S_CHAR BoundaryId;
 
75
#endif
 
76
 
 
77
    typedef U_CHAR ElementType;
 
78
 
 
79
    typedef ALBERTA FE_SPACE DofSpace;
 
80
 
 
81
 
 
82
 
 
83
    // Memory Manipulation Functions
 
84
    // -----------------------------
 
85
 
 
86
    template< class Data >
 
87
    inline Data *memAlloc ( size_t size )
 
88
    {
 
89
      return MEM_ALLOC( size, Data );
 
90
    }
 
91
 
 
92
    template< class Data >
 
93
    inline Data *memCAlloc ( size_t size )
 
94
    {
 
95
      return MEM_CALLOC( size, Data );
 
96
    }
 
97
 
 
98
    template< class Data >
 
99
    inline Data *memReAlloc ( Data *ptr, size_t oldSize, size_t newSize )
 
100
    {
 
101
      return MEM_REALLOC( ptr, oldSize, newSize, Data );
 
102
    }
 
103
 
 
104
    template< class Data >
 
105
    inline void memFree ( Data *ptr, size_t size )
 
106
    {
 
107
      return MEM_FREE( ptr, size, Data );
 
108
    }
 
109
 
 
110
 
 
111
 
 
112
    // GlobalSpace
 
113
    // -----------
 
114
 
 
115
    class GlobalSpace
 
116
    {
 
117
      typedef GlobalSpace This;
 
118
 
 
119
    public:
 
120
      typedef GlobalMatrix Matrix;
 
121
      typedef GlobalVector Vector;
 
122
 
 
123
    private:
 
124
      Matrix identityMatrix_;
 
125
      Vector nullVector_;
 
126
 
 
127
      GlobalSpace ()
 
128
      {
 
129
        for( int i = 0; i < dimWorld; ++i )
 
130
        {
 
131
          for( int j = 0; j < dimWorld; ++j )
 
132
            identityMatrix_[ i ][ j ] = Real( 0 );
 
133
          identityMatrix_[ i ][ i ] = Real( 1 );
 
134
          nullVector_[ i ] = Real( 0 );
 
135
        }
 
136
      }
 
137
 
 
138
      static This &instance ()
 
139
      {
 
140
        static This theInstance;
 
141
        return theInstance;
 
142
      }
 
143
 
 
144
    public:
 
145
      static const Matrix &identityMatrix ()
 
146
      {
 
147
        return instance().identityMatrix_;
 
148
      }
 
149
 
 
150
      static const Vector &nullVector ()
 
151
      {
 
152
        return instance().nullVector_;
 
153
      }
 
154
    };
 
155
 
 
156
 
 
157
 
 
158
    // NumSubEntities
 
159
    // --------------
 
160
 
 
161
    template< int dim, int codim >
 
162
    struct NumSubEntities;
 
163
 
 
164
    template< int dim >
 
165
    struct NumSubEntities< dim, 0 >
 
166
    {
 
167
      static const int value = 1;
 
168
    };
 
169
 
 
170
    template< int dim >
 
171
    struct NumSubEntities< dim, dim >
 
172
    {
 
173
      static const int value = dim+1;
 
174
    };
 
175
 
 
176
    template<>
 
177
    struct NumSubEntities< 0, 0 >
 
178
    {
 
179
      static const int value = 1;
 
180
    };
 
181
 
 
182
    template<>
 
183
    struct NumSubEntities< 2, 1 >
 
184
    {
 
185
      static const int value = 3;
 
186
    };
 
187
 
 
188
    template<>
 
189
    struct NumSubEntities< 3, 1 >
 
190
    {
 
191
      static const int value = 4;
 
192
    };
 
193
 
 
194
    template<>
 
195
    struct NumSubEntities< 3, 2 >
 
196
    {
 
197
      static const int value = 6;
 
198
    };
 
199
 
 
200
 
 
201
 
 
202
    // CodimType
 
203
    // ---------
 
204
 
 
205
    template< int dim, int codim >
 
206
    struct CodimType;
 
207
 
 
208
    template< int dim >
 
209
    struct CodimType< dim, 0 >
 
210
    {
 
211
      static const int value = CENTER;
 
212
    };
 
213
 
 
214
    template< int dim >
 
215
    struct CodimType< dim, dim >
 
216
    {
 
217
      static const int value = VERTEX;
 
218
    };
 
219
 
 
220
    template<>
 
221
    struct CodimType< 2, 1 >
 
222
    {
 
223
      static const int value = EDGE;
 
224
    };
 
225
 
 
226
    template<>
 
227
    struct CodimType< 3, 1 >
 
228
    {
 
229
      static const int value = FACE;
 
230
    };
 
231
 
 
232
    template<>
 
233
    struct CodimType< 3, 2 >
 
234
    {
 
235
      static const int value = EDGE;
 
236
    };
 
237
 
 
238
 
 
239
 
 
240
    // FillFlags
 
241
    // ---------
 
242
 
 
243
    template< int dim >
 
244
    struct FillFlags
 
245
    {
 
246
      typedef ALBERTA FLAGS Flags;
 
247
 
 
248
      static const Flags nothing = FILL_NOTHING;
 
249
 
 
250
      static const Flags coords = FILL_COORDS;
 
251
 
 
252
      static const Flags neighbor = FILL_NEIGH;
 
253
 
 
254
      static const Flags orientation = (dim == 3 ? FILL_ORIENTATION : FILL_NOTHING);
 
255
 
 
256
      static const Flags projection = FILL_PROJECTION;
 
257
 
 
258
#if DUNE_ALBERTA_VERSION >= 0x300
 
259
      static const Flags elementType = FILL_NOTHING;
 
260
#else
 
261
      static const Flags elementType = (dim == 3 ? FILL_EL_TYPE : FILL_NOTHING);
 
262
#endif
 
263
 
 
264
#if DUNE_ALBERTA_VERSION >= 0x300
 
265
      static const Flags boundaryId = FILL_MACRO_WALLS;
 
266
#else
 
267
      static const Flags boundaryId = FILL_BOUND;
 
268
#endif
 
269
 
 
270
#if DUNE_ALBERTA_VERSION >= 0x300
 
271
      static const Flags nonPeriodic = FILL_NON_PERIODIC;
 
272
#else
 
273
      static const Flags nonPeriodic = FILL_NOTHING;
 
274
#endif
 
275
 
 
276
      static const Flags all = coords | neighbor | boundaryId | nonPeriodic
 
277
                               | orientation | projection | elementType;
 
278
 
 
279
#if DUNE_ALBERTA_VERSION >= 0x300
 
280
      static const Flags standardWithCoords = all & ~nonPeriodic & ~projection;
 
281
#else
 
282
      static const Flags standardWithCoords = all;
 
283
#endif
 
284
 
 
285
#if DUNE_ALBERTA_CACHE_COORDINATES
 
286
      static const Flags standard = standardWithCoords & ~coords;
 
287
#else
 
288
      static const Flags standard = standardWithCoords;
 
289
#endif
 
290
    };
 
291
 
 
292
 
 
293
 
 
294
    // RefinementEdge
 
295
    // --------------
 
296
 
 
297
    template< int dim >
 
298
    struct RefinementEdge
 
299
    {
 
300
      static const int value = 0;
 
301
    };
 
302
 
 
303
    template<>
 
304
    struct RefinementEdge< 2 >
 
305
    {
 
306
      static const int value = 2;
 
307
    };
 
308
 
 
309
 
 
310
 
 
311
    // Dune2AlbertaNumbering
 
312
    // ---------------------
 
313
 
 
314
    template< int dim, int codim >
 
315
    struct Dune2AlbertaNumbering
 
316
    {
 
317
      static int apply ( const int i )
 
318
      {
 
319
        assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
 
320
        return i;
 
321
      }
 
322
    };
 
323
 
 
324
    template<>
 
325
    struct Dune2AlbertaNumbering< 3, 2 >
 
326
    {
 
327
      static const int numSubEntities = NumSubEntities< 3, 2 >::value;
 
328
 
 
329
      static int apply ( const int i )
 
330
      {
 
331
        assert( (i >= 0) && (i < numSubEntities) );
 
332
        static int dune2alberta[ numSubEntities ] = { 0, 3, 1, 2, 4, 5 };
 
333
        return dune2alberta[ i ];
 
334
      }
 
335
    };
 
336
 
 
337
 
 
338
 
 
339
    // Generic2AlbertaNumbering
 
340
    // ------------------------
 
341
 
 
342
    template< int dim, int codim >
 
343
    struct Generic2AlbertaNumbering
 
344
    {
 
345
      static int apply ( const int i )
 
346
      {
 
347
        assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
 
348
        return i;
 
349
      }
 
350
    };
 
351
 
 
352
    template< int dim >
 
353
    struct Generic2AlbertaNumbering< dim, 1 >
 
354
    {
 
355
      static int apply ( const int i )
 
356
      {
 
357
        assert( (i >= 0) && (i < NumSubEntities< dim, 1 >::value) );
 
358
        return dim - i;
 
359
      }
 
360
    };
 
361
 
 
362
    template<>
 
363
    struct Generic2AlbertaNumbering< 1, 1 >
 
364
    {
 
365
      static int apply ( const int i )
 
366
      {
 
367
        assert( (i >= 0) && (i < NumSubEntities< 1, 1 >::value) );
 
368
        return i;
 
369
      }
 
370
    };
 
371
 
 
372
    template<>
 
373
    struct Generic2AlbertaNumbering< 3, 2 >
 
374
    {
 
375
      static const int numSubEntities = NumSubEntities< 3, 2 >::value;
 
376
 
 
377
      static int apply ( const int i )
 
378
      {
 
379
        assert( (i >= 0) && (i < numSubEntities) );
 
380
        static int generic2alberta[ numSubEntities ] = { 0, 1, 3, 2, 4, 5 };
 
381
        return generic2alberta[ i ];
 
382
      }
 
383
    };
 
384
 
 
385
 
 
386
 
 
387
    // NumberingMap
 
388
    // ------------
 
389
 
 
390
    template< int dim, template< int, int > class Numbering = Generic2AlbertaNumbering >
 
391
    class NumberingMap
 
392
    {
 
393
      typedef NumberingMap< dim, Numbering > This;
 
394
 
 
395
      template< int codim >
 
396
      struct Initialize;
 
397
 
 
398
      int *dune2alberta_[ dim+1 ];
 
399
      int *alberta2dune_[ dim+1 ];
 
400
      int numSubEntities_[ dim+1 ];
 
401
 
 
402
      NumberingMap ( const This & );
 
403
      This &operator= ( const This & );
 
404
 
 
405
    public:
 
406
      NumberingMap ()
 
407
      {
 
408
        ForLoop< Initialize, 0, dim >::apply( *this );
 
409
      }
 
410
 
 
411
      ~NumberingMap ()
 
412
      {
 
413
        for( int codim = 0; codim <= dim; ++codim )
 
414
        {
 
415
          delete[]( dune2alberta_[ codim ] );
 
416
          delete[]( alberta2dune_[ codim ] );
 
417
        }
 
418
      }
 
419
 
 
420
      int dune2alberta ( int codim, int i ) const
 
421
      {
 
422
        assert( (codim >= 0) && (codim <= dim) );
 
423
        assert( (i >= 0) && (i < numSubEntities( codim )) );
 
424
        return dune2alberta_[ codim ][ i ];
 
425
      }
 
426
 
 
427
      int alberta2dune ( int codim, int i ) const
 
428
      {
 
429
        assert( (codim >= 0) && (codim <= dim) );
 
430
        assert( (i >= 0) && (i < numSubEntities( codim )) );
 
431
        return alberta2dune_[ codim ][ i ];
 
432
      }
 
433
 
 
434
      int numSubEntities ( int codim ) const
 
435
      {
 
436
        assert( (codim >= 0) && (codim <= dim) );
 
437
        return numSubEntities_[ codim ];
 
438
      }
 
439
    };
 
440
 
 
441
 
 
442
 
 
443
    // NumberingMap::Initialize
 
444
    // ------------------------
 
445
 
 
446
    template< int dim, template< int, int > class Numbering >
 
447
    template< int codim >
 
448
    struct NumberingMap< dim, Numbering >::Initialize
 
449
    {
 
450
      static const int numSubEntities = NumSubEntities< dim, codim >::value;
 
451
 
 
452
      static void apply ( NumberingMap< dim, Numbering > &map )
 
453
      {
 
454
        map.numSubEntities_[ codim ] = numSubEntities;
 
455
        map.dune2alberta_[ codim ] = new int[ numSubEntities ];
 
456
        map.alberta2dune_[ codim ] = new int[ numSubEntities ];
 
457
 
 
458
        for( int i = 0; i < numSubEntities; ++i )
 
459
        {
 
460
          const int j = Numbering< dim, codim >::apply( i );
 
461
          map.dune2alberta_[ codim ][ i ] = j;
 
462
          map.alberta2dune_[ codim ][ j ] = i;
 
463
        }
 
464
      }
 
465
    };
 
466
 
 
467
 
 
468
 
 
469
    // MapVertices
 
470
    // -----------
 
471
 
 
472
    template< int dim, int codim >
 
473
    struct MapVertices;
 
474
 
 
475
    template< int dim >
 
476
    struct MapVertices< dim, 0 >
 
477
    {
 
478
      static int apply ( int subEntity, int vertex )
 
479
      {
 
480
        assert( subEntity == 0 );
 
481
        assert( (vertex >= 0) && (vertex <= NumSubEntities< dim, dim >::value) );
 
482
        return vertex;
 
483
      }
 
484
    };
 
485
 
 
486
    template<>
 
487
    struct MapVertices< 2, 1 >
 
488
    {
 
489
      static int apply ( int subEntity, int vertex )
 
490
      {
 
491
        assert( (subEntity >= 0) && (subEntity < 3) );
 
492
        assert( (vertex >= 0) && (vertex < 2) );
 
493
        //static const int map[ 3 ][ 2 ] = { {1,2}, {2,0}, {0,1} };
 
494
        static const int map[ 3 ][ 2 ] = { {1,2}, {0,2}, {0,1} };
 
495
        return map[ subEntity ][ vertex ];
 
496
      }
 
497
    };
 
498
 
 
499
    template<>
 
500
    struct MapVertices< 3, 1 >
 
501
    {
 
502
      static int apply ( int subEntity, int vertex )
 
503
      {
 
504
        assert( (subEntity >= 0) && (subEntity < 4) );
 
505
        assert( (vertex >= 0) && (vertex < 3) );
 
506
        //static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,3,2}, {0,1,3}, {0,2,1} };
 
507
        static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,2,3}, {0,1,3}, {0,1,2} };
 
508
        return map[ subEntity ][ vertex ];
 
509
      }
 
510
    };
 
511
 
 
512
    template<>
 
513
    struct MapVertices< 3, 2 >
 
514
    {
 
515
      static int apply ( int subEntity, int vertex )
 
516
      {
 
517
        assert( (subEntity >= 0) && (subEntity < 6) );
 
518
        assert( (vertex >= 0) && (vertex < 2) );
 
519
        static const int map[ 6 ][ 2 ] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
 
520
        return map[ subEntity ][ vertex ];
 
521
      }
 
522
    };
 
523
 
 
524
    template< int dim >
 
525
    struct MapVertices< dim, dim >
 
526
    {
 
527
      static int apply ( int subEntity, int vertex )
 
528
      {
 
529
        assert( (subEntity >= 0) && (subEntity < NumSubEntities< dim, 1 >::value) );
 
530
        assert( vertex == 0 );
 
531
        return subEntity;
 
532
      }
 
533
    };
 
534
 
 
535
 
 
536
 
 
537
    // Twist
 
538
    // -----
 
539
 
 
540
    // ******************************************************************
 
541
    // Meaning of the twist (same as in ALU)
 
542
    // -------------------------------------
 
543
    //
 
544
    // Consider a fixed ordering of the vertices v_1, ... v_n of a face
 
545
    // (here, we assume their indices to be increasing). Denote by k the
 
546
    // local number of a vertex v within the element and by t the twist.
 
547
    // Then, v = v_j, where j is computed by the following formula:
 
548
    //
 
549
    //        / (2n + 1 - k + t) % n, if t < 0
 
550
    //   j = <
 
551
    //        \ (k + t) % n,          if t >= 0
 
552
    //
 
553
    //  Note: We use the order of the 0-th vertex dof to assign the twist.
 
554
    //        This is ok for two reasons:
 
555
    //        - ALBERTA preserves the relative order of the dofs during
 
556
    //          dof compression.
 
557
    //        - ALBERTA enforces the first vertex dof admin to be periodic.
 
558
    // ******************************************************************
 
559
 
 
560
    template< int dim, int subdim >
 
561
    struct Twist
 
562
    {
 
563
      static const int numSubEntities = NumSubEntities< dim, dim-subdim >::value;
 
564
 
 
565
      static const int minTwist = 0;
 
566
      static const int maxTwist = 0;
 
567
 
 
568
      static int twist ( const Element *element, int subEntity )
 
569
      {
 
570
        assert( (subEntity >= 0) && (subEntity < numSubEntities) );
 
571
        return 0;
 
572
      }
 
573
    };
 
574
 
 
575
    template< int dim >
 
576
    struct Twist< dim, 1 >
 
577
    {
 
578
      static const int numSubEntities = NumSubEntities< dim, dim-1 >::value;
 
579
 
 
580
      static const int minTwist = 0;
 
581
      static const int maxTwist = 1;
 
582
 
 
583
      static int twist ( const Element *element, int subEntity )
 
584
      {
 
585
        assert( (subEntity >= 0) && (subEntity < numSubEntities) );
 
586
        const int numVertices = NumSubEntities< 1, 1 >::value;
 
587
        int dof[ numVertices ];
 
588
        for( int i = 0; i < numVertices; ++i )
 
589
        {
 
590
          const int j = MapVertices< dim, dim-1 >::apply( subEntity, i );
 
591
          dof[ i ] = element->dof[ j ][ 0 ];
 
592
        }
 
593
        return (dof[ 0 ] < dof[ 1 ] ? 0 : 1);
 
594
      }
 
595
    };
 
596
 
 
597
 
 
598
    template<>
 
599
    struct Twist< 1, 1 >
 
600
    {
 
601
      static const int minTwist = 0;
 
602
      static const int maxTwist = 0;
 
603
 
 
604
      static int twist ( const Element *element, int subEntity )
 
605
      {
 
606
        assert( subEntity == 0 );
 
607
        return 0;
 
608
      }
 
609
    };
 
610
 
 
611
 
 
612
    template< int dim >
 
613
    struct Twist< dim, 2 >
 
614
    {
 
615
      static const int numSubEntities = NumSubEntities< dim, dim-2 >::value;
 
616
 
 
617
      static const int minTwist = -3;
 
618
      static const int maxTwist = 2;
 
619
 
 
620
      static int twist ( const Element *element, int subEntity )
 
621
      {
 
622
        assert( (subEntity >= 0) && (subEntity < numSubEntities) );
 
623
        const int numVertices = NumSubEntities< 2, 2 >::value;
 
624
        int dof[ numVertices ];
 
625
        for( int i = 0; i < numVertices; ++i )
 
626
        {
 
627
          const int j = MapVertices< dim, dim-2 >::apply( subEntity, i );
 
628
          dof[ i ] = element->dof[ j ][ 0 ];
 
629
        }
 
630
 
 
631
        const int twist[ 8 ] = { -2, 1, 666, -1, 2, 666, -3, 0 };
 
632
        const int k = int( dof[ 0 ] < dof[ 1 ] )
 
633
                      | (int( dof[ 0 ] < dof[ 2 ] ) << 1)
 
634
                      | (int( dof[ 1 ] < dof[ 2 ] ) << 2);
 
635
        assert( twist[ k ] != 666 );
 
636
        return twist[ k ];
 
637
      }
 
638
    };
 
639
 
 
640
 
 
641
    template<>
 
642
    struct Twist< 2, 2 >
 
643
    {
 
644
      static const int minTwist = 0;
 
645
      static const int maxTwist = 0;
 
646
 
 
647
      static int twist ( const Element *element, int subEntity )
 
648
      {
 
649
        assert( subEntity == 0 );
 
650
        return 0;
 
651
      }
 
652
    };
 
653
 
 
654
 
 
655
 
 
656
    template< int dim >
 
657
    inline int applyTwist ( int twist, int i )
 
658
    {
 
659
      const int numCorners = NumSubEntities< dim, dim >::value;
 
660
      return (twist < 0 ? (2*numCorners + 1 - i + twist) : i + twist) % numCorners;
 
661
    }
 
662
 
 
663
    template< int dim >
 
664
    inline int applyInverseTwist ( int twist, int i )
 
665
    {
 
666
      const int numCorners = NumSubEntities< dim, dim >::value;
 
667
      return (twist < 0 ? (2*numCorners + 1 - i + twist) : numCorners + i - twist) % numCorners;
 
668
    }
 
669
 
 
670
  }
 
671
 
 
672
}
 
673
 
 
674
#endif // #if HAVE_ALBERTA
 
675
 
 
676
#endif // #ifndef DUNE_ALBERTA_MISC_HH