finiteelements/mechanics/materials/interface.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 <ikarus/finiteelements/mechanics/materials/materialhelpers.hh>
#include <ikarus/finiteelements/mechanics/materials/numericalmaterialinversion.hh>
#include <ikarus/finiteelements/mechanics/materials/strainconversions.hh>
#include <ikarus/finiteelements/mechanics/materials/tags.hh>
#include <ikarus/finiteelements/physicshelper.hh>
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/linearalgebrahelper.hh>
#include <ikarus/utils/traits.hh>
 
namespace Ikarus::Materials {
 
#ifndef DOXYGEN
template <class MImpl>
struct Material;
 
template <auto stressIndexPair, typename MImpl>
struct VanishingStress;
 
template <auto strainIndexPair, typename MImpl>
struct VanishingStrain;
#endif
 
template <typename MAT, typename S>
consteval bool hasCorrectSize() {
  if constexpr (Concepts::EigenVector6<S> or Concepts::EigenMatrix33<S>)
    return true;
  if constexpr (MAT::isReduced and Concepts::EigenVector<S>) {
    return S::RowsAtCompileTime == MAT::freeStrains;
  } else
    return false;
}
 
template <typename MAT, typename S>
concept CorrectStrainSize = hasCorrectSize<MAT, S>();
 
template <typename MAT, bool isReduced = MAT::isReduced>
struct UnderlyingMaterial;
 
template <typename MAT>
struct UnderlyingMaterial<MAT, true>
{
  using type = MAT::Underlying;
  static auto material(const MAT& mat) { return mat.underlying(); }
};
 
template <typename MAT>
struct UnderlyingMaterial<MAT, false>
{
  using type = MAT;
  static auto material(const MAT& mat) { return mat; }
};
 
template <typename MAT>
using UnderlyingMaterial_t = typename UnderlyingMaterial<MAT>::type;
 
template <typename MAT>
auto underlyingMaterial(const MAT& mat) {
  return UnderlyingMaterial<MAT>::material(mat);
}
 
template <class MI>
struct Material
{
  using MaterialImpl = MI; 
 
  static constexpr bool isStrainVanished =
      traits::isSpecializationNonTypeAndTypes<VanishingStrain, MaterialImpl>::value;
  static constexpr bool isStressVanished =
      traits::isSpecializationNonTypeAndTypes<VanishingStress, MaterialImpl>::value;
 
  static constexpr bool isReduced = isStrainVanished or isStressVanished;
 
  static constexpr bool isLinear = MI::strainTag == StrainTags::linear;
 
  static constexpr bool isHyperelastic = [] {
    if constexpr (isReduced) {
      return requires {
        typename MI::Underlying::DeviatoricType;
        typename MI::Underlying::VolumetricType;
      };
    } else {
      return requires {
        typename MI::DeviatoricType;
        typename MI::VolumetricType;
      };
    }
  }();
 
  constexpr const MaterialImpl& impl() const { return static_cast<const MaterialImpl&>(*this); }
 
  constexpr MaterialImpl& impl() { return static_cast<MaterialImpl&>(*this); }
 
  [[nodiscard]] constexpr static std::string name() { return MI::nameImpl(); }
 
  [[nodiscard]] auto materialParameters() const { return impl().materialParametersImpl(); }
 
  static constexpr double derivativeFactor = MI::derivativeFactorImpl;
 
  template <StrainTags tag, typename Derived>
  requires CorrectStrainSize<MaterialImpl, Derived>
  [[nodiscard]] auto storedEnergy(const Eigen::MatrixBase<Derived>& Eraw) const {
    decltype(auto) Ev = enlargeIfReduced<Material>(Eraw);
    decltype(auto) E  = transformStrain<tag, MaterialImpl::strainTag>(Ev);
 
    if constexpr (Concepts::EigenVector<Derived>) { // receiving vector means voigt notation
      if constexpr (MaterialImpl::energyAcceptsVoigt)
        return impl().storedEnergyImpl(toVoigt(E));
      else
        return impl().storedEnergyImpl(E);
    } else
      return impl().storedEnergyImpl(E);
  }
 
