[local.functionals] update intersection functionals according to integrands

dune/gdt/functionals/l2.hh

@@ -67,10 +67,8 @@ public:
                            const ElementFilterType& filter = ApplyOnAllElements())
     : BaseType(assembly_grid_view, source_spc, vec)
   {
-    this->append(LocalFunctionalType(local_binary_to_unary_element_integrand(inducing_function, LocalIntegrandType(1)),
-                                     over_integrate),
-                 param,
-                 filter);
+    this->append(
+        LocalFunctionalType(LocalIntegrandType(1).with_ansatz(inducing_function), over_integrate), param, filter);
   }
 
   L2VolumeVectorFunctional(AssemblyGridView assembly_grid_view,
@@ -81,10 +79,8 @@ public:
                            const ElementFilterType& filter = ApplyOnAllElements())
     : BaseType(assembly_grid_view, source_spc)
   {
-    this->append(LocalFunctionalType(local_binary_to_unary_element_integrand(inducing_function, LocalIntegrandType(1)),
-                                     over_integrate),
-                 param,
-                 filter);
+    this->append(
+        LocalFunctionalType(LocalIntegrandType(1).with_ansatz(inducing_function), over_integrate), param, filter);
   }
 }; // class L2VolumeVectorFunctional
 

dune/gdt/functionals/vector-based.hh

@@ -158,10 +158,11 @@ public:
   using typename FunctionalBaseType::SourceVectorType;
 
   using typename WalkerBaseType::ElementType;
+  using typename WalkerBaseType::IntersectionType;
   using ElementFilterType = XT::Grid::ElementFilter<AssemblyGridViewType>;
   using IntersectionFilterType = XT::Grid::IntersectionFilter<AssemblyGridViewType>;
   using ApplyOnAllElements = XT::Grid::ApplyOn::AllElements<AssemblyGridViewType>;
