defaultglobalassembler.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE 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_ASSEMBLER_GLOBALASSEMBLER_HH
#define DUNE_ASSEMBLER_GLOBALASSEMBLER_HH
 
#include <dune/assembler/assemblerconcepts.hh>
#include <dune/assembler/localviewpair.hh>
#include <dune/assembler/backends/flatarray.hh>
#include <dune/assembler/backends/flatmatrix.hh>
 
#include <dune/common/concept.hh>
#include <dune/common/typeutilities.hh>
#include <dune/common/referencehelper.hh>
 
#include <dune/functions/common/typeerasure.hh>
#include <dune/functions/functionspacebases/concepts.hh>
 
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include <dune/assembler/sequentialexecutor.hh>
 
namespace Dune::Assembler {
 
 
template <Concepts::GlobalBasis RowBasis, Concepts::GlobalBasis ColBasis = RowBasis, class Executor = SequentialExecutor>
struct Assembler
{
  using RowLocalView = typename RowBasis::LocalView;
  using ColLocalView = typename ColBasis::LocalView;
 
  template <class LA>
  static constexpr bool isLocalElementMatrixAssembler =
    Concepts::LocalElementMatrixAssembler<LA,RowLocalView,ColLocalView>;
 
  template <class LA>
  static constexpr bool isLocalElementVectorAssembler =
    Concepts::LocalElementVectorAssembler<LA,RowLocalView>;
 
  template <class LA>
  static constexpr bool isLocalInteriorIntersectionMatrixAssembler =
    Concepts::LocalInteriorIntersectionMatrixAssembler<LA,RowLocalView,ColLocalView>;
 
  template <class LA>
  static constexpr bool isLocalInteriorIntersectionVectorAssembler =
    Concepts::LocalInteriorIntersectionVectorAssembler<LA,RowLocalView>;
 
  template <class LA>
  static constexpr bool isLocalBoundaryIntersectionMatrixAssembler =
    Concepts::LocalBoundaryIntersectionMatrixAssembler<LA,RowLocalView,ColLocalView>;
 
  template <class LA>
  static constexpr bool isLocalBoundaryIntersectionVectorAssembler =
    Concepts::LocalBoundaryIntersectionVectorAssembler<LA,RowLocalView>;
 
private:
  const RowBasis& rowBasis_;
  const ColBasis& colBasis_;
  Executor executor_;
 
public:
 
  Assembler (const RowBasis& rowBasis, const ColBasis& colBasis, Executor executor)
    : rowBasis_(rowBasis)
    , colBasis_(colBasis)
    , executor_(executor)
  {
    if constexpr (Dune::ResolveRef_t<Executor>::isParallel)
    {
      constexpr bool threadSafeViews =
        Dune::Capabilities::viewThreadSafe<typename RowBasis::GridView::Grid>::v &&
        Dune::Capabilities::viewThreadSafe<typename ColBasis::GridView::Grid>::v;
      static_assert(threadSafeViews, "Trying to use Assembler with parallel Executor, but grid view is not viewThreadSafe.");
    }
  }
 
  Assembler (const RowBasis& rowBasis, const ColBasis& colBasis)
  requires std::is_default_constructible_v<Executor>
    : Assembler(rowBasis, colBasis, Executor())
  {}
 
  Assembler (const RowBasis& rowBasis)
  requires std::same_as<RowBasis, ColBasis> && std::is_default_constructible_v<Executor>
    : Assembler(rowBasis, rowBasis, Executor())
  {}
 
