common/geometry.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (C) 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_GRID_COMMON_GEOMETRY_HH
#define DUNE_GRID_COMMON_GEOMETRY_HH
 
#include <cassert>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/transpose.hh>
#include <dune/common/std/type_traits.hh>
 
#include <dune/geometry/referenceelements.hh>
 
namespace Dune
{
 
  // External Forward Declarations
  // -----------------------------
 
  template< int dim, int dimworld, class ct, class GridFamily >
  class GridDefaultImplementation;
 
 
 
  //*****************************************************************************
  //
  // Geometry
  // forwards the interface to the implementation
  //
  //*****************************************************************************
 
  template< int mydim, int cdim, class GridImp, template< int, int, class > class GeometryImp >
  class Geometry
  {
  public:
    typedef GeometryImp< mydim, cdim, GridImp > Implementation;
 
    Implementation &impl () { return realGeometry; }
    const Implementation &impl () const { return realGeometry; }
 
    constexpr static int mydimension = mydim;
 
    constexpr static int coorddimension = cdim;
 
    typedef typename GridImp::ctype ctype;
 
    typedef FieldVector<ctype, mydim> LocalCoordinate;
 
    typedef FieldVector< ctype, cdim > GlobalCoordinate;
 
    typedef decltype(std::declval<Implementation>().volume()) Volume;
 
    typedef typename Implementation::JacobianInverseTransposed JacobianInverseTransposed;
 
    typedef typename Implementation::JacobianTransposed JacobianTransposed;
 
  private:
 
    template<class Implementation_T>
    using JacobianInverseOfImplementation = decltype(typename Implementation_T::JacobianInverse{std::declval<Implementation_T>().jacobianInverse(std::declval<LocalCoordinate>())});
 
    using JacobianInverseDefault = decltype(transpose(std::declval<JacobianInverseTransposed>()));
 
    template<class Implementation_T>
    using JacobianOfImplementation = decltype(typename Implementation_T::Jacobian{std::declval<Implementation_T>().jacobian(std::declval<LocalCoordinate>())});
 
    using JacobianDefault = decltype(transpose(std::declval<JacobianTransposed>()));
 
    auto jacobianImpl ( const LocalCoordinate &local, std::true_type /*implDetected*/ ) const
    {
      return impl().jacobian(local);
    }
 
    [[deprecated("Geometry implementatons are required to provide a jacobian(local) method. The default implementation is deprecated and will be removed after release 2.9")]]
    auto jacobianImpl ( const LocalCoordinate &local, std::false_type /*implNotDetected*/ ) const
    {
      return transpose(jacobianTransposed(local));
    }
 
    auto jacobianInverseImpl ( const LocalCoordinate &local, std::true_type /*implDetected*/ ) const
    {
      return impl().jacobianInverse(local);
    }
 
    [[deprecated("Geometry implementatons are required to provide a jacobianInverse(local) method. The default implementation is deprecated and will be removed after release 2.9")]]
    auto jacobianInverseImpl ( const LocalCoordinate &local, std::false_type /*implNotDetected*/ ) const
    {
      return transpose(jacobianInverseTransposed(local));
    }
 
  public:
 
    using JacobianInverse = Std::detected_or_t<JacobianInverseDefault, JacobianInverseOfImplementation, Implementation>;
 
    using Jacobian = Std::detected_or_t<JacobianDefault, JacobianOfImplementation, Implementation>;
 
    GeometryType type () const { return impl().type(); }
 
    bool affine() const { return impl().affine(); }
 
    int corners () const { return impl().corners(); }
 
    GlobalCoordinate corner ( int i ) const
    {
      return impl().corner( i );
    }
 
    GlobalCoordinate global (const LocalCoordinate& local) const
    {
      return impl().global( local );
    }
 
    LocalCoordinate local (const GlobalCoordinate& global) const
    {
      return impl().local( global );
    }
 
    Volume integrationElement (const LocalCoordinate& local) const
    {
      return impl().integrationElement( local );
    }
 
    Volume volume () const
    {
      return impl().volume();
    }
 
    GlobalCoordinate center () const
    {
      return impl().center();
    }
 
    JacobianTransposed jacobianTransposed ( const LocalCoordinate& local ) const
    {
      return impl().jacobianTransposed( local );
    }
 
    JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate &local ) const
    {
      return impl().jacobianInverseTransposed(local);
    }
 
    Jacobian jacobian ( const LocalCoordinate& local ) const
    {
      auto implDetected = Std::is_detected<JacobianOfImplementation, Implementation>{};
      return jacobianImpl(local, implDetected);
    }
 
    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
    {
      auto implDetected = Std::is_detected<JacobianInverseOfImplementation, Implementation>{};
      return jacobianInverseImpl(local, implDetected);
    }
 
    //===========================================================
    //===========================================================
 
    explicit Geometry ( const Implementation &impl )
      : realGeometry( impl )
    {}
 
 
  protected:
 
    Implementation realGeometry;
  };
 
 
 
  //************************************************************************
  // GEOMETRY Default Implementations
  //*************************************************************************
  //
  // --GeometryDefault
  //
  template<int mydim, int cdim, class GridImp, template<int,int,class> class GeometryImp>
  class GeometryDefaultImplementation
  {
  public:
    static const int mydimension = mydim;
    static const int coorddimension = cdim;
 
    // save typing
    typedef typename GridImp::ctype ctype;
 
    typedef FieldVector< ctype, mydim > LocalCoordinate;
    typedef FieldVector< ctype, cdim > GlobalCoordinate;
 
    typedef ctype Volume;
 
    typedef FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed;
 
    typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
 
    typedef FieldMatrix< ctype, mydim, cdim > JacobianInverse;
 
    typedef FieldMatrix< ctype, cdim, mydim > Jacobian;
 
    Volume volume () const
    {
      GeometryType type = asImp().type();
 
      // get corresponding reference element
      auto refElement = referenceElement< ctype , mydim >(type);
 
      LocalCoordinate localBaryCenter ( 0 );
      // calculate local bary center
      const int corners = refElement.size(0,0,mydim);
      for(int i=0; i<corners; ++i) localBaryCenter += refElement.position(i,mydim);
      localBaryCenter *= (ctype) (1.0/corners);
 
      // volume is volume of reference element times integrationElement
      return refElement.volume() * asImp().integrationElement(localBaryCenter);
    }
 
    GlobalCoordinate center () const
    {
      GeometryType type = asImp().type();
 
      // get corresponding reference element
      auto refElement = referenceElement< ctype , mydim >(type);
 
      // center is (for now) the centroid of the reference element mapped to
      // this geometry.
      return asImp().global(refElement.position(0,0));
    }
 
    Jacobian jacobian ( const LocalCoordinate& local ) const
    {
      return asImp().jacobianTransposed(local).transposed();
    }
 
    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
    {
      return asImp().jacobianInverseTransposed(local).transposed();
    }
 
  private:
    GeometryImp<mydim,cdim,GridImp>& asImp () {return static_cast<GeometryImp<mydim,cdim,GridImp>&>(*this);}
    const GeometryImp<mydim,cdim,GridImp>& asImp () const {return static_cast<const GeometryImp<mydim,cdim,GridImp>&>(*this);}
  }; // end GeometryDefault
 
  template<int cdim, class GridImp, template<int,int,class> class GeometryImp>
  class GeometryDefaultImplementation<0,cdim,GridImp,GeometryImp>
  {
    // my dimension
    constexpr static int mydim = 0;
 
  public:
    static const int mydimension = mydim;
    static const int coorddimension = cdim;
 
    // save typing
    typedef typename GridImp::ctype ctype;
 
    typedef FieldVector< ctype, mydim > LocalCoordinate;
    typedef FieldVector< ctype, cdim > GlobalCoordinate;
    typedef ctype Volume;
 
    typedef FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed;
 
    typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
 
    typedef FieldMatrix< ctype, mydim, cdim > JacobianInverse;
 
    typedef FieldMatrix< ctype, cdim, mydim > Jacobian;
 
    FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const
    {
      return asImp().corner(0);
    }
 
    FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& ) const
    {
      return FieldVector<ctype, mydim>();
    }
 
    Volume volume () const
    {
      return Volume(1.0);
    }
 
    FieldVector<ctype, cdim> center () const
    {
      return asImp().corner(0);
    }
 
    Jacobian jacobian ( const LocalCoordinate& local ) const
    {
      return asImp().jacobianTransposed(local).transposed();
    }
 
    JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
    {
      return asImp().jacobianInverseTransposed(local).transposed();
    }
 
  private:
    // Barton-Nackman trick
    GeometryImp<mydim,cdim,GridImp>& asImp () {return static_cast<GeometryImp<mydim,cdim,GridImp>&>(*this);}
    const GeometryImp<mydim,cdim,GridImp>& asImp () const {return static_cast<const GeometryImp<mydim,cdim,GridImp>&>(*this);}
  }; // end GeometryDefault
 
 
 
  // referenceElement
  // ----------------
 
  template< int mydim, int cdim, class GridImp, template< int, int, class > class GeometryImp>
  auto referenceElement(const Geometry<mydim,cdim,GridImp,GeometryImp>& geo)
    -> decltype(referenceElement(geo,geo.impl()))
  {
    return referenceElement(geo,geo.impl());
  }
 
 
  template< int mydim, int cdim, class GridImp, template< int, int, class > class GeometryImp, typename Impl>
  auto referenceElement(const Geometry<mydim,cdim,GridImp,GeometryImp>& geo, const Impl&)
    -> decltype(referenceElement<typename GridImp::ctype,mydim>(geo.type()))
  {
    using Geo = Geometry<mydim,cdim,GridImp,GeometryImp>;
    return referenceElement<typename Geo::ctype,Geo::mydimension>(geo.type());
  }
 
} // namespace Dune
 
#endif // DUNE_GRID_COMMON_GEOMETRY_HH