DiscreteFunction.hpp

#pragma once

#include <optional>
#include <vector>

#include <dune/functions/common/defaultderivativetraits.hh>
#include <dune/functions/functionspacebases/defaultnodetorangemap.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/typetree/FiniteElementType.hpp>
#include <amdis/typetree/RangeType.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS
{

  template <class Coeff, class GB, class TreePath>
  class DiscreteFunction;


  template <class Coeff, class GB, class TreePath>
  class DiscreteFunction
      : public DiscreteFunction<Coeff const,GB,TreePath>
  {
    using Self = DiscreteFunction<Coeff,GB,TreePath>;
    using Super = DiscreteFunction<Coeff const,GB,TreePath>;

    using Coefficients = Coeff;
    using GlobalBasis = GB;
    using ValueType = typename Coeff::value_type;

  public:
    template <class... Path>
    DiscreteFunction(Coefficients& dofVector, GlobalBasis const& basis, Path... path)
      : Super(dofVector, basis, path...)
      , mutableCoeff_(&dofVector)
    {}

    template <class DV, class... Path,
      class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
      class GB_    = TYPEOF(*std::declval<DV>().basis())>
    DiscreteFunction(DV&& dofVector, Path... path)
      : DiscreteFunction(dofVector.coefficients(), *dofVector.basis(), path...)
    {}

  public:

    template <class Expr, class Tag = tag::average>
    void interpolate_noalias(Expr&& expr, Tag strategy = {});


    template <class Expr, class Tag = tag::average>
    void interpolate(Expr&& expr, Tag strategy = {});

    template <class Expr>
    Self& operator<<(Expr&& expr)
    {
      interpolate(FWD(expr));
      return *this;
    }

    template <class Expr>
    Self& operator+=(Expr&& expr)
    {
      interpolate((*this) + expr);
      return *this;
    }

    template <class Expr>
    Self& operator-=(Expr&& expr)
    {
      interpolate((*this) - expr);
      return *this;
    }

    Coefficients& coefficients()
    {
      return *mutableCoeff_;
    }

    using Super::coefficients;

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

    using Super::child;

  protected:
    Coefficients* mutableCoeff_;
  };

  template <class Coeff, class GB, class TreePath>
  class DiscreteFunction<Coeff const,GB,TreePath>
  {
  private:
    using Coefficients = std::remove_const_t<Coeff>;
    using GlobalBasis = GB;
    using GridView = typename GlobalBasis::GridView;
    using ValueType = typename Coefficients::value_type;

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

    using TreeCache = typename GlobalBasis::LocalView::TreeCache;
    using SubTreeCache = typename Dune::TypeTree::ChildForTreePath<TreeCache, TreePath>;

    using NodeToRangeEntry = Dune::Functions::DefaultNodeToRangeMap<SubTree>;

  public:
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;

    using Domain = typename EntitySet::GlobalCoordinate;

    using Range = RangeType_t<SubTree,ValueType>;

    enum { hasDerivative = false };

  public:
    template <class Type>
    class DerivativeLocalFunctionBase;

    class GradientLocalFunction;
    class PartialLocalFunction;
    class DivergenceLocalFunction;

    class LocalFunction;

  public:
    template <class... Path>
    DiscreteFunction(Coefficients const& coefficients, GlobalBasis const& basis, Path... path)
      : coefficients_(&coefficients)
      , basis_(&basis)
      , treePath_(makeTreePath(path...))
      , entitySet_(basis_->gridView())
      , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*basis_, treePath_))
    {}

    template <class DV, class... Path,
      class Coeff_ = TYPEOF(std::declval<DV>().coefficients()),
      class GB_    = TYPEOF(*std::declval<DV>().basis())>
    DiscreteFunction(DV const& dofVector, Path... path)
      : DiscreteFunction(dofVector.coefficients(), *dofVector.basis(), path...)
    {}


    Range operator()(Domain const& x) const;

    LocalFunction makeLocalFunction() const
    {
      return LocalFunction{*this};
    }

    EntitySet const& entitySet() const
    {
      return entitySet_;
    }

    GlobalBasis const& basis() const
    {
      return *basis_;
    }

    TreePath const& treePath() const
    {
      return treePath_;
    }

    Coefficients const& coefficients() const
    {
      return *coefficients_;
    }

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

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

  // deduction guides

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

  template <class DV, class... Path,
    class Coeff = decltype(std::declval<DV>().coefficients()),
    class GB    = decltype(*std::declval<DV>().basis()),
    class TP    = TYPEOF(makeTreePath(std::declval<Path>()...))>
  DiscreteFunction(DV&, Path...)
    -> DiscreteFunction<std::remove_reference_t<Coeff>,std::decay_t<GB>,TP>;


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

  template <class Coeff, class GB, class... Path, class... Indices>
  auto valueOf(DiscreteFunction<Coeff,GB,Path...>& df, Indices... ii)
  {
    return df.child(ii...);
  }

  template <class Coeff, class GB, class... Path, class... Indices>
  auto valueOf(DiscreteFunction<Coeff,GB,Path...> const& df, Indices... ii)
  {
    return df.child(ii...);
  }

  template <class DV, class... Indices,
    class = decltype(std::declval<DV>().coefficients()),
    class = decltype(std::declval<DV>().basis())>
  auto valueOf(DV& dofVec, Indices... ii)
  {
    return DiscreteFunction{dofVec.coefficients(), *dofVec.basis(), makeTreePath(ii...)};
  }

} // end namespace AMDiS

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