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_LEVEL_HH
4 #define DUNE_ALBERTA_LEVEL_HH
5 
6 #include <cassert>
7 #include <cstdlib>
8 
9 #include <dune/grid/albertagrid/meshpointer.hh>
10 #include <dune/grid/albertagrid/dofadmin.hh>
11 #include <dune/grid/albertagrid/dofvector.hh>
12 
13 #if HAVE_ALBERTA
14 
15 namespace Dune
16 {
17 
18   // AlbertaGridLevelProvider
19   // ------------------------
20 
21   template< int dim >
22   class AlbertaGridLevelProvider
23   {
24     typedef AlbertaGridLevelProvider< dim > This;
25 
26     typedef unsigned char Level;
27 
28     typedef Alberta::DofVectorPointer< Level > DofVectorPointer;
29     typedef Alberta::DofAccess< dim, 0 > DofAccess;
30 
31     typedef Alberta::FillFlags< dim > FillFlags;
32 
33     static const Level isNewFlag = (1 << 7);
34     static const Level levelMask = (1 << 7) - 1;
35 
36     class SetLocal;
37     class CalcMaxLevel;
38 
39     template< Level flags >
40     struct ClearFlags;
41 
42     struct Interpolation;
43 
44   public:
45     typedef Alberta::ElementInfo< dim > ElementInfo;
46     typedef Alberta::MeshPointer< dim > MeshPointer;
47     typedef Alberta::HierarchyDofNumbering< dim > DofNumbering;
48 
operator ()(const Alberta::Element * element) const49     Level operator() ( const Alberta::Element *element ) const
50     {
51       const Level *array = (Level *)level_;
52       return array[ dofAccess_( element, 0 ) ] & levelMask;
53     }
54 
operator ()(const ElementInfo & elementInfo) const55     Level operator() ( const ElementInfo &elementInfo ) const
56     {
57       return (*this)( elementInfo.el() );
58     }
59 
isNew(const Alberta::Element * element) const60     bool isNew ( const Alberta::Element *element ) const
61     {
62       const Level *array = (Level *)level_;
63       return ((array[ dofAccess_( element, 0 ) ] & isNewFlag) != 0);
64     }
65 
isNew(const ElementInfo & elementInfo) const66     bool isNew ( const ElementInfo &elementInfo ) const
67     {
68       return isNew( elementInfo.el() );
69     }
70 
maxLevel() const71     Level maxLevel () const
72     {
73       CalcMaxLevel calcFromCache;
74       level_.forEach( calcFromCache );
75 #ifndef NDEBUG
76       CalcMaxLevel calcFromGrid;
77       mesh().leafTraverse( calcFromGrid, FillFlags::nothing );
78       assert( calcFromCache.maxLevel() == calcFromGrid.maxLevel() );
79 #endif
80       return calcFromCache.maxLevel();;
81     }
82 
mesh() const83     MeshPointer mesh () const
84     {
85       return MeshPointer( level_.dofSpace()->mesh );
86     }
87 
markAllOld()88     void markAllOld ()
89     {
90       ClearFlags< isNewFlag > clearIsNew;
91       level_.forEach( clearIsNew );
92     }
93 
create(const DofNumbering & dofNumbering)94     void create ( const DofNumbering &dofNumbering )
95     {
96       const Alberta::DofSpace *const dofSpace = dofNumbering.dofSpace( 0 );
97       dofAccess_ = DofAccess( dofSpace );
98 
99       level_.create( dofSpace, "Element level" );
100       assert( level_ );
101       level_.template setupInterpolation< Interpolation >();
102 
103       SetLocal setLocal( level_ );
104       mesh().hierarchicTraverse( setLocal, FillFlags::nothing );
105     }
106 
release()107     void release ()
108     {
109       level_.release();
110       dofAccess_ = DofAccess();
111     }
112 
113   private:
114     DofVectorPointer level_;
115     DofAccess dofAccess_;
116   };
117 
118 
119 
120   // AlbertaGridLevelProvider::SetLocal
121   // ----------------------------------
122 
123   template< int dim >
124   class AlbertaGridLevelProvider< dim >::SetLocal
125   {
126     DofVectorPointer level_;
127     DofAccess dofAccess_;
128 
129   public:
SetLocal(const DofVectorPointer & level)130     explicit SetLocal ( const DofVectorPointer &level )
131       : level_( level ),
132         dofAccess_( level.dofSpace() )
133     {}
134 
operator ()(const Alberta::ElementInfo<dim> & elementInfo) const135     void operator() ( const Alberta::ElementInfo< dim > &elementInfo ) const
136     {
137       Level *const array = (Level *)level_;
138       array[ dofAccess_( elementInfo, 0 ) ] = elementInfo.level();
139     }
140   };
141 
142 
143 
144   // AlbertaGridLevelProvider::CalcMaxLevel
145   // --------------------------------------
146 
147   template< int dim >
148   class AlbertaGridLevelProvider< dim >::CalcMaxLevel
149   {
150     Level maxLevel_;
151 
152   public:
CalcMaxLevel()153     CalcMaxLevel ()
154       : maxLevel_( 0 )
155     {}
156 
operator ()(const Level & dof)157     void operator() ( const Level &dof )
158     {
159       maxLevel_ = std::max( maxLevel_, Level( dof & levelMask ) );
160     }
161 
operator ()(const Alberta::ElementInfo<dim> & elementInfo)162     void operator() ( const Alberta::ElementInfo< dim > &elementInfo )
163     {
164       maxLevel_ = std::max( maxLevel_, Level( elementInfo.level() ) );
165     }
166 
maxLevel() const167     Level maxLevel () const
168     {
169       return maxLevel_;
170     }
171   };
172 
173 
174 
175   // AlbertaGridLevelProvider::ClearFlags
176   // ------------------------------------
177 
178   template< int dim >
179   template< typename AlbertaGridLevelProvider< dim >::Level flags >
180   struct AlbertaGridLevelProvider< dim >::ClearFlags
181   {
operator ()Dune::AlbertaGridLevelProvider::ClearFlags182     void operator() ( Level &dof ) const
183     {
184       dof &= ~flags;
185     }
186   };
187 
188 
189 
190   // AlbertaGridLevelProvider::Interpolation
191   // ---------------------------------------
192 
193   template< int dim >
194   struct AlbertaGridLevelProvider< dim >::Interpolation
195   {
196     static const int dimension = dim;
197 
198     typedef Alberta::Patch< dimension > Patch;
199 
interpolateVectorDune::AlbertaGridLevelProvider::Interpolation200     static void interpolateVector ( const DofVectorPointer &dofVector,
201                                     const Patch &patch )
202     {
203       const DofAccess dofAccess( dofVector.dofSpace() );
204       Level *array = (Level *)dofVector;
205 
206       for( int i = 0; i < patch.count(); ++i )
207       {
208         const Alberta::Element *const father = patch[ i ];
209         assert( (array[ dofAccess( father, 0 ) ] & levelMask) < levelMask );
210         const Level childLevel = (array[ dofAccess( father, 0 ) ] + 1) | isNewFlag;
211         for( int i = 0; i < 2; ++i )
212         {
213           const Alberta::Element *child = father->child[ i ];
214           array[ dofAccess( child, 0 ) ] = childLevel;
215         }
216       }
217     }
218   };
219 
220 }
221 
222 #endif // #if HAVE_ALBERTA
223 
224 #endif
225