integratedboundarybilinearform.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_INTEGRATEDBOUNDARYBILINEARFORM_HH
#define DUNE_FUFEM_FORMS_INTEGRATEDBOUNDARYBILINEARFORM_HH
 
#include <cstddef>
#include <type_traits>
#include <utility>
 
#include <dune/common/indices.hh>
#include <dune/common/rangeutilities.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 Domain, class QuadratureHint=Dune::Fufem::Forms::QuadratureHints::KeepOrder>
  class IntegratedBoundaryBilinearForm
  {
    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 FacetKey = typename MultipleQuadratureCacheManager<double, BilinearOperator::Element::dimension>::FacetKey;
 
  public:
 
    using Element = typename LocalOperator::Element;
    using Intersection = typename LocalOperator::Intersection;
 
    IntegratedBoundaryBilinearForm(const BilinearOperator& sumOperator, const Domain& domain, QuadratureHint quadratureHint)
      : sumOperator_(sumOperator)
      , sumLocalOperator_(localOperator(sumOperator_))
      , domain_(domain)
      , quadratureHint_(quadratureHint)
    {}
 
    IntegratedBoundaryBilinearForm(const BilinearOperator& sumOperator, const Domain& domain)
    requires (std::is_default_constructible_v<QuadratureHint>)
      : IntegratedBoundaryBilinearForm(sumOperator, domain, QuadratureHint{})
    {}
 
    const BilinearOperator& integrandOperator() const
    {
      return sumOperator_;
    }
 
    // Dune::Assembler interface
 
    void bindLocalViews (const TestRootLocalView& testLocalView, const AnsatzRootLocalView& ansatzLocalView)
    {
      cacheManager_.clear();
      sumLocalOperator_.registerLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperator_.registerCaches(cacheManager_.prototype());
      insidePtr_ = &testLocalView.element();
    }
 
    void bindElement (const Element& element)
    {}
 
    template <class LocalPattern>
    void assembleBoundaryIntersectionMatrixPattern(const Intersection& intersection, LocalPattern& localPattern)
    {
      if (not domain_.contains(intersection))
        return;
      localPattern.addAll();
    }
 
    template<class LocalMatrix>
    void assembleBoundaryIntersectionMatrix(const Intersection& intersection, LocalMatrix& localMatrix)
    {
      using namespace Dune::Indices;
      using namespace Dune::Fufem::Forms::Impl::Tensor;
 
      if (not domain_.contains(intersection))
        return;
 
      sumLocalOperator_.bind(intersection, *insidePtr_, *insidePtr_);
 
      auto facet = intersection.indexInInside();
 
      const auto& geometry = intersection.geometry();
 
      Impl::forEachTupleEntry(sumLocalOperator_.operators(), [&](auto& op) {
        cacheManager_[FacetKey(quadratureHint_(op.quadratureRuleKey()), facet)].invalidate();
      });
 
      // Lift rank by 2 to prepend zero indices for inside
      auto localTensor = LiftRank(TensorView(localMatrix, _2), _2);
 
      Impl::forEachTupleEntry(sumLocalOperator_.operators(), [&](auto& op) {
        auto& cacheForRule = cacheManager_[FacetKey(quadratureHint_(op.quadratureRuleKey()), facet)];
        op.bindToCaches(cacheForRule);
        const auto& quadRule = cacheForRule.rule();
        for (auto k : Dune::range(quadRule.size()))
        {
          const auto& quadPoint = quadRule[k];
          auto quadPointPositionInFacet = intersection.geometryInInside().local(quadPoint.position());
          const auto integrationWeight = quadPoint.weight() * geometry.integrationElement(quadPointPositionInFacet);
          axpy(integrationWeight, op(k), localTensor);
        }
      });
    }
 
    // Dune::Fufem interface
 
    template<class TestLocalView, class AnsatzLocalView>
    void preprocess(const TestLocalView& testLocalView, const AnsatzLocalView& ansatzLocalView)
    {
      cacheManager_.clear();
      sumLocalOperator_.registerLocalViews(testLocalView.rootLocalView(), ansatzLocalView.rootLocalView());
      sumLocalOperator_.registerCaches(cacheManager_.prototype());
      insidePtr_ = &testLocalView.element();
    }
 
    template<class LocalMatrix, class TestLocalView, class AnsatzLocalView>
    void operator()(const Element& element, LocalMatrix& localMatrix, const TestLocalView& testSubspaceLocalView, const AnsatzLocalView& ansatzSubspaceLocalView)
    {
      bool processElement = true;
      if constexpr (requires{domain_.containsFaceOf(element);})
        processElement = domain_.containsFaceOf(element);
      if (processElement)
      {
        bindElement(element);
        const auto& gridView = testSubspaceLocalView.globalBasis().rootBasis().gridView();
        for(const auto& intersection : intersections(gridView, element))
          if (intersection.boundary() && !intersection.neighbor())
            assembleBoundaryIntersectionMatrix(intersection, localMatrix);
      }
    }
 
  private:
    const BilinearOperator sumOperator_;
    mutable LocalOperator sumLocalOperator_;
    const Element* insidePtr_ = nullptr;
    MultipleQuadratureCacheManager<double, Element::dimension> cacheManager_;
    const Domain& domain_;
    QuadratureHint quadratureHint_;
  };
 
  template<class BilinearOperator, class Domain>
  IntegratedBoundaryBilinearForm(const BilinearOperator&, const Domain&) -> IntegratedBoundaryBilinearForm<BilinearOperator, Domain>;
 
 
 
  template<class BilinearOperator, class Domain, class QuadratureHint>
  struct IsLocalAssembler<IntegratedBoundaryBilinearForm<BilinearOperator, Domain, QuadratureHint>> : public std::true_type {};
 
 
 
} // namespace Dune::Fufem::Forms
 
 
#endif // DUNE_FUFEM_FORMS_INTEGRATEDBOUNDARYBILINEARFORM_HH