utils/concepts.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 <concepts>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
 
#include <Eigen/Dense>
#include <Eigen/Sparse>
 
#include <autodiff/forward/dual/dual.hpp>
 
#include "ikarus/assembler/dirichletbcenforcement.hh"
#include "ikarus/finiteelements/mechanics/materials/tags.hh"
#include <ikarus/utils/traits.hh>
 
namespace Eigen {
template <typename Derived>
struct EigenBase;
}
 
namespace Ikarus {
 
struct ControlInformation;
struct NonLinearSolverInformation;
 
template <typename Derived>
auto transpose(const Eigen::EigenBase<Derived>& A);
namespace Concepts {
 
  template <typename T>
  concept EigenType = std::is_base_of_v<Eigen::MatrixBase<std::decay_t<T>>, std::decay_t<T>>;
 
  template <typename Basis>
  concept FlatInterLeavedBasis = requires {
    std::is_same_v<typename Basis::PreBasis::IndexMergingStrategy, Dune::Functions::BasisFactory::FlatInterleaved>;
  };
 
  namespace Impl {
    template <template <typename, int, typename> class U, typename T>
    struct LagrangeNodeHelper : std::false_type
    {
    };
    template <template <typename, int, typename> class U, typename GV, int k, typename R>
    struct LagrangeNodeHelper<U, U<GV, k, R>> : std::true_type
    {
    };
 
    template <template <typename, int, typename> class U, typename T, int k>
    struct LagrangeNodeHelperOfOrder : std::false_type
    {
    };
    template <template <typename, int, typename> class U, typename GV, int k, typename R>
    struct LagrangeNodeHelperOfOrder<U, U<GV, k, R>, k> : std::true_type
    {
    };
  } // namespace Impl
 
  template <typename N>
  concept LagrangeNode = Impl::LagrangeNodeHelper<Dune::Functions::LagrangeNode, N>::value;
 
  template <typename N, int order>
  concept LagrangeNodeOfOrder = Impl::LagrangeNodeHelperOfOrder<Dune::Functions::LagrangeNode, N, order>::value;
 
  template <typename B>
  concept FlatLexicographicBasis = requires {
    std::is_same_v<typename B::PreBasis::IndexMergingStrategy, Dune::Functions::BasisFactory::FlatLexicographic>;
  };
 
  template <typename B>
  concept FlatIndexBasis = FlatLexicographicBasis<B> or FlatInterLeavedBasis<B>;
 
  template <typename Basis>
  concept BlockedInterLeavedBasis = requires {
    std::is_same_v<typename Basis::PreBasis::IndexMergingStrategy, Dune::Functions::BasisFactory::BlockedInterleaved>;
  };
 
  template <typename B>
  concept BlockedLexicographicBasis = requires {
    std::is_same_v<typename B::PreBasis::IndexMergingStrategy, Dune::Functions::BasisFactory::BlockedLexicographic>;
  };
 
  template <typename DLB>
  concept DuneLocalBasis = requires(DLB& duneLocalBasis) {
    typename DLB::Traits::RangeType;
    typename DLB::Traits::JacobianType;
    DLB::Traits::dimDomain;
    typename DLB::Traits::DomainType;
 
    typename DLB::Traits::DomainFieldType;
    typename DLB::Traits::RangeFieldType;
 
    duneLocalBasis.evaluateFunction(std::declval<typename DLB::Traits::DomainType>(),
                                    std::declval<std::vector<typename DLB::Traits::RangeType>&>());
    duneLocalBasis.evaluateJacobian(std::declval<typename DLB::Traits::DomainType>(),
                                    std::declval<std::vector<typename DLB::Traits::JacobianType>&>());
  };
 
  template <typename B>
  concept BlockedIndexBasis = BlockedLexicographicBasis<B> or BlockedInterLeavedBasis<B>;
 
  template <typename PF, typename F, typename SA>
  concept PathFollowingStrategy = requires(PF pft, F nop, SA args, typename F::Domain req) {
    { pft(args) } -> std::same_as<void>;
    { pft.initialPrediction(req, nop, args) } -> std::same_as<void>;
    { pft.intermediatePrediction(req, nop, args) } -> std::same_as<void>;
  };
 
  template <typename ASS, typename NLSI, typename SA, typename DifferentiableFunction>
  concept AdaptiveStepSizingStrategy =
      requires(ASS adaptiveStepSizing, NLSI info, SA args, DifferentiableFunction nop) {
        { adaptiveStepSizing(info, args, nop) } -> std::same_as<void>;
        { adaptiveStepSizing.targetIterations() } -> std::same_as<int>;
        { adaptiveStepSizing.setTargetIterations(std::declval<int>()) } -> std::same_as<void>;
      };
 
