From 4a7bfa1c580294dc789953f128c522b841becaa7 Mon Sep 17 00:00:00 2001
From: Felix Schindler <felix.schindler@wwu.de>
Date: Sat, 14 Mar 2020 21:26:12 +0100
Subject: [PATCH] [local.integrands] some major updates | WILL BREAK YOUR CODE

* unary intersection integrands are applied to one test basis on one side
  of an intersection, meant for functionals
* binary intersection integrands are applied to one test and one ansatz
  basis on only one side of the intersection, meant for bilinear forms
* integrands can derive from multiple interfaces, we need to provide an
  additional combined interface for each combination for operator +
* one can now use with_ansatz(u) to obtain integrands depending only on
  the test basis
---
 dune/gdt/local/integrands/abs.hh           |   2 +-
 dune/gdt/local/integrands/combined.hh      | 519 +++++++++++----------
 dune/gdt/local/integrands/conversion.hh    | 132 ++++--
 dune/gdt/local/integrands/div.hh           |   6 +-
 dune/gdt/local/integrands/elliptic-ipdg.hh |  18 +-
 dune/gdt/local/integrands/elliptic.hh      |   2 +-
 dune/gdt/local/integrands/generic.hh       |  64 ++-
 dune/gdt/local/integrands/identity.hh      |   2 +-
 dune/gdt/local/integrands/interfaces.hh    | 362 ++++++++++++--
 dune/gdt/local/integrands/ipdg.hh          |  56 +--
 dune/gdt/local/integrands/laplace-ipdg.hh  | 129 +++--
 dune/gdt/local/integrands/laplace.hh       |   2 +-
 dune/gdt/local/integrands/product.hh       | 126 ++++-
 13 files changed, 1006 insertions(+), 414 deletions(-)

diff --git a/dune/gdt/local/integrands/abs.hh b/dune/gdt/local/integrands/abs.hh
index 970c2c3ae..c35c9e2c7 100644
--- a/dune/gdt/local/integrands/abs.hh
+++ b/dune/gdt/local/integrands/abs.hh
@@ -30,7 +30,7 @@ public:
   using typename BaseType::DomainType;
   using typename BaseType::LocalBasisType;
 
-  std::unique_ptr<BaseType> copy() const
+  std::unique_ptr<BaseType> copy_as_unary_element_integrand() const
   {
     return std::make_unique<ThisType>();
   }
