neohooke.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <ikarus/finiteelements/mechanics/materials/interface.hh>
#include <ikarus/utils/tensorutils.hh>
 
namespace Ikarus {
 
template <typename ST>
struct NeoHookeT : public Material<NeoHookeT<ST>>
{
  using ScalarType                    = ST;
  static constexpr int worldDimension = 3;
  using StrainMatrix                  = Eigen::Matrix<ScalarType, worldDimension, worldDimension>;
  using StressMatrix                  = StrainMatrix;
  using MaterialParameters            = LamesFirstParameterAndShearModulus;
 
  static constexpr auto strainTag              = StrainTags::rightCauchyGreenTensor;
  static constexpr auto stressTag              = StressTags::PK2;
  static constexpr auto tangentModuliTag       = TangentModuliTags::Material;
  static constexpr bool energyAcceptsVoigt     = false;
  static constexpr bool stressToVoigt          = false;
  static constexpr bool stressAcceptsVoigt     = false;
  static constexpr bool moduliToVoigt          = false;
  static constexpr bool moduliAcceptsVoigt     = false;
  static constexpr double derivativeFactorImpl = 2;
 
  [[nodiscard]] constexpr static std::string nameImpl() noexcept { return "NeoHooke"; }
 
  explicit NeoHookeT(const MaterialParameters& mpt)
      : materialParameter_{mpt} {}
 
  MaterialParameters materialParametersImpl() const { return materialParameter_; }
 
  template <typename Derived>
  ScalarType storedEnergyImpl(const Eigen::MatrixBase<Derived>& C) const {
    static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
    if constexpr (!Concepts::EigenVector<Derived>) {
      const auto traceC = C.trace();
      const auto detC   = C.determinant();
      checkPositiveDetC(detC);
      const auto logdetF = log(sqrt(detC));
      return materialParameter_.mu / 2.0 * (traceC - 3 - 2 * logdetF) +
             materialParameter_.lambda / 2.0 * logdetF * logdetF;
    } else
      static_assert(!Concepts::EigenVector<Derived>,
                    "NeoHooke energy can only be called with a matrix and not a vector in Voigt notation");
  }
 
  template <bool voigt, typename Derived>
  auto stressesImpl(const Eigen::MatrixBase<Derived>& C) const {
    static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
    if constexpr (!voigt) {
      if constexpr (!Concepts::EigenVector<Derived>) {
        const auto detC = C.determinant();
        checkPositiveDetC(detC);
        const auto logdetF = log(sqrt(detC));
        const auto invC    = C.inverse().eval();
        return (materialParameter_.mu * (StrainMatrix::Identity() - invC) + materialParameter_.lambda * logdetF * invC)
            .eval();
      } else
        static_assert(!Concepts::EigenVector<Derived>,
                      "NeoHooke can only be called with a matrix and not a vector in Voigt notation");
    } else
      static_assert(voigt == false, "NeoHooke does not support returning stresses in Voigt notation");
  }
 
  template <bool voigt, typename Derived>
  auto tangentModuliImpl(const Eigen::MatrixBase<Derived>& C) const {
    static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
    if constexpr (!voigt) {
      const auto invC = C.inverse().eval();
      const auto detC = C.determinant();
      checkPositiveDetC(detC);
      const auto logdetF = log(sqrt(detC));
      const auto CTinv   = tensorView(invC, std::array<Eigen::Index, 2>({3, 3}));
      static_assert(Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3>>::NumIndices == 2);
      Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>> moduli =
          (materialParameter_.lambda * dyadic(CTinv, CTinv) +
           2 * (materialParameter_.mu - materialParameter_.lambda * logdetF) *
               symTwoSlots(fourthOrderIKJL(invC, invC), {2, 3}))
              .eval();
      return moduli;
    } else
      static_assert(voigt == false, "NeoHooke does not support returning tangent moduli in Voigt notation");
  }
 
  template <typename STO>
  auto rebind() const {
    return NeoHookeT<STO>(materialParameter_);
  }
 
private:
  MaterialParameters materialParameter_;
 
  void checkPositiveDetC(ScalarType detC) const {
    if (Dune::FloatCmp::le(static_cast<double>(detC), 0.0, 1e-10))
      DUNE_THROW(Dune::InvalidStateException,
                 "Determinant of right Cauchy Green tensor C must be greater than zero. detC = " +
                     std::to_string(static_cast<double>(detC)));
  }
};
 
using NeoHooke = NeoHookeT<double>;
 
} // namespace Ikarus