alu3dgridfactory.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_ALU3DGRID_FACTORY_HH
#define DUNE_ALU3DGRID_FACTORY_HH
 
#include <map>
#include <vector>
 
#if HAVE_ALUGRID
 
#include <dune/common/array.hh>
#include <dune/common/parallel/mpihelper.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/common/boundaryprojection.hh>
 
#include <dune/grid/alugrid/common/transformation.hh>
#include <dune/grid/alugrid/3d/alugrid.hh>
 
namespace Dune
{
 
  template< class ALUGrid >
  class ALU3dGridFactory
    : public GridFactoryInterface< ALUGrid >
  {
    typedef ALU3dGridFactory< ALUGrid > ThisType;
    typedef GridFactoryInterface< ALUGrid > BaseType;
 
  public:
    typedef ALUGrid Grid;
 
    typedef typename Grid::ctype ctype;
 
    static const ALU3dGridElementType elementType = Grid::elementType;
 
    static const unsigned int dimension = Grid::dimension;
    static const unsigned int dimensionworld = Grid::dimensionworld;
 
    typedef typename Grid::MPICommunicatorType MPICommunicatorType;
 
    typedef DuneBoundaryProjection< 3 >  DuneBoundaryProjectionType;
 
    template< int codim >
    struct Codim
    {
      typedef typename Grid::template Codim< codim >::Entity Entity;
    };
 
    typedef unsigned int VertexId;
 
    typedef ALUGridTransformation< ctype, dimensionworld > Transformation;
 
    typedef typename Transformation::WorldVector WorldVector;
    typedef typename Transformation::WorldMatrix WorldMatrix;
 
  private:
    dune_static_assert( (elementType == tetra || elementType == hexa),
                        "ALU3dGridFactory supports only grids containing "
                        "tetrahedrons or hexahedrons exclusively." );
 
    typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapperType;
 
    static const unsigned int numCorners = EntityCount< elementType >::numVertices;
    static const unsigned int numFaces = EntityCount< elementType >::numFaces;
    static const unsigned int numFaceCorners = EntityCount< elementType >::numVerticesPerFace;
 
    typedef ElementTopologyMapping< elementType > ElementTopologyMappingType;
    typedef FaceTopologyMapping< elementType > FaceTopologyMappingType;
 
    typedef FieldVector< ctype, dimensionworld > VertexType;
    typedef std::vector< unsigned int > ElementType;
    typedef array< unsigned int, numFaceCorners > FaceType;
 
    struct FaceLess;
 
    typedef std::vector< std::pair< VertexType, size_t > > VertexVector;
    typedef std::vector< ElementType > ElementVector;
    typedef std::pair< FaceType, int > BndPair ;
    typedef std::map< FaceType,  int > BoundaryIdMap;
    typedef std::vector< std::pair< BndPair, BndPair > > PeriodicBoundaryVector;
    typedef std::pair< unsigned int, int > SubEntity;
    typedef std::map< FaceType, SubEntity, FaceLess > FaceMap;
 
    typedef std::map< FaceType, const DuneBoundaryProjectionType* > BoundaryProjectionMap;
    typedef std::vector< const DuneBoundaryProjectionType* > BoundaryProjectionVector;
 
    typedef std::vector< Transformation > FaceTransformationVector;
 
    // copy vertex numbers and store smalled #dimension ones
    void copyAndSort ( const std::vector< unsigned int > &vertices, FaceType &faceId ) const
    {
      std::vector<unsigned int> tmp( vertices );
      std::sort( tmp.begin(), tmp.end() );
 
      // copy only the first dimension vertices (enough for key)
      for( size_t i = 0; i < faceId.size(); ++i ) faceId[ i ] = tmp[ i ];
    }
 
  private:
    // return grid object
    virtual Grid* createGridObj( BoundaryProjectionVector* bndProjections, const std::string& name ) const
    {
      return ( allowGridGeneration_ ) ?
             new Grid( communicator_, globalProjection_, bndProjections , name, realGrid_ ) :
             new Grid( communicator_ );
    }
 
  public:
    explicit ALU3dGridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator(),
                                bool removeGeneratedFile = true );
 
