geometry.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALU3DGRIDGEOMETRY_HH
#define DUNE_ALU3DGRIDGEOMETRY_HH
 
// System includes
 
// Dune includes
#include <dune/common/power.hh>
#include <dune/grid/common/grid.hh>
 
// Local includes
#include "alu3dinclude.hh"
#include "topology.hh"
#include "mappings.hh"
#include <dune/grid/alugrid/common/objectfactory.hh>
 
namespace Dune
{
 
  // Forward declarations
  template<int cd, int dim, class GridImp>
  class ALU3dGridEntity;
  template<int cd, class GridImp >
  class ALU3dGridEntityPointer;
  template<int mydim, int coorddim, class GridImp>
  class ALU3dGridGeometry;
  template< ALU3dGridElementType, class >
  class ALU3dGrid;
  class BilinearSurfaceMapping;
  class TrilinearMapping;
 
  template< class GridImp >
  class ALU3dGridIntersectionIterator;
 
  template <int cdim>
  class MyALUGridGeometryImplementation
  {
  public:
    typedef FieldVector<alu3d_ctype, cdim> CoordinateVectorType;
 
    static const signed char invalid      = -1; // means geometry is not meaningful
    static const signed char updated      =  0; // means the point values have been set
    static const signed char buildmapping =  1; // means updated and mapping was build
 
    template <int dim, int corners, class Mapping>
    class GeometryImplBase
    {
    private:
      // prohibited due to reference counting
      GeometryImplBase( const GeometryImplBase& );
 
    protected:
      static const int corners_ = corners ;
 
      typedef FieldMatrix<alu3d_ctype, corners , cdim>  CoordinateMatrixType;
 
      template <int dummy, int dimused>
      struct CoordTypeExtractorType
      {
        typedef CoordinateMatrixType Type;
      };
 
      template <int dummy>
      struct CoordTypeExtractorType< dummy, 3 >
      {
        typedef CoordinateMatrixType* Type;
      };
 
      typedef typename CoordTypeExtractorType< 0, dim > :: Type CoordinateStorageType ;
 
      typedef Mapping MappingType;
 
      CoordinateStorageType coord_ ;
 
      MappingType map_;
 
      double volume_ ;
 
      mutable unsigned int refCount_;
 
      signed char status_ ;
    public:
      GeometryImplBase()
        : coord_( 0 ),
          map_(),
          volume_( 1.0 )
      {
        reset();
      }
 
      void reset()
      {
        // reset reference counter
        refCount_ = 1;
        // reset status
        status_   = invalid ;
      }
 
      void operator ++ () { ++ refCount_; }
 
      void operator -- () { assert( refCount_ > 0 ); --refCount_; }
 
      bool operator ! () const { return refCount_ == 0; }
 
      bool stillUsed () const { return refCount_ > 1 ; }
 
      // copy coordinate vector from field vector or alu3d_ctype[cdim]
      template <class CoordPtrType>
      static inline void copy(const CoordPtrType& p,
                              CoordinateVectorType& c)
      {
        assert( cdim == 3 );
        c[0] = p[0];
        c[1] = p[1];
        c[2] = p[2];
      }
 
      template <class CoordPtrType>
      void update(const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType& ) const
      {
        DUNE_THROW(InvalidStateException,"This method should not be called!");
      }
 
      template <class CoordPtrType>
      void update(const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType& ) const
      {
        DUNE_THROW(InvalidStateException,"This method should not be called!");
      }
 
      template <class CoordPtrType>
      void update(const CoordPtrType&,
                  const CoordPtrType&,
                  const CoordPtrType& ) const
      {
        DUNE_THROW(InvalidStateException,"This method should not be called!");
      }
 
      // set status to invalid
      void invalidate () { status_ = invalid ; }
 
      // return true if geometry is valid
      bool valid () const { return status_ != invalid ; }
 
      // set volume
      void setVolume( const double volume ) { volume_ = volume ; }
 
      // return volume
      double volume() const { return volume_; }
    };
 
    template <int dummy, int dim,
        ALU3dGridElementType eltype> class GeometryImpl;
  public:
    // geometry implementation for edges and vertices
    template <int dummy, int dim, ALU3dGridElementType eltype>
    class GeometryImpl : public GeometryImplBase< dim, dim+1, LinearMapping<cdim, dim> >
    {
      typedef GeometryImplBase< dim, dim+1, LinearMapping<cdim, dim> > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType ;
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      // return coordinate vector
      inline const CoordinateVectorType& operator [] (const int i) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        return coord_[i];
      }
 
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( coord_[0] );
        status_ = buildmapping ;
        return map_;
      }
 
