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_DOFADMIN_HH
4 #define DUNE_ALBERTA_DOFADMIN_HH
5 
6 #include <utility>
7 
8 #include <dune/common/hybridutilities.hh>
9 
10 #include <dune/grid/albertagrid/misc.hh>
11 #include <dune/grid/albertagrid/elementinfo.hh>
12 
13 #if HAVE_ALBERTA
14 
15 namespace Dune
16 {
17 
18   namespace Alberta
19   {
20 
21     // External Forward Declarations
22     // -----------------------------
23 
24     template< int dim >
25     class MeshPointer;
26 
27 
28 
29     // DofAccess
30     // ---------
31 
32     template< int dim, int codim >
33     class DofAccess
34     {
35       static const int codimtype = CodimType< dim, codim >::value;
36 
37     public:
38       static const int numSubEntities = NumSubEntities< dim, codim >::value;
39 
40       static const int dimension = dim;
41       static const int codimension = codim;
42 
43       typedef Alberta::ElementInfo< dimension > ElementInfo;
44 
DofAccess()45       DofAccess ()
46         : node_( -1 )
47       {}
48 
DofAccess(const DofSpace * dofSpace)49       explicit DofAccess ( const DofSpace *dofSpace )
50       {
51         assert( dofSpace );
52         node_ = dofSpace->admin->mesh->node[ codimtype ];
53         index_ = dofSpace->admin->n0_dof[ codimtype ];
54       }
55 
operator ()(const Element * element,int subEntity,int i) const56       int operator() ( const Element *element, int subEntity, int i ) const
57       {
58         assert( element );
59         assert( node_ != -1 );
60         assert( subEntity < numSubEntities );
61         return element->dof[ node_ + subEntity ][ index_ + i ];
62       }
63 
operator ()(const Element * element,int subEntity) const64       int operator() ( const Element *element, int subEntity ) const
65       {
66         return (*this)( element, subEntity, 0 );
67       }
68 
operator ()(const ElementInfo & elementInfo,int subEntity,int i) const69       int operator() ( const ElementInfo &elementInfo, int subEntity, int i ) const
70       {
71         return (*this)( elementInfo.el(), subEntity, i );
72       }
73 
operator ()(const ElementInfo & elementInfo,int subEntity) const74       int operator() ( const ElementInfo &elementInfo, int subEntity ) const
75       {
76         return (*this)( elementInfo.el(), subEntity );
77       }
78 
79     private:
80       int node_;
81       int index_;
82     };
83 
84 
85 
86     // HierarchyDofNumbering
87     // ---------------------
88 
89     template< int dim >
90     class HierarchyDofNumbering
91     {
92       typedef HierarchyDofNumbering< dim > This;
93 
94     public:
95       static const int dimension = dim;
96 
97       typedef Alberta::MeshPointer< dimension > MeshPointer;
98       typedef Alberta::ElementInfo< dimension > ElementInfo;
99 
100     private:
101       static const int nNodeTypes = N_NODE_TYPES;
102 
103       template< int codim >
104       struct CreateDofSpace;
105 
106       template< int codim >
107       struct CacheDofSpace;
108 
109       typedef std::pair< int, int > Cache;
110 
111     public:
HierarchyDofNumbering()112       HierarchyDofNumbering ()
113       {}
114 
115     private:
116       HierarchyDofNumbering ( const This & );
117       This &operator= ( const This & );
118 
119     public:
~HierarchyDofNumbering()120       ~HierarchyDofNumbering ()
121       {
122         release();
123       }
124 
operator ()(const Element * element,int codim,unsigned int subEntity) const125       int operator() ( const Element *element, int codim, unsigned int subEntity ) const
126       {
127         assert( !(*this) == false );
128         assert( (codim >= 0) && (codim <= dimension) );
129         const Cache &cache = cache_[ codim ];
130         return element->dof[ cache.first + subEntity ][ cache.second ];
131       }
132 
operator ()(const ElementInfo & element,int codim,unsigned int subEntity) const133       int operator() ( const ElementInfo &element, int codim, unsigned int subEntity ) const
134       {
135         return (*this)( element.el(), codim, subEntity );
136       }
137 
operator bool() const138       explicit operator bool () const
139       {
140         return (bool)mesh_;
141       }
142 
dofSpace(int codim) const143       const DofSpace *dofSpace ( int codim ) const
144       {
145         assert( *this );
146         assert( (codim >= 0) && (codim <= dimension) );
147         return dofSpace_[ codim ];
148       }
149 
emptyDofSpace() const150       const DofSpace *emptyDofSpace () const
151       {
152         assert( *this );
153         return emptySpace_;
154       }
155 
mesh() const156       const MeshPointer &mesh () const
157       {
158         return mesh_;
159       }
160 
size(int codim) const161       int size ( int codim ) const
162       {
163         return dofSpace( codim )->admin->size;
164       }
165 
166       void create ( const MeshPointer &mesh );
167 
release()168       void release ()
169       {
170         if( *this )
171         {
172           for( int codim = 0; codim <= dimension; ++codim )
173             freeDofSpace( dofSpace_[ codim ] );
174           freeDofSpace( emptySpace_ );
175           mesh_ = MeshPointer();
176         }
177       }
178 
179     private:
180       static const DofSpace *createEmptyDofSpace ( const MeshPointer &mesh );
181       static const DofSpace *createDofSpace ( const MeshPointer &mesh,
182                                               const std::string &name,
183                                               const int (&ndof)[ nNodeTypes ],
184                                               const bool periodic = false );
185       static void freeDofSpace ( const DofSpace *dofSpace );
186 
187       MeshPointer mesh_;
188       const DofSpace *emptySpace_;
189       const DofSpace *dofSpace_[ dimension+1 ];
190       Cache cache_[ dimension+1 ];
191     };
192 
193 
194 
195     template< int dim >
196     inline void
create(const MeshPointer & mesh)197     HierarchyDofNumbering< dim >::create ( const MeshPointer &mesh )
198     {
199       release();
200 
201       if( !mesh )
202         return;
203 
204       mesh_ = mesh;
205 
206       Hybrid::forEach( std::make_index_sequence< dimension+1 >{}, [ & ]( auto i ){ CreateDofSpace< i >::apply( mesh_, dofSpace_ ); } );
207       Hybrid::forEach( std::make_index_sequence< dimension+1 >{}, [ & ]( auto i ){ CacheDofSpace< i >::apply( dofSpace_, cache_ ); } );
208 
209       emptySpace_ = createEmptyDofSpace( mesh_ );
210       for( int i = 0; i < nNodeTypes; ++i )
211         assert( emptySpace_->admin->n_dof[ i ] == 0 );
212     }
213 
214 
215 
216     template< int dim >
217     inline const DofSpace *
createEmptyDofSpace(const MeshPointer & mesh)218     HierarchyDofNumbering< dim >::createEmptyDofSpace ( const MeshPointer &mesh )
219     {
220       int ndof[ nNodeTypes ];
221       for( int i = 0; i < nNodeTypes; ++i )
222         ndof[ i ] = 0;
223       std::string name = "Empty";
224       return createDofSpace( mesh, name, ndof );
225     }
226 
227 
228     template< int dim >
229     inline const DofSpace *
createDofSpace(const MeshPointer & mesh,const std::string & name,const int (& ndof)[nNodeTypes],const bool periodic)230     HierarchyDofNumbering< dim >::createDofSpace ( const MeshPointer &mesh,
231                                                    const std::string &name,
232                                                    const int (&ndof)[ nNodeTypes ],
233                                                    const bool periodic )
234     {
235       const ALBERTA FLAGS flags
236         = ADM_PRESERVE_COARSE_DOFS | (periodic ? ADM_PERIODIC : 0);
237       return ALBERTA get_dof_space ( mesh, name.c_str(), ndof, flags );
238     }
239 
240 
241     template< int dim >
242     inline void
freeDofSpace(const DofSpace * dofSpace)243     HierarchyDofNumbering< dim >::freeDofSpace ( const DofSpace *dofSpace )
244     {
245       ALBERTA free_fe_space( dofSpace );
246     }
247 
248 
249 
250     // HierarchyDofNumbering::CreateDofSpace
251     // -------------------------------------
252 
253     template< int dim >
254     template< int codim >
255     struct HierarchyDofNumbering< dim >::CreateDofSpace
256     {
applyDune::Alberta::HierarchyDofNumbering::CreateDofSpace257       static void apply ( const MeshPointer &mesh, const DofSpace *(&dofSpace)[ dim+1 ] )
258       {
259         int ndof[ nNodeTypes ];
260         for( int i = 0; i < nNodeTypes; ++i )
261           ndof[ i ] = 0;
262         ndof[ CodimType< dim, codim >::value ] = 1;
263 
264         std::string name = "Codimension ";
265         name += (char)(codim + '0');
266 
267         dofSpace[ codim ] = createDofSpace( mesh, name, ndof );
268         assert( dofSpace[ codim ] );
269       }
270     };
271 
272 
273 
274     // HierarchyDofNumbering::CacheDofSpace
275     // ------------------------------------
276 
277     template< int dim >
278     template< int codim >
279     struct HierarchyDofNumbering< dim >::CacheDofSpace
280     {
applyDune::Alberta::HierarchyDofNumbering::CacheDofSpace281       static void apply ( const DofSpace *(&dofSpace)[ dim+1 ], Cache (&cache)[ dim+1 ] )
282       {
283         assert( dofSpace[ codim ] );
284         const int codimtype = CodimType< dim, codim >::value;
285         cache[ codim ].first = dofSpace[ codim ]->mesh->node[ codimtype ];
286         cache[ codim ].second = dofSpace[ codim ]->admin->n0_dof[ codimtype ];
287       }
288     };
289 
290   } // namespace Alberta
291 
292 } // namespace Dune
293 
294 #endif // #if HAVE_ALBERTA
295 
296 #endif // #ifndef DUNE_ALBERTA_DOFADMIN_HH
297