istlmatrixbackend.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_BACKENDS_ISTLMATRIXBACKEND_HH
#define DUNE_FUFEM_BACKENDS_ISTLMATRIXBACKEND_HH
 
#ifndef DUNE_FUFEM_DISABLE_HEADER_DEPRECATION
#warning This header is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.
#endif
 
#include <type_traits>
#include <utility>
 
#include <dune/common/deprecated.hh>
#include <dune/common/indices.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/hybridutilities.hh>
 
#include <dune/functions/common/indexaccess.hh>
 
#include <dune/assembler/backends/singlerowcolmatrix.hh>
 
#include <dune/fufem/backends/matrixbuilder.hh>
 
 
 
namespace Dune::Fufem::Concept {
 
// Concept for a ConstMatrixBackend
// This class is deprecated and will be removed after 2.11.
template<class RowBasis, class ColBasis=RowBasis>
struct
ConstMatrixBackend
{
  using RowMultiIndex = typename RowBasis::MultiIndex;
  using ColMultiIndex = typename ColBasis::MultiIndex;
 
  template<class V>
  auto require(const V& v) -> decltype(
    v(std::declval<RowMultiIndex>(), std::declval<ColMultiIndex>())
  );
};
 
// Concept for a MatrixBackend
template<class RowBasis, class ColBasis=RowBasis>
struct
[[deprecated("This class is deprecated and will be removed after 2.11.")]]
MatrixBackend : Dune::Concept::Refines< ConstMatrixBackend<RowBasis, ColBasis> >
{
  using RowMultiIndex = typename RowBasis::MultiIndex;
  using ColMultiIndex = typename ColBasis::MultiIndex;
 
  template<class V>
  auto require(const V& v) -> decltype(
    const_cast<V&>(v).patternBuilder(),
    const_cast<V&>(v)(std::declval<RowMultiIndex>(), std::declval<ColMultiIndex>()) = v(std::declval<RowMultiIndex>(), std::declval<ColMultiIndex>())
  );
};
 
} // namespace Dune::Fufem::Concept
 
 
 
namespace Dune::Fufem {
 
 
 
namespace Impl {
 
 
 
  // The following is essentially a copy of Dune::Functions::hybridIndexAccess
  // But since matrices do not provide size(), we have to use N() instead.
  template<class C, class I, class F,
    typename std::enable_if< Dune::models<Dune::Functions::Imp::Concept::HasDynamicIndexAccess<I>, C>(), int>::type = 0>
  static auto hybridRowIndexAccess(C&& c, const I& i, F&& f)
    -> decltype(f(c[i]))
  {
    return f(c[i]);
  }
 
  template<class C, class I, class F,
    typename std::enable_if< not Dune::models<Dune::Functions::Imp::Concept::HasDynamicIndexAccess<I>, C>(), int>::type = 0>
  static decltype(auto) hybridRowIndexAccess(C&& c, const I& i, F&& f)
  {
    return Hybrid::switchCases(std::make_index_sequence<std::decay_t<C>::N()>(), i,
        [&](const auto& ii) -> decltype(auto){
          return f(c[ii]);
        }, [&]() -> decltype(auto){
          return f(c[Dune::Indices::_0]);
        });
  }
 
 
 
  // Call f(matrix[i][j]) and return the result.
  // This works with dynamic indices i,j, even for a multi-type matrix.
  // However, it requires that f() has a unique return type.
  template<class Matrix, class RowIndex, class ColIndex, class F>
  decltype(auto) visitMatrixEntry(Matrix&& matrix, const RowIndex& i, const ColIndex& j, F&& f)
  {
    using namespace Dune::Functions;
    return hybridRowIndexAccess(matrix, i, [&](auto&& matrix_i) -> decltype(auto) {
      return hybridIndexAccess(matrix_i, j, f);
    });
  }
 
 
 