diff --git a/dune/gdt/local/integrands/combined.hh b/dune/gdt/local/integrands/combined.hh
index 3c48f068b..5a73a9afd 100644
--- a/dune/gdt/local/integrands/combined.hh
+++ b/dune/gdt/local/integrands/combined.hh
@@ -18,244 +18,277 @@ namespace Dune {
 namespace GDT {
 
 
-/// \todo add LocalUnaryElementIntegrandSum
-/// \todo add operator+ to LocalUnaryElementIntegrandInterface
-/// \sa LocalQuaternaryIntersectionIntegrandSum
-// template <class E, size_t r, size_t rC, class R, class F>
-// class LocalUnaryElementIntegrandSum
-//  : public XT::Common::ParametricInterface
-//  , public XT::Grid::ElementBoundObject<Element>
-//{
-//  static_assert(XT::Grid::is_entity<Element>::value, "");
-
-//  using ThisType = LocalUnaryElementIntegrandInterface<Element, range_dim, range_dim_cols, RangeField, Field>;
-
-// public:
-//  using E = Element;
-//  using D = typename Element::Geometry::ctype;
-//  static const constexpr size_t d = E::dimension;
-//  using F = Field;
-
-//  using R = RangeField;
-//  static const constexpr size_t r = range_dim;
-//  static const constexpr size_t rC = range_dim_cols;
-
-//  using typename XT::Grid::ElementBoundObject<Element>::ElementType;
-//  using DomainType = FieldVector<D, d>;
-//  using LocalBasisType = XT::Functions::ElementFunctionSetInterface<E, r, rC, R>;
-
-//  LocalUnaryElementIntegrandInterface(const XT::Common::ParameterType& param_type = {})
-//    : XT::Common::ParametricInterface(param_type)
-//  {}
-
-//  virtual ~LocalUnaryElementIntegrandInterface() = default;
-
-//  virtual std::unique_ptr<ThisType> copy() const = 0;
-
-// protected:
-//  virtual void post_bind(const IntersectionType& intrsctn) = 0;
-
-// public:
-
-//  /**
-//   * Returns the polynomial order of the integrand, given the basis.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual int order(const LocalBasisType& basis, const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * Computes the evaluation of this integrand at the given point for each function in the basis.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual void evaluate(const LocalBasisType& basis,
-//                        const DomainType& point_in_reference_element,
-//                        DynamicVector<F>& result,
-//                        const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * This method is provided for convenience and should not be used within library code.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual DynamicVector<F> evaluate(const LocalBasisType& basis,
-//                                    const DomainType& point_in_reference_element,
-//                                    const XT::Common::Parameter& param = {}) const
-//  {
-//    DynamicVector<F> result(basis.size(param));
-//    evaluate(basis, point_in_reference_element, result, param);
-//    return result;
-//  }
-//}; // class LocalUnaryElementIntegrandSum
-
-
-/// \todo add LocalBinaryElementIntegrandSum
-/// \todo add operator+ to LocalBinaryElementIntegrandInterface
-/// \sa LocalQuaternaryIntersectionIntegrandSum
-// template <class E, size_t t_r, size_t t_RC, class TF, class F, size_t a_r, size_t a_rC, class AF>
-// class LocalBinaryElementIntegrandSum
-//  : public XT::Common::ParametricInterface
-//  , public XT::Grid::ElementBoundObject<Element>
-//{
-//  static_assert(XT::Grid::is_entity<Element>::value, "");
-
-//  using ThisType = LocalBinaryElementIntegrandInterface<Element,
-//                                                        test_range_dim,
-//                                                        test_range_dim_cols,
-//                                                        TestRangeField,
-//                                                        Field,
-//                                                        ansatz_range_dim,
-//                                                        ansatz_range_dim_cols,
-//                                                        AnsatzRangeField>;
-
-// public:
-//  using E = Element;
-//  using D = typename Element::Geometry::ctype;
-//  static const constexpr size_t d = E::dimension;
-//  using F = Field;
-
-//  using TR = TestRangeField;
-//  static const constexpr size_t t_r = test_range_dim;
-//  static const constexpr size_t t_rC = test_range_dim_cols;
-
-//  using AR = AnsatzRangeField;
-//  static const constexpr size_t a_r = ansatz_range_dim;
-//  static const constexpr size_t a_rC = ansatz_range_dim_cols;
-
-//  using typename XT::Grid::ElementBoundObject<Element>::ElementType;
-//  using DomainType = FieldVector<D, d>;
-//  using LocalTestBasisType = XT::Functions::ElementFunctionSetInterface<E, t_r, t_rC, TR>;
-//  using LocalAnsatzBasisType = XT::Functions::ElementFunctionSetInterface<E, a_r, a_rC, AR>;
-
-//  LocalBinaryElementIntegrandInterface(const XT::Common::ParameterType& param_type = {})
-//    : XT::Common::ParametricInterface(param_type)
-//  {}
-
-//  virtual ~LocalBinaryElementIntegrandInterface() = default;
-
-//  virtual std::unique_ptr<ThisType> copy() const = 0;
-
-// protected:
-//  virtual void post_bind(const IntersectionType& intrsctn) = 0;
-
-// public:
-//  /**
-//   * Returns the polynomial order of the integrand, given the bases.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual int order(const LocalTestBasisType& test_basis,
-//                    const LocalAnsatzBasisType& ansatz_basis,
-//                    const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * Computes the evaluation of this integrand at the given point for each combination of functions from the two
-//   bases.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual void evaluate(const LocalTestBasisType& test_basis,
-//                        const LocalAnsatzBasisType& ansatz_basis,
-//                        const DomainType& point_in_reference_element,
-//                        DynamicMatrix<F>& result,
-//                        const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * This method is provided for convenience and should not be used within library code.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual DynamicMatrix<F> evaluate(const LocalTestBasisType& test_basis,
-//                                    const LocalAnsatzBasisType& ansatz_basis,
-//                                    const DomainType& point_in_reference_element,
-//                                    const XT::Common::Parameter& param = {}) const
-//  {
-//    DynamicMatrix<F> result(test_basis.size(param), ansatz_basis.size(param), 0);
-//    evaluate(test_basis, ansatz_basis, point_in_reference_element, result, param);
-//    return result;
-//  }
-//}; // class LocalBinaryElementIntegrandInterface
-
-
-/// \todo add LocalBinaryIntersectionIntegrandSum
-/// \todo add operator+ to LocalBinaryIntersectionIntegrandInterface
-/// \sa LocalQuaternaryIntersectionIntegrandSum
-// template <class I, size_t r, size_t rC, class R, class F>
-// class LocalBinaryIntersectionIntegrandSum
-//  : public XT::Common::ParametricInterface
-//  , public XT::Grid::IntersectionBoundObject<Intersection>
-//{
-//  static_assert(XT::Grid::is_intersection<Intersection>::value, "");
-
-//  using ThisType =
-//      LocalBinaryIntersectionIntegrandInterface<Intersection, range_dim, range_dim_cols, RangeField, Field>;
-
-// public:
-//  using typename XT::Grid::IntersectionBoundObject<Intersection>::IntersectionType;
-//  using ElementType = XT::Grid::extract_inside_element_t<Intersection>;
-
-//  using I = Intersection;
-//  using E = ElementType;
-//  using D = typename ElementType::Geometry::ctype;
-//  static const constexpr size_t d = E::dimension;
-//  using F = Field;
-
-//  using RF = RangeField;
-//  static const constexpr size_t r = range_dim;
-//  static const constexpr size_t rC = range_dim_cols;
-
-//  using DomainType = FieldVector<D, d - 1>;
-//  using LocalBasisType = XT::Functions::ElementFunctionSetInterface<E, r, rC, RF>;
-
-//  LocalBinaryIntersectionIntegrandInterface(const XT::Common::ParameterType& param_type = {})
-//    : XT::Common::ParametricInterface(param_type)
-//  {}
-
-//  virtual ~LocalBinaryIntersectionIntegrandInterface() = default;
-
-//  virtual std::unique_ptr<ThisType> copy() const = 0;
-
-// protected:
-//  virtual void post_bind(const IntersectionType& intrsctn) = 0;
-
-// public:
-//  /**
-//   * Returns the polynomial order of the integrand, given the bases.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual int order(const LocalBasisType& inside_basis,
-//                    const LocalBasisType& outside_basis,
-//                    const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * Computes the evaluation of this integrand at the given point for each combination of functions from the two
-//   bases.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual void evaluate(const LocalBasisType& inside_basis,
-//                        const LocalBasisType& outside_basis,
-//                        const DomainType& point_in_reference_element,
-//                        DynamicMatrix<F>& result,
-//                        const XT::Common::Parameter& param = {}) const = 0;
-
-//  /**
-//   * This method is provided for convenience and should not be used within library code.
-//   *
-//   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
-//   **/
-//  virtual DynamicMatrix<F> evaluate(const LocalBasisType& inside_basis,
-//                                    const LocalBasisType& outside_basis,
-//                                    const DomainType& point_in_reference_element,
-//                                    const XT::Common::Parameter& param = {}) const
-//  {
-//    DynamicMatrix<F> result(inside_basis.size(param), outside_basis.size(param), 0);
-//    evaluate(inside_basis, inside_basis, point_in_reference_element, result, param);
-//    return result;
-//  }
-//}; // class LocalBinaryIntersectionIntegrandSum
+template <class I, size_t r, size_t rC, class R, class F>
+class LocalUnaryIntersectionIntegrandSum : public LocalUnaryIntersectionIntegrandInterface<I, r, rC, R, F>
+{
+  using ThisType = LocalUnaryIntersectionIntegrandSum<I, r, rC, R, F>;
+  using BaseType = LocalUnaryIntersectionIntegrandInterface<I, r, rC, R, F>;
+
+public:
+  using typename BaseType::DomainType;
+  using typename BaseType::IntersectionType;
+  using typename BaseType::LocalBasisType;
+
+  LocalUnaryIntersectionIntegrandSum(const BaseType& left, const BaseType& right)
+    : BaseType(left.parameter_type() + right.parameter_type())
+    , left_(left.copy_as_unary_intersection_integrand().release())
+    , right_(right.copy_as_unary_intersection_integrand().release())
+  {
+    DUNE_THROW_IF(left.inside() != right.inside(),
+                  Exceptions::integrand_error,
+                  "left.inside() = " << left.inside() << "\n   right.inside() = " << right.inside());
+  }
+
+  LocalUnaryIntersectionIntegrandSum(const ThisType& other)
+    : BaseType(other.parameter_type())
+    , left_(other.left_.access().copy_as_unary_intersection_integrand().release())
+    , right_(other.right_.access().copy_as_unary_intersection_integrand().release())
+  {}
+
+  LocalUnaryIntersectionIntegrandSum(ThisType&& source) = default;
+
+  std::unique_ptr<BaseType> copy_as_unary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+protected:
+  void post_bind(const IntersectionType& intrsctn) override final
+  {
+    left_.access().bind(intrsctn);
+    right_.access().bind(intrsctn);
+  }
+
+public:
+  bool inside() const override final
+  {
+    return left_.access().inside();
+  }
+
+  int order(const LocalBasisType& basis, const XT::Common::Parameter& param = {}) const override final
+  {
+    return std::max(left_.access().order(basis, param), right_.access().order(basis, param));
+  }
+
+  using BaseType::evaluate;
+
+  void evaluate(const LocalBasisType& basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicVector<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // Each integrand clears its storage, so we let the left one write into ...
+    left_.access().evaluate(basis, point_in_reference_intersection, result, param);
+    // ..., the right one into ...
+    right_.access().evaluate(basis, point_in_reference_intersection, right_result_, param);
+    // ... and simply add them up (cannot use += here, vectors might have different sizes).
+    const size_t size = basis.size(param);
+    for (size_t ii = 0; ii < size; ++ii)
+      result[ii] += right_result_[ii];
+  } // ... evaluate(...)
+
+private:
+  XT::Common::StorageProvider<BaseType> left_;
+  XT::Common::StorageProvider<BaseType> right_;
+  mutable DynamicVector<F> right_result_;
+}; // class LocalUnaryIntersectionIntegrandSum
+
+
+template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+class LocalBinaryIntersectionIntegrandSum
+  : public LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>
+{
+  using BaseType = LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+  using ThisType = LocalBinaryIntersectionIntegrandSum<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+
+public:
+  using typename BaseType::DomainType;
+  using typename BaseType::IntersectionType;
+  using typename BaseType::LocalAnsatzBasisType;
+  using typename BaseType::LocalTestBasisType;
+
+  LocalBinaryIntersectionIntegrandSum(const BaseType& left, const BaseType& right)
+    : BaseType(left.parameter_type() + right.parameter_type())
+    , left_(left.copy_as_binary_intersection_integrand().release())
+    , right_(right.copy_as_binary_intersection_integrand().release())
+  {
+    DUNE_THROW_IF(left.inside() != right.inside(),
+                  Exceptions::integrand_error,
+                  "left.inside() = " << left.inside() << "\n   right.inside() = " << right.inside());
+  }
+
+  LocalBinaryIntersectionIntegrandSum(const ThisType& other)
+    : BaseType(other.parameter_type())
+    , left_(other.left_.access().copy_as_binary_intersection_integrand().release())
+    , right_(other.right_.access().copy_as_binary_intersection_integrand().release())
+  {}
+
+  LocalBinaryIntersectionIntegrandSum(ThisType&& source) = default;
+
+  std::unique_ptr<BaseType> copy_as_binary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+protected:
+  void post_bind(const IntersectionType& intrsctn) override final
+  {
+    left_.access().bind(intrsctn);
+    right_.access().bind(intrsctn);
+  }
+
+public:
+  bool inside() const override final
+  {
+    return left_.access().inside();
+  }
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
+            const XT::Common::Parameter& param = {}) const override final
+  {
+    return std::max(left_.access().order(test_basis, ansatz_basis, param),
+                    right_.access().order(test_basis, ansatz_basis, param));
+  }
+
+  using BaseType::evaluate;
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicMatrix<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // Each integrand clears its storage, so we let the left one write into ...
+    left_.access().evaluate(test_basis, ansatz_basis, point_in_reference_intersection, result, param);
+    // ..., the right one into ...
+    right_.access().evaluate(test_basis, ansatz_basis, point_in_reference_intersection, right_result_, param);
+    // ... and simply add them up (cannot use += here, matrices might have different sizes).
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    for (size_t ii = 0; ii < rows; ++ii)
+      for (size_t jj = 0; jj < cols; ++jj)
+        result[ii][jj] += right_result_[ii][jj];
+  } // ... evaluate(...)
+
+private:
+  XT::Common::StorageProvider<BaseType> left_;
+  XT::Common::StorageProvider<BaseType> right_;
+  mutable DynamicMatrix<F> right_result_;
+}; // class LocalBinaryIntersectionIntegrandSum
+
+
+template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+class LocalUnaryAndBinaryIntersectionIntegrandSum
+  : public LocalUnaryAndBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>
+{
+  using BaseType = LocalUnaryAndBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+  using ThisType = LocalUnaryAndBinaryIntersectionIntegrandSum<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+
+public:
+  using typename BaseType::DomainType;
+  using typename BaseType::IntersectionType;
+  using typename BaseType::LocalAnsatzBasisType;
+  using typename BaseType::LocalTestBasisType;
+
+  LocalUnaryAndBinaryIntersectionIntegrandSum(const BaseType& left, const BaseType& right)
+    : BaseType(left.parameter_type() + right.parameter_type())
+    , left_(left.copy_as_binary_intersection_integrand().release())
+    , right_(right.copy_as_binary_intersection_integrand().release())
+  {
+    DUNE_THROW_IF(left.inside() != right.inside(),
+                  Exceptions::integrand_error,
+                  "left.inside() = " << left.inside() << "\n   right.inside() = " << right.inside());
+  }
+
+  LocalUnaryAndBinaryIntersectionIntegrandSum(const ThisType& other)
+    : BaseType(other.parameter_type())
+    , left_(other.left_.access().copy_as_binary_intersection_integrand().release())
+    , right_(other.right_.access().copy_as_binary_intersection_integrand().release())
+  {}
+
+  LocalUnaryAndBinaryIntersectionIntegrandSum(ThisType&& source) = default;
+
+  std::unique_ptr<BaseType> copy_as_unary_and_binary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+protected:
+  void post_bind(const IntersectionType& intrsctn) override final
+  {
+    left_.access().bind(intrsctn);
+    right_.access().bind(intrsctn);
+  }
+
+public:
+  bool inside() const override final
+  {
+    return left_.access().inside();
+  }
+
+  using BaseType::evaluate;
+
+  /// \name Required for LocalUnaryIntersectionIntegrandInterface
+  /// \{
+
+  int order(const LocalTestBasisType& test_basis, const XT::Common::Parameter& param = {}) const override final
+  {
+    return std::max(left_.access().order(test_basis, param), right_.access().order(test_basis, param));
+  }
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicVector<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // Each integrand clears its storage, so we let the left one write into ...
+    left_.access().evaluate(test_basis, point_in_reference_intersection, result, param);
+    // ..., the right one into ...
+    right_.access().evaluate(test_basis, point_in_reference_intersection, right_unary_result_, param);
+    // ... and simply add them up (cannot use += here, vectors might have different sizes).
+    const size_t size = test_basis.size(param);
+    for (size_t ii = 0; ii < size; ++ii)
+      result[ii] += right_unary_result_[ii];
+  } // ... evaluate(...)
+
+  /// \}
+  /// \name Required for LocalBinaryIntersectionIntegrandInterface
+  /// \{
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
+            const XT::Common::Parameter& param = {}) const override final
+  {
+    return std::max(left_.access().order(test_basis, ansatz_basis, param),
+                    right_.access().order(test_basis, ansatz_basis, param));
+  }
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicMatrix<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // Each integrand clears its storage, so we let the left one write into ...
+    left_.access().evaluate(test_basis, ansatz_basis, point_in_reference_intersection, result, param);
+    // ..., the right one into ...
+    right_.access().evaluate(test_basis, ansatz_basis, point_in_reference_intersection, right_binary_result_, param);
+    // ... and simply add them up (cannot use += here, matrices might have different sizes).
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    for (size_t ii = 0; ii < rows; ++ii)
+      for (size_t jj = 0; jj < cols; ++jj)
+        result[ii][jj] += right_binary_result_[ii][jj];
+  } // ... evaluate(...)
+
+  /// \}
+private:
+  XT::Common::StorageProvider<BaseType> left_;
+  XT::Common::StorageProvider<BaseType> right_;
+  mutable DynamicVector<F> right_unary_result_;
+  mutable DynamicMatrix<F> right_binary_result_;
+}; // class LocalUnaryAndBinaryIntersectionIntegrandSum
 
 
 template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
