HMHStorakers.cpp

Go to the documentation of this file.

/**
 * @file HMHStorakers.cpp
 * @brief Storakers material model
 *
 * The model can either approximate the full natural stretch H = log(U), or
 * recover U locally from the Biot stress in the no-stretch formulation.
 */
 
#include <Lie.hpp>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdio>
#include <iomanip>
#include <limits>
#include <regex>
#include <sstream>
#include <utility>
 
namespace EshelbianPlasticity {
 
struct HMHStorakers : public PhysicalEquations {
  struct MaterialParameters {
    double eta;
    double mu;
    double beta;
  };
 
  struct DerivedParameters {
    double bulkModulus;
    double poissonRatio;
    double lameLambda;
  };
 
  HMHStorakers(MoFEM::Interface &m_field, const double eta, const double mu,
               const double beta)
      : mField(m_field), defaultParameters{eta, mu, beta} {
    CHK_THROW_MESSAGE(getOptions(), "get options failed");
    CHK_THROW_MESSAGE(extractBlockData(Sev::inform),
                      "extract block data failed");
 
    if (EshelbianCore::isNoStretch()) {
      EshelbianCore::stretchHandling = EshelbianCore::NO_STREACH_NONLINEAR;
    }
    if (EshelbianCore::stretchSelector != StretchSelector::LOG) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Storakers requires natural logarithmic stretch");
    }
  }
 
  static DerivedParameters getDerivedParameters(const double mu,
                                                const double beta) {
    return {
        2 * mu * (beta + 1. / 3.),
        beta / (1 + 2 * beta),
        2 * beta * mu,
    };
  }
 
  static inline double getMu1(const MaterialParameters &parameters) {
    return (1 - parameters.eta) * parameters.mu;
  }
 
  static inline double getMu2(const MaterialParameters &parameters) {
    return parameters.eta * parameters.mu;
  }
 
  template <typename T>
  static inline double getLogJacobian(T &principal_hencky_stretches) {
    FTENSOR_INDEX(SPACE_DIM, i);
    constexpr auto t_one = FTensor::One<>();
    return principal_hencky_stretches(i) * t_one(i);
  }
 
  template <typename T>
  static inline double getJacobianPower(
      const MaterialParameters &parameters,
      T &principal_hencky_stretches) {
    return std::exp(-parameters.beta *
                    getLogJacobian(principal_hencky_stretches));
  }
 
  static inline double getGamma(const MaterialParameters &parameters,
                                const double jacobian_power) {
    const double mu1 = getMu1(parameters);
    const double mu2 = getMu2(parameters);
    return 2 * mu1 * jacobian_power +
           mu2 * jacobian_power * jacobian_power;
  }
 
  static inline double getGammaY(const MaterialParameters &parameters,
                                 const double jacobian_power) {
    const double mu1 = getMu1(parameters);
    const double mu2 = getMu2(parameters);
    return 2 * mu1 * jacobian_power +
           2 * mu2 * jacobian_power * jacobian_power;
  }
 
  static inline double getPrincipalCoordinateStress(
      const MaterialParameters &parameters, const double gamma,
      const double hencky_stretch) {
    const double stretch = std::exp(hencky_stretch);
    return 2 * getMu1(parameters) * stretch +
           getMu2(parameters) * stretch * stretch - gamma;
  }
 
  static inline double getPrincipalCoordinateStressDerivativeAtFixedGamma(
      const MaterialParameters &parameters, const double hencky_stretch) {
    const double stretch = std::exp(hencky_stretch);
    return 2 * getMu1(parameters) * stretch +
           2 * getMu2(parameters) * stretch * stretch;
  }
 
