ComposerVectorGridFunction.hpp

#pragma once
 
#include <tuple>
#include <type_traits>
 
#include <amdis/Operations.hpp>
#include <amdis/algorithm/Map.hpp>
#include <amdis/common/ForEach.hpp>
#include <amdis/common/Logical.hpp>
#include <amdis/common/Order.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/Derivative.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
 
namespace AMDiS
{
#ifndef DOXYGEN
  template <class Signatur, class Element, class Functor, class LocalFunctions>
  class ComposerVectorLocalFunction;
#endif
 
  // implementation
  template <class R, class D, class E, class Functor, class LocalFct>
  class ComposerVectorLocalFunction<R(D), E, Functor, LocalFct>
  {
  public:
    using Range = R;
    using Domain = D;
    using Element = E;
 
    enum { hasDerivative = true };
 
  private:
    using LocalFctRange = TYPEOF( std::declval<LocalFct>()(std::declval<Domain>()) );
 
  public:
    ComposerVectorLocalFunction(Functor fct, std::vector<LocalFct> localFcts)
      : fct_(std::move(fct))
      , localFcts_(std::move(localFcts))
    {}
 
    void bind(Element const& element)
    {
      for (auto& lf : localFcts_)
        lf.bind(element);
    }
 
    void unbind()
    {
      for (auto& lf : localFcts_)
        lf.unbind();
    }
 
    bool bound() const
    {
      return localFcts_[0].bound();
    }
 
    Range operator()(Domain const& x) const
    {
      return fct_(Recursive::map([&](auto const& lf) { return lf(x); }, localFcts_));
    }
 
    Element const& localContext() const
    {
      return localFcts_[0].localContext();
    }
 
 
  public:
    Functor const& fct() const
    {
      return fct_;
    }
 
    auto const& localFunctions() const
    {
      return localFcts_;
    }
 
  private:
    Functor fct_;
    std::vector<LocalFct> localFcts_;
  };
 
  template <class Element, class Functor, class LocalFct>
  auto makeComposerVectorLocalFunction(Functor const& f, std::vector<LocalFct> lfs)
  {
    using D = typename Element::Geometry::LocalCoordinate;
    using LocalFctRange = TYPEOF(std::declval<LocalFct>()(std::declval<D>()));
    using R = TYPEOF(f(std::vector<LocalFctRange>{}));
    return ComposerLocalFunction<R(D), Element, Functor, LocalFct>{f, std::move(lfs)};
  }
 
  namespace Operation
  {
    struct PlusVector
    {
      template <class T>
      T operator()(std::vector<T> const& ts) const
      {
        return std::accumulate(ts.begin(),ts.end(),T(0));
      }
    };
 
    inline int order(PlusVector const&, std::vector<int> const& orders)
    {
      return *std::max_element(orders.begin(), orders.end());
    }
 
    inline auto partial(PlusVector const&, std::size_t i)
    {
      return One{};
    }
 
 
  } // end namespace Operation
 
 
  template <class Sig, class E, class F, class LF, class Type,
    REQUIRES(Concepts::HasPartial<F>)>
  auto derivativeOf(ComposerVectorLocalFunction<Sig,E,F,LF> const& composed, Type const& type)
  {
    // d_i(f)[lgfs...] * lgfs_i
    auto term_i = [&](std::size_t i)
    {
      auto di_f = makeComposerVectorLocalFunction<E>(partial(composed.fct(), i),
        composed.localFunctions());
      auto df = makeComposerLocalFunction<E>(Operation::Multiplies{}, di_f,
        derivativeOf(composed.localFunctions()[i], type));
      if (composed.bound())
        df.bind(composed.localContext());
      return df;
    };
 
    // sum_i [ d_i(f)[lgfs...] * derivativeOf(lgfs_i)
    std::vector<TYPEOF(term_i(0))> terms;
    for (std::size_t i = 0; i < composed.localFunctions().size(); ++i)
      terms.emplace_back(term_i(i));
    return makeComposerVectorLocalFunction<E>(Operation::PlusVector{}, terms);
  }
 
 
 
