dunefunctionsoperatorassembler.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_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
 
#ifndef DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
#warning This header is deprecated and will be removed after 2.11. Use the global assemblers from dune/assembler/defaultglobalassembler.hh instead.
#endif
 
#include <cstddef>
#include <type_traits>
 
#include <dune/common/deprecated.hh>
#include <dune/common/std/type_traits.hh>
 
#include <dune/istl/matrix.hh>
 
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include <dune/fufem/parallel/parallelalgorithm.hh>
#include <dune/fufem/backends/istlmatrixbackend.hh>
 
 
namespace Dune {
namespace Fufem {
 
 
 
template <class TestBasis, class AnsatzBasis>
class
[[deprecated("This class will be removed after 2.11. Use Dune::Assembler::Assembler together with Dune::Fufem::Forms based local assemblers instead.")]]
DuneFunctionsOperatorAssembler
{
  using GridView = typename TestBasis::GridView;
 
  // Zero-initialize local entries contained in the tree.
  // In case of a subspace basis, this does not touch
  // entries not contained in the subspace.
  template<class TestLocalView, class AnsatzLocalView, class LocalMatrix>
  static void zeroInitializeLocal(const TestLocalView& testLocalView, const AnsatzLocalView& ansatzLocalView, LocalMatrix& localMatrix)
  {
    for (size_t i=0; i<testLocalView.tree().size(); ++i)
    {
      auto localRow = testLocalView.tree().localIndex(i);
      for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
      {
        auto localCol = ansatzLocalView.tree().localIndex(j);
        localMatrix[localRow][localCol] = 0;
      }
    }
  }
 
  // Distribute local matrix pattern to global pattern builder
  // In case of a subspace basis, this does not touch
  // entries not contained in the subspace.
  template<class TestLocalView, class AnsatzLocalView, class PatternBuilder>
  static void distributeLocalPattern(const TestLocalView& testLocalView, const AnsatzLocalView& ansatzLocalView, PatternBuilder& patternBuilder)
  {
    for (size_t i=0; i<testLocalView.tree().size(); ++i)
    {
      auto localRow = testLocalView.tree().localIndex(i);
      auto row = testLocalView.index(localRow);
      for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
      {
        auto localCol = ansatzLocalView.tree().localIndex(j);
        auto col = ansatzLocalView.index(localCol);
        patternBuilder.insertEntry(row,col);
      }
    }
  }
 
  // Distribute local matrix entries to global matrix
  // In case of a subspace basis, this does not touch
  // entries not contained in the subspace.
  template<class TestLocalView, class AnsatzLocalView, class LocalMatrix, class MatrixBackend>
  static void distributeLocalEntries(const TestLocalView& testLocalView, const AnsatzLocalView& ansatzLocalView, const LocalMatrix& localMatrix, MatrixBackend& matrixBackend)
  {
    for (size_t i=0; i<testLocalView.tree().size(); ++i)
    {
      auto localRow = testLocalView.tree().localIndex(i);
      auto row = testLocalView.index(localRow);
      for (size_t j=0; j<ansatzLocalView.tree().size(); ++j)
      {
        auto localCol = ansatzLocalView.tree().localIndex(j);
        auto col = ansatzLocalView.index(localCol);
        matrixBackend(row,col) += localMatrix[localRow][localCol];
      }
    }
  }
 
  // SFINAE expression check for preprocess() method of local assembler
  template <class LocalAssembler>
  using LocalAssemblerPreprocess = decltype(std::declval<LocalAssembler>().preprocess(std::declval<typename TestBasis::LocalView>(), std::declval<typename AnsatzBasis::LocalView>()));
 
public:
  [[deprecated("This class will be removed after 2.11. Use Dune::Assembler::Assembler together with Dune::Fufem::Forms based local assemblers instead.")]]
  DuneFunctionsOperatorAssembler(const TestBasis& tBasis, const AnsatzBasis& aBasis) :
      testBasis_(tBasis),
      ansatzBasis_(aBasis)
  {}
 
 
 
  template <class PatternBuilder>
  void assembleBulkPattern(PatternBuilder& patternBuilder) const
  {
    patternBuilder.resize(testBasis_, ansatzBasis_);
 
    // create two localViews but use only one if bases are the same
    auto ansatzLocalView = ansatzBasis_.localView();
    auto separateTestLocalView = testBasis_.localView();
    auto& testLocalView = selectTestLocalView(ansatzLocalView, separateTestLocalView);
 
    for (const auto& element : elements(testBasis_.gridView()))
    {
      // bind the localViews to the element
      bind(ansatzLocalView, testLocalView, element);
 
      distributeLocalPattern(testLocalView, ansatzLocalView, patternBuilder);
    }
  }
 
  template <class PatternBuilder, class ElementPartition>
  void assembleBulkPattern(PatternBuilder& patternBuilder, const ElementPartition& elementPartition, std::size_t threadCount) const
  {
    static_assert(Dune::Capabilities::viewThreadSafe<typename GridView::Grid>::v, "Trying to use thread parallel assembler but grid is not viewThreadSafe.");
 
    patternBuilder.resize(testBasis_, ansatzBasis_);
 
    parallelAlgorithm(threadCount, [&](auto parallel) mutable {
      // create two localViews but use only one if bases are the same
      auto ansatzLocalView = ansatzBasis_.localView();
      auto separateTestLocalView = testBasis_.localView();
      auto& testLocalView = selectTestLocalView(ansatzLocalView, separateTestLocalView);
 
      parallel.coloredForEach(elementPartition, [&](const auto& element) {
        bind(ansatzLocalView, testLocalView, element);
        distributeLocalPattern(testLocalView, ansatzLocalView, patternBuilder);
      });
    });
  }
 
 
  template <class PatternBuilder>
  void assembleSkeletonPattern(PatternBuilder& patternBuilder) const
  {
    patternBuilder.resize(testBasis_, ansatzBasis_);
 
    auto insideTestLocalView        = testBasis_.localView();
 
    auto insideAnsatzLocalView      = ansatzBasis_.localView();
 
    auto outsideAnsatzLocalView     = ansatzBasis_.localView();
 
    for (const auto& element : elements(testBasis_.gridView()))
    {
      insideTestLocalView.bind(element);
 
      insideAnsatzLocalView.bind(element);
 
      /* Coupling on the same element */
      distributeLocalPattern(insideTestLocalView, insideAnsatzLocalView, patternBuilder);
 
      for (const auto& is : intersections(testBasis_.gridView(), element))
      {
        /* Elements on the boundary have been treated above, iterate along the inner edges */
        if (is.neighbor())
        {
          // Current outside element
          const auto& outsideElement = is.outside();
 
          // Bind the outer parts to the outer element
          outsideAnsatzLocalView.bind(outsideElement);
 
          // We assume that all basis functions of the inner element couple with all basis functions from the outer one
          distributeLocalPattern(insideTestLocalView, outsideAnsatzLocalView, patternBuilder);
        }
      }
    }
  }
 
 
 
  template <class Matrix, class LocalAssembler>
  void assembleBulkEntries(Matrix&& matrix, LocalAssembler&& localAssembler) const
  {
DUNE_NO_DEPRECATED_BEGIN
    auto&& matrixBackend = toMatrixBackend<TestBasis, AnsatzBasis>(matrix);
DUNE_NO_DEPRECATED_END
 
    // create two localViews but use only one if bases are the same
    auto ansatzLocalView = ansatzBasis_.localView();
    auto separateTestLocalView = testBasis_.localView();
    auto& testLocalView = selectTestLocalView(ansatzLocalView, separateTestLocalView);
 
    if constexpr(Dune::Std::is_detected_v<LocalAssemblerPreprocess, LocalAssembler>)
      localAssembler.preprocess(testLocalView, ansatzLocalView);
 
    using Field = std::decay_t<decltype(matrixBackend(testLocalView.index(0), ansatzLocalView.index(0)))>;
    using LocalMatrix = Dune::Matrix<Field>;
 
    auto localMatrix = LocalMatrix();
 
    localMatrix.setSize(testLocalView.maxSize(), ansatzLocalView.maxSize());
 
    for (const auto& element : elements(testBasis_.gridView()))
    {
      // bind the localViews to the element
      bind(ansatzLocalView, testLocalView, element);
 
      zeroInitializeLocal(testLocalView, ansatzLocalView, localMatrix);
 
      localAssembler(element, localMatrix, testLocalView, ansatzLocalView);
 
      // Add element stiffness matrix onto the global stiffness matrix
      distributeLocalEntries(testLocalView, ansatzLocalView, localMatrix, matrixBackend);
    }
  }
 
 
  template <class Matrix, class LocalAssembler, class ElementPartition>
  void assembleBulkEntries(Matrix&& matrix, LocalAssembler&& localAssembler, const ElementPartition& elementPartition, std::size_t threadCount) const
  {
    static_assert(Dune::Capabilities::viewThreadSafe<typename GridView::Grid>::v, "Trying to use thread parallel assembler but grid is not viewThreadSafe.");
 
DUNE_NO_DEPRECATED_BEGIN
    auto&& matrixBackend = toMatrixBackend<TestBasis, AnsatzBasis>(matrix);
DUNE_NO_DEPRECATED_END
 
    parallelAlgorithm(threadCount, [&, localAssembler](auto parallel) mutable {
      // create two localViews but use only one if bases are the same
      auto ansatzLocalView = ansatzBasis_.localView();
      auto separateTestLocalView = testBasis_.localView();
      auto& testLocalView = selectTestLocalView(ansatzLocalView, separateTestLocalView);
 
      if constexpr(Dune::Std::is_detected_v<LocalAssemblerPreprocess, std::decay_t<LocalAssembler>>)
        localAssembler.preprocess(testLocalView, ansatzLocalView);
 
      using Field = std::decay_t<decltype(matrixBackend(testLocalView.index(0), ansatzLocalView.index(0)))>;
      using LocalMatrix = Dune::Matrix<Field>;
 
      auto localMatrix = LocalMatrix();
 
      localMatrix.setSize(testLocalView.maxSize(), ansatzLocalView.maxSize());
 
      parallel.coloredForEach(elementPartition, [&](const auto& element) {
        // bind the localViews to the element
        bind(ansatzLocalView, testLocalView, element);
 
        zeroInitializeLocal(testLocalView, ansatzLocalView, localMatrix);
 
        localAssembler(element, localMatrix, testLocalView, ansatzLocalView);
 
        // Add element stiffness matrix onto the global stiffness matrix
        distributeLocalEntries(testLocalView, ansatzLocalView, localMatrix, matrixBackend);
      });
    });
  }
 
 
 
 
  /* This variant of the IntersectionOperatorAssembler that works
   * with dune-functions bases.
   *
   * Currently not supported are constrained bases and lumping.
   *
   * Note that this was written (and tested) having discontinuous Galerkin bases in mind.
   * There may be cases where things do not work as expected for standard (continuous) Galerkin spaces.
   *
   * If passed matrix does not satisfy the MatrixBackend concept,
   * it is wrapped using `istlMatrixBackend()`.
   */
  template <class Matrix, class LocalAssembler, class LocalBoundaryAssembler>
  void assembleSkeletonEntries(Matrix&& matrix, LocalAssembler&& localAssembler, LocalBoundaryAssembler&& localBoundaryAssembler) const
  {
DUNE_NO_DEPRECATED_BEGIN
    auto&& matrixBackend = toMatrixBackend<TestBasis, AnsatzBasis>(matrix);
DUNE_NO_DEPRECATED_END
 
    Dune::MultipleCodimMultipleGeomTypeMapper<GridView> faceMapper(testBasis_.gridView(), mcmgElementLayout());
 
    auto insideTestLocalView        = testBasis_.localView();
 
    auto insideAnsatzLocalView      = ansatzBasis_.localView();
 
    auto outsideTestLocalView       = testBasis_.localView();
 
    auto outsideAnsatzLocalView     = ansatzBasis_.localView();
 
    using Field = std::decay_t<decltype(matrixBackend(insideTestLocalView.index(0), insideAnsatzLocalView.index(0)))>;
    using LocalMatrix = Dune::Matrix<Field>;
    using MatrixContainer = Dune::Matrix<LocalMatrix>;
    auto matrixContrainer = MatrixContainer(2,2);
 
    matrixContrainer[0][0].setSize(insideTestLocalView.maxSize(),  insideAnsatzLocalView.maxSize());
    matrixContrainer[0][1].setSize(insideTestLocalView.maxSize(),  outsideAnsatzLocalView.maxSize());
    matrixContrainer[1][0].setSize(outsideTestLocalView.maxSize(), insideAnsatzLocalView.maxSize());
    matrixContrainer[1][1].setSize(outsideTestLocalView.maxSize(), outsideAnsatzLocalView.maxSize());
 
    for (const auto& element : elements(testBasis_.gridView()))
    {
      insideTestLocalView.bind(element);
 
      insideAnsatzLocalView.bind(element);
 
      for (const auto& is : intersections(testBasis_.gridView(), element))
      {
 
        /* Assemble (depending on whether we have a boundary edge or not) */
        if (is.boundary())
        {
          auto& localMatrix = matrixContrainer[0][0];
 
          zeroInitializeLocal(insideTestLocalView, insideAnsatzLocalView, localMatrix);
 
          localBoundaryAssembler(is, localMatrix, insideTestLocalView, insideAnsatzLocalView);
 
          distributeLocalEntries(insideTestLocalView, insideAnsatzLocalView, localMatrix, matrixBackend);
 
        }
        else if (is.neighbor())
        {
          /* Only handle every intersection once; we choose the case where the inner element has the lower index
           *
           * This should also work with grids where elements can have different geometries. TODO: Actually test this somewhere!*/
          const auto& outsideElement = is.outside();
          if (faceMapper.index(element) > faceMapper.index(outsideElement))
            continue;
 
 
          // Bind the outer parts to the outer element
          outsideTestLocalView.bind(outsideElement);
 
          outsideAnsatzLocalView.bind(outsideElement);
 
          zeroInitializeLocal(insideTestLocalView,  insideAnsatzLocalView,  matrixContrainer[0][0]);
          zeroInitializeLocal(insideTestLocalView,  outsideAnsatzLocalView, matrixContrainer[0][1]);
          zeroInitializeLocal(outsideTestLocalView, insideAnsatzLocalView,  matrixContrainer[1][0]);
          zeroInitializeLocal(outsideTestLocalView, outsideAnsatzLocalView, matrixContrainer[1][1]);
 
          localAssembler(is, matrixContrainer, insideTestLocalView, insideAnsatzLocalView, outsideTestLocalView, outsideAnsatzLocalView);
 
          distributeLocalEntries(insideTestLocalView,  insideAnsatzLocalView,  matrixContrainer[0][0], matrixBackend);
          distributeLocalEntries(insideTestLocalView,  outsideAnsatzLocalView, matrixContrainer[0][1], matrixBackend);
          distributeLocalEntries(outsideTestLocalView, insideAnsatzLocalView,  matrixContrainer[1][0], matrixBackend);
          distributeLocalEntries(outsideTestLocalView, outsideAnsatzLocalView, matrixContrainer[1][1], matrixBackend);
        }
      }
    }
  }
 
 
  template <class Matrix, class LocalAssembler>
  void assembleBulk(Matrix&& matrix, LocalAssembler&& localAssembler) const
  {
DUNE_NO_DEPRECATED_BEGIN
    auto&& matrixBackend = toMatrixBackend<TestBasis, AnsatzBasis>(matrix);
DUNE_NO_DEPRECATED_END
 
    auto patternBuilder = matrixBackend.patternBuilder();
    assembleBulkPattern(patternBuilder);
    patternBuilder.setupMatrix();
    matrixBackend.assign(0);
    assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler));
  }
 