      // update vertex
      template <class CoordPtrType>
      inline void update(const CoordPtrType& p0)
      {
        assert( corners_ == 1 );
        copy( p0, coord_[0] );
        // we need to update the mapping
        status_ = updated ;
      }
    };
 
    // geometry implementation for edges and vertices
    template <int dummy, ALU3dGridElementType eltype>
    class GeometryImpl<dummy,1,eltype>
      : public GeometryImplBase< 1, 2, LinearMapping<cdim, 1> >
    {
      enum { dim = 1 };
      typedef GeometryImplBase< dim, dim+1, LinearMapping<cdim, dim> > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType;
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      // return coordinate vector
      inline const CoordinateVectorType& operator [] (const int i) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        return coord_[i];
      }
 
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( coord_[0], coord_[1] );
        status_ = buildmapping ;
        return map_;
      }
 
      // update edge
      template <class CoordPtrType>
      inline void update(const CoordPtrType& p0,
                         const CoordPtrType& p1)
      {
        assert( corners_ == 2 );
        copy( p0, coord_[0] );
        copy( p1, coord_[1] );
        status_ = updated;
      }
    };
 
    // geom impl for simplex faces (triangles)
    template <int dummy>
    class GeometryImpl<dummy, 2, tetra>
      : public GeometryImplBase< 2, 3, LinearMapping<cdim, 2> >
    {
      // dim = 2, corners = 3
      typedef GeometryImplBase< 2, 3, LinearMapping<cdim, 2> > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType ;
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      // return coordinate vector
      inline const CoordinateVectorType& operator [] (const int i) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        return coord_[i];
      }
 
      // update geometry coordinates
      template <class CoordPtrType>
      inline void update(const CoordPtrType& p0,
                         const CoordPtrType& p1,
                         const CoordPtrType& p2)
      {
        copy(p0, coord_[0] );
        copy(p1, coord_[1] );
        copy(p2, coord_[2] );
        status_ = updated;
      }
 
      // return mapping (always up2date)
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( coord_[0], coord_[1], coord_[2] );
        status_ = buildmapping ;
        return map_;
      }
    };
 
    //
    //  hexa specializations
    //
 
    // geom impl for cube faces (quadrilaterals)
    template <int dummy>
    class GeometryImpl<dummy, 2, hexa>
      : public GeometryImplBase< 2, 4, BilinearSurfaceMapping >
    {
      // dim = 2, corners = 4
      typedef GeometryImplBase< 2, 4, BilinearSurfaceMapping > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType ;
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      // return coordinate vector
      inline const CoordinateVectorType& operator [] (const int i) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        return coord_[i];
      }
 
      // update geometry coordinates
      template <class CoordPtrType>
      inline void update(const CoordPtrType& p0,
                         const CoordPtrType& p1,
                         const CoordPtrType& p2,
                         const CoordPtrType& p3)
      {
        copy(p0, coord_[0] );
        copy(p1, coord_[1] );
        copy(p2, coord_[2] );
        copy(p3, coord_[3] );
        status_ = updated;
      }
 
      // return mapping (always up2date)
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( coord_[0], coord_[1], coord_[2], coord_[3] );
        status_ = buildmapping ;
        return map_;
      }
    };
 
    // geometry impl for hexahedrons
    template <int dummy>
    class GeometryImpl<dummy,3, hexa>
      : public GeometryImplBase< 3, 8, TrilinearMapping >
    {
      // dim = 3, corners = 8
      typedef GeometryImplBase< 3, 8, TrilinearMapping > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType ;
      typedef typename BaseType :: CoordinateMatrixType CoordinateMatrixType;
 
      typedef alu3d_ctype CoordPtrType[cdim];
 
      // coordinate pointer vector
      const alu3d_ctype* coordPtr_[ corners_ ];
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      GeometryImpl() : BaseType()
      {
        // set initialize coord pointers
        for( int i=0; i<corners_; ++i )
          coordPtr_[ i ] = 0;
      }
 
      // desctructor
      ~GeometryImpl()
      {
        if( coord_ ) delete coord_;
      }
 
      const alu3d_ctype* point( const int i ) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        assert( coordPtr_[i] );
        return coordPtr_[ i ];
      }
 
      // return coordinates
      inline CoordinateVectorType operator [] (const int i) const
      {
        CoordinateVectorType coord ;
        copy( point( i ), coord );
        return coord ;
      }
 
      // update geometry coordinates
      inline void update(const CoordPtrType& p0,
                         const CoordPtrType& p1,
                         const CoordPtrType& p2,
                         const CoordPtrType& p3,
                         const CoordPtrType& p4,
                         const CoordPtrType& p5,
                         const CoordPtrType& p6,
                         const CoordPtrType& p7)
      {
        coordPtr_[0] = &p0[ 0 ];
        coordPtr_[1] = &p1[ 0 ];
        coordPtr_[2] = &p2[ 0 ];
        coordPtr_[3] = &p3[ 0 ];
        coordPtr_[4] = &p4[ 0 ];
        coordPtr_[5] = &p5[ 0 ];
        coordPtr_[6] = &p6[ 0 ];
        coordPtr_[7] = &p7[ 0 ];
        status_ = updated;
      }
 
      // update geometry in father coordinates
      template <class GeometryImp>
      inline void updateInFather(const GeometryImp &fatherGeom ,
                                 const GeometryImp &myGeom)
      {
        if( coord_ == 0 )
        {
          coord_ = new CoordinateMatrixType();
        }
 
        CoordinateMatrixType& coord = *coord_;
        // compute the local coordinates in father refelem
        for(int i=0; i < myGeom.corners() ; ++i)
        {
          // calculate coordinate
          coord[i] = fatherGeom.local( myGeom.corner( i ) );
 
          // set pointer
          coordPtr_[i] = (&(coord[i][0]));
 
          // to avoid rounding errors
          for(int j=0; j<cdim; ++j)
          {
            if ( coord[i][j] < 1e-16) coord[i][j] = 0.0;
          }
        }
 
        status_ = updated ;
      }
 
      // return mapping (always up2date)
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( point( 0 ), point( 1 ), point( 2 ), point( 3 ),
                           point( 4 ), point( 5 ), point( 6 ), point( 7 ) );
 
        status_ = buildmapping;
        return map_;
      }
 
      // set status to invalid
      void invalidate () { status_ = invalid ; }
 
      // return true if geometry is valid
      bool valid () const { return status_ != invalid ; }
    };
 
 
    // geometry impl for hexahedrons
    template <int dummy>
    class GeometryImpl<dummy,3, tetra>
      : public GeometryImplBase< 3, 4, LinearMapping<cdim, cdim> >
    {
      // dim = 3, corners = 8
      typedef GeometryImplBase< 3, 4, LinearMapping<cdim, cdim> > BaseType;
 
      using BaseType :: corners_ ;
      using BaseType :: copy ;
      using BaseType :: coord_ ;
      using BaseType :: map_ ;
      using BaseType :: status_ ;
 
      typedef typename BaseType :: MappingType MappingType ;
      typedef typename BaseType :: CoordinateMatrixType CoordinateMatrixType;
 
      typedef alu3d_ctype CoordPtrType[cdim];
 
      // coordinate pointer vector
      const alu3d_ctype* coordPtr_[ corners_ ];
    public:
      using BaseType :: update ;
      using BaseType :: valid ;
 
      // default constructor
      GeometryImpl() : BaseType()
      {
        // set initialize coord pointers
        for( int i=0; i<corners_; ++i )
          coordPtr_[ i ] = 0;
      }
 
      // destructor
      ~GeometryImpl()
      {
        if( coord_ ) delete coord_;
      }
 
      const alu3d_ctype* point( const int i ) const
      {
        assert( valid() );
        assert( i>=0 && i<corners_ );
        assert( coordPtr_[ i ] );
        return coordPtr_[ i ];
      }
 
      // return coordinate vector
      inline CoordinateVectorType operator [] (const int i) const
      {
        CoordinateVectorType coord ;
        copy( point( i ), coord );
        return coord ;
      }
 
      // update geometry coordinates
      inline void update(const CoordPtrType& p0,
                         const CoordPtrType& p1,
                         const CoordPtrType& p2,
                         const CoordPtrType& p3)
      {
        coordPtr_[0] = &p0[ 0 ];
        coordPtr_[1] = &p1[ 0 ];
        coordPtr_[2] = &p2[ 0 ];
        coordPtr_[3] = &p3[ 0 ];
        status_ = updated;
      }
 
      // update geometry in father coordinates
      template <class GeometryImp>
      inline void updateInFather(const GeometryImp &fatherGeom ,
                                 const GeometryImp & myGeom)
      {
        if( coord_ == 0 )
        {
          coord_ = new CoordinateMatrixType();
        }
 
        CoordinateMatrixType& coord = *coord_;
        // compute the local coordinates in father refelem
        for(int i=0; i < myGeom.corners() ; ++i)
        {
          // calculate coordinate
          coord[i] = fatherGeom.local( myGeom.corner( i ) );
 
          // set pointer
          coordPtr_[i] = (&(coord[i][0]));
 
          // to avoid rounding errors
          for(int j=0; j<cdim; ++j)
          {
            if ( coord[i][j] < 1e-16) coord[i][j] = 0.0;
          }
        }
 
        status_ = updated;
      }
 
      // return mapping (always up2date)
      inline MappingType& mapping()
      {
        assert( valid() );
        if( status_ == buildmapping ) return map_;
 
        map_.buildMapping( point( 0 ), point( 1 ), point( 2 ), point( 3 ) );
 
        status_ = buildmapping;
        return map_;
      }
    };
  }; // end of class ALUGridGeometryImplementation
 
  template <int mydim, int cdim, class GridImp>
  class ALU3dGridGeometry :
    public GeometryDefaultImplementation<mydim, cdim, GridImp, ALU3dGridGeometry>
  {
    static const ALU3dGridElementType elementType = GridImp::elementType;
 
    typedef typename GridImp::MPICommunicatorType Comm;
 
    friend class ALU3dGridIntersectionIterator<GridImp>;
 
    typedef typename ALU3dImplTraits< elementType, Comm >::IMPLElementType IMPLElementType;
    typedef typename ALU3dImplTraits< elementType, Comm >::GEOFaceType GEOFaceType;
    typedef typename ALU3dImplTraits< elementType, Comm >::GEOEdgeType GEOEdgeType;
    typedef typename ALU3dImplTraits< elementType, Comm >::GEOVertexType GEOVertexType;
 
    // interface types
    typedef typename ALU3dImplTraits< elementType, Comm >::HFaceType HFaceType;
    typedef typename ALU3dImplTraits< elementType, Comm >::HEdgeType HEdgeType;
    typedef typename ALU3dImplTraits< elementType, Comm >::VertexType VertexType;
 
    typedef ElementTopologyMapping<elementType> ElementTopo;
    typedef FaceTopologyMapping<elementType> FaceTopo;
 
    enum { corners_      = (elementType == hexa) ? StaticPower<2,mydim>::power : mydim+1 };
 
    // type of specialized geometry implementation
    typedef typename MyALUGridGeometryImplementation<cdim> ::
    template GeometryImpl<0, mydim, elementType > GeometryImplType;
 
  public:
    typedef typename GridImp :: ctype ctype;
 
    typedef FieldVector<ctype, mydim> LocalCoordinate;
 
    typedef FieldVector<ctype, cdim > GlobalCoordinate;
 
    typedef FieldMatrix<ctype,cdim,mydim> JacobianInverseTransposed;
 
    typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
 
    // type of coordinate matrix for faces
    typedef FieldMatrix<ctype,
        EntityCount< elementType > :: numVerticesPerFace , 3> FaceCoordinatesType;
 
    ALU3dGridGeometry();
 
    ALU3dGridGeometry( const ALU3dGridGeometry& );
 
    ALU3dGridGeometry& operator = ( const ALU3dGridGeometry& );
 
    ~ALU3dGridGeometry( );
 
    GeometryType type () const;
 
    int corners () const;
 
    GlobalCoordinate corner (int i) const;
 
    GlobalCoordinate global (const LocalCoordinate& local) const;
 
    LocalCoordinate local (const GlobalCoordinate& global) const;
 
    ctype integrationElement (const LocalCoordinate& local) const;
 
    const JacobianInverseTransposed &jacobianInverseTransposed (const LocalCoordinate& local) const;
 
    const JacobianTransposed& jacobianTransposed (const LocalCoordinate& local) const;
 
    inline bool affine () const;
 
    ctype volume () const;
 
    //***********************************************************************
    //***********************************************************************
    bool buildGeom(const IMPLElementType & item);
    bool buildGeom(const HFaceType & item, int twist, int faceNum);
    bool buildGeom(const HEdgeType & item, int twist, int);
    bool buildGeom(const VertexType & item, int twist, int);
 
    // this method is used by the intersection iterator
    bool buildGeom(const FaceCoordinatesType& coords);
 
    // this method is used by the intersection iterator
    template <class coord_t>
    bool buildGeom(const coord_t& p0,
                   const coord_t& p1,
                   const coord_t& p2,
                   const coord_t& p3);
 
    // this method is used by the intersection iterator
    template <class coord_t>
    bool buildGeom(const coord_t& p0,
                   const coord_t& p1,
                   const coord_t& p2);
 
    template <class GeometryType>
    bool buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom);
 
    void print (std::ostream& ss) const;
 
    void invalidate () ;
 
    bool valid () const ;
 
  protected:
    void assign( const ALU3dGridGeometry& other );
    void removeObj();
    void getObject();
 
    // type of object provider
    typedef ALUMemoryProvider< GeometryImplType > GeometryProviderType ;
 
    static GeometryProviderType& geoProvider()
    {
#ifdef USE_SMP_PARALLEL
      typedef ALUGridObjectFactory< GridImp >  GridObjectFactoryType;
      static std::vector< GeometryProviderType > storage( GridObjectFactoryType :: maxThreads() );
      return storage[ GridObjectFactoryType :: threadNumber () ];
#else
      static GeometryProviderType storage;
      return storage;
#endif
    }
 
    // return reference to geometry implementation
    GeometryImplType& geoImpl() const
    {
      assert( geoImpl_ );
      return *geoImpl_;
    }
 
    // implementation of the coordinates and mapping
    GeometryImplType* geoImpl_;
  };
 
} // end namespace Dune
 
#include "geometry_imp.cc"
 
#endif