  // ....
  template <Concepts::LocalMatrixAssembler<RowLocalView,ColLocalView> LocalAssembler,
            class Matrix>
  void assembleMatrixEntries (LocalAssembler& la, Matrix& matrix)
  {
    executor_([&, la](auto context) mutable
    {
      auto localViews = LocalViewPair{rowBasis_, colBasis_};
      la.bindLocalViews(get<0>(localViews), get<1>(localViews));
 
      auto outsideLocalViews = LocalViewPair{rowBasis_, colBasis_};
      if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler>)
        la.bindOutsideLocalViews(get<0>(outsideLocalViews), get<1>(outsideLocalViews));
 
      // local matrices for local assembling
      using LocalMatrix = typename Matrix::LocalMatrix;
      FlatArray<LocalMatrix,2,2> localMatrices{LocalMatrix(get<0>(localViews).maxSize(), get<1>(localViews).maxSize())};
 
      using LocalPattern = typename Matrix::LocalPattern;
      FlatArray<LocalPattern,2,2> localPattern{LocalPattern(get<0>(localViews).maxSize(), get<1>(localViews).maxSize())};
 
      auto mapper = MultipleCodimMultipleGeomTypeMapper{rowBasis_.gridView(), mcmgElementLayout()};
      context.forEach(elements(rowBasis_.gridView()), [&](const auto& e)
      {
        const auto insideIndex = mapper.index(e);
        localViews.bind(e);
 
        // initialize the matrix to zero
        localMatrices[0][0].resize(get<0>(localViews).size(), get<1>(localViews).size());
        localPattern[0][0].resize(get<0>(localViews).size(), get<1>(localViews).size());
        localMatrices[0][0] = 0;
        localPattern[0][0].clear();
 
        // Whenever any contribution was computed on this
        // element, we need to scatter the local element data.
        bool elementContributionsComputed = false;
 
        const bool hasBoundaryIntersections = e.hasBoundaryIntersections();
 
        // only bind to an element if necessary
        if (isLocalElementMatrixAssembler<LocalAssembler> ||
            isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> || (isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler> && hasBoundaryIntersections))
          la.bindElement(e);
 
        if constexpr (isLocalElementMatrixAssembler<LocalAssembler>)
        {
          la.assembleElementMatrix(e, localMatrices[0][0]);
          la.assembleElementMatrixPattern(e, localPattern[0][0]);
          elementContributionsComputed = true;
        }
 
        // Only assemble on inner or boundary intersections if the assembler
        // supports this.
        if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> ||
                      isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler>)
        {
          // We only need to traverse the intersections if the assembler
          // supports interior intersections or if there are boundary intersections.
          if (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> ||
              hasBoundaryIntersections)
          {
            for (auto const& is : intersections(rowBasis_.gridView(), e))
            {
              // the intersection is on a domain boundary
              if (is.boundary() && !is.neighbor())
              {
                if constexpr (isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler>)
                {
                  la.assembleBoundaryIntersectionMatrix(is, localMatrices[0][0]);
                  la.assembleBoundaryIntersectionMatrixPattern(is, localPattern[0][0]);
                  elementContributionsComputed = true;
                }
              }
 
              // the intersection is a regular inner boundary
              if (!is.boundary() && is.neighbor())
              {
                if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler>)
                {
                  const auto outsideIndex = mapper.index(is.outside());
                  // visit each intersection only once
                  if (insideIndex < outsideIndex)
                  {
                    outsideLocalViews.bind(is.outside());
 
                    // update the sizes of the local matrices and patterns
                    localMatrices[0][1].resize(get<0>(localViews).size(), get<1>(outsideLocalViews).size());
                    localMatrices[1][0].resize(get<0>(outsideLocalViews).size(), get<1>(localViews).size());
                    localMatrices[1][1].resize(get<0>(outsideLocalViews).size(), get<1>(outsideLocalViews).size());
                    localPattern[0][1].resize(get<0>(localViews).size(), get<1>(outsideLocalViews).size());
                    localPattern[1][0].resize(get<0>(outsideLocalViews).size(), get<1>(localViews).size());
                    localPattern[1][1].resize(get<0>(outsideLocalViews).size(), get<1>(outsideLocalViews).size());
 
                    // initialize the coupling matrices to zero
                    localMatrices[0][1] = 0;
                    localMatrices[1][0] = 0;
                    localMatrices[1][1] = 0;
                    localPattern[0][1].clear();
                    localPattern[1][0].clear();
                    localPattern[1][1].clear();
 
                    la.assembleInteriorIntersectionMatrix(is, localMatrices);
                    la.assembleInteriorIntersectionMatrixPattern(is, localPattern);
                    elementContributionsComputed = true;
 
                    matrix.scatter(get<0>(localViews),get<1>(outsideLocalViews), localMatrices[0][1], localPattern[0][1]);
                    matrix.scatter(get<0>(outsideLocalViews),get<1>(localViews), localMatrices[1][0], localPattern[1][0]);
                    matrix.scatter(get<0>(outsideLocalViews),get<1>(outsideLocalViews), localMatrices[1][1], localPattern[1][1]);
                  }
                }
              }
            }
          }
        }
        if (elementContributionsComputed)
          matrix.scatter(get<0>(localViews),get<1>(localViews), localMatrices[0][0], localPattern[0][0]);
      });
    });
  }
 
  template <Concepts::LocalMatrixAssembler<RowLocalView,ColLocalView> LocalAssembler,
            class MatrixPattern>
  void assembleMatrixPattern (LocalAssembler& la, MatrixPattern& matrixPattern)
  {
    executor_([&, la](auto context) mutable
    {
      auto localViews = LocalViewPair{rowBasis_, colBasis_};
      la.bindLocalViews(get<0>(localViews), get<1>(localViews));
 
      auto outsideLocalViews = LocalViewPair{rowBasis_, colBasis_};
      if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler>)
        la.bindOutsideLocalViews(get<0>(outsideLocalViews), get<1>(outsideLocalViews));
 
      using LocalPattern = typename MatrixPattern::LocalPattern;
      FlatArray<LocalPattern,2,2> localPattern{LocalPattern(get<0>(localViews).maxSize(), get<1>(localViews).maxSize())};
 
      auto mapper = MultipleCodimMultipleGeomTypeMapper{rowBasis_.gridView(), mcmgElementLayout()};
      context.forEach(elements(rowBasis_.gridView()), [&](const auto& e)
      {
        const auto insideIndex = mapper.index(e);
        localViews.bind(e);
 
        localPattern[0][0].resize(get<0>(localViews).size(), get<1>(localViews).size());
        localPattern[0][0].clear();
 
        // Whenever any contribution was computed on this
        // element, we need to scatter the local element data.
        bool elementContributionsComputed = false;
 
        const bool hasBoundaryIntersections = e.hasBoundaryIntersections();
 
        // only bind to an element if necessary
        if (isLocalElementMatrixAssembler<LocalAssembler> ||
            isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> || (isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler> && hasBoundaryIntersections))
          la.bindElement(e);
 
        if constexpr (isLocalElementMatrixAssembler<LocalAssembler>)
        {
          la.assembleElementMatrixPattern(e, localPattern[0][0]);
          elementContributionsComputed = true;
        }
 
        // Only assemble on inner or boundary intersections if the assembler
        // supports this.
        if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> ||
                      isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler>)
        {
          // We only need to traverse the intersections if the assembler
          // supports interior intersections or if there are boundary intersections.
          if (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler> ||
              hasBoundaryIntersections)
          {
            for (auto const& is : intersections(rowBasis_.gridView(), e))
            {
              // the intersection is on a domain boundary
              if (is.boundary() && !is.neighbor())
              {
                if constexpr (isLocalBoundaryIntersectionMatrixAssembler<LocalAssembler>)
                {
                  la.assembleBoundaryIntersectionMatrixPattern(is, localPattern[0][0]);
                  elementContributionsComputed = true;
                }
              }
 
              // the intersection is a regular inner boundary
              if (!is.boundary() && is.neighbor())
              {
                if constexpr (isLocalInteriorIntersectionMatrixAssembler<LocalAssembler>)
                {
                  const auto outsideIndex = mapper.index(is.outside());
                  // visit each intersection only once
                  if (insideIndex < outsideIndex)
                  {
                    outsideLocalViews.bind(is.outside());
 
                    localPattern[0][1].resize(get<0>(localViews).size(), get<1>(outsideLocalViews).size());
                    localPattern[1][0].resize(get<0>(outsideLocalViews).size(), get<1>(localViews).size());
                    localPattern[1][1].resize(get<0>(outsideLocalViews).size(), get<1>(outsideLocalViews).size());
 
                    localPattern[0][1].clear();
                    localPattern[1][0].clear();
                    localPattern[1][1].clear();
 
                    la.assembleInteriorIntersectionMatrixPattern(is, localPattern);
                    elementContributionsComputed = true;
 
                    matrixPattern.scatter(get<0>(localViews),get<1>(outsideLocalViews), localPattern[0][1]);
                    matrixPattern.scatter(get<0>(outsideLocalViews),get<1>(localViews), localPattern[1][0]);
                    matrixPattern.scatter(get<0>(outsideLocalViews),get<1>(outsideLocalViews), localPattern[1][1]);
                  }
                }
              }
            }
          }
        }
        if (elementContributionsComputed)
          matrixPattern.scatter(get<0>(localViews),get<1>(localViews), localPattern[0][0]);
      });
    });
  }
 
  template <Concepts::LocalVectorAssembler<typename RowBasis::LocalView> LocalAssembler, class Vector>
  void assembleVectorEntries (LocalAssembler& la, Vector& vector)
  {
    executor_([&, la](auto context) mutable
    {
      auto localView = rowBasis_.localView();
      la.bindLocalView(localView);
 
      auto outsideLocalView = rowBasis_.localView();
      if constexpr (isLocalInteriorIntersectionVectorAssembler<LocalAssembler>)
        la.bindOutsideLocalView(outsideLocalView);
 
      // local vectors for local assembling
      using LocalVector = typename Vector::LocalVector;
      std::array<LocalVector,2> localVectors{LocalVector(localView.maxSize())};
 
      auto mapper = MultipleCodimMultipleGeomTypeMapper{rowBasis_.gridView(), mcmgElementLayout()};
      context.forEach(elements(rowBasis_.gridView()), [&](const auto& e)
      {
        const auto insideIndex = mapper.index(e);
        localView.bind(e);
 
        // initialize the matrix to zero
        localVectors[0].resize(localView.size());
        localVectors[0] = 0;
 
        // Whenever any contribution was computed on this
        // element, we need to scatter the local element data.
        bool elementContributionsComputed = false;
 
        const bool hasBoundaryIntersections = e.hasBoundaryIntersections();
 
        // only bind to an element if necessary
        if (isLocalElementVectorAssembler<LocalAssembler> ||
            isLocalInteriorIntersectionVectorAssembler<LocalAssembler> || (isLocalBoundaryIntersectionVectorAssembler<LocalAssembler> && hasBoundaryIntersections))
          la.bindElement(e);
 
        if constexpr (isLocalElementVectorAssembler<LocalAssembler>)
        {
          la.assembleElementVector(e, localVectors[0]);
          elementContributionsComputed = true;
        }
 
        // Only assemble on inner or boundary intersections if the assembler
        // supports this.
        if constexpr (isLocalInteriorIntersectionVectorAssembler<LocalAssembler> ||
                      isLocalBoundaryIntersectionVectorAssembler<LocalAssembler>)
        {
          // We only need to traverse the intersections if the assembler
          // supports interior intersections or if there are boundary intersections.
          if (isLocalInteriorIntersectionVectorAssembler<LocalAssembler> ||
              hasBoundaryIntersections)
          {
            for (auto const& is : intersections(rowBasis_.gridView(), e))
            {
              // the intersection is on a domain boundary
              if (is.boundary() && !is.neighbor())
              {
                if constexpr (isLocalBoundaryIntersectionVectorAssembler<LocalAssembler>)
                {
                  la.assembleBoundaryIntersectionVector(is, localVectors[0]);
                  elementContributionsComputed = true;
                }
              }
 
              // the intersection is a regular inner boundary
              if (!is.boundary() && is.neighbor())
              {
                if constexpr (isLocalInteriorIntersectionVectorAssembler<LocalAssembler>)
                {
                  const auto outsideIndex = mapper.index(is.outside());
                  // visit each intersection only once
                  if (insideIndex < outsideIndex)
                  {
                    outsideLocalView.bind(is.outside());
 
                    // initialize the coupling matrices to zero
                    localVectors[1].resize(outsideLocalView.size());
                    localVectors[1] = 0;
 
                    la.assembleInteriorIntersectionVector(is, localVectors);
                    elementContributionsComputed = true;
 
                    vector.scatter(outsideLocalView, localVectors[1]);
                  }
                }
              }
            }
          }
        }
        if (elementContributionsComputed)
          vector.scatter(localView, localVectors[0]);
      });
    });
  }
 
};
 
 
 
} // end namespace Dune::Assembler
 
#endif // DUNE_ASSEMBLER_GLOBALASSEMBLER_HH