lineartransformedlocalfiniteelement.hh

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

#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LINEARTRANSFORMEDLOCALFINITEELEMENT_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LINEARTRANSFORMEDLOCALFINITEELEMENT_HH
#include <dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh>
#include <dune/istl/scaledidmatrix.hh>

namespace Dune
{
  namespace Functions::Impl
  {
    /**
     * @brief This set of classes models a global valued Finite Element similar to the classes in
     * dune/functions/functionspacebases/globalvaluedlocalfiniteelement.hh
     *
     * However, here we restrict the admissable transformations to linear transformations of the basisfunctions,
     * as is the case for various Elements like Hermite, Morley and Argyris. Additionally the Transformator classes can implemenent caching of the transformation and thus cannot be purely static anymore. This leads to some changes in the interface in constrast to range-space transformations like Piola-Transformations.
     *
     * In particular, the interfaces differ from the "GlobalValued"-classes in the following points:
     *    LinearTransformator interface:
     *      - The LinearTransformator classes implement a linear transformation from one set of basisfunction
     *        onto another, that means, the transformation is invariant to differentiation, and the need
     *        for methods like "applyJacobian" vanishes.
     *      - The LinearTransformator classes are expected to implement a method bind(Element)
     *        which sets up the transformation, for example by filling a sparse Matrix with values.
     *      - The LinearTransformator classes are not static.
     *      - The inner class LocalValuedFunction cannot implement the inverse of the transformation in the current
     *        implementation, as this would require changes to LinearTransformedInterpolation and also the inverse of
     *        the transformation (sparse Matrix inversion!). Instead, the LocalValuedFunction is the object, that,
     *        when evaluated by the LocalInterpolation class, gives the coefficients as if it was evaluated by the
     *        GlobalValued Interpolation. It can (and should) therefore implement methods used by the, for
     *        example normalDerivative(size_t) and should make use of the fact that the dune-functions interface
     *        mandates the derivative of a LocalFunction to be global-valued.
     *      - Additionally the inner class LocalValuedFunction should be differentiable as often as needed by the Node
     *        set.
     *
     *    LinearTransformedLocalBasis:
     *      - holds an instance of the Transformator class
     *      - calls transformator.prepare() whenever bound to an element
     *      - higher derivatives are implemented
     *
     *    LinearTransformedLocalInterpolation (defined only as typedef):
     *      - same as GlobalValuedInterpolation
     *
     *    LinearTransformedLocalCoefficients (class does not exist):
     *      - same as in the LocalCoefficients of the localvalued finite element
     *
     *    LinearTransformedLocalFiniteElement:
     *      - copy and paste from GlobalValuedFiniteElement with adapted typedefs, TODO maybe merge those classes
     */
    /** \brief Example of LinearTransformation Interface
     *  \tparam F Field type
     *  \tparam size number of basisfunctions
     *
     */
    template <class F, unsigned int size>
    struct IdentityTransformation
    {
      IdentityTransformation() : m(1) {}

      template <class Element>
      void bind(Element e)
      {
        // Adapt values
      }

      template <typename Values, typename LocalCoordinate, typename Geometry>
      void apply(Values &values, const LocalCoordinate &xi, const Geometry &geometry)
      {
        Values tmp = values;
        m.mv(tmp, values);
      }

      template <class Function, class LocalCoordinate, class Element>
      class LocalValuedFunction
      {
        const Function &f_;
        const Element &element_;

      public:
        LocalValuedFunction(const Function &f, const Element &e)
            : f_(f), element_(e) {}

        auto operator()(const LocalCoordinate &xi) const
        {
          auto &&f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
          return f(xi); // no transformation for functionvalues needed
        }

        friend auto derivative(LocalValuedFunction const &t)
        {
          return derivative(t.f_);
        }
      };

      ScaledIdMatrix<F, size> m;
    };

    /** \brief Implementation of a dune-localfunctions LocalBasis that applies a transformation
     *
     * \tparam Transformator The transformation that is to be applied
     * \tparam LocalValuedLocalBasis The local-valued LocalBasis that is getting transformed
     * \tparam Element The element that the global-valued FE lives on
     */
    template <class LinearTransformator, class LocalValuedLocalBasis, class Element>
    class LinearTransformedLocalBasis
    {
    public:
      using Traits = typename LocalValuedLocalBasis::Traits;

      /** \brief Bind the local basis to a particular grid element
       */
      void bind(const LocalValuedLocalBasis &localValuedLocalBasis, const Element &element)
      {
        localValuedLocalBasis_ = &localValuedLocalBasis;
        element_ = &element;
        transformator.bind(element_);
      }

      /** \brief Number of shape functions
       */
      auto size() const
      {
        return localValuedLocalBasis_->size();
      }

      //! \brief Evaluate all shape functions
      void evaluateFunction(const typename Traits::DomainType &x,
                            std::vector<typename Traits::RangeType> &out) const
      {
        localValuedLocalBasis_->evaluateFunction(x, out);

        transformator.apply(out, x, element_->geometry());
      }

      /** \brief Evaluate Jacobian of all shape functions
       *
       * \param x Point in the reference element where to evaluation the Jacobians
       * \param[out] out The Jacobians of all shape functions at the point x
       */
      void evaluateJacobian(const typename Traits::DomainType &x,
                            std::vector<typename Traits::JacobianType> &out) const
      {
        localValuedLocalBasis_->evaluateJacobian(x, out);

        transformator.apply(out, x, element_->geometry());
      }

      /** \brief Evaluate partial derivatives of any order of all shape functions
       *
       * \param order Order of the partial derivatives, in the classic multi-index notation
       * \param in Position where to evaluate the derivatives
       * \param[out] out The desired partial derivatives
       */
      void partial(const std::array<unsigned int, Traits::dimDomain> &order,
                   const typename Traits::DomainType &x,
                   std::vector<typename Traits::RangeType> &out) const
      {
        localValuedLocalBasis_->partial(x, out);
        transformator.apply(out, x, element_->geometry());
      }

      //! \brief Polynomial order of the shape functions
      auto order() const
      {
        return localValuedLocalBasis_->order();
      }

      const LocalValuedLocalBasis *localValuedLocalBasis_;
      const Element *element_;
      LinearTransformator transformator;
    };

    /** \brief Implementation of a dune-localfunctions LocalInterpolation
     *    that accepts global-valued functions
     *
     * \tparam Transformator The transformation (e.g., Piola) that transforms from local to global values
     * \tparam LocalValuedLocalInterpolation The local-valued LocalInterpolation
     *    that is used for the actual interpolation
     * \tparam Element The element that the global-valued FE lives on
     */

    template <class LinearTransformator, class LocalValuedLocalInterpolation, class Element>
    using LinearTransformedLocalInterpolatioin = GlobalValuedLocalInterpolation<LinearTransformator, LocalValuedLocalInterpolation, Element>;

    /** \brief LocalFiniteElement implementation that uses values defined wrt particular grid elements
     *
     * \tparam Transformator Class implementing range-space transformations (like the Piola transforms)
     * \tparam LocalValuedLFE LocalFiniteElement implementation whose values are to be transformed
     * \tparam Element Element where to transform the FE values to
     */
    template <class LinearTransformator, class LocalValuedLFE, class Element>
    class LinearTransformedLocalFiniteElement
    {
      using LocalBasis = LinearTransformedLocalBasis<
          Transformator, typename LocalValuedLFE::Traits::LocalBasisType, Element>;
      using LocalInterpolation = LinearTransformedLocalInterpolation<
          Transformator, typename LocalValuedLFE::Traits::LocalInterpolationType, Element>;

    public:
      /** \brief Export number types, dimensions, etc.
       */
      using Traits = LocalFiniteElementTraits<
          LocalBasis, typename LocalValuedLFE::Traits::LocalCoefficientsType, LocalInterpolation>;

      LinearTransformedLocalFiniteElement() {}

      void bind(const LocalValuedLFE &localValuedLFE, const Element &element)
      {
        globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
        globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
        localValuedLFE_ = &localValuedLFE;
      }

      /** \brief Returns the local basis, i.e., the set of shape functions
       */
      const typename Traits::LocalBasisType &localBasis() const
      {
        return globalValuedLocalBasis_;
      }

      /** \brief Returns the assignment of the degrees of freedom to the element subentities
       */
      const typename Traits::LocalCoefficientsType &localCoefficients() const
      {
        return localValuedLFE_->localCoefficients();
      }

      /** \brief Returns object that evaluates degrees of freedom
       */
      const typename Traits::LocalInterpolationType &localInterpolation() const
      {
        return globalValuedLocalInterpolation_;
      }

      /** \brief The number of shape functions */
      std::size_t size() const
      {
        return localValuedLFE_->size();
      }

      /** \brief The reference element that the local finite element is defined on
       */
      GeometryType type() const
      {
        return localValuedLFE_->type();
      }

    private:
      typename Traits::LocalBasisType globalValuedLocalBasis_;
      typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
      const LocalValuedLFE *localValuedLFE_;
    };
  } // namespace Functions::Impl
} // namespace Dune
#endif