Add ComposerVectorGridFunction a composer of a vector of gridfunctions

amdis/gridfunctions/CMakeLists.txt

@@ -4,6 +4,7 @@ install(FILES
     AnalyticGridFunction.hpp
     CoarsenedGridFunction.hpp
     ComposerGridFunction.hpp
+    ComposerVectorGridFunction.hpp
     ConstantGridFunction.hpp
     CoordsGridFunction.hpp
     Derivative.hpp
@@ -15,3 +16,5 @@ install(FILES
     GridFunction.hpp
     Operations.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/gridfunctions)
+
+add_subdirectory(test)

amdis/gridfunctions/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:
+    /// Constructor. Stores copies of the functor and localFunction(gridfunction)s.
+    ComposerVectorLocalFunction(Functor fct, std::vector<LocalFct> localFcts)
+      : fct_(std::move(fct))
+      , localFcts_(std::move(localFcts))
+    {}
+
+    /// Calls \ref bind for all localFunctions
+    void bind(Element const& element)
+    {
+      for (auto& lf : localFcts_)
+        lf.bind(element);
+    }
+
+    /// Calls \ref unbind for all localFunctions
+    void unbind()
+    {
+      for (auto& lf : localFcts_)
+        lf.unbind();
+    }
+
+    /// \brief Check whether the LocalFunction is bound to an element
+    bool bound() const
+    {
+      return localFcts_[0].bound();
+    }
+
+    /// Applies the functor \ref fct_ to the evaluated localFunctions
+    Range operator()(Domain const& x) const
+    {
+      return fct_(Recursive::map([&](auto const& lf) { return lf(x); }, localFcts_));
+    }
+
+    /// Get the element this localfunction (and all inner localfunctions) are bound to
+    Element const& localContext() const
+    {
+      return localFcts_[0].localContext();
+    }
+
+
+  public:
+    /// Return the stored functor
+    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
+  {
+    /** \addtogroup operations
+     *  @{
+     **/
+
+    /// Functor that represents A+B
+    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)
+    {
+      assert(i < 2 && "Derivatives of `PlusVector` only defined for the binary case.");
+      return One{};
+    }
+
+    /// @}
+
+  } // end namespace Operation
+
+  /// \fn derivative
+  /// \brief Derivative of the LocalFunction of a ComposerGridFunction, utilizing
+  /// the chain-rule. Only available if the functor provides partial derivatives.
+  /**
+   * \f$ d_x(f(lf(x)...)) = \sum_i d_i(f)[lf(x)...] * derivativeOf(lf[i]) \f$
+   *
+   * **Requirements:**
+   * - The Functor `F` must model `Concepts::HasPartial`
+   **/
+  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);
+  }
+
+
+  /// \fn order
+  /// \brief Calculate the polynomial order of a functor `F` that provides a free
+  /// function order(f, [degs...]), where degs are the orders of the LocalFunctions.
+  /**
+   * **Requirements:**
+   * - The functor `F` must model `Concepts::HasFunctorVectorOrder`
+   * - All localFunctions `LF` must model `Concepts::Polynomial`
+   **/
+  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()));
+  }
+
+
+  /// \class ComposerGridFunction
+  /// \brief A Gridfunction that applies a functor to the evaluated Gridfunctions
+  /**
+   * \ingroup GridFunctions
+   * Composition of GridFunctions `g_i` by applying a functor `f` locally, i.e.
+   * locally it is evaluated
+   * \f$ f([g_0(x), g_1(x), ...]) \f$
+   *
+   * \tparam Sig        Signature of the GridFunction
+   * \tparam EntitySet  The EntitySet this GridFunction is defined on
+   * \tparam Functor    The type of the outer functor `f`
+   * \tparam GridFct    Type of the GridFunction `g_i`
+   *
+   * Requirements:
+   * - `arity(f) == 1`
+   * - `arity(g_i) == arity(g_j) for i != j`
+   * - `g_i` models concept \ref GridFunction
+   **/
+  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>;
+
+    /// \brief Constructor. Stores copies of the functor and gridfunctions.
+    ComposerVectorGridFunction(EntitySet const& entitySet, Functor fct, std::vector<GridFct> gridFcts)
+      : entitySet_(entitySet)
+      , fct_(std::move(fct))
+      , gridFcts_(std::move(gridFcts))
+    {}
+
+    /// Applies the functor to the evaluated gridfunctions
+    Range operator()(Domain const& x) const
+    {
+      return fct_(Recursive::map([&](auto const& gf) { return gf(x); }, gridFcts_));
+    }
+
+    /// Return the stored \ref EntitySet of the first GridFunction
+    EntitySet const& entitySet() const
+    {
+      return entitySet_;
+    }
+
+    /// Create the localFunction by composition of the inner localFunctions
+    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
+
+
+  /// \brief Generator function for ComposerVectorPreGridFunction.
+  /// \relates ComposerVectorPreGridFunction
+  /**
+   * \ingroup GridFunctions
+   * Applies the functor `f` to the grid-functions `[gridFcts...]`.
+   * See \ref ComposerVectorPreGridFunction.
+   *
+   * **Examples:**
+   * - `invokeAtQP(Operation::PlusVector{}, std::vector{X(0), X(1)});`
+   **/
+  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

amdis/gridfunctions/test/CMakeLists.txt

+# additional compiler option for CMAKE_BUILD_TYPE=RelWithDebInfo only for gcc
+set(NoVarTrackingAssignments $<$<AND:$<CONFIG:RelWithDebInfo>,$<CXX_COMPILER_ID:GNU>>:-fno-var-tracking-assignments>)
+
+dune_add_test(SOURCES CoordGridFunction.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES ComposerVectorGridFunction.cpp
+  LINK_LIBRARIES amdis)

