utils/dirichletvalues.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 <algorithm>
#include <concepts>
#include <cstddef>
#include <functional>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/common/float_cmp.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/indices.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
 
#include <Eigen/Core>
 
#include <autodiff/forward/real/real.hpp>
 
#include <ikarus/utils/concepts.hh>
 
namespace Dune {
template <class K, int N>
class FieldVector;
}
 
namespace Ikarus {
 
template <class>
struct DeriveSizeType;
 
template <Concepts::EigenVector T>
struct DeriveSizeType<T>
{
  using SizeType = Eigen::Index;
};
 
template <>
struct DeriveSizeType<std::vector<bool>>
{
  using SizeType = std::vector<bool>::size_type;
};
 
namespace Impl {
  template <typename Tree>
  struct PreBasisInfo
  {
  };
 
  template <typename Tree>
  requires(Tree::isLeaf)
  struct PreBasisInfo<Tree>
  {
    static constexpr std::size_t size = 0;
    using NodalSolutionType           = double;
  };
 
  template <typename Tree>
  requires(Tree::isPower)
  struct PreBasisInfo<Tree>
  {
    static constexpr std::size_t size = Tree::degree();
    using NodalSolutionType           = Eigen::Vector<double, size>;
  };
 
  template <typename Tree>
  requires(Tree::isComposite)
  struct PreBasisInfo<Tree>
  {
    using ChildTreeType = Tree::template Child<0>::Type;
    static_assert(not ChildTreeType::isComposite,
                  "DirichletValues is not implemented to handle a composite basis within a composite basis.");
 
    static constexpr std::size_t size = PreBasisInfo<ChildTreeType>::size;
    using NodalSolutionType           = PreBasisInfo<ChildTreeType>::NodalSolutionType;
  };
} // namespace Impl
 
template <typename B, typename FC = std::vector<bool>>
class DirichletValues
{
public:
  using Basis                         = std::remove_cvref_t<B>;
  using FlagsType                     = FC;
  using BackendType                   = decltype(Dune::Functions::istlVectorBackend(std::declval<FlagsType&>()));
  using SizeType                      = typename DeriveSizeType<FlagsType>::SizeType;
  static constexpr int worldDimension = Basis::GridView::dimensionworld;
  using LocalView                     = Basis::LocalView;
  using Tree                          = LocalView::Tree;
  using NodalSolutionType             = typename Impl::PreBasisInfo<Tree>::NodalSolutionType;
  static constexpr std::size_t numberOfChildrenAtNode = Impl::PreBasisInfo<Tree>::size;
  explicit DirichletValues(const B& basis)
      : basis_{basis},
        dirichletFlagsBackend_{dirichletFlags_} {
    dirichletFlagsBackend_.resize(basis_);
    std::fill(dirichletFlags_.begin(), dirichletFlags_.end(), false);
    if constexpr (Tree::isComposite) {
      Dune::Hybrid::forEach(
          Dune::Hybrid::integralRange(Dune::index_constant<std::tuple_size_v<typename Basis::PreBasis::SubPreBases>>()),
          [&](const auto i) {
            static_assert(not Tree::template Child<i>::Type::isComposite,
                          "DirichletValues is not implemented to handle a composite basis within a composite basis.");
          });
    }
  }
 