@@ -273,19 +306,19 @@ public:
 
   LocalQuaternaryIntersectionIntegrandSum(const BaseType& left, const BaseType& right)
     : BaseType(left.parameter_type() + right.parameter_type())
-    , left_(left.copy().release())
-    , right_(right.copy().release())
+    , left_(left.copy_as_quaternary_intersection_integrand().release())
+    , right_(right.copy_as_quaternary_intersection_integrand().release())
   {}
 
   LocalQuaternaryIntersectionIntegrandSum(const ThisType& other)
     : BaseType(other)
-    , left_(other.left_.access().copy().release())
-    , right_(other.right_.access().copy().release())
+    , left_(other.left_.access().copy_as_quaternary_intersection_integrand().release())
+    , right_(other.right_.access().copy_as_quaternary_intersection_integrand().release())
   {}
 
   LocalQuaternaryIntersectionIntegrandSum(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
diff --git a/dune/gdt/local/integrands/conversion.hh b/dune/gdt/local/integrands/conversion.hh
index 2f1ec208c..2b77c96c1 100644
--- a/dune/gdt/local/integrands/conversion.hh
+++ b/dune/gdt/local/integrands/conversion.hh
@@ -12,6 +12,7 @@
 #define DUNE_GDT_LOCAL_INTEGRANDS_CONVERSION_HH
 
 #include <dune/xt/functions/interfaces/grid-function.hh>
+#include <dune/xt/functions/generic/element-function.hh>
 
 #include "interfaces.hh"
 
@@ -20,14 +21,15 @@ namespace GDT {
 
 
 /**
- * Given a function f and a binary element integrand bi(\cdot, cdot), models the unary element integrand bi(f, \cdot).
+ * Given a function f and a binary element integrand bi(\cdot, \cdot), models the unary element integrand bi(f, \cdot).
+ *
+ * Most likey you do not want to use this class directly, but LocalUnaryElementIntegrandInterface::with_ansatz.
  *
  * See also LocalUnaryElementIntegrandInterface for a description of the template arguments.
  *
- * \sa local_binary_to_unary_element_integrand
  * \sa LocalUnaryElementIntegrandInterface
  * \sa LocalBinaryElementIntegrandInterface
- * \sa XT::Functions::GridFunctionInterface
+ * \sa XT::Functions::GridFunction
  */
 template <class E,
           size_t t_r = 1,
@@ -37,34 +39,33 @@ template <class E,
           size_t a_r = t_r,
           size_t a_rC = t_rC,
           class AF = TF>
-class LocalBinaryToUnaryElementIntegrand : public LocalUnaryElementIntegrandInterface<E, a_r, a_rC, AF, F>
+class LocalBinaryToUnaryElementIntegrand : public LocalUnaryElementIntegrandInterface<E, t_r, t_rC, TF, F>
 {
   using ThisType = LocalBinaryToUnaryElementIntegrand<E, t_r, t_rC, TF, F, a_r, a_rC, AF>;
-  using BaseType = LocalUnaryElementIntegrandInterface<E, a_r, a_rC, AF, F>;
+  using BaseType = LocalUnaryElementIntegrandInterface<E, t_r, t_rC, TF, F>;
 
 public:
   using typename BaseType::DomainType;
   using typename BaseType::ElementType;
   using typename BaseType::LocalBasisType;
 
-  using LocalizableFunctionType = XT::Functions::GridFunctionInterface<E, t_r, t_rC, TF>;
   using LocalBinaryElementIntegrandType = LocalBinaryElementIntegrandInterface<E, t_r, t_rC, TF, F, a_r, a_rC, AF>;
 
-  LocalBinaryToUnaryElementIntegrand(const LocalizableFunctionType& inducing_function_as_test_basis,
-                                     const LocalBinaryElementIntegrandType& local_binary_integrand)
-    : inducing_function_as_test_basis_(inducing_function_as_test_basis)
-    , local_function_(inducing_function_as_test_basis_.local_function())
-    , local_binary_integrand_(local_binary_integrand.copy())
+  LocalBinaryToUnaryElementIntegrand(const LocalBinaryElementIntegrandType& local_binary_integrand,
+                                     XT::Functions::GridFunction<E, a_r, a_rC, AF> inducing_function_as_ansatz_basis)
+    : inducing_function_as_ansatz_basis_(inducing_function_as_ansatz_basis)
+    , local_function_(inducing_function_as_ansatz_basis_.local_function())
+    , local_binary_integrand_(local_binary_integrand.copy_as_binary_element_integrand())
   {}
 
   LocalBinaryToUnaryElementIntegrand(const ThisType& other)
     : BaseType()
-    , inducing_function_as_test_basis_(other.inducing_function_as_test_basis_)
-    , local_function_(inducing_function_as_test_basis_.local_function())
-    , local_binary_integrand_(other.local_binary_integrand_->copy())
+    , inducing_function_as_ansatz_basis_(other.inducing_function_as_ansatz_basis_)
+    , local_function_(inducing_function_as_ansatz_basis_.local_function())
+    , local_binary_integrand_(other.local_binary_integrand_->copy_as_binary_element_integrand())
   {}
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_unary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -79,7 +80,7 @@ protected:
 public:
   int order(const LocalBasisType& basis, const XT::Common::Parameter& param = {}) const override final
   {
-    return local_binary_integrand_->order(*local_function_, basis, param);
+    return local_binary_integrand_->order(basis, *local_function_, param);
   }
 
   using BaseType::evaluate;
@@ -95,14 +96,14 @@ public:
       result.resize(size);
     // evaluate
     local_binary_integrand_->evaluate(
-        *local_function_, basis, point_in_reference_element, local_binary_integrand_result_, param);
+        basis, *local_function_, point_in_reference_element, local_binary_integrand_result_, param);
     // extract result
     result = local_binary_integrand_result_[0];
   } // ... evaluate(...)
 
 private:
-  const LocalizableFunctionType& inducing_function_as_test_basis_;
-  std::unique_ptr<typename LocalizableFunctionType::LocalFunctionType> local_function_;
+  const XT::Functions::GridFunction<E, a_r, a_rC, AF> inducing_function_as_ansatz_basis_;
+  std::unique_ptr<typename XT::Functions::GridFunction<E, a_r, a_rC, AF>::LocalFunctionType> local_function_;
   std::unique_ptr<LocalBinaryElementIntegrandType> local_binary_integrand_;
   mutable DynamicMatrix<F> local_binary_integrand_result_;
 }; // class LocalBinaryToUnaryElementIntegrand
@@ -110,15 +111,94 @@ private:
 
 /**
  * \sa LocalBinaryToUnaryElementIntegrand
+ * \sa LocalUnaryIntersectionIntegrandInterface
+ * \sa LocalBinaryIntersectionIntegrandInterface
+ * \sa XT::Functions::GridFunction
  */
-template <class E, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
-LocalBinaryToUnaryElementIntegrand<E, t_r, t_rC, TF, F, a_r, a_rC, AF> local_binary_to_unary_element_integrand(
-    const LocalBinaryElementIntegrandInterface<E, t_r, t_rC, TF, F, a_r, a_rC, AF>& local_binary_element_integrand,
-    const XT::Functions::GridFunctionInterface<E, t_r, t_rC, TF>& inducing_function_as_test_basis)
+template <class I,
+          size_t t_r = 1,
+          size_t t_rC = 1,
+          class TF = double,
+          class F = double,
+          size_t a_r = t_r,
+          size_t a_rC = t_rC,
+          class AF = TF>
+class LocalBinaryToUnaryIntersectionIntegrand : public LocalUnaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F>
 {
-  return LocalBinaryToUnaryElementIntegrand<E, t_r, t_rC, TF, F, a_r, a_rC, AF>(inducing_function_as_test_basis,
-                                                                                local_binary_element_integrand);
-}
+  using ThisType = LocalBinaryToUnaryIntersectionIntegrand<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+  using BaseType = LocalUnaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F>;
+
+public:
+  using typename BaseType::DomainType;
+  using typename BaseType::E;
+  using typename BaseType::IntersectionType;
+  using typename BaseType::LocalBasisType;
+
+  using LocalBinaryIntersectionIntegrandType =
+      LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+
+  LocalBinaryToUnaryIntersectionIntegrand(
+      const LocalBinaryIntersectionIntegrandType& local_binary_integrand,
+      XT::Functions::GridFunction<E, a_r, a_rC, AF> inducing_function_as_ansatz_basis)
+    : inducing_function_as_ansatz_basis_(inducing_function_as_ansatz_basis)
+    , local_function_(inducing_function_as_ansatz_basis_.local_function())
+    , local_binary_integrand_(local_binary_integrand.copy_as_binary_intersection_integrand())
+  {}
+
+  LocalBinaryToUnaryIntersectionIntegrand(const ThisType& other)
+    : BaseType()
+    , inducing_function_as_ansatz_basis_(other.inducing_function_as_ansatz_basis_)
+    , local_function_(inducing_function_as_ansatz_basis_.local_function())
+    , local_binary_integrand_(other.local_binary_integrand_->copy_as_binary_intersection_integrand())
+  {}
+
+  std::unique_ptr<BaseType> copy_as_unary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+  bool inside() const override final
+  {
+    return local_binary_integrand_->inside();
+  }
+
+protected:
+  void post_bind(const IntersectionType& intersctn) override final
+  {
+    local_function_->bind(this->inside() ? intersctn.inside() : intersctn.outside());
+    local_binary_integrand_->bind(intersctn);
+  }
+
+public:
+  int order(const LocalBasisType& test_basis, const XT::Common::Parameter& param = {}) const override final
+  {
+    return local_binary_integrand_->order(test_basis, *local_function_, param);
+  }
+
+  using BaseType::evaluate;
+
+  void evaluate(const LocalBasisType& test_basis,
+                const DomainType& point_in_reference_Intersection,
+                DynamicVector<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // prepare storage
+    const auto size = test_basis.size(param);
+    if (result.size() < size)
+      result.resize(size);
+    // evaluate
+    local_binary_integrand_->evaluate(
+        test_basis, *local_function_, point_in_reference_Intersection, local_binary_integrand_result_, param);
+    // extract result
+    result = local_binary_integrand_result_[0];
+  } // ... evaluate(...)
+
+private:
+  const XT::Functions::GridFunction<E, a_r, a_rC, AF> inducing_function_as_ansatz_basis_;
+  std::unique_ptr<typename XT::Functions::GridFunction<E, a_r, a_rC, AF>::LocalFunctionType> local_function_;
+  std::unique_ptr<LocalBinaryIntersectionIntegrandType> local_binary_integrand_;
+  mutable DynamicMatrix<F> local_binary_integrand_result_;
+}; // class LocalBinaryToUnaryIntersectionIntegrand
 
 
 } // namespace GDT
