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");
264 dgf::PeriodicFaceTransformationBlock trafoBlock( input, dimworld );
265 if( trafoBlock.isactive() )
267 Dune::GridPtr< Grid > grid( input, mpiComm );
271 Dune::dgf::IntervalBlock intervalBlock( input );
272 if( !intervalBlock.isactive() )
274 std::cerr <<
"No interval block found, using default DGF method to create ALUGrid!" << std::endl;
275 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
278 if( intervalBlock.numIntervals() != 1 )
280 std::cerr <<
"ALUGrid creation from YaspGrid can only handle 1 interval block, using default DGF method to create ALUGrid!" << std::endl;
281 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
284 if( intervalBlock.dimw() != dim )
286 std::cerr <<
"ALUGrid creation from YaspGrid only works for dim == dimworld, using default DGF method to create ALUGrid!" << std::endl;
287 return SharedPtrType( GridPtr< Grid > (input, mpiComm ).release());
290 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
296 std::array<int, dim> dims;
297 FieldVector<ctype, dimworld> lowerLeft;
298 FieldVector<ctype, dimworld> upperRight;
299 for(
int i=0; i<dim; ++i )
301 dims[ i ] = interval.n[ i ] ;
302 lowerLeft[ i ] = interval.p[ 0 ][ i ];
303 upperRight[ i ] = interval.p[ 1 ][ i ];
307 comm.broadcast( &dims[ 0 ], dim, 0 );
308 comm.broadcast( &lowerLeft [ 0 ], dim, 0 );
309 comm.broadcast( &upperRight[ 0 ], dim, 0 );
311 std::string nameYasp( name );
312 nameYasp +=
" via YaspGrid";
314 return SGF :: createCubeGridImpl( lowerLeft, upperRight, dims, comm, nameYasp );
317 template <
class int_t >
320 const FieldVector<ctype,dimworld>& upperRight,
321 const std::array< int_t, dim>& elements,
325 std::stringstream dgfstream;
326 dgfstream <<
"DGF" << std::endl;
327 dgfstream <<
"Interval" << std::endl;
328 dgfstream << lowerLeft << std::endl;
329 dgfstream << upperRight << std::endl;
330 for(
int i=0; i<dim; ++ i)
331 dgfstream << elements[ i ] <<
" ";
332 dgfstream << std::endl;
333 dgfstream <<
"#" << std::endl;
334 dgfstream <<
"Cube" << std::endl;
335 dgfstream <<
"#" << std::endl;
336 dgfstream <<
"Simplex" << std::endl;
337 dgfstream <<
"#" << std::endl;
339 std::cout << dgfstream.str() << std::endl;
341 Dune::GridPtr< Grid > grid( dgfstream, mpiComm );
345 template <
class int_t >
348 const FieldVector<ctype,dimworld>& upperRight,
349 const std::array< int_t, dim>& elements,
353 std::string name(
"Cartesian ALUGrid via YaspGrid" );
354 return createCubeGridImpl( lowerLeft, upperRight, elements, comm, name );
358 template <
int codim,
class Entity>
361 return entity.subEntities( codim );
364 template <
class int_t >
367 const FieldVector<ctype,dimworld>& upperRight,
368 const std::array< int_t, dim>& elements,
370 const std::string& name )
372 const int myrank = comm.rank();
374 typedef YaspGrid< dimworld, EquidistantOffsetCoordinates<double,dimworld> > CartesianGridType ;
375 std::array< int, dim > dims;
376 for(
int i=0; i<dim; ++i ) dims[ i ] = elements[ i ];
380 CartesianGridType sgrid( lowerLeft, upperRight, dims, std::bitset<dim>(0ULL), 1, commSelf );
382 typedef typename CartesianGridType :: LeafGridView GridView ;
383 typedef typename GridView :: IndexSet IndexSet ;
384 typedef typename IndexSet :: IndexType IndexType ;
385 typedef typename GridView :: template Codim< 0 > :: Iterator ElementIterator ;
386 typedef typename ElementIterator::Entity Entity ;
387 typedef typename GridView :: IntersectionIterator IntersectionIterator ;
388 typedef typename IntersectionIterator :: Intersection Intersection ;
390 GridView gridView = sgrid.leafGridView();
391 const IndexSet &indexSet = gridView.indexSet();
394 SimplePartitioner< GridView, InteriorBorder_Partition > partitioner( gridView, comm, lowerLeft, upperRight );
397 GridFactory< Grid > factory;
400 std::map< IndexType, unsigned int > vtxMap;
401 std::map< double, const Entity > sortedElementList;
403 const int numVertices = (1 << dim);
404 std::vector< unsigned int > vertices( numVertices );
406 const ElementIterator end = gridView.template end< 0 >();
407 for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
409 const Entity &entity = *it;
412 if( partitioner.rank( entity ) != myrank )
415 const double elIndex = partitioner.index( entity );
416 assert( sortedElementList.find( elIndex ) == sortedElementList.end() );
417 sortedElementList.insert( std::make_pair( elIndex, entity ) );
420 int nextElementIndex = 0;
421 const auto endSorted = sortedElementList.end();
422 for(
auto it = sortedElementList.begin(); it != endSorted; ++it )
424 const Entity &entity = (*it).second;
427 const typename Entity::Geometry geo = entity.geometry();
429 for(
int i = 0; i < numVertices; ++i )
431 const IndexType vtxId = indexSet.subIndex( entity, i, dim );
433 std::pair< typename std::map< IndexType, unsigned int >::iterator,
bool > result
434 = vtxMap.insert( std::make_pair( vtxId, vtxMap.size() ) );
436 factory.insertVertex( geo.corner( i ), vtxId );
437 vertices[ i ] = result.first->second;
440 factory.insertElement( entity.type(), vertices );
441 const int elementIndex = nextElementIndex++;
445 const IntersectionIterator iend = gridView.iend( entity );
446 for( IntersectionIterator iit = gridView.ibegin( entity ); iit != iend; ++iit )
448 const Intersection &isec = *iit;
449 const int faceNumber = isec.indexInInside();
451 if( isec.boundary() )
452 factory.insertBoundary( elementIndex, faceNumber );
454 if( isec.neighbor() && (partitioner.rank( isec.outside() ) != myrank) )
455 factory.insertProcessBorder( elementIndex, faceNumber );
461 factory.setLongestEdgeFlag(
false);
464 return SharedPtrType( factory.createGrid(
true,
true, name ) );
#define alugrid_assert(EX)
Definition: alugrid_assert.hh:20
Definition: alu3dinclude.hh:33
Definition: alu3dinclude.hh:63
unstructured parallel implementation of the DUNE grid interface
Definition: alugrid.hh:31
BaseType::ctype ctype
Definition: alugrid.hh:46
Definition: structuredgridfactory.hh:34
Dune ::CollectiveCommunication< MPICommunicatorType > CollectiveCommunication
Definition: structuredgridfactory.hh:239
MPIHelper::MPICommunicator MPICommunicatorType
Definition: structuredgridfactory.hh:235
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
ALUGrid< dim, dimworld, eltype, refineType, Comm > Grid
Definition: structuredgridfactory.hh:45
int subEntities(const Entity &entity) const
Definition: structuredgridfactory.hh:359
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:319
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:347
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:366
Grid::ctype ctype
Definition: structuredgridfactory.hh:234
Dune::GridPtr< Grid >::mygrid_ptr SharedPtrType
Definition: structuredgridfactory.hh:50