  template <class Sig, class E, class F, class LF,
    REQUIRES(Concepts::HasFunctorVectorOrder<F> && Concepts::Polynomial<LF>)>
  int order(ComposerVectorLocalFunction<Sig,E,F,LF> const& composed)
  {
    return order(composed.fct(),
      Recursive::map([](auto const& lf) { return order(lf); }, composed.localFunctions()));
  }
 
 
 
  template <class Sig, class EntitySet, class Functor, class GridFct>
  class ComposerVectorGridFunction;
 
  template <class R, class D, class ES, class Functor, class GridFct>
  class ComposerVectorGridFunction<R(D), ES, Functor, GridFct>
  {
  public:
    using Range = R;
    using Domain = D;
    using EntitySet = ES;
 
    enum { hasDerivative = false };
 
  private:
    using LocalFct = TYPEOF( localFunction(std::declval<GridFct const&>()) );
    using LocalDomain = typename EntitySet::LocalCoordinate;
    using Element = typename EntitySet::Element;
 
  public:
    using LocalFunction
      = ComposerVectorLocalFunction<Range(LocalDomain), Element, Functor, LocalFct>;
 
    ComposerVectorGridFunction(EntitySet const& entitySet, Functor fct, std::vector<GridFct> gridFcts)
      : entitySet_(entitySet)
      , fct_(std::move(fct))
      , gridFcts_(std::move(gridFcts))
    {}
 
    Range operator()(Domain const& x) const
    {
      return fct_(Recursive::map([&](auto const& gf) { return gf(x); }, gridFcts_));
    }
 
    EntitySet const& entitySet() const
    {
      return entitySet_;
    }
 
    LocalFunction makeLocalFunction() const
    {
      return LocalFunction{fct_,
        Recursive::map([](auto const& gf) { return localFunction(gf); }, gridFcts_)};
    }
 
  private:
    EntitySet entitySet_;
    Functor fct_;
    std::vector<GridFct> gridFcts_;
  };
 
 
  // Generator function for ComposerGridFunction expressions
  template <class Functor, class GridView, class GridFct>
  auto makeComposerVectorGridFunction(Functor const& f, GridView const& gridView,
                              std::vector<GridFct> gridFcts)
  {
    static_assert((Concepts::GridFunction<GridFct>),
      "All passed parameters must be GridFunctions.");
 
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
    using Domain = typename EntitySet::GlobalCoordinate;
    using GridFctRange = TYPEOF( std::declval<GridFct>()(std::declval<Domain>()) );
 
    static_assert(Concepts::Callable<Functor, std::vector<GridFctRange>>,
      "Range types of grid functions are not compatible with the functor.");
    using Range = TYPEOF(f(std::vector<GridFctRange>{}));
 
    using FGF = ComposerVectorGridFunction<Range(Domain), EntitySet, Functor, GridFct>;
    return FGF{EntitySet{gridView},f, std::move(gridFcts)};
  }
 
 
#ifndef DOXYGEN
  // PreGridFunction related to ComposerGridFunction.
  template <class Functor, class PreGridFct>
  struct ComposerVectorPreGridFunction
  {
    using Self = ComposerVectorPreGridFunction;
 
    struct Creator
    {
      template <class GridView>
      static auto create(Self const& self, GridView const& gridView)
      {
        return makeComposerVectorGridFunction(self.fct_, gridView,
          Recursive::map([&](auto const& pgf) { return makeGridFunction(pgf, gridView); },
            self.preGridFcts_));
      }
    };
 
    ComposerVectorPreGridFunction(Functor fct, std::vector<PreGridFct> pgfs)
      : fct_(std::move(fct))
      , preGridFcts_(std::move(pgfs))
    {}
 
  private:
    Functor fct_;
    std::vector<PreGridFct> preGridFcts_;
  };
 
  namespace Traits
  {
    template <class Functor, class PreGridFct>
    struct IsPreGridFunction<ComposerVectorPreGridFunction<Functor, PreGridFct>>
      : std::true_type {};
  }
#endif
 
 
 
  template <class Functor, class PreGridFct>
  auto invokeVectorAtQP(Functor f, std::vector<PreGridFct> gridFcts)
  {
    return ComposerVectorPreGridFunction<Functor, PreGridFct>{std::move(f), std::move(gridFcts)};
  }
 
} // end namespace AMDiS