1 #ifndef DUNE_ALUGRID_COMMON_TWISTS_HH 2 #define DUNE_ALUGRID_COMMON_TWISTS_HH 3 4 #include <iterator> 5 6 #include <dune/common/fvector.hh> 7 8 #include <dune/geometry/affinegeometry.hh> 9 #include <dune/geometry/referenceelements.hh> 10 11 #include <dune/alugrid/common/capabilities.hh> 12 13 namespace Dune 14 { 15 16 // Internal Forward Declarations 17 // ----------------------------- 18 19 template< int corners, int dim > 20 class ALUTwist; 21 22 template< int corners, int dim > 23 class ALUTwists; 24 25 26 27 // ALUTwistIterator 28 // ---------------- 29 30 template< class Twist > 31 struct ALUTwistIterator 32 : public std::iterator< std::input_iterator_tag, Twist, int > 33 { ALUTwistIteratorDune::ALUTwistIterator34 explicit ALUTwistIterator ( Twist twist ) : twist_( twist ) {} 35 operator *Dune::ALUTwistIterator36 const Twist &operator* () const { return twist_; } operator ->Dune::ALUTwistIterator37 const Twist *operator-> () const { return &twist_; } 38 operator ==Dune::ALUTwistIterator39 bool operator== ( const ALUTwistIterator &other ) const { return (twist_ == other.twist_); } operator !=Dune::ALUTwistIterator40 bool operator!= ( const ALUTwistIterator &other ) const { return (twist_ != other.twist_); } 41 operator ++Dune::ALUTwistIterator42 ALUTwistIterator &operator++ () { ++twist_.aluTwist_; return *this; } operator ++Dune::ALUTwistIterator43 ALUTwistIterator operator++ ( int ) { ALUTwistIterator other( *this ); ++(*this); return other; } 44 45 private: 46 Twist twist_; 47 }; 48 49 50 51 // ALUTwist for dimension 2 52 // ------------------------ 53 54 template< int corners > 55 class ALUTwist< corners, 2 > 56 { 57 typedef ALUTwist< corners, 2 > Twist; 58 59 friend struct ALUTwistIterator< Twist >; 60 61 template< class ctype > 62 struct CoordVector 63 { 64 // type of container for reference elements 65 typedef ReferenceElements< ctype, 2 > ReferenceElementContainerType; 66 67 // type of reference element 68 typedef std::decay_t< decltype( ReferenceElementContainerType::general( std::declval< const Dune::GeometryType & >() ) ) > ReferenceElementType; 69 CoordVectorDune::ALUTwist::CoordVector70 CoordVector ( const Twist &twist ) 71 : twist_( twist ), 72 refElement_( ReferenceElements< ctype, 2 >::general( twist.type() ) ) 73 {} 74 operator []Dune::ALUTwist::CoordVector75 const FieldVector< ctype, 2 > operator[] ( int i ) const { return refElement().position( twist_( i ), 2 ); } 76 refElementDune::ALUTwist::CoordVector77 const ReferenceElementType& refElement () const { return refElement_; } 78 79 private: 80 const Twist &twist_; 81 const ReferenceElementType& refElement_; 82 }; 83 84 public: 85 /** \brief dimension */ 86 static const int dimension = 2; 87 88 /** 89 * \name Construction 90 * \{ 91 */ 92 93 /** \brief default constructor; results in identity twist */ ALUTwist()94 ALUTwist () : aluTwist_( 0 ) {} 95 ALUTwist(GeometryType)96 explicit ALUTwist ( GeometryType ) : aluTwist_( 0 ) {} 97 ALUTwist(int aluTwist)98 explicit ALUTwist ( int aluTwist ) : aluTwist_( aluTwist ) {} 99 ALUTwist(int origin,bool positive)100 ALUTwist ( int origin, bool positive ) 101 : aluTwist_( positive ? origin : (origin + corners - 1) % corners - corners ) 102 {} 103 104 /** \} */ 105 106 /** 107 * \name Copying and Assignment 108 * \{ 109 */ 110 111 /** \brief copy constructor */ ALUTwist(const ALUTwist & other)112 ALUTwist ( const ALUTwist &other ) : aluTwist_( other.aluTwist_ ) {} 113 114 /** \brief assignment operator */ operator =(const ALUTwist & other)115 ALUTwist &operator= ( const ALUTwist &other ) { aluTwist_ = other.aluTwist_; return *this; } 116 117 /** \} */ 118 119 /** 120 * \name Group Operations 121 * \{ 122 */ 123 124 /** \brief composition */ operator *(const ALUTwist & other) const125 ALUTwist operator* ( const ALUTwist &other ) const 126 { 127 return ALUTwist( apply( other.apply( 0 ) ), !(positive() ^ other.positive()) ); 128 } 129 130 /** \brief composition with inverse */ operator /(const ALUTwist & other) const131 ALUTwist operator/ ( const ALUTwist &other ) const { return (*this * other.inverse()); } 132 133 /** \brief composition */ operator *=(const ALUTwist & other)134 ALUTwist &operator*= ( const ALUTwist &other ) { return (*this = (*this) * other); } 135 136 /** \brief composition with inverse */ operator /=(const ALUTwist & other)137 ALUTwist &operator/= ( const ALUTwist &other ) { return (*this = (*this) / other); } 138 139 /** \brief obtain inverse */ inverse() const140 ALUTwist inverse () const { return ALUTwist( positive() ? (corners - aluTwist_) % corners : aluTwist_ ); } 141 142 /** \} */ 143 144 /** 145 * \brief Comparison 146 * \{ 147 */ 148 149 /** \brief check for equality */ operator ==(const ALUTwist & other) const150 bool operator== ( const ALUTwist &other ) const { return (aluTwist_ == other.aluTwist_); } 151 152 /** \brief check for inequality */ operator !=(const ALUTwist & other) const153 bool operator!= ( const ALUTwist &other ) const { return (aluTwist_ != other.aluTwist_); } 154 155 /** \} */ 156 157 /** 158 * \brief Topological Corner Mapping 159 * \{ 160 */ 161 162 /** \brief obtain type of reference element */ type() const163 GeometryType type () const 164 { 165 if( corners == 3 ) 166 return GeometryTypes::simplex(dimension); 167 else 168 return GeometryTypes::cube(dimension); 169 } 170 171 /** 172 * \brief map i-th corner 173 * 174 * \param[in] i corner index 175 * 176 * \returns mapped index 177 */ operator ()(int i) const178 int operator() ( int i ) const { return aluVertex2duneVertex( apply( duneVertex2aluVertex( i ) ) ); } 179 operator ()(int i,int codim) const180 int operator() ( int i, int codim ) const { return alu2dune( apply( dune2alu( i, codim ), codim ), codim ); } 181 182 /** \} */ 183 184 /** 185 * \brief Geometric Equivalent 186 * \{ 187 */ 188 189 /** \brief cast into affine geometry \f$F\f$ */ 190 template< class ctype > operator AffineGeometry<ctype,dimension,dimension>() const191 operator AffineGeometry< ctype, dimension, dimension > () const 192 { 193 const CoordVector< ctype > coordVector( *this ); 194 return AffineGeometry< ctype, dimension, dimension >( coordVector.refElement(), coordVector ); 195 } 196 197 /** \brief equivalent to \f$|DF| > 0\f$ */ positive() const198 bool positive () const { return (aluTwist_ >= 0); } 199 200 /** \} */ 201 202 // non-interface methods 203 operator int() const204 operator int () const { return aluTwist_; } 205 206 // apply twist in ALU numbering apply(int i) const207 int apply ( int i ) const { return ((positive() ? i : 2*corners + 1 - i) + aluTwist_) % corners; } apply(int i,int codim) const208 int apply ( int i, int codim ) const { return (codim == 0 ? i : (codim == 1 ? applyEdge( i ) : apply( i ))); } 209 210 private: applyEdge(int i) const211 int applyEdge ( int i ) const { return ((positive() ? i : 2*corners - i) + aluTwist_) % corners; } 212 213 aluEdge2duneEdge(int i)214 static int aluEdge2duneEdge ( int i ) { return ((corners == 3) ? (3 - i) % 3 : (6 - swap23( i )) % 4); } duneEdge2aluEdge(int i)215 static int duneEdge2aluEdge ( int i ) { return ((corners == 3) ? (3 - i) % 3 : swap23( (6 - i) % 4 )); } 216 aluVertex2duneVertex(int i)217 static int aluVertex2duneVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); } duneVertex2aluVertex(int i)218 static int duneVertex2aluVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); } 219 alu2dune(int i,int codim)220 static int alu2dune ( int i, int codim ) { return (codim == 0 ? i : (codim == 1 ? aluEdge2duneEdge( i ) : aluVertex2duneVertex( i ))); } dune2alu(int i,int codim)221 static int dune2alu ( int i, int codim ) { return (codim == 0 ? i : (codim == 1 ? duneEdge2aluEdge( i ) : aluVertex2duneVertex( i ))); } 222 swap23(int i)223 static int swap23 ( int i ) { return i ^ (i >> 1); } 224 225 int aluTwist_; 226 }; 227 228 229 230 // ALUTwist for dimension 1 231 // ------------------------ 232 233 template<> 234 class ALUTwist< 2, 1 > 235 { 236 typedef ALUTwist< 2, 1 > Twist; 237 238 friend struct ALUTwistIterator< Twist >; 239 240 template< class ctype > 241 struct CoordVector 242 { CoordVectorDune::ALUTwist::CoordVector243 CoordVector ( const Twist &twist ) : twist_( twist ) {} 244 operator []Dune::ALUTwist::CoordVector245 FieldVector< ctype, 1 > operator[] ( int i ) const { return FieldVector< ctype, 1 >( ctype( twist_( i ) ) ); } 246 247 private: 248 const Twist &twist_; 249 }; 250 251 public: 252 /** \brief dimension */ 253 static const int dimension = 1; 254 255 /** 256 * \name Construction 257 * \{ 258 */ 259 260 /** \brief default constructor; results in identity twist */ ALUTwist()261 ALUTwist () : aluTwist_( 0 ) {} 262 ALUTwist(GeometryType)263 explicit ALUTwist ( GeometryType ) : aluTwist_( 0 ) {} 264 ALUTwist(int aluTwist)265 explicit ALUTwist ( int aluTwist ) : aluTwist_( aluTwist ) {} 266 ALUTwist(bool positive)267 explicit ALUTwist ( bool positive ) : aluTwist_( positive ) {} 268 269 /** \} */ 270 271 /** 272 * \name Copying and Assignment 273 * \{ 274 */ 275 276 /** \brief copy constructor */ ALUTwist(const ALUTwist & other)277 ALUTwist ( const ALUTwist &other ) : aluTwist_( other.aluTwist_ ) {} 278 279 /** \brief assignment operator */ operator =(const ALUTwist & other)280 ALUTwist &operator= ( const ALUTwist &other ) { aluTwist_ = other.aluTwist_; return *this; } 281 282 /** \} */ 283 284 /** 285 * \name Group Operations 286 * \{ 287 */ 288 289 /** \brief composition */ operator *(const ALUTwist & other) const290 ALUTwist operator* ( const ALUTwist &other ) const { return ALUTwist( aluTwist_ ^ other.aluTwist_ ); } 291 292 /** \brief composition with inverse */ operator /(const ALUTwist & other) const293 ALUTwist operator/ ( const ALUTwist &other ) const { return (*this * other.inverse()); } 294 295 /** \brief composition */ operator *=(const ALUTwist & other)296 ALUTwist &operator*= ( const ALUTwist &other ) { return (*this = (*this) * other); } 297 298 /** \brief composition with inverse */ operator /=(const ALUTwist & other)299 ALUTwist &operator/= ( const ALUTwist &other ) { return (*this = (*this) / other); } 300 301 /** \brief obtain inverse */ inverse() const302 ALUTwist inverse () const { return *this; } 303 304 /** \} */ 305 306 /** 307 * \brief Comparison 308 * \{ 309 */ 310 311 /** \brief check for equality */ operator ==(const ALUTwist & other) const312 bool operator== ( const ALUTwist &other ) const { return (aluTwist_ == other.aluTwist_); } 313 314 /** \brief check for inequality */ operator !=(const ALUTwist & other) const315 bool operator!= ( const ALUTwist &other ) const { return (aluTwist_ != other.aluTwist_); } 316 317 /** \} */ 318 319 /** 320 * \brief Topological Corner Mapping 321 * \{ 322 */ 323 324 /** \brief obtain type of reference element */ type() const325 GeometryType type () const { return GeometryTypes::cube(dimension); } 326 327 /** 328 * \brief map i-th corner 329 * 330 * \param[in] i corner index 331 * 332 * \returns mapped index 333 */ operator ()(int i) const334 int operator() ( int i ) const { return (i ^ aluTwist_); } 335 operator ()(int i,int codim) const336 int operator() ( int i, int codim ) const { return (codim == 0 ? i : (*this)( i )); } 337 338 /** \} */ 339 340 /** 341 * \brief Geometric Equivalent 342 * \{ 343 */ 344 345 /** \brief cast into affine geometry \f$F\f$ */ 346 template< class ctype > operator AffineGeometry<ctype,dimension,dimension>() const347 operator AffineGeometry< ctype, dimension, dimension > () const 348 { 349 const CoordVector< ctype > coordVector( *this ); 350 return AffineGeometry< ctype, dimension, dimension >( type(), coordVector ); 351 } 352 353 /** \brief equivalent to \f$|DF| > 0\f$ */ positive() const354 bool positive () const { return (aluTwist_ == 0); } 355 356 /** \} */ 357 operator int() const358 operator int () const { return aluTwist_; } 359 360 private: 361 int aluTwist_; 362 }; 363 364 365 366 367 // ALUTwists for dimension 2 368 // ------------------------- 369 370 template< int corners > 371 class ALUTwists< corners, 2 > 372 { 373 public: 374 /** \brief dimension */ 375 static const int dimension = 2; 376 377 typedef ALUTwist< corners, 2 > Twist; 378 379 typedef ALUTwistIterator< Twist > Iterator; 380 381 /** \brief obtain type of reference element */ type() const382 GeometryType type () const 383 { 384 if( corners == 3 ) 385 return GeometryTypes::simplex(dimension); 386 else 387 return GeometryTypes::cube(dimension); 388 } 389 390 /** \brief obtain number of twists in the group */ size() const391 std::size_t size () const { return 2*corners; } 392 393 /** \brief obtain index of a given twist */ index(const Twist & twist) const394 std::size_t index ( const Twist &twist ) const { return static_cast< int >( twist ) + corners; } 395 396 /** 397 * \name Iterators 398 * \{ 399 */ 400 401 /** \brief obtain begin iterator */ begin() const402 Iterator begin () const { return Iterator( Twist( -corners ) ); } 403 404 /** \brief obtain end iterator */ end() const405 Iterator end () const { return Iterator( Twist( corners ) ); } 406 407 template< class Permutation > find(const Permutation & permutation) const408 Iterator find ( const Permutation &permutation ) const 409 { 410 // calculate twist (if permutation is a valid twist) 411 const int origin = duneVertex2aluVertex( permutation( aluVertex2duneVertex( 0 ) ) ); 412 const int next = duneVertex2aluVertex( permutation( aluVertex2duneVertex( 1 ) ) ); 413 const Twist twist( origin, (origin + 1) % corners == next ); 414 415 // check twist equals permutation (i.e., permutation is valid) 416 for( int i = 0; i < corners; ++i ) 417 { 418 if( twist( i ) != permutation( i ) ) 419 return end(); 420 } 421 return Iterator( twist ); 422 } 423 424 /** \} */ 425 426 private: aluVertex2duneVertex(int i)427 static int aluVertex2duneVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); } duneVertex2aluVertex(int i)428 static int duneVertex2aluVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); } 429 swap23(int i)430 static int swap23 ( int i ) { return i ^ (i >> 1); } 431 }; 432 433 434 435 // ALUTwists for dimension 1 436 // ------------------------- 437 438 template<> 439 class ALUTwists< 2, 1 > 440 { 441 public: 442 /** \brief dimension */ 443 static const int dimension = 1; 444 445 typedef ALUTwist< 2, 1 > Twist; 446 447 typedef ALUTwistIterator< Twist > Iterator; 448 449 /** \brief obtain type of reference element */ type() const450 GeometryType type () const { return GeometryTypes::cube(dimension); } 451 452 /** \brief obtain number of twists in the group */ size() const453 std::size_t size () const { return 2; } 454 455 /** \brief obtain index of a given twist */ index(const Twist & twist) const456 std::size_t index ( const Twist &twist ) const { return static_cast< int >( twist ); } 457 458 /** 459 * \name Iterators 460 * \{ 461 */ 462 463 /** \brief obtain begin iterator */ begin() const464 Iterator begin () const { return Iterator( Twist( 0 ) ); } 465 466 /** \brief obtain end iterator */ end() const467 Iterator end () const { return Iterator( Twist( int( size() ) ) ); } 468 469 template< class Permutation > find(const Permutation & permutation) const470 Iterator find ( const Permutation &permutation ) const { return Iterator( Twist( permutation( 0 ) ) ); } 471 472 /** \} */ 473 }; 474 475 476 477 // TrivialTwist 478 // ------------ 479 480 template< unsigned int topologyId, int dim > 481 struct TrivialTwist 482 { 483 /** \brief dimension */ 484 static const int dimension = dim; 485 486 /** 487 * \name Construction 488 * \{ 489 */ 490 491 /** \brief default constructor; results in identity twist */ TrivialTwistDune::TrivialTwist492 TrivialTwist () {} 493 TrivialTwistDune::TrivialTwist494 explicit TrivialTwist ( GeometryType ) {} 495 496 /** \} */ 497 498 /** 499 * \name Copying and Assignment 500 * \{ 501 */ 502 503 /** \brief copy constructor */ TrivialTwistDune::TrivialTwist504 TrivialTwist ( const TrivialTwist & ) {} 505 506 /** \brief assignment operator */ operator =Dune::TrivialTwist507 TrivialTwist &operator= ( const TrivialTwist & ) { return *this; } 508 509 /** \} */ 510 511 /** 512 * \name Group Operations 513 * \{ 514 */ 515 516 /** \brief composition */ operator *Dune::TrivialTwist517 TrivialTwist operator* ( const TrivialTwist & ) const { return *this; } 518 519 /** \brief composition with inverse */ operator /Dune::TrivialTwist520 TrivialTwist operator/ ( const TrivialTwist & ) const { return *this; } 521 522 /** \brief composition */ operator *=Dune::TrivialTwist523 TrivialTwist &operator*= ( const TrivialTwist & ) const { return *this; } 524 525 /** \brief composition with inverse */ operator /=Dune::TrivialTwist526 TrivialTwist &operator/= ( const TrivialTwist & ) const { return *this; } 527 528 /** \brief obtain inverse */ inverseDune::TrivialTwist529 TrivialTwist inverse () const { return *this; } 530 531 /** \} */ 532 533 /** 534 * \brief Comparison 535 * \{ 536 */ 537 538 /** \brief check for equality */ operator ==Dune::TrivialTwist539 bool operator== ( const TrivialTwist & ) const { return true; } 540 541 /** \brief check for inequality */ operator !=Dune::TrivialTwist542 bool operator!= ( const TrivialTwist & ) const { return false; } 543 544 /** \} */ 545 546 /** 547 * \brief Topological Corner Mapping 548 * \{ 549 */ 550 551 /** \brief obtain type of reference element */ typeDune::TrivialTwist552 GeometryType type () const { return GeometryType( topologyId, dim ); } 553 554 /** 555 * \brief map i-th corner 556 * 557 * \param[in] i corner index 558 * 559 * \returns mapped index 560 */ operator ()Dune::TrivialTwist561 int operator() ( int i ) const { return i; } 562 operator ()Dune::TrivialTwist563 int operator() ( int i, int codim ) const { return i; } 564 565 /** \} */ 566 567 /** 568 * \brief Geometric Equivalent 569 * \{ 570 */ 571 572 /** \brief cast into affine geometry \f$F\f$ */ 573 template< class ctype > operator AffineGeometry<ctype,dimension,dimension>Dune::TrivialTwist574 operator AffineGeometry< ctype, dimension, dimension > () const 575 { 576 return ReferenceElements< ctype, dimension >::general( type() ).template geometry< 0 >( 0 ); 577 } 578 579 /** \brief equivalent to \f$|DF| > 0\f$ */ positiveDune::TrivialTwist580 bool positive () const { return true; } 581 582 /** \} */ 583 operator intDune::TrivialTwist584 operator int () const { return 0; } 585 }; 586 587 588 589 // TrivialTwists 590 // ------------- 591 592 template< unsigned int topologyId, int dim > 593 struct TrivialTwists 594 : private TrivialTwist< topologyId, dim > 595 { 596 /** \brief dimension */ 597 static const int dimension = dim; 598 599 typedef TrivialTwist< topologyId, dim > Twist; 600 601 typedef const Twist *Iterator; 602 TrivialTwistsDune::TrivialTwists603 TrivialTwists () {} TrivialTwistsDune::TrivialTwists604 explicit TrivialTwists ( GeometryType type ) {} 605 606 /** \brief obtain type of reference element */ typeDune::TrivialTwists607 GeometryType type () const { return twist().type(); } 608 609 /** \brief obtain number of twists in the group */ sizeDune::TrivialTwists610 static std::size_t size () { return 1; } 611 612 /** \brief obtain index of a given twist */ indexDune::TrivialTwists613 static std::size_t index ( const Twist & ) { return 0; } 614 615 /** 616 * \name Iterators 617 * \{ 618 */ 619 620 /** \brief obtain begin iterator */ beginDune::TrivialTwists621 Iterator begin () const { return &twist(); } 622 623 /** \brief obtain end iterator */ endDune::TrivialTwists624 Iterator end () const { return begin() + size(); } 625 626 template< class Permutation > findDune::TrivialTwists627 Iterator find ( const Permutation &permutation ) const // noexcept 628 { 629 // check whether permutation is the identity (i.e., permutation is valid) 630 const int corners = ReferenceElements< double, dimension >::general( type() ).size( dimension ); 631 for( int i = 0; i < corners; ++i ) 632 { 633 if( i != permutation( i ) ) 634 return end(); 635 } 636 return begin(); 637 } 638 639 /** \} */ 640 641 private: twistDune::TrivialTwists642 const Twist &twist () const { return *this; } 643 }; 644 645 } // namespace Dune 646 647 #endif // #ifndef DUNE_ALUGRID_COMMON_TWISTS_HH 648