integratedskeletonbilinearform.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_FORMS_INTEGRATEDSKELETONBILINEARFORM_HH
#define DUNE_FUFEM_FORMS_INTEGRATEDSKELETONBILINEARFORM_HH
 
#include <cstddef>
#include <type_traits>
#include <utility>
 
#include <dune/common/indices.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/stringutility.hh>
 
#include <dune/fufem/forms/shapefunctioncache.hh>
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/localsumassembler.hh>
#include <dune/fufem/forms/quadraturehints.hh>
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class BilinearOperator, class Patch, class QuadratureHint=Dune::Fufem::Forms::QuadratureHints::KeepOrder>
  class IntegratedSkeletonBilinearForm
  {
    using TestRootBasis = std::decay_t<decltype(std::get<0>(std::declval<BilinearOperator>().basis()))>;
    using TestRootLocalView = typename TestRootBasis::LocalView;
    using TestRootTree = typename TestRootLocalView::Tree;
 
    using AnsatzRootBasis = std::decay_t<decltype(std::get<1>(std::declval<BilinearOperator>().basis()))>;
    using AnsatzRootLocalView = typename AnsatzRootBasis::LocalView;
    using AnsatzRootTree = typename AnsatzRootLocalView::Tree;
 
    using LocalOperator = decltype(localOperator(std::declval<BilinearOperator>()));
 
    using CacheManager = typename Dune::Fufem::Forms::CacheManager<double,  BilinearOperator::Element::dimension>;
    using QuadratureRule = typename CacheManager::QuadratureRule;
    using QuadraturePoint = typename QuadratureRule::value_type;
 
  public:
    using Element = typename BilinearOperator::Element;
    using Intersection = typename TestRootBasis::GridView::Intersection;
 
    IntegratedSkeletonBilinearForm(const BilinearOperator& sumOperator, const Patch& patch, QuadratureHint quadratureHint)
      : sumOperator_(sumOperator)
      , sumLocalOperator_(localOperator(sumOperator_))
      , sumLocalOperatorFlipped_(localOperator(sumOperator_))
      , patch_(patch)
      , quadratureHint_(quadratureHint)
    {}
 
    IntegratedSkeletonBilinearForm(const BilinearOperator& sumOperator, const Patch& patch)
    requires (std::is_default_constructible_v<QuadratureHint>)
      : IntegratedSkeletonBilinearForm(sumOperator, patch, QuadratureHint{})
    {}
 
    const BilinearOperator& integrandOperator() const
    {
      return sumOperator_;
    }
 
    // Dune::Assembler interface
 
    void bindLocalViews (const TestRootLocalView& testLocalView, const AnsatzRootLocalView& ansatzLocalView)
    {
      sumLocalOperator_.registerLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperator_.registerCaches(insideCache_);
      sumLocalOperatorFlipped_.registerOutsideLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperatorFlipped_.registerOutsideCaches(insideCache_);
      insidePtr_ = &testLocalView.element();
    }
 
    void bindOutsideLocalViews (const TestRootLocalView& testLocalView, const AnsatzRootLocalView& ansatzLocalView)
    {
      sumLocalOperator_.registerOutsideLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperator_.registerOutsideCaches(outsideCache_);
      sumLocalOperatorFlipped_.registerLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperatorFlipped_.registerCaches(outsideCache_);
      outsidePtr_ = &testLocalView.element();
    }
 
    void bindElement (const Element& element)
    {}
 
    template <class LocalPatterns>
    void assembleInteriorIntersectionMatrixPattern (const Intersection& intersection, LocalPatterns& localPatterns)
    {
      if (not patch_.contains(intersection))
        return;
      localPatterns[0][0].addAll();
      localPatterns[0][1].addAll();
      localPatterns[1][0].addAll();
      localPatterns[1][1].addAll();
    }
 
    template <class LocalMatrices>
    void assembleInteriorIntersectionMatrix (const Intersection& intersection, LocalMatrices& localMatrices)
    {
      using namespace Dune::Indices;
      using namespace Dune::Fufem::Forms::Impl::Tensor;
 
      if (not patch_.contains(intersection))
        return;
 
      bool isOriented = true;
      if constexpr (requires { patch_.isOriented(intersection); })
        isOriented = patch_.isOriented(intersection);
 
      // LocalMatrices only supports localMatrices[i][j][k][l] and localMatrices[i][j](k,l).
      // The TensorView provides the unified access via localTensor(i,j,k,l).
      auto localTensor = TensorView(localMatrices, _4);
 
      if (isOriented)
      {
        sumLocalOperator_.bind(intersection, *insidePtr_, *outsidePtr_);
 
        const auto& geometry = intersection.geometry();
        const auto& geometryInInside = intersection.geometryInInside();
        const auto& geometryInOutside = intersection.geometryInOutside();
 
        // Since we will often have many boundary terms (e.g. for DG)
        // while evaluations are currently not cached across elements,
        // we use a single quadrature rule for all addends despite the
        // fact that a lower order might also be sufficient for some of them.
        auto key = quadratureHint_(sumLocalOperator_.quadratureRuleKey());
        key.setGeometryType(geometry.type());
        const auto& rule = QuadratureRuleCache<double, Element::dimension-1>::rule(key);
 
        // Notice that we cannot easily cache evaluation by identifying quadrature rules
        // on intersections across different elements. To guarantee that quadrature
        // points match up on inside and outside, we have to compute the element coordinates
        // of intersection quadrature points using the intersection. This has the effect
        // that the e.g. points on face 2 may be different on two elements because the
        // intersections are not oriented consistently. The latter problem could be solved
        // by computing the element coordinates using the reference element embedding, but
        // then inside and outside quadrature points would not match up.
        //
        // Hence we take the intersection approach (which is also robust wrt. non-conforming
        // intersections) and do not cache across different intersections so far.
        // Since the operator expects a cache, we use caches that are reset for any intersection.
        insideRule_.resize(rule.size(), QuadraturePoint(typename QuadraturePoint::Vector(), 0));
        for (auto k : Dune::range(rule.size()))
          insideRule_[k] = QuadraturePoint(geometryInInside.global(rule[k].position()), 0);
        insideCache_.setRule(insideRule_);
 
        outsideRule_.resize(rule.size(), QuadraturePoint(typename QuadraturePoint::Vector(), 0));
        for (auto k : Dune::range(rule.size()))
          outsideRule_[k] = QuadraturePoint(geometryInOutside.global(rule[k].position()), 0);
        outsideCache_.setRule(outsideRule_);
 
        sumLocalOperator_.bindToCaches(insideCache_, outsideCache_);
 
        for (auto k : Dune::range(rule.size()))
        {
          Impl::forEachTupleEntry(sumLocalOperator_.operators(), [&](auto& op) {
            const auto integrationWeight = rule[k].weight() * geometry.integrationElement(rule[k].position());
            axpy(integrationWeight, op(k), localTensor);
          });
        }
      }
      else
      {
        // If the intersection is not oriented, we need to flip the two leading indices
        auto flippedLocalTensor = FlipIndices(localTensor, _2);
 
        sumLocalOperatorFlipped_.bind(intersection, *outsidePtr_, *insidePtr_);
 
        const auto& geometry = intersection.geometry();
        const auto& geometryInInside = intersection.geometryInInside();
        const auto& geometryInOutside = intersection.geometryInOutside();
 
        // Since we will often have many boundary terms (e.g. for DG)
        // while evaluations are currently not cached across elements,
        // we use a single quadrature rule for all addends despite the
        // fact that a lower order might also be sufficient for some of them.
        auto key = quadratureHint_(sumLocalOperatorFlipped_.quadratureRuleKey());
        key.setGeometryType(geometry.type());
        const auto& rule = QuadratureRuleCache<double, Element::dimension-1>::rule(key);
 
        // Notice that we cannot easily cache evaluation by identifying quadrature rules
        // on intersections across different elements. To guarantee that quadrature
        // points match up on inside and outside, we have to compute the element coordinates
        // of intersection quadrature points using the intersection. This has the effect
        // that the e.g. points on face 2 may be different on two elements because the
        // intersections are not oriented consistently. The latter problem could be solved
        // by computing the element coordinates using the reference element embedding, but
        // then inside and outside quadrature points would not match up.
        //
        // Hence we take the intersection approach (which is also robust wrt. non-conforming
        // intersections) and do not cache across different intersections so far.
        // Since the operator expects a cache, we use caches that are reset for any intersection.
        insideRule_.resize(rule.size(), QuadraturePoint(typename QuadraturePoint::Vector(), 0));
        for (auto k : Dune::range(rule.size()))
          insideRule_[k] = QuadraturePoint(geometryInInside.global(rule[k].position()), 0);
        insideCache_.setRule(insideRule_);
 
        outsideRule_.resize(rule.size(), QuadraturePoint(typename QuadraturePoint::Vector(), 0));
        for (auto k : Dune::range(rule.size()))
          outsideRule_[k] = QuadraturePoint(geometryInOutside.global(rule[k].position()), 0);
        outsideCache_.setRule(outsideRule_);
 
        sumLocalOperatorFlipped_.bindToCaches(outsideCache_, insideCache_);
 
        for (auto k : Dune::range(rule.size()))
        {
          Impl::forEachTupleEntry(sumLocalOperatorFlipped_.operators(), [&](auto& op) {
            const auto integrationWeight = rule[k].weight() * geometry.integrationElement(rule[k].position());
            axpy(integrationWeight, op(k), flippedLocalTensor);
          });
        }
      }
 
    }
 
  private:
    const BilinearOperator sumOperator_;
    mutable LocalOperator sumLocalOperator_;
    mutable LocalOperator sumLocalOperatorFlipped_;
    const Element* insidePtr_ = nullptr;
    const Element* outsidePtr_ = nullptr;
    mutable QuadratureRule insideRule_;
    mutable QuadratureRule outsideRule_;
    CacheManager insideCache_;
    CacheManager outsideCache_;
    const Patch& patch_;
    QuadratureHint quadratureHint_;
  };
 
  template<class BilinearOperator, class Patch>
  IntegratedSkeletonBilinearForm(const BilinearOperator&, const Patch&) -> IntegratedSkeletonBilinearForm<BilinearOperator, Patch>;
 
 
 
  template<class BilinearOperator, class Patch, class QuadratureHint>
  struct IsLocalAssembler<IntegratedSkeletonBilinearForm<BilinearOperator, Patch, QuadratureHint>> : public std::true_type {};
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_INTEGRATEDSKELETONBILINEARFORM_HH