  MaterialParameters getMaterialParameters(const EntityHandle ent) const {
    for (const auto &block : blockData) {
      if (block.blockEnts.find(ent) != block.blockEnts.end()) {
        return {block.eta, block.mu, block.beta};
      }
    }
 
    if (!blockData.empty()) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Storakers material block not found for entity. If "
                        "blocks are used, assign all material entities");
    }
 
    return defaultParameters;
  }
 
  struct OpJacobian : public EshelbianPlasticity::OpJacobian {
    using EshelbianPlasticity::OpJacobian::OpJacobian;
 
    MoFEMErrorCode evaluateRhs(EntData &) override { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &) override { return 0; }
  };
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) override {
    return new OpJacobian(eval_rhs, eval_lhs, data_ptr, physics_ptr);
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "storakers_", "Storakers material",
                      "none");
    CHKERR PetscOptionsScalar("-eta", "Blend fraction eta", "",
                              defaultParameters.eta, &defaultParameters.eta,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-mu", "Infinitesimal shear modulus mu", "",
                              defaultParameters.mu, &defaultParameters.mu,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-beta", "Compressibility parameter beta", "",
                              defaultParameters.beta, &defaultParameters.beta,
                              PETSC_NULLPTR);
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    CHKERR validateParameters(defaultParameters, "default options");
 
    const auto derived =
        getDerivedParameters(defaultParameters.mu, defaultParameters.beta);
    MOFEM_LOG("EP", Sev::inform)
        << "MatBlock Storakers (default): "
        << "eta = " << defaultParameters.eta
        << ", mu = " << defaultParameters.mu
        << ", beta = " << defaultParameters.beta
        << ", K = " << derived.bulkModulus << ", nu = " << derived.poissonRatio
        << ", lambda_L = " << derived.lameLambda
        << ", grad alpha u = " << alphaGradU;
    MOFEM_LOG("EP", Sev::inform)
        << "Storakers conversions: mu_1 = (1 - eta)*mu; mu_2 = eta*mu; "
           "lambda_L = 2*beta*mu; K = 2*mu*(beta + 1/3); "
           "nu = beta/(1 + 2*beta); beta = nu/(1 - 2*nu)";
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(const Sev sev) {
    return extractBlockData(
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % "MAT_STORAKES").str())),
        sev);
  }
 
  MoFEMErrorCode
  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
                   const Sev sev) {
    MoFEMFunctionBegin;
 
    for (const auto meshset : meshset_vec_ptr) {
      MOFEM_LOG("EP", sev) << *meshset;
 
      std::vector<double> attributes;
      CHKERR meshset->getAttributes(attributes);
 
      const auto json_parameters =
          mField.getInterface<JsonConfigManager>()->getParamsFromBlockset(
              "MAT_STORAKES", meshset->getMeshsetId());
 
      MaterialParameters parameters{};
      if (!json_parameters.empty()) {
        if (json_parameters.size() != 3 ||
            json_parameters.find("eta") == json_parameters.end() ||
            json_parameters.find("mu") == json_parameters.end() ||
            json_parameters.find("beta") == json_parameters.end()) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "MAT_STORAKES JSON block must have exactly three "
                  "attributes: eta, mu, beta");
        }
        parameters = {json_parameters.at("eta"), json_parameters.at("mu"),
                      json_parameters.at("beta")};
      } else {
        if (attributes.size() != 3) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "MAT_STORAKES block must have exactly three attributes: "
                  "eta, mu, beta");
        }
        parameters = {attributes[0], attributes[1], attributes[2]};
      }
      CHKERR validateParameters(parameters, "MAT_STORAKES block");
 
      Range entities;
      CHKERR mField.get_moab().get_entities_by_handle(meshset->meshset,
                                                      entities, true);
      blockData.push_back(
          {parameters.eta, parameters.mu, parameters.beta, entities});
 
      const auto derived = getDerivedParameters(parameters.mu, parameters.beta);
      MOFEM_LOG("EP", sev) << "MatBlock Storakers: eta = " << parameters.eta
                           << ", mu = " << parameters.mu
                           << ", beta = " << parameters.beta
                           << ", K = " << derived.bulkModulus
                           << ", nu = " << derived.poissonRatio
                           << ", lambda_L = " << derived.lameLambda
                           << ", nb. entities = " << entities.size();
    }
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * @brief Assemble the full-order equation conjugate to H = log(U)
   *
   * The material part of the coordinate stress is
   * \f[
   * \boldsymbol{\tau}_H = \mathbf U\mathbf T
   * = 2\mu_1\mathbf U + \mu_2\mathbf U^2 - \gamma\mathbf I.
   * \f]
   */
  struct OpSpatialPhysical : public OpAssembleVolume {
    OpSpatialPhysical(const std::string &field_name,
                      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                      const double alpha_u);
 
    MoFEMErrorCode integrate(EntData &data) override;
 
  private:
    const double alphaU;
  };
 
  VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) override {
    return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
  }
 
  /**
   * @brief Assemble the consistent full-order H--H tangent
   *
   * \f[
   * \mathbb C_H = D_H(2\mu_1e^H + \mu_2e^{2H})
   * + \beta\gamma_Y\mathbf I\otimes\mathbf I.
   * \f]
   */
  struct OpSpatialPhysical_du_du : public OpAssembleVolumePositiveDefine {
    OpSpatialPhysical_du_du(std::string row_field, std::string col_field,
                            boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                            const double alpha_u);
 
    MoFEMErrorCode getOptions();
    MoFEMErrorCode integrate(EntData &row_data, EntData &col_data) override;
 
  private:
    const double alphaU;
    double minimEigenValue = 0;
  };
 
  VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const double alpha_u) override {
    return new OpSpatialPhysical_du_du(
        std::move(row_field), std::move(col_field), std::move(data_ptr),
        alpha_u);
  }
 
  /**
   * @brief Recover the right stretch from a symmetric Biot stress.
   *
   * This class inverts the Storakers constitutive equation
   *
   * \f[
   * \mathbf{T} = 2\mu_1\mathbf{I} + \mu_2\mathbf{U}
   *               - (\gamma+q)\mathbf{U}^{-1}, \qquad
   * \mu_1 = (1-\eta)\mu, \quad \mu_2 = \eta\mu,
   * \f]
   *
   * by diagonalising
   * \f$\mathbf{T}=\mathbf{Q}\operatorname{diag}(t_i)\mathbf{Q}^T\f$.
   * For a trial \f$y=\log z\f$, where \f$z=J^{-\beta}\f$, it evaluates
   *
   * \f[
   * \gamma = 2\mu_1 z + \mu_2 z^2, \qquad
   * \gamma_{\mathrm{eff}}=\gamma+q, \qquad
   * u_i(y) =
   * \frac{b_i+\sqrt{b_i^2+4\mu_2\gamma_{\mathrm{eff}}}}
   *      {2\mu_2}, \qquad
   * b_i=t_i-2\mu_1.
   * \f]
   * The positive-tangent \f$+\f$ root is used. If
   * \f$\gamma_{\mathrm{eff}}<0\f$, it is admissible only while
   * \f$b_i>0\f$ and
   * \f$b_i^2+4\mu_2\gamma_{\mathrm{eff}}>0\f$ for every principal value.
   *
   * The only nonlinear unknown is obtained with safeguarded Newton--Brent
   * iterations from
   *
   * \f[
   * \phi(y)=\frac{y}{\beta}+\sum_i\log u_i(y)=0.
   * \f]
   *
   * Finally,
   * \f$\mathbf{U}=\mathbf{Q}\operatorname{diag}(u_i)\mathbf{Q}^T\f$ and
   * \f$J=\exp(-y/\beta)\f$. The spelling "Strech" is retained to match the
   * existing material interface.
   */
  struct StressToStrech {
    struct ScalarEquationState {
      double y = 0;
      double z = 1;
      double gamma = 0;
      double gammaY = 0;
      double q = 0;
      double gammaEffective = 0;
      int gammaEffectiveSign = 0;
      double logAbsGammaEffective =
          -std::numeric_limits<double>::infinity();
      double phi = 0;
      double dPhi = 0;
      FTensor::Tensor1<double, SPACE_DIM> tS{0., 0., 0.};
      FTensor::Tensor1<double, SPACE_DIM> tU{0., 0., 0.};
      FTensor::Tensor1<double, SPACE_DIM> tLogU{0., 0., 0.};
    };
 
    explicit StressToStrech(const MaterialParameters &parameters)
        : materialParameters(parameters),
          mu1((1 - parameters.eta) * parameters.mu),
          mu2(parameters.eta * parameters.mu) {
      CHK_THROW_MESSAGE(validateParameters(parameters, "StressToStrech"),
                        "invalid Storakers parameters");
    }
 
    /**
     * @brief Diagonalise sym(stress) and calculate b_i = t_i - 2 mu_1
     */
    template <typename T> MoFEMErrorCode stressToSpectralForm(T &t_stress) {
      return stressToSpectralForm(t_stress, 0.);
    }
 
    template <typename T>
    MoFEMErrorCode stressToSpectralForm(T &t_stress, const double q) {
      MoFEMFunctionBegin;
 
      if (!std::isfinite(q)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Storakers external pressure q must be finite");
      }
      externalPressure = q;
 
      FTENSOR_INDEXES(SPACE_DIM, i, j);
      tEigenVectors(i, j) =
          0.5 * (t_stress(i, j) + t_stress(j, i));
 
      for (int ii = 0; ii != SPACE_DIM; ++ii) {
        for (int jj = 0; jj != SPACE_DIM; ++jj) {
          if (!std::isfinite(tEigenVectors(ii, jj))) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                    "Storakers Biot stress contains a non-finite value");
          }
        }
      }
 
      CHKERR computeEigenValuesSymmetric(tEigenVectors, tPrincipalStress);
      nbUniqueEigenvalues = getUniqNb<SPACE_DIM>(tPrincipalStress);
      if (nbUniqueEigenvalues < SPACE_DIM) {
        CHKERR sortEigenVals<SPACE_DIM>(tPrincipalStress, tEigenVectors);
      }
      tB(i) = tPrincipalStress(i) - 2 * mu1;
 
      hasSpectralForm = true;
      hasScalarSolution = false;
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Evaluate u_i(y), phi(y), and the exact Newton derivative
     */
    MoFEMErrorCode
    calculatePrincipalStretches(const double y,
                                ScalarEquationState &state) const {
      MoFEMFunctionBegin;
 
      if (!hasSpectralForm) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Call stressToSpectralForm before evaluating the scalar "
                "equation");
      }
      if (!std::isfinite(y)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Non-finite Storakers scalar iterate");
      }
 
      state = {};
      state.y = y;
      FTENSOR_INDEX(SPACE_DIM, i);
      constexpr auto t_one = FTensor::One<>();
      const double log_minimum =
          std::log(std::numeric_limits<double>::min());
      const double log_maximum =
          std::log(std::numeric_limits<double>::max());
      const double negative_infinity =
          -std::numeric_limits<double>::infinity();
      const double log_linear_gamma =
          mu1 > 0 ? std::log(2 * mu1) + y : negative_infinity;
      const double log_quadratic_gamma = std::log(mu2) + 2 * y;
      const double log_gamma =
          logSumExp(log_linear_gamma, log_quadratic_gamma);
      const double log_gamma_y = logSumExp(
          log_linear_gamma, std::log(2 * mu2) + 2 * y);
 
      int gamma_effective_sign = 1;
      double log_abs_gamma_effective = log_gamma;
      if (externalPressure > 0) {
        log_abs_gamma_effective =
            logSumExp(log_gamma, std::log(externalPressure));
      } else if (externalPressure < 0) {
        const double log_minus_q = std::log(-externalPressure);
        if (log_gamma > log_minus_q) {
          log_abs_gamma_effective =
              log_gamma +
              std::log(-std::expm1(log_minus_q - log_gamma));
        } else if (log_gamma < log_minus_q) {
          gamma_effective_sign = -1;
          log_abs_gamma_effective =
              log_minus_q +
              std::log(-std::expm1(log_gamma - log_minus_q));
        } else {
          gamma_effective_sign = 0;
          log_abs_gamma_effective = negative_infinity;
        }
      }
 
      if (y > log_maximum || log_gamma > log_maximum ||
          log_abs_gamma_effective > log_maximum ||
          log_gamma_y > log_maximum) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Storakers scalar iterate is outside the floating-point "
                "range");
      }
 
      state.z = y >= log_minimum ? std::exp(y) : 0;
      state.gamma = log_gamma >= log_minimum ? std::exp(log_gamma) : 0;
      state.gammaY =
          log_gamma_y >= log_minimum ? std::exp(log_gamma_y) : 0;
      state.q = externalPressure;
      state.gammaEffectiveSign = gamma_effective_sign;
      state.logAbsGammaEffective = log_abs_gamma_effective;
      if (gamma_effective_sign &&
          log_abs_gamma_effective >= log_minimum) {
        state.gammaEffective =
            gamma_effective_sign * std::exp(log_abs_gamma_effective);
      }
 
      const double log_two_mu2 = std::log(2.) + std::log(mu2);
      double positive_gamma_term = 0;
      if (gamma_effective_sign > 0) {
        const double log_gamma_term =
            std::log(2.) +
            0.5 * (std::log(mu2) + log_abs_gamma_effective);
        if (log_gamma_term > log_maximum) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "Storakers principal-stretch discriminant is outside the "
                  "floating-point range");
        }
        positive_gamma_term =
            log_gamma_term >= log_minimum ? std::exp(log_gamma_term) : 0;
      }
 
      double log_derivative_sum = negative_infinity;
      for (int ii = 0; ii != SPACE_DIM; ++ii) {
        double log_s = negative_infinity;
        if (gamma_effective_sign > 0) {
          state.tS(ii) = std::hypot(tB(ii), positive_gamma_term);
          log_s = std::log(state.tS(ii));
        } else if (gamma_effective_sign == 0) {
          if (tB(ii) <= 0) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                    "No positive stable Storakers stretch at gamma_eff=0");
          }
          state.tS(ii) = tB(ii);
          log_s = std::log(tB(ii));
        } else {
          if (tB(ii) <= 0) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                    "Negative Storakers gamma_eff requires positive shifted "
                    "principal stresses");
          }
          const double log_b = std::log(tB(ii));
          const double log_negative_term =
              std::log(4.) + std::log(mu2) + log_abs_gamma_effective;
          const double log_b_squared = 2 * log_b;
          if (!(log_negative_term < log_b_squared)) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                    "Storakers state is outside the stable principal-stretch "
                    "branch");
          }
          log_s =
              log_b +
              0.5 * std::log(-std::expm1(log_negative_term -
                                          log_b_squared));
          state.tS(ii) = log_s >= log_minimum ? std::exp(log_s) : 0;
        }
        if (!std::isfinite(state.tS(ii)) || state.tS(ii) <= 0) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "Invalid Storakers principal-stretch discriminant");
        }
 
        if (tB(ii) >= 0) {
          const double log_b =
              tB(ii) > 0 ? std::log(tB(ii)) : negative_infinity;
          state.tLogU(ii) =
              logSumExp(log_b, log_s) - log_two_mu2;
        } else {
          const double log_denominator =
              logSumExp(log_s, std::log(-tB(ii)));
          state.tLogU(ii) =
              std::log(2.) + log_abs_gamma_effective - log_denominator;
        }
 
        if (!std::isfinite(state.tLogU(ii)) ||
            state.tLogU(ii) < log_minimum ||
            state.tLogU(ii) > log_maximum) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "Invalid Storakers principal stretch");
        }
        state.tU(ii) = std::exp(state.tLogU(ii));
 
        const double log_derivative_term =
            -state.tLogU(ii) - log_s;
        log_derivative_sum =
            logSumExp(log_derivative_sum, log_derivative_term);
      }
 
      state.phi =
          y / materialParameters.beta + state.tLogU(i) * t_one(i);
      const double log_volumetric_derivative =
          log_gamma_y + log_derivative_sum;
      if (log_volumetric_derivative > log_maximum) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Invalid Storakers scalar derivative term");
      }
      const double volumetric_derivative =
          log_volumetric_derivative >= log_minimum
              ? std::exp(log_volumetric_derivative)
              : 0;
      state.dPhi = 1 / materialParameters.beta + volumetric_derivative;
      if (!std::isfinite(state.phi) || !std::isfinite(state.dPhi) ||
          state.dPhi <= 0) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Invalid Storakers scalar residual or derivative");
      }
 
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Solve phi(y)=0 with safeguarded Newton--Brent iterations
     *
     * The last converged value is used as the initial guess, which allows one
     * instance to be reused over a sequence of nearby quadrature-point states.
     */
    MoFEMErrorCode solveScalarEquation() {
      MoFEMFunctionBegin;
      const auto [safe_lower_y, safe_upper_y] = getSafeYBounds();
      if (!(safe_lower_y < safe_upper_y)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "No admissible Storakers scalar interval for q=%.16g",
                externalPressure);
      }
      CHKERR solveScalarEquation(
          std::clamp(scalarState.y, safe_lower_y, safe_upper_y));
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode solveScalarEquation(const double initial_y) {
      MoFEMFunctionBegin;
 
      if (!hasSpectralForm) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Call stressToSpectralForm before solveScalarEquation");
      }
 
      const auto [safe_lower_y, safe_upper_y] = getSafeYBounds();
      if (!(safe_lower_y < safe_upper_y)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "No admissible Storakers scalar interval for q=%.16g",
                externalPressure);
      }
      if (!std::isfinite(initial_y) || initial_y < safe_lower_y ||
          initial_y > safe_upper_y) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Initial Storakers scalar iterate %.16g is outside the safe "
                "interval [%.16g, %.16g]",
                initial_y, safe_lower_y, safe_upper_y);
      }
 
      ScalarEquationState initial;
      CHKERR calculatePrincipalStretches(initial_y, initial);
      if (std::abs(initial.phi) <= absoluteTolerance) {
        scalarState = initial;
        hasScalarSolution = true;
        MoFEMFunctionReturnHot(0);
      }
 
      double lower_y = initial_y;
      double upper_y = initial_y;
      ScalarEquationState lower = initial;
      ScalarEquationState upper = initial;
      double bracket_step = 1;
      bool bracketed = false;
      double bracket_endpoint_y = initial_y;
      constexpr double bracket_growth_factor = 1.5;
 
      for (int iteration = 0; iteration < maxBracketIterations; ++iteration) {
        const bool search_lower = initial.phi > 0;
        const double available_distance =
            search_lower ? bracket_endpoint_y - safe_lower_y
                         : safe_upper_y - bracket_endpoint_y;
        if (available_distance <= 0) {
          break;
        }
 
        const double increment = std::min(bracket_step, available_distance);
        const double candidate_y =
            bracket_endpoint_y + (search_lower ? -increment : increment);
        if (candidate_y == bracket_endpoint_y) {
          break;
        }
 
        ScalarEquationState candidate;
        CHKERR calculatePrincipalStretches(candidate_y, candidate);
        if (initial.phi > 0) {
          lower_y = candidate_y;
          lower = candidate;
          upper_y = initial_y;
          upper = initial;
          bracketed = lower.phi <= 0;
        } else {
          lower_y = initial_y;
          lower = initial;
          upper_y = candidate_y;
          upper = candidate;
          bracketed = upper.phi >= 0;
        }
        if (bracketed) {
          break;
        }
        bracket_endpoint_y = candidate_y;
        bracket_step *= bracket_growth_factor;
      }
 
      if (!bracketed) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Could not bracket the Storakers scalar root inside the safe "
                "interval [%.16g, %.16g]",
                safe_lower_y, safe_upper_y);
      }
 
      double current_y =
          std::abs(lower.phi) < std::abs(upper.phi) ? lower_y : upper_y;
      ScalarEquationState current =
          std::abs(lower.phi) < std::abs(upper.phi) ? lower : upper;
      double history_y = std::numeric_limits<double>::quiet_NaN();
      double history_phi = std::numeric_limits<double>::quiet_NaN();
 
      for (int iteration = 0; iteration < maxIterations; ++iteration) {
        if (hasConverged(current, lower_y, upper_y)) {
          scalarState = current;
          hasScalarSolution = true;
          MoFEMFunctionReturnHot(0);
        }
 
        const double width = upper_y - lower_y;
        const double guard = 0.05 * width;
        auto is_safely_bracketed = [&](const double candidate) {
          return std::isfinite(candidate) && candidate > lower_y + guard &&
                 candidate < upper_y - guard;
        };
 
        // First use the exact derivative from Section 3.2.
        double candidate_y = current_y - current.phi / current.dPhi;
        bool accepted = is_safely_bracketed(candidate_y) &&
                        std::abs(candidate_y - current_y) <= 0.5 * width;
 
        // Brent inverse quadratic interpolation, when three distinct residual
        // values are available.
        if (!accepted && std::isfinite(history_phi) && lower.phi != upper.phi &&
            lower.phi != history_phi && upper.phi != history_phi) {
          candidate_y =
              lower_y * upper.phi * history_phi /
                  ((lower.phi - upper.phi) * (lower.phi - history_phi)) +
              upper_y * lower.phi * history_phi /
                  ((upper.phi - lower.phi) * (upper.phi - history_phi)) +
              history_y * lower.phi * upper.phi /
                  ((history_phi - lower.phi) * (history_phi - upper.phi));
          accepted = is_safely_bracketed(candidate_y);
        }
 
        // Brent secant interpolation, followed by the bisection safeguard.
        if (!accepted && lower.phi != upper.phi) {
          candidate_y = upper_y - upper.phi * (upper_y - lower_y) /
                                      (upper.phi - lower.phi);
          accepted = is_safely_bracketed(candidate_y);
        }
        if (!accepted) {
          candidate_y = 0.5 * (lower_y + upper_y);
        }
 
        ScalarEquationState candidate;
        CHKERR calculatePrincipalStretches(candidate_y, candidate);
        history_y = current_y;
        history_phi = current.phi;
 
        if (candidate.phi <= 0) {
          lower_y = candidate_y;
          lower = candidate;
        } else {
          upper_y = candidate_y;
          upper = candidate;
        }
 
        current_y = candidate_y;
        current = candidate;
      }
 
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Storakers Newton--Brent solver did not converge in %d "
              "iterations",
              maxIterations);
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Reconstruct U = Q diag(u_i) Q^T in the working basis
     */
    template <typename T>
    MoFEMErrorCode spectralFormToStrech(T &t_stretch) const {
      MoFEMFunctionBegin;
 
      if (!hasScalarSolution) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Call solveScalarEquation before spectralFormToStrech");
      }
 
      FTENSOR_INDEXES(SPACE_DIM, i, j);
      auto identity = [](const double u) { return u; };
      auto t_reconstructed_stretch =
          EigenMatrix::getMat(scalarState.tU, tEigenVectors, identity);
      t_stretch(i, j) = t_reconstructed_stretch(i, j);
 
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Calculate C = dT/dU in the working/reference coordinates
     *
     * EigenMatrix::getDiffMat evaluates the derivative at fixed
     * gamma_eff=gamma+q and rotates it from the principal basis with the stored
     * eigenvectors. The rank-one term adds the variation of J, z and material
     * gamma; prescribed q has zero derivative.
     */
    template <typename T>
    MoFEMErrorCode calculateStiffness(T &t_stiffness) const {
      MoFEMFunctionBegin;
 
      if (!hasScalarSolution) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Call solveScalarEquation before calculateStiffness");
      }
 
      FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
      const double gamma = scalarState.gammaEffective;
      auto principal_stress = [this, gamma](const double u) {
        return 2 * mu1 + mu2 * u - gamma / u;
      };
      auto principal_stiffness = [this, gamma](const double u) {
        return mu2 + (gamma / u) / u;
      };
      auto t_frozen_gamma_stiffness = EigenMatrix::getDiffMat(
          scalarState.tU, tEigenVectors, principal_stress,
          principal_stiffness, nbUniqueEigenvalues);
 
      auto inverse_stretch = [](const double u) { return 1 / u; };
      auto t_inverse_stretch = EigenMatrix::getMat(
          scalarState.tU, tEigenVectors, inverse_stretch);
      const double alpha = materialParameters.beta * scalarState.gammaY;
 
      t_stiffness(i, j, k, l) =
          t_frozen_gamma_stiffness(i, j, k, l) +
          alpha * t_inverse_stretch(i, j) * t_inverse_stretch(k, l);
 
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Calculate S = dU/dT in the working/reference coordinates
     *
     * getDiffMat applied to the inverse principal stress law gives L^{-1}
     * already rotated from the principal basis. The second term is the exact
     * Sherman--Morrison correction.
     */
    template <typename T>
    MoFEMErrorCode calculateCompliance(T &t_compliance) const {
      MoFEMFunctionBegin;
 
      if (!hasScalarSolution) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Call solveScalarEquation before calculateCompliance");
      }
 
      FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
      const int gamma_sign = scalarState.gammaEffectiveSign;
      const double log_abs_gamma = scalarState.logAbsGammaEffective;
      auto inverse_state = [this, gamma_sign,
                            log_abs_gamma](const double stress) {
        const double b = stress - 2 * mu1;
        double log_s;
        if (gamma_sign > 0) {
          const double log_gamma_term =
              std::log(2.) +
              0.5 * (std::log(mu2) + log_abs_gamma);
          const double gamma_term = std::exp(log_gamma_term);
          const double s = std::hypot(b, gamma_term);
          log_s = std::log(s);
          if (b < 0) {
            const double log_denominator =
                logSumExp(log_s, std::log(-b));
            const double log_u =
                std::log(2.) + log_abs_gamma - log_denominator;
            return std::pair{std::exp(log_u), s};
          }
        } else if (gamma_sign == 0) {
          log_s = std::log(b);
        } else {
          const double log_b = std::log(b);
          const double log_ratio = std::log(4.) + std::log(mu2) +
                                   log_abs_gamma - 2 * log_b;
          log_s = log_b + 0.5 * std::log(-std::expm1(log_ratio));
        }
 
        const double log_b = std::log(b);
        const double log_u =
            logSumExp(log_b, log_s) - std::log(2.) - std::log(mu2);
        return std::pair{std::exp(log_u), std::exp(log_s)};
      };
      auto inverse_stress_law = [inverse_state](const double stress) {
        return inverse_state(stress).first;
      };
      auto inverse_stress_law_derivative =
          [inverse_state](const double stress) {
            const auto [u, s] = inverse_state(stress);
            return u / s;
      };
      auto t_inverse_l = EigenMatrix::getDiffMat(
          tPrincipalStress, tEigenVectors, inverse_stress_law,
          inverse_stress_law_derivative, nbUniqueEigenvalues);
 
      auto inverse_discriminant = [](const double s) { return 1 / s; };
      auto t_y = EigenMatrix::getMat(scalarState.tS, tEigenVectors,
                                     inverse_discriminant);
 
      const double alpha = materialParameters.beta * scalarState.gammaY;
      FTensor::Tensor1<double, SPACE_DIM> t_inverse_stretch_dot_y_terms;
      for (int aa = 0; aa != SPACE_DIM; ++aa) {
        t_inverse_stretch_dot_y_terms(aa) =
            1 / (scalarState.tU(aa) * scalarState.tS(aa));
      }
      constexpr auto t_one = FTensor::One<>();
      const double inverse_stretch_dot_y =
          t_inverse_stretch_dot_y_terms(i) * t_one(i);
      const double denominator = 1 + alpha * inverse_stretch_dot_y;
      if (!std::isfinite(denominator) || denominator <= 0) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Invalid Storakers compliance denominator");
      }
      const double correction = alpha / denominator;
 
      t_compliance(i, j, k, l) =
          t_inverse_l(i, j, k, l) - correction * t_y(i, j) * t_y(k, l);
 
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Calculate stiffness and compliance in reference coordinates
     *
     * The output arguments can be owning Ddg tensors or MoFEM Ddg views over
     * matD and matInvD storage.
     */
    template <typename TStiffness, typename TCompliance>
    MoFEMErrorCode calculateBiotTangent(TStiffness &t_stiffness,
                                        TCompliance &t_compliance) const {
      MoFEMFunctionBegin;
      CHKERR calculateStiffness(t_stiffness);
      CHKERR calculateCompliance(t_compliance);
      MoFEMFunctionReturn(0);
    }
 
    template <typename TStress, typename TStretch>
    MoFEMErrorCode calculate(TStress &t_stress, TStretch &t_stretch) {
      return calculate(t_stress, t_stretch, 0.);
    }
 
    template <typename TStress, typename TStretch>
    MoFEMErrorCode calculate(TStress &t_stress, TStretch &t_stretch,
                             const double q) {
      MoFEMFunctionBegin;
      CHKERR stressToSpectralForm(t_stress, q);
      CHKERR solveScalarEquation();
      CHKERR spectralFormToStrech(t_stretch);
      MoFEMFunctionReturn(0);
    }
 
    /**
     * @brief Calculate U(T), dT/dU, and dU/dT in reference coordinates
     */
    template <typename TStress, typename TStretch, typename TStiffness,
              typename TCompliance>
    MoFEMErrorCode calculate(TStress &t_stress, TStretch &t_stretch,
                             TStiffness &t_stiffness,
                             TCompliance &t_compliance) {
      return calculate(t_stress, t_stretch, t_stiffness, t_compliance, 0.);
    }
 
    template <typename TStress, typename TStretch, typename TStiffness,
              typename TCompliance>
    MoFEMErrorCode calculate(TStress &t_stress, TStretch &t_stretch,
                             TStiffness &t_stiffness,
                             TCompliance &t_compliance, const double q) {
      MoFEMFunctionBegin;
      CHKERR calculate(t_stress, t_stretch, q);
      CHKERR calculateBiotTangent(t_stiffness, t_compliance);
      MoFEMFunctionReturn(0);
    }
 
    const ScalarEquationState &getScalarEquationState() const {
      return scalarState;
    }
 
    double getJacobian() const {
      return std::exp(-scalarState.y / materialParameters.beta);
    }
 
    const FTensor::Tensor1<double, SPACE_DIM> &getPrincipalStress() const {
      return tPrincipalStress;
    }
 
    const FTensor::Tensor1<double, SPACE_DIM> &
    getShiftedPrincipalStress() const {
      return tB;
    }
 
    const FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> &
    getEigenVectors() const {
      return tEigenVectors;
    }
 
    int maxIterations = 100;
    int maxBracketIterations = 64;
    double absoluteTolerance = 1e-12;
    double relativeTolerance = 1e-12;
 
  private:
    static double logSumExp(const double first, const double second) {
      if (first == -std::numeric_limits<double>::infinity()) {
        return second;
      }
      if (second == -std::numeric_limits<double>::infinity()) {
        return first;
      }
      const double maximum = std::max(first, second);
      const double minimum = std::min(first, second);
      return maximum + std::log1p(std::exp(minimum - maximum));
    }
 
    std::pair<double, double> getSafeYBounds() const {
      const double safety_margin = std::log(16.);
      const double log_maximum =
          std::log(std::numeric_limits<double>::max()) - safety_margin;
 
      // gamma, gammaY, and z may safely underflow to zero because the scalar
      // residual is evaluated from their logarithms. The lower limit instead
      // keeps J = exp(-y/beta) representable.
      double lower_y = -materialParameters.beta * log_maximum;
      if (externalPressure < 0) {
        const double negative_infinity =
            -std::numeric_limits<double>::infinity();
        const double log_minus_q = std::log(-externalPressure);
        bool all_shifted_stresses_positive = true;
        double minimum_shifted_stress = std::numeric_limits<double>::max();
        for (int ii = 0; ii != SPACE_DIM; ++ii) {
          if (tB(ii) <= 0) {
            all_shifted_stresses_positive = false;
            break;
          }
          minimum_shifted_stress =
              std::min(minimum_shifted_stress, tB(ii));
        }
 
        double log_required_gamma = log_minus_q;
        if (all_shifted_stresses_positive) {
          const double log_spinodal_magnitude =
              2 * std::log(minimum_shifted_stress) -
              std::log(4.) - std::log(mu2);
          if (log_spinodal_magnitude > log_minus_q) {
            log_required_gamma = negative_infinity;
          } else if (log_spinodal_magnitude == log_minus_q) {
            log_required_gamma =
                log_minus_q +
                std::log(64 * std::numeric_limits<double>::epsilon());
          } else {
            log_required_gamma =
                log_minus_q +
                std::log(-std::expm1(log_spinodal_magnitude - log_minus_q));
          }
        }
 
        if (std::isfinite(log_required_gamma)) {
          const double log_mu1 =
              mu1 > 0 ? std::log(mu1) : negative_infinity;
          const double log_sqrt_discriminant =
              0.5 * logSumExp(2 * log_mu1,
                              std::log(mu2) + log_required_gamma);
          const double log_denominator =
              logSumExp(log_mu1, log_sqrt_discriminant);
          const double boundary_y =
              log_required_gamma - log_denominator;
          const double boundary_margin =
              64 * std::numeric_limits<double>::epsilon() *
              std::max(1., std::abs(boundary_y));
          lower_y = std::max(lower_y, boundary_y + boundary_margin);
        }
      }
 
      double upper_y = log_maximum;
      const int active_terms = mu1 > 0 ? 2 : 1;
      const double log_term_maximum =
          log_maximum - std::log(static_cast<double>(active_terms));
 
      // Bound gammaY, whose quadratic coefficient is the largest coefficient
      // used by gamma and gammaY.
      if (mu1 > 0) {
        upper_y = std::min(
            upper_y, log_term_maximum - std::log(2 * mu1));
      }
      upper_y = std::min(
          upper_y,
          0.5 * (log_term_maximum - std::log(2 * mu2)));
 
      // Also reserve enough range for 2*sqrt(mu2*gamma), which enters the
      // principal-stretch discriminant.
      const double log_gamma_effective_limit = std::min(
          log_maximum,
          2 * (log_maximum - std::log(2.)) - std::log(mu2));
      double log_gamma_limit = log_gamma_effective_limit;
      if (externalPressure > 0) {
        const double log_q = std::log(externalPressure);
        if (log_q >= log_gamma_effective_limit) {
          return {1., 0.};
        }
        log_gamma_limit =
            log_gamma_effective_limit +
            std::log(-std::expm1(log_q - log_gamma_effective_limit));
      }
      const double log_gamma_term_limit =
          log_gamma_limit - std::log(static_cast<double>(active_terms));
      if (mu1 > 0) {
        upper_y = std::min(
            upper_y, log_gamma_term_limit - std::log(2 * mu1));
      }
      upper_y = std::min(
          upper_y,
          0.5 * (log_gamma_term_limit - std::log(mu2)));
 
      return {lower_y, upper_y};
    }
 
    bool hasConverged(const ScalarEquationState &state, const double lower_y,
                      const double upper_y) const {
      const double y_scale = std::max(1., std::abs(state.y));
      const double y_tolerance =
          absoluteTolerance + relativeTolerance * y_scale;
      return std::abs(state.phi) <= absoluteTolerance ||
             std::abs(upper_y - lower_y) <= y_tolerance;
    }
 
    MaterialParameters materialParameters;
    double mu1;
    double mu2;
    double externalPressure = 0;
    FTensor::Tensor1<double, SPACE_DIM> tPrincipalStress;
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> tEigenVectors;
    FTensor::Tensor1<double, SPACE_DIM> tB;
    ScalarEquationState scalarState;
    int nbUniqueEigenvalues = 0;
    bool hasSpectralForm = false;
    bool hasScalarSolution = false;
  };
 
  /** @brief Assemble the isotropic external-strain contribution */
  struct OpSpatialPhysicalExternalStrain : public OpAssembleVolume {
    OpSpatialPhysicalExternalStrain(
        const std::string &field_name,
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
        std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
    MoFEMErrorCode integrate(EntData &data) override;
 
  private:
    boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
    std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
  };
 
  VolUserDataOperator *returnOpSpatialPhysicalExternalStrain(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) override {
    return new OpSpatialPhysicalExternalStrain(
        field_name, data_ptr, external_strain_vec_ptr, std::move(smv));
  }
 
  VolUserDataOperator *returnOpCalculateExternalPressure(
      VectorPtr external_pressure_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) override {
    return new OpCalculateExternalPressure(
        std::move(external_pressure_ptr),
        std::move(external_strain_vec_ptr), std::move(smv));
  }
 
  /**
   * @brief Evaluate Storakers stretch and its tangent from the supplied stress
   */
  struct OpCalculateStretchFromStress : public VolUserDataOperator {
    OpCalculateStretchFromStress(
        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
        boost::shared_ptr<MatrixDouble> strain_ptr,
        boost::shared_ptr<MatrixDouble> stress_ptr,
        boost::shared_ptr<HMHStorakers> storakers_ptr,
        VectorPtr external_pressure_ptr);
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override;
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<MatrixDouble> strainPtr;
    boost::shared_ptr<MatrixDouble> stressPtr;
    boost::shared_ptr<HMHStorakers> storakersPtr;
    VectorPtr externalPressurePtr;
  };
 
  VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr) override {
    return returnOpCalculateStretchFromStress(
        std::move(data_ptr), std::move(physics_ptr), std::move(strain_ptr),
        nullptr);
  }
 
  VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr,
      boost::shared_ptr<MatrixDouble> strain_ptr,
      VectorPtr external_pressure_ptr) override {
    auto storakers_ptr = boost::dynamic_pointer_cast<HMHStorakers>(physics_ptr);
    if (!storakers_ptr) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Pointer to HMHStorakers is null");
    }
 
    return new OpCalculateStretchFromStress(
        data_ptr,
        strain_ptr ? strain_ptr : data_ptr->getLogStretchTensorAtPts(),
        data_ptr->getApproxPAtPts(), storakers_ptr,
        std::move(external_pressure_ptr));
  }
 
