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