Dune Core Modules (2.4.1)

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