From 6f902b1cf0330f2a1984ec3ec1bf51920dbc9423 Mon Sep 17 00:00:00 2001
From: Felix Schindler <felix.schindler@wwu.de>
Date: Tue, 23 Feb 2021 09:44:23 +0100
Subject: [PATCH] [local.integrands] add
LocalIntersectionNormalComponentProductIntegrand
---
dune/gdt/local/integrands/product.hh | 113 +++++++++++++++++++++++++++
1 file changed, 113 insertions(+)
diff --git a/dune/gdt/local/integrands/product.hh b/dune/gdt/local/integrands/product.hh
index ab8a9fc58..fe209f629 100644
--- a/dune/gdt/local/integrands/product.hh
+++ b/dune/gdt/local/integrands/product.hh
@@ -386,6 +386,119 @@ private:
}; // class LocalIntersectionProductIntegrand
+/**
+ * Given an inducing vector-valued function f, computes `(f * n) * phi * psi` for all combinations of phi and psi in the bases.
+ *
+ * \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 LocalIntersectionNormalComponentProductIntegrand : public LocalBinaryIntersectionIntegrandInterface<I, r, 1, TR, F, r, 1, AR>
+{
+ using ThisType = LocalIntersectionNormalComponentProductIntegrand;
+ 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, d, 1, F>;
+
+ LocalIntersectionNormalComponentProductIntegrand(XT::Functions::GridFunction<E, d, 1, F> weight,
+ const bool use_inside_bases = true,
+ const std::string& logging_prefix = "")
+ : BaseType({},
+ logging_prefix.empty() ? "LocalIntersectionNormalComponentProductIntegrand" : logging_prefix,
+ /*logging_disabled=*/logging_prefix.empty())
+ , weight_(weight.copy_as_grid_function())
+ , local_weight_(weight_->local_function())
+ , inside_(use_inside_bases)
+ {}
+
+ LocalIntersectionNormalComponentProductIntegrand(const ThisType& other)
+ : BaseType(other)
+ , weight_(other.weight_->copy_as_grid_function())
+ , local_weight_(weight_->local_function())
+ , inside_(other.inside_)
+ {}
+
+ LocalIntersectionNormalComponentProductIntegrand(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
+ F weight;
+ const auto normal = this->intersection().unitOuterNormal(point_in_reference_intersection);
+ 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) * normal);
+ } 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) * normal);
+ }
+ // 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:
+ const std::unique_ptr<XT::Functions::GridFunctionInterface<E, d, 1, F>> weight_;
+ std::unique_ptr<typename XT::Functions::GridFunctionInterface<E, d, 1, 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 LocalIntersectionNormalComponentProductIntegrand
+
+
template <class E_or_I, size_t r = 1, class TR = double, class F = double, class AR = TR>
class LocalProductIntegrand
: public std::conditional_t<XT::Grid::is_entity<E_or_I>::value,
--
GitLab