fehelper.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 <functional>
#include <optional>
 
#include <dune/localfefunctions/manifolds/realTuple.hh>
 
#include <ikarus/utils/traversal.hh>
 
namespace Ikarus::FEHelper {
template <typename Traits, typename Vector, typename ST>
auto localSolutionBlockVector(
    const Vector& x, const typename Traits::LocalView& localView,
    const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& dx = std::nullopt) {
  constexpr int worldDim = Traits::worlddim;
  const auto& fe         = localView.tree().child(0).finiteElement();
  Dune::BlockVector<Dune::RealTuple<ST, worldDim>> localX(fe.size());
  if (dx) {
    for (auto i = 0U; i < localX.size(); ++i)
      for (auto j = 0U; j < worldDim; ++j)
        localX[i][j] =
            dx.value().get()[i * worldDim + j] + x[localView.index(localView.tree().child(j).localIndex(i))[0]];
  } else {
    for (auto i = 0U; i < localX.size(); ++i)
      for (auto j = 0U; j < worldDim; ++j)
        localX[i][j] = x[localView.index(localView.tree().child(j).localIndex(i))[0]];
  }
 
  return localX;
}
 
template <typename Traits, typename Vector, typename ST>
auto localSolutionBlockVectorScalar(
    const Vector& x, const typename Traits::LocalView& localView, auto&& treePath,
    const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& dx = std::nullopt) {
  constexpr int worldDim = Traits::worlddim;
  const auto& fe         = localView.tree().child(treePath).finiteElement();
  Dune::BlockVector<Dune::RealTuple<ST, 1>> localX(fe.size());
  if (dx) {
    for (auto i = 0U; i < localX.size(); ++i) {
      auto iIdx    = localView.index(localView.tree().child(treePath).localIndex(i))[0];
      localX[i][0] = dx.value().get()[iIdx] + x[iIdx];
    }
  } else {
    for (auto i = 0U; i < localX.size(); ++i)
      localX[i][0] = x[localView.index(localView.tree().child(treePath).localIndex(i))[0]];
  }
 
  return localX;
}
 
template <typename Traits, typename Vector, typename ST>
auto localSolutionBlockVectorComposite(
    const Vector& x, const typename Traits::LocalView& localView, auto&& treePath,
    const std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>& dx = std::nullopt) {
  constexpr int worldDim = Traits::worlddim;
  const auto& fe         = localView.tree().child(treePath).child(0).finiteElement();
  Dune::BlockVector<Dune::RealTuple<ST, worldDim>> localX(fe.size());
  if (dx) {
    for (auto i = 0U; i < localX.size(); ++i)
      for (auto j = 0U; j < worldDim; ++j) {
        auto ijIndx  = localView.index(localView.tree().child(treePath).child(j).localIndex(i))[0];
        localX[i][j] = dx.value().get()[ijIndx] + x[ijIndx];
      }
  } else {
    for (auto i = 0U; i < localX.size(); ++i)
      for (auto j = 0U; j < worldDim; ++j)
        localX[i][j] = x[localView.index(localView.tree().child(treePath).child(j).localIndex(i))[0]];
  }
 
  return localX;
}
 
namespace Impl {
  template <typename LocalView, typename Node>
  void leafNodeIndices(const LocalView& localView, const Node& node,
                       std::vector<typename LocalView::MultiIndex>& globalIndices) {
    const auto& fe = node.finiteElement();
    for (size_t i = 0; i < fe.size(); ++i)
      globalIndices.push_back(localView.index(node.localIndex(i)));
  }
 
  template <typename LocalView, typename Node>
  void powerNodeIndices(const LocalView& localView, const Node& node,
                        std::vector<typename LocalView::MultiIndex>& globalIndices) {
    const auto& fe         = node.child(0).finiteElement();
    const int childrenSize = node.degree();
    for (size_t i = 0; i < fe.size(); ++i)
      for (int j = 0; j < childrenSize; ++j)
        globalIndices.push_back(localView.index(node.child(j).localIndex(i)));
  }
} // namespace Impl
 
template <typename LocalView>
void globalIndicesFromLocalView(const LocalView& localView,
                                std::vector<typename LocalView::MultiIndex>& globalIndices) {
  assert(localView.bound());
  globalIndices.clear();
  using namespace Dune::Indices;
  using namespace FEHelper::Impl;
 
  auto leafOpFunc = [&](auto&& node, [[maybe_unused]] auto&& treePath) {
    leafNodeIndices(localView, node, globalIndices);
  };
 
  auto powerOpFunc = [&](auto&& node, [[maybe_unused]] auto&& treePath) {
    powerNodeIndices(localView, node, globalIndices);
  };
 
  utils::forEachLeafOrPowerLeafNode(localView.tree(), Dune::TypeTree::hybridTreePath(), powerOpFunc, leafOpFunc);
}
 
template <typename FiniteElement>
void globalIndices(const FiniteElement& fe, std::vector<typename FiniteElement::LocalView::MultiIndex>& globalIndices) {
  globalIndicesFromLocalView(fe.localView(), globalIndices);
}
} // namespace Ikarus::FEHelper