    explicit ALU3dGridFactory ( const std::string &filename,
                                const MPICommunicatorType &communicator = Grid::defaultCommunicator() );
 
    explicit ALU3dGridFactory ( const bool verbose, const MPICommunicatorType &communicator );
 
    virtual ~ALU3dGridFactory ();
 
    virtual void insertVertex ( const VertexType &pos );
 
    // for testing parallel GridFactory
    VertexId insertVertex ( const VertexType &pos, const size_t globalId );
 
    virtual void
    insertElement ( const GeometryType &geometry,
                    const std::vector< VertexId > &vertices );
 
    virtual void
    insertBoundary ( const GeometryType &geometry,
                     const std::vector< VertexId > &faceVertices,
                     const int id );
 
    virtual void insertBoundary ( const int element, const int face, const int id );
 
    // for testing parallel GridFactory
    void insertProcessBorder ( const int element, const int face )
    {
      insertBoundary( element, face, ALU3DSPACE ProcessorBoundary_t );
    }
 
    virtual void
    insertBoundaryProjection ( const GeometryType &type,
                               const std::vector< VertexId > &vertices,
                               const DuneBoundaryProjectionType *projection );
 
    virtual void
    insertBoundarySegment ( const std::vector< VertexId >& vertices ) ;
 
    virtual void
    insertBoundarySegment ( const std::vector< VertexId >& vertices,
                            const shared_ptr<BoundarySegment<3,3> >& boundarySegment ) ;
 
    virtual void insertBoundaryProjection ( const DuneBoundaryProjectionType& bndProjection );
 
    void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
 
    Grid *createGrid ();
 
    Grid *createGrid ( const bool addMissingBoundaries, const std::string dgfName = "" );
 
    Grid *createGrid ( const bool addMissingBoundaries, bool temporary, const std::string dgfName = "" );
 
    virtual unsigned int
    insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
    {
      return Grid::getRealImplementation( entity ).getIndex();
    }
    virtual unsigned int
    insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
    {
      return Grid::getRealImplementation( entity ).getIndex();
    }
    virtual unsigned int
    insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
    {
      return intersection.boundarySegmentIndex();
    }
    virtual bool
    wasInserted ( const typename Grid::LeafIntersection &intersection ) const
    {
      return intersection.boundary() &&
             ( insertionIndex(intersection) < numFacesInserted_ );
    }
 
  private:
    size_t globalId ( const VertexId &id ) const
    {
      assert( id < vertices_.size() );
      return vertices_[ id ].second;
    }
 
    const VertexType &position ( const VertexId &id ) const
    {
      assert( id < vertices_.size() );
      return vertices_[ id ].first;
    }
 
    void assertGeometryType( const GeometryType &geometry );
    static void generateFace ( const ElementType &element, const int f, FaceType &face );
    void generateFace ( const SubEntity &subEntity, FaceType &face ) const;
    void correctElementOrientation ();
    bool identifyFaces ( const Transformation &transformation, const FaceType &key1, const FaceType &key2, const int defaultId );
    void searchPeriodicNeighbor ( FaceMap &faceMap, const typename FaceMap::iterator &pos, const int defaultId  );
    void reinsertBoundary ( const FaceMap &faceMap, const typename FaceMap::const_iterator &pos, const int id );
    void recreateBoundaryIds ( const int defaultId = 1 );
 
    int rank_;
 
    VertexVector vertices_;
    ElementVector elements_;
    BoundaryIdMap boundaryIds_;
    PeriodicBoundaryVector periodicBoundaries_;
    const DuneBoundaryProjectionType* globalProjection_ ;
    BoundaryProjectionMap boundaryProjections_;
    FaceTransformationVector faceTransformations_;
    unsigned int numFacesInserted_;
    bool realGrid_;
    const bool allowGridGeneration_;
 
