geometry.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH
#define DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH

#include <cassert>
#include <functional>
#include <iterator>
#include <limits>
#include <vector>

#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/optional.hh>
#include <dune/common/std/type_traits.hh>

#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/curvedgeometry/localfunctiongeometry.hh>

#include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>

#include "localgeometrywrapper.hh"

namespace Dune
{
  namespace CGeo
  {
    /// \brief Curved geometry implementation based on local basis function parametrization
    /**
     *  Parametrization of the geometry by a differentiable localFunction
     *
     *  \tparam  LF      localFunction parametrizing the geometry
     *  \tparam  LG      localGeometry for coordinate transform from subEntity to element,
     *                   see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface
     *  \tparam  ORDER   Polynomial order of lagrange parametrization. If order < 0 use localFunction.
     */
    template <class ct, int mydim, int cdim, class LF, class LG, int ORDER = -1>
    class Geometry
        : public LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>
    {
      using Super = LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>;

      //! type of the localFunction representation the geometry parametrization
      using LocalFunction = LF;

      //! type of coordinate transformation for subEntities to codim=0 entities
      using LocalGeometry = LG;

      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

    public:
      /// \brief constructor from localFunction to parametrize the geometry
      /**
       *  \param[in]  refElement     reference element for the geometry
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       *
       */
      template <class... Args>
      Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
        : Super(refElement, [localFunction, localGeometry=LocalGeometry{std::forward<Args>(args)...}](auto const& local)
          {
            return localFunction(localGeometry.global(local));
          })
      {}

      /// \brief constructor, forwarding to the other constructor that take a reference-element
      /**
       *  \param[in]  gt             geometry type
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       */
      template <class... Args>
      Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
        : Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...)
      {}
    };


    // Specialization for the case that ORDER < 0: Use LocalFunction `LF` directly as parametrization
    template <class ct, int mydim, int cdim, class LF, class LG>
    class Geometry<ct,mydim,cdim,LF,LG,-1>
        : public LocalFunctionGeometry<LF,LG>
    {
      using Super = LocalFunctionGeometry<LF,LG>;

      //! type of the localFunction representation the geometry parametrization
      using LocalFunction = LF;

      //! type of coordinate transformation for subEntities to codim=0 entities
      using LocalGeometry = LG;

      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

    public:
      /// \brief constructor from localFunction to parametrize the geometry
      /**
       *  \param[in]  refElement     reference element for the geometry
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       *
       */
      template <class... Args>
      Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
        : Super(refElement, localFunction, LocalGeometry{std::forward<Args>(args)...})
      {}

      /// \brief constructor, forwarding to the other constructor that take a reference-element
      /**
       *  \param[in]  gt             geometry type
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       */
      template <class... Args>
      Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
        : Super(gt, localFunction, LocalGeometry{std::forward<Args>(args)...})
      {}
    };


  #ifndef DOXYGEN

    // Specialization for vertex geometries
    template <class ct, int cdim, class LF, class LG>
    class VertexGeometry
        : public AffineGeometry<ct,0,cdim>
    {
      using Super = AffineGeometry<ct,0,cdim>;

      using LocalFunction = LF;
      using LocalGeometry = LG;

      //! tolerance to numerical algorithms
      static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }

      struct Construct {};

    public:
      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, 0>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

      /// type of local coordinates
      using LocalCoordinate = FieldVector<ct, 0>;

      /// type of global coordinates
      using GlobalCoordinate = FieldVector<ct, cdim>;

    protected:
      VertexGeometry (Construct, const ReferenceElement& refElement, const LocalFunction& localFunction,
                      const LocalGeometry& lg)
        : Super(refElement, std::vector<GlobalCoordinate>{localFunction(lg.global(refElement.position(0,0)))})
      {}

    public:
      template <class... Args>
      VertexGeometry (const ReferenceElement& refElement, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, refElement, lf, LocalGeometry(std::forward<Args>(args)...))
      {}

      template <class... Args>
      VertexGeometry (Dune::GeometryType gt, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, ReferenceElements::general(gt), lf,
                   LocalGeometry(std::forward<Args>(args)...))
      {}

      Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
      {
        auto localCoord = Super::local(globalCoord);
        if ((globalCoord - Super::global(localCoord)).two_norm2() > tolerance())
          return {};

        return localCoord;
      }

      GlobalCoordinate normal (const LocalCoordinate& local) const
      {
        DUNE_THROW(Dune::NotImplemented,
          "ERROR: normal() method only defined for edges in 2D and faces in 3D");
        return GlobalCoordinate(0);
      }

    private:
      VertexGeometry const& flatGeometry () const { return *this; }
    };


    // Specialization for vertex geometries
    template <class ct, int cdim, class LF, class LG, int order>
    class Geometry<ct,0,cdim,LF,LG,order>
        : public VertexGeometry<ct,cdim,LF,LG>
    {
      using Super = VertexGeometry<ct,cdim,LF,LG>;
    public:
      using Super::Super;
    };

    template <class ct, int cdim, class LF, class LG>
    class Geometry<ct,0,cdim,LF,LG,-1>
        : public VertexGeometry<ct,cdim,LF,LG>
    {
      using Super = VertexGeometry<ct,cdim,LF,LG>;
    public:
      using Super::Super;
    };

#endif // DOXYGEN

  } // end namespace CGeo
} // namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH