traversal.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2026 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include <ikarus/utils/concepts.hh>
 
namespace Ikarus::utils {
 
template <class T, class TreePath, class PowerFunc, class LeafFunc>
void forEachLeafOrPowerLeafNode(T&& tree, TreePath&& treePath, PowerFunc&& powerFunc, LeafFunc&& leafFunc) {
  using Tree = std::decay_t<T>;
  if constexpr (Tree::isLeaf) {
    leafFunc(tree, treePath);
  } else if constexpr (Tree::isPower) {
    if constexpr (Tree::template Child<Dune::Indices::_0>::Type::isLeaf) {
      powerFunc(tree, treePath);
    } else {
      for (std::size_t i = 0; i < tree.degree(); ++i) {
        auto childTreePath = Dune::TypeTree::push_back(treePath, i);
        forEachLeafOrPowerLeafNode(tree.child(i), childTreePath, powerFunc, leafFunc);
      }
    }
  } else {
    auto indices = std::make_index_sequence<Tree::degree()>{};
    Dune::Hybrid::forEach(indices, [&](auto i) {
      auto childTreePath = Dune::TypeTree::push_back(treePath, i);
      forEachLeafOrPowerLeafNode(tree.child(i), childTreePath, powerFunc, leafFunc);
    });
  }
}
 
namespace Impl {
  template <bool isScalar, typename NT>
  void checkLagrangeNode(const NT& node) {
    if constexpr (isScalar) {
      static_assert(Concepts::LagrangeNode<std::remove_cvref_t<decltype(node)>>, "Only supported for Lagrange basis");
    } else {
      static_assert(Concepts::LagrangeNode<std::remove_cvref_t<decltype(node.child(0))>>,
                    "Only supported for Lagrange basis");
    }
  }
} // namespace Impl
 
template <typename LV, typename F, int myDim = LV::Element::mydimension>
requires(std::convertible_to<F, std::function<bool(int, Dune::FieldVector<double, myDim> &&)>>)
void forEachLagrangeNodePosition(const LV& localView, F&& f) {
  using Tree                     = LV::Tree;
  static constexpr bool isScalar = Tree::isLeaf;
  static constexpr bool isPower  = Tree::isPower;
  const auto& node               = [&]() {
    if constexpr (isScalar or isPower)
      return localView.tree();
    else
      return localView.tree().template child<0>();
  }();
 
  Impl::checkLagrangeNode<isScalar>(node);
  assert(localView.bound() && "The local view must be bound to an element");
  const auto& localFE = [&]() {
    if constexpr (isScalar)
      return node.finiteElement();
    else
      return node.child(0).finiteElement();
  }();
  std::vector<Dune::FieldVector<double, myDim>> lagrangeNodeCoords;
  lagrangeNodeCoords.resize(localFE.size());
  std::vector<double> out;
  for (int i = 0; i < myDim; i++) {
    auto ithCoord = [&i](const Dune::FieldVector<double, myDim>& x) { return x[i]; };
    localFE.localInterpolation().interpolate(ithCoord, out);
    for (std::size_t j = 0; j < out.size(); j++)
      lagrangeNodeCoords[j][i] = out[j];
  }
  for (int nodeNumber = 0; auto& nCoord : lagrangeNodeCoords)
    if (f(nodeNumber++, std::move(nCoord)))
      break;
}
} // namespace Ikarus::utils