HMHNeohookean.cpp

Go to the documentation of this file.

/**
 * @file HMHNeohookean.cpp
 * @brief Abaqus-compatible compressible Neo-Hookean material in principal
 * stretches
 * @date 2024-08-31
 *
 * @copyright Copyright (c) 2024
 *
 * Constitutive synopsis
 * ---------------------
 * The deformation gradient is split as F = R U, where R is a rigid rotation
 * and U is the symmetric positive-definite stretch tensor.  If lambda_a are
 * the eigenvalues (principal stretches) of U, this implementation evaluates
 *
 *   W(U) = c10 (J^{-2/3} tr(U^2) - 3) + K/2 (J - 1)^2,
 *   J = det(U),  mu_0 = 2 c10,  D_1 = 2/K.
 *
 * The internal scalar variable is the natural Hencky stretch
 * h_a=log(lambda_a), so lambda_a=exp(h_a).
 *
 * This is the same reduced-polynomial N=1 potential used by
 * mofem/src/materials/impl/MatNeohookean.cpp. Require c10 > 0 and K > 0. The
 * energy is objective, so rigid rotations do not change it: zero rotational
 * modes in a full F-based Hessian are expected and are not a loss of stretch
 * stability.
 */
 
#include <Lie.hpp>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <iomanip>
#include <limits>
#include <sstream>
#include <vector>
 
namespace EshelbianPlasticity {
 
struct HMHNeohookean : public PhysicalEquations {
 
  // Abaqus-compatible compressible Neo-Hookean model
  //
  //   W(U) = c10 (J^(-2/3) tr(U^2) - 3) + K/2 (J - 1)^2,
  //   mu_0 = 2 c10, D_1 = 2/K, J = det(U).
 
  static inline double getShearModulus(const double c10) { return 2. * c10; }
 
  static inline double getAbaqusD1(const double K) { return 2. / K; }
 
  struct PrincipalState {
    double jacobian;
    double jacobianToMinusTwoThirds;
    double firstInvariant;
    double volumetricStress;
    double volumetricTangent;
  };
 
  static MoFEMErrorCode
  getCoordinateStretchFromStretch(const double stretch,
                                  double &coordinate_stretch) {
    MoFEMFunctionBegin;
    if (!std::isfinite(stretch) || stretch <= 0.) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Abaqus Neo-Hookean stretch must be finite and positive");
    }
    coordinate_stretch = std::log(stretch);
    if (!std::isfinite(coordinate_stretch)) {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Non-finite Abaqus Neo-Hookean coordinate stretch");
    }
    MoFEMFunctionReturn(0);
  }
 
  template <typename T>
  static inline double getLogJacobian(T &principal_coordinate_stretches) {
    FTENSOR_INDEX(SPACE_DIM, i);
    constexpr auto t_one = FTensor::One<>();
    return principal_coordinate_stretches(i) * t_one(i);
  }
 
  template <typename T>
  static inline PrincipalState
  getPrincipalState(const double K, T &principal_coordinate_stretches) {
    FTENSOR_INDEX(SPACE_DIM, i);
    constexpr auto t_one = FTensor::One<>();
    const double log_jacobian = getLogJacobian(principal_coordinate_stretches);
    const double jacobian = std::exp(log_jacobian);
    const double jacobian_to_minus_two_thirds =
        std::exp((-2. / 3.) * log_jacobian);
    FTensor::Tensor1<double, SPACE_DIM> t_squared_stretches;
    for (int aa = 0; aa != SPACE_DIM; ++aa) {
      t_squared_stretches(aa) =
          std::exp(2. * principal_coordinate_stretches(aa));
    }
    const double first_invariant = t_squared_stretches(i) * t_one(i);
    return {jacobian, jacobian_to_minus_two_thirds, first_invariant,
            K * jacobian * (jacobian - 1.),
            K * jacobian * (2. * jacobian - 1.)};
  }
 
  static MoFEMErrorCode validatePrincipalState(const PrincipalState &state,
                                               const char *source) {
    MoFEMFunctionBegin;
    if (!std::isfinite(state.jacobian) || state.jacobian <= 0. ||
        !std::isfinite(state.jacobianToMinusTwoThirds) ||
        state.jacobianToMinusTwoThirds <= 0. ||
        !std::isfinite(state.firstInvariant) || state.firstInvariant <= 0. ||
        !std::isfinite(state.volumetricStress) ||
        !std::isfinite(state.volumetricTangent)) {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
              "Non-finite Abaqus Neo-Hookean principal state in %s", source);
    }
    MoFEMFunctionReturn(0);
  }
 
  static inline double getPrincipalCoordinateStress(const double c10,
                                                    const PrincipalState &state,
                                                    const double v) {
    const double squared_stretch = std::exp(2. * v);
    return getShearModulus(c10) * state.jacobianToMinusTwoThirds *
               (squared_stretch - state.firstInvariant / 3.) +
           state.volumetricStress;
  }
 
  static inline double getPrincipalCoordinateStressDerivativeAtFixedInvariants(
      const double c10, const PrincipalState &state, const double v) {
    return 2. * getShearModulus(c10) * state.jacobianToMinusTwoThirds *
           std::exp(2. * v);
  }
 
  template <typename T>
  static inline double getStrainEnergy(const double c10, const double K,
                                       T &principal_coordinate_stretches) {
    const auto state = getPrincipalState(K, principal_coordinate_stretches);
    const double jacobian_minus_one = state.jacobian - 1.;
    return c10 * (state.jacobianToMinusTwoThirds * state.firstInvariant - 3.) +
           0.5 * K * jacobian_minus_one * jacobian_minus_one;
  }
 
  static MoFEMErrorCode validateMaterialParameters(const double c10,
                                                   const double K,
                                                   const char *source) {
    MoFEMFunctionBegin;
    if (!std::isfinite(c10) || c10 <= 0.) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Neo-Hookean c10 must be finite and positive in %s", source);
    }
    if (!std::isfinite(K) || K <= 0.) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Abaqus Neo-Hookean K must be finite and positive in %s "
              "(received c10=%.16g, K=%.16g)",
              source, c10, K);
    }
    MoFEMFunctionReturn(0);
  }
 
  HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
      : PhysicalEquations(), mField(m_field), c10_default(c10), K_default(K) {
 
    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,
                        "Neo-Hookean requires natural logarithmic stretch");
    }
  }
 
  auto getMaterialParameters(EntityHandle ent) {
    for (auto &b : blockData) {
      if (b.blockEnts.find(ent) != b.blockEnts.end()) {
        return std::make_pair(b.c10, b.K);
      }
    }
    if (blockData.size() != 0)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Block not found for entity handle. If you mat set "
                        "block, set it to all elements");
    return std::make_pair(c10_default, K_default);
  }
 
  struct OpJacobian : public EshelbianPlasticity::OpJacobian {
    using EshelbianPlasticity::OpJacobian::OpJacobian;
    MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
  };
 
  UserDataOperator *
  returnOpJacobian(const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
    return (new OpJacobian(eval_rhs, eval_lhs, data_ptr, physics_ptr));
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
 
    CHKERR PetscOptionsScalar("-c10", "C10", "", c10_default, &c10_default,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K_default, &K_default,
                              PETSC_NULLPTR);
 
    alphaGradU = 0;
    CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
                              alphaGradU, &alphaGradU, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    CHKERR validateMaterialParameters(c10_default, K_default,
                                      "default options");
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock Neo-Hookean (default)")
        << " c10 = " << c10_default << " K = " << K_default
        << " mu = " << getShearModulus(c10_default)
        << " Abaqus D1 = " << getAbaqusD1(K_default)
        << " grad alpha u = " << alphaGradU;
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(Sev sev) {
    return extractBlockData(
 
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
            (boost::format("%s(.*)") % "MAT_NEOHOOKEAN").str()
 
                )),
 
        sev);
  }
 
  MoFEMErrorCode
  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
                   Sev sev) {
    MoFEMFunctionBegin;
 
    for (auto m : meshset_vec_ptr) {
      MOFEM_LOG("EP", sev) << *m;
      std::vector<double> block_data;
      CHKERR m->getAttributes(block_data);
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
        return ents;
      };
 
      const auto json_parameters =
          mField.getInterface<JsonConfigManager>()->getParamsFromBlockset(
              "MAT_NEOHOOKEAN", m->getMeshsetId());
      double c10;
      double K;
      if (!json_parameters.empty()) {
        if (json_parameters.size() != 2 ||
            json_parameters.find("c10") == json_parameters.end() ||
            json_parameters.find("k") == json_parameters.end()) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "MAT_NEOHOOKEAN JSON block must have exactly two "
                  "attributes: c10, k");
        }
        c10 = json_parameters.at("c10");
        K = json_parameters.at("k");
      } else {
        if (block_data.size() < 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "MAT_NEOHOOKEAN block must have at least two attributes: "
                  "c10, K");
        }
        c10 = block_data[0];
        K = block_data[1];
      }
 
      CHKERR validateMaterialParameters(c10, K, "MAT_NEOHOOKEAN block");
 
      blockData.push_back({c10, K, get_block_ents()});
 
      MOFEM_LOG("EP", sev) << "MatBlock Neo-Hookean c10 = "
                           << blockData.back().c10
                           << " K = " << blockData.back().K
                           << " mu = " << getShearModulus(c10)
                           << " Abaqus D1 = " << getAbaqusD1(K)
                           << " nb ents. = "
                           << blockData.back().blockEnts.size();
    }
    MoFEMFunctionReturn(0);
  }
 
  struct OpSpatialPhysical : public OpAssembleVolume {
 
    OpSpatialPhysical(const std::string &field_name,
                      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                      const double alpha_u);
 
    MoFEMErrorCode integrate(EntData &data);
 
  private:
    const double alphaU;
  };
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u) {
    return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
  }
 
  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);
 
  private:
    boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
    std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
  };
 
  virtual 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) {
    return new OpSpatialPhysicalExternalStrain(field_name, data_ptr,
                                               external_strain_vec_ptr, 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));
  }
 
  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);
    MoFEMErrorCode getOptions();
    MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
  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) {
    return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
  }
 
  /**
   * @brief Recover the stretch tensor from a prescribed Biot stress.
   *
   * After diagonalising the Biot stress, the principal Hencky stretches
   * \f$h_a = \log(\lambda_a)\f$ are found from
   *
   * \f[
   * \mathcal{R}_a(\boldsymbol{h}) =
   * s_a e^{h_a} - \tau_a + q = 0, \qquad
   * \tau_a = 2c_{10}J^{-2/3}
   * \left(e^{2h_a} - \frac{I_1}{3}\right) + KJ(J-1),
   * \quad I_1 = \sum_b e^{2h_b},
   * \f]
   *
   * where \f$s_a\f$ is a principal Biot stress and
   * \f$q=3K_{\rm ext}\varepsilon_{\rm ext}\f$ is prescribed at the
   * integration point. The recovered stretch is
   * \f$\mathbf{U} = \mathbf{Q}\,\mathrm{diag}(e^{h_a})\mathbf{Q}^T\f$.
   */
  template <typename T_Biota, typename T_Stretch>
  struct CalculateStretchFromStress {
 
    MoFEMErrorCode evaluateRhs();
    MoFEMErrorCode evaluateLhs();
    MoFEMErrorCode evaluateFullLhs();
    double evaluateComplementaryPotential() const;
    MoFEMErrorCode calculateStretch();
    MoFEMErrorCode calculateBiotStretch();
    MoFEMErrorCode calculateDHenckyStreachDBiotStress();
    MoFEMErrorCode logSnesFailure(const PetscErrorCode snes_solve_error,
                                  const PetscErrorCode accepted_state_error,
                                  const PetscErrorCode reason_query_error);
    MoFEMErrorCode setUPSnes();
    static MoFEMErrorCode snesObjective(SNES snes, Vec x, PetscReal *objective,
                                        void *ctx);
    static MoFEMErrorCode snesRhs(SNES snes, Vec x, Vec r, void *ctx);
    static MoFEMErrorCode snesLhs(SNES snes, Vec x, Mat A, Mat B, void *ctx);
 
    FTensor::Tensor2_symmetric<T_Biota, 3> tBiotStress;
    FTensor::Tensor2_symmetric<T_Stretch, 3> tBiotStreach;
    FTensor::Tensor2_symmetric<double, 3> tHenckyStreach;
    FTensor::Ddg<double, 3, 3> tDRe;
    FTensor::Ddg<double, 3, 3> tDHenckyStreachDBiotStress;
    FTensor::Tensor1<double, 3> tBiotStressEigenVals;
    FTensor::Tensor1<double, 3> tHenckyStreachEigenVals;
    FTensor::Tensor2<double, 3, 3> tBiotStressEigenVecs;
    FTensor::Tensor1<double, SPACE_DIM> tReVec;
    MatrixDouble tPrincipalDReMat;
    MatrixDouble tDReMat;
    MatrixDouble tDHenckyStreachDBiotStressMat;
    /// SNES unknown: principal Hencky stretches.
    FTensor::Tensor1<double, SPACE_DIM> tStretchVec;
    SmartPetscObj<Mat> A;
    SmartPetscObj<Vec> R;
    SmartPetscObj<Vec> Chi;
    SmartPetscObj<SNES> sNes;
    double c10;
    double K;
    double externalPressure = 0;
    int max_iter = 100;
    double tol = 1e-12;
    double minimEigenValue = 0.;
    bool allowNonConverged = true;
    double minimumComplementaryHessianEigenvalue =
        std::numeric_limits<double>::infinity();
    double complementaryHessianSpectralRadius = 0.;
  };
 
  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<HMHNeohookean> neohookean_ptr,
        VectorPtr external_pressure_ptr);
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
    boost::shared_ptr<MatrixDouble> strainPtr;
    boost::shared_ptr<MatrixDouble> stressPtr;
    boost::shared_ptr<HMHNeohookean> neohookeanPtr;
    VectorPtr externalPressurePtr;
    CalculateStretchFromStress<double, double> stretchFromStress;
  };
 
  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 neohookean_ptr =
        boost::dynamic_pointer_cast<HMHNeohookean>(physics_ptr);
    if (!neohookean_ptr) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Pointer to HMHNeohookean is null");
    }
    return new OpCalculateStretchFromStress(
        data_ptr,
        strain_ptr ? strain_ptr : data_ptr->getLogStretchTensorAtPts(),
        data_ptr->getApproxPAtPts(), neohookean_ptr,
        std::move(external_pressure_ptr));
  }
 