  template <typename LS, typename M, typename V>
  concept LinearSolverCheck = requires(LS linearSolver, M A, V vec) {
    linearSolver.analyzePattern(A);
    linearSolver.factorize(A);
    linearSolver.solve(vec, vec);
  };
 
  template <typename NLS>
  concept NonLinearSolverCheckForPathFollowing = requires {
    not(std::is_same_v<typename NLS::DifferentiableFunction::Domain, double> and
        ((traits::isSpecializationTypeAndNonTypes<
             Eigen::Matrix, typename NLS::DifferentiableFunction::Traits::template Range<1>>::value) or
         (traits::isSpecializationTypeNonTypeAndType<
             Eigen::SparseMatrix, typename NLS::DifferentiableFunction::Traits::template Range<1>>::value)));
  };
 
  template <typename L, typename R>
  concept MultiplyAble = requires(L x, R y) { x* y; };
 
  template <typename L, typename R>
  concept AddAble = requires(L x, R y) { x + y; };
 
  template <typename L, typename R>
  concept SubstractAble = requires(L x, R y) { x - y; };
 
  template <typename L, typename R>
  concept MultiplyAssignAble = requires(L x, R y) { x *= y; };
 
  template <typename L, typename R>
  concept DivideAssignAble = requires(L x, R y) { x /= y; };
 
  template <typename L, typename R>
  concept AddAssignAble = requires(L x, R y) { x += y; };
 
  template <typename L, typename R>
  concept SubstractAssignAble = requires(L x, R y) { x -= y; };
 
  template <typename L, typename R>
  concept DivideAble = requires(L x, R y) { x / y; };
 
  template <typename L>
  concept NegateAble = requires(L x) { -x; };
 
  template <typename L>
  concept TransposeAble = requires(L x) { transpose(x); };
 
  template <typename Op, typename... Args>
  concept IsFunctorWithArgs = requires(Op op, Args... args) { op(args...); };
 
  template <typename V>
  concept EigenVector = static_cast<bool>(V::IsVectorAtCompileTime);
 
  template <typename M>
  concept EigenMatrix = traits::isSpecializationTypeAndNonTypes<Eigen::Matrix, M>::value;
 
  template <typename M>
  concept SparseEigenMatrix = traits::isSpecializationTypeNonTypeAndType<Eigen::SparseMatrix, M>::value;
 
  template <typename M>
  concept DenseOrSparseEigenMatrix = SparseEigenMatrix<M> or EigenMatrix<M>;
 
#define MAKE_EIGEN_FIXED_VECTOR_CONCEPT(Size) \
  template <typename V>                       \
  concept EigenVector##Size =                 \
      static_cast<bool>(V::IsVectorAtCompileTime) and static_cast<bool>(V::SizeAtCompileTime == Size);
 
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(1)
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(2)
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(3)
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(4)
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(5)
  MAKE_EIGEN_FIXED_VECTOR_CONCEPT(6)
 
#define MAKE_EIGEN_FIXED_MATRIX_CONCEPT(Size1, Size2)                                                           \
  template <typename M>                                                                                         \
  concept EigenMatrix##Size1##Size2 = static_cast<bool>(std::remove_cvref_t<M>::RowsAtCompileTime == Size1) and \
                                      static_cast<bool>(std::remove_cvref_t<M>::ColsAtCompileTime == Size2);
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 1)
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 2)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(2, 2)
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 3)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(2, 3)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(3, 2)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(3, 3)
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 4)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(2, 4)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(3, 4)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(4, 2)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(4, 3)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(4, 4)
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 5)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(2, 5)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(3, 5)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(4, 5)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(5, 2)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(5, 3)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(5, 4)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(5, 5)
 
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(1, 6)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(2, 6)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(3, 6)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(4, 6)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(5, 6)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(6, 2)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(6, 3)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(6, 4)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(6, 5)
  MAKE_EIGEN_FIXED_MATRIX_CONCEPT(6, 6)
 
