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