1 #ifndef DUNE_ALUGRID_STRUCTUREDGRIDFACTORY_HH
2 #define DUNE_ALUGRID_STRUCTUREDGRIDFACTORY_HH
8 #include <dune/common/version.hh>
10 #include <dune/common/classname.hh>
11 #include <dune/common/exceptions.hh>
12 #include <dune/common/fvector.hh>
13 #include <dune/common/shared_ptr.hh>
15 #include <dune/grid/common/gridfactory.hh>
16 #include <dune/grid/common/exceptions.hh>
23 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
25 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
28 #include <dune/grid/io/file/dgfparser/dgfs.hh>
37 template<
class Gr
id >
45 template<
int dim,
int dimworld, ALUGr
idElementType eltype, ALUGr
idRefinementType refineType,
class Comm >
59 template<
class GV, PartitionIteratorType pitype,
class IS =
typename GV::IndexSet >
60 class SimplePartitioner
62 typedef SimplePartitioner< GV, pitype, IS >
This;
66 typedef typename GridView::Grid
Grid;
71 typedef typename IndexSet::IndexType IndexType;
73 static const int dimension = Grid::dimension;
75 typedef typename Grid::template Codim< 0 >::Entity Element;
77 typedef typename Element::Geometry::GlobalCoordinate VertexType;
80 typedef Dune :: CollectiveCommunication< typename MPIHelper :: MPICommunicator >
83 #ifdef USE_ALUGRID_SFC_ORDERING
89 const VertexType& lowerLeft,
const VertexType& upperRight )
91 gridView_( gridView ),
92 indexSet_( gridView_.indexSet() ),
93 pSize_( comm_.size() ),
94 elementCuts_( pSize_, -1 ),
95 #ifdef USE_ALUGRID_SFC_ORDERING
96 sfc_( lowerLeft, upperRight, comm_ ),
98 maxIndex_(
double(indexSet_.size(0)-1) )
101 calculateElementCuts();
105 template<
class Entity >
106 double index(
const Entity &entity )
const
109 #ifdef USE_ALUGRID_SFC_ORDERING
111 VertexType center = entity.geometry().center();
113 return sfc_.index( center );
115 return double(indexSet_.index( entity ));
119 template<
class Entity >
120 int rank(
const Entity &entity )
const
123 #ifdef USE_ALUGRID_SFC_ORDERING
125 VertexType center = entity.geometry().center();
127 const double hidx = sfc_.index( center );
129 const long int index = (hidx * maxIndex_);
131 const long int index = indexSet_.index( entity );
133 return rank( index );
137 int rank(
long int index )
const
139 if( index < elementCuts_[ 0 ] )
return 0;
140 for(
int p=1; p<pSize_; ++p )
142 if( index >= elementCuts_[ p-1 ] && index < elementCuts_[ p ] )
148 void calculateElementCuts()
150 const size_t nElements = indexSet_.size( 0 );
153 const int nRanks = pSize_;
156 const size_t minPerProc = (double(nElements) / double( nRanks ));
157 size_t maxPerProc = minPerProc ;
158 if( nElements % nRanks != 0 )
163 double percentage = (double(nElements) / double( nRanks ));
164 percentage -= minPerProc ;
165 percentage *= nRanks ;
168 size_t elementCount = maxPerProc ;
169 size_t elementNumber = 0;
170 size_t localElementNumber = 0;
171 const int lastRank = nRanks - 1;
173 const size_t size = indexSet_.size( 0 );
174 for(
size_t i=0; i<size; ++i )
176 if( localElementNumber >= elementCount )
178 elementCuts_[ rank ] = i ;
181 if( rank < lastRank ) ++ rank;
184 localElementNumber = 0;
187 if( elementCount == maxPerProc && rank >= percentage )
188 elementCount = minPerProc ;
193 ++localElementNumber;
197 elementCuts_[ lastRank ] = size ;
205 const GridView& gridView_;
206 const IndexSet &indexSet_;
209 std::vector< long int > elementCuts_ ;
211 #ifdef USE_ALUGRID_SFC_ORDERING
215 const double maxIndex_ ;
223 typedef Dune :: CollectiveCommunication< MPICommunicatorType >
230 std::ifstream file( filename.c_str() );
233 DUNE_THROW(InvalidStateException,
"file not found " << filename );
235 return createCubeGrid( file, filename, mpiComm );
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 );
294 template <
class int_t >
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 );
318 template <
class int_t >
321 const FieldVector<ctype,dimworld>& upperRight,
322 const array< int_t, dim>& elements,
326 std::string name(
"Cartesian ALUGrid via SGrid" );
327 return createCubeGridImpl( lowerLeft, upperRight, elements, comm, name );
331 template <
int codim,
class Entity>
334 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
335 return entity.subEntities( codim );
337 return entity.template count< codim > ();
341 template <
class int_t >
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 ) );
456 #endif // #ifndef DUNE_ALUGRID_STRUCTUREDGRIDFACTORY_HH