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>
14 #include <dune/common/version.hh>
16 #include <dune/grid/common/gridfactory.hh>
17 #include <dune/grid/common/exceptions.hh>
25 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
33 template<
class Gr
id >
41 template<
int dim,
int dimworld, ALUGr
idElementType eltype, ALUGr
idRefinementType refineType,
class Comm >
46 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6)
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_( SpaceFillingCurveOrderingType::ZCurve, lowerLeft, upperRight, comm_ ),
100 #ifdef USE_ALUGRID_SFC_ORDERING
101 const auto end = gridView_.template end< 0 > ();
102 for(
auto it = gridView_.template begin< 0 > (); it != end; ++it )
104 VertexType center = (*it).geometry().center();
106 const double hidx = sfc_.index( center );
107 maxIndex_ = std::max( maxIndex_, hidx );
111 maxIndex_ /= indexSet_.size( 0 );
115 calculateElementCuts();
119 template<
class Entity >
120 double index(
const Entity &entity )
const
123 #ifdef USE_ALUGRID_SFC_ORDERING
125 VertexType center = entity.geometry().center();
127 return sfc_.index( center );
129 return double(indexSet_.index( entity ));
133 template<
class Entity >
134 int rank(
const Entity &entity )
const
137 #ifdef USE_ALUGRID_SFC_ORDERING
139 VertexType center = entity.geometry().center();
141 const double hidx = sfc_.index( center );
143 const long int index = (hidx / maxIndex_);
146 const long int index = indexSet_.index( entity );
148 return rank( index );
152 int rank(
long int index )
const
154 if( index < elementCuts_[ 0 ] )
return 0;
155 for(
int p=1; p<pSize_; ++p )
157 if( index >= elementCuts_[ p-1 ] && index < elementCuts_[ p ] )
163 void calculateElementCuts()
165 const size_t nElements = indexSet_.size( 0 );
168 const int nRanks = pSize_;
171 const size_t minPerProc = (double(nElements) / double( nRanks ));
172 size_t maxPerProc = minPerProc ;
173 if( nElements % nRanks != 0 )
178 double percentage = (double(nElements) / double( nRanks ));
179 percentage -= minPerProc ;
180 percentage *= nRanks ;
183 size_t elementCount = maxPerProc ;
184 size_t elementNumber = 0;
185 size_t localElementNumber = 0;
186 const int lastRank = nRanks - 1;
188 const size_t size = indexSet_.size( 0 );
189 for(
size_t i=0; i<size; ++i )
191 if( localElementNumber >= elementCount )
193 elementCuts_[ rank ] = i ;
196 if( rank < lastRank ) ++ rank;
199 localElementNumber = 0;
202 if( elementCount == maxPerProc && rank >= percentage )
203 elementCount = minPerProc ;
208 ++localElementNumber;
212 elementCuts_[ lastRank ] = size ;
220 const GridView& gridView_;
221 const IndexSet &indexSet_;
224 std::vector< long int > elementCuts_ ;
226 #ifdef USE_ALUGRID_SFC_ORDERING
238 typedef Dune :: CollectiveCommunication< MPICommunicatorType >
245 std::ifstream file( filename.c_str() );
248 DUNE_THROW(InvalidStateException,
"file not found " << filename );
250 return createCubeGrid( file, filename, mpiComm );
255 const std::string& name,
259 static_assert( dim == dimworld,
"YaspGrid is used for creation of the structured grid which only supports dim == dimworld");
261 Dune::dgf::IntervalBlock intervalBlock( input );
262 if( !intervalBlock.isactive() )
264 std::cerr <<
"No interval block found, using default DGF method to create ALUGrid!" << std::endl;
265 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
268 if( intervalBlock.numIntervals() != 1 )
270 std::cerr <<
"ALUGrid creation from YaspGrid can only handle 1 interval block, using default DGF method to create ALUGrid!" << std::endl;
271 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
274 if( intervalBlock.dimw() != dim )
276 std::cerr <<
"ALUGrid creation from YaspGrid only works for dim == dimworld, using default DGF method to create ALUGrid!" << std::endl;
277 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
280 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
286 std::array<int, dim> dims;
287 FieldVector<ctype, dimworld> lowerLeft;
288 FieldVector<ctype, dimworld> upperRight;
289 for(
int i=0; i<dim; ++i )
291 dims[ i ] = interval.n[ i ] ;
292 lowerLeft[ i ] = interval.p[ 0 ][ i ];
293 upperRight[ i ] = interval.p[ 1 ][ i ];
297 comm.broadcast( &dims[ 0 ], dim, 0 );
298 comm.broadcast( &lowerLeft [ 0 ], dim, 0 );
299 comm.broadcast( &upperRight[ 0 ], dim, 0 );
301 std::string nameYasp( name );
302 nameYasp +=
" via YaspGrid";
304 return SGF :: createCubeGridImpl( lowerLeft, upperRight, dims, comm, nameYasp );
307 template <
class int_t >
310 const FieldVector<ctype,dimworld>& upperRight,
311 const std::array< int_t, dim>& elements,
315 std::stringstream dgfstream;
316 dgfstream <<
"DGF" << std::endl;
317 dgfstream <<
"Interval" << std::endl;
318 dgfstream << lowerLeft << std::endl;
319 dgfstream << upperRight << std::endl;
320 for(
int i=0; i<dim; ++ i)
321 dgfstream << elements[ i ] <<
" ";
322 dgfstream << std::endl;
323 dgfstream <<
"#" << std::endl;
324 dgfstream <<
"Cube" << std::endl;
325 dgfstream <<
"#" << std::endl;
326 dgfstream <<
"Simplex" << std::endl;
327 dgfstream <<
"#" << std::endl;
329 std::cout << dgfstream.str() << std::endl;
331 Dune::GridPtr< Grid > grid( dgfstream, mpiComm );
335 template <
class int_t >
338 const FieldVector<ctype,dimworld>& upperRight,
339 const std::array< int_t, dim>& elements,
343 std::string name(
"Cartesian ALUGrid via YaspGrid" );
344 return createCubeGridImpl( lowerLeft, upperRight, elements, comm, name );
348 template <
int codim,
class Entity>
351 return entity.subEntities( codim );
354 template <
class int_t >
357 const FieldVector<ctype,dimworld>& upperRight,
358 const std::array< int_t, dim>& elements,
360 const std::string& name )
362 const int myrank = comm.rank();
364 typedef YaspGrid< dimworld, EquidistantOffsetCoordinates<double,dimworld> > CartesianGridType ;
365 std::array< int, dim > dims;
366 for(
int i=0; i<dim; ++i ) dims[ i ] = elements[ i ];
370 CartesianGridType sgrid( lowerLeft, upperRight, dims, std::bitset<dim>(0ULL), 1, commSelf );
372 typedef typename CartesianGridType :: LeafGridView GridView ;
373 typedef typename GridView :: IndexSet IndexSet ;
374 typedef typename IndexSet :: IndexType IndexType ;
375 typedef typename GridView :: template Codim< 0 > :: Iterator ElementIterator ;
376 typedef typename ElementIterator::Entity Entity ;
377 typedef typename GridView :: IntersectionIterator IntersectionIterator ;
378 typedef typename IntersectionIterator :: Intersection Intersection ;
380 GridView gridView = sgrid.leafGridView();
381 const IndexSet &indexSet = gridView.indexSet();
384 SimplePartitioner< GridView, InteriorBorder_Partition > partitioner( gridView, comm, lowerLeft, upperRight );
387 GridFactory< Grid > factory;
390 std::map< IndexType, unsigned int > vtxMap;
391 std::map< double, const Entity > sortedElementList;
393 const int numVertices = (1 << dim);
394 std::vector< unsigned int > vertices( numVertices );
396 const ElementIterator end = gridView.template end< 0 >();
397 for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
399 const Entity &entity = *it;
402 if( partitioner.rank( entity ) != myrank )
405 const double elIndex = partitioner.index( entity );
406 assert( sortedElementList.find( elIndex ) == sortedElementList.end() );
407 sortedElementList.insert( std::make_pair( elIndex, entity ) );
410 int nextElementIndex = 0;
411 const auto endSorted = sortedElementList.end();
412 for(
auto it = sortedElementList.begin(); it != endSorted; ++it )
414 const Entity &entity = (*it).second;
417 const typename Entity::Geometry geo = entity.geometry();
419 for(
int i = 0; i < numVertices; ++i )
421 const IndexType vtxId = indexSet.subIndex( entity, i, dim );
423 std::pair< typename std::map< IndexType, unsigned int >::iterator,
bool > result
424 = vtxMap.insert( std::make_pair( vtxId, vtxMap.size() ) );
426 factory.insertVertex( geo.corner( i ), vtxId );
427 vertices[ i ] = result.first->second;
430 factory.insertElement( entity.type(), vertices );
431 const int elementIndex = nextElementIndex++;
435 const IntersectionIterator iend = gridView.iend( entity );
436 for( IntersectionIterator iit = gridView.ibegin( entity ); iit != iend; ++iit )
438 const Intersection &isec = *iit;
439 const int faceNumber = isec.indexInInside();
441 if( isec.boundary() )
442 factory.insertBoundary( elementIndex, faceNumber );
444 if( isec.neighbor() && (partitioner.rank( isec.outside() ) != myrank) )
445 factory.insertProcessBorder( elementIndex, faceNumber );
450 return SharedPtrType( factory.createGrid(
true,
true, name ) );
#define alugrid_assert(EX)
Definition: alugrid_assert.hh:20
Definition: alu3dinclude.hh:50
Definition: alu3dinclude.hh:80
unstructured parallel implementation of the DUNE grid interface
Definition: alugrid.hh:31
BaseType::ctype ctype
Definition: alugrid.hh:48
Definition: structuredgridfactory.hh:34
Dune ::CollectiveCommunication< MPICommunicatorType > CollectiveCommunication
Definition: structuredgridfactory.hh:239
static SharedPtrType createCubeGrid(const std::string &filename, MPICommunicatorType mpiComm=MPIHelper ::getCommunicator())
Definition: structuredgridfactory.hh:242
static SharedPtrType createCubeGrid(std::istream &input, const std::string &name, MPICommunicatorType mpiComm=MPIHelper ::getCommunicator())
Definition: structuredgridfactory.hh:254
MPIHelper ::MPICommunicator MPICommunicatorType
Definition: structuredgridfactory.hh:235
ALUGrid< dim, dimworld, eltype, refineType, Comm > Grid
Definition: structuredgridfactory.hh:45
int subEntities(const Entity &entity) const
Definition: structuredgridfactory.hh:349
Dune::GridPtr< Grid >::mygrid_ptr SharedPtrType
Definition: structuredgridfactory.hh:50
StructuredGridFactory< Grid > This
Definition: structuredgridfactory.hh:54
static SharedPtrType createSimplexGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const std::array< int_t, dim > &elements, MPICommunicatorType mpiComm=MPIHelper ::getCommunicator())
Definition: structuredgridfactory.hh:309
static SharedPtrType createCubeGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const std::array< int_t, dim > &elements, MPICommunicatorType mpiComm=MPIHelper ::getCommunicator())
Definition: structuredgridfactory.hh:337
static SharedPtrType createCubeGridImpl(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const std::array< int_t, dim > &elements, const CollectiveCommunication &comm, const std::string &name)
Definition: structuredgridfactory.hh:356
Grid::ctype ctype
Definition: structuredgridfactory.hh:234