    MPICommunicatorType communicator_;
  };
 
 
 
  template< class ALUGrid >
  struct ALU3dGridFactory< ALUGrid >::FaceLess
    : public std::binary_function< FaceType, FaceType, bool >
  {
    bool operator() ( const FaceType &a, const FaceType &b ) const
    {
      for( unsigned int i = 0; i < numFaceCorners; ++i )
      {
        if( a[ i ] != b[ i ] )
          return (a[ i ] < b[ i ]);
      }
      return false;
    }
  };
 
 
  template< class ALUGrid >
  inline void ALU3dGridFactory< ALUGrid >
  ::assertGeometryType( const GeometryType &geometry )
  {
    if( elementType == tetra )
    {
      if( !geometry.isSimplex() )
        DUNE_THROW( GridError, "Only simplex geometries can be inserted into "
                    "ALUSimplexGrid< 3, 3 >." );
    }
    else
    {
      if( !geometry.isCube() )
        DUNE_THROW( GridError, "Only cube geometries can be inserted into "
                    "ALUCubeGrid< 3, 3 >." );
    }
  }
 
 
  template< int dimw >
  class GridFactory< ALUSimplexGrid< 3, dimw > >
    : public ALU3dGridFactory< ALUSimplexGrid< 3, dimw > >
  {
    typedef GridFactory< ALUSimplexGrid< 3, dimw > > ThisType;
    typedef ALU3dGridFactory< ALUSimplexGrid< 3, dimw > > BaseType;
 
  public:
    typedef typename BaseType::Grid Grid;
 
    typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
 
    explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( communicator )
    {}
 
    GridFactory ( const std::string &filename,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( filename, communicator )
    {}
 
  protected:
    template< class, class, int > friend class ALULocalGeometryStorage;
    GridFactory ( const bool realGrid,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( realGrid, communicator )
    {}
  };
 
 
 
  template< int dimw >
  class GridFactory< ALUCubeGrid< 3, dimw > >
    : public ALU3dGridFactory< ALUCubeGrid< 3, dimw > >
  {
    typedef GridFactory< ALUCubeGrid< 3, dimw > > ThisType;
    typedef ALU3dGridFactory< ALUCubeGrid< 3, dimw > > BaseType;
 
  public:
    typedef typename BaseType::Grid Grid;
 
    typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
 
    explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( communicator )
    {}
 
    GridFactory ( const std::string &filename,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( filename, communicator )
    {}
 
  protected:
    template< class, class, int > friend class ALULocalGeometryStorage;
    GridFactory ( const bool realGrid,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( realGrid, communicator )
    {}
  };
 
 
  template<ALUGridElementType eltype, ALUGridRefinementType refinementtype , class Comm >
  class GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
    : public ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
  {
    typedef GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > ThisType;
    typedef ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > BaseType;
 
  public:
    typedef typename BaseType::Grid Grid;
 
    typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
 
    explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( communicator )
    {}
 
    GridFactory ( const std::string &filename,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( filename, communicator )
    {}
 
  protected:
    template< class, class, int > friend class ALULocalGeometryStorage;
    GridFactory ( const bool realGrid,
                  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
      : BaseType( realGrid, communicator )
    {}
  };
 
 
 
  // Implementation of ALU3dGridFactory
  // ----------------------------------
 
  template< class ALUGrid >
  inline
  ALU3dGridFactory< ALUGrid >
  :: ALU3dGridFactory ( const MPICommunicatorType &communicator,
                        bool removeGeneratedFile )
    : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
      globalProjection_ ( 0 ),
      numFacesInserted_ ( 0 ),
      realGrid_( true ),
      allowGridGeneration_( rank_ == 0 ),
      communicator_( communicator )
  {}
 
  template< class ALUGrid >
  inline
  ALU3dGridFactory< ALUGrid >
  :: ALU3dGridFactory ( const std::string &filename,
                        const MPICommunicatorType &communicator )
    : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
      globalProjection_ ( 0 ),
      numFacesInserted_ ( 0 ),
      realGrid_( true ),
      allowGridGeneration_( rank_ == 0 ),
      communicator_( communicator )
  {}
 
  template< class ALUGrid >
  inline
  ALU3dGridFactory< ALUGrid >
  :: ALU3dGridFactory ( const bool realGrid,
                        const MPICommunicatorType &communicator )
    : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
      globalProjection_ ( 0 ),
      numFacesInserted_ ( 0 ),
      realGrid_( realGrid ),
      allowGridGeneration_( true ),
      communicator_( communicator )
  {}
 
  template< class ALUGrid >
  inline void ALU3dGridFactory< ALUGrid > ::
  insertBoundarySegment ( const std::vector< unsigned int >& vertices )
  {
    FaceType faceId;
    copyAndSort( vertices, faceId );
 
    if( vertices.size() != numFaceCorners )
      DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
 
    if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
      DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
    //  DUNE_THROW( NotImplemented, "insertBoundarySegment with a single argument" );
 
    boundaryProjections_[ faceId ] = 0;
 
    BndPair boundaryId;
    for( unsigned int i = 0; i < numFaceCorners; ++i )
    {
      const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
      boundaryId.first[ j ] = vertices[ i ];
    }
    boundaryId.second = 1;
    boundaryIds_.insert( boundaryId );
  }
 
  template< class ALUGrid >
  inline void ALU3dGridFactory< ALUGrid > ::
  insertBoundarySegment ( const std::vector< unsigned int >& vertices,
                          const shared_ptr<BoundarySegment<3,3> >& boundarySegment )
  {
    FaceType faceId;
    copyAndSort( vertices, faceId );
 
    if( vertices.size() != numFaceCorners )
      DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
 
    if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
      DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
 
    const size_t numVx = vertices.size();
    GeometryType type;
    if( numVx == 3 )
      type.makeSimplex( dimension-1 );
    else
      type.makeCube( dimension-1 );
 
    // we need double here because of the structure of BoundarySegment
    // and BoundarySegmentWrapper which have double as coordinate type
    typedef FieldVector< double, dimensionworld > CoordType;
    std::vector< CoordType > coords( numVx );
    for( size_t i = 0; i < numVx; ++i )
    {
      // if this assertions is thrown vertices were not inserted at first
      assert( vertices_.size() > vertices[ i ] );
 
      // get global coordinate and copy it
      const VertexType &x = position( vertices[ i ] );
      for( unsigned int j = 0; j < dimensionworld; ++j )
        coords[ i ][ j ] = x[ j ];
    }
 
    BoundarySegmentWrapperType* prj
      = new BoundarySegmentWrapperType( type, coords, boundarySegment );
    boundaryProjections_[ faceId ] = prj;
#ifndef NDEBUG
    // consistency check
    for( size_t i = 0; i < numVx; ++i )
    {
      CoordType global = (*prj)( coords [ i ] );
      if( (global - coords[ i ]).two_norm() > 1e-6 )
        DUNE_THROW(GridError,"BoundarySegment does not map face vertices to face vertices.");
    }
#endif
 
    BndPair boundaryId;
    for( unsigned int i = 0; i < numFaceCorners; ++i )
    {
      const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
      boundaryId.first[ j ] = vertices[ i ];
    }
    boundaryId.second = 1;
    boundaryIds_.insert( boundaryId );
  }
 
 
  template< class ALUGrid >
  inline void ALU3dGridFactory< ALUGrid >
  ::generateFace ( const SubEntity &subEntity, FaceType &face ) const
  {
    generateFace( elements_[ subEntity.first ], subEntity.second, face );
  }
 
} // end namespace Dune
 
#endif // #if HAVE_ALUGRID
 
#if COMPILE_ALUGRID_INLINE
  #include "alu3dgridfactory.cc"
#endif
#endif