private:
  static MoFEMErrorCode validateParameters(const MaterialParameters &parameters,
                                           const char *source) {
    MoFEMFunctionBegin;
 
    if (!std::isfinite(parameters.eta) || parameters.eta <= 0 ||
        parameters.eta > 1) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Storakers eta must satisfy 0 < eta <= 1 in %s", source);
    }
    if (!std::isfinite(parameters.mu) || parameters.mu <= 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Storakers mu must be positive in %s", source);
    }
    if (!std::isfinite(parameters.beta) || parameters.beta <= 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Storakers beta must be positive in %s", source);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  struct BlockData {
    double eta;
    double mu;
    double beta;
    Range blockEnts;
  };
 
  MoFEM::Interface &mField;
  MaterialParameters defaultParameters;
  double alphaGradU = 0;
  std::vector<BlockData> blockData;
};
 
HMHStorakers::OpSpatialPhysical::OpSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
    : OpAssembleVolume(field_name, std::move(data_ptr), OPROW),
      alphaU(alpha_u) {}
 
MoFEMErrorCode HMHStorakers::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto storakers_ptr =
      boost::dynamic_pointer_cast<HMHStorakers>(dataAtPts->physicsPtr);
  if (!storakers_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHStorakers is null");
  }
 
  const auto parameters =
      storakers_ptr->getMaterialParameters(getFEEntityHandle());
  const double alpha_u = alphaU;
  const double alpha_grad_u = storakers_ptr->alphaGradU;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  FTENSOR_INDEX(size_symm, L);
 
  const auto t_L = FTensor::SymmLTensor<double, SPACE_DIM>();
  const int nb_dofs = data.getIndices().size();
  const int nb_integration_pts = data.getN().size1();
  const int nb_base_functions = data.getN().size2();
  const double volume = getVolume();
 
  auto t_w = getFTensor0IntegrationWeight();
  auto t_adjoint = dataAtPts->getFTensorAdjointPdU(nb_integration_pts);
  auto t_dot_log_u =
      dataAtPts->getFTensorLogStretchDot(nb_integration_pts);
  auto t_grad_dot_log_u =
      dataAtPts->getFTensorGradLogStretchDot(nb_integration_pts);
  auto t_eigen_vals = dataAtPts->getFTensorEigenVals(nb_integration_pts);
  auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(nb_integration_pts);
 
  auto t_row_base_fun = data.getFTensor0N();
  auto t_row_grad_fun = data.getFTensor1DiffN<SPACE_DIM>();
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double a = volume * t_w;
      ++t_w;
 
      const double jacobian_power =
          getJacobianPower(parameters, t_eigen_vals);
      const double gamma = getGamma(parameters, jacobian_power);
      if (!std::isfinite(jacobian_power) || jacobian_power < 0 ||
          !std::isfinite(gamma)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
                "Non-finite Storakers volumetric response in spatial "
                "residual");
      }
 
      auto coordinate_stress = [parameters, gamma](const double h) {
        return getPrincipalCoordinateStress(parameters, gamma, h);
      };
      auto t_material_stress = EigenMatrix::getMat(
          t_eigen_vals, t_eigen_vecs, coordinate_stress);
 
      FTensor::Tensor1<double, size_symm> t_material;
      t_material(L) = t_L(i, j, L) * t_material_stress(i, j);
 
      ++t_eigen_vals;
      ++t_eigen_vecs;
 
      FTensor::Tensor1<double, size_symm> t_viscous;
      t_viscous(L) =
          alpha_u * (t_L(i, j, L) * t_dot_log_u(i, j));
 
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = t_adjoint(L) - t_material(L) - t_viscous(L);
      t_residual(L) *= a;
 
      FTensor::Tensor2<double, size_symm, SPACE_DIM> t_grad_residual;
      t_grad_residual(L, i) = alpha_grad_u * t_grad_dot_log_u(L, i);
      t_grad_residual(L, i) *= a;
 
      ++t_adjoint;
      ++t_dot_log_u;
      ++t_grad_dot_log_u;
 
      auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
      int bb = 0;
      for (; bb != nb_dofs / size_symm; ++bb) {
        t_nf(L) -= t_row_base_fun * t_residual(L);
        t_nf(L) += t_row_grad_fun(i) * t_grad_residual(L, i);
        ++t_nf;
        ++t_row_base_fun;
        ++t_row_grad_fun;
      }
      for (; bb != nb_base_functions; ++bb) {
        ++t_row_base_fun;
        ++t_row_grad_fun;
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto not_implemented = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Storakers full-order operator is implemented only for no_h1");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
    CHKERR not_implemented();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
  }
 
  MoFEMFunctionReturn(0);
}
 
