albertagrid/geometry.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALBERTA_GEOMETRY_HH
#define DUNE_ALBERTA_GEOMETRY_HH

#include <dune/grid/common/geometry.hh>
#include <dune/grid/albertagrid/misc.hh>
#include <dune/grid/albertagrid/elementinfo.hh>

#if HAVE_ALBERTA

namespace Dune
{

  // Forward Declarations
  // --------------------

  template< int dim, int dimworld >
  class AlbertaGrid;



  // AlbertaGridCoordinateReader
  // ---------------------------

  template< int codim, class GridImp >
  struct AlbertaGridCoordinateReader
  {
    typedef typename std::remove_const< GridImp >::type Grid;

    static constexpr int dimension = Grid::dimension;
    static constexpr int codimension = codim;
    static constexpr int mydimension = dimension - codimension;
    static constexpr int coorddimension = Grid::dimensionworld;

    typedef Alberta::Real ctype;

    typedef Alberta::ElementInfo< dimension > ElementInfo;
    typedef FieldVector< ctype, coorddimension > Coordinate;

    AlbertaGridCoordinateReader ( const GridImp &grid,
                                  const ElementInfo &elementInfo,
                                  int subEntity )
      : grid_( grid ),
        elementInfo_( elementInfo ),
        subEntity_( subEntity )
    {}

    const ElementInfo &elementInfo () const
    {
      return elementInfo_;
    }

    void coordinate ( int i, Coordinate &x ) const
    {
      assert( !elementInfo_ == false );
      assert( (i >= 0) && (i <= mydimension) );

      const int k = mapVertices( subEntity_, i );
      const Alberta::GlobalVector &coord = grid_.getCoord( elementInfo_, k );
      for( int j = 0; j < coorddimension; ++j )
        x[ j ] = coord[ j ];
    }

    bool hasDeterminant () const
    {
      return false;
    }

    ctype determinant () const
    {
      assert( hasDeterminant() );
      return ctype( 0 );
    }

  private:
    static int mapVertices ( int subEntity, int i )
    {
      return Alberta::MapVertices< dimension, codimension >::apply( subEntity, i );
    }

    const Grid &grid_;
    const ElementInfo &elementInfo_;
    const int subEntity_;
  };



  // AlbertaGridGeometry
  // -------------------

  template< int mydim, int cdim, class GridImp >
  class AlbertaGridGeometry
  {
    typedef AlbertaGridGeometry< mydim, cdim, GridImp > This;

    // remember type of the grid
    typedef GridImp Grid;

    // dimension of barycentric coordinates
    static constexpr int dimbary = mydim + 1;

  public:
    typedef Alberta::Real ctype;

    static constexpr int dimension = Grid :: dimension;
    static constexpr int mydimension = mydim;
    static constexpr int codimension = dimension - mydimension;
    static constexpr int coorddimension = cdim;

    typedef FieldVector< ctype, mydimension > LocalCoordinate;
    typedef FieldVector< ctype, coorddimension > GlobalCoordinate;

    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
    typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;

    typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse;

  private:
    static constexpr int numCorners = mydimension + 1;

    typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;

  public:
    AlbertaGridGeometry ()
    {
      invalidate();
    }

    template< class CoordReader >
    AlbertaGridGeometry ( const CoordReader &coordReader )
    {
      build( coordReader );
    }

    GeometryType type () const
    {
      return GeometryTypes::simplex( mydimension );
    }

    bool affine () const { return true; }

    int corners () const
    {
      return numCorners;
    }

    GlobalCoordinate corner ( const int i ) const
    {
      assert( (i >= 0) && (i < corners()) );
      return coord_[ i ];
    }

    GlobalCoordinate center () const
    {
      return centroid_;
    }

    GlobalCoordinate global ( const LocalCoordinate &local ) const;

    LocalCoordinate local ( const GlobalCoordinate &global ) const;

    ctype integrationElement () const
    {
      assert( calcedDet_ );
      return elDet_;
    }

    ctype integrationElement ( [[maybe_unused]] const LocalCoordinate &local ) const
    {
      return integrationElement();
    }

    ctype volume () const
    {
      return integrationElement() / ctype( factorial(mydimension) );
    }

    const JacobianTransposed &jacobianTransposed () const;

    const JacobianTransposed &
    jacobianTransposed ( [[maybe_unused]] const LocalCoordinate &local ) const
    {
      return jacobianTransposed();
    }