  template <class Matrix, class LocalAssembler, class ElementPartition>
  void assembleBulk(Matrix&& matrix, LocalAssembler&& localAssembler, const ElementPartition& elementPartition, std::size_t threadCount) const
  {
    static_assert(Dune::Capabilities::viewThreadSafe<typename GridView::Grid>::v, "Trying to use thread parallel assembler but grid is not viewThreadSafe.");
 
DUNE_NO_DEPRECATED_BEGIN
    auto&& matrixBackend = toMatrixBackend<TestBasis, AnsatzBasis>(matrix);
DUNE_NO_DEPRECATED_END
 
    auto patternBuilder = matrixBackend.patternBuilder();
    assembleBulkPattern(patternBuilder, elementPartition, threadCount);
    patternBuilder.setupMatrix();
    matrixBackend.assign(0);
    assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler), elementPartition, threadCount);
  }
 
private:
 
template<class AnsatzLocalView, class TestLocalView>
TestLocalView& selectTestLocalView(AnsatzLocalView& ansatzLocalView, TestLocalView& testLocalView) const
{
  if constexpr (std::is_same<TestBasis,AnsatzBasis>::value)
  {
    if (&testBasis_ == &ansatzBasis_)
      return ansatzLocalView;
    else
      return testLocalView;
  }
  else
    return testLocalView;
}
 
template<class AnsatzLocalView, class TestLocalView, class E>
void bind(AnsatzLocalView& ansatzLocalView, TestLocalView& testLocalView, const E& e) const
{
  ansatzLocalView.bind(e);
  if constexpr (std::is_same<TestBasis,AnsatzBasis>::value)
    if (&testLocalView == &ansatzLocalView)
      return;
  // localViews differ: bind testLocalView too
  testLocalView.bind(e);
}
 
protected:
  const TestBasis& testBasis_;
  const AnsatzBasis& ansatzBasis_;
};
 
 
template <class TestBasis, class AnsatzBasis>
// [[deprecated("This function will be removed after 2.11. Use Dune::Assembler::Assembler together with Dune::Fufem::Forms based local assemblers instead.")]]
auto duneFunctionsOperatorAssembler(const TestBasis& testBasis, const AnsatzBasis& ansatzBasis)
{
DUNE_NO_DEPRECATED_BEGIN
  return DuneFunctionsOperatorAssembler<TestBasis, AnsatzBasis>(testBasis, ansatzBasis);
DUNE_NO_DEPRECATED_END
}
 
 
 
} // namespace Fufem
} // namespace Dune
 
 
#endif // DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH