DiscreteFunction.hpp

#pragma once

#include <optional>
#include <vector>

#include <dune/common/ftraits.hh>
#include <dune/common/typeutilities.hh>
#include <dune/functions/common/defaultderivativetraits.hh>
#include <dune/functions/functionspacebases/flatvectorview.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/typetree/childextraction.hh>

#include <amdis/common/Tags.hpp>
#include <amdis/functions/NodeCache.hpp>
#include <amdis/linearalgebra/Access.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/typetree/RangeType.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS
{
  /// \class DiscreteFunction
  /// \brief A view on a subspace of a \ref DOFVector
  /**
  * \ingroup GridFunctions
  *
  * \tparam Coeff     Const or mutable coefficient vector
  * \tparam GB        The type of the global basis
  * \tparam TreePath  A realization of \ref Dune::TypeTree::HybridTreePath
  * \tparam Range     The range type for th evaluation of the discrete function
  *
  * **Requirements:**
  * - GB models \ref Concepts::GlobalBasis
  **/
  template <class Coeff, class GB, class TreePath, class Range = void>
  class DiscreteFunction;


  /// A mutable view on the subspace of a DOFVector, \relates DiscreteFunction
  template <class Coeff, class GB, class TreePath, class R>
  class DiscreteFunction
      : public DiscreteFunction<Coeff const,GB,TreePath,R>
  {
    using Self = DiscreteFunction<Coeff,GB,TreePath,R>;
    using Super = DiscreteFunction<Coeff const,GB,TreePath,R>;

    using Coefficients = Coeff;
    using GlobalBasis = GB;

  public:
    /// Constructor. Stores a pointer to the mutable `dofvector`.
    template <class... Path>
    DiscreteFunction(Coefficients& coefficients, GlobalBasis const& basis, Path... path)
      : Super(coefficients, basis, path...)
      , mutableCoeff_(&coefficients)
    {}

    template <class... Path>
    DiscreteFunction(Coefficients& coefficients, std::shared_ptr<GlobalBasis const> const& basis, Path... path)
      : DiscreteFunction(coefficients, *basis, path...)
    {}

    /// Construct a DiscreteFunction directly from a DOFVector
    template <class DV, class... Path,
      Dune::disableCopyMove<DiscreteFunction,DV> = 0,
      class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
      class GB_    = TYPEOF(std::declval<DV>().basis())>
    explicit DiscreteFunction(DV&& dofVector, Path... path)
      : DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
    {}

  public:
    /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
    /// reference to this DOFVector in the expression.
    /**
     * **Example:**
     * ```
     * auto v = valueOf(prob.solutionVector(),0);
     * v.interpolate_noalias([](auto const& x) { return x[0]; });
     * ```
     **/
    template <class Expr, class Tag = tag::average>
    void interpolate_noalias(Expr&& expr, Tag strategy = {});

    /// \brief Interpolation of GridFunction to DOFVector
    /**
     * **Example:**
     * ```
     * auto v = valueOf(prob.solutionVector(),0);
     * v.interpolate(v + [](auto const& x) { return x[0]; });
     * ```
     * Allows to have a reference to the DOFVector in the expression, e.g. as
     * \ref DiscreteFunction or \ref gradientOf() of a DiscreteFunction.
     **/
    template <class Expr, class Tag = tag::average>
    void interpolate(Expr&& expr, Tag strategy = {});

    /// \brief Interpolation of GridFunction to DOFVector, alias to \ref interpolate()
    template <class Expr>
    Self& operator<<(Expr&& expr)
    {
      interpolate(FWD(expr));
      return *this;
    }

    /// \brief interpolate `(*this) + expr` to DOFVector
    template <class Expr>
    Self& operator+=(Expr&& expr)
    {
      interpolate((*this) + expr);
      return *this;
    }

    /// \brief interpolate `(*this) - expr` to DOFVector
    template <class Expr>
    Self& operator-=(Expr&& expr)
    {
      interpolate((*this) - expr);
      return *this;
    }

    /// Return the mutable DOFVector
    Coefficients& coefficients()
    {
      return *mutableCoeff_;
    }

    /// Return the const DOFVector
    using Super::coefficients;

    template <class Range = void, class... Indices>
    auto child(Indices... ii)
    {
      auto tp = cat(this->treePath_, makeTreePath(ii...));
      return DiscreteFunction<Coeff, GB, TYPEOF(makeTreePath(tp)), Range>{*mutableCoeff_, this->basis(), tp};
    }

    using Super::child;

  protected:
    Coefficients* mutableCoeff_;
  };

  /// A Const DiscreteFunction
  template <class Coeff, class GB, class TreePath, class R>
  class DiscreteFunction<Coeff const,GB,TreePath,R>
  {
  private:
    using Coefficients = std::remove_const_t<Coeff>;
    using GlobalBasis = GB;
    using LocalView = typename GlobalBasis::LocalView;
    using GridView = typename GlobalBasis::GridView;

    using Coefficient = TYPEOF(std::declval<Traits::Access>()(std::declval<Coeff const>(), std::declval<typename GB::MultiIndex>()));

    using Tree = typename GlobalBasis::LocalView::Tree;
    using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;

  public:
    /// Set of entities the DiscreteFunction is defined on
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;

    /// Global coordinates of the EntitySet
    using Domain = typename EntitySet::GlobalCoordinate;

    /// Range type of this DiscreteFunction. If R=void deduce the Range automatically, using RangeType_t
    using Range = std::conditional_t<std::is_same_v<R,void>, RangeType_t<SubTree,Coefficient>, R>;

    /// \brief This GridFunction has no derivative function, it can be created
    /// by \ref DiscreteGridFunction.
    enum { hasDerivative = false };

  private:
    /// A LocalFunction representing the localfunction and derivative of the DOFVector on a bound element
    template <class Type>
    class LocalFunction;

  public:
    /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
    template <class... Path>
    DiscreteFunction(Coefficients const& coefficients, GlobalBasis const& basis, Path... path)
      : coefficients_(&coefficients)
      , basis_(&basis)
      , treePath_(makeTreePath(path...))
      , entitySet_(basis_->gridView())
    {}

    template <class... Path>
    DiscreteFunction(Coefficients const& coefficients, std::shared_ptr<GlobalBasis const> const& basisPtr, Path... path)
      : DiscreteFunction(coefficients, *basisPtr, path...)
    {}

    /// Construct a DiscreteFunction directly from a DOFVector
    template <class DV, class... Path,
      Dune::disableCopyMove<DiscreteFunction,DV> = 0,
      class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
      class GB_    = TYPEOF(std::declval<DV>().basis())>
    explicit DiscreteFunction(DV const& dofVector, Path... path)
      : DiscreteFunction(dofVector.coefficients(), dofVector.basis(), path...)
    {}

    /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
    Range operator()(Domain const& x) const;

    /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
    LocalFunction<tag::value> makeLocalFunction() const
    {
      return {std::make_shared<LocalView>(basis().localView()), treePath(), coefficients(), tag::value{}};
    }

    /// \brief Return a \ref Dune::Functions::GridViewEntitySet
    EntitySet const& entitySet() const
    {
      return entitySet_;
    }

    /// \brief Return global basis bound to the DOFVector
    GlobalBasis const& basis() const
    {
      return *basis_;
    }

    /// \brief Return treePath associated with this discrete function
    TreePath const& treePath() const
    {
      return treePath_;
    }

    /// \brief Return const coefficient vector
    Coefficients const& coefficients() const
    {
      return *coefficients_;
    }

    template <class Range = void, class... Indices>
    auto child(Indices... ii) const
    {
      auto tp = cat(this->treePath_, makeTreePath(ii...));
      return DiscreteFunction<Coeff const, GB, TYPEOF(makeTreePath(tp)), Range>{*coefficients_, *basis_, tp};
    }

  protected:
    Coefficients const* coefficients_;
    GlobalBasis const* basis_;
    TreePath treePath_;
    EntitySet entitySet_;
  };


  // deduction guides

  template <class C, class Basis, class... Indices,
    class GB = Underlying_t<Basis>,
    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...)),
    REQUIRES(Concepts::GlobalBasis<GB>)>
  DiscreteFunction(C&, Basis const&, Indices...)
    -> DiscreteFunction<C,GB,TP>;

  template <class DV, class... Indices,
    class C  = decltype(std::declval<DV>().coefficients()),
    class GB = decltype(std::declval<DV>().basis()),
    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
  DiscreteFunction(DV&, Indices...)
    -> DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP>;


  // grid functions representing the DOFVector
  // -----------------------------------------

  /// A Generator for the childs of a mutable \ref DiscreteFunction
  template <class Range = void, class C, class GB, class TP, class R, class... Indices>
  auto valueOf(DiscreteFunction<C,GB,TP,R>& df, Indices... ii)
  {
    return df.template child<Range>(ii...);
  }

  /// A Generator for the childs of a const \ref DiscreteFunction
  template <class Range = void, class C, class GB, class TP, class R, class... Indices>
  auto valueOf(DiscreteFunction<C,GB,TP,R> const& df, Indices... ii)
  {
    return df.template child<Range>(ii...);
  }

  /// A Generator to transform a DOFVector into a \ref DiscreteFunction
  template <class Range = void, class DV, class... Indices,
    class C  = decltype(std::declval<DV>().coefficients()),
    class GB = decltype(std::declval<DV>().basis()),
    class TP = TYPEOF(makeTreePath(std::declval<Indices>()...))>
  auto valueOf(DV& dofVec, Indices... ii)
  {
    using DF = DiscreteFunction<std::remove_reference_t<C>,Underlying_t<GB>,TP,Range>;
    return DF{dofVec.coefficients(), dofVec.basis(), makeTreePath(ii...)};
  }

} // end namespace AMDiS

#include "DiscreteLocalFunction.inc.hpp"
#include "DiscreteFunction.inc.hpp"