  template <typename F, typename TreePath = Dune::TypeTree::HybridTreePath<>>
  void fixBoundaryDOFs(F&& f, TreePath&& treePath = {}) {
    using namespace Dune::Functions;
    using SubSpaceLocalView =
        typename std::remove_cvref_t<decltype(subspaceBasis(basis_, std::forward<TreePath>(treePath)))>::LocalView;
 
    if constexpr (Concepts::IsFunctorWithArgs<F, BackendType, typename Basis::MultiIndex>) {
      auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); };
      Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
    } else if constexpr (Concepts::IsFunctorWithArgs<F, BackendType, int, SubSpaceLocalView>) {
      auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); };
      Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
    } else if constexpr (Concepts::IsFunctorWithArgs<F, BackendType, int, SubSpaceLocalView,
                                                     typename Basis::GridView::Intersection>) {
      auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) {
        f(dirichletFlagsBackend_, localIndex, localView, intersection);
      };
      Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
    } else {
      static_assert(Dune::AlwaysFalse<F>(), "fixBoundaryDOFs: A function with this signature is not supported");
    }
  }
 
  template <typename F>
  void fixDOFs(F&& f) {
    f(basis_, dirichletFlagsBackend_);
  }
 
  template <typename MultiIndex>
  requires(not std::integral<MultiIndex>)
  void setSingleDOF(const MultiIndex i, bool flag) {
    dirichletFlagsBackend_[i] = flag;
  }
 
  void setSingleDOF(std::size_t i, bool flag)
  requires(std::same_as<typename Basis::MultiIndex, Dune::Functions::FlatMultiIndex<size_t>>)
  {
    dirichletFlags_[i] = flag;
  }
 
  void reset() { std::fill(dirichletFlags_.begin(), dirichletFlags_.end(), false); }
 
  /* \brief Returns the local basis object */
  const auto& basis() const { return basis_; }
 
  /* \brief Returns a boolean values, if the give multiIndex is constrained */
  template <typename MultiIndex>
  requires(not std::integral<MultiIndex>)
  [[nodiscard]] bool isConstrained(const MultiIndex& multiIndex) const {
    return dirichletFlagsBackend_[multiIndex];
  }
 
  /* \brief Returns a boolean values, if the i-th degree of freedom is constrained */
  [[nodiscard]] bool isConstrained(std::size_t i) const
  requires(std::same_as<typename Basis::MultiIndex, Dune::Functions::FlatMultiIndex<size_t>>)
  {
    return dirichletFlags_[i];
  }
 
  /* \brief Returns how many degrees of freedoms are fixed */
  auto fixedDOFsize() const { return std::ranges::count(dirichletFlags_, true); }
 
  /* \brief Returns how many degrees of freedom there are */
  auto size() const { return dirichletFlags_.size(); }
 
  /* \brief Returns the underlying dof flag container */
  auto& container() const { return dirichletFlags_; }
 
  template <typename F>
  void storeInhomogeneousBoundaryCondition(F&& f, double lambda = 1.0) {
    auto derivativeLambda = [&](const auto& globalCoord, const double& lambda) {
      autodiff::real lambdaDual = lambda;
      lambdaDual[1]             = 1; // Setting the derivative in lambda direction to 1
      return derivative(f(globalCoord, lambdaDual));
    };
    dirichletFunctions_.push_back({f, derivativeLambda});
    setInhomogeneousBoundaryConditionFlag(lambda);
  }
 
  void setZeroAtConstrainedDofs(Eigen::VectorXd& xIh) const {
    for (Eigen::Index i = 0; i < xIh.size(); ++i)
      if (this->isConstrained(i))
        xIh[i] = 0.0;
  }
 
  void evaluateInhomogeneousBoundaryCondition(Eigen::VectorXd& xIh, const double& lambda) const {
    Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
    inhomogeneousBoundaryVectorDummy.setZero(this->size());
    xIh.resizeLike(inhomogeneousBoundaryVectorDummy);
    xIh.setZero();
    auto runInterpolateImpl = [&](const auto& basis) {
      for (auto& f : dirichletFunctions_) {
        interpolate(basis, inhomogeneousBoundaryVectorDummy,
                    [&](const auto& globalCoord) { return f.value(globalCoord, lambda); });
        xIh += inhomogeneousBoundaryVectorDummy;
      }
    };
    runInterpolate(runInterpolateImpl);
  }
 
  void evaluateInhomogeneousBoundaryConditionDerivative(Eigen::VectorXd& xIh, const double& lambda) const {
    Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
    inhomogeneousBoundaryVectorDummy.setZero(this->size());
    xIh.resizeLike(inhomogeneousBoundaryVectorDummy);
    xIh.setZero();
    auto runInterpolateImpl = [&](const auto& basis) {
      for (auto& f : dirichletFunctions_) {
        interpolate(basis, inhomogeneousBoundaryVectorDummy,
                    [&](const auto& globalCoord) { return f.derivative(globalCoord, lambda); });
        xIh += inhomogeneousBoundaryVectorDummy;
      }
    };
    runInterpolate(runInterpolateImpl);
  }
 
private:
  Basis basis_;
  FlagsType dirichletFlags_;
  BackendType dirichletFlagsBackend_;
  struct DirichletFunctions
  {
    using Signature = std::function<NodalSolutionType(const Dune::FieldVector<double, worldDimension>&, const double&)>;
    Signature value;
    Signature derivative;
  };
  std::vector<DirichletFunctions> dirichletFunctions_;
 
  void setInhomogeneousBoundaryConditionFlag(double lambda) {
    Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
    inhomogeneousBoundaryVectorDummy.setZero(this->size());
    this->evaluateInhomogeneousBoundaryCondition(inhomogeneousBoundaryVectorDummy, lambda);
    for (const std::size_t i : Dune::range(this->size()))
      if (Dune::FloatCmp::ne(inhomogeneousBoundaryVectorDummy[i], 0.0))
        this->setSingleDOF(i, true);
  }
 
  template <typename F>
  void runInterpolate(F&& f) const {
    if constexpr (Tree::isComposite)
      f(subspaceBasis(basis_, Dune::Indices::_0));
    else
      f(basis_);
  }
};
 
} // namespace Ikarus