  template <StrainTags tag, bool voigt = true, typename Derived>
  requires CorrectStrainSize<MaterialImpl, Derived>
  [[nodiscard]] auto stresses(const Eigen::MatrixBase<Derived>& Eraw) const {
    decltype(auto) Ev = enlargeIfReduced<Material>(Eraw);
    decltype(auto) E  = transformStrain<tag, MaterialImpl::strainTag>(Ev);
    if constexpr (voigt and MaterialImpl::stressToVoigt == false)
      // user requests a Voigt shaped return but material is not able to. Therefore, we transform it here.
      return toVoigt(stressesMaybeTransformInputToVoigt<false>(E), false);
    else
      return stressesMaybeTransformInputToVoigt<voigt>(E);
  }
 
  template <StrainTags tag, bool voigt = true, typename Derived>
  requires CorrectStrainSize<MaterialImpl, Derived>
  [[nodiscard]] auto tangentModuli(const Eigen::MatrixBase<Derived>& Eraw) const {
    decltype(auto) Ev = enlargeIfReduced<Material>(Eraw);
    decltype(auto) E  = transformStrain<tag, MaterialImpl::strainTag>(Ev);
    if constexpr (voigt and MaterialImpl::moduliToVoigt == false)
      // user request a Voigt shaped return but material is not able to. Therefore, we transform it here.
      return toVoigt(tangentModuliMaybeTransformInputToVoigt<false>(E));
    else
      return tangentModuliMaybeTransformInputToVoigt<voigt>(E);
  }
 
  template <StrainTags tag, bool voigt = true, bool useNumeric = false, typename Derived>
  requires CorrectStrainSize<MaterialImpl, Derived>
  [[nodiscard]] auto materialInversion(const Eigen::MatrixBase<Derived>& Sraw,
                                       const Eigen::MatrixBase<Derived>& EstartRaw = Derived::Zero().eval(),
                                       const double tol = 1e-12, const int maxIter = 20) const {
    const auto S = Impl::maybeFromVoigt(Sraw.derived(), false).eval();
 
    auto [D, Eraw] = [&]() {
      if constexpr (requires { impl().materialInversionImpl(S); } and not useNumeric)
        return impl().materialInversionImpl(S);
      else {
        const auto Estart = Impl::maybeFromVoigt(EstartRaw.derived(), true).eval();
        return numericalMaterialInversion(impl(), S, Estart, tol, maxIter);
      }
    }();
 
    const auto E = transformStrain<MaterialImpl::strainTag, tag>(Eraw).eval();
    if constexpr (voigt)
      return std::make_pair(D, toVoigt(E));
    else
      return std::make_pair(fromVoigt(D), E);
  }
 
  template <typename STO>
  auto rebind() const {
    return impl().template rebind<STO>();
  }
 
private:
  template <bool voigt = true, typename Derived>
  auto stressesMaybeTransformInputToVoigt(const Eigen::MatrixBase<Derived>& E) const {
    if constexpr (Concepts::EigenVector<Derived>) { // receiving vector means Voigt notation
      if constexpr (MaterialImpl::stressAcceptsVoigt)
        return impl().template stressesImpl<voigt>(E);
      else // material is not able to accept Voigt shaped Input. Therefore, we transform it before.
        return impl().template stressesImpl<voigt>(fromVoigt(E.derived()));
    } else
      return impl().template stressesImpl<voigt>(E.derived());
  }
 
  template <bool voigt = true, typename Derived>
  auto tangentModuliMaybeTransformInputToVoigt(const Eigen::MatrixBase<Derived>& E) const {
    if constexpr (Concepts::EigenVector<Derived>) { // receiving vector means voigt notation
      if constexpr (MaterialImpl::moduliAcceptsVoigt)
        return impl().template tangentModuliImpl<voigt>(E);
      else
        return impl().template tangentModuliImpl<voigt>(fromVoigt(E.derived()));
    } else
      return impl().template tangentModuliImpl<voigt>(E.derived());
  }
};
 
} // namespace Ikarus::Materials