private:
  MoFEM::Interface &mField;
 
  double c10_default;
  double K_default;
  double alphaGradU;
  double lambdaMinU = 1e-4;
  double lambdaMaxU = 1e4;
  double betaStretchBoxMin = 1e-6;
  double betaStretchBoxMax = 1e-8;
  struct BlockData {
    double c10;
    double K;
    Range blockEnts;
  };
  std::vector<BlockData> blockData;
};
 
HMHNeohookean::OpSpatialPhysical::OpSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
    : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {}
 
extern "C" {
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
                        double circumcenter[3], double *xi, double *eta,
                        double *zeta);
}
 
MoFEMErrorCode HMHNeohookean::OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
  auto [def_c10, def_K] =
      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  const double c10 = def_c10;
  const double alpha_u = alphaU;
  const double bulk_modulus = def_K;
  const double alpha_grad_u = neohookean_ptr->alphaGradU;
 
  FTENSOR_INDEX(size_symm, L);
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjoint_log_du =
      dataAtPts->getFTensorAdjointPdU(nb_integration_pts);
  auto t_u = dataAtPts->getFTensorStretch(nb_integration_pts);
  auto t_dot_log_u = dataAtPts->getFTensorLogStretchDot(nb_integration_pts);
  auto t_diff_u = dataAtPts->getFTensorDiffStretch(nb_integration_pts);
  auto t_grad_log_u =
      dataAtPts->getFTensorGradLogStretchDot(nb_integration_pts);
  auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(nb_integration_pts);
 
  auto t_eigen_vals = dataAtPts->getFTensorEigenVals(nb_integration_pts);
  auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(nb_integration_pts);
  auto &nbUniq = dataAtPts->nbUniq;
  auto t_nb_uniq =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
 
  auto t_diff = FTensor::DiffTensor<double>();
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
 
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, size_symm>, size_symm>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  auto t_grad_base_fun = data.getFTensor1DiffN<3>();
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      const auto principal_state =
          getPrincipalState(bulk_modulus, t_eigen_vals);
      CHKERR validatePrincipalState(principal_state, "spatial residual");
      auto coordinate_stress = [c10, principal_state](const double v) {
        return getPrincipalCoordinateStress(c10, principal_state, v);
      };
      auto t_neohookean_hencky =
          EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, coordinate_stress);
 
      FTensor::Tensor1<double, size_symm> t_P;
      t_P(L) = t_L(i, j, L) * t_neohookean_hencky(i, j);
 
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_nb_uniq;
 
      FTensor::Tensor1<double, size_symm> t_viscous_P;
      t_viscous_P(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_approx_P_adjoint_log_du(L) - t_P(L) - t_viscous_P(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_log_u(L, i);
      t_grad_residual(L, i) *= a;
 
      ++t_approx_P_adjoint_log_du;
      ++t_dot_log_u;
      ++t_grad_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_grad_base_fun(i) * t_grad_residual(L, i);
        ++t_nf;
        ++t_row_base_fun;
        ++t_grad_base_fun;
      }
      for (; bb != nb_base_functions; ++bb) {
        ++t_row_base_fun;
        ++t_grad_base_fun;
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto large = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Not implemented for Neo-Hookean (used ADOL-C)");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
    CHKERR large();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
  };
 
  MoFEMFunctionReturn(0);
}
 
HMHNeohookean::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(external_strain_vec_ptr), scalingMethodsMap{smv} {}
 
MoFEMErrorCode
HMHNeohookean::OpSpatialPhysicalExternalStrain::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
 
  double time = OpAssembleVolume::getFEMethod()->ts_t;
  if (EshelbianCore::physicalTimeFlg) {
    time = EshelbianCore::currentPhysicalTime;
  }
  // get entity of tet
  EntityHandle fe_ent = OpAssembleVolume::getFEEntityHandle();
  // iterate over all block data
 
  for (auto &ext_strain_block : (*externalStrainVecPtr)) {
    auto block_name = "(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(ext_strain_block.blockName, reg_name)) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Analytical external strain not implemented for Neo-Hookean "
              "material.");
    }
    // check if finite element entity is part of the EXTERNALSTRAIN block
    if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
      double scale = 1;
      if (scalingMethodsMap.find(ext_strain_block.blockName) !=
          scalingMethodsMap.end()) {
        scale *=
            scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
      } else {
        MOFEM_LOG("SELF", Sev::warning)
            << "No scaling method found for " << ext_strain_block.blockName;
      }
 
      // get ExternalStrain block data
      double external_strain_val = scale * ext_strain_block.val;
      double K = ext_strain_block.bulkModulusK;
 
      FTENSOR_INDEX(size_symm, L);
      auto t_L = FTensor::SymmLTensor<double, 3>();
      constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
      int nb_dofs = data.getIndices().size();
      int nb_integration_pts = data.getN().size1();
      auto vol = getVolume();
      auto t_w = getFTensor0IntegrationWeight();
 
      FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
 
      int nb_base_functions = data.getN().size2();
      auto t_row_base_fun = data.getFTensor0N();
 
      const double bulk_modulus = K;
      const double diag_val = external_strain_val;
      // Preserve the block's original meaning as an independent linear
      // hydrostatic eigenstress.
      const double sigma_J = bulk_modulus * 3. * diag_val;
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
        double a = vol * t_w;
        ++t_w;
 
        FTensor::Tensor1<double, size_symm> t_residual;
        t_residual(L) = (t_L(i, j, L) * t_kd(i, j)) * sigma_J;
        t_residual(L) *= a;
 
        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);
}
 
HMHNeohookean::OpSpatialPhysical_du_du::OpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
    : OpAssembleVolumePositiveDefine(row_field, col_field, data_ptr, OPROWCOL,
                                     false),
      alphaU(alpha) {
 
  CHK_THROW_MESSAGE(getOptions(), "get options failed");
 
  sYmm = false;
}
 