HMHStorakers::OpSpatialPhysical_du_du::OpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
    : OpAssembleVolumePositiveDefine(std::move(row_field),
                                     std::move(col_field),
                                     std::move(data_ptr), OPROWCOL, false),
      alphaU(alpha_u) {
  CHK_THROW_MESSAGE(getOptions(), "get options failed");
  sYmm = false;
}
 
MoFEMErrorCode HMHStorakers::OpSpatialPhysical_du_du::getOptions() {
  MoFEMFunctionBegin;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "storakers_", "Storakers material",
                    "none");
  CHKERR PetscOptionsScalar("-min_eigen_value", "Minimum eigenvalue", "",
                            minimEigenValue, &minimEigenValue, PETSC_NULLPTR);
  PetscOptionsEnd();
  if (!std::isfinite(minimEigenValue) || minimEigenValue < 0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Storakers min_eigen_value must be finite and non-negative");
  }
  MOFEM_LOG("EP", Sev::inform)
      << "Storakers min_eigen_value = " << minimEigenValue;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode HMHStorakers::OpSpatialPhysical_du_du::integrate(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto storakers_ptr =
      boost::dynamic_pointer_cast<HMHStorakers>(dataAtPts->physicsPtr);
  if (!storakers_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHStorakers is null");
  }
 
  const auto parameters =
      storakers_ptr->getMaterialParameters(getFEEntityHandle());
  const double alpha_u = alphaU;
  const double alpha_grad_u = storakers_ptr->alphaGradU;
  const double bulk_modulus =
      getDerivedParameters(parameters.mu, parameters.beta).bulkModulus;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
  FTENSOR_INDEXES(size_symm, L, J);
 
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
  const auto t_L = FTensor::SymmLTensor<double, SPACE_DIM>();
  const auto t_diff = FTensor::DiffTensor<double>();
 
  const int nb_integration_pts = row_data.getN().size1();
  const int row_nb_dofs = row_data.getIndices().size();
  const int col_nb_dofs = col_data.getIndices().size();
  const int row_nb_base_functions = row_data.getN().size2();
 
  auto get_ftensor2 = [](MatrixDouble &matrix, const int row,
                         const int col) {
    std::array<double *, size_symm * size_symm> pointers;
    for (int rr = 0; rr != size_symm; ++rr) {
      for (int cc = 0; cc != size_symm; ++cc) {
        pointers[size_symm * rr + cc] = &matrix(row + rr, col + cc);
      }
    }
    return FTensor::Tensor2<FTensor::PackPtr<double *, size_symm>, size_symm,
                            size_symm>(pointers);
  };
 
  const double volume = getVolume();
  const double ts_a = getTSa();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_row_grad_fun = row_data.getFTensor1DiffN<SPACE_DIM>();
 
  auto t_adjoint_dstretch =
      dataAtPts->getFTensorAdjointPdstretch(nb_integration_pts);
  auto t_eigen_vals = dataAtPts->getFTensorEigenVals(nb_integration_pts);
  auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(nb_integration_pts);
  auto t_nb_uniq = FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(
      dataAtPts->nbUniq.data().data());
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      const double a = volume * t_w;
      ++t_w;
 
      const double jacobian_power =
          getJacobianPower(parameters, t_eigen_vals);
      const double gamma = getGamma(parameters, jacobian_power);
      const double gamma_y = getGammaY(parameters, jacobian_power);
      const double volumetric_tangent = parameters.beta * gamma_y;
      if (!std::isfinite(jacobian_power) || jacobian_power < 0 ||
          !std::isfinite(gamma) || !std::isfinite(gamma_y) ||
          !std::isfinite(volumetric_tangent)) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
                "Non-finite Storakers volumetric response in spatial "
                "tangent");
      }
 
      auto coordinate_stress = [parameters, gamma](const double h) {
        return getPrincipalCoordinateStress(parameters, gamma, h);
      };
      auto coordinate_stress_derivative =
          [parameters](const double h) {
            return getPrincipalCoordinateStressDerivativeAtFixedGamma(
                parameters, h);
          };
      auto t_frozen_gamma_tangent = EigenMatrix::getDiffMat(
          t_eigen_vals, t_eigen_vecs, coordinate_stress,
          coordinate_stress_derivative, t_nb_uniq);
      auto identity = [](const double) { return 1.; };
      auto t_log_jacobian_gradient =
          EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, identity);
 
      FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_material_tangent;
      t_material_tangent(i, j, k, l) =
          t_frozen_gamma_tangent(i, j, k, l) +
          volumetric_tangent * t_log_jacobian_gradient(i, j) *
              t_log_jacobian_gradient(k, l);
 
      FTensor::Tensor2<double, size_symm, size_symm> t_dP;
      t_dP(L, J) =
          t_L(i, j, L) *
          (t_material_tangent(i, j, k, l) * t_L(k, l, J));
      t_dP(L, J) +=
          (alpha_u * ts_a) *
          (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_deltaP;
      t_deltaP(i, j) =
          (t_adjoint_dstretch(i, j) || t_adjoint_dstretch(j, i)) / 2.;
      auto t_diff2_uP = EigenMatrix::getDiffDiffMat(
          t_eigen_vals, t_eigen_vecs,
          static_cast<double (*)(double)>(std::exp),
          static_cast<double (*)(double)>(std::exp),
          static_cast<double (*)(double)>(std::exp), t_deltaP, t_nb_uniq);
      t_dP(L, J) -=
          t_L(i, j, L) * (t_diff2_uP(i, j, k, l) * t_L(k, l, J));
 
      ++t_adjoint_dstretch;
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_nb_uniq;
 
      t_dP(L, J) *= a;
      if (minimEigenValue > 0) {
        FTensor::Tensor2<double, size_symm, size_symm> t_hessian_eig_vecs;
        t_hessian_eig_vecs(L, J) =
            0.5 * (t_dP(L, J) + t_dP(J, L));
        FTensor::Tensor1<double, size_symm> t_hessian_eig_vals;
        CHKERR computeEigenValuesSymmetric(t_hessian_eig_vecs,
                                           t_hessian_eig_vals);
 
        const double eigenvalue_floor = bulk_modulus * minimEigenValue;
        bool project_hessian = false;
        for (int aa = 0; aa != size_symm; ++aa) {
          project_hessian =
              project_hessian || t_hessian_eig_vals(aa) < eigenvalue_floor;
        }
        if (project_hessian) {
          auto min_eig_val = [eigenvalue_floor](const double value) {
            return 0.5 *
                   (value + eigenvalue_floor +
                    std::abs(value - eigenvalue_floor));
          };
          auto t_dP_min_eig = EigenMatrix::getMat(
              t_hessian_eig_vals, t_hessian_eig_vecs, min_eig_val);
          t_dP(L, J) = t_dP_min_eig(L, J);
        }
      }
 
      int rr = 0;
      for (; rr != row_nb_dofs / size_symm; ++rr) {
        auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
        auto t_col_grad_fun =
            col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
        auto t_m = get_ftensor2(K, size_symm * rr, 0);
 
        for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
          t_m(L, J) +=
              (t_row_base_fun * t_col_base_fun) * t_dP(L, J);
          t_m(L, J) +=
              (a * alpha_grad_u * ts_a) *
              (t_row_grad_fun(i) * t_col_grad_fun(i)) * t_kd_sym(L, J);
          ++t_m;
          ++t_col_base_fun;
          ++t_col_grad_fun;
        }
        ++t_row_base_fun;
        ++t_row_grad_fun;
      }
      for (; rr != row_nb_base_functions; ++rr) {
        ++t_row_base_fun;
        ++t_row_grad_fun;
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto not_implemented = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Storakers full-order tangent is implemented only for no_h1");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
    CHKERR not_implemented();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
  }
 
  MoFEMFunctionReturn(0);
}
 
