1 #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH 2 #define DUNE_ALUGRIDGEOMETRYSTORAGE_HH 3 4 #include <array> 5 #include <memory> 6 7 #include <dune/common/exceptions.hh> 8 9 #include <dune/grid/common/grid.hh> 10 #include <dune/grid/common/gridfactory.hh> 11 12 #include <dune/alugrid/common/declaration.hh> 13 #include <dune/alugrid/3d/alu3dinclude.hh> 14 15 namespace Dune 16 { 17 template< class Grid > 18 class ReferenceGridFactory; 19 20 21 template< class GridImp, class GeometryImpl, int nChild > 22 class ALULocalGeometryStorage 23 { 24 typedef ALULocalGeometryStorage< GridImp, GeometryImpl, nChild > ThisType; 25 26 // array with pointers to the geometries 27 std::array< GeometryImpl, nChild > geoms_; 28 29 // count local geometry creation 30 int count_; 31 32 // true if geoms have been initialized 33 bool initialized_; 34 35 // type of grid impl 36 typedef typename GridImp :: ctype ctype; 37 enum{ dimension = GridImp :: dimension }; 38 enum{ dimensionworld = GridImp :: dimensionworld }; 39 40 template <int dummy, int dim, int dimworld, int > 41 struct CreateGeometries; 42 43 template <int dummy, int dimworld> 44 struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE triangle > 45 { 46 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries47 static inline void createGeometries(Storage& storage, 48 const GeometryType& type, 49 const bool nonConform ) 50 { 51 if( nonConform ) 52 { 53 typedef ALUGrid< 2, dimworld, simplex, nonconforming, ALUGridNoComm > Grid; 54 storage.template createGeometries< Grid > (type); 55 } 56 else 57 { 58 typedef ALUGrid< 2, dimworld, simplex, conforming, ALUGridNoComm > Grid; 59 storage.template createGeometries< Grid > (type); 60 } 61 } 62 }; 63 64 template <int dummy, int dimworld> 65 struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE tetra > 66 { 67 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries68 static inline void createGeometries(Storage& storage, 69 const GeometryType& type, 70 const bool nonConform ) 71 { 72 if( nonConform ) 73 { 74 typedef ALUGrid< 2, dimworld, simplex, nonconforming, ALUGridNoComm > Grid; 75 storage.template createGeometries< Grid > (type); 76 } 77 else 78 { 79 typedef ALUGrid< 2, dimworld, simplex, conforming, ALUGridNoComm > Grid; 80 storage.template createGeometries< Grid > (type); 81 } 82 } 83 }; 84 85 template <int dummy> 86 struct CreateGeometries<dummy, 3, 3, ALU3DSPACE tetra > 87 { 88 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries89 static inline void createGeometries(Storage& storage, 90 const GeometryType& type, 91 const bool nonConform ) 92 { 93 alugrid_assert ( nonConform ) ; 94 { 95 typedef ALUGrid< 3, 3, simplex, nonconforming, ALUGridNoComm > Grid; 96 storage.template createGeometries< Grid > (type); 97 } 98 /* 99 // TODO, implement this for refinement of all edges (conforming) 100 else 101 { 102 typedef ALUGrid< 3, 3, simplex, conforming, ALUGridNoComm > Grid; 103 storage.template createGeometries< Grid > (type); 104 } 105 */ 106 } 107 }; 108 109 template <int dummy, int dimworld> 110 struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE quadrilateral > 111 { 112 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries113 static void createGeometries(Storage& storage, 114 const GeometryType& type, 115 const bool nonConform ) 116 { 117 alugrid_assert ( nonConform ) ; 118 { 119 typedef ALUGrid< 2, dimworld, cube, nonconforming, ALUGridNoComm > Grid; 120 storage.template createGeometries< Grid > (type); 121 } 122 } 123 }; 124 125 template <int dummy, int dimworld> 126 struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE hexa > 127 { 128 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries129 static void createGeometries(Storage& storage, 130 const GeometryType& type, 131 const bool nonConform ) 132 { 133 alugrid_assert ( nonConform ) ; 134 { 135 typedef ALUGrid< 2, dimworld, cube, nonconforming, ALUGridNoComm > Grid; 136 storage.template createGeometries< Grid > (type); 137 } 138 } 139 }; 140 141 template <int dummy> 142 struct CreateGeometries<dummy, 3, 3, ALU3DSPACE hexa > 143 { 144 template <class Storage> createGeometriesDune::ALULocalGeometryStorage::CreateGeometries145 static void createGeometries(Storage& storage, 146 const GeometryType& type, 147 const bool nonConform ) 148 { 149 alugrid_assert ( nonConform ); 150 { 151 typedef ALUGrid< 3, 3, cube, nonconforming, ALUGridNoComm > Grid; 152 storage.template createGeometries< Grid > (type); 153 } 154 } 155 }; 156 157 public: 158 // create empty storage ALULocalGeometryStorage(const GeometryType type,const bool nonConform)159 inline ALULocalGeometryStorage ( const GeometryType type, const bool nonConform ) 160 : count_( 0 ), initialized_( false ) 161 { 162 // initialize geometries 163 initialize( type, nonConform ); 164 } 165 166 // create empty storage ALULocalGeometryStorage()167 inline ALULocalGeometryStorage () 168 : count_( 0 ), initialized_( false ) 169 { 170 } 171 172 // return reference to local geometry operator [](int child) const173 inline const GeometryImpl& operator [] (int child) const 174 { 175 alugrid_assert ( geomCreated(child) ); 176 // this method is not thread safe yet 177 return geoms_[child]; 178 } 179 180 //! access local geometry storage storage(const GeometryType type,const bool nonConforming)181 static inline const ThisType& storage( const GeometryType type, const bool nonConforming ) 182 { 183 if( type.isSimplex() ) 184 { 185 // create static variable for this thread 186 static ThisType simplexGeoms( type, nonConforming ); 187 return simplexGeoms ; 188 } 189 else 190 { 191 // should be a cube geometry a this point 192 alugrid_assert( type.isCube() ); 193 194 // create static variable 195 static ThisType cubeGeoms( type, nonConforming ); 196 return cubeGeoms ; 197 } 198 } 199 200 protected: 201 // check if geometry has been created geomCreated(int child) const202 inline bool geomCreated(int child) const { return geoms_[child].valid(); } 203 204 //! initialize local geometries initialize(const GeometryType type,const bool nonConform)205 inline bool initialize( const GeometryType type, const bool nonConform ) 206 { 207 if( ! initialized_ ) 208 { 209 // first set flag, because this method might be called again during 210 // creation of local geometries and then result in an infinite loop 211 initialized_ = true ; 212 213 // the idea is to create a grid containing the reference element, 214 // refine once and the store the father - child relations 215 CreateGeometries<0, dimension, dimensionworld, GridImp :: elementType > 216 ::createGeometries(*this, type, nonConform); 217 return true; 218 } 219 return false; 220 } 221 222 template < class Grid > createGeometries(const GeometryType & type)223 inline void createGeometries(const GeometryType& type) 224 { 225 static bool firstCall = true ; 226 if( firstCall ) 227 { 228 firstCall = false ; 229 230 // create factory for the reference element grid 231 ReferenceGridFactory< Grid > factory; 232 233 const auto& refElem = 234 Dune::ReferenceElements< ctype, dimension >::general( type ); 235 236 // insert vertices 237 FieldVector<ctype, dimensionworld> pos( 0 ); 238 const int vxSize = refElem.size(dimension); 239 for(int i=0; i<vxSize; ++i) 240 { 241 FieldVector<ctype, dimension> position = refElem.position(i, dimension ); 242 // copy position 243 for(int d = 0; d<dimension; ++d ) 244 pos[ d ] = position[ d ]; 245 246 factory.insertVertex( pos ); 247 } 248 249 std::vector< unsigned int > vertices( vxSize ); 250 // create grid with reference element 251 for(size_t i=0; i<vertices.size(); ++i) vertices[ i ] = i; 252 factory.insertElement(type, vertices); 253 254 std::unique_ptr< Grid > gridPtr( factory.createGrid() ); 255 Grid& grid = *gridPtr; 256 257 // refine once to get children in the reference element 258 const int level = 1; 259 grid.globalRefine( level ); 260 261 { 262 typedef typename Grid :: MacroGridView MacroGridView; 263 MacroGridView macroView = grid.template macroGridView< All_Partition > (); 264 typedef typename MacroGridView :: template Codim< 0 > :: Iterator Iterator; 265 266 Iterator it = macroView.template begin<0> (); 267 268 if( it == macroView.template end<0>() ) 269 DUNE_THROW(InvalidStateException,"Empty Grid, should contain at least 1 element"); 270 271 typedef typename Iterator :: Entity EntityType; 272 273 const EntityType& entity = *it; 274 const typename EntityType :: Geometry& geo = entity.geometry(); 275 typedef typename EntityType :: HierarchicIterator HierarchicIteratorType; 276 const HierarchicIteratorType end = entity.hend( level ); 277 278 int childNum = 0; 279 for( HierarchicIteratorType child = entity.hbegin( level ); 280 child != end; ++child, ++childNum ) 281 { 282 create( geo, child->geometry(), childNum ); 283 } 284 } 285 } 286 } 287 288 // create local geometry 289 template< class Geometry > create(const Geometry & father,const Geometry & son,const int child)290 inline void create ( const Geometry &father, 291 const Geometry &son, 292 const int child ) 293 { 294 alugrid_assert ( !geomCreated( child ) ); 295 alugrid_assert ( (child >= 0) && (child < nChild) ); 296 297 alugrid_assert ( (count_ < nChild) ); 298 ++count_; 299 300 geoms_[ child ].buildGeomInFather( father, son ); 301 } 302 303 }; 304 305 } // namespace Dune 306 307 #endif // #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH 308