#define MAKE_EIGEN_FIXED_MATRIX_OR_VOIGT_CONCEPT(Size1, Size2) \
  template <typename M>                                        \
  concept EigenMatrixOrVoigtNotation##Size1 = EigenMatrix##Size1##Size1<M> or EigenVector##Size2<M>;
 
  MAKE_EIGEN_FIXED_MATRIX_OR_VOIGT_CONCEPT(1, 1)
  MAKE_EIGEN_FIXED_MATRIX_OR_VOIGT_CONCEPT(2, 3)
  MAKE_EIGEN_FIXED_MATRIX_OR_VOIGT_CONCEPT(3, 6)
 
  namespace Impl {
    template <template <typename...> class MaterialToCheck, typename Material>
    consteval bool isMaterial() {
      if constexpr (traits::isSpecialization<MaterialToCheck, Material>::value)
        return true;
 
      if constexpr (Material::isReduced) {
        if constexpr (traits::isSpecialization<MaterialToCheck, typename Material::Underlying>::value) {
          return true;
        } else {
          return false;
        }
      } else
        return false;
    }
  } // namespace Impl
 
  template <template <typename...> class MaterialToCheck, typename Material>
  concept IsMaterial = Impl::isMaterial<MaterialToCheck, Material>();
 
  namespace Impl {
    template <typename T>
    concept ResultType = requires(T t) {
      typename T::type;                             // The nested type 'type'
      typename T::Vectorizer;                       // The nested type 'Vectorizer'
      typename T::Matricizer;                       // The nested type 'Matricizer'
      { toString(t) } -> std::same_as<std::string>; // The toString function
    };
  } // namespace Impl
 
  template <template <typename, int, int> typename RT>
  concept ResultType =
      Impl::ResultType<RT<double, 1, 1>> or Impl::ResultType<RT<double, 1, 2>> or Impl::ResultType<RT<double, 1, 3>> or
      Impl::ResultType<RT<double, 2, 3>> or Impl::ResultType<RT<double, 3, 3>>;
 
  template <typename T>
  concept FlatAssembler = requires(T t, const typename T::FERequirement& req,
                                   typename T::AffordanceCollectionType affordance, DBCOption dbcOption) {
    { t.requirement() } -> std::convertible_to<typename T::FERequirement&>;
    { t.affordanceCollection() } -> std::convertible_to<typename T::AffordanceCollectionType>;
    { t.dBCOption() } -> std::convertible_to<DBCOption>;
 
    { t.bind(req, affordance, dbcOption) } -> std::same_as<void>;
    { t.bind(req) } -> std::same_as<void>;
    { t.bind(affordance) } -> std::same_as<void>;
    { t.bind(dbcOption) } -> std::same_as<void>;
 
    { t.bound() } -> std::convertible_to<bool>;
    { t.boundToRequirement() } -> std::convertible_to<bool>;
    { t.boundToAffordanceCollection() } -> std::convertible_to<bool>;
    { t.boundToDBCOption() } -> std::convertible_to<bool>;
    { t.estimateOfConnectivity() } -> std::convertible_to<size_t>;
 
    { t.createFullVector(std::declval<Eigen::Ref<const Eigen::VectorXd>>()) } -> std::convertible_to<Eigen::VectorXd>;
    { t.constraintsBelow(std::declval<size_t>()) } -> std::convertible_to<size_t>;
    { t.isConstrained(std::declval<size_t>()) } -> std::convertible_to<bool>;
    { t.size() } -> std::convertible_to<size_t>;
    { t.reducedSize() } -> std::convertible_to<size_t>;
  };
 
  template <typename T>
  concept ScalarFlatAssembler =
      Concepts::FlatAssembler<T> and requires(T t, const typename T::FERequirement& req,
                                              typename T::AffordanceCollectionType affordance, DBCOption dbcOption) {
        { t.scalar(req, affordance.scalarAffordance()) } -> std::convertible_to<const double&>;
        { t.scalar() } -> std::convertible_to<const double&>;
      };
 
  template <typename T>
  concept VectorFlatAssembler = Concepts::ScalarFlatAssembler<T> and
                                requires(T t, const typename T::FERequirement& req,
                                         typename T::AffordanceCollectionType affordance, DBCOption dbcOption) {
                                  {
                                    t.vector(req, affordance.vectorAffordance(), dbcOption)
                                  } -> std::convertible_to<const Eigen::VectorXd&>;
                                  { t.vector(dbcOption) } -> std::convertible_to<const Eigen::VectorXd&>;
                                  { t.vector() } -> std::convertible_to<const Eigen::VectorXd&>;
                                };
 
  template <typename T>
  concept MatrixFlatAssembler = Concepts::VectorFlatAssembler<T> and
                                requires(T t, const typename T::FERequirement& req,
                                         typename T::AffordanceCollectionType affordance, DBCOption dbcOption) {
                                  { t.matrix(req, affordance.matrixAffordance(), dbcOption) };
                                  { t.matrix(dbcOption) };
                                  { t.matrix() };
                                };
 
  // adapted from /dune/dune-vtk/dune/vtk/utility/concepts.hh
  template <class DC>
  concept DataCollector = requires(DC dc) {
    typename DC::GridView;
    { dc.update() } -> std::same_as<void>;
    { dc.numPoints() } -> std::convertible_to<std::uint64_t>;
    { dc.numCells() } -> std::convertible_to<std::uint64_t>;
    { dc.gridView() } -> std::same_as<const typename DC::GridView&>;
  };
 
  template <class GV>
  concept GridView = requires(GV g) {
    typename GV::Grid;
    GV::dimension;
    GV::dimensionworld;
 
    { g.grid() };
  };
 
  namespace Impl {
    template <typename T>
    struct is_dual : std::false_type
    {
    };
 
    // Specialization for Dual<T, U>: this will be true for Dual types
    template <typename T, typename U>
    struct is_dual<autodiff::detail::Dual<T, U>> : std::true_type
    {
    };
  } // namespace Impl
 
  template <typename T>
  concept AutodiffScalar = Impl::is_dual<T>::value;
 
  template <typename T>
  concept SmartPointer = traits::isSharedPtr<T>::value || traits::isUniquePtr<T>::value;
 
  template <typename T>
  concept PointerOrSmartPointer = SmartPointer<T> || traits::Pointer<T>;
 
  template <typename S>
  concept ControlRoutineState = requires(S s) {
    typename S::Domain;
 
    { s.information } -> std::same_as<const ControlInformation&>;
    { s.loadStep } -> std::convertible_to<int>;
    { s.stepSize } -> std::convertible_to<double>;
    { s.domain } -> std::same_as<const typename S::Domain&>;
  };
 
  template <typename S>
  concept NonLinearSolverState = requires(S s) {
    typename S::Domain;
    typename S::CorrectionType;
 
    { s.information } -> std::same_as<const NonLinearSolverInformation&>;
    { s.domain } -> std::same_as<const typename S::Domain&>;
    { s.correction } -> std::same_as<const typename S::CorrectionType&>;
  };
 
  template <StrainTags tag>
  concept ReferenceConfiguraionStrain = (tag == StrainTags::greenLagrangian) or
                                        (tag == StrainTags::rightCauchyGreenTensor) or (tag == StrainTags::linear);
 
  template <StressTags tag>
  concept ReferenceConfiguraionStress = (tag == StressTags::PK2) or (tag == StressTags::linear);
 
  namespace Formulations {
    template <StrainTags T1, StressTags T2>
    concept TotalLagrangian = (T1 == StrainTags::linear and T2 == StressTags::linear) or
                              (T1 == StrainTags::greenLagrangian and T2 == StressTags::PK2);
 
    template <StrainTags T1, StressTags T2>
    concept TwoPoint =
        (T1 == StrainTags::deformationGradient or T1 == StrainTags::displacementGradient) and T2 == StressTags::PK1;
  } // namespace Formulations
 
  template <typename M>
  concept Material = requires(M m, Eigen::Matrix3d C) {
    M::strainTag;
    M::stressTag;
    M::tangentModuliTag;
    M::derivativeFactor;
    M::isReduced;
 
    typename M::ScalarType;
 
    { m.name() } -> std::convertible_to<std::string>;
    { m.materialParameters() } -> std::same_as<typename M::MaterialParameters>;
    { m.template storedEnergy<M::strainTag>(C) } -> std::same_as<typename M::ScalarType>;
    { m.template stresses<M::strainTag>(C) };
    { m.template tangentModuli<M::strainTag>(C) };
  };
 
  template <typename MAT>
  concept GeometricallyLinearMaterial = Material<MAT> && (MAT::strainTag == StrainTags::linear);
 
  template <typename ES>
  concept EigenValueSolver = requires(ES es) {
    typename ES::MatrixType;
    typename ES::ScalarType;
    { es.compute() } -> std::same_as<bool>;
    { es.eigenvalues() } -> std::convertible_to<Eigen::VectorX<typename ES::ScalarType>>;
    { es.eigenvectors() } -> std::convertible_to<Eigen::MatrixX<typename ES::ScalarType>>;
    { es.nev() } -> std::convertible_to<int>;
  };
 
} // namespace Concepts
 
namespace traits {
#ifndef DOXYGEN
 
  // Primary template: for non-pointer types.
  template <typename T, bool = Concepts::PointerOrSmartPointer<T>>
  struct MaybeDereference
  {
    using type = T;
  };
 
  // Specialization: for pointer (or smart pointer) types.
  template <typename T>
  struct MaybeDereference<T, true>
  {
    // We know T is pointer-like, so *std::declval<T&>() is well-formed.
    using type = std::remove_reference_t<decltype(*std::declval<T&>())>;
  };
#endif
 
  template <typename T>
  using MaybeDereferencedType = typename MaybeDereference<T>::type;
 
} // namespace traits
 
} // namespace Ikarus