ogden.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@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/tags.hh>
#include <ikarus/utils/tensorutils.hh>
 
namespace Ikarus::Materials {
 
template <typename ST_, int n, PrincipalStretchTags tag>
struct OgdenT
{
  using ScalarType = ST_;
 
  template <typename ST = ScalarType>
  using PrincipalStretches = Eigen::Vector<ST, 3>;
 
  static constexpr PrincipalStretchTags stretchTag = tag;
  static constexpr int numMatParameters            = n;
  static constexpr int dim                         = 3;
  static constexpr bool usesDeviatoricStretches    = stretchTag == PrincipalStretchTags::deviatoric;
 
  template <typename ST = ScalarType>
  using FirstDerivative = Eigen::Vector<ST, dim>;
  template <typename ST = ScalarType>
  using SecondDerivative = Eigen::Matrix<ST, dim, dim>;
 
  using MaterialParameters = std::array<double, numMatParameters>;
  using MaterialExponents  = std::array<double, numMatParameters>;
 
  [[nodiscard]] constexpr static std::string name() noexcept {
    return "Ogden (n = " + std::to_string(numMatParameters) + ", stretch type = " + toString(tag) + ")";
  }
 
  explicit OgdenT(const MaterialParameters& mpt, const MaterialExponents& mex)
      : materialParameters_{mpt},
        materialExponents_{mex} {}
 
  MaterialParameters materialParametersImpl() const { return materialParameters_; }
 
  const MaterialExponents& materialExponents() const { return materialExponents_; }
 
  template <typename ST>
  auto storedEnergyImpl(const PrincipalStretches<ST>& lambda) const {
    auto& mu    = materialParameters_;
    auto& alpha = materialExponents_;
 
    ST energy{};
 
    if constexpr (usesDeviatoricStretches) {
      auto lambdaBar = Impl::deviatoricStretches(lambda);
      for (auto i : parameterRange())
        energy += mu[i] / alpha[i] * (lambdaBar.array().pow(alpha[i]).sum() - 3);
    } else {
      auto J    = lambda[0] * lambda[1] * lambda[2];
      auto logJ = log(J);
      for (auto i : parameterRange())
        energy += mu[i] / alpha[i] * (lambda.array().pow(alpha[i]).sum() - 3) - mu[i] * logJ;
    }
    return energy;
  }
 
  template <typename ST>
  auto firstDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
    auto& mu       = materialParameters_;
    auto& alpha    = materialExponents_;
    auto dWdLambda = FirstDerivative<ST>::Zero().eval();
 
    if constexpr (usesDeviatoricStretches) {
      auto lambdaBar    = Impl::deviatoricStretches(lambda);
      auto dWdLambdaBar = Eigen::Array<ST, dim, 1>::Zero().eval();
 
      for (const auto j : parameterRange())
        dWdLambdaBar += mu[j] * lambdaBar.array().pow(alpha[j] - 1);
 
      const ST sumLambdaBar = (lambdaBar.array() * dWdLambdaBar).sum();
      dWdLambda             = (lambdaBar.array() * dWdLambdaBar - (1.0 / 3.0) * sumLambdaBar) / lambda.array();
    } else
      for (const auto j : parameterRange())
        dWdLambda.array() += (mu[j] * (lambda.array().pow(alpha[j]) - 1)) / lambda.array();
 
    return dWdLambda;
  }
 
  template <typename ST>
  auto secondDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
    auto& mu    = materialParameters_;
    auto& alpha = materialExponents_;
    auto dS     = SecondDerivative<ST>::Zero().eval();
 
    if constexpr (usesDeviatoricStretches) {
      const auto lambdaBar = Impl::deviatoricStretches(lambda);
      const auto dWdLambda = firstDerivativeImpl(lambda);
 
      auto lambdaBarPowSum = Eigen::Array<ST, dim, 1>::Zero().eval();
      for (const auto p : parameterRange())
        lambdaBarPowSum[p] = lambdaBar.array().pow(alpha[p]).sum();
 
      for (const auto a : dimensionRange()) {
        for (const auto b : dimensionRange()) {
          if (a == b)
            for (const auto p : parameterRange())
              dS(a, b) += mu[p] * alpha[p] * (1.0 / 3.0 * pow(lambdaBar[a], alpha[p]) + 1.0 / 9.0 * lambdaBarPowSum[p]);
 
          else
            for (const auto p : parameterRange())
              dS(a, b) += mu[p] * alpha[p] *
                          (-(1.0 / 3.0) * (pow(lambdaBar[a], alpha[p]) + pow(lambdaBar[b], alpha[p])) +
                           1.0 / 9.0 * lambdaBarPowSum[p]);
 
          dS(a, b) *= 1.0 / (lambda[a] * lambda[b]);
          if (a == b)
            dS(a, b) -= (2.0 / lambda[a]) * dWdLambda[a];
        }
      }
    } else {
      for (const auto j : parameterRange()) {
        dS.diagonal().array() +=
            (-2 * mu[j] * (lambda.array().pow(alpha[j]) - 1) + mu[j] * lambda.array().pow(alpha[j]) * alpha[j]) /
            lambda.array().square();
      }
    }
    return dS;
  }
 
  template <typename STO>
  auto rebind() const {
    return OgdenT<STO, numMatParameters, stretchTag>(materialParameters_, materialExponents_);
  }
 
private:
  MaterialParameters materialParameters_;
  MaterialExponents materialExponents_;
 
  inline static constexpr auto parameterRange() { return Dune::range(numMatParameters); }
  inline static constexpr auto dimensionRange() { return Dune::range(dim); }
};
 
template <int n, PrincipalStretchTags tag>
using Ogden = OgdenT<double, n, tag>;
 
} // namespace Ikarus::Materials