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