-  using ApplyOnAllIntersection = XT::Grid::ApplyOn::AllIntersections<AssemblyGridViewType>;
+  using ApplyOnAllIntersections = XT::Grid::ApplyOn::AllIntersections<AssemblyGridViewType>;
 
   /**
    * \name Ctors which accept an existing vector into which to assemble.
@@ -230,7 +231,23 @@ public:
         local_functional, param, XT::Grid::ApplyOn::GenericFilteredElements<AssemblyGridViewType>(filter_lambda));
   }
 
-  // similar append for LocalIntersectionFunctionalInterface ...
+  ThisType& append(const LocalIntersectionFunctionalInterface<IntersectionType, r, rC, F, DofFieldType>& local_functional,
+                   const XT::Common::Parameter& param = {},
+                   const IntersectionFilterType& filter = ApplyOnAllIntersections())
+  {
+    LocalIntersectionFunctionalAssembler<V, AssemblyGridViewType, r, rC, F, GV> tmp(
+        this->source_space(), local_functional, VectorStorage::access(), param);
+    this->append(tmp, filter);
+    return *this;
+  }
+
+  ThisType& append(const LocalIntersectionFunctionalInterface<IntersectionType, r, rC, F, DofFieldType>& local_functional,
+                   const XT::Common::Parameter& param,
+                   std::function<bool(const AssemblyGridViewType&, const IntersectionType&)> filter_lambda)
+  {
+    return append(
+        local_functional, param, XT::Grid::ApplyOn::GenericFilteredIntersections<AssemblyGridViewType>(filter_lambda));
+  }
 
   void assemble(const bool use_tbb = false) override final
   {

dune/gdt/local/assembler/functional-assemblers.hh

@@ -96,6 +96,85 @@ private:
 }; // class LocalElementFunctionalAssembler
 
 
+template <class Vector, class GridView, size_t r, size_t rC, class R = double, class SpaceGridView = GridView>
+class LocalIntersectionFunctionalAssembler : public XT::Grid::IntersectionFunctor<GridView>
+{
+  static_assert(XT::LA::is_vector<Vector>::value, "");
+  static_assert(XT::Grid::is_view<GridView>::value, "");
+  static_assert(XT::Grid::is_view<SpaceGridView>::value, "");
+
+  using ThisType = LocalIntersectionFunctionalAssembler<Vector, GridView, r, rC, R, SpaceGridView>;
+  using BaseType = XT::Grid::IntersectionFunctor<GridView>;
+
+public:
+  using typename BaseType::ElementType;
+  using typename BaseType::IntersectionType;
+  using VectorType = Vector;
+  using FieldType = typename VectorType::ScalarType;
+  using SpaceType = SpaceInterface<SpaceGridView, r, rC, R>;
+  using LocalFunctionalType = LocalIntersectionFunctionalInterface<IntersectionType, r, rC, R, FieldType>;
+
+  LocalIntersectionFunctionalAssembler(const SpaceType& space,
+                                       const LocalFunctionalType& local_functional,
+                                       VectorType& global_vector,
+                                       const XT::Common::Parameter& param = {})
+    : space_(space)
+    , local_functional_(local_functional.copy())
+    , global_vector_(global_vector)
+    , param_(param)
+    , local_vector_(space_.mapper().max_local_size())
+    , global_indices_(space_.mapper().max_local_size())
+    , basis_(space_.basis().localize())
+  {
+    DUNE_THROW_IF(global_vector_.size() != space_.mapper().size(),
+                  XT::Common::Exceptions::shapes_do_not_match,
+                  "global_vector.size() = " << global_vector_.size() << "\n  "
+                                            << "space.mapper().size()" << space_.mapper().size());
+  }
+
+  LocalIntersectionFunctionalAssembler(const ThisType& other)
+    : BaseType()
+    , space_(other.space_)
+    , local_functional_(other.local_functional_->copy())
+    , global_vector_(other.global_vector_)
+    , param_(other.param_)
+    , local_vector_(other.local_vector_)
+    , global_indices_(other.global_indices_)
+    , basis_(space_.basis().localize())
+  {}
+
+  LocalIntersectionFunctionalAssembler(ThisType&& source) = default;
+
+  BaseType* copy() override final
+  {
+    return new ThisType(*this);
+  }
+
+  void apply_local(const IntersectionType& intersection,
+                   const ElementType& inside_element,
+                   const ElementType& outside_element) override final
+  {
+    const auto& element = local_functional_->inside() ? inside_element : outside_element;
+    // apply functional
+    basis_->bind(element);
+    local_functional_->apply(intersection, *basis_, local_vector_, param_);
+    // copy local vector to global
+    space_.mapper().global_indices(element, global_indices_);
+    for (size_t jj = 0; jj < basis_->size(param_); ++jj)
+      global_vector_.add_to_entry(global_indices_[jj], local_vector_[jj]);
+  }
+
+private:
+  const SpaceType& space_;
+  std::unique_ptr<LocalFunctionalType> local_functional_;
+  VectorType& global_vector_;
+  XT::Common::Parameter param_;
+  DynamicVector<FieldType> local_vector_;
+  DynamicVector<size_t> global_indices_;
+  mutable std::unique_ptr<typename SpaceType::GlobalBasisType::LocalizedType> basis_;
+}; // class LocalIntersectionFunctionalAssembler
+
+
 } // namespace GDT
 } // namespace Dune
 

dune/gdt/local/functionals/integrals.hh

@@ -38,7 +38,7 @@ public:
 
   LocalElementIntegralFunctional(const IntegrandType& integrand, const int over_integrate = 0)
     : BaseType(integrand.parameter_type())
-    , integrand_(integrand.copy())
+    , integrand_(integrand.copy_as_unary_element_integrand())
     , over_integrate_(over_integrate)
   {}
 
@@ -49,13 +49,13 @@ public:
                                  const XT::Common::ParameterType& param_type = {},
                                  const int over_integrate = 0)
     : BaseType(param_type)
-    , integrand_(GenericIntegrand(order_function, evaluate_function, post_bind_function).copy())
+    , integrand_(GenericIntegrand(order_function, evaluate_function, post_bind_function).copy_as_unary_element_integrand())
     , over_integrate_(over_integrate)
   {}
 
   LocalElementIntegralFunctional(const ThisType& other)
     : BaseType(other.parameter_type())
-    , integrand_(other.integrand_->copy())
+    , integrand_(other.integrand_->copy_as_unary_element_integrand())
     , over_integrate_(other.over_integrate_)
   {}
 
@@ -114,29 +114,27 @@ public:
   using typename BaseType::D;
   using typename BaseType::IntersectionType;
   using typename BaseType::LocalBasisType;
-  using IntegrandType = LocalBinaryIntersectionIntegrandInterface<I, r, rC, R>;
-  using GenericIntegrand = GenericLocalBinaryIntersectionIntegrand<I, r, rC, R>;
+  using IntegrandType = LocalUnaryIntersectionIntegrandInterface<I, r, rC, R>;
+  //  using GenericIntegrand = GenericLocalUnaryIntersectionIntegrand<I, r, rC, R>;
 
   LocalIntersectionIntegralFunctional(const IntegrandType& integrand, const int over_integrate = 0)
     : BaseType(integrand.parameter_type())
-    , integrand_(integrand.copy())
+    , integrand_(integrand.copy_as_unary_intersection_integrand())
     , over_integrate_(over_integrate)
   {}
 
-  LocalIntersectionIntegralFunctional(typename GenericIntegrand::GenericOrderFunctionType order_function,
-                                      typename GenericIntegrand::GenericEvaluateFunctionType evaluate_function,
-                                      typename GenericIntegrand::GenericPostBindFunctionType post_bind_function =
-                                          [](const I&) {},
-                                      const XT::Common::ParameterType& param_type = {},
-                                      const int over_integrate = 0)
-    : BaseType(param_type)
-    , integrand_(GenericIntegrand(order_function, evaluate_function, post_bind_function).copy())
-    , over_integrate_(over_integrate)
-  {}
+  //  LocalIntersectionIntegralFunctional(typename GenericIntegrand::GenericOrderFunctionType order_function,
+  //                                      typename GenericIntegrand::GenericEvalauteFunctionType evaluate_function,
+  //                                      const XT::Common::ParameterType& param_type = {},
+  //                                      const int over_integrate = 0)
+  //    : BaseType(param_type)
+  //    , integrand_(GenericIntegrand(order_function, evaluate_function).copy())
+  //    , over_integrate_(over_integrate)
+  //  {}
 
   LocalIntersectionIntegralFunctional(const ThisType& other)
     : BaseType(other.parameter_type())
-    , integrand_(other.integrand_->copy())
+    , integrand_(other.integrand_->copy_as_unary_intersection_integrand())
     , over_integrate_(other.over_integrate_)
   {}
 
@@ -147,44 +145,45 @@ public:
     return std::make_unique<ThisType>(*this);
   }
 
+  bool inside() const override final
+  {
+    return integrand_->inside();
+  }
+
   using BaseType::apply;
 
   void apply(const IntersectionType& intersection,
-             const LocalBasisType& inside_basis,
-             const LocalBasisType& outside_basis,
-             DynamicMatrix<F>& result,
+             const LocalBasisType& test_basis,
+             DynamicVector<F>& result,
              const XT::Common::Parameter& param = {}) const override final
   {
     // prepare integand
     integrand_->bind(intersection);
     // prepare storage
-    const size_t rows = inside_basis.size(param);
-    const size_t cols = outside_basis.size(param);
-    if (result.rows() < rows || result.cols() < cols)
-      result.resize(rows, cols);
+    const size_t size = test_basis.size(param);
+    if (result.size() < size)
+      result.resize(size);
     result *= 0;
     // loop over all quadrature points
-    const auto integrand_order = integrand_->order(inside_basis, outside_basis, param) + over_integrate_;
+    const auto integrand_order = integrand_->order(test_basis, param) + over_integrate_;
     for (const auto& quadrature_point : QuadratureRules<D, d - 1>::rule(intersection.type(), integrand_order)) {
       const auto point_in_reference_intersection = quadrature_point.position();
       // integration factors
       const auto integration_factor = intersection.geometry().integrationElement(point_in_reference_intersection);
       const auto quadrature_weight = quadrature_point.weight();
       // evaluate the integrand
-      integrand_->evaluate(inside_basis, outside_basis, point_in_reference_intersection, integrand_values_, param);
-      assert(integrand_values_.rows() >= rows && "This must not happen!");
-      assert(integrand_values_.cols() >= cols && "This must not happen!");
+      integrand_->evaluate(test_basis, point_in_reference_intersection, integrand_values_, param);
+      assert(integrand_values_.size() >= size && "This must not happen!");
       // compute integral
-      for (size_t ii = 0; ii < rows; ++ii)
-        for (size_t jj = 0; jj < cols; ++jj)
-          result[ii][jj] += integrand_values_[ii][jj] * integration_factor * quadrature_weight;
+      for (size_t ii = 0; ii < size; ++ii)
+        result[ii] += integrand_values_[ii] * integration_factor * quadrature_weight;
     } // loop over all quadrature points
   } // ... apply(...)
 
 private:
   mutable std::unique_ptr<IntegrandType> integrand_;
   const int over_integrate_;
-  mutable DynamicMatrix<F> integrand_values_;
+  mutable DynamicVector<F> integrand_values_;
 }; // class LocalIntersectionIntegralFunctional
 
 

dune/gdt/local/functionals/interfaces.hh

@@ -126,24 +126,30 @@ public:
   virtual std::unique_ptr<ThisType> copy() const = 0;
 
   /**
-   * Computes the application of this functional for each combination of basis functions.
+   * Flag to document which element the basis is expected to be bound to.
+   */
+  virtual bool inside() const
+  {
+    return true;
+  }
+
+  /**
+   * Computes the application of this functional for each basis function.
    */
   virtual void apply(const IntersectionType& intersection,
-                     const LocalBasisType& inside_basis,
-                     const LocalBasisType& outside_basis,
-                     DynamicMatrix<F>& result,
+                     const LocalBasisType& test_basis,
+                     DynamicVector<F>& result,
                      const XT::Common::Parameter& param = {}) const = 0;
 
   /**
    * This method is provided for convenience and should not be used within library code.
    */
-  virtual DynamicMatrix<F> apply(const IntersectionType& intersection,
-                                 const LocalBasisType& inside_basis,
-                                 const LocalBasisType& outside_basis,
+  virtual DynamicVector<F> apply(const IntersectionType& intersection,
+                                 const LocalBasisType& test_basis,
                                  const XT::Common::Parameter& param = {}) const
   {
-    DynamicMatrix<F> ret(inside_basis.size(param), outside_basis.size(param), 0);
-    this->apply(intersection, inside_basis, outside_basis, ret, param);
+    DynamicVector<F> ret(test_basis.size(param), 0);
+    this->apply(intersection, test_basis, ret, param);
     return ret;
   }
 }; // class LocalIntersectionFunctionalInterface