    const JacobianInverseTransposed &jacobianInverseTransposed () const;

    const JacobianInverseTransposed &
    jacobianInverseTransposed ( [[maybe_unused]] const LocalCoordinate &local ) const
    {
      return jacobianInverseTransposed();
    }

    Jacobian jacobian ( const LocalCoordinate &local ) const
    {
      return jacobianTransposed(local).transposed();
    }

    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
    {
      return jacobianInverseTransposed(local).transposed();
    }

    //***********************************************************************
    //  Methods that not belong to the Interface, but have to be public
    //***********************************************************************

    void invalidate ()
    {
      builtJT_ = false;
      builtJTInv_ = false;
      calcedDet_ = false;
    }

    template< class CoordReader >
    void build ( const CoordReader &coordReader );

    void print ( std::ostream &out ) const;

  private:
    // calculates the volume of the element
    ctype elDeterminant () const
    {
      return std::abs( Alberta::determinant( jacobianTransposed() ) );
    }

    CoordMatrix coord_;

    GlobalCoordinate centroid_;

    // storage for the transposed of the jacobian
    mutable JacobianTransposed jT_;

    // storage for the transposed inverse of the jacboian
    mutable JacobianInverseTransposed jTInv_;

    // has jT_ been computed, yet?
    mutable bool builtJT_;
    // has jTInv_ been computed, yet?
    mutable bool builtJTInv_;

    mutable bool calcedDet_; 
    mutable ctype elDet_; 
  };



  // AlbertaGridGlobalGeometry
  // -------------------------

  template< int mydim, int cdim, class GridImp >
  class AlbertaGridGlobalGeometry
    : public AlbertaGridGeometry< mydim, cdim, GridImp >
  {
    typedef AlbertaGridGlobalGeometry< mydim, cdim, GridImp > This;
    typedef AlbertaGridGeometry< mydim, cdim, GridImp > Base;

  public:
    AlbertaGridGlobalGeometry ()
      : Base()
    {}

    template< class CoordReader >
    AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
      : Base( coordReader )
    {}
  };


#if !DUNE_ALBERTA_CACHE_COORDINATES
  template< int dim, int cdim >
  class AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > >
  {
    typedef AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > > This;

    // remember type of the grid
    typedef AlbertaGrid< dim, cdim > Grid;

    // dimension of barycentric coordinates
    static constexpr int dimbary = dim + 1;

    typedef Alberta::ElementInfo< dim > ElementInfo;

  public:
    typedef Alberta::Real ctype;

    static constexpr int dimension = Grid::dimension;
    static constexpr int mydimension = dimension;
    static constexpr int codimension = dimension - mydimension;
    static constexpr int coorddimension = cdim;

    typedef FieldVector< ctype, mydimension > LocalCoordinate;
    typedef FieldVector< ctype, coorddimension > GlobalCoordinate;

    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
    typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;

    typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse;

  private:
    static constexpr int numCorners = mydimension + 1;

    typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;

  public:
    AlbertaGridGlobalGeometry ()
    {
      invalidate();
    }

    template< class CoordReader >
    AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
    {
      build( coordReader );
    }

    GeometryType type () const
    {
      return GeometryTypes::simplex( mydimension );
    }

    int corners () const
    {
      return numCorners;
    }

    GlobalCoordinate corner ( const int i ) const
    {
      assert( (i >= 0) && (i < corners()) );
      const Alberta::GlobalCoordinate &x = elementInfo_.coordinate( i );
      GlobalCoordinate y;
      for( int j = 0; j < coorddimension; ++j )
        y[ j ] = x[ j ];
      return y;
    }

    GlobalCoordinate center () const
    {
      GlobalCoordinate centroid_ = corner( 0 );
      for( int i = 1; i < numCorners; ++i )
        centroid_ += corner( i );
      centroid_ *= ctype( 1 ) / ctype( numCorners );
      return centroid_;
    }

    GlobalCoordinate global ( const LocalCoordinate &local ) const;

    LocalCoordinate local ( const GlobalCoordinate &global ) const;

    ctype integrationElement () const
    {
      return elementInfo_.geometryCache().integrationElement();
    }

    ctype integrationElement ( const LocalCoordinate &local ) const
    {
      return integrationElement();
    }