diff --git a/dune/gdt/local/integrands/div.hh b/dune/gdt/local/integrands/div.hh
index 8c1629f6c..acf8dcd14 100644
--- a/dune/gdt/local/integrands/div.hh
+++ b/dune/gdt/local/integrands/div.hh
@@ -72,8 +72,6 @@ protected:
  * bases.
  *
  * \note Note that f can also be given as a scalar value or omitted.
- *
- * \sa local_binary_to_unary_element_integrand
  */
 template <class E, class TR = double, class F = double, class AR = TR>
 class LocalElementAnsatzValueTestDivProductIntegrand
@@ -123,7 +121,7 @@ public:
 
   LocalElementAnsatzValueTestDivProductIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -218,7 +216,7 @@ public:
 
   LocalElementAnsatzDivTestValueProductIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
diff --git a/dune/gdt/local/integrands/elliptic-ipdg.hh b/dune/gdt/local/integrands/elliptic-ipdg.hh
index 1308b9ea7..bad91885a 100644
--- a/dune/gdt/local/integrands/elliptic-ipdg.hh
+++ b/dune/gdt/local/integrands/elliptic-ipdg.hh
@@ -81,7 +81,7 @@ public:
 
   InnerCoupling(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -253,7 +253,7 @@ public:
 
   DirichletCoupling(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -558,7 +558,7 @@ public:
 
   Inner(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -1008,7 +1008,7 @@ private:
  */
 template <class I, class F = double, Method method = default_method>
 class DXT_DEPRECATED_MSG(
-    "Use LocalLaplaceIPDGIntegrands::DirichletCoupling + LocalIPDGIntegrands::boundaryPenalty instead (10.08.2019)!")
+    "Use LocalLaplaceIPDGIntegrands::DirichletCoupling + LocalIPDGIntegrands::Penalty instead (10.08.2019)!")
     DirichletBoundaryLhs : public LocalQuaternaryIntersectionIntegrandInterface<I, 1, 1, F, F, 1, 1, F>
 {
   using BaseType = LocalQuaternaryIntersectionIntegrandInterface<I, 1, 1, F, F, 1, 1, F>;
@@ -1047,7 +1047,7 @@ public:
 
   DirichletBoundaryLhs(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -1492,7 +1492,7 @@ public:
 
   InnerOnlyPenalty(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -1901,8 +1901,8 @@ private:
  * \sa [Epshteyn, Riviere, 2007] for the meaning of beta
  */
 template <class I, class F = double, Method method = default_method>
-class DXT_DEPRECATED_MSG("Use LocalIPDGIntegrands::BoundaryPenalty instead (05.08.2019)!")
-    DirichletBoundaryLhsOnlyPenalty : public LocalQuaternaryIntersectionIntegrandInterface<I, 1, 1, F, F, 1, 1, F>
+class DXT_DEPRECATED_MSG("Use LocalIPDGIntegrands::Penalty instead (05.08.2019)!") DirichletBoundaryLhsOnlyPenalty
+  : public LocalQuaternaryIntersectionIntegrandInterface<I, 1, 1, F, F, 1, 1, F>
 {
   using BaseType = LocalQuaternaryIntersectionIntegrandInterface<I, 1, 1, F, F, 1, 1, F>;
   using ThisType = DirichletBoundaryLhsOnlyPenalty<I, F, method>;
@@ -1940,7 +1940,7 @@ public:
 
   DirichletBoundaryLhsOnlyPenalty(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
diff --git a/dune/gdt/local/integrands/elliptic.hh b/dune/gdt/local/integrands/elliptic.hh
index 9cf927b79..1d80582d8 100644
--- a/dune/gdt/local/integrands/elliptic.hh
+++ b/dune/gdt/local/integrands/elliptic.hh
@@ -89,7 +89,7 @@ public:
 
   LocalEllipticIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
diff --git a/dune/gdt/local/integrands/generic.hh b/dune/gdt/local/integrands/generic.hh
index d6597958a..2e27e644c 100644
--- a/dune/gdt/local/integrands/generic.hh
+++ b/dune/gdt/local/integrands/generic.hh
@@ -52,7 +52,7 @@ public:
     , evaluate_(other.evaluate_)
   {}
 
-  std::unique_ptr<BaseType> copy() const
+  std::unique_ptr<BaseType> copy_as_unary_element_integrand() const
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -129,7 +129,7 @@ public:
     , evaluate_(other.evaluate_)
   {}
 
-  std::unique_ptr<BaseType> copy() const
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -170,76 +170,94 @@ private:
 }; // class GenericLocalBinaryElementIntegrand
 
 
-template <class I, size_t r = 1, size_t rC = 1, class RF = double, class F = double>
-class GenericLocalBinaryIntersectionIntegrand : public LocalBinaryIntersectionIntegrandInterface<I, r, rC, RF, F>
+template <class I,
+          size_t t_r = 1,
+          size_t t_rC = 1,
+          class TF = double,
+          class F = double,
+          size_t a_r = t_r,
+          size_t a_rC = t_rC,
+          class AF = TF>
+class GenericLocalBinaryIntersectionIntegrand
+  : public LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>
 {
-  using BaseType = LocalBinaryIntersectionIntegrandInterface<I, r, rC, RF, F>;
-  using ThisType = GenericLocalBinaryIntersectionIntegrand<I, r, rC, RF, F>;
+  using BaseType = LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
+  using ThisType = GenericLocalBinaryIntersectionIntegrand<I, t_r, t_rC, TF, F, a_r, a_rC, AF>;
 
 public:
   using typename BaseType::DomainType;
-  using typename BaseType::LocalBasisType;
+  using typename BaseType::LocalAnsatzBasisType;
+  using typename BaseType::LocalTestBasisType;
 
-  using GenericOrderFunctionType = std::function<int(const LocalBasisType& /*inside_basis*/,
-                                                     const LocalBasisType& /*outside_basis*/,
+  using GenericOrderFunctionType = std::function<int(const LocalTestBasisType& /*test_basis*/,
+                                                     const LocalAnsatzBasisType& /*ansatz_basis*/,
                                                      const XT::Common::Parameter& /*param*/)>;
-  using GenericEvalauteFunctionType = std::function<void(const LocalBasisType& /*inside_basis*/,
-                                                         const LocalBasisType& /*outside_basis*/,
+  using GenericEvalauteFunctionType = std::function<void(const LocalTestBasisType& /*test_basis*/,
+                                                         const LocalAnsatzBasisType& /*ansatz_basis*/,
                                                          const DomainType& /*point_in_reference_intersection*/,
                                                          DynamicMatrix<F>& /*result*/,
                                                          const XT::Common::Parameter& /*param*/)>;
 
   GenericLocalBinaryIntersectionIntegrand(GenericOrderFunctionType order_function,
                                           GenericEvalauteFunctionType evaluate_function,
+                                          bool use_inside_bases = true,
                                           const XT::Common::ParameterType& param_type = {})
     : BaseType(param_type)
     , order_(order_function)
     , evaluate_(evaluate_function)
+    , inside_(use_inside_bases)
   {}
 
   GenericLocalBinaryIntersectionIntegrand(const ThisType& other)
     : BaseType(other.parameter_type())
     , order_(other.order_)
     , evaluate_(other.evaluate_)
+    , inside_(other.inside_)
   {}
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
 
-  int order(const LocalBasisType& inside_basis,
-            const LocalBasisType& outside_basis,
+  bool inside() const override final
+  {
+    return inside_;
+  }
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
             const XT::Common::Parameter& param = {}) const override final
   {
-    return order_(inside_basis, outside_basis, this->parse_parameter(param));
+    return order_(test_basis, ansatz_basis, this->parse_parameter(param));
   }
 
   using BaseType::evaluate;
 
-  void evaluate(const LocalBasisType& inside_basis,
-                const LocalBasisType& outside_basis,
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
                 const DomainType& point_in_reference_intersection,
                 DynamicMatrix<F>& result,
                 const XT::Common::Parameter& param = {}) const override final
   { // prepare storage
-    const size_t rows = inside_basis.size(param);
-    const size_t cols = outside_basis.size(param);
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
     if (result.rows() < rows || result.cols() < cols)
       result.resize(rows, cols);
     // evaluate
-    evaluate_(inside_basis, outside_basis, point_in_reference_intersection, result, this->parse_parameter(param));
+    evaluate_(test_basis, ansatz_basis, point_in_reference_intersection, result, this->parse_parameter(param));
     // check
     DUNE_THROW_IF(result.rows() < rows || result.cols() < cols,
                   Exceptions::integrand_error,
-                  "inside_basis.size(param) = " << rows << "\n   result.rows() = " << result.rows()
-                                                << "outside_basis.size(param) = " << cols
-                                                << "\n   result.cols() = " << result.cols());
+                  "test_basis.size(param) = " << rows << "\n   result.rows() = " << result.rows()
+                                              << "ansatz_basis.size(param) = " << cols
+                                              << "\n   result.cols() = " << result.cols());
   } // ... evaluate(...)
 
 private:
   const GenericOrderFunctionType order_;
   const GenericEvalauteFunctionType evaluate_;
+  const bool inside_;
 }; // class GenericLocalBinaryIntersectionIntegrand
 
 
diff --git a/dune/gdt/local/integrands/identity.hh b/dune/gdt/local/integrands/identity.hh
index 663822eba..83107230c 100644
--- a/dune/gdt/local/integrands/identity.hh
+++ b/dune/gdt/local/integrands/identity.hh
@@ -31,7 +31,7 @@ public:
   using typename BaseType::DomainType;
   using typename BaseType::LocalBasisType;
 
-  std::unique_ptr<BaseType> copy() const
+  std::unique_ptr<BaseType> copy_as_unary_element_integrand() const
   {
     return std::make_unique<ThisType>();
   }
diff --git a/dune/gdt/local/integrands/interfaces.hh b/dune/gdt/local/integrands/interfaces.hh
index 96798d086..92a9fc9b6 100644
--- a/dune/gdt/local/integrands/interfaces.hh
+++ b/dune/gdt/local/integrands/interfaces.hh
@@ -19,6 +19,7 @@
 #include <dune/xt/grid/bound-object.hh>
 #include <dune/xt/grid/type_traits.hh>
 #include <dune/xt/functions/interfaces/element-functions.hh>
+#include <dune/xt/functions/grid-function.hh>
 
 #include <dune/gdt/exceptions.hh>
 
@@ -26,20 +27,34 @@ namespace Dune {
 namespace GDT {
 
 
-// forwards (required for operator+, include are below)
-template <class E, size_t r, size_t rC, class R, class F>
-class LocalUnaryElementIntegrandSum;
+// forwards (required for operator+), includes are below
+// template <class E, size_t r, size_t rC, class R, class F>
+// class LocalUnaryElementIntegrandSum;
 
-template <class E, size_t t_r, size_t t_RC, class TF, class F, size_t a_r, size_t a_rC, class AF>
-class LocalBinaryElementIntegrandSum;
+// template <class E, size_t t_r, size_t t_RC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+// class LocalBinaryElementIntegrandSum;
 
 template <class I, size_t r, size_t rC, class R, class F>
+class LocalUnaryIntersectionIntegrandSum;
+
+template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
 class LocalBinaryIntersectionIntegrandSum;
 
+template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+class LocalUnaryAndBinaryIntersectionIntegrandSum;
+
 template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
 class LocalQuaternaryIntersectionIntegrandSum;
 
 
+// forwards (required for with_ansatz), includes are below
+template <class E, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+class LocalBinaryToUnaryElementIntegrand;
+
+template <class I, size_t t_r, size_t t_rC, class TF, class F, size_t a_r, size_t a_rC, class AF>
+class LocalBinaryToUnaryIntersectionIntegrand;
+
+
 /**
  * Interface for integrands in integrals over grid elements, which depend on one argument only (usually the test basis
  * in an integral-based functional).
@@ -81,7 +96,7 @@ public:
 
   virtual ~LocalUnaryElementIntegrandInterface() = default;
 
-  virtual std::unique_ptr<ThisType> copy() const = 0;
+  virtual std::unique_ptr<ThisType> copy_as_unary_element_integrand() const = 0;
 
   /**
    * Returns the polynomial order of the integrand, given the basis.
@@ -172,7 +187,13 @@ public:
 
   virtual ~LocalBinaryElementIntegrandInterface() = default;
 
-  virtual std::unique_ptr<ThisType> copy() const = 0;
+  template <class... Args>
+  LocalBinaryToUnaryElementIntegrand<E, t_r, t_rC, TR, F, a_r, a_rC, AR> with_ansatz(Args&&... args) const
+  {
+    return LocalBinaryToUnaryElementIntegrand<E, t_r, t_rC, TR, F, a_r, a_rC, AR>(*this, std::forward<Args>(args)...);
+  }
+
+  virtual std::unique_ptr<ThisType> copy_as_binary_element_integrand() const = 0;
 
   /**
    * Returns the polynomial order of the integrand, given the bases.
@@ -208,12 +229,85 @@ public:
     evaluate(test_basis, ansatz_basis, point_in_reference_element, result, param);
     return result;
   }
+
+protected:
+  void ensure_size_and_clear_results(const LocalTestBasisType& test_basis,
+                                     const LocalAnsatzBasisType& ansatz_basis,
+                                     DynamicMatrix<F>& result,
+                                     const XT::Common::Parameter& param) const
+  {
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    if (result.rows() < rows || result.cols() < cols)
+      result.resize(rows, cols);
+    result *= 0;
+  } // ... ensure_size_and_clear_results(...)
 }; // class LocalBinaryElementIntegrandInterface
 
 
+template <class I,
+          size_t t_r = 1,
+          size_t t_rC = 1,
+          class TR = double,
+          class F_ = double,
+          size_t a_r = t_r,
+          size_t a_rC = t_rC,
+          class AR = TR>
+class LocalUnaryAndBinaryElementIntegrandInterface
+  : public LocalUnaryElementIntegrandInterface<I, t_r, t_rC, TR, F_>
+  , public LocalBinaryElementIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>
+{
+
+  using ThisType = LocalUnaryAndBinaryElementIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>;
+
+protected:
+  using UnaryBaseType = LocalUnaryElementIntegrandInterface<I, t_r, t_rC, TR, F_>;
+  using BinaryBaseType = LocalBinaryElementIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>;
+
+public:
+  /// \name Members and typedefs required for disambiuation.
+  /// \{
+
+  using typename UnaryBaseType::D;
+  using typename UnaryBaseType::DomainType;
+  using typename UnaryBaseType::E;
+  using typename UnaryBaseType::ElementType;
+  using typename UnaryBaseType::F;
+  using UnaryBaseType::d;
+
+  /// \}
+
+  LocalUnaryAndBinaryElementIntegrandInterface(const XT::Common::ParameterType& param_type = {})
+    : UnaryBaseType(param_type)
+    , BinaryBaseType(param_type)
+  {}
+
+  virtual ~LocalUnaryAndBinaryElementIntegrandInterface() = default;
+
+  virtual std::unique_ptr<ThisType> copy_as_unary_and_binary_element_integrand() const = 0;
+
+  //  using UnaryBaseType::operator+;
+  //  using BinaryBaseType::operator+;
+
+  //  LocalUnaryAndBinaryElementIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR> operator+(const ThisType& other) const
+  //  {
+  //    return LocalUnaryAndBinaryElementIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR>(*this, other);
+  //  }
+
+  /// \name Methods required for disambiuation.
+  /// \{
+
+  using UnaryBaseType::is_parametric;
+  using UnaryBaseType::parameter_type;
+  using UnaryBaseType::parse_parameter;
+
+  /// \}
+}; // class LocalUnaryAndBinaryElementIntegrandInterface
+
+
 /**
- * Interface for integrands in integrals over grid intersections, which depend on two arguments (usually the bases on
- * the inside and outside the in an integral-based functional).
+ * Interface for integrands in integrals over grid intersections, which depend on one argument (usually the test basis
+ * on the inside of the intersection in an integral-based functional).
  *
  * \note Regarding SMP: the integrand is copied for each thread, so
  *       - no shared mutable state between copies to be thread safe, but
@@ -224,14 +318,13 @@ template <class Intersection,
           size_t range_dim_cols = 1,
           class RangeField = double,
           class Field = double>
-class LocalBinaryIntersectionIntegrandInterface
-  : public XT::Common::ParametricInterface
-  , public XT::Grid::IntersectionBoundObject<Intersection>
+class LocalUnaryIntersectionIntegrandInterface
+  : virtual public XT::Common::ParametricInterface
+  , virtual public XT::Grid::IntersectionBoundObject<Intersection>
 {
   static_assert(XT::Grid::is_intersection<Intersection>::value, "");
 
-  using ThisType =
-      LocalBinaryIntersectionIntegrandInterface<Intersection, range_dim, range_dim_cols, RangeField, Field>;
+  using ThisType = LocalUnaryIntersectionIntegrandInterface<Intersection, range_dim, range_dim_cols, RangeField, Field>;
 
 public:
   using typename XT::Grid::IntersectionBoundObject<Intersection>::IntersectionType;
@@ -250,21 +343,159 @@ public:
   using DomainType = FieldVector<D, d - 1>;
   using LocalBasisType = XT::Functions::ElementFunctionSetInterface<E, r, rC, RF>;
 
+  LocalUnaryIntersectionIntegrandInterface(const XT::Common::ParameterType& param_type = {})
+    : XT::Common::ParametricInterface(param_type)
+  {}
+
+  virtual ~LocalUnaryIntersectionIntegrandInterface() = default;
+
+  virtual std::unique_ptr<ThisType> copy_as_unary_intersection_integrand() const = 0;
+
+  LocalUnaryIntersectionIntegrandSum<I, r, rC, RF, F> operator+(const ThisType& other) const
+  {
+    return LocalUnaryIntersectionIntegrandSum<I, r, rC, RF, F>(*this, other);
+  }
+
+  /**
+   * Flag to document which element the basis is expected to be bound to.
+   */
+  virtual bool inside() const
+  {
+    return true;
+  }
+
+  /**
+   * Returns the polynomial order of the integrand, given the bases.
+   *
+   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
+   **/
+  virtual int order(const LocalBasisType& basis, const XT::Common::Parameter& param = {}) const = 0;
+
+  /**
+   * Computes the evaluation of this integrand at the given point for each basis function.
+   *
+   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
+   **/
+  virtual void evaluate(const LocalBasisType& basis,
+                        const DomainType& point_in_reference_intersection,
+                        DynamicVector<F>& result,
+                        const XT::Common::Parameter& param = {}) const = 0;
+
+  /**
+   * This method is provided for convenience and should not be used within library code.
+   *
+   * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
+   **/
+  virtual DynamicVector<F> evaluate(const LocalBasisType& basis,
+                                    const DomainType& point_in_reference_intersection,
+                                    const XT::Common::Parameter& param = {}) const
+  {
+    DynamicVector<F> result(basis.size(param), 0);
+    evaluate(basis, point_in_reference_intersection, result, param);
+    return result;
+  }
+
+protected:
+  void ensure_size_and_clear_results(const LocalBasisType& basis,
+                                     DynamicVector<F>& result,
+                                     const XT::Common::Parameter& param) const
+  {
+    const size_t size = basis.size(param);
+    if (result.size() < size)
+      result.resize(size);
+    result *= 0;
+  }
+}; // class LocalUnaryIntersectionIntegrandInterface
+
+
+/**
+ * Interface for integrands in integrals over grid intersections, which depend on two arguments (usually the test and
+ * ansatz bases on the inside of the intersection in an integral-based functional).
+ *
+ * \note Regarding SMP: the integrand is copied for each thread, so
+ *       - no shared mutable state between copies to be thread safe, but
+ *       - local mutable state is ok.
+ */
+template <class Intersection,
+          size_t test_range_dim = 1,
+          size_t test_range_dim_cols = 1,
+          class TestRangeField = double,
+          class Field = double,
+          size_t ansatz_range_dim = test_range_dim,
+          size_t ansatz_range_dim_cols = test_range_dim_cols,
+          class AnsatzRangeField = TestRangeField>
+class LocalBinaryIntersectionIntegrandInterface
+  : virtual public XT::Common::ParametricInterface
+  , virtual public XT::Grid::IntersectionBoundObject<Intersection>
+{
+  static_assert(XT::Grid::is_intersection<Intersection>::value, "");
+
+  using ThisType = LocalBinaryIntersectionIntegrandInterface<Intersection,
+                                                             test_range_dim,
+                                                             test_range_dim_cols,
+                                                             TestRangeField,
+                                                             Field,
+                                                             ansatz_range_dim,
+                                                             ansatz_range_dim_cols,
+                                                             AnsatzRangeField>;
+
+public:
+  using typename XT::Grid::IntersectionBoundObject<Intersection>::IntersectionType;
+  using ElementType = XT::Grid::extract_inside_element_t<Intersection>;
+
+  using I = Intersection;
+  using E = ElementType;
+  using D = typename ElementType::Geometry::ctype;
+  static const constexpr size_t d = E::dimension;
+  using F = Field;
+
+  using TR = TestRangeField;
+  static const constexpr size_t t_r = test_range_dim;
+  static const constexpr size_t t_rC = test_range_dim_cols;
+
+  using AR = AnsatzRangeField;
+  static const constexpr size_t a_r = ansatz_range_dim;
+  static const constexpr size_t a_rC = ansatz_range_dim_cols;
+
+  using DomainType = FieldVector<D, d - 1>;
+  using LocalTestBasisType = XT::Functions::ElementFunctionSetInterface<E, t_r, t_rC, TR>;
+  using LocalAnsatzBasisType = XT::Functions::ElementFunctionSetInterface<E, a_r, a_rC, AR>;
+
   LocalBinaryIntersectionIntegrandInterface(const XT::Common::ParameterType& param_type = {})
     : XT::Common::ParametricInterface(param_type)
   {}
 
   virtual ~LocalBinaryIntersectionIntegrandInterface() = default;
 
-  virtual std::unique_ptr<ThisType> copy() const = 0;
+  template <class... Args>
+  LocalBinaryToUnaryIntersectionIntegrand<I, t_r, t_rC, TR, F, a_r, a_rC, AR> with_ansatz(Args&&... args) const
+  {
+    return LocalBinaryToUnaryIntersectionIntegrand<I, t_r, t_rC, TR, F, a_r, a_rC, AR>(*this,
+                                                                                       std::forward<Args>(args)...);
+  }
+
+  virtual std::unique_ptr<ThisType> copy_as_binary_intersection_integrand() const = 0;
+
+  LocalBinaryIntersectionIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR> operator+(const ThisType& other) const
+  {
+    return LocalBinaryIntersectionIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR>(*this, other);
+  }
+
+  /**
+   * Flag to document which element the bases are expected to be bound to.
+   */
+  virtual bool inside() const
+  {
+    return true;
+  }
 
   /**
    * Returns the polynomial order of the integrand, given the bases.
    *
    * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
    **/
-  virtual int order(const LocalBasisType& inside_basis,
-                    const LocalBasisType& outside_basis,
+  virtual int order(const LocalTestBasisType& test_basis,
+                    const LocalAnsatzBasisType& ansatz_basis,
                     const XT::Common::Parameter& param = {}) const = 0;
 
   /**
@@ -272,9 +503,9 @@ public:
    *
    * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
    **/
-  virtual void evaluate(const LocalBasisType& inside_basis,
-                        const LocalBasisType& outside_basis,
-                        const DomainType& point_in_reference_element,
+  virtual void evaluate(const LocalTestBasisType& test_basis,
+                        const LocalAnsatzBasisType& ansatz_basis,
+                        const DomainType& point_in_reference_intersection,
                         DynamicMatrix<F>& result,
                         const XT::Common::Parameter& param = {}) const = 0;
 
@@ -283,18 +514,97 @@ public:
    *
    * \note Will throw Exceptions::not_bound_to_an_element_yet error if not bound yet!
    **/
-  virtual DynamicMatrix<F> evaluate(const LocalBasisType& inside_basis,
-                                    const LocalBasisType& outside_basis,
-                                    const DomainType& point_in_reference_element,
+  virtual DynamicMatrix<F> evaluate(const LocalTestBasisType& test_basis,
+                                    const LocalAnsatzBasisType& ansatz_basis,
+                                    const DomainType& point_in_reference_intersection,
                                     const XT::Common::Parameter& param = {}) const
   {
-    DynamicMatrix<F> result(inside_basis.size(param), outside_basis.size(param), 0);
-    evaluate(inside_basis, inside_basis, point_in_reference_element, result, param);
+    DynamicMatrix<F> result(test_basis.size(param), ansatz_basis.size(param), 0);
+    evaluate(test_basis, ansatz_basis, point_in_reference_intersection, result, param);
     return result;
   }
+
+protected:
+  void ensure_size_and_clear_results(const LocalTestBasisType& test_basis,
+                                     const LocalAnsatzBasisType& ansatz_basis,
+                                     DynamicMatrix<F>& result,
+                                     const XT::Common::Parameter& param) const
+  {
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    if (result.rows() < rows || result.cols() < cols)
+      result.resize(rows, cols);
+    result *= 0;
+  } // ... ensure_size_and_clear_results(...)
 }; // class LocalBinaryIntersectionIntegrandInterface
 
 
+/**
+ * \attention We do not handle the case when the parametric nature of the integrand as a unary one differs from the
+ *            integrand as a binary one!
+ */
+template <class I,
+          size_t t_r = 1,
+          size_t t_rC = 1,
+          class TR = double,
+          class F_ = double,
+          size_t a_r = t_r,
+          size_t a_rC = t_rC,
+          class AR = TR>
+class LocalUnaryAndBinaryIntersectionIntegrandInterface
+  : public LocalUnaryIntersectionIntegrandInterface<I, t_r, t_rC, TR, F_>
+  , public LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>
+{
+
+  using ThisType = LocalUnaryAndBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>;
+
+protected:
+  using UnaryBaseType = LocalUnaryIntersectionIntegrandInterface<I, t_r, t_rC, TR, F_>;
+  using BinaryBaseType = LocalBinaryIntersectionIntegrandInterface<I, t_r, t_rC, TR, F_, a_r, a_rC, AR>;
+
+public:
+  /// \name Members and typedefs required for disambiuation.
+  /// \{
+
+  using typename UnaryBaseType::D;
+  using typename UnaryBaseType::DomainType;
+  using typename UnaryBaseType::E;
+  using typename UnaryBaseType::ElementType;
+  using typename UnaryBaseType::F;
+  using typename UnaryBaseType::IntersectionType;
+  using UnaryBaseType::d;
+
+  /// \}
+
+  LocalUnaryAndBinaryIntersectionIntegrandInterface(const XT::Common::ParameterType& param_type = {})
+    : UnaryBaseType(param_type)
+    , BinaryBaseType(param_type)
+  {}
+
+  virtual ~LocalUnaryAndBinaryIntersectionIntegrandInterface() = default;
+
+  virtual std::unique_ptr<ThisType> copy_as_unary_and_binary_intersection_integrand() const = 0;
+
+  using UnaryBaseType::operator+;
+  using BinaryBaseType::operator+;
+
+  LocalUnaryAndBinaryIntersectionIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR> operator+(const ThisType& other) const
+  {
+    return LocalUnaryAndBinaryIntersectionIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR>(*this, other);
+  }
+
+  /// \name Methods required for disambiuation.
+  /// \{
+
+  using UnaryBaseType::bind;
+  using UnaryBaseType::intersection;
+  using UnaryBaseType::parameter_type;
+  using UnaryBaseType::parse_parameter;
+
+  /// \}
+}; // class LocalUnaryAndBinaryIntersectionIntegrandInterface
+
+
 template <class Intersection,
           size_t test_range_dim = 1,
           size_t test_range_dim_cols = 1,
@@ -346,7 +656,7 @@ public:
 
   virtual ~LocalQuaternaryIntersectionIntegrandInterface() = default;
 
-  virtual std::unique_ptr<ThisType> copy() const = 0;
+  virtual std::unique_ptr<ThisType> copy_as_quaternary_intersection_integrand() const = 0;
 
   LocalQuaternaryIntersectionIntegrandSum<I, t_r, t_rC, TR, F, a_r, a_rC, AR> operator+(const ThisType& other) const
   {
diff --git a/dune/gdt/local/integrands/ipdg.hh b/dune/gdt/local/integrands/ipdg.hh
index c65453f9b..ac8e538e9 100644
--- a/dune/gdt/local/integrands/ipdg.hh
+++ b/dune/gdt/local/integrands/ipdg.hh
@@ -93,7 +93,7 @@ public:
 
   InnerPenalty(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -193,9 +193,9 @@ private:
 
 
 template <class I>
-class BoundaryPenalty : public LocalQuaternaryIntersectionIntegrandInterface<I>
+class BoundaryPenalty : public LocalBinaryIntersectionIntegrandInterface<I>
 {
-  using BaseType = LocalQuaternaryIntersectionIntegrandInterface<I>;
+  using BaseType = LocalBinaryIntersectionIntegrandInterface<I>;
   using ThisType = BoundaryPenalty<I>;
 
 public:
@@ -228,7 +228,7 @@ public:
 
   BoundaryPenalty(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -240,54 +240,48 @@ protected:
   }
 
 public:
-  int order(const LocalTestBasisType& test_basis_inside,
-            const LocalAnsatzBasisType& ansatz_basis_inside,
-            const LocalTestBasisType& /*test_basis_outside*/,
-            const LocalAnsatzBasisType& /*ansatz_basis_outside*/,
+  bool inside() const override final
+  {
+    return true; // We expect the bases to be bound to the inside (see evaluate).
+  }
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
             const XT::Common::Parameter& param = {}) const override final
   {
-    return local_weight_->order(param) + test_basis_inside.order(param) + ansatz_basis_inside.order(param);
+    return local_weight_->order(param) + test_basis.order(param) + ansatz_basis.order(param);
   }
 
-  void evaluate(const LocalTestBasisType& test_basis_inside,
-                const LocalAnsatzBasisType& ansatz_basis_inside,
-                const LocalTestBasisType& test_basis_outside,
-                const LocalAnsatzBasisType& ansatz_basis_outside,
+  using BaseType::evaluate;
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
                 const DomainType& point_in_reference_intersection,
-                DynamicMatrix<F>& result_in_in,
-                DynamicMatrix<F>& result_in_out,
-                DynamicMatrix<F>& result_out_in,
-                DynamicMatrix<F>& result_out_out,
+                DynamicMatrix<F>& result,
                 const XT::Common::Parameter& param = {}) const override final
   {
     // Prepare sotrage, ...
-    this->ensure_size_and_clear_results(test_basis_inside,
-                                        ansatz_basis_inside,
-                                        test_basis_outside,
-                                        ansatz_basis_outside,
-                                        result_in_in,
-                                        result_in_out,
-                                        result_out_in,
-                                        result_out_out,
-                                        param);
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    if (result.rows() < rows || result.cols() < cols)
+      result.resize(rows, cols);
+    result *= 0;
     // evaluate ...
     const auto point_in_inside_reference_element =
         this->intersection().geometryInInside().global(point_in_reference_intersection);
     const auto normal = this->intersection().unitOuterNormal(point_in_reference_intersection);
     // ... basis functions ...
-    test_basis_inside.evaluate(point_in_inside_reference_element, test_basis_values_, param);
-    ansatz_basis_inside.evaluate(point_in_inside_reference_element, ansatz_basis_values_, param);
+    test_basis.evaluate(point_in_inside_reference_element, test_basis_values_, param);
+    ansatz_basis.evaluate(point_in_inside_reference_element, ansatz_basis_values_, param);
     // ... and data functions, ....
     const auto weight = local_weight_->evaluate(point_in_inside_reference_element, param);
     // compute the weighted penalty ...
     const auto h = intersection_diameter_(this->intersection());
     const auto penalty = (penalty_ * (normal * (weight * normal))) / h;
     // and finally compute integrand.
-    const size_t rows = test_basis_inside.size(param);
-    const size_t cols = ansatz_basis_inside.size(param);
     for (size_t ii = 0; ii < rows; ++ii)
       for (size_t jj = 0; jj < cols; ++jj)
-        result_in_in[ii][jj] += penalty * ansatz_basis_values_[jj] * test_basis_values_[ii];
+        result[ii][jj] += penalty * ansatz_basis_values_[jj] * test_basis_values_[ii];
   } // ... evaluate(...)
 
 private:
diff --git a/dune/gdt/local/integrands/laplace-ipdg.hh b/dune/gdt/local/integrands/laplace-ipdg.hh
index 23e258511..7cd2bb424 100644
--- a/dune/gdt/local/integrands/laplace-ipdg.hh
+++ b/dune/gdt/local/integrands/laplace-ipdg.hh
@@ -70,7 +70,7 @@ public:
 
   InnerCoupling(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -211,9 +211,9 @@ private:
  *       * symmetry_prefactor = 1 && weight_function = diffusion => SWIPDG
  */
 template <class I>
-class DirichletCoupling : public LocalQuaternaryIntersectionIntegrandInterface<I>
+class DirichletCoupling : public LocalUnaryAndBinaryIntersectionIntegrandInterface<I>
 {
-  using BaseType = LocalQuaternaryIntersectionIntegrandInterface<I>;
+  using BaseType = LocalUnaryAndBinaryIntersectionIntegrandInterface<I>;
   using ThisType = DirichletCoupling<I>;
 
 public:
@@ -225,23 +225,42 @@ public:
   using typename BaseType::LocalAnsatzBasisType;
   using typename BaseType::LocalTestBasisType;
 
-  DirichletCoupling(const double& symmetry_prefactor, XT::Functions::GridFunction<E, d, d> diffusion)
-    : BaseType(diffusion.parameter_type())
+  /**
+   * \note dirichlet_data is only required if used as a unary integrand, i.e. for the right hand side
+   */
+  DirichletCoupling(const double& symmetry_prefactor,
+                    XT::Functions::GridFunction<E, d, d> diffusion,
+                    XT::Functions::GridFunction<E> dirichlet_data = 0.)
+    : BaseType(diffusion.parameter_type() + dirichlet_data.parameter_type())
     , symmetry_prefactor_(symmetry_prefactor)
     , diffusion_(diffusion)
+    , dirichlet_data_(dirichlet_data)
     , local_diffusion_(diffusion_.local_function())
+    , local_dirichlet_data_(dirichlet_data_.local_function())
   {}
 
   DirichletCoupling(const ThisType& other)
     : BaseType(other.parameter_type())
     , symmetry_prefactor_(other.symmetry_prefactor_)
     , diffusion_(other.diffusion_)
+    , dirichlet_data_(other.dirichlet_data_)
     , local_diffusion_(diffusion_.local_function())
+    , local_dirichlet_data_(dirichlet_data_.local_function())
   {}
 
   DirichletCoupling(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<typename BaseType::UnaryBaseType> copy_as_unary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+  std::unique_ptr<typename BaseType::BinaryBaseType> copy_as_binary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+  std::unique_ptr<BaseType> copy_as_unary_and_binary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -251,70 +270,102 @@ protected:
   {
     const auto inside_element = intersection.inside();
     local_diffusion_->bind(inside_element);
+    local_dirichlet_data_->bind(inside_element);
   }
 
 public:
-  int order(const LocalTestBasisType& test_basis_inside,
-            const LocalAnsatzBasisType& ansatz_basis_inside,
-            const LocalTestBasisType& /*test_basis_outside*/,
-            const LocalAnsatzBasisType& /*ansatz_basis_outside*/,
+  bool inside() const override final
+  {
+    return true; // We expect the bases to be bound to the inside (see evaluate and post_bind).
+  }
+
+  /// \name Required by LocalUnaryIntersectionIntegrandInterface.
+  /// \{
+
+  int order(const LocalTestBasisType& test_basis, const XT::Common::Parameter& param = {}) const override final
+  {
+    return local_dirichlet_data_->order(param) + local_diffusion_->order(param)
+           + std::max(test_basis.order(param) - 1, 0);
+  }
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicVector<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    using Base = typename BaseType::UnaryBaseType;
+    // Prepare sotrage, ...
+    Base::ensure_size_and_clear_results(test_basis, result, param);
+    // evaluate ...
+    const auto point_in_inside_reference_element =
+        Base::intersection().geometryInInside().global(point_in_reference_intersection);
+    const auto normal = Base::intersection().unitOuterNormal(point_in_reference_intersection);
+    // ... basis functions and ...
+    test_basis.jacobians(point_in_inside_reference_element, test_basis_grads_, param);
+    // ... data functions, ...
+    const auto diffusion = local_diffusion_->evaluate(point_in_inside_reference_element, param);
+    const auto dirichlet_data = local_dirichlet_data_->evaluate(point_in_inside_reference_element, param);
+    // ... and finally compute the integrand.
+    const size_t size = test_basis.size(param);
+    for (size_t ii = 0; ii < size; ++ii)
+      result[ii] += -1.0 * symmetry_prefactor_ * dirichlet_data * ((diffusion * test_basis_grads_[ii][0]) * normal);
+  } // ... evaluate(...)
+
+  /// \}
+  /// \name Required by LocalBinaryIntersectionIntegrandInterface.
+  /// \{
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
             const XT::Common::Parameter& param = {}) const override final
   {
-    return local_diffusion_->order(param) + test_basis_inside.order(param) + ansatz_basis_inside.order(param);
+    return local_diffusion_->order(param) + test_basis.order(param) + ansatz_basis.order(param);
   }
 
-  void evaluate(const LocalTestBasisType& test_basis_inside,
-                const LocalAnsatzBasisType& ansatz_basis_inside,
-                const LocalTestBasisType& test_basis_outside,
-                const LocalAnsatzBasisType& ansatz_basis_outside,
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
                 const DomainType& point_in_reference_intersection,
-                DynamicMatrix<F>& result_in_in,
-                DynamicMatrix<F>& result_in_out,
-                DynamicMatrix<F>& result_out_in,
-                DynamicMatrix<F>& result_out_out,
+                DynamicMatrix<F>& result,
                 const XT::Common::Parameter& param = {}) const override final
   {
+    using Base = typename BaseType::BinaryBaseType;
     // Prepare sotrage, ...
-    this->ensure_size_and_clear_results(test_basis_inside,
-                                        ansatz_basis_inside,
-                                        test_basis_outside,
-                                        ansatz_basis_outside,
-                                        result_in_in,
-                                        result_in_out,
-                                        result_out_in,
-                                        result_out_out,
-                                        param);
+    Base::ensure_size_and_clear_results(test_basis, ansatz_basis, result, param);
     // evaluate ...
     const auto point_in_inside_reference_element =
-        this->intersection().geometryInInside().global(point_in_reference_intersection);
-    const auto normal = this->intersection().unitOuterNormal(point_in_reference_intersection);
+        Base::intersection().geometryInInside().global(point_in_reference_intersection);
+    const auto normal = Base::intersection().unitOuterNormal(point_in_reference_intersection);
     // ... basis functions and ...
-    test_basis_inside.evaluate(point_in_inside_reference_element, test_basis_values_, param);
-    test_basis_inside.jacobians(point_in_inside_reference_element, test_basis_grads_, param);
-    ansatz_basis_inside.evaluate(point_in_inside_reference_element, ansatz_basis_values_, param);
-    ansatz_basis_inside.jacobians(point_in_inside_reference_element, ansatz_basis_grads_, param);
+    test_basis.evaluate(point_in_inside_reference_element, test_basis_values_, param);
+    test_basis.jacobians(point_in_inside_reference_element, test_basis_grads_, param);
+    ansatz_basis.evaluate(point_in_inside_reference_element, ansatz_basis_values_, param);
+    ansatz_basis.jacobians(point_in_inside_reference_element, ansatz_basis_grads_, param);
     // ... data functions, ...
     const auto diffusion = local_diffusion_->evaluate(point_in_inside_reference_element, param);
     // ... and finally compute the integrand.
-    const size_t rows = test_basis_inside.size(param);
-    const size_t cols = ansatz_basis_inside.size(param);
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
     for (size_t ii = 0; ii < rows; ++ii)
       for (size_t jj = 0; jj < cols; ++jj) {
-        result_in_in[ii][jj] += -1.0 * ((diffusion * ansatz_basis_grads_[jj][0]) * normal) * test_basis_values_[ii];
-        result_in_in[ii][jj] +=
+        result[ii][jj] += -1.0 * ((diffusion * ansatz_basis_grads_[jj][0]) * normal) * test_basis_values_[ii];
+        result[ii][jj] +=
             -1.0 * symmetry_prefactor_ * ansatz_basis_values_[jj] * ((diffusion * test_basis_grads_[ii][0]) * normal);
       }
   } // ... evaluate(...)
 
+  /// \}
+
 private:
   const double symmetry_prefactor_;
   XT::Functions::GridFunction<E, d, d> diffusion_;
+  XT::Functions::GridFunction<E> dirichlet_data_;
   std::unique_ptr<typename XT::Functions::GridFunctionInterface<E, d, d>::LocalFunctionType> local_diffusion_;
+  std::unique_ptr<typename XT::Functions::GridFunctionInterface<E>::LocalFunctionType> local_dirichlet_data_;
   mutable std::vector<typename LocalTestBasisType::RangeType> test_basis_values_;
   mutable std::vector<typename LocalTestBasisType::DerivativeRangeType> test_basis_grads_;
   mutable std::vector<typename LocalAnsatzBasisType::RangeType> ansatz_basis_values_;
   mutable std::vector<typename LocalAnsatzBasisType::DerivativeRangeType> ansatz_basis_grads_;
-}; // DirichletCoupling
+}; // class DirichletCoupling
 
 
 } // namespace LocalLaplaceIPDGIntegrands
diff --git a/dune/gdt/local/integrands/laplace.hh b/dune/gdt/local/integrands/laplace.hh
index 95b480b0c..54a74dff6 100644
--- a/dune/gdt/local/integrands/laplace.hh
+++ b/dune/gdt/local/integrands/laplace.hh
@@ -52,7 +52,7 @@ public:
 
   LocalLaplaceIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
diff --git a/dune/gdt/local/integrands/product.hh b/dune/gdt/local/integrands/product.hh
index 78a561aec..98ac3bd5e 100644
--- a/dune/gdt/local/integrands/product.hh
+++ b/dune/gdt/local/integrands/product.hh
@@ -31,8 +31,6 @@ namespace GDT {
  * Given an inducing function f, computes `f(x) * phi(x) * psi(x)` for all combinations of phi and psi in the bases.
  *
  * \note Note that f can also be given as a scalar value or omitted.
- *
- * \sa local_binary_to_unary_element_integrand
  */
 template <class E, size_t r = 1, class TR = double, class F = double, class AR = TR>
 class LocalElementProductIntegrand : public LocalBinaryElementIntegrandInterface<E, r, 1, TR, F, r, 1, AR>
@@ -61,7 +59,7 @@ public:
 
   LocalElementProductIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_binary_element_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -118,9 +116,10 @@ private:
  * \note Note that f can also be given as a scalar value or omitted.
  */
 template <class I, size_t r = 1, class TR = double, class F = double, class AR = TR>
-class LocalIntersectionProductIntegrand : public LocalQuaternaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>
+class LocalCouplingIntersectionProductIntegrand
+  : public LocalQuaternaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>
 {
-  using ThisType = LocalIntersectionProductIntegrand<I, r, TR, F, AR>;
+  using ThisType = LocalCouplingIntersectionProductIntegrand<I, r, TR, F, AR>;
   using BaseType = LocalQuaternaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>;
 
 public:
@@ -133,23 +132,23 @@ public:
 
   using GridFunctionType = XT::Functions::GridFunctionInterface<E, r, r, F>;
 
-  LocalIntersectionProductIntegrand(XT::Functions::GridFunction<E, r, r, F> weight = {1.})
+  LocalCouplingIntersectionProductIntegrand(XT::Functions::GridFunction<E, r, r, F> weight = {1.})
     : BaseType()
     , weight_(weight)
     , local_weight_in_(weight_.local_function())
     , local_weight_out_(weight_.local_function())
   {}
 
-  LocalIntersectionProductIntegrand(const ThisType& other)
+  LocalCouplingIntersectionProductIntegrand(const ThisType& other)
     : BaseType(other.parameter_type())
     , weight_(other.weight_)
     , local_weight_in_(weight_.local_function())
     , local_weight_out_(weight_.local_function())
   {}
 
-  LocalIntersectionProductIntegrand(ThisType&& source) = default;
+  LocalCouplingIntersectionProductIntegrand(ThisType&& source) = default;
 
-  std::unique_ptr<BaseType> copy() const override final
+  std::unique_ptr<BaseType> copy_as_quaternary_intersection_integrand() const override final
   {
     return std::make_unique<ThisType>(*this);
   }
@@ -238,6 +237,115 @@ private:
   mutable std::vector<typename LocalTestBasisType::RangeType> test_basis_out_values_;
   mutable std::vector<typename LocalAnsatzBasisType::RangeType> ansatz_basis_in_values_;
   mutable std::vector<typename LocalAnsatzBasisType::RangeType> ansatz_basis_out_values_;
+}; // class LocalCouplingIntersectionProductIntegrand
+
+
+/**
+ * Given an inducing function f, computes `<f> * phi * psi` for all combinations of phi and psi in the bases, where
+ * `<f>` denotes the average of f evaluated on the inside and evaluated on the outside.
+ *
+ * \note Note that f can also be given as a scalar value or omitted.
+ */
+template <class I, size_t r = 1, class TR = double, class F = double, class AR = TR>
+class LocalIntersectionProductIntegrand : public LocalBinaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>
+{
+  using ThisType = LocalIntersectionProductIntegrand<I, r, TR, F, AR>;
+  using BaseType = LocalBinaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>;
+
+public:
+  using BaseType::d;
+  using typename BaseType::DomainType;
+  using typename BaseType::E;
+  using typename BaseType::IntersectionType;
+  using typename BaseType::LocalAnsatzBasisType;
+  using typename BaseType::LocalTestBasisType;
+
+  using GridFunctionType = XT::Functions::GridFunctionInterface<E, r, r, F>;
+
+  LocalIntersectionProductIntegrand(XT::Functions::GridFunction<E, r, r, F> weight = {1.}, bool use_inside_bases = true)
+    : BaseType()
+    , weight_(weight)
+    , local_weight_(weight_.local_function())
+    , inside_(use_inside_bases)
+  {}
+
+  LocalIntersectionProductIntegrand(const ThisType& other)
+    : BaseType(other.parameter_type())
+    , weight_(other.weight_)
+    , local_weight_(weight_.local_function())
+    , inside_(other.inside_)
+  {}
+
+  LocalIntersectionProductIntegrand(ThisType&& source) = default;
+
+  std::unique_ptr<BaseType> copy_as_binary_intersection_integrand() const override final
+  {
+    return std::make_unique<ThisType>(*this);
+  }
+
+protected:
+  void post_bind(const IntersectionType& intersct) override final
+  {
+    if (inside_)
+      local_weight_->bind(intersct.inside());
+    else {
+      DUNE_THROW_IF(!intersct.neighbor(), Exceptions::integrand_error, "");
+      local_weight_->bind(intersct.outside());
+    }
+  } // ... post_bind(...)
+
+public:
+  bool inside() const override final
+  {
+    return inside_;
+  }
+
+  int order(const LocalTestBasisType& test_basis,
+            const LocalAnsatzBasisType& ansatz_basis,
+            const XT::Common::Parameter& param = {}) const override final
+  {
+    return local_weight_->order(param) + test_basis.order(param) + ansatz_basis.order(param);
+  }
+
+  using BaseType::evaluate;
+
+  void evaluate(const LocalTestBasisType& test_basis,
+                const LocalAnsatzBasisType& ansatz_basis,
+                const DomainType& point_in_reference_intersection,
+                DynamicMatrix<F>& result,
+                const XT::Common::Parameter& param = {}) const override final
+  {
+    // prepare sotrage
+    this->ensure_size_and_clear_results(test_basis, ansatz_basis, result, param);
+    // evaluate
+    typename XT::Functions::GridFunction<E, r, r, F>::LocalFunctionType::RangeReturnType weight;
+    if (inside_) {
+      const auto point_in_inside_reference_element =
+          this->intersection().geometryInInside().global(point_in_reference_intersection);
+      test_basis.evaluate(point_in_inside_reference_element, test_basis_values_, param);
+      ansatz_basis.evaluate(point_in_inside_reference_element, ansatz_basis_values_, param);
+      weight = local_weight_->evaluate(point_in_inside_reference_element, param);
+    } else {
+      const auto point_in_outside_reference_element =
+          this->intersection().geometryInOutside().global(point_in_reference_intersection);
+      test_basis.evaluate(point_in_outside_reference_element, test_basis_values_, param);
+      ansatz_basis.evaluate(point_in_outside_reference_element, ansatz_basis_values_, param);
+      weight = local_weight_->evaluate(point_in_outside_reference_element, param);
+    }
+    // compute integrand
+    const size_t rows = test_basis.size(param);
+    const size_t cols = ansatz_basis.size(param);
+    for (size_t ii = 0; ii < rows; ++ii)
+      for (size_t jj = 0; jj < cols; ++jj)
+        result[ii][jj] = (weight * ansatz_basis_values_[jj]) * test_basis_values_[ii];
+  } // ... evaluate(...)
+
+private:
+  XT::Functions::GridFunction<E, r, r, F> weight_;
+  std::unique_ptr<typename XT::Functions::GridFunction<E, r, r, F>::LocalFunctionType> local_weight_;
+  const bool inside_;
+  mutable std::vector<typename LocalTestBasisType::RangeType> test_basis_values_;
+  mutable std::vector<typename LocalAnsatzBasisType::RangeType> ansatz_basis_values_;
 }; // class LocalIntersectionProductIntegrand
 
 
-- 
GitLab