functionhelper.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@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 worldDim, typename Basis>
auto globalIndexFromGlobalPosition(const Basis& basis, const Dune::FieldVector<double, worldDim>& 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;
  using Element        = typename LocalView::Element;
  constexpr int myDim  = Element::mydimension;
  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, worldDim>> globalIndices;
 
  auto fT = [&](int nodeNumber, Dune::FieldVector<double, myDim>&& localCoordinate) {
    if (Dune::FloatCmp::eq(geo.global(localCoordinate), pos, tol)) {
      globalIndices.emplace();
      for (int j = 0; j < worldDim; 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);
}
 
template <typename T>
decltype(auto) maybeDeref(T& t) {
  if constexpr (Concepts::PointerOrSmartPointer<T>)
    return *t;
  else
    return t;
}
 
template <typename A>
requires(std::is_pointer_v<std::remove_cvref_t<A>> || traits::isSharedPtr<std::remove_cvref_t<A>>::value)
auto obtainForcesDueToIDBC(const A& assembler) {
  const auto& K = assembler->matrix(DBCOption::Raw);
  auto& dv      = assembler->dirichletValues();
  auto newInc   = Eigen::VectorXd::Zero(dv.size()).eval();
  dv.evaluateInhomogeneousBoundaryConditionDerivative(newInc, 1.0);
 
  Eigen::VectorXd F_dirichlet;
  F_dirichlet = K * newInc;
 
  if (assembler->dBCOption() == DBCOption::Full)
    assembler->dirichletValues().setZeroAtConstrainedDofs(F_dirichlet);
  else
    F_dirichlet = assembler->createReducedVector(F_dirichlet);
  return F_dirichlet;
}
 
} // namespace Ikarus::utils