MoFEMErrorCode HMHNeohookean::OpSpatialPhysical_du_du::getOptions() {
  MoFEMFunctionBegin;
  PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
  CHKERR PetscOptionsScalar("-min_eigen_value", "Minimum eigenvalue", "",
                            minimEigenValue, &minimEigenValue, PETSC_NULLPTR);
  PetscOptionsEnd();
  MOFEM_LOG("EP", Sev::inform)
      << "Neo-Hookean min_eigen_value = " << minimEigenValue;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
HMHNeohookean::OpSpatialPhysical_du_du::integrate(EntData &row_data,
                                                  EntData &col_data) {
  MoFEMFunctionBegin;
 
  auto neohookean_ptr =
      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
  if (!neohookean_ptr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Pointer to HMHNeohookean is null");
  }
  auto [def_c10, def_K] =
      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  const double c10 = def_c10;
  const double alpha_u = alphaU;
  const double bulk_modulus = def_K;
  const double alpha_grad_u = neohookean_ptr->alphaGradU;
 
  FTENSOR_INDEXES(size_symm, L, J);
 
  constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
 
  auto t_L = FTensor::SymmLTensor<double, 3>();
  auto t_diff = FTensor::DiffTensor<double>();
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
 
  auto v = getVolume();
  auto ts_a = getTSa();
  auto t_w = getFTensor0IntegrationWeight();
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
  auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
 
  auto t_diff_u = dataAtPts->getFTensorDiffStretch(nb_integration_pts);
  auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(nb_integration_pts);
  auto t_u = dataAtPts->getFTensorStretch(nb_integration_pts);
  auto t_approx_P_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 &nbUniq = dataAtPts->nbUniq;
  auto t_nb_uniq =
      FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
 
  auto no_h1 = [&]() {
    MoFEMFunctionBegin;
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      double a = v * t_w;
      ++t_w;
 
      const auto principal_state =
          getPrincipalState(bulk_modulus, t_eigen_vals);
      CHKERR validatePrincipalState(principal_state, "spatial tangent");
      auto coordinate_stress = [c10, principal_state](const double v) {
        return getPrincipalCoordinateStress(c10, principal_state, v);
      };
      auto coordinate_stress_derivative = [c10,
                                           principal_state](const double v) {
        return getPrincipalCoordinateStressDerivativeAtFixedInvariants(
            c10, principal_state, v);
      };
      auto squared_stretch = [](const double v) { return std::exp(2. * v); };
      auto identity = [](const double) { return 1.; };
      auto isochoric_coordinate_stress = [c10,
                                          principal_state](const double v) {
        return getShearModulus(c10) * principal_state.jacobianToMinusTwoThirds *
               (std::exp(2. * v) - principal_state.firstInvariant / 3.);
      };
      auto t_diff_neohookean =
          EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, coordinate_stress,
                                  coordinate_stress_derivative, t_nb_uniq);
      auto t_squared_stretch =
          EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, squared_stretch);
      auto t_identity =
          EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, identity);
      auto t_isochoric_coordinate_stress = EigenMatrix::getMat(
          t_eigen_vals, t_eigen_vecs, isochoric_coordinate_stress);
      const double scaled_shear_modulus =
          getShearModulus(c10) * principal_state.jacobianToMinusTwoThirds;
      FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_material_tangent;
      t_material_tangent(i, j, k, l) =
          t_diff_neohookean(i, j, k, l) -
          (2. / 3.) * t_isochoric_coordinate_stress(i, j) * t_identity(k, l) -
          (2. / 3.) * scaled_shear_modulus * t_identity(i, j) *
              t_squared_stretch(k, l) +
          principal_state.volumetricTangent * t_identity(i, j) *
              t_identity(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, 3> t_deltaP;
      t_deltaP(i, j) = (t_approx_P_adjoint_dstretch(i, j) ||
                        t_approx_P_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_approx_P_adjoint_dstretch;
      ++t_eigen_vals;
      ++t_eigen_vecs;
      ++t_nb_uniq;
 
      t_dP(L, J) *= a;
      if (minimEigenValue > 0) {
        // Symmetrize tangesnt stiffness matrix, and add calulate eigen values,
        // tham add minimum eigen value to the tangent stiffness matrix
        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](double v) {
            return (v + eigenvalue_floor + std::abs(v - eigenvalue_floor)) / 2.;
          };
          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<3>(gg, 0);
 
        auto t_m = get_ftensor2(K, 6 * rr, 0);
        for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
          double b = t_row_base_fun * t_col_base_fun;
          t_m(L, J) += b * t_dP(L, J);
          double c = (a * alpha_grad_u * ts_a) *
                     (t_row_grad_fun(i) * t_col_grad_fun(i));
          t_m(L, J) += c * 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 large = [&]() {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Not implemented for Neo-Hookean (used ADOL-C)");
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    CHKERR no_h1();
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
    CHKERR large();
    break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "gradApproximator not handled");
  };
 
  MoFEMFunctionReturn(0);
}
 
HMHNeohookean::OpCalculateStretchFromStress::OpCalculateStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<MatrixDouble> strain_ptr,
    boost::shared_ptr<MatrixDouble> stress_ptr,
    boost::shared_ptr<HMHNeohookean> neohookean_ptr,
    VectorPtr external_pressure_ptr)
    : VolUserDataOperator(H1, OPLAST), dataAtPts(data_ptr),
      strainPtr(strain_ptr), stressPtr(stress_ptr),
      neohookeanPtr(neohookean_ptr),
      externalPressurePtr(std::move(external_pressure_ptr)) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBVERTEX] = true;
  CHK_THROW_MESSAGE(stretchFromStress.setUPSnes(),
                    "setUPSnes failed for Neo-Hookean stretch solve");
}
 
MoFEMErrorCode
HMHNeohookean::OpCalculateStretchFromStress::doWork(int side, EntityType type,
                                                    EntData &data) {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l, m, n);
  FTENSOR_INDEXES(size_symm, L, J);
 
  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 integration points");
  }
#endif // NDEBUG
  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 [def_c10, def_K] =
      neohookeanPtr->getMaterialParameters(getFEEntityHandle());
  const double c10 = def_c10;
  const double bulk_modulus = def_K;
  CHKERR validateMaterialParameters(c10, bulk_modulus,
                                    "no-stretch integration point");
 
  if (EshelbianCore::rotSelector != LARGE_ROT) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Rotation selector not handled by Abaqus Neo-Hookean");
  }
 
  auto &stretch_from_stress = stretchFromStress;
  stretch_from_stress.c10 = c10;
  stretch_from_stress.K = bulk_modulus;
  stretch_from_stress.allowNonConverged = true;
  constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
 
  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);
    stretch_from_stress.tBiotStress(i, j) =
        (t_biot_stress(i, j) || t_biot_stress(j, i)) / 2.;
    stretch_from_stress.tHenckyStreach(i, j) = t_strain(i, j);
    stretch_from_stress.externalPressure =
        externalPressurePtr ? (*externalPressurePtr)[gg] : 0.;
 
    CHKERR stretch_from_stress.calculateStretch();
    CHKERR stretch_from_stress.calculateDHenckyStreachDBiotStress();
 
    t_strain(i, j) = stretch_from_stress.tHenckyStreach(i, j);
    t_log_stretch_total(i, j) = stretch_from_stress.tHenckyStreach(i, j);
    t_log_u2_h1(i, j) = 0;
    t_biot_stretch(i, j) = stretch_from_stress.tBiotStreach(i, j);
 
    t_eigen_vals(i) = stretch_from_stress.tHenckyStreachEigenVals(i);
    t_eigen_vecs(i, j) = stretch_from_stress.tBiotStressEigenVecs(i, j);
    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);
 
    MatrixDouble d_biot_stretch_d_biot_stress_mat;
    auto get_d_biot_stretch_d_biot_stress = MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
        DL>::size(d_biot_stretch_d_biot_stress_mat, 1);
    auto t_d_biot_stretch_d_biot_stress = get_d_biot_stretch_d_biot_stress();
 
    t_d_biot_stretch_d_biot_stress(i, j, k, l) =
        t_diff_stretch(i, j, m, n) *
        stretch_from_stress.tDHenckyStreachDBiotStress(m, n, k, l);
 
    // The tensor chain above is the major-symmetric compliance C = dU/dS in
    // DDg component storage. Positive definiteness belongs to C, not to its
    // generally nonsymmetric packed action matrix Q = C W.
    d_biot_stretch_d_biot_stress_mat.resize(size_symm, size_symm, false);
    auto t_packed_compliance = getFTensor2FromMat<size_symm, size_symm>(
        d_biot_stretch_d_biot_stress_mat);
    FTensor::Tensor2<double, size_symm, size_symm> t_compliance_eig_vecs;
    t_compliance_eig_vecs(L, J) =
        0.5 * (t_packed_compliance(L, J) + t_packed_compliance(J, L));
    FTensor::Tensor1<double, size_symm> t_compliance_eig_vals;
    CHKERR computeEigenValuesSymmetric(t_compliance_eig_vecs,
                                       t_compliance_eig_vals);
    const double numerical_relative_floor =
        std::sqrt(std::numeric_limits<double>::epsilon());
    double compliance_spectral_radius = 0.;
    for (int aa = 0; aa != size_symm; ++aa) {
      compliance_spectral_radius = std::max(
          compliance_spectral_radius, std::abs(t_compliance_eig_vals(aa)));
    }
    // The material floor alone does not control the condition number. If one
    // compliance eigenvalue becomes very large, reconstructing V diag(lambda)
    // V^T can round a fixed small eigenvalue to zero. Tie the numerical floor
    // to the spectral radius so every retained mode remains representable.
    const double material_compliance_eigenvalue_floor =
        stretch_from_stress.minimEigenValue /
        std::max(std::abs(stretch_from_stress.K),
                 std::numeric_limits<double>::min());
    const double numerical_compliance_eigenvalue_floor =
        numerical_relative_floor * compliance_spectral_radius;
    const double compliance_eigenvalue_floor =
        std::max(material_compliance_eigenvalue_floor,
                 numerical_compliance_eigenvalue_floor);
    bool project_compliance = false;
    for (int aa = 0; aa != size_symm; ++aa) {
      project_compliance =
          project_compliance ||
          t_compliance_eig_vals(aa) < compliance_eigenvalue_floor;
    }
    if (project_compliance) {
      auto floor_compliance_eigenvalue =
          [compliance_eigenvalue_floor](const double v) {
            return (v + compliance_eigenvalue_floor +
                    std::abs(v - compliance_eigenvalue_floor)) /
                   2.;
          };
      auto t_projected_compliance =
          EigenMatrix::getMat(t_compliance_eig_vals, t_compliance_eig_vecs,
                              floor_compliance_eigenvalue);
      t_packed_compliance(L, J) = t_projected_compliance(L, J);
    }
 
    // Store the projected compliance directly in DDg component storage.
    auto d_biot_stress_d_hencky_mat = d_biot_stretch_d_biot_stress_mat;
    d_biot_stretch_d_biot_stress_mat.resize(1, size_symm * size_symm, false);
    t_mat_inv_d(i, j, k, l) = t_d_biot_stretch_d_biot_stress(i, j, k, l);
 
    // Tensor contraction acts through Q = C W. Invert Q, then convert Q^-1
    // back to DDg component storage Q^-1 W^-1.
    for (int row = 0; row != size_symm; ++row) {
      for (const int col : {1, 2, 4}) {
        d_biot_stress_d_hencky_mat(row, col) *= 2.;
      }
    }
    CHKERR computeMatrixInverse(d_biot_stress_d_hencky_mat);
    for (int row = 0; row != size_symm; ++row) {
      for (const int col : {1, 2, 4}) {
        d_biot_stress_d_hencky_mat(row, col) *= 0.5;
      }
    }
    auto get_d_biot_stress_d_hencky = MatrixSizeHelper<
        GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
        DL>::size(d_biot_stress_d_hencky_mat, 1);
    auto t_d_biot_stress_d_hencky = get_d_biot_stress_d_hencky();
    t_mat_d(i, j, k, l) = t_d_biot_stress_d_hencky(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);
}
 