amdis/gridfunctions/test/ComposerVectorGridFunction.cpp

+#include <cassert>
+#include <limits>
+#include <vector>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/gridfunctions/ComposerVectorGridFunction.hpp>
+#include <amdis/gridfunctions/CoordsGridFunction.hpp>
+
+#include <dune/grid/yaspgrid.hh>
+
+template <class T>
+bool almost_equal(T a, T b, int ulp = 2)
+{
+  using std::abs;
+  using F = typename Dune::FieldTraits<T>::field_type;
+  return abs(a-b) <= abs(a+b) * std::numeric_limits<F>::epsilon() * ulp ||
+         abs(a-b) < std::numeric_limits<F>::min();
+}
+
+using namespace AMDiS;
+int main(int argc, char** argv)
+{
+  Environment env{argc, argv};
+
+  Dune::YaspGrid<2> grid({1.0,1.0},{4u,4u});
+  auto gridView = grid.leafGridView();
+
+  auto pgf0 = X(0);
+  auto pgf1 = X(1);
+  std::vector pgfs = {pgf0, pgf1};
+
+  auto f = [](std::vector<double> x) { return x[0]+x[1]; };
+
+  auto vec_pgf = invokeVectorAtQP(f, pgfs);
+  static_assert(Concepts::AnyGridFunction<TYPEOF(vec_pgf)>);
+
+  auto vec_pgf_gf = makeGridFunction(vec_pgf, gridView);
+  static_assert(Concepts::GridFunction<TYPEOF(vec_pgf_gf)>);
+
+  auto gf0 = makeGridFunction(pgf0, gridView);
+  auto gf1 = makeGridFunction(pgf1, gridView);
+  auto vec_gf = makeComposerVectorGridFunction(f, gridView, std::vector{gf0,gf1});
+  static_assert(Concepts::GridFunction<TYPEOF(vec_gf)>);
+
+  FieldVector<double,2> x = {0.5, 0.5};
+  assert(almost_equal(gf0(x)+gf1(x), vec_gf(x)));
+
+  auto lf0 = localFunction(gf0);
+  auto lf1 = localFunction(gf1);
+  auto vec_lf = localFunction(vec_gf);
+  for (auto const& e : elements(gridView))
+  {
+    lf0.bind(e);
+    lf1.bind(e);
+    vec_lf.bind(e);
+
+    auto refElem = referenceElement(e);
+    auto local = refElem.position(0,0);
+    assert(almost_equal(lf0(local)+lf1(local), vec_lf(local)));
+
+    lf0.unbind();
+    lf1.unbind();
+    vec_lf.unbind();
+  }
+}

amdis/gridfunctions/test/CoordGridFunction.cpp

+#include <cassert>
+#include <limits>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/gridfunctions/CoordsGridFunction.hpp>
+
+#include <dune/grid/yaspgrid.hh>
+
+template <class T>
+bool almost_equal(T a, T b, int ulp = 2)
+{
+  using std::abs;
+  using F = typename Dune::FieldTraits<T>::field_type;
+  return abs(a-b) <= abs(a+b) * std::numeric_limits<F>::epsilon() * ulp ||
+         abs(a-b) < std::numeric_limits<F>::min();
+}
+
+using namespace AMDiS;
+int main(int argc, char** argv)
+{
+  Environment env{argc, argv};
+
+  Dune::YaspGrid<2> grid({1.0,1.0},{4u,4u});
+  auto gridView = grid.leafGridView();
+
+  auto pgf0 = X(0);
+  auto pgf1 = X(1);
+  static_assert(Concepts::AnyGridFunction<TYPEOF(pgf0)>);
+  static_assert(Concepts::AnyGridFunction<TYPEOF(pgf1)>);
+  static_assert(std::is_same_v<TYPEOF(pgf0),TYPEOF(pgf1)>);
+
+  auto gf0 = makeGridFunction(pgf0, gridView);
+  auto gf1 = makeGridFunction(pgf1, gridView);
+
+  FieldVector<double,2> x = {0.2, 0.6};
+  assert(almost_equal(gf0(x), x[0]));
+  assert(almost_equal(gf1(x), x[1]));
+
+  auto lf0 = localFunction(gf0);
+  auto lf1 = localFunction(gf1);
+  for (auto const& e : elements(gridView))
+  {
+    lf0.bind(e);
+    lf1.bind(e);
+
+    auto refElem = referenceElement(e);
+    auto geo = e.geometry();
+    auto local = refElem.position(0,0);
+    auto global = geo.global(local);
+    assert(almost_equal(lf0(local), global[0]));
+    assert(almost_equal(lf1(local), global[1]));
+
+    lf0.unbind();
+    lf1.unbind();
+  }
+}

amdis/operations/Basic.hpp

@@ -17,11 +17,20 @@ namespace AMDiS
         template <class F, int... I>
         auto require(F&& f, std::integer_sequence<int,I...>) -> decltype( order(f, I...) );
       };
+
+      struct HasFunctorVectorOrder
+      {
+        template <class F>
+        auto require(F&& f, std::vector<int> o = {}) -> decltype( order(f, o) );
+      };
     }
 
     template <class F, int N>
     constexpr bool HasFunctorOrder = models<Definition::HasFunctorOrder(F, std::make_integer_sequence<int,N>)>;
 
+    template <class F>
+    constexpr bool HasFunctorVectorOrder = models<Definition::HasFunctorVectorOrder(F)>;
+
     template <class F, int N>
     using HasFunctorOrder_t = bool_t<HasFunctorOrder<F,N>>;
   }