1 #ifndef DUNE_ALU3DGRID_FACTORY_HH
2 #define DUNE_ALU3DGRID_FACTORY_HH
8 #include <dune/common/shared_ptr.hh>
9 #include <dune/common/parallel/mpihelper.hh>
10 #include <dune/common/version.hh>
12 #include <dune/geometry/referenceelements.hh>
14 #include <dune/grid/common/gridfactory.hh>
15 #include <dune/grid/common/boundaryprojection.hh>
25 template<
class ALUGr
id >
26 class ALU3dGridFactory
27 :
public GridFactoryInterface< ALUGrid >
29 typedef ALU3dGridFactory< ALUGrid > ThisType;
30 typedef GridFactoryInterface< ALUGrid > BaseType;
51 #if !DUNE_VERSION_NEWER(DUNE_GRID,3,0)
53 #endif // #if !DUNE_VERSION_NEWER(DUNE_GRID,3,0)
68 "ALU3dGridFactory supports only grids containing "
69 "tetrahedrons or hexahedrons exclusively." );
71 typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapperType;
77 static const unsigned int boundaryId2d = 87;
83 typedef FieldVector< ctype, dimensionworld > VertexInputType;
87 typedef FieldVector< ctype, 3 > VertexType;
89 typedef std::vector< unsigned int > ElementType;
90 typedef std::array< unsigned int, numFaceCorners > FaceType;
94 typedef std::vector< std::pair< VertexType, GlobalIdType > > VertexVector;
95 typedef std::vector< ElementType > ElementVector;
96 typedef std::pair< FaceType, int > BndPair ;
97 typedef std::map< FaceType, int > BoundaryIdMap;
98 typedef std::vector< std::pair< BndPair, BndPair > > PeriodicBoundaryVector;
99 typedef std::pair< unsigned int, int > SubEntity;
100 typedef std::map< FaceType, SubEntity, FaceLess > FaceMap;
102 typedef std::map< FaceType, const DuneBoundaryProjectionType* > BoundaryProjectionMap;
103 typedef std::vector< const DuneBoundaryProjectionType* > BoundaryProjectionVector;
105 typedef std::vector< Transformation > FaceTransformationVector;
107 static void copy (
const std::initializer_list< unsigned int > &vertices, FaceType &faceId )
109 std::copy_n( vertices.begin(), faceId.size(), faceId.begin() );
112 static FaceType makeFace (
const std::vector< unsigned int > &vertices )
114 if( vertices.size() != (
dimension == 2 ? 2 : numFaceCorners) )
115 DUNE_THROW( GridError,
"Wrong number of face vertices passed: " << vertices.size() <<
"." );
121 copy( { 0, vertices[ 1 ]+1, vertices[ 0 ]+1 }, faceId );
123 copy( { 2*vertices[ 0 ], 2*vertices[ 1 ], 2*vertices[ 0 ]+1, 2*vertices[ 1 ]+1 }, faceId );
126 std::copy_n( vertices.begin(), numFaceCorners, faceId.begin() );
130 static BndPair makeBndPair (
const FaceType &face,
int id )
133 for(
unsigned int i = 0; i < numFaceCorners; ++i )
136 bndPair.first[ j ] = face[ i ];
144 virtual Grid* createGridObj( BoundaryProjectionVector* bndProjections,
const std::string& name )
const
146 return new Grid( communicator_, globalProjection_, bndProjections , name, realGrid_ );
156 bool removeGeneratedFile =
true );
169 virtual void insertVertex (
const VertexInputType &pos );
188 const std::vector< VertexId > &vertices );
202 insertBoundary (
const GeometryType &geometry,
const std::vector< VertexId > &faceVertices,
int boundaryId = 1 );
212 virtual void insertBoundary (
int element,
int face,
int boundaryId = 1 );
230 const std::vector< VertexId > &vertices,
250 const shared_ptr<BoundarySegment<dimension,dimensionworld> >& boundarySegment ) ;
275 Grid *
createGrid (
const bool addMissingBoundaries,
const std::string dgfName =
"" );
277 Grid *
createGrid (
const bool addMissingBoundaries,
bool temporary,
const std::string dgfName =
"" );
282 alugrid_assert( Grid::getRealImplementation( entity ).getIndex() <
int(ordering_.size()) );
283 return ordering_[ Grid::getRealImplementation( entity ).getIndex() ];
291 return Grid::getRealImplementation( entity ).getIndex()/2;
294 return Grid::getRealImplementation( entity ).getIndex() - 1;
296 return Grid::getRealImplementation( entity ).getIndex();
302 std::vector< unsigned int > vertices;
303 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
310 const Dune::ReferenceElement< double, dimension > &refElem =
311 Dune::ReferenceElements< double, dimension >::general( in.type() );
312 int faceNr = intersection.indexInInside();
313 const int vxSize = refElem.size( faceNr, 1,
dimension );
314 for (
int i=0;i<vxSize;++i)
316 int vxIdx = refElem.subEntity( faceNr, 1 , i ,
dimension);
317 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
318 vertices.push_back(
insertionIndex( in.template subEntity<dimension>(vxIdx) ) );
320 vertices.push_back(
insertionIndex( *(in.template subEntity<dimension>(vxIdx) ) ) );
324 FaceType faceId = makeFace( vertices );
325 std::sort( faceId.begin(), faceId.end() );
327 typename BoundaryIdMap::const_iterator pos = insertionOrder_.find( faceId );
328 if( pos != insertionOrder_.end() )
331 return std::numeric_limits<unsigned int>::max();
335 wasInserted (
const typename Grid::LeafIntersection &intersection )
const
337 return intersection.boundary() &&
338 (
insertionIndex(intersection) < std::numeric_limits<unsigned int>::max());
341 const std::vector<unsigned int>&
ordering ()
const {
return ordering_; }
344 void doInsertVertex (
const VertexInputType &pos,
const GlobalIdType globalId );
345 void doInsertBoundary (
int element,
int face,
int boundaryId );
350 return vertices_[ id ].second;
353 const VertexType &position (
const VertexId &
id )
const
356 return vertices_[ id ].first;
359 const VertexInputType inputPosition (
const VertexId &
id )
const
362 VertexType vertex = vertices_[ id ].first;
363 VertexInputType iVtx(0.);
369 void assertGeometryType(
const GeometryType &geometry );
370 static void generateFace (
const ElementType &element,
const int f, FaceType &face );
371 void generateFace (
const SubEntity &subEntity, FaceType &face )
const;
372 void correctElementOrientation ();
373 bool identifyFaces (
const Transformation &transformation,
const FaceType &key1,
const FaceType &key2,
const int defaultId );
374 void searchPeriodicNeighbor ( FaceMap &faceMap,
const typename FaceMap::iterator &pos,
const int defaultId );
375 void reinsertBoundary (
const FaceMap &faceMap,
const typename FaceMap::const_iterator &pos,
const int id );
376 void recreateBoundaryIds (
const int defaultId = 1 );
379 void sortElements(
const VertexVector& vertices,
const ElementVector& elements, std::vector< unsigned int >&
ordering );
383 VertexVector vertices_;
384 ElementVector elements_;
385 BoundaryIdMap boundaryIds_,insertionOrder_;
386 PeriodicBoundaryVector periodicBoundaries_;
388 BoundaryProjectionMap boundaryProjections_;
389 FaceTransformationVector faceTransformations_;
390 unsigned int numFacesInserted_;
392 const bool allowGridGeneration_;
393 bool foundGlobalIndex_ ;
398 std::vector< unsigned int > ordering_;
403 template<
class ALUGr
id >
405 :
public std::binary_function< FaceType, FaceType, bool >
407 bool operator() (
const FaceType &a,
const FaceType &b )
const
409 for(
unsigned int i = 0; i < numFaceCorners; ++i )
411 if( a[ i ] != b[ i ] )
412 return (a[ i ] < b[ i ]);
419 template<
class ALUGr
id >
423 if( elementType ==
tetra )
425 if( !geometry.isSimplex() )
426 DUNE_THROW( GridError,
"Only simplex geometries can be inserted into "
427 "ALUGrid< 3, 3, simplex, refrule >." );
431 if( !geometry.isCube() )
432 DUNE_THROW( GridError,
"Only cube geometries can be inserted into "
433 "ALUGrid< 3, 3, cube, refrule >." );
440 template<
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype ,
class Comm >
441 class GridFactory<
ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
442 :
public ALU3dGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
444 typedef GridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
ThisType;
460 :
BaseType( filename, communicator )
464 template<
class Gr
id >
468 template<
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype ,
class Comm >
470 :
public ALU3dGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
491 template<
class ALUGr
id >
495 bool removeGeneratedFile )
497 globalProjection_ ( 0 ),
498 numFacesInserted_ ( 0 ),
500 allowGridGeneration_( rank_ == 0 ),
501 foundGlobalIndex_( false ),
502 communicator_( communicator ),
506 template<
class ALUGr
id >
512 globalProjection_ ( 0 ),
513 numFacesInserted_ ( 0 ),
515 allowGridGeneration_( rank_ == 0 ),
516 foundGlobalIndex_( false ),
517 communicator_( communicator ),
521 template<
class ALUGr
id >
527 globalProjection_ ( 0 ),
528 numFacesInserted_ ( 0 ),
529 realGrid_( realGrid ),
530 allowGridGeneration_( true ),
531 foundGlobalIndex_( false ),
532 communicator_( communicator ),
536 template<
class ALUGr
id >
540 FaceType faceId = makeFace( vertices );
542 boundaryIds_.insert( makeBndPair( faceId, 1 ) );
544 std::sort( faceId.begin(), faceId.end() );
545 if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
546 DUNE_THROW( GridError,
"Only one boundary projection can be attached to a face." );
548 boundaryProjections_[ faceId ] =
nullptr;
550 insertionOrder_.insert( std::make_pair( faceId, insertionOrder_.size() ) );
553 template<
class ALUGr
id >
557 FaceType faceId = makeFace( vertices );
561 std::sort( faceId.begin(), faceId.end() );
562 boundaryProjections_[ faceId ] =
nullptr;
565 template<
class ALUGr
id >
568 const shared_ptr<BoundarySegment<dimension,dimensionworld> >& boundarySegment )
570 FaceType faceId = makeFace( vertices );
572 const size_t numVx = vertices.size();
574 if( elementType ==
tetra )
575 type.makeSimplex( dimension-1 );
577 type.makeCube( dimension-1 );
581 typedef FieldVector< double, dimensionworld > CoordType;
582 std::vector< CoordType > coords( numVx );
583 for(
size_t i = 0; i < numVx; ++i )
589 const VertexType &x = position( vertices[ i ] );
590 for(
unsigned int j = 0; j < dimensionworld; ++j )
591 coords[ i ][ j ] = x[ j ];
594 boundaryIds_.insert( makeBndPair( faceId, 1 ) );
596 std::sort( faceId.begin(), faceId.end() );
597 if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
598 DUNE_THROW( GridError,
"Only one boundary projection can be attached to a face." );
600 BoundarySegmentWrapperType* prj
601 =
new BoundarySegmentWrapperType( type, coords, boundarySegment );
602 boundaryProjections_[ faceId ] = prj;
605 for(
size_t i = 0; i < numVx; ++i )
607 CoordType global = (*prj)( coords [ i ] );
608 if( (global - coords[ i ]).two_norm() > 1e-6 )
609 DUNE_THROW(GridError,
"BoundarySegment does not map face vertices to face vertices.");
613 insertionOrder_.insert( std::make_pair( faceId, insertionOrder_.size() ) );
617 template<
class ALUGr
id >
619 ::generateFace (
const SubEntity &subEntity, FaceType &face )
const
621 generateFace( elements_[ subEntity.first ], subEntity.second, face );
626 #if COMPILE_ALUGRID_INLINE