dune-alugrid  2.6-git
twists.hh
Go to the documentation of this file.
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 
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 >
32  : public std::iterator< std::input_iterator_tag, Twist, int >
33  {
34  explicit ALUTwistIterator ( Twist twist ) : twist_( twist ) {}
35 
36  const Twist &operator* () const { return twist_; }
37  const Twist *operator-> () const { return &twist_; }
38 
39  bool operator== ( const ALUTwistIterator &other ) const { return (twist_ == other.twist_); }
40  bool operator!= ( const ALUTwistIterator &other ) const { return (twist_ != other.twist_); }
41 
42  ALUTwistIterator &operator++ () { ++twist_.aluTwist_; return *this; }
43  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  {
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 
70  CoordVector ( const Twist &twist )
71  : twist_( twist ),
72  refElement_( ReferenceElements< ctype, 2 >::general( twist.type() ) )
73  {}
74 
75  const FieldVector< ctype, 2 > operator[] ( int i ) const { return refElement().position( twist_( i ), 2 ); }
76 
77  const ReferenceElementType& refElement () const { return refElement_; }
78 
79  private:
80  const Twist &twist_;
81  const ReferenceElementType& refElement_;
82  };
83 
84  public:
86  static const int dimension = 2;
87 
94  ALUTwist () : aluTwist_( 0 ) {}
95 
96  explicit ALUTwist ( GeometryType ) : aluTwist_( 0 ) {}
97 
98  explicit ALUTwist ( int aluTwist ) : aluTwist_( aluTwist ) {}
99 
100  ALUTwist ( int origin, bool positive )
101  : aluTwist_( positive ? origin : (origin + corners - 1) % corners - corners )
102  {}
103 
112  ALUTwist ( const ALUTwist &other ) : aluTwist_( other.aluTwist_ ) {}
113 
115  ALUTwist &operator= ( const ALUTwist &other ) { aluTwist_ = other.aluTwist_; return *this; }
116 
125  ALUTwist operator* ( const ALUTwist &other ) const
126  {
127  return ALUTwist( apply( other.apply( 0 ) ), !(positive() ^ other.positive()) );
128  }
129 
131  ALUTwist operator/ ( const ALUTwist &other ) const { return (*this * other.inverse()); }
132 
134  ALUTwist &operator*= ( const ALUTwist &other ) { return (*this = (*this) * other); }
135 
137  ALUTwist &operator/= ( const ALUTwist &other ) { return (*this = (*this) / other); }
138 
140  ALUTwist inverse () const { return ALUTwist( positive() ? (corners - aluTwist_) % corners : aluTwist_ ); }
141 
150  bool operator== ( const ALUTwist &other ) const { return (aluTwist_ == other.aluTwist_); }
151 
153  bool operator!= ( const ALUTwist &other ) const { return (aluTwist_ != other.aluTwist_); }
154 
163  GeometryType type () const
164  {
165  if( corners == 3 )
166  return GeometryType( typename Impl::SimplexTopology< dimension >::type() );
167  else
168  return GeometryType( typename Impl::CubeTopology< dimension >::type() );
169  }
170 
178  int operator() ( int i ) const { return aluVertex2duneVertex( apply( duneVertex2aluVertex( i ) ) ); }
179 
180  int operator() ( int i, int codim ) const { return alu2dune( apply( dune2alu( i, codim ), codim ), codim ); }
181 
190  template< class ctype >
191  operator AffineGeometry< ctype, dimension, dimension > () const
192  {
193  const CoordVector< ctype > coordVector( *this );
194  return AffineGeometry< ctype, dimension, dimension >( coordVector.refElement(), coordVector );
195  }
196 
198  bool positive () const { return (aluTwist_ >= 0); }
199 
202  // non-interface methods
203 
204  operator int () const { return aluTwist_; }
205 
206  // apply twist in ALU numbering
207  int apply ( int i ) const { return ((positive() ? i : 2*corners + 1 - i) + aluTwist_) % corners; }
208  int apply ( int i, int codim ) const { return (codim == 0 ? i : (codim == 1 ? applyEdge( i ) : apply( i ))); }
209 
210  private:
211  int applyEdge ( int i ) const { return ((positive() ? i : 2*corners - i) + aluTwist_) % corners; }
212 
213 
214  static int aluEdge2duneEdge ( int i ) { return ((corners == 3) ? (3 - i) % 3 : (6 - swap23( i )) % 4); }
215  static int duneEdge2aluEdge ( int i ) { return ((corners == 3) ? (3 - i) % 3 : swap23( (6 - i) % 4 )); }
216 
217  static int aluVertex2duneVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); }
218  static int duneVertex2aluVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); }
219 
220  static int alu2dune ( int i, int codim ) { return (codim == 0 ? i : (codim == 1 ? aluEdge2duneEdge( i ) : aluVertex2duneVertex( i ))); }
221  static int dune2alu ( int i, int codim ) { return (codim == 0 ? i : (codim == 1 ? duneEdge2aluEdge( i ) : aluVertex2duneVertex( i ))); }
222 
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  {
237 
238  friend struct ALUTwistIterator< Twist >;
239 
240  template< class ctype >
241  struct CoordVector
242  {
243  CoordVector ( const Twist &twist ) : twist_( twist ) {}
244 
245  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:
253  static const int dimension = 1;
254 
261  ALUTwist () : aluTwist_( 0 ) {}
262 
263  explicit ALUTwist ( GeometryType ) : aluTwist_( 0 ) {}
264 
265  explicit ALUTwist ( int aluTwist ) : aluTwist_( aluTwist ) {}
266 
267  explicit ALUTwist ( bool positive ) : aluTwist_( positive ) {}
268 
277  ALUTwist ( const ALUTwist &other ) : aluTwist_( other.aluTwist_ ) {}
278 
280  ALUTwist &operator= ( const ALUTwist &other ) { aluTwist_ = other.aluTwist_; return *this; }
281 
290  ALUTwist operator* ( const ALUTwist &other ) const { return ALUTwist( aluTwist_ ^ other.aluTwist_ ); }
291 
293  ALUTwist operator/ ( const ALUTwist &other ) const { return (*this * other.inverse()); }
294 
296  ALUTwist &operator*= ( const ALUTwist &other ) { return (*this = (*this) * other); }
297 
299  ALUTwist &operator/= ( const ALUTwist &other ) { return (*this = (*this) / other); }
300 
302  ALUTwist inverse () const { return *this; }
303 
312  bool operator== ( const ALUTwist &other ) const { return (aluTwist_ == other.aluTwist_); }
313 
315  bool operator!= ( const ALUTwist &other ) const { return (aluTwist_ != other.aluTwist_); }
316 
325  GeometryType type () const { return GeometryType( Impl::CubeTopology< dimension >::type() ); }
326 
334  int operator() ( int i ) const { return (i ^ aluTwist_); }
335 
336  int operator() ( int i, int codim ) const { return (codim == 0 ? i : (*this)( i )); }
337 
346  template< class ctype >
347  operator AffineGeometry< ctype, dimension, dimension > () const
348  {
349  const CoordVector< ctype > coordVector( *this );
350  return AffineGeometry< ctype, dimension, dimension >( type(), coordVector );
351  }
352 
354  bool positive () const { return (aluTwist_ == 0); }
355 
358  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:
375  static const int dimension = 2;
376 
378 
380 
382  GeometryType type () const
383  {
384  if( corners == 3 )
385  return GeometryType( typename Impl::SimplexTopology< dimension >::type() );
386  else
387  return GeometryType( typename Impl::CubeTopology< dimension >::type() );
388  }
389 
391  std::size_t size () const { return 2*corners; }
392 
394  std::size_t index ( const Twist &twist ) const { return static_cast< int >( twist ) + corners; }
395 
402  Iterator begin () const { return Iterator( Twist( -corners ) ); }
403 
405  Iterator end () const { return Iterator( Twist( corners ) ); }
406 
407  template< class Permutation >
408  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 
426  private:
427  static int aluVertex2duneVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); }
428  static int duneVertex2aluVertex ( int i ) { return ((corners == 3) ? i : swap23( i )); }
429 
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:
443  static const int dimension = 1;
444 
446 
448 
450  GeometryType type () const { return GeometryType( Impl::CubeTopology< dimension >::type() ); }
451 
453  std::size_t size () const { return 2; }
454 
456  std::size_t index ( const Twist &twist ) const { return static_cast< int >( twist ); }
457 
464  Iterator begin () const { return Iterator( Twist( 0 ) ); }
465 
467  Iterator end () const { return Iterator( Twist( int( size() ) ) ); }
468 
469  template< class Permutation >
470  Iterator find ( const Permutation &permutation ) const { return Iterator( Twist( permutation( 0 ) ) ); }
471 
473  };
474 
475 
476 
477  // TrivialTwist
478  // ------------
479 
480  template< unsigned int topologyId, int dim >
482  {
484  static const int dimension = dim;
485 
493 
494  explicit TrivialTwist ( GeometryType ) {}
495 
504  TrivialTwist ( const TrivialTwist & ) {}
505 
507  TrivialTwist &operator= ( const TrivialTwist & ) { return *this; }
508 
517  TrivialTwist operator* ( const TrivialTwist & ) const { return *this; }
518 
520  TrivialTwist operator/ ( const TrivialTwist & ) const { return *this; }
521 
523  TrivialTwist &operator*= ( const TrivialTwist & ) const { return *this; }
524 
526  TrivialTwist &operator/= ( const TrivialTwist & ) const { return *this; }
527 
529  TrivialTwist inverse () const { return *this; }
530 
539  bool operator== ( const TrivialTwist & ) const { return true; }
540 
542  bool operator!= ( const TrivialTwist & ) const { return false; }
543 
552  GeometryType type () const { return GeometryType( topologyId, dim ); }
553 
561  int operator() ( int i ) const { return i; }
562 
563  int operator() ( int i, int codim ) const { return i; }
564 
573  template< class ctype >
574  operator AffineGeometry< ctype, dimension, dimension > () const
575  {
576  return ReferenceElements< ctype, dimension >::general( type() ).template geometry< 0 >( 0 );
577  }
578 
580  bool positive () const { return true; }
581 
584  operator int () const { return 0; }
585  };
586 
587 
588 
589  // TrivialTwists
590  // -------------
591 
592  template< unsigned int topologyId, int dim >
594  : private TrivialTwist< topologyId, dim >
595  {
597  static const int dimension = dim;
598 
600 
601  typedef const Twist *Iterator;
602 
604  explicit TrivialTwists ( GeometryType type ) {}
605 
607  GeometryType type () const { return twist().type(); }
608 
610  static std::size_t size () { return 1; }
611 
613  static std::size_t index ( const Twist & ) { return 0; }
614 
621  Iterator begin () const { return &twist(); }
622 
624  Iterator end () const { return begin() + size(); }
625 
626  template< class Permutation >
627  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 
641  private:
642  const Twist &twist () const { return *this; }
643  };
644 
645 } // namespace Dune
646 
647 #endif // #ifndef DUNE_ALUGRID_COMMON_TWISTS_HH
Capabilities for ALUGrid.
Definition: alu3dinclude.hh:80
Definition: twists.hh:20
Definition: twists.hh:23
Definition: twists.hh:33
ALUTwistIterator & operator++()
Definition: twists.hh:42
bool operator==(const ALUTwistIterator &other) const
Definition: twists.hh:39
const Twist * operator->() const
Definition: twists.hh:37
const Twist & operator*() const
Definition: twists.hh:36
bool operator!=(const ALUTwistIterator &other) const
Definition: twists.hh:40
ALUTwistIterator(Twist twist)
Definition: twists.hh:34
Definition: twists.hh:56
ALUTwist(int origin, bool positive)
Definition: twists.hh:100
ALUTwist()
default constructor; results in identity twist
Definition: twists.hh:94
ALUTwist(const ALUTwist &other)
copy constructor
Definition: twists.hh:112
int apply(int i) const
Definition: twists.hh:207
ALUTwist(GeometryType)
Definition: twists.hh:96
bool positive() const
equivalent to
Definition: twists.hh:198
GeometryType type() const
Topological Corner Mapping.
Definition: twists.hh:163
ALUTwist inverse() const
obtain inverse
Definition: twists.hh:140
int apply(int i, int codim) const
Definition: twists.hh:208
ALUTwist(int aluTwist)
Definition: twists.hh:98
Definition: twists.hh:235
ALUTwist()
default constructor; results in identity twist
Definition: twists.hh:261
ALUTwist(GeometryType)
Definition: twists.hh:263
ALUTwist(bool positive)
Definition: twists.hh:267
ALUTwist(const ALUTwist &other)
copy constructor
Definition: twists.hh:277
ALUTwist(int aluTwist)
Definition: twists.hh:265
bool positive() const
equivalent to
Definition: twists.hh:354
GeometryType type() const
Topological Corner Mapping.
Definition: twists.hh:325
ALUTwist inverse() const
obtain inverse
Definition: twists.hh:302
std::size_t index(const Twist &twist) const
obtain index of a given twist
Definition: twists.hh:394
Iterator find(const Permutation &permutation) const
Definition: twists.hh:408
Iterator begin() const
obtain begin iterator
Definition: twists.hh:402
ALUTwistIterator< Twist > Iterator
Definition: twists.hh:379
Iterator end() const
obtain end iterator
Definition: twists.hh:405
ALUTwist< corners, 2 > Twist
Definition: twists.hh:377
GeometryType type() const
obtain type of reference element
Definition: twists.hh:382
std::size_t size() const
obtain number of twists in the group
Definition: twists.hh:391
Iterator find(const Permutation &permutation) const
Definition: twists.hh:470
std::size_t index(const Twist &twist) const
obtain index of a given twist
Definition: twists.hh:456
Iterator end() const
obtain end iterator
Definition: twists.hh:467
GeometryType type() const
obtain type of reference element
Definition: twists.hh:450
Iterator begin() const
obtain begin iterator
Definition: twists.hh:464
ALUTwistIterator< Twist > Iterator
Definition: twists.hh:447
std::size_t size() const
obtain number of twists in the group
Definition: twists.hh:453
ALUTwist< 2, 1 > Twist
Definition: twists.hh:445
Definition: twists.hh:482
TrivialTwist & operator=(const TrivialTwist &)
assignment operator
Definition: twists.hh:507
int operator()(int i) const
map i-th corner
Definition: twists.hh:561
TrivialTwist operator*(const TrivialTwist &) const
composition
Definition: twists.hh:517
TrivialTwist(GeometryType)
Definition: twists.hh:494
TrivialTwist inverse() const
obtain inverse
Definition: twists.hh:529
TrivialTwist & operator/=(const TrivialTwist &) const
composition with inverse
Definition: twists.hh:526
static const int dimension
dimension
Definition: twists.hh:484
TrivialTwist & operator*=(const TrivialTwist &) const
composition
Definition: twists.hh:523
bool operator==(const TrivialTwist &) const
Comparison.
Definition: twists.hh:539
bool positive() const
equivalent to
Definition: twists.hh:580
TrivialTwist()
default constructor; results in identity twist
Definition: twists.hh:492
TrivialTwist(const TrivialTwist &)
copy constructor
Definition: twists.hh:504
GeometryType type() const
Topological Corner Mapping.
Definition: twists.hh:552
TrivialTwist operator/(const TrivialTwist &) const
composition with inverse
Definition: twists.hh:520
bool operator!=(const TrivialTwist &) const
check for inequality
Definition: twists.hh:542
Definition: twists.hh:595
Iterator find(const Permutation &permutation) const
Definition: twists.hh:627
static std::size_t index(const Twist &)
obtain index of a given twist
Definition: twists.hh:613
TrivialTwist< topologyId, dim > Twist
Definition: twists.hh:599
static std::size_t size()
obtain number of twists in the group
Definition: twists.hh:610
TrivialTwists(GeometryType type)
Definition: twists.hh:604
Iterator end() const
obtain end iterator
Definition: twists.hh:624
TrivialTwists()
Definition: twists.hh:603
GeometryType type() const
obtain type of reference element
Definition: twists.hh:607
const Twist * Iterator
Definition: twists.hh:601
Iterator begin() const
obtain begin iterator
Definition: twists.hh:621
static const int dimension
dimension
Definition: twists.hh:597