DUNE-FEM (unstable)

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