HMHStorakers::OpSpatialPhysicalExternalStrain::OpSpatialPhysicalExternalStrain(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
    std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
    : OpAssembleVolume(field_name, data_ptr, OPROW),
      externalStrainVecPtr(std::move(external_strain_vec_ptr)),
      scalingMethodsMap(std::move(smv)) {}
 
MoFEMErrorCode
HMHStorakers::OpSpatialPhysicalExternalStrain::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto storakers_ptr =
      boost::dynamic_pointer_cast<HMHStorakers>(dataAtPts->physicsPtr);
  if (!storakers_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHStorakers is null");
  }
  if (!externalStrainVecPtr || externalStrainVecPtr->empty()) {
    MoFEMFunctionReturnHot(0);
  }
 
  double time = OpAssembleVolume::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
  const EntityHandle fe_ent = OpAssembleVolume::getFEEntityHandle();
 
  for (const auto &external_strain_block : *externalStrainVecPtr) {
    if (std::regex_match(external_strain_block.blockName,
                         std::regex("(.*)ANALYTICAL_EXTERNALSTRAIN(.*)"))) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Analytical external strain not implemented for Storakers "
              "material");
    }
    if (external_strain_block.ents.find(fe_ent) ==
        external_strain_block.ents.end()) {
      continue;
    }
 
    double scale = 1;
    const auto scaling_method =
        scalingMethodsMap.find(external_strain_block.blockName);
    if (scaling_method != scalingMethodsMap.end()) {
      if (!scaling_method->second) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Null scaling method for %s",
                external_strain_block.blockName.c_str());
      }
      scale *= scaling_method->second->getScale(time);
    } else {
      MOFEM_LOG("EP", Sev::warning)
          << "No scaling method found for "
          << external_strain_block.blockName;
    }
 
    const double q = 3 * external_strain_block.bulkModulusK * scale *
                     external_strain_block.val;
    if (!std::isfinite(q)) {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Non-finite Storakers external pressure for %s",
              external_strain_block.blockName.c_str());
    }
 
    FTENSOR_INDEXES(SPACE_DIM, i, j);
    FTENSOR_INDEX(size_symm, L);
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
    const auto t_L = FTensor::SymmLTensor<double, SPACE_DIM>();
 
    const int nb_dofs = data.getIndices().size();
    const int nb_integration_pts = data.getN().size1();
    const int nb_base_functions = data.getN().size2();
    const double volume = getVolume();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base_fun = data.getFTensor0N();
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, size_symm> t_residual;
      t_residual(L) = volume * t_w * q *
                      (t_L(i, j, L) * t_kd(i, j));
      ++t_w;
 
      auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
      int bb = 0;
      for (; bb != nb_dofs / size_symm; ++bb) {
        t_nf(L) -= t_row_base_fun * t_residual(L);
        ++t_nf;
        ++t_row_base_fun;
      }
      for (; bb != nb_base_functions; ++bb) {
        ++t_row_base_fun;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
HMHStorakers::OpCalculateStretchFromStress::OpCalculateStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<HMHStorakers> storakers_ptr,
    VectorPtr external_pressure_ptr)
    : VolUserDataOperator(H1, OPLAST), dataAtPts(std::move(data_ptr)),
      strainPtr(std::move(strain_ptr)), stressPtr(std::move(stress_ptr)),
      storakersPtr(std::move(storakers_ptr)),
      externalPressurePtr(std::move(external_pressure_ptr)) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBVERTEX] = true;
}
 
