240 const std::string& name,
244#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
245 static_assert( dim == dimworld,
"YaspGrid is used for creation of the structured grid which only supports dim == dimworld");
248 Dune::dgf::IntervalBlock intervalBlock( input );
249 if( !intervalBlock.isactive() )
251 std::cerr <<
"No interval block found, using default DGF method to create ALUGrid!" << std::endl;
252 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
255 if( intervalBlock.numIntervals() != 1 )
257 std::cerr <<
"ALUGrid creation from YaspGrid can only handle 1 interval block, using default DGF method to create ALUGrid!" << std::endl;
258 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
261 if( intervalBlock.dimw() != dim )
263 std::cerr <<
"ALUGrid creation from YaspGrid only works for dim == dimworld, using default DGF method to create ALUGrid!" << std::endl;
264 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
267 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
273 std::array<int, dim> dims;
274 FieldVector<ctype, dimworld> lowerLeft;
275 FieldVector<ctype, dimworld> upperRight;
276 for(
int i=0; i<dim; ++i )
278 dims[ i ] = interval.n[ i ] ;
279 lowerLeft[ i ] = interval.p[ 0 ][ i ];
280 upperRight[ i ] = interval.p[ 1 ][ i ];
284 comm.broadcast( &dims[ 0 ], dim, 0 );
285 comm.broadcast( &lowerLeft [ 0 ], dim, 0 );
286 comm.broadcast( &upperRight[ 0 ], dim, 0 );
288 std::string nameYasp( name );
289 nameYasp +=
" via YaspGrid";
291 return SGF :: createCubeGridImpl( lowerLeft, upperRight, dims, comm, nameYasp );
297 const FieldVector<ctype,dimworld>& upperRight,
298 const array< int_t, dim>& elements,
302 std::stringstream dgfstream;
303 dgfstream <<
"DGF" << std::endl;
304 dgfstream <<
"Interval" << std::endl;
305 dgfstream << lowerLeft << std::endl;
306 dgfstream << upperRight << std::endl;
307 for(
int i=0; i<dim; ++ i)
308 dgfstream << elements[ i ] <<
" ";
309 dgfstream << std::endl;
310 dgfstream <<
"#" << std::endl;
311 dgfstream <<
"Simplex" << std::endl;
312 dgfstream <<
"#" << std::endl;
314 Dune::GridPtr< Grid > grid( dgfstream, mpiComm );
344 const FieldVector<ctype,dimworld>& upperRight,
345 const array< int_t, dim>& elements,
347 const std::string& name )
349 const int myrank = comm.rank();
351#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
352 typedef YaspGrid< dimworld, EquidistantOffsetCoordinates<double,dimworld> > CartesianGridType ;
353 std::array< int, dim > dims;
354 for(
int i=0; i<dim; ++i ) dims[ i ] = elements[ i ];
358 CartesianGridType sgrid( lowerLeft, upperRight, dims, std::bitset<dim>(0ULL), 1, commSelf );
360 typedef SGrid< dim, dimworld, double > CartesianGridType ;
361 FieldVector< int, dim > dims;
362 for(
int i=0; i<dim; ++i ) dims[ i ] = elements[ i ];
365 CartesianGridType sgrid( dims, lowerLeft, upperRight );
368 typedef typename CartesianGridType :: LeafGridView GridView ;
369 typedef typename GridView :: IndexSet IndexSet ;
370 typedef typename IndexSet :: IndexType IndexType ;
371 typedef typename GridView :: template Codim< 0 > :: Iterator ElementIterator ;
372 typedef typename ElementIterator::Entity Entity ;
373 typedef typename GridView :: IntersectionIterator IntersectionIterator ;
374 typedef typename IntersectionIterator :: Intersection Intersection ;
376 GridView gridView = sgrid.leafGridView();
377 const IndexSet &indexSet = gridView.indexSet();
380 SimplePartitioner< GridView, InteriorBorder_Partition > partitioner( gridView, comm, lowerLeft, upperRight );
383 GridFactory< Grid > factory;
386 std::map< IndexType, unsigned int > vtxMap;
387 std::map< double, const Entity > sortedElementList;
389 const int numVertices = (1 << dim);
390 std::vector< unsigned int > vertices( numVertices );
392 const ElementIterator end = gridView.template end< 0 >();
393 for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
395 const Entity &entity = *it;
398 if( partitioner.rank( entity ) != myrank )
401 const double elIndex = partitioner.index( entity );
402 assert( sortedElementList.find( elIndex ) == sortedElementList.end() );
403 sortedElementList.insert( std::make_pair( elIndex, entity ) );
406 int nextElementIndex = 0;
407 const auto endSorted = sortedElementList.end();
408 for(
auto it = sortedElementList.begin(); it != endSorted; ++it )
410 const Entity &entity = (*it).second;
413 const typename Entity::Geometry geo = entity.geometry();
415 for(
int i = 0; i < numVertices; ++i )
417 const IndexType vtxId = indexSet.subIndex( entity, i, dim );
419 std::pair< typename std::map< IndexType, unsigned int >::iterator,
bool > result
420 = vtxMap.insert( std::make_pair( vtxId, vtxMap.size() ) );
422 factory.insertVertex( geo.corner( i ), vtxId );
423 vertices[ i ] = result.first->second;
426 factory.insertElement( entity.type(), vertices );
427 const int elementIndex = nextElementIndex++;
431 const IntersectionIterator iend = gridView.iend( entity );
432 for( IntersectionIterator iit = gridView.ibegin( entity ); iit != iend; ++iit )
434 const Intersection &isec = *iit;
435 const int faceNumber = isec.indexInInside();
437 if( isec.boundary() )
438 factory.insertBoundary( elementIndex, faceNumber );
440#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
441 if( isec.neighbor() && (partitioner.rank( isec.outside() ) != myrank) )
443 if( isec.neighbor() && (partitioner.rank( *isec.outside() ) != myrank) )
445 factory.insertProcessBorder( elementIndex, faceNumber );
450 return SharedPtrType( factory.createGrid(
true,
true, name ) );