    ctype volume () const
    {
      return integrationElement() / ctype( factorial(mydimension) );
    }

    const JacobianTransposed &jacobianTransposed () const
    {
      return elementInfo_.geometryCache().jacobianTransposed();
    }

    const JacobianTransposed &
    jacobianTransposed ( const LocalCoordinate &local ) const
    {
      return jacobianTransposed();
    }

    const JacobianInverseTransposed &jacobianInverseTransposed () const
    {
      return elementInfo_.geometryCache().jacobianInverseTransposed();
    }

    const JacobianInverseTransposed &
    jacobianInverseTransposed ( const LocalCoordinate &local ) const
    {
      return jacobianInverseTransposed();
    }

    Jacobian jacobian ( const LocalCoordinate &local ) const
    {
      return jacobianTransposed(local).transposed();
    }

    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
    {
      return jacobianInverseTransposed(local).transposed();
    }

    //***********************************************************************
    //  Methods that not belong to the Interface, but have to be public
    //***********************************************************************

    void invalidate ()
    {
      elementInfo_ = ElementInfo();
    }

    template< class CoordReader >
    void build ( const CoordReader &coordReader )
    {
      elementInfo_ = coordReader.elementInfo();
    }

  private:
    ElementInfo elementInfo_;
  };
#endif // #if !DUNE_ALBERTA_CACHE_COORDINATES



  // AlbertaGridLocalGeometryProvider
  // --------------------------------

  template< class Grid >
  class AlbertaGridLocalGeometryProvider
  {
    typedef AlbertaGridLocalGeometryProvider< Grid > This;

  public:
    typedef typename Grid::ctype ctype;

    static constexpr int dimension = Grid::dimension;

    template< int codim >
    struct Codim
    {
      typedef AlbertaGridGeometry< dimension-codim, dimension, Grid > LocalGeometry;
    };

    typedef typename Codim< 0 >::LocalGeometry LocalElementGeometry;
    typedef typename Codim< 1 >::LocalGeometry LocalFaceGeometry;

    static constexpr int numChildren = 2;
    static constexpr int numFaces = dimension + 1;

    static constexpr int minFaceTwist = Alberta::Twist< dimension, dimension-1 >::minTwist;
    static constexpr int maxFaceTwist = Alberta::Twist< dimension, dimension-1 >::maxTwist;
    static constexpr int numFaceTwists = maxFaceTwist - minFaceTwist + 1;

  private:
    struct GeoInFatherCoordReader;
    struct FaceCoordReader;

    AlbertaGridLocalGeometryProvider ()
    {
      buildGeometryInFather();
      buildFaceGeometry();
    }

    ~AlbertaGridLocalGeometryProvider ()
    {
      for( int child = 0; child < numChildren; ++child )
      {
        delete geometryInFather_[ child ][ 0 ];
        delete geometryInFather_[ child ][ 1 ];
      }

      for( int i = 0; i < numFaces; ++i )
      {
        for( int j = 0; j < numFaceTwists; ++j )
          delete faceGeometry_[ i ][ j ];
      }
    }

    void buildGeometryInFather();
    void buildFaceGeometry();

  public:
    const LocalElementGeometry &
    geometryInFather ( int child, const int orientation = 1 ) const
    {
      assert( (child >= 0) && (child < numChildren) );
      assert( (orientation == 1) || (orientation == -1) );
      return *geometryInFather_[ child ][ (orientation + 1) / 2 ];
    }

    const LocalFaceGeometry &
    faceGeometry ( int face, int twist = 0 ) const
    {
      assert( (face >= 0) && (face < numFaces) );
      assert( (twist >= minFaceTwist) && (twist <= maxFaceTwist) );
      return *faceGeometry_[ face ][ twist - minFaceTwist ];
    }

    static const This &instance ()
    {
      static This theInstance;
      return theInstance;
    }

  private:
    template< int codim >
    static int mapVertices ( int subEntity, int i )
    {
      return Alberta::MapVertices< dimension, codim >::apply( subEntity, i );
    }

    const LocalElementGeometry *geometryInFather_[ numChildren ][ 2 ];
    const LocalFaceGeometry *faceGeometry_[ numFaces ][ numFaceTwists ];
  };

} // namespace Dune

#endif // #if HAVE_ALBERTA

#endif // #ifndef DUNE_ALBERTA_GEOMETRY_HH