MoFEMErrorCode
HMHStorakers::OpCalculateStretchFromStress::doWork(int, EntityType, EntData &) {
  MoFEMFunctionBegin;
 
  if (!dataAtPts || !strainPtr || !stressPtr || !storakersPtr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Invalid data passed to Storakers stretch-from-stress operator");
  }
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
 
  const int nb_integration_pts = stressPtr->size1();
#ifndef NDEBUG
  if (nb_integration_pts != getGaussPts().size2()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Inconsistent number of Storakers integration points");
  }
#endif
  if (externalPressurePtr &&
      externalPressurePtr->size() != nb_integration_pts) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Inconsistent number of external-pressure integration points");
  }
 
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                   DL>::size(*strainPtr, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                   DL>::size(*dataAtPts->getStretchTensorAtPts(),
                             nb_integration_pts);
  MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(*dataAtPts->getDiffStretchTensorAtPts(),
                             nb_integration_pts);
  MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
      *dataAtPts->getEigenVals(), nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(*dataAtPts->getEigenVecs(), nb_integration_pts);
  dataAtPts->nbUniq.resize(nb_integration_pts, false);
  MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(dataAtPts->matD, nb_integration_pts);
  MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(dataAtPts->matInvD, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                   DL>::size(dataAtPts->logStretch2H1AtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2SymmetricFromMatType<SPACE_DIM, -1, DL>,
                   DL>::size(dataAtPts->logStretchTotalTensorAtPts,
                             nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(dataAtPts->rotMatAtPts, nb_integration_pts);
  MatrixSizeHelper<GetFTensor2FromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                   DL>::size(*dataAtPts->getAdjointPdstretchAtPts(),
                             nb_integration_pts);
 
  strainPtr->clear();
  dataAtPts->getStretchTensorAtPts()->clear();
  dataAtPts->getDiffStretchTensorAtPts()->clear();
  dataAtPts->getEigenVals()->clear();
  dataAtPts->getEigenVecs()->clear();
  dataAtPts->nbUniq.clear();
  dataAtPts->matD.clear();
  dataAtPts->matInvD.clear();
  dataAtPts->logStretch2H1AtPts.clear();
  dataAtPts->logStretchTotalTensorAtPts.clear();
 
  auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM, -1, DL>(*strainPtr);
  auto t_biot_stretch = dataAtPts->getFTensorStretch(nb_integration_pts);
  auto t_diff_stretch = dataAtPts->getFTensorDiffStretch(nb_integration_pts);
  auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM, -1, DL>(*stressPtr);
  auto t_omega = dataAtPts->getFTensorRotAxis(nb_integration_pts);
  auto t_R = dataAtPts->getFTensorRotMat(nb_integration_pts);
  auto t_biot_stress =
      dataAtPts->getFTensorAdjointPdstretch(nb_integration_pts);
  auto t_eigen_vals = dataAtPts->getFTensorEigenVals(nb_integration_pts);
  auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(nb_integration_pts);
  auto t_mat_d =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, -1, DL>(dataAtPts->matD);
  auto t_mat_inv_d =
      getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, -1, DL>(dataAtPts->matInvD);
  auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(nb_integration_pts);
  auto t_log_stretch_total =
      dataAtPts->getFTensorLogStretchTotal(nb_integration_pts);
  auto t_nb_uniq = FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(
      dataAtPts->nbUniq.data().data());
 
  const auto parameters =
      storakersPtr->getMaterialParameters(getFEEntityHandle());
  StressToStrech stress_to_stretch(parameters);
  constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
  constexpr auto identity = [](const double value) { return value; };
 
  if (EshelbianCore::rotSelector != LARGE_ROT) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Rotation selector not handled by Storakers");
  }
 
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
 
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_stress;
    t_rotated_stress(i, j) = t_R(k, i) * t_stress(k, j);
    t_biot_stress(i, j) = t_diff_sym(i, j, k, l) * t_rotated_stress(k, l);
 
    const double initial_y = stress_to_stretch.getScalarEquationState().y;
    const double q = externalPressurePtr ? (*externalPressurePtr)[gg] : 0.;
    const auto calculate_error = stress_to_stretch.calculate(
        t_biot_stress, t_biot_stretch, t_mat_d, t_mat_inv_d, q);
    if (calculate_error) {
      // Keep the production diagnostic to one copy-paste line. The atom prints
      // the material and spectral state, plus the scalar state on convergence.
      std::ostringstream atom_reproducer;
      atom_reproducer
          << std::scientific
          << std::setprecision(std::numeric_limits<double>::max_digits10)
          << "./storakers_stress_to_stretch_atom"
          << " -storakers_atom_eta " << parameters.eta
          << " -storakers_atom_mu " << parameters.mu
          << " -storakers_atom_beta " << parameters.beta
          << " -storakers_atom_q " << q
          << " -storakers_atom_initial_y " << initial_y
          << " -storakers_atom_solver_tol "
          << stress_to_stretch.absoluteTolerance
          << " -storakers_atom_biot_stress " << t_biot_stress(0, 0) << ","
          << t_biot_stress(0, 1) << "," << t_biot_stress(0, 2) << ","
          << t_biot_stress(1, 1) << "," << t_biot_stress(1, 2) << ","
          << t_biot_stress(2, 2) << " -log_no_color";
 
      PetscMPIInt rank = -1;
      (void)MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
      std::ostringstream atom_reproducer_log;
      atom_reproducer_log << "[" << rank
                          << "] <warning> atom_reproducer_command="
                          << atom_reproducer.str() << '\n';
      const auto atom_reproducer_log_string = atom_reproducer_log.str();
      std::fwrite(atom_reproducer_log_string.data(), 1,
                  atom_reproducer_log_string.size(), stderr);
      std::fflush(stderr);
    }
    CHKERR calculate_error;
 
    const auto &state = stress_to_stretch.getScalarEquationState();
    const auto &eigenvectors = stress_to_stretch.getEigenVectors();
    t_eigen_vals(i) = state.tLogU(i);
    t_eigen_vecs(i, j) = eigenvectors(i, j);
    for (int aa = 0; aa != SPACE_DIM; ++aa) {
      if (!std::isfinite(t_eigen_vals(aa))) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "Non-finite Storakers strain eigenvalue");
      }
    }
 
    auto t_log_stretch =
        EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, identity);
    t_strain(i, j) = t_log_stretch(i, j);
    t_log_stretch_total(i, j) = t_strain(i, j);
    t_log_u2_h1(i, j) = 0;
 
    t_nb_uniq =
        getUniqNb<SPACE_DIM>(getVectorAdaptor(&t_eigen_vals(0), SPACE_DIM),
                             FTensor::Number<SPACE_DIM>{});
    if (t_nb_uniq < SPACE_DIM) {
      CHKERR sortEigenVals<SPACE_DIM>(
          getVectorAdaptor(&t_eigen_vals(0), SPACE_DIM),
          getMatrixAdaptor(&t_eigen_vecs(0, 0), SPACE_DIM, SPACE_DIM),
          FTensor::Number<SPACE_DIM>{});
    }
    auto t_diff_stretch_mat = EigenMatrix::getDiffMat(
        t_eigen_vals, t_eigen_vecs,
        static_cast<double (*)(double)>(std::exp),
        static_cast<double (*)(double)>(std::exp), t_nb_uniq);
    t_diff_stretch(i, j, k, l) = t_diff_stretch_mat(i, j, k, l);
 
    ++t_strain;
    ++t_biot_stretch;
    ++t_diff_stretch;
    ++t_stress;
    ++t_omega;
    ++t_R;
    ++t_biot_stress;
    ++t_eigen_vals;
    ++t_eigen_vecs;
    ++t_mat_d;
    ++t_mat_inv_d;
    ++t_log_u2_h1;
    ++t_log_stretch_total;
    ++t_nb_uniq;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testHMHStorakersStressToStretch() {
  MoFEMFunctionBegin;
 
  using MaterialParameters = HMHStorakers::MaterialParameters;
  using StressToStrech = HMHStorakers::StressToStrech;
 
  PetscReal finite_difference_step = 1e-6;
  PetscReal stretch_tolerance = 1e-10;
  PetscReal compliance_tolerance = 1e-5;
  PetscReal solver_tolerance = 1e-14;
  PetscReal atom_eta = 0.37;
  PetscReal atom_mu = 2.4;
  PetscReal atom_beta = 0.8;
  PetscReal atom_q = 0;
  PetscReal atom_initial_y = 0;
  constexpr PetscInt symmetric_tensor_size = SPACE_DIM * (SPACE_DIM + 1) / 2;
  std::array<PetscReal, symmetric_tensor_size + 1> atom_biot_stress_input{};
  std::array<PetscReal, symmetric_tensor_size> atom_biot_stress{};
  PetscInt atom_biot_stress_size = atom_biot_stress_input.size();
  PetscBool atom_biot_stress_set = PETSC_FALSE;
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_fd_step", &finite_difference_step,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_stretch_tol", &stretch_tolerance,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_compliance_tol",
                             &compliance_tolerance, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_solver_tol", &solver_tolerance,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_eta", &atom_eta, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_mu", &atom_mu, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_beta", &atom_beta, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_q", &atom_q, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-storakers_atom_initial_y", &atom_initial_y,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetRealArray(
      PETSC_NULLPTR, PETSC_NULLPTR, "-storakers_atom_biot_stress",
      atom_biot_stress_input.data(), &atom_biot_stress_size,
      &atom_biot_stress_set);
  if (!std::isfinite(finite_difference_step) || finite_difference_step <= 0 ||
      !std::isfinite(stretch_tolerance) || stretch_tolerance <= 0 ||
      !std::isfinite(compliance_tolerance) || compliance_tolerance <= 0 ||
      !std::isfinite(solver_tolerance) || solver_tolerance <= 0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Invalid Storakers atom-test tolerance");
  }
  if (atom_biot_stress_set &&
      atom_biot_stress_size != symmetric_tensor_size) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "-storakers_atom_biot_stress needs exactly six values in order "
            "s00,s01,s02,s11,s12,s22");
  }
  if (atom_biot_stress_set) {
    if (!std::isfinite(atom_eta) || !std::isfinite(atom_mu) ||
        !std::isfinite(atom_beta) || !std::isfinite(atom_q) ||
        !std::isfinite(atom_initial_y)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Storakers reproduction parameters must be finite");
    }
    for (int rr = 0; rr != symmetric_tensor_size; ++rr) {
      if (!std::isfinite(atom_biot_stress_input[rr])) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "-storakers_atom_biot_stress values must be finite");
      }
      atom_biot_stress[rr] = atom_biot_stress_input[rr];
    }
  }
 
  struct StretchCase {
    const char *name;
    MaterialParameters parameters;
    std::array<double, SPACE_DIM> principalStretch;
    bool rotate;
    double q;
  };
 
  const MaterialParameters parameters{0.37, 2.4, 0.8};
  const MaterialParameters near_incompressible_parameters{0.5, 0.0080194,
                                                           248.};
  // eta=1 is the power-law Neo-Hookean specialization used by
  // HMHNeohookean: mu=2*c10 and beta=p/2.
  const MaterialParameters power_neohookean_parameters{1., 3.4, 11. / 12.};
  const std::array<StretchCase, 7> stretch_cases{{
      {"reference", parameters, {1., 1., 1.}, false, 0.},
      {"rotated_distinct_q", parameters, {0.71, 1.08, 1.43}, true, 0.65},
      {"rotated_repeated_negative_q", parameters, {0.82, 0.82, 1.36}, true,
       -0.4},
      {"negative_gamma_eff", parameters, {1., 1., 1.}, false, -4.},
      {"large_beta_dilated", near_incompressible_parameters, {4., 4., 4.},
       false, 0.},
      {"power_neohookean_distinct", power_neohookean_parameters,
       {0.62, 1.12, 1.51}, true, 0.},
      {"power_neohookean_repeated_q", power_neohookean_parameters,
       {0.55, 0.55, 0.9}, true, 0.75},
  }};
  const std::array<std::array<int, 2>, 6> mandel_directions{{
      {0, 0},
      {1, 1},
      {2, 2},
      {1, 2},
      {0, 2},
      {0, 1},
  }};
 
  auto set_solver_tolerances = [solver_tolerance](auto &solver) {
    solver.absoluteTolerance = solver_tolerance;
    solver.relativeTolerance = solver_tolerance;
  };
 
  if (atom_biot_stress_set) {
    const MaterialParameters atom_parameters{atom_eta, atom_mu, atom_beta};
    StressToStrech response(atom_parameters);
    set_solver_tolerances(response);
 
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress;
    t_stress(0, 0) = atom_biot_stress[0];
    t_stress(0, 1) = atom_biot_stress[1];
    t_stress(0, 2) = atom_biot_stress[2];
    t_stress(1, 1) = atom_biot_stress[3];
    t_stress(1, 2) = atom_biot_stress[4];
    t_stress(2, 2) = atom_biot_stress[5];
 
    MOFEM_LOG("EP", Sev::inform)
        << std::scientific
        << std::setprecision(std::numeric_limits<double>::max_digits10)
        << "Storakers reproduction input: eta=" << atom_eta
        << " mu=" << atom_mu << " beta=" << atom_beta
        << " q=" << atom_q
        << " initial_y=" << atom_initial_y
        << " solver_tol=" << solver_tolerance
        << " biot_stress=[" << atom_biot_stress[0] << ","
        << atom_biot_stress[1] << "," << atom_biot_stress[2] << ","
        << atom_biot_stress[3] << "," << atom_biot_stress[4] << ","
        << atom_biot_stress[5] << "]";
 
    CHKERR response.stressToSpectralForm(t_stress, atom_q);
    const auto &principal_stress = response.getPrincipalStress();
    const auto &shifted_principal_stress =
        response.getShiftedPrincipalStress();
    MOFEM_LOG("EP", Sev::inform)
        << std::scientific
        << std::setprecision(std::numeric_limits<double>::max_digits10)
        << "Storakers reproduction spectrum: principal_stress=["
        << principal_stress(0) << "," << principal_stress(1) << ","
        << principal_stress(2) << "] shifted_principal_stress=["
        << shifted_principal_stress(0) << ","
        << shifted_principal_stress(1) << ","
        << shifted_principal_stress(2) << "]";
 
    CHKERR response.solveScalarEquation(atom_initial_y);
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_recovered_stretch;
    FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_stiffness;
    FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_compliance;
    CHKERR response.spectralFormToStrech(t_recovered_stretch);
    CHKERR response.calculateBiotTangent(t_stiffness, t_compliance);
 
    const auto &final_state = response.getScalarEquationState();
    MOFEM_LOG("EP", Sev::inform)
        << std::scientific
        << std::setprecision(std::numeric_limits<double>::max_digits10)
        << "Storakers reproduction converged scalar state: y=" << final_state.y
        << " z=" << final_state.z << " gamma=" << final_state.gamma
        << " q=" << final_state.q
        << " gamma_eff=" << final_state.gammaEffective
        << " gamma_y=" << final_state.gammaY << " phi=" << final_state.phi
        << " d_phi=" << final_state.dPhi << " u=[" << final_state.tU(0)
        << "," << final_state.tU(1) << "," << final_state.tU(2) << "]"
        << " jacobian=" << response.getJacobian();
    MOFEM_LOG("EP", Sev::inform)
        << "Storakers reproduction recovered stretch:\n"
        << t_recovered_stretch;
    MoFEMFunctionReturnHot(0);
  }
 
  for (const auto &stretch_case : stretch_cases) {
    const auto &p = stretch_case.parameters;
    const double mu1 = (1 - p.eta) * p.mu;
    const double mu2 = p.eta * p.mu;
 
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigenvectors;
    for (int ii = 0; ii != SPACE_DIM; ++ii) {
      for (int jj = 0; jj != SPACE_DIM; ++jj) {
        t_eigenvectors(ii, jj) = ii == jj ? 1. : 0.;
      }
    }
    if (stretch_case.rotate) {
      // Rows are orthonormal eigenvectors in the MoFEM convention.
      t_eigenvectors(0, 0) = 2. / 3.;
      t_eigenvectors(0, 1) = -2. / 3.;
      t_eigenvectors(0, 2) = 1. / 3.;
      t_eigenvectors(1, 0) = 2. / 3.;
      t_eigenvectors(1, 1) = 1. / 3.;
      t_eigenvectors(1, 2) = -2. / 3.;
      t_eigenvectors(2, 0) = 1. / 3.;
      t_eigenvectors(2, 1) = 2. / 3.;
      t_eigenvectors(2, 2) = 2. / 3.;
    }
 
    const double jacobian = stretch_case.principalStretch[0] *
                            stretch_case.principalStretch[1] *
                            stretch_case.principalStretch[2];
    const double z = std::pow(jacobian, -p.beta);
    const double gamma = 2 * mu1 * z + mu2 * z * z;
 
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_expected_stretch;
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_inverse_stretch;
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress;
    for (int ii = 0; ii != SPACE_DIM; ++ii) {
      for (int jj = ii; jj != SPACE_DIM; ++jj) {
        t_expected_stretch(ii, jj) = 0;
        t_inverse_stretch(ii, jj) = 0;
        for (int aa = 0; aa != SPACE_DIM; ++aa) {
          const double projector =
              t_eigenvectors(aa, ii) * t_eigenvectors(aa, jj);
          t_expected_stretch(ii, jj) +=
              projector * stretch_case.principalStretch[aa];
          t_inverse_stretch(ii, jj) +=
              projector / stretch_case.principalStretch[aa];
        }
        t_stress(ii, jj) = (ii == jj ? 2 * mu1 : 0.) +
                           mu2 * t_expected_stretch(ii, jj) -
                           (gamma + stretch_case.q) *
                               t_inverse_stretch(ii, jj);
      }
    }
 
    StressToStrech base_response(p);
    set_solver_tolerances(base_response);
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_recovered_stretch;
    FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_stiffness;
    FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_compliance;
    CHKERR base_response.calculate(t_stress, t_recovered_stretch, t_stiffness,
                                   t_compliance, stretch_case.q);
 
    double maximum_stretch_error = 0;
    double stretch_scale = 1;
    double stress_norm_squared = 0;
    for (int ii = 0; ii != SPACE_DIM; ++ii) {
      for (int jj = 0; jj != SPACE_DIM; ++jj) {
        if (!std::isfinite(t_recovered_stretch(ii, jj)) ||
            !std::isfinite(t_expected_stretch(ii, jj)) ||
            !std::isfinite(t_stress(ii, jj))) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Non-finite Storakers stretch-recovery result");
        }
        maximum_stretch_error = std::max(
            maximum_stretch_error,
            std::abs(t_recovered_stretch(ii, jj) - t_expected_stretch(ii, jj)));
        stretch_scale =
            std::max(stretch_scale, std::abs(t_expected_stretch(ii, jj)));
        stress_norm_squared += t_stress(ii, jj) * t_stress(ii, jj);
      }
    }
    const double scaled_stretch_error = maximum_stretch_error / stretch_scale;
    MOFEM_LOG("EP", Sev::inform)
        << "Storakers stretch recovery: case=" << stretch_case.name
        << " q=" << stretch_case.q
        << " scaled_error=" << scaled_stretch_error;
    if (!std::isfinite(scaled_stretch_error) ||
        scaled_stretch_error > stretch_tolerance) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Storakers stress-to-stretch recovery failed");
    }
 
    const double perturbation =
        finite_difference_step *
        std::max({1., p.mu, std::sqrt(stress_norm_squared)});
    double maximum_compliance_error = 0;
    int worst_direction = -1;
 
    for (int column = 0; column != static_cast<int>(mandel_directions.size());
         ++column) {
      const int kk = mandel_directions[column][0];
      const int ll = mandel_directions[column][1];
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction;
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress_plus;
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress_minus;
      for (int ii = 0; ii != SPACE_DIM; ++ii) {
        for (int jj = ii; jj != SPACE_DIM; ++jj) {
          t_direction(ii, jj) = 0;
        }
      }
      t_direction(kk, ll) = kk == ll ? 1. : 1. / std::sqrt(2.);
      for (int ii = 0; ii != SPACE_DIM; ++ii) {
        for (int jj = ii; jj != SPACE_DIM; ++jj) {
          t_stress_plus(ii, jj) =
              t_stress(ii, jj) + perturbation * t_direction(ii, jj);
          t_stress_minus(ii, jj) =
              t_stress(ii, jj) - perturbation * t_direction(ii, jj);
        }
      }
 
      StressToStrech plus_response(p);
      StressToStrech minus_response(p);
      set_solver_tolerances(plus_response);
      set_solver_tolerances(minus_response);
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stretch_plus;
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stretch_minus;
      CHKERR plus_response.calculate(t_stress_plus, t_stretch_plus,
                                     stretch_case.q);
      CHKERR minus_response.calculate(t_stress_minus, t_stretch_minus,
                                      stretch_case.q);
 
      double error_norm_squared = 0;
      double analytical_norm_squared = 0;
      double numerical_norm_squared = 0;
      for (int ii = 0; ii != SPACE_DIM; ++ii) {
        for (int jj = 0; jj != SPACE_DIM; ++jj) {
          double analytical = 0;
          for (int mm = 0; mm != SPACE_DIM; ++mm) {
            for (int nn = 0; nn != SPACE_DIM; ++nn) {
              analytical += t_compliance(ii, jj, mm, nn) * t_direction(mm, nn);
            }
          }
          const double numerical =
              (t_stretch_plus(ii, jj) - t_stretch_minus(ii, jj)) /
              (2 * perturbation);
          const double error = analytical - numerical;
          if (!std::isfinite(analytical) || !std::isfinite(numerical) ||
              !std::isfinite(error)) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Non-finite Storakers compliance result");
          }
          error_norm_squared += error * error;
          analytical_norm_squared += analytical * analytical;
          numerical_norm_squared += numerical * numerical;
        }
      }
      const double scaled_compliance_error =
          std::sqrt(error_norm_squared) /
          std::max({1., std::sqrt(analytical_norm_squared),
                    std::sqrt(numerical_norm_squared)});
      if (!std::isfinite(scaled_compliance_error)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Non-finite Storakers compliance error");
      }
      if (scaled_compliance_error > maximum_compliance_error) {
        maximum_compliance_error = scaled_compliance_error;
        worst_direction = column;
      }
    }
 
    MOFEM_LOG("EP", Sev::inform)
        << "Storakers compliance finite difference: case=" << stretch_case.name
        << " max_scaled_error=" << maximum_compliance_error
        << " direction=" << worst_direction;
    if (!std::isfinite(maximum_compliance_error) ||
        maximum_compliance_error > compliance_tolerance) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Storakers compliance finite-difference check failed");
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity