common/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_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 implementations 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 implementations 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