traversal.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_TRAVERSAL_HH
#define DUNE_FUFEM_BACKENDS_TRAVERSAL_HH
 
#include <type_traits>
#include <utility>
#include <functional>
#include <variant>
 
#include <dune/common/typetraits.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/hybridutilities.hh>
 
#include <dune/istl/multitypeblockmatrix.hh>
 
#include <dune/assembler/backends/singlerowcolmatrix.hh>
 
#include <dune/functions/common/multiindex.hh>
 
 
 
namespace Dune::Fufem::Impl {
 
  // Shallow traversal implementation
 
  template<class Vector, class F>
  static void forEachVectorEntry(Vector&& vector, F&& f)
  {
    Dune::Hybrid::forEach(Dune::range(Dune::Hybrid::size(vector)), [&](auto i) {
      f(vector[i], i);
    });
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(Matrix&& matrix, F&& f)
  {
    for(auto i : Dune::range(matrix.N()))
      for(auto&& [matrix_ij, j] : Dune::sparseRange(matrix[i]))
        f(matrix_ij, i, j);
  }
 
  template<class... MatrixRows, class F>
  static void forEachMatrixEntry(const Dune::MultiTypeBlockMatrix<MatrixRows...>& matrix, F&& f)
  {
    using Matrix = Dune::MultiTypeBlockMatrix<MatrixRows...>;
    Dune::Hybrid::forEach(std::make_index_sequence<Matrix::N()>(), [&](auto i) {
      Dune::Hybrid::forEach(std::make_index_sequence<Matrix::M()>(), [&](auto j) {
        f(matrix[i][j], i, j);
      });
    });
  }
 
  template<class... MatrixRows, class F>
  static void forEachMatrixEntry(Dune::MultiTypeBlockMatrix<MatrixRows...>& matrix, F&& f)
  {
    using Matrix = Dune::MultiTypeBlockMatrix<MatrixRows...>;
    Dune::Hybrid::forEach(std::make_index_sequence<Matrix::N()>(), [&](auto i) {
      Dune::Hybrid::forEach(std::make_index_sequence<Matrix::M()>(), [&](auto j) {
        f(matrix[i][j], i, j);
      });
    });
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(const Dune::Assembler::SingleRowMatrix<Matrix>& matrix, F&& f)
  {
    for(auto&& [matrix_0j, j] : Dune::sparseRange(matrix[0]))
      f(matrix_0j, std::monostate(), j);
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(Dune::Assembler::SingleRowMatrix<Matrix>& matrix, F&& f)
  {
    for(auto&& [matrix_0j, j] : Dune::sparseRange(matrix[0]))
      f(matrix_0j, std::monostate(), j);
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(const Dune::Assembler::SingleColumnMatrix<Matrix>& matrix, F&& f)
  {
    for(auto i : Dune::range(matrix.N()))
      for(auto&& [matrix_i0, j] : Dune::sparseRange(matrix[i]))
        f(matrix_i0, i, std::monostate());
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(Dune::Assembler::SingleColumnMatrix<Matrix>& matrix, F&& f)
  {
    for(auto i : Dune::range(matrix.N()))
      for(auto&& [matrix_i0, j] : Dune::sparseRange(matrix[i]))
        f(matrix_i0, i, std::monostate());
  }
 
  // Utilities for deep traversal
 
  template<class size_type, std::size_t n, class Tail>
  auto extendMultiIndex(const Dune::Functions::StaticMultiIndex<size_type, n>& head, Tail tail)
  {
    if constexpr (std::is_same_v<Tail, std::monostate>)
      return head;
    else
      return std::apply([&](const auto&... headEntries) {
        return Dune::Functions::StaticMultiIndex<size_type, n+1>{headEntries..., tail};
      }, head);
  }
 
  template<class F, class MultiIndex>
  class RecursiveVectorVisitor
  {
  public:
 
    RecursiveVectorVisitor(F f, MultiIndex index)
      : f_(f)
      , index_(index)
    {}
 
    template<class Entry, class FlatIndex>
    void operator()(Entry&& vector_i, FlatIndex i)
    {
      if constexpr (Dune::IsNumber<typename Dune::AutonomousValueType<std::decay_t<Entry>>::type>::value)
        f_(vector_i, extendMultiIndex(index_, i));
      else
        Impl::forEachVectorEntry(vector_i, Dune::Fufem::Impl::RecursiveVectorVisitor(f_, extendMultiIndex(index_, i)));
    }
 
  private:
    F f_;
    const MultiIndex index_;
  };
 
 
  template<class F, class RowMultiIndex, class ColMultiIndex>
  class RecursiveMatrixVisitor
  {
  public:
 
    RecursiveMatrixVisitor(F f, RowMultiIndex rowIndex, ColMultiIndex colIndex)
      : f_(f)
      , rowIndex_(rowIndex)
      , colIndex_(colIndex)
    {}
 
    template<class Entry, class FlatRowIndex, class FlatColIndex>
    void operator()(Entry&& matrix_ij, FlatRowIndex i, FlatColIndex j)
    {
      if constexpr (Dune::IsNumber<typename Dune::AutonomousValueType<std::decay_t<Entry>>::type>::value)
        f_(matrix_ij, extendMultiIndex(rowIndex_, i), extendMultiIndex(colIndex_, j));
      else
        Impl::forEachMatrixEntry(matrix_ij, Dune::Fufem::Impl::RecursiveMatrixVisitor(f_, extendMultiIndex(rowIndex_, i), extendMultiIndex(colIndex_, j)));
    }
 
  private:
    F f_;
    const RowMultiIndex rowIndex_;
    const RowMultiIndex colIndex_;
  };
 
 
} // namespace Dune::Fufem::Impl
 
 
 
namespace Dune::Fufem {
 
  template<class Vector, class F>
  static void forEachVectorEntry(Vector&& vector, F&& f)
  {
    Impl::forEachVectorEntry(vector, f);
  }
 
  template<class Matrix, class F>
  static void forEachMatrixEntry(Matrix&& matrix, F&& f)
  {
    Impl::forEachMatrixEntry(matrix, f);
  }
 
  template<class Vector, class F>
  static void recursiveForEachVectorEntry(Vector&& vector, F&& f)
  {
    using EmptyIndex = Dune::Functions::StaticMultiIndex<std::size_t, 0>;
    auto recursive_f = Dune::Fufem::Impl::RecursiveVectorVisitor(std::ref(f), EmptyIndex{});
    Impl::forEachVectorEntry(vector, recursive_f);
  }
 
  template<class Matrix, class F>
  static void recursiveForEachMatrixEntry(Matrix&& matrix, F&& f)
  {
    using EmptyIndex = Dune::Functions::StaticMultiIndex<std::size_t, 0>;
    auto recursive_f = Dune::Fufem::Impl::RecursiveMatrixVisitor(std::ref(f), EmptyIndex{}, EmptyIndex{});
    Impl::forEachMatrixEntry(matrix, recursive_f);
  }
 
} // namespace Dune::Fufem
 
 
 
#endif // DUNE_FUFEM_BACKENDS_TRAVERSAL_HH