template <typename T_Biota, typename T_Stretch>
double HMHNeohookean::CalculateStretchFromStress<
    T_Biota, T_Stretch>::evaluateComplementaryPotential() const {
  double complementary_potential = getStrainEnergy(c10, K, tStretchVec);
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    const double h = tStretchVec(aa);
    complementary_potential -= tBiotStressEigenVals(aa) * (std::exp(h) - 1.);
    complementary_potential -= externalPressure * h;
  }
  return complementary_potential;
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::evaluateRhs() {
  MoFEMFunctionBegin;
 
  const auto principal_state = getPrincipalState(K, tStretchVec);
  CHKERR validatePrincipalState(principal_state, "principal residual");
 
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    const double h = tStretchVec(aa);
    tReVec(aa) = std::exp(h) * tBiotStressEigenVals(aa) -
                 getPrincipalCoordinateStress(c10, principal_state, h) +
                 externalPressure;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::evaluateLhs() {
  MoFEMFunctionBegin;
 
  const auto principal_state = getPrincipalState(K, tStretchVec);
  CHKERR validatePrincipalState(principal_state, "principal Hessian");
  const double scaled_shear_modulus =
      getShearModulus(c10) * principal_state.jacobianToMinusTwoThirds;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  auto t_principal_d_re =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(tPrincipalDReMat);
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    const double h_a = tStretchVec(aa);
    const double stretch_a = std::exp(h_a);
    const double squared_stretch_a = stretch_a * stretch_a;
    const double isochoric_deviator_a =
        squared_stretch_a - principal_state.firstInvariant / 3.;
    for (int bb = 0; bb != SPACE_DIM; ++bb) {
      const double squared_stretch_b = std::exp(2. * tStretchVec(bb));
      const double isochoric_tangent =
          scaled_shear_modulus *
          ((aa == bb ? 2. * squared_stretch_a : 0.) -
           (2. / 3.) * squared_stretch_b - (2. / 3.) * isochoric_deviator_a);
      const double coordinate_stress_tangent =
          isochoric_tangent + principal_state.volumetricTangent;
      tPrincipalDReMat(aa, bb) = -coordinate_stress_tangent;
      if (aa == bb) {
        tPrincipalDReMat(aa, bb) += stretch_a * tBiotStressEigenVals(aa);
      }
    }
  }
 
  // R = d(S:U - W)/dh, so the complementary Hessian is -dR/dh.  Inspect the
  // unprojected Hessian even when optional Newton regularisation is enabled.
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_hessian_eig_vecs;
  t_hessian_eig_vecs(i, j) =
      -0.5 * (t_principal_d_re(i, j) + t_principal_d_re(j, i));
  FTensor::Tensor1<double, SPACE_DIM> t_hessian_eig_vals;
  CHKERR computeEigenValuesSymmetric(t_hessian_eig_vecs, t_hessian_eig_vals);
  minimumComplementaryHessianEigenvalue =
      std::numeric_limits<double>::infinity();
  complementaryHessianSpectralRadius = 0.;
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    minimumComplementaryHessianEigenvalue =
        std::min(minimumComplementaryHessianEigenvalue, t_hessian_eig_vals(aa));
    complementaryHessianSpectralRadius = std::max(
        complementaryHessianSpectralRadius, std::abs(t_hessian_eig_vals(aa)));
  }
 
  if (minimEigenValue > 0) {
    // Project the Newton matrix to positive definite and restore the residual
    // Jacobian sign. The residual and the stability measurement above remain
    // unmodified.
    const double eigenvalue_floor = K * minimEigenValue;
    bool project_hessian = false;
    for (int aa = 0; aa != SPACE_DIM; ++aa) {
      project_hessian =
          project_hessian || t_hessian_eig_vals(aa) < eigenvalue_floor;
    }
    if (project_hessian) {
      auto floor_eigenvalue = [eigenvalue_floor](const double v) {
        return (v + eigenvalue_floor + std::abs(v - eigenvalue_floor)) / 2.;
      };
      auto t_projected_hessian = EigenMatrix::getMat(
          t_hessian_eig_vals, t_hessian_eig_vecs, floor_eigenvalue);
      t_principal_d_re(i, j) = -t_projected_hessian(i, j);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota,
                                          T_Stretch>::evaluateFullLhs() {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
  FTENSOR_INDEXES(size_symm, L, J);
 
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  FTensor::Tensor1<double, SPACE_DIM> t_eigen_vals;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vecs;
  t_eigen_vecs(i, j) = tHenckyStreach(i, j);
 
  CHKERR computeEigenValuesSymmetric(t_eigen_vecs, t_eigen_vals);
 
  const int nb_uniq = getUniqNb<SPACE_DIM>(t_eigen_vals);
  if (nb_uniq < SPACE_DIM) {
    CHKERR sortEigenVals<SPACE_DIM>(t_eigen_vals, t_eigen_vecs);
  }
 
  const auto principal_state = getPrincipalState(K, t_eigen_vals);
  CHKERR validatePrincipalState(principal_state, "full Hessian");
  auto coordinate_stress = [this, principal_state](const double v) {
    return getPrincipalCoordinateStress(c10, principal_state, v);
  };
  auto coordinate_stress_derivative = [this, principal_state](const double v) {
    return getPrincipalCoordinateStressDerivativeAtFixedInvariants(
        c10, principal_state, v);
  };
  auto squared_stretch = [](const double v) { return std::exp(2. * v); };
  auto identity = [](const double) { return 1.; };
  auto isochoric_coordinate_stress = [this, principal_state](const double v) {
    return getShearModulus(c10) * principal_state.jacobianToMinusTwoThirds *
           (std::exp(2. * v) - principal_state.firstInvariant / 3.);
  };
 
  auto t_diff_neohookean =
      EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs, coordinate_stress,
                              coordinate_stress_derivative, nb_uniq);
  auto t_squared_stretch =
      EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, squared_stretch);
  auto t_identity = EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, identity);
  auto t_isochoric_coordinate_stress = EigenMatrix::getMat(
      t_eigen_vals, t_eigen_vecs, isochoric_coordinate_stress);
  const double scaled_shear_modulus =
      getShearModulus(c10) * principal_state.jacobianToMinusTwoThirds;
 
  FTensor::Ddg<double, 3, 3> t_model_log_stress_diff;
  t_model_log_stress_diff(i, j, k, l) =
      t_diff_neohookean(i, j, k, l) -
      (2. / 3.) * t_isochoric_coordinate_stress(i, j) * t_identity(k, l) -
      (2. / 3.) * scaled_shear_modulus * t_identity(i, j) *
          t_squared_stretch(k, l) +
      principal_state.volumetricTangent * t_identity(i, j) * t_identity(k, l);
 
  auto t_diff2_stretch_stress = 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), tBiotStress, nb_uniq);
 
  tDRe(i, j, k, l) =
      t_diff2_stretch_stress(i, j, k, l) - t_model_log_stress_diff(i, j, k, l);
 
  if (tDReMat.size1() != size_symm || tDReMat.size2() != size_symm) {
    tDReMat.resize(size_symm, size_symm, false);
  }
  auto t_d_re_mat = getFTensor2FromMat<size_symm, size_symm>(tDReMat);
  t_d_re_mat(L, J) = t_L(i, j, L) * (tDRe(i, j, k, l) * t_L(k, l, J));
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota,
                                          T_Stretch>::calculateBiotStretch() {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  tBiotStreach(i, j) =
      EigenMatrix::getMat(tHenckyStreachEigenVals, tBiotStressEigenVecs,
                          static_cast<double (*)(double)>(std::exp))(i, j);
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode HMHNeohookean::CalculateStretchFromStress<
    T_Biota, T_Stretch>::calculateDHenckyStreachDBiotStress() {
  MoFEMFunctionBegin;
 
  CHKERR evaluateFullLhs();
 
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
  FTENSOR_INDEXES(size_symm, L, J, M);
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  FTensor::Tensor1<double, SPACE_DIM> t_eigen_vals;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vecs;
  t_eigen_vecs(i, j) = tHenckyStreach(i, j);
 
  CHKERR computeEigenValuesSymmetric(t_eigen_vecs, t_eigen_vals);
 
  const int nb_uniq = getUniqNb<SPACE_DIM>(t_eigen_vals);
  if (nb_uniq < SPACE_DIM) {
    CHKERR sortEigenVals<SPACE_DIM>(t_eigen_vals, t_eigen_vecs);
  }
 
  auto t_diff_stretch = EigenMatrix::getDiffMat(
      t_eigen_vals, t_eigen_vecs, static_cast<double (*)(double)>(std::exp),
      static_cast<double (*)(double)>(std::exp), nb_uniq);
 
  MatrixDouble d_re_d_biot_stress_mat;
  d_re_d_biot_stress_mat.resize(size_symm, size_symm, false);
  tDHenckyStreachDBiotStressMat.resize(size_symm, size_symm, false);
 
  auto t_d_re_d_biot_stress =
      getFTensor2FromMat<size_symm, size_symm>(d_re_d_biot_stress_mat);
  t_d_re_d_biot_stress(L, J) =
      t_L(i, j, L) * (t_diff_stretch(k, l, i, j) * t_L(k, l, J));
 
  auto inv_d_re_d_hencky_mat = tDReMat;
  CHKERR computeMatrixInverse(inv_d_re_d_hencky_mat);
 
  auto t_inv_d_re_d_hencky =
      getFTensor2FromMat<size_symm, size_symm>(inv_d_re_d_hencky_mat);
  auto t_d_hencky_d_biot_stress =
      getFTensor2FromMat<size_symm, size_symm>(tDHenckyStreachDBiotStressMat);
 
  t_d_hencky_d_biot_stress(L, J) =
      -t_inv_d_re_d_hencky(L, M) * t_d_re_d_biot_stress(M, J);
 
  // tDHenckyStreachDBiotStressMat is the packed symmetric operator Q.
  // DDg storage is Q W^-1, so remove the repeated off-diagonal contraction
  // weights before exposing the same memory as a fourth-order tensor.
  auto d_hencky_d_biot_stress_ddg_mat = tDHenckyStreachDBiotStressMat;
  for (int row = 0; row != size_symm; ++row) {
    for (const int col : {1, 2, 4}) {
      d_hencky_d_biot_stress_ddg_mat(row, col) *= 0.5;
    }
  }
  auto get_d_hencky_d_biot_stress_ddg =
      MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, -1, DL>,
                       DL>::size(d_hencky_d_biot_stress_ddg_mat, 1);
  auto t_d_hencky_d_biot_stress_ddg = get_d_hencky_d_biot_stress_ddg();
  tDHenckyStreachDBiotStress(i, j, k, l) =
      t_d_hencky_d_biot_stress_ddg(i, j, k, l);
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::setUPSnes() {
  MoFEMFunctionBegin;
 
  if (sNes) {
    MoFEMFunctionReturnHot(0);
  }
 
  tPrincipalDReMat.resize(SPACE_DIM, SPACE_DIM, false);
  tDReMat.resize(size_symm, size_symm, false);
  tDHenckyStreachDBiotStressMat.resize(size_symm, size_symm, false);
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  FTENSOR_INDEXES(size_symm, L, J);
  tHenckyStreach(i, j) = 0;
  tReVec(i) = 0;
  tStretchVec(i) = 0;
  tBiotStressEigenVals(i) = 0;
  tHenckyStreachEigenVals(i) = 0;
  tBiotStressEigenVecs(i, j) = 0;
  auto t_principal_d_re =
      getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(tPrincipalDReMat);
  auto t_d_re = getFTensor2FromMat<size_symm, size_symm>(tDReMat);
  auto t_d_hencky_d_biot_stress =
      getFTensor2FromMat<size_symm, size_symm>(tDHenckyStreachDBiotStressMat);
  t_principal_d_re(i, j) = 0;
  t_d_re(L, J) = 0;
  t_d_hencky_d_biot_stress(L, J) = 0;
 
  Mat a;
  Vec r;
  Vec chi;
  CHKERR MatCreateSeqDense(PETSC_COMM_SELF, SPACE_DIM, SPACE_DIM, PETSC_NULLPTR,
                           &a);
  CHKERR VecCreateSeq(PETSC_COMM_SELF, SPACE_DIM, &r);
  CHKERR VecCreateSeq(PETSC_COMM_SELF, SPACE_DIM, &chi);
 
  A = SmartPetscObj<Mat>(a);
  R = SmartPetscObj<Vec>(r);
  Chi = SmartPetscObj<Vec>(chi);
 
  sNes = createSNES(PETSC_COMM_SELF);
  CHKERR SNESSetType(sNes, SNESNEWTONLS);
  CHKERR SNESSetObjective(sNes, snesObjective, this);
  CHKERR SNESSetFunction(sNes, R, snesRhs, this);
  CHKERR SNESSetJacobian(sNes, A, A, snesLhs, this);
 
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-neo_hookean_min_eigen_value", &minimEigenValue,
                             PETSC_NULLPTR);
  if (!std::isfinite(minimEigenValue) || minimEigenValue < 0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "-neo_hookean_min_eigen_value must be finite and non-negative");
  }
 
  KSP ksp;
  PC pc;
  SNESLineSearch line_search;
  CHKERR SNESGetKSP(sNes, &ksp);
  CHKERR KSPGetPC(ksp, &pc);
  CHKERR KSPSetType(ksp, KSPPREONLY);
  CHKERR PCSetType(pc, PCLU);
  CHKERR SNESSetTolerances(sNes, tol, tol, tol, max_iter, PETSC_DEFAULT);
  CHKERR SNESGetLineSearch(sNes, &line_search);
  CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHL2);
  PetscReal min_lambda;
  PetscReal max_step;
  PetscReal line_search_rtol;
  PetscReal line_search_atol;
  PetscReal line_search_ltol;
  PetscInt line_search_max_it;
  CHKERR SNESLineSearchGetTolerances(line_search, &min_lambda, &max_step,
                                     &line_search_rtol, &line_search_atol,
                                     &line_search_ltol, &line_search_max_it);
  CHKERR SNESLineSearchSetTolerances(line_search, min_lambda, 1.,
                                     line_search_rtol, line_search_atol,
                                     line_search_ltol, line_search_max_it);
  CHKERR SNESAppendOptionsPrefix(sNes, "nh_stretch_");
  CHKERR SNESSetFromOptions(sNes);
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::snesObjective(
    SNES, Vec x, PetscReal *objective, void *ctx) {
  MoFEMFunctionBegin;
 
  auto stretch_from_stress =
      static_cast<CalculateStretchFromStress<T_Biota, T_Stretch> *>(ctx);
 
  const double *x_array;
  CHKERR VecGetArrayRead(x, &x_array);
  for (int ii = 0; ii != SPACE_DIM; ++ii) {
    stretch_from_stress->tStretchVec(ii) = x_array[ii];
  }
  CHKERR VecRestoreArrayRead(x, &x_array);
 
  *objective = stretch_from_stress->evaluateComplementaryPotential();
  if (!std::isfinite(*objective)) {
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
            "Non-finite Neo-Hookean complementary potential");
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::snesRhs(
    SNES, Vec x, Vec r, void *ctx) {
  MoFEMFunctionBegin;
 
  auto stretch_from_stress =
      static_cast<CalculateStretchFromStress<T_Biota, T_Stretch> *>(ctx);
 
  const double *x_array;
  CHKERR VecGetArrayRead(x, &x_array);
  for (int ii = 0; ii != SPACE_DIM; ++ii) {
    stretch_from_stress->tStretchVec(ii) = x_array[ii];
  }
  CHKERR VecRestoreArrayRead(x, &x_array);
 
  CHKERR stretch_from_stress->evaluateRhs();
 
  double *r_array;
  CHKERR VecGetArray(r, &r_array);
  for (int ii = 0; ii != SPACE_DIM; ++ii) {
    // R = d(S:U-W)/dh, therefore -R is the gradient of W-S:U.
    r_array[ii] = -stretch_from_stress->tReVec(ii);
  }
  CHKERR VecRestoreArray(r, &r_array);
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::snesLhs(
    SNES, Vec x, Mat A, Mat B, void *ctx) {
  MoFEMFunctionBegin;
 
  auto stretch_from_stress =
      static_cast<CalculateStretchFromStress<T_Biota, T_Stretch> *>(ctx);
 
  const double *x_array;
  CHKERR VecGetArrayRead(x, &x_array);
  for (int ii = 0; ii != SPACE_DIM; ++ii) {
    stretch_from_stress->tStretchVec(ii) = x_array[ii];
  }
  CHKERR VecRestoreArrayRead(x, &x_array);
 
  CHKERR stretch_from_stress->evaluateLhs();
 
  PetscScalar *a_array;
  CHKERR MatDenseGetArray(A, &a_array);
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    for (int bb = 0; bb != SPACE_DIM; ++bb) {
      // Use the symmetric complementary Hessian. evaluateLhs() optionally
      // projects this Hessian before restoring the residual-Jacobian sign.
      a_array[aa + SPACE_DIM * bb] =
          -0.5 * (stretch_from_stress->tPrincipalDReMat(aa, bb) +
                  stretch_from_stress->tPrincipalDReMat(bb, aa));
    }
  }
  CHKERR MatDenseRestoreArray(A, &a_array);
 
  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  if (A != B) {
    CHKERR MatCopy(A, B, SAME_NONZERO_PATTERN);
  }
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota, T_Stretch>::logSnesFailure(
    const PetscErrorCode snes_solve_error,
    const PetscErrorCode accepted_state_error,
    const PetscErrorCode reason_query_error) {
  MoFEMFunctionBegin;
 
  // Failure reporting must never replace the error which triggered it. Keep
  // every query best-effort and report the first diagnostic error.
  PetscErrorCode diagnostic_error = 0;
  auto record_diagnostic_error = [&](const PetscErrorCode error) {
    if (error && !diagnostic_error) {
      diagnostic_error = error;
    }
    return error == 0;
  };
 
  SNESConvergedReason snes_reason = SNES_CONVERGED_ITERATING;
  PetscInt snes_iterations = -1;
  PetscReal function_norm = std::numeric_limits<PetscReal>::quiet_NaN();
  record_diagnostic_error(SNESGetConvergedReason(sNes, &snes_reason));
  record_diagnostic_error(SNESGetIterationNumber(sNes, &snes_iterations));
  record_diagnostic_error(SNESGetFunctionNorm(sNes, &function_norm));
 
  PetscMPIInt rank = -1;
  record_diagnostic_error(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 
  std::ostringstream atom_reproducer;
  atom_reproducer << std::scientific
                  << std::setprecision(
                         std::numeric_limits<double>::max_digits10)
                  << "./neohookean_stretch_gradient_atom"
                  << " -nh_stretch_atom_c10 " << c10 << " -nh_stretch_atom_k "
                  << K << " -nh_stretch_atom_q " << externalPressure
                  << " -nh_stretch_atom_max_it " << max_iter
                  << " -nh_stretch_atom_tol " << tol
                  << " -nh_stretch_atom_biot_stress " << tBiotStress(0, 0)
                  << "," << tBiotStress(0, 1) << "," << tBiotStress(0, 2) << ","
                  << tBiotStress(1, 1) << "," << tBiotStress(1, 2) << ","
                  << tBiotStress(2, 2);
 
  // setUPSnes() gives the atom the same defaults as production. Append only
  // explicitly supplied options which can change this local solve.
  const char *const reproduced_options[] = {
      "-neo_hookean_min_eigen_value",
      "-nh_stretch_test_zero_initial_guess",
      "-nh_stretch_snes_type",
      "-nh_stretch_snes_linesearch_type",
      "-nh_stretch_snes_linesearch_minlambda",
      "-nh_stretch_snes_linesearch_maxstep",
      "-nh_stretch_snes_linesearch_rtol",
      "-nh_stretch_snes_linesearch_atol",
      "-nh_stretch_snes_linesearch_ltol",
      "-nh_stretch_snes_linesearch_max_it",
      "-nh_stretch_snes_linesearch_damping",
      "-nh_stretch_snes_linesearch_order",
      "-nh_stretch_snes_linesearch_alpha",
      "-nh_stretch_ksp_type",
      "-nh_stretch_ksp_rtol",
      "-nh_stretch_ksp_atol",
      "-nh_stretch_ksp_divtol",
      "-nh_stretch_ksp_max_it",
      "-nh_stretch_ksp_gmres_restart",
      "-nh_stretch_pc_type",
      "-nh_stretch_snes_atol",
      "-nh_stretch_snes_rtol",
      "-nh_stretch_snes_stol",
      "-nh_stretch_snes_max_it",
      "-nh_stretch_snes_max_funcs",
      "-nh_stretch_snes_error_if_not_converged"};
  for (const auto option_name : reproduced_options) {
    char option_value[PETSC_MAX_PATH_LEN] = {};
    PetscBool option_set = PETSC_FALSE;
    if (record_diagnostic_error(PetscOptionsGetString(
            PETSC_NULLPTR, PETSC_NULLPTR, option_name, option_value,
            sizeof(option_value), &option_set)) &&
        option_set) {
      atom_reproducer << " " << option_name;
      if (option_value[0]) {
        atom_reproducer << " " << option_value;
      }
    }
  }
  atom_reproducer << " -log_no_color";
 
  // MoFEM's PETSc log formatter has a fixed 1024-byte buffer. Bypass it for
  // this one potentially long record so the copy-paste command is never cut.
  const auto atom_reproducer_string = atom_reproducer.str();
  std::ostringstream atom_reproducer_log;
  atom_reproducer_log << "[" << rank << "] <warning> atom_reproducer_command="
                      << atom_reproducer_string << '\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);
 
  const char *reason_string = "unavailable";
  record_diagnostic_error(SNESGetConvergedReasonString(sNes, &reason_string));
  MOFEM_LOG("SELF", Sev::warning)
      << std::scientific
      << std::setprecision(std::numeric_limits<double>::max_digits10)
      << "Neo-Hookean principal SNES failed: reason=" << reason_string << " ("
      << static_cast<int>(snes_reason) << ")"
      << " iterations=" << snes_iterations << " fnorm=" << function_norm
      << " errors=[" << snes_solve_error << "," << accepted_state_error << ","
      << reason_query_error << "," << diagnostic_error << "]";
 
  MoFEMFunctionReturn(0);
}
 
template <typename T_Biota, typename T_Stretch>
MoFEMErrorCode
HMHNeohookean::CalculateStretchFromStress<T_Biota,
                                          T_Stretch>::calculateStretch() {
  MoFEMFunctionBegin;
 
  if (!sNes) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Neo-Hookean stretch-from-stress SNES is not set up");
  }
  CHKERR validateMaterialParameters(c10, K, "stretch-from-stress solve");
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  CHKERR computeEigenValuesSymmetric(tBiotStress, tBiotStressEigenVals,
                                     tBiotStressEigenVecs);
 
  FTensor::Tensor1<double, SPACE_DIM> t_target_biot_stress_eigenvalues;
  t_target_biot_stress_eigenvalues(i) = tBiotStressEigenVals(i);
 
  constexpr auto t_one = FTensor::One<>();
  const double tr_biot_stress =
      tBiotStressEigenVals(i) * t_one(i) + SPACE_DIM * externalPressure;
  auto set_linear_initial_guess = [&]() {
    tStretchVec(i) =
        (tBiotStressEigenVals(i) + externalPressure - tr_biot_stress / 3.) /
            (4 * c10) +
        tr_biot_stress / (9 * K);
  };
  set_linear_initial_guess();
  PetscBool test_zero_initial_guess = PETSC_FALSE;
  CHKERR PetscOptionsHasName(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-nh_stretch_test_zero_initial_guess",
                             &test_zero_initial_guess);
  if (test_zero_initial_guess) {
    tStretchVec(i) = 0;
  }
  struct LocalSolveState {
    PetscErrorCode solveError = 0;
    PetscErrorCode acceptedStateError = 0;
    PetscErrorCode reasonQueryError = 0;
    PetscErrorCode residualError = 0;
    PetscErrorCode hessianError = 0;
    SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
    double residualNorm = std::numeric_limits<double>::infinity();
    bool finiteIterate = false;
    bool boundedIterate = false;
    bool stableHessian = false;
    bool acceptable = false;
  };
 
  auto solve_from_current_guess = [&]() {
    LocalSolveState state;
    double initial_residual_norm = std::numeric_limits<double>::infinity();
    if (!evaluateRhs()) {
      initial_residual_norm = tReVec.l2();
    }
 
    double *chi_array = PETSC_NULLPTR;
    state.acceptedStateError = VecGetArray(Chi, &chi_array);
    if (!state.acceptedStateError) {
      for (int aa = 0; aa != SPACE_DIM; ++aa) {
        chi_array[aa] = tStretchVec(aa);
      }
      state.acceptedStateError = VecRestoreArray(Chi, &chi_array);
    }
    if (!state.acceptedStateError) {
      state.acceptedStateError = VecAssemblyBegin(Chi);
    }
    if (!state.acceptedStateError) {
      state.acceptedStateError = VecAssemblyEnd(Chi);
    }
    if (state.acceptedStateError) {
      return state;
    }
 
    state.solveError = SNESSolve(sNes, PETSC_NULLPTR, Chi);
 
    // The callbacks can finish at a rejected trial point. Copy Chi and
    // explicitly re-evaluate the residual and unprojected Hessian there.
    const double *accepted_chi_array = PETSC_NULLPTR;
    state.acceptedStateError = VecGetArrayRead(Chi, &accepted_chi_array);
    if (!state.acceptedStateError) {
      for (int aa = 0; aa != SPACE_DIM; ++aa) {
        tStretchVec(aa) = accepted_chi_array[aa];
      }
      state.acceptedStateError = VecRestoreArrayRead(Chi, &accepted_chi_array);
    }
    state.reasonQueryError = SNESGetConvergedReason(sNes, &state.reason);
    if (!state.acceptedStateError) {
      state.residualError = evaluateRhs();
      state.hessianError = evaluateLhs();
    }
 
    state.finiteIterate = !state.acceptedStateError;
    double maximum_absolute_natural_hencky = 0.;
    for (int aa = 0; aa != SPACE_DIM && state.finiteIterate; ++aa) {
      state.finiteIterate = std::isfinite(tStretchVec(aa));
      maximum_absolute_natural_hencky =
          std::max(maximum_absolute_natural_hencky, std::abs(tStretchVec(aa)));
    }
    constexpr double catastrophic_hencky_limit = 10. * EshelbianCore::v_max;
    state.boundedIterate =
        state.finiteIterate &&
        maximum_absolute_natural_hencky <= catastrophic_hencky_limit;
 
    if (!state.residualError) {
      state.residualNorm = tReVec.l2();
    }
 
    const double stability_tolerance =
        std::sqrt(std::numeric_limits<double>::epsilon()) *
        std::max({1., std::abs(K), complementaryHessianSpectralRadius});
    state.stableHessian =
        !state.hessianError &&
        std::isfinite(minimumComplementaryHessianEigenvalue) &&
        minimumComplementaryHessianEigenvalue > stability_tolerance;
 
    PetscReal absolute_tolerance = tol;
    PetscReal relative_tolerance = tol;
    PetscReal step_tolerance = PETSC_DEFAULT;
    PetscInt maximum_iterations = max_iter;
    PetscInt maximum_function_evaluations = PETSC_DEFAULT;
    const PetscErrorCode tolerance_error = SNESGetTolerances(
        sNes, &absolute_tolerance, &relative_tolerance, &step_tolerance,
        &maximum_iterations, &maximum_function_evaluations);
    const double acceptable_residual =
        tolerance_error
            ? 10. * tol
            : 10. * std::max(static_cast<double>(absolute_tolerance),
                             static_cast<double>(relative_tolerance) *
                                 initial_residual_norm);
    const bool converged =
        !state.reasonQueryError && state.reason > SNES_CONVERGED_ITERATING;
    const bool mildly_nonconverged = allowNonConverged &&
                                     !state.reasonQueryError &&
                                     state.reason == SNES_DIVERGED_MAX_IT &&
                                     state.residualNorm <= acceptable_residual;
    const bool admissible_error =
        !state.solveError || state.solveError == PETSC_ERR_NOT_CONVERGED;
    state.acceptable =
        admissible_error && !state.acceptedStateError && !state.residualError &&
        !state.hessianError && state.finiteIterate && state.boundedIterate &&
        state.stableHessian && (converged || mildly_nonconverged);
    return state;
  };
 
  struct ContinuationState {
    LocalSolveState finalState;
    double acceptedFraction = 0.;
    double rejectedFraction = 0.;
    FTensor::Tensor1<double, SPACE_DIM> tRejectedHencky{0., 0., 0.};
  };
 
  auto run_stress_continuation = [&]() {
    ContinuationState continuation;
    FTensor::Tensor1<double, SPACE_DIM> t_accepted_hencky{0., 0., 0.};
    double fraction_step = 0.25;
    constexpr double minimum_fraction_step = 1.e-6;
    constexpr int maximum_continuation_attempts = 200;
 
    for (int attempt = 0; attempt != maximum_continuation_attempts &&
                          continuation.acceptedFraction < 1.;
         ++attempt) {
      const double trial_fraction =
          std::min(1., continuation.acceptedFraction + fraction_step);
      tBiotStressEigenVals(i) =
          trial_fraction * t_target_biot_stress_eigenvalues(i);
      tStretchVec(i) = t_accepted_hencky(i);
 
      continuation.finalState = solve_from_current_guess();
      if (continuation.finalState.acceptable) {
        continuation.acceptedFraction = trial_fraction;
        t_accepted_hencky(i) = tStretchVec(i);
        fraction_step =
            std::min(1. - continuation.acceptedFraction, 1.5 * fraction_step);
      } else {
        continuation.rejectedFraction = trial_fraction;
        continuation.tRejectedHencky(i) = tStretchVec(i);
        fraction_step *= 0.5;
        tStretchVec(i) = t_accepted_hencky(i);
      }
 
      if (continuation.acceptedFraction < 1. &&
          fraction_step < minimum_fraction_step) {
        break;
      }
    }
 
    tBiotStressEigenVals(i) = t_target_biot_stress_eigenvalues(i);
    return continuation;
  };
 
  auto direct_state = solve_from_current_guess();
  auto final_state = direct_state;
  bool used_continuation = false;
 
  if (!direct_state.acceptable) {
    used_continuation = true;
    const auto original_continuation = run_stress_continuation();
    final_state = original_continuation.finalState;
 
    if (original_continuation.acceptedFraction < 1.) {
      const double failure_minimum_hessian_eigenvalue =
          minimumComplementaryHessianEigenvalue;
      (void)logSnesFailure(final_state.solveError,
                           final_state.acceptedStateError,
                           final_state.reasonQueryError);
      const double reached_fraction = original_continuation.acceptedFraction;
      MOFEM_LOG("EP", Sev::warning)
          << "Neo-Hookean stress continuation stopped at accepted fraction "
          << reached_fraction << " after rejecting fraction "
          << original_continuation.rejectedFraction
          << "; residual = " << final_state.residualNorm
          << "; minimum complementary Hessian eigenvalue = "
          << failure_minimum_hessian_eigenvalue
          << "; rejected principal Hencky = ["
          << original_continuation.tRejectedHencky(0) << ","
          << original_continuation.tRejectedHencky(1) << ","
          << original_continuation.tRejectedHencky(2) << "]";
      MOFEM_LOG("EP", Sev::error)
          << "Abaqus Neo-Hookean stretch-from-stress solve failed to reach a "
             "globally unique constitutive target; this indicates a numerical "
             "range or nonlinear-solver failure; aborting all MPI ranks";
      // This failure is detected by one integration-point owner only. Returning
      // a rank-local PETSc error lets that rank unwind distributed solvers
      // while the other ranks remain in residual assembly, producing mismatched
      // MPI collectives (for example MUMPS MPI_Bcast versus Vec MPI_Allreduce).
      // Abort before rank-asymmetric destruction of distributed PETSc objects.
      (void)MPI_Abort(PETSC_COMM_WORLD, MOFEM_OPERATION_UNSUCCESSFUL);
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "Abaqus Neo-Hookean stretch-from-stress numerical solve failed");
    }
  }
 
  if (used_continuation) {
    MOFEM_LOG("EP", Sev::verbose)
        << "Neo-Hookean stretch-from-stress recovered by adaptive stress "
           "continuation";
  } else if (final_state.reason == SNES_DIVERGED_MAX_IT) {
    MOFEM_LOG("EP", Sev::warning)
        << "Neo-Hookean stretch-from-stress accepted a bounded, stable "
           "maximum-iteration state with residual "
        << final_state.residualNorm;
  }
 
  for (int aa = 0; aa != SPACE_DIM; ++aa) {
    if (aa != 0 &&
        isEq(tBiotStressEigenVals(aa), tBiotStressEigenVals(aa - 1))) {
      continue;
    }
    double group_sum = 0;
    int group_size = 0;
    for (int bb = aa; bb != SPACE_DIM; ++bb) {
      if (isEq(tBiotStressEigenVals(aa), tBiotStressEigenVals(bb))) {
        group_sum += tStretchVec(bb);
        ++group_size;
      }
    }
    const double group_average = group_sum / group_size;
    for (int bb = aa; bb != SPACE_DIM; ++bb) {
      if (isEq(tBiotStressEigenVals(aa), tBiotStressEigenVals(bb))) {
        tStretchVec(bb) = group_average;
      }
    }
  }
 
  tHenckyStreachEigenVals(i) = tStretchVec(i);
  auto identity = [](const double v) { return v; };
  tHenckyStreach(i, j) = EigenMatrix::getMat(
      tHenckyStreachEigenVals, tBiotStressEigenVecs, identity)(i, j);
  CHKERR calculateBiotStretch();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testHMHNeohookeanStretchGradient() {
  MoFEMFunctionBegin;
 
  using StretchFromStress =
      HMHNeohookean::CalculateStretchFromStress<double, double>;
 
  PetscReal atom_c10 = 1.7;
  PetscReal atom_K = 8.5;
  PetscReal atom_q = 0.;
  PetscReal atom_tol = 1e-13;
  PetscInt atom_max_iter = 100;
  std::array<PetscReal, size_symm + 1> atom_biot_stress_input{};
  std::array<PetscReal, size_symm> 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,
                             "-nh_stretch_atom_c10", &atom_c10, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-nh_stretch_atom_k",
                             &atom_K, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-nh_stretch_atom_q",
                             &atom_q, PETSC_NULLPTR);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR,
                            "-nh_stretch_atom_max_it", &atom_max_iter,
                            PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-nh_stretch_atom_tol", &atom_tol, PETSC_NULLPTR);
  CHKERR PetscOptionsGetRealArray(
      PETSC_NULLPTR, PETSC_NULLPTR, "-nh_stretch_atom_biot_stress",
      atom_biot_stress_input.data(), &atom_biot_stress_size,
      &atom_biot_stress_set);
  if (!std::isfinite(atom_c10) || atom_c10 <= 0. || !std::isfinite(atom_K) ||
      !std::isfinite(atom_q) || !std::isfinite(atom_tol) || atom_tol <= 0. ||
      atom_max_iter < 0 ||
      atom_max_iter > static_cast<PetscInt>(std::numeric_limits<int>::max())) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Invalid Neo-Hookean atom material or SNES parameters");
  }
  CHKERR HMHNeohookean::validateMaterialParameters(atom_c10, atom_K,
                                                   "Neo-Hookean atom test");
  if (atom_biot_stress_set && atom_biot_stress_size != size_symm) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "-nh_stretch_atom_biot_stress needs exactly six values in "
            "order s00,s01,s02,s11,s12,s22");
  }
  if (atom_biot_stress_set) {
    for (int rr = 0; rr != size_symm; ++rr) {
      const auto value = atom_biot_stress_input[rr];
      if (!std::isfinite(value)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                "-nh_stretch_atom_biot_stress values must be finite");
      }
      atom_biot_stress[rr] = value;
    }
  }
 
  EshelbianCore::stretchSelector = StretchSelector::LOG;
 
  for (const double coordinate_stretch :
       {-2., 0., EshelbianCore::v_max, EshelbianCore::v_max + 3.}) {
    const double stretch = std::exp(coordinate_stretch);
    double recovered_coordinate_stretch = 0.;
    CHKERR HMHNeohookean::getCoordinateStretchFromStretch(
        stretch, recovered_coordinate_stretch);
    const double inverse_error =
        std::abs(recovered_coordinate_stretch - coordinate_stretch);
    const double inverse_tolerance =
        1e-12 * std::max(1., std::abs(coordinate_stretch));
    if (!std::isfinite(inverse_error) || inverse_error > inverse_tolerance) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "HMHNeohookean stretch-coordinate inverse check failed");
    }
  }
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  FTENSOR_INDEX(size_symm, L);
  auto t_L = FTensor::SymmLTensor<double, 3>();
 
  auto set_symm_tensor = [&](auto &t_symm, VectorDouble &values) {
    auto t_values = getFTensor1FromPtr<size_symm>(&*values.data().begin());
    t_symm(i, j) = t_L(i, j, L) * t_values(L);
  };
 
  auto get_symm_tensor = [&](auto &t_symm, VectorDouble &values) {
    values.resize(size_symm, false);
    values[0] = t_symm(0, 0);
    values[1] = t_symm(0, 1);
    values[2] = t_symm(0, 2);
    values[3] = t_symm(1, 1);
    values[4] = t_symm(1, 2);
    values[5] = t_symm(2, 2);
  };
 
  auto set_material = [&](StretchFromStress &stretch_from_stress) {
    stretch_from_stress.c10 = atom_c10;
    stretch_from_stress.K = atom_K;
    stretch_from_stress.externalPressure = atom_q;
    stretch_from_stress.max_iter = atom_max_iter;
    stretch_from_stress.tol = atom_tol;
  };
 
  auto check_snes_solve =
      [&](StretchFromStress &stretch_from_stress, const char *case_name,
          const int col, const double perturbation, const bool log_values) {
        MoFEMFunctionBeginHot;
 
        SNESConvergedReason reason;
        const char *reason_string = "unknown";
        PetscInt iterations;
        PetscReal fnorm;
        CHKERR SNESGetConvergedReason(stretch_from_stress.sNes, &reason);
        CHKERR SNESGetConvergedReasonString(stretch_from_stress.sNes,
                                            &reason_string);
        CHKERR SNESGetIterationNumber(stretch_from_stress.sNes, &iterations);
        CHKERR SNESGetFunctionNorm(stretch_from_stress.sNes, &fnorm);
 
        const bool exceeded_built_in_iteration_limit =
            !atom_biot_stress_set && iterations > 25;
        const bool inaccurate_built_in_solution =
            !atom_biot_stress_set && fnorm > 1e-10;
        if (reason <= 0 || !std::isfinite(fnorm) ||
            inaccurate_built_in_solution || exceeded_built_in_iteration_limit) {
          MOFEM_LOG("EP", Sev::error)
              << "HMHNeohookean principal SNES failed: case=" << case_name
              << " col=" << col << " perturbation=" << perturbation
              << " iterations=" << iterations << " fnorm=" << fnorm
              << " reason=" << static_cast<int>(reason) << " (" << reason_string
              << ")";
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "HMHNeohookean principal stretch SNES did not converge");
        }
 
        if (log_values) {
          MOFEM_LOG("EP", Sev::inform)
              << "HMHNeohookean principal SNES solved: case=" << case_name
              << " iterations=" << iterations << " fnorm=" << fnorm
              << " reason=" << static_cast<int>(reason) << " (" << reason_string
              << ")\nbiot_stress=" << stretch_from_stress.tBiotStress
              << "\nhencky_stretch=" << stretch_from_stress.tHenckyStreach
              << "\nbiot_stretch=" << stretch_from_stress.tBiotStreach;
        }
 
        MoFEMFunctionReturnHot(0);
      };
 
  auto solve_for_hencky = [&](VectorDouble &stress_vec,
                              VectorDouble &hencky_vec, const char *case_name,
                              const int col, const double perturbation) {
    MoFEMFunctionBeginHot;
 
    StretchFromStress stretch_from_stress;
    set_material(stretch_from_stress);
    CHKERR stretch_from_stress.setUPSnes();
    set_symm_tensor(stretch_from_stress.tBiotStress, stress_vec);
    CHKERR stretch_from_stress.calculateStretch();
    CHKERR check_snes_solve(stretch_from_stress, case_name, col, perturbation,
                            false);
    get_symm_tensor(stretch_from_stress.tHenckyStreach, hencky_vec);
 
    MoFEMFunctionReturnHot(0);
  };
 
  struct StretchCase {
    const char *name;
    std::array<double, size_symm> stress;
    std::array<double, size_symm> expected_hencky;
    bool check_expected;
    bool check_zero_tangent;
    double tangent_tolerance;
  };
 
  const std::array<StretchCase, 8> built_in_stretch_cases{{
      {"generic",
       {0.18, -0.035, 0.024, -0.095, 0.041, 0.13},
       {},
       false,
       false,
       5e-5},
      {"zero", {0., 0., 0., 0., 0., 0.}, {}, true, true, 5e-5},
      {"hydrostatic",
       {3.6322073565896487, 0., 0., 3.6322073565896487, 0., 3.6322073565896487},
       {0.1, 0., 0., 0.1, 0., 0.1},
       true,
       false,
       5e-5},
      {"rotated_repeated",
       {0.09259999671274205, -0.73713612813466134, -0.36856806406733089,
        0.092599996712742161, 0.368568064067331, -0.4602520993882544},
       {0.02666666666666665, -0.1066666666666666, -0.05333333333333333,
        0.02666666666666666, 0.05333333333333333, -0.05333333333333334},
       true,
       false,
       5e-5},
      {"rotated_near_repeated",
       {0.092600176677642032, -0.73713594355526935, -0.36856807329630148,
        0.092600112074853747, 0.36856820250187705, -0.46025239471529633},
       {0.02666669266666667, -0.10666664, -0.05333333466666666,
        0.02666668333333332, 0.05333335333333333, -0.05333337600000001},
       true,
       false,
       1e-4},
      {"rotated_distinct",
       {0.8407464902540408, 0.82227899319182596, 0.28074304918139947,
        0.75776693280844076, -0.1147839342902007, 0.66914993941471534},
       {0.044, 0.14, 0.048, 0.03, -0.02, 0.016},
       true,
       false,
       5e-5},
      {"tensile_distortional",
       {3.037800922108802, 0., 0., 1.1493777171719499, 0., 1.7405883229709418},
       {0.3, 0., 0., -0.12, 0., 0.02},
       true,
       false,
       5e-5},
      {"compressive_distortional",
       {-3.6858407112404157, 0., 0., 0.067487441652110958, 0.,
        -1.1528550088055771},
       {-0.3, 0., 0., 0.12, 0., -0.02},
       true,
       false,
       5e-5},
  }};
 
  std::vector<StretchCase> stretch_cases;
  if (atom_biot_stress_set) {
    StretchCase command_line_case{"command_line", {}, {}, false, false, 5e-5};
    for (int rr = 0; rr != size_symm; ++rr) {
      command_line_case.stress[rr] = atom_biot_stress[rr];
    }
    stretch_cases.push_back(command_line_case);
  } else {
    stretch_cases.assign(built_in_stretch_cases.begin(),
                         built_in_stretch_cases.end());
  }
  const bool use_reference_configuration =
      atom_c10 == 1.7 && atom_K == 8.5 && atom_q == 0.;
 
  auto make_vector = [](const std::array<double, size_symm> &values) {
    VectorDouble vector;
    vector.resize(size_symm, false);
    for (int rr = 0; rr != size_symm; ++rr) {
      vector[rr] = values[rr];
    }
    return vector;
  };
 
  auto commutator_error = [](auto &t_a, auto &t_b) {
    double commutator_norm_squared = 0.;
    double a_norm_squared = 0.;
    double b_norm_squared = 0.;
    for (int rr = 0; rr != SPACE_DIM; ++rr) {
      for (int cc = 0; cc != SPACE_DIM; ++cc) {
        double commutator = 0.;
        for (int kk = 0; kk != SPACE_DIM; ++kk) {
          commutator += t_a(rr, kk) * t_b(kk, cc) - t_b(rr, kk) * t_a(kk, cc);
        }
        commutator_norm_squared += commutator * commutator;
        a_norm_squared += t_a(rr, cc) * t_a(rr, cc);
        b_norm_squared += t_b(rr, cc) * t_b(rr, cc);
      }
    }
    return std::sqrt(commutator_norm_squared) /
           std::max(1., std::sqrt(a_norm_squared * b_norm_squared));
  };
 
  constexpr double eps = 1e-6;
  for (const auto &stretch_case : stretch_cases) {
    auto stress_vec = make_vector(stretch_case.stress);
 
    StretchFromStress base_stretch_from_stress;
    set_material(base_stretch_from_stress);
    CHKERR base_stretch_from_stress.setUPSnes();
    set_symm_tensor(base_stretch_from_stress.tBiotStress, stress_vec);
    CHKERR base_stretch_from_stress.calculateStretch();
    CHKERR check_snes_solve(base_stretch_from_stress, stretch_case.name, -1, 0.,
                            true);
    CHKERR base_stretch_from_stress.calculateDHenckyStreachDBiotStress();
 
    VectorDouble base_hencky_vec;
    get_symm_tensor(base_stretch_from_stress.tHenckyStreach, base_hencky_vec);
 
    for (int rr = 0; rr != SPACE_DIM; ++rr) {
      for (int cc = rr; cc != SPACE_DIM; ++cc) {
        if (!std::isfinite(base_stretch_from_stress.tBiotStress(rr, cc)) ||
            !std::isfinite(base_stretch_from_stress.tHenckyStreach(rr, cc)) ||
            !std::isfinite(base_stretch_from_stress.tBiotStreach(rr, cc))) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "HMHNeohookean spectral reconstruction is not finite");
        }
      }
    }
 
    if (atom_q == 0. && stretch_case.check_expected &&
        (use_reference_configuration || stretch_case.check_zero_tangent)) {
      double max_expected_error = 0.;
      for (int rr = 0; rr != size_symm; ++rr) {
        const double expected_error =
            std::abs(base_hencky_vec[rr] - stretch_case.expected_hencky[rr]);
        if (!std::isfinite(expected_error)) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "HMHNeohookean spectral reconstruction is not finite");
        }
        max_expected_error = std::max(max_expected_error, expected_error);
      }
      MOFEM_LOG("EP", Sev::inform)
          << "HMHNeohookean spectral reconstruction check: case="
          << stretch_case.name << " max_abs=" << max_expected_error;
      if (max_expected_error > 1e-9) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "HMHNeohookean spectral reconstruction check failed");
      }
    }
 
    const double stress_hencky_commutator =
        commutator_error(base_stretch_from_stress.tBiotStress,
                         base_stretch_from_stress.tHenckyStreach);
    const double stress_biot_commutator =
        commutator_error(base_stretch_from_stress.tBiotStress,
                         base_stretch_from_stress.tBiotStreach);
    const double hencky_biot_commutator =
        commutator_error(base_stretch_from_stress.tHenckyStreach,
                         base_stretch_from_stress.tBiotStreach);
    const double max_commutator =
        std::max(stress_hencky_commutator,
                 std::max(stress_biot_commutator, hencky_biot_commutator));
    if (!std::isfinite(max_commutator) || max_commutator > 1e-10) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "HMHNeohookean coaxial tensor check failed");
    }
 
    FTensor::Tensor1<double, SPACE_DIM> t_hencky_eigen_vals;
    FTensor::Tensor1<double, SPACE_DIM> t_biot_eigen_vals;
    FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vecs;
    t_eigen_vecs(i, j) = base_stretch_from_stress.tHenckyStreach(i, j);
    CHKERR computeEigenValuesSymmetric(t_eigen_vecs, t_hencky_eigen_vals);
    auto t_expected_biot_stretch =
        EigenMatrix::getMat(t_hencky_eigen_vals, t_eigen_vecs,
                            static_cast<double (*)(double)>(std::exp));
    for (int rr = 0; rr != SPACE_DIM; ++rr) {
      for (int cc = rr; cc != SPACE_DIM; ++cc) {
        const double biot_stretch_error =
            std::abs(base_stretch_from_stress.tBiotStreach(rr, cc) -
                     t_expected_biot_stretch(rr, cc));
        if (!std::isfinite(biot_stretch_error) || biot_stretch_error > 1e-10) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "HMHNeohookean Biot stretch does not match the selected "
                  "stretch map");
        }
      }
    }
    t_eigen_vecs(i, j) = base_stretch_from_stress.tBiotStreach(i, j);
    CHKERR computeEigenValuesSymmetric(t_eigen_vecs, t_biot_eigen_vals);
    for (int aa = 0; aa != SPACE_DIM; ++aa) {
      if (!std::isfinite(t_biot_eigen_vals(aa)) ||
          t_biot_eigen_vals(aa) <= 0.) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "HMHNeohookean Biot stretch is not positive definite");
      }
    }
 
    MatrixDouble fd_gradient;
    fd_gradient.resize(size_symm, size_symm, false);
    for (int col = 0; col != size_symm; ++col) {
      auto stress_plus_vec = stress_vec;
      auto stress_minus_vec = stress_vec;
      stress_plus_vec[col] += eps;
      stress_minus_vec[col] -= eps;
 
      VectorDouble plus_hencky_vec;
      VectorDouble minus_hencky_vec;
      CHKERR solve_for_hencky(stress_plus_vec, plus_hencky_vec,
                              stretch_case.name, col, eps);
      CHKERR solve_for_hencky(stress_minus_vec, minus_hencky_vec,
                              stretch_case.name, col, -eps);
 
      for (int row = 0; row != size_symm; ++row) {
        fd_gradient(row, col) =
            (plus_hencky_vec[row] - minus_hencky_vec[row]) / (2. * eps);
      }
    }
 
    double max_abs_error = 0.;
    double max_scaled_error = 0.;
    double max_analytical = 0.;
    double max_numerical = 0.;
    int max_row = 0;
    int max_col = 0;
    for (int row = 0; row != size_symm; ++row) {
      for (int col = 0; col != size_symm; ++col) {
        const double analytical =
            base_stretch_from_stress.tDHenckyStreachDBiotStressMat(row, col);
        const double numerical = fd_gradient(row, col);
        if (!std::isfinite(analytical) || !std::isfinite(numerical)) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "HMHNeohookean full tangent is not finite");
        }
        const double abs_error = std::abs(analytical - numerical);
        const double scale =
            std::max(1., std::max(std::abs(analytical), std::abs(numerical)));
        const double scaled_error = abs_error / scale;
 
        if (scaled_error > max_scaled_error) {
          max_scaled_error = scaled_error;
          max_abs_error = abs_error;
          max_analytical = analytical;
          max_numerical = numerical;
          max_row = row;
          max_col = col;
        }
      }
    }
 
    MOFEM_LOG("EP", Sev::inform)
        << "HMHNeohookean full tangent finite difference check: case="
        << stretch_case.name << " max_abs=" << max_abs_error
        << " max_scaled=" << max_scaled_error
        << " analytical=" << max_analytical << " numerical=" << max_numerical
        << " at (" << max_row << "," << max_col << ")";
 
    if (max_scaled_error > stretch_case.tangent_tolerance) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "HMHNeohookean full tangent finite difference check failed");
    }
 
    if (atom_q == 0. && stretch_case.check_zero_tangent) {
      const double a = 1. / (4. * atom_c10);
      const double q = 1. / (9. * atom_K) - a / 3.;
      MatrixDouble expected_tangent;
      expected_tangent.resize(size_symm, size_symm, false);
      expected_tangent.clear();
      for (const int row : {0, 3, 5}) {
        for (const int col : {0, 3, 5}) {
          expected_tangent(row, col) = q;
        }
        expected_tangent(row, row) += a;
      }
      for (const int shear : {1, 2, 4}) {
        expected_tangent(shear, shear) = a;
      }
 
      double max_zero_tangent_error = 0.;
      for (int row = 0; row != size_symm; ++row) {
        for (int col = 0; col != size_symm; ++col) {
          max_zero_tangent_error = std::max(
              max_zero_tangent_error,
              std::abs(base_stretch_from_stress.tDHenckyStreachDBiotStressMat(
                           row, col) -
                       expected_tangent(row, col)));
        }
      }
      if (!std::isfinite(max_zero_tangent_error) ||
          max_zero_tangent_error > 1e-10) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "HMHNeohookean zero-stress analytical tangent check failed");
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity