dune-grid
2.2.0
|
00001 #ifndef DUNE_ALBERTA_REFINEMENT_HH 00002 #define DUNE_ALBERTA_REFINEMENT_HH 00003 00010 #include <cassert> 00011 00012 #include <dune/grid/albertagrid/misc.hh> 00013 #include <dune/grid/albertagrid/elementinfo.hh> 00014 00015 #if HAVE_ALBERTA 00016 00017 namespace Dune 00018 { 00019 00020 namespace Alberta 00021 { 00022 00023 // Internal Forward Declarations 00024 // ----------------------------- 00025 00026 template< int dim, int codim > 00027 struct ForEachInteriorSubChild; 00028 00029 00030 00031 // Patch 00032 // ----- 00033 00034 template< int dim > 00035 class Patch 00036 { 00037 typedef Patch< dim > This; 00038 00039 dune_static_assert( ((dim >= 1) && (dim <= 3)), 00040 "Alberta supports only dimensions 1, 2, 3" ); 00041 00042 public: 00043 static const int dimension = dim; 00044 00045 typedef Alberta::ElementInfo< dimension > ElementInfo; 00046 00047 typedef ALBERTA RC_LIST_EL ElementList; 00048 00049 private: 00050 ElementList *list_; 00051 int count_; 00052 00053 public: 00054 Patch ( ElementList *list, int count ) 00055 : list_( list ), 00056 count_( count ) 00057 { 00058 assert( count > 0 ); 00059 } 00060 00061 Element *operator[] ( int i ) const; 00062 00063 int count () const 00064 { 00065 return count_; 00066 } 00067 00068 template< class LevelProvider > 00069 ElementInfo elementInfo ( int i, const LevelProvider &levelProvider ) const; 00070 00071 int elementType ( int i ) const; 00072 bool hasNeighbor ( int i, int neighbor ) const; 00073 int neighborIndex ( int i, int neighbor ) const; 00074 00075 template< class Functor > 00076 void forEach ( Functor &functor ) const 00077 { 00078 for( int i = 0; i < count(); ++i ) 00079 functor( (*this)[ i ] ); 00080 } 00081 00082 template< int codim, class Functor > 00083 void forEachInteriorSubChild ( Functor &functor ) const 00084 { 00085 ForEachInteriorSubChild< dim, codim >::apply( functor, *this ); 00086 } 00087 }; 00088 00089 00090 template< int dim > 00091 inline Element *Patch< dim >::operator[] ( int i ) const 00092 { 00093 assert( (i >= 0) && (i < count()) ); 00094 return list_[ i ].el_info.el; 00095 } 00096 00097 00098 template< int dim > 00099 template< class LevelProvider > 00100 inline typename Patch< dim >::ElementInfo 00101 Patch< dim >::elementInfo ( int i, const LevelProvider &levelProvider ) const 00102 { 00103 assert( (i >= 0) && (i < count()) ); 00104 return ElementInfo::createFake( list_[ i ].el_info ); 00105 } 00106 00107 template<> 00108 template< class LevelProvider > 00109 inline typename Patch< 2 >::ElementInfo 00110 Patch< 2 >::elementInfo ( int i, const LevelProvider &levelProvider ) const 00111 { 00112 assert( (i >= 0) && (i < count()) ); 00113 const MeshPointer< 2 > &mesh = levelProvider.mesh(); 00114 const Element *element = (*this)[ i ]; 00115 const int level = levelProvider( element ); 00116 return ElementInfo::createFake( mesh, element, level ); 00117 } 00118 00119 00120 template< int dim > 00121 inline int Patch< dim >::elementType ( int i ) const 00122 { 00123 assert( (i >= 0) && (i < count()) ); 00124 return list_[ i ].el_info.el_type; 00125 } 00126 00127 00128 template< int dim > 00129 inline bool Patch< dim >::hasNeighbor ( int i, int neighbor ) const 00130 { 00131 return (list_[ i ].neigh[ neighbor ] != NULL); 00132 } 00133 00134 template< int dim > 00135 inline int Patch< dim >::neighborIndex ( int i, int neighbor ) const 00136 { 00137 assert( hasNeighbor( i, neighbor ) ); 00138 return (list_[ i ].neigh[ neighbor ]->no); 00139 } 00140 00141 00142 00143 // ForEachInteriorSubChild 00144 // ----------------------- 00145 00146 template< int dim > 00147 struct ForEachInteriorSubChild< dim, 0 > 00148 { 00149 template< class Functor > 00150 static void apply ( Functor &functor, const Patch< dim > &patch ) 00151 { 00152 for( int i = 0; i < patch.count(); ++i ) 00153 { 00154 Element *const father = patch[ i ]; 00155 functor( father->child[ 0 ], 0 ); 00156 functor( father->child[ 1 ], 0 ); 00157 } 00158 } 00159 }; 00160 00161 template< int dim > 00162 struct ForEachInteriorSubChild< dim, dim > 00163 { 00164 template< class Functor > 00165 static void apply ( Functor &functor, const Patch< dim > &patch ) 00166 { 00167 functor( patch[ 0 ]->child[ 0 ], dim ); 00168 } 00169 }; 00170 00171 template<> 00172 struct ForEachInteriorSubChild< 2, 1 > 00173 { 00174 template< class Functor > 00175 static void apply ( Functor &functor, const Patch< 2 > &patch ) 00176 { 00177 // see alberta/src/2d/lagrange_2_2d.c for details 00178 Element *const firstFather = patch[ 0 ]; 00179 00180 Element *const firstChild = firstFather->child[ 0 ]; 00181 functor( firstChild, 0 ); 00182 functor( firstChild, 1 ); 00183 00184 functor( firstFather->child[ 1 ], 1 ); 00185 00186 if( patch.count() > 1 ) 00187 { 00188 Element *const father = patch[ 1 ]; 00189 functor( father->child[ 0 ], 1 ); 00190 } 00191 } 00192 }; 00193 00194 template<> 00195 struct ForEachInteriorSubChild< 3, 1 > 00196 { 00197 template< class Functor > 00198 static void apply ( Functor &functor, const Patch< 3 > &patch ) 00199 { 00200 // see alberta/src/3d/lagrange_3_3d.c for details 00201 Element *const firstFather = patch[ 0 ]; 00202 00203 Element *const firstChild = firstFather->child[ 0 ]; 00204 functor( firstChild, 0 ); 00205 functor( firstChild, 1 ); 00206 functor( firstChild, 2 ); 00207 00208 Element *const secondChild = firstFather->child[ 1 ]; 00209 functor( secondChild, 1 ); 00210 functor( secondChild, 2 ); 00211 00212 for( int i = 1; i < patch.count(); ++i ) 00213 { 00214 Element *const father = patch[ i ]; 00215 const int type = patch.elementType( i ); 00216 00217 int lr_set = 0; 00218 if( patch.hasNeighbor( i, 0 ) && (patch.neighborIndex( i, 0 ) < i) ) 00219 lr_set |= 1; 00220 if( patch.hasNeighbor( i, 1 ) && (patch.neighborIndex( i, 1 ) < i) ) 00221 lr_set |= 2; 00222 assert( lr_set != 0 ); 00223 00224 functor( father->child[ 0 ], 0 ); 00225 switch( lr_set ) 00226 { 00227 case 1: 00228 functor( father->child[ 0 ], 2 ); 00229 functor( father->child[ 1 ], (type == 0 ? 1 : 2) ); 00230 break; 00231 00232 case 2: 00233 functor( father->child[ 0 ], 1 ); 00234 functor( father->child[ 1 ], (type == 0 ? 2 : 1) ); 00235 break; 00236 } 00237 } 00238 } 00239 }; 00240 00241 template<> 00242 struct ForEachInteriorSubChild< 3, 2 > 00243 { 00244 template< class Functor > 00245 static void apply ( Functor &functor, const Patch< 3 > &patch ) 00246 { 00247 // see alberta/src/3d/lagrange_2_3d.c for details 00248 Element *const firstFather = patch[ 0 ]; 00249 00250 Element *const firstChild = firstFather->child[ 0 ]; 00251 functor( firstChild, 2 ); 00252 functor( firstChild, 4 ); 00253 functor( firstChild, 5 ); 00254 00255 functor( firstFather->child[ 1 ], 2 ); 00256 00257 for( int i = 1; i < patch.count(); ++i ) 00258 { 00259 Element *const father = patch[ i ]; 00260 00261 int lr_set = 0; 00262 if( patch.hasNeighbor( i, 0 ) && (patch.neighborIndex( i, 0 ) < i) ) 00263 lr_set = 1; 00264 if( patch.hasNeighbor( i, 1 ) && (patch.neighborIndex( i, 1 ) < i) ) 00265 lr_set += 2; 00266 assert( lr_set != 0 ); 00267 00268 switch( lr_set ) 00269 { 00270 case 1: 00271 functor( father->child[ 0 ], 4 ); 00272 break; 00273 00274 case 2: 00275 functor( father->child[ 0 ], 5 ); 00276 break; 00277 } 00278 } 00279 } 00280 }; 00281 00282 00283 00284 // GeometryInFather 00285 // ---------------- 00286 00287 template< int dim > 00288 struct GeometryInFather; 00289 00290 template<> 00291 struct GeometryInFather< 1 > 00292 { 00293 static const int dim = 1; 00294 00295 typedef Real LocalVector[ dim ]; 00296 00297 static const LocalVector & 00298 coordinate ( int child, int orientation, int i ) 00299 { 00300 static const Real coords[ 2 ][ dim+1 ][ dim ] 00301 = { { {0.0}, {0.5} }, { {0.5}, {1.0} } }; 00302 assert( (i >= 0) && (i <= dim) ); 00303 return coords[ child ][ i ]; 00304 } 00305 }; 00306 00307 template<> 00308 struct GeometryInFather< 2 > 00309 { 00310 static const int dim = 2; 00311 00312 typedef Real LocalVector[ dim ]; 00313 00314 static const LocalVector & 00315 coordinate ( int child, int orientation, int i ) 00316 { 00317 static const Real coords[ 2 ][ dim+1 ][ dim ] 00318 = { { {0.0, 1.0}, {0.0, 0.0}, {0.5, 0.0} }, 00319 { {1.0, 0.0}, {0.0, 1.0}, {0.5, 0.0} } }; 00320 assert( (i >= 0) && (i <= dim) ); 00321 return coords[ child ][ i ]; 00322 } 00323 }; 00324 00325 template<> 00326 struct GeometryInFather< 3 > 00327 { 00328 static const int dim = 3; 00329 00330 typedef Real LocalVector[ dim ]; 00331 00332 static const LocalVector & 00333 coordinate ( int child, int orientation, int i ) 00334 { 00335 static const Real coords[ 2 ][ dim+1 ][ dim ] 00336 = { { {0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {0.5, 0.0, 0.0} }, 00337 { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {0.5, 0.0, 0.0} } }; 00338 static const int flip[ 2 ][ 2 ][ dim+1 ] 00339 = { { {0, 1, 2, 3}, {0, 1, 2, 3} }, { {0, 2, 1, 3}, {0, 1, 2, 3} } }; 00340 assert( (i >= 0) && (i <= dim) ); 00341 i = flip[ child ][ orientation ][ i ]; 00342 return coords[ child ][ i ]; 00343 } 00344 }; 00345 00346 } 00347 00348 } 00349 00350 #endif // #if HAVE_ALBERTA 00351 00352 #endif