Recover the stretch tensor from a prescribed Biot stress. More...

Inheritance diagram for EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >:

Collaboration diagram for EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >:

Public Member Functions
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 Public Member Functions
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)

Public Attributes
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

FTensor::Tensor1< double, SPACE_DIM >	tStretchVec
	SNES unknown: principal Hencky stretches.

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

double	complementaryHessianSpectralRadius = 0.

Detailed Description

template<typename T_Biota, typename T_Stretch>
struct EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >

Recover the stretch tensor from a prescribed Biot stress.

After diagonalising the Biot stress, the principal Hencky stretches \(h_a = \log(\lambda_a)\) are found from

\[ \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}, \]

where \(s_a\) is a principal Biot stress and \(q=3K_{\rm ext}\varepsilon_{\rm ext}\) is prescribed at the integration point. The recovered stretch is \(\mathbf{U} = \mathbf{Q}\,\mathrm{diag}(e^{h_a})\mathbf{Q}^T\).

Definition at line 383 of file HMHNeohookean.cpp.

Member Function Documentation

◆ calculateBiotStretch()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::calculateBiotStretch ( )

Definition at line 1408 of file HMHNeohookean.cpp.

                                                                           {
  MoFEMFunctionBegin;
 
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  tBiotStreach(i, j) =
      EigenMatrix::getMat(tHenckyStreachEigenVals, tBiotStressEigenVecs,
                          static_cast<double (*)(double)>(std::exp))(i, j);
 
  MoFEMFunctionReturn(0);
}

◆ calculateDHenckyStreachDBiotStress()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::calculateDHenckyStreachDBiotStress ( )

Definition at line 1421 of file HMHNeohookean.cpp.

                                                            {
  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);
}

◆ calculateStretch()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::calculateStretch ( )

Definition at line 1775 of file HMHNeohookean.cpp.

                                                                       {
  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);
}

◆ evaluateComplementaryPotential()

template<typename T_Biota , typename T_Stretch >

double EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::evaluateComplementaryPotential ( ) const

Definition at line 1229 of file HMHNeohookean.cpp.

                                                              {
  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;
}

◆ evaluateFullLhs()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::evaluateFullLhs ( )

Definition at line 1334 of file HMHNeohookean.cpp.

                                                                      {
  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);
}

◆ evaluateLhs()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::evaluateLhs ( )

Definition at line 1259 of file HMHNeohookean.cpp.

                                                                         {
  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);
}

◆ evaluateRhs()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::evaluateRhs ( )

Definition at line 1241 of file HMHNeohookean.cpp.

                                                                         {
  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);
}

◆ logSnesFailure()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::logSnesFailure	(	const PetscErrorCode	snes_solve_error,
		const PetscErrorCode	accepted_state_error,
		const PetscErrorCode	reason_query_error
	)

Definition at line 1663 of file HMHNeohookean.cpp.

                                             {
  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);
}

◆ setUPSnes()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::setUPSnes ( )

Definition at line 1486 of file HMHNeohookean.cpp.

                                                                       {
  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);
}

◆ snesLhs()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::snesLhs	(	SNES	snes,
		Vec	x,
		Mat	A,
		Mat	B,
		void *	ctx
	)

static

Definition at line 1623 of file HMHNeohookean.cpp.

                                          {
  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);
}

◆ snesObjective()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::snesObjective	(	SNES	snes,
		Vec	x,
		PetscReal *	objective,
		void *	ctx
	)

static

Definition at line 1569 of file HMHNeohookean.cpp.

                                                  {
  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);
}

◆ snesRhs()

template<typename T_Biota , typename T_Stretch >

MoFEMErrorCode EshelbianPlasticity::HMHNeohookean::CalculateStretchFromStress< T_Biota, T_Stretch >::snesRhs	(	SNES	snes,
		Vec	x,
		Vec	r,
		void *	ctx
	)

static

Definition at line 1594 of file HMHNeohookean.cpp.

                                   {
  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);
}