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