From 255b3b2a312aa5fc5b99cc4a8ddfb9c62d0dc344 Mon Sep 17 00:00:00 2001
From: Felix Schindler <felix.schindler@wwu.de>
Date: Sat, 14 Mar 2020 21:28:41 +0100
Subject: [PATCH] [local.functionals] update intersection functionals according
to integrands
---
dune/gdt/functionals/l2.hh | 12 +--
dune/gdt/functionals/vector-based.hh | 21 ++++-
.../local/assembler/functional-assemblers.hh | 79 +++++++++++++++++++
dune/gdt/local/functionals/integrals.hh | 63 ++++++++-------
dune/gdt/local/functionals/interfaces.hh | 24 +++---
5 files changed, 148 insertions(+), 51 deletions(-)
diff --git a/dune/gdt/functionals/l2.hh b/dune/gdt/functionals/l2.hh
index e71cd509b..87e6b44b8 100644
--- a/dune/gdt/functionals/l2.hh
+++ b/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
diff --git a/dune/gdt/functionals/vector-based.hh b/dune/gdt/functionals/vector-based.hh
index 0841e094f..8f6c037b8 100644
--- a/dune/gdt/functionals/vector-based.hh
+++ b/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
{
diff --git a/dune/gdt/local/assembler/functional-assemblers.hh b/dune/gdt/local/assembler/functional-assemblers.hh
index 841fb6ade..8827fed64 100644
--- a/dune/gdt/local/assembler/functional-assemblers.hh
+++ b/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
diff --git a/dune/gdt/local/functionals/integrals.hh b/dune/gdt/local/functionals/integrals.hh
index a27b88a75..dff5877ce 100644
--- a/dune/gdt/local/functionals/integrals.hh
+++ b/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
diff --git a/dune/gdt/local/functionals/interfaces.hh b/dune/gdt/local/functionals/interfaces.hh
index a9cb21a82..6424b4bda 100644
--- a/dune/gdt/local/functionals/interfaces.hh
+++ b/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
--
GitLab