  // Call f(matrix[i0][j0]...[in][jm]), by recursively resolving row- and column-multi-indices.
  // Whenenver a SingleRowMatrix or SingleColumnMatrix is encountered, a zero row- or column-index
  // is inserted, respectively. The recursion ends, whenever f(matrixEntry) can be called.
  // This works with dynamic indices i,j, even for a multi-type matrix.
  // However, it requires that f() has a unique return type.
  template<class Matrix, class RowIndex, class ColIndex, class F>
  static decltype(auto) visitMatrixEntryRecursive(Matrix& matrix, const RowIndex& iii, const ColIndex& jjj, F&& f)
  {
    using namespace Dune::Indices;
    using namespace Dune::Functions::Imp;
    static constexpr bool isSingleRow = Dune::Assembler::IsSingleRowMatrix<Matrix>::value;
    static constexpr bool isSingleCol = Dune::Assembler::IsSingleColumnMatrix<Matrix>::value;
    auto splitIndex = [] (auto&& multiIndex) { return std::make_tuple(multiIndex[_0], shiftedDynamicMultiIndex<1>(multiIndex)) ;};
    if constexpr (std::is_invocable_v<F, Matrix&>)
    {
      return f(matrix);
    }
    else if constexpr (isSingleRow)
    {
      assert(jjj.size()>0);
      auto [j, jj] = splitIndex(jjj);
      // We have to capture jj explicitly by value, because capturing structured bindings
      // by reference is not allowed before C++20
      return visitMatrixEntry(matrix, _0, j, [&, jj=jj](auto&& matrix_0j) -> decltype(auto) {
        return visitMatrixEntryRecursive(matrix_0j, iii, jj, f);
      });
    }
    else if constexpr (isSingleCol)
    {
      assert(iii.size()>0);
      auto [i, ii] = splitIndex(iii);
      // We have to capture ii explicitly by value, because capturing structured bindings
      // by reference is not allowed before C++20
      return visitMatrixEntry(matrix, i, _0, [&, ii=ii](auto&& matrix_i0) -> decltype(auto) {
        return visitMatrixEntryRecursive(matrix_i0, ii, jjj, f);
      });
    }
    else
    {
      assert(iii.size()>0);
      assert(jjj.size()>0);
      auto [i, ii] = splitIndex(iii);
      auto [j, jj] = splitIndex(jjj);
      // We have to capture ii and jj explicitly by value, because capturing structured bindings
      // by reference is not allowed before C++20
      return visitMatrixEntry(matrix, i, j, [&, ii=ii, jj=jj](auto&& matrix_ij) -> decltype(auto) {
        return visitMatrixEntryRecursive(matrix_ij, ii, jj, f);
      });
    }
  }
 
 
 
} // namespace Impl
 
 
 
template<class M, class E=typename M::field_type>
class
[[deprecated("This class is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.")]]
ISTLMatrixBackend
{
 
  template<class Result>
  struct ToScalar;
 
  template<class R>
  struct ToScalar<R&>
  {
    template<class Matrix,
      std::enable_if_t<std::is_convertible_v<Matrix&, R&>, int> = 0>
    R& operator()(Matrix& matrix) {
      return matrix;
    }
 
    R& operator()(Dune::FieldMatrix<R, 1, 1>& matrix) {
      return matrix[0][0];
    }
  };
 
  template<class R>
  struct ToScalar<const R&>
  {
    template<class Matrix,
      std::enable_if_t<std::is_convertible_v<Matrix&, const R&>, int> = 0>
    const R& operator()(Matrix& matrix) {
      return matrix;
    }
 
    const R& operator()(const Dune::FieldMatrix<R, 1, 1>& matrix) {
      return matrix[0][0];
    }
  };
 
 
 
public:
 
  using Matrix = M;
  using Entry = E;
  using value_type = Entry;
 
  ISTLMatrixBackend(Matrix& matrix) :
    matrix_(&matrix)
  {}
 
  MatrixBuilder<Matrix> patternBuilder()
  {
    return {*matrix_};
  }
 
  template<class RowMultiIndex, class ColMultiIndex>
  const Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col) const
  {
    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<const Entry&>());
  }
 
  template<class RowMultiIndex, class ColMultiIndex>
  Entry& operator()(const RowMultiIndex& row, const ColMultiIndex& col)
  {
    return Impl::visitMatrixEntryRecursive(*matrix_, row, col, ToScalar<Entry&>());
  }
 
  const Matrix& matrix() const
  {
    return *matrix_;
  }
 
  Matrix& matrix()
  {
    return *matrix_;
  }
 
  template<class Value>
  void assign(const Value& value)
  {
    matrix() = value;
  }
 
protected:
 
  Matrix* matrix_;
};
 
 
 
template<class Entry, class Matrix>
[[deprecated("This function is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.")]]
auto istlMatrixBackend(Matrix& matrix)
{
DUNE_NO_DEPRECATED_BEGIN
  return ISTLMatrixBackend<Matrix, Entry>(matrix);
DUNE_NO_DEPRECATED_END
}
 
template<class Matrix>
[[deprecated("This function is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.")]]
auto istlMatrixBackend(Matrix& matrix)
{
DUNE_NO_DEPRECATED_BEGIN
  return ISTLMatrixBackend<Matrix, typename Matrix::field_type>(matrix);
DUNE_NO_DEPRECATED_END
}
 
 
 
template<class RowBasis, class ColBasis, class M>
[[deprecated("This function is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.")]]
decltype(auto) toConstMatrixBackend(M& m) {
DUNE_NO_DEPRECATED_BEGIN
  if constexpr (Dune::models<Dune::Fufem::Concept::ConstMatrixBackend<RowBasis,ColBasis>, M>()) {
    return m;
  } else {
    return Dune::Fufem::istlMatrixBackend(m);
  }
DUNE_NO_DEPRECATED_END
}
 
 
 
template<class RowBasis, class ColBasis, class M>
[[deprecated("This function is deprecated and will be removed after 2.11. Use the ISTLMatrixBackend from dune-assembler instead.")]]
decltype(auto) toMatrixBackend(M& m) {
DUNE_NO_DEPRECATED_BEGIN
  if constexpr (Dune::models<Dune::Fufem::Concept::MatrixBackend<RowBasis, ColBasis>, M>()) {
    return m;
  } else {
    return Dune::Fufem::istlMatrixBackend(m);
  }
DUNE_NO_DEPRECATED_END
}
 
 
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_BACKENDS_ISTLMATRIXBACKEND_HH