functionhelper.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <ranges>
 
#include <dune/common/float_cmp.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
 
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/traversal.hh>
 
namespace Ikarus::utils {
template <int size, typename LV>
void obtainLagrangeGlobalNodePositions(const LV& localView,
                                       std::vector<Dune::FieldVector<double, size>>& lagrangeNodeGlobalCoords) {
  auto fT = [&]([[maybe_unused]] int nodeNumber, Dune::FieldVector<double, size>&& localCoordinate) {
    lagrangeNodeGlobalCoords.emplace_back(localView.element().geometry().global(localCoordinate));
    return false;
  };
  forEachLagrangeNodePosition(localView, fT);
}
 
template <int size, typename Basis>
auto globalIndexFromGlobalPosition(const Basis& basis, const Dune::FieldVector<double, size>& pos) {
  static_assert(Concepts::LagrangeNode<std::remove_cvref_t<decltype(basis.localView().tree().child(0))>>,
                "globalIndexFromGlobalPosition is only supported for Lagrange basis");
  constexpr double tol = 1e-8;
  using LocalView      = std::remove_cvref_t<decltype(basis.localView())>;
  using MultiIndex     = typename LocalView::MultiIndex;
  Dune::HierarchicSearch hSearch(basis.gridView().grid(), basis.gridView().indexSet());
  const auto& ele = hSearch.findEntity(pos);
  auto localView  = basis.localView();
  localView.bind(ele);
  const auto geo   = localView.element().geometry();
  const auto& node = localView.tree();
  std::optional<std::array<MultiIndex, size>> globalIndices;
 
  auto fT = [&](int nodeNumber, Dune::FieldVector<double, size>&& localCoordinate) {
    if (Dune::FloatCmp::eq(geo.global(localCoordinate), pos, tol)) {
      globalIndices.emplace();
      for (int j = 0; j < size; j++)
        globalIndices.value()[j] = localView.index(node.child(j).localIndex(nodeNumber));
      return true;
    }
    return false;
  };
  forEachLagrangeNodePosition(localView, fT);
  if (not globalIndices.has_value())
    DUNE_THROW(Dune::InvalidStateException, "No Lagrange node found at the given position in the grid.");
  return globalIndices.value();
}
 
template <typename FE>
auto referenceElementSubEntityPositions(FE& fe, int codim) {
  constexpr int dim            = FE::Traits::mydim;
  const auto& element          = fe.gridElement();
  const auto& referenceElement = Dune::referenceElement<double, dim>(element.type());
  const int numberOfVertices   = referenceElement.size(codim);
 
  auto getPosition = [=](const int i) { return referenceElement.position(i, codim); };
  return std::views::transform(std::views::iota(0, numberOfVertices), getPosition);
}
 
template <typename FE>
auto referenceElementVertexPositions(FE& fe) {
  return referenceElementSubEntityPositions(fe, FE::Traits::mydim);
}
 
} // namespace Ikarus::utils