neohookean_auxiliary_log_stress_atom.cpp

Go to the documentation of this file.

/**
 * @file neohookean_auxiliary_log_stress_atom.cpp
 * @brief Verify the five-coordinate Neo-Hookean material map.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <NeoHookeanLogarithmicMaterial.hpp>
using Material = EshelbianPlasticity::NeoHookeanLogarithmicMaterial;
 
namespace {
 
MoFEMErrorCode checkNear(const double actual, const double expected,
                         const double tolerance, const char *label) {
  MoFEMFunctionBegin;
  if (!std::isfinite(actual) ||
      std::abs(actual - expected) >
          tolerance * std::max(1., std::abs(expected)))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "%s: expected %.16e, obtained %.16e", label, expected, actual);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkCoordinates(const Material::Coordinates &actual,
                                const Material::Coordinates &expected,
                                const double tolerance, const char *label) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  double scale = 1.;
  for (int aa = 0; aa != 5; ++aa)
    scale = std::max(scale, std::abs(expected(aa)));
  Material::Coordinates t_error;
  t_error(A) = (actual(A) - expected(A)) / scale;
  CHKERR checkNear(std::sqrt(t_error(A) * t_error(A)), 0., tolerance, label);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkPositiveSymmetric(const Material::Tangent &t_matrix,
                                      const double tolerance,
                                      const char *label) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
  Material::Tangent t_error;
  t_error(A, B) = t_matrix(A, B) - t_matrix(B, A);
  CHKERR checkNear(std::sqrt(t_error(A, B) * t_error(A, B)), 0.,
                   tolerance *
                       std::max(1., std::sqrt(t_matrix(A, B) * t_matrix(A, B))),
                   label);
  Material::Tangent t_eigen_vectors;
  t_eigen_vectors(A, B) = t_matrix(A, B);
  Material::Coordinates t_eigen_values;
  CHKERR computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values);
  for (int aa = 0; aa != 5; ++aa)
    if (!std::isfinite(t_eigen_values(aa)) || t_eigen_values(aa) <= 0.)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "%s: eigenvalue %d is %.16e", label, aa, t_eigen_values(aa));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkReference(const Material::Parameters &parameters,
                              const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  Material::Coordinates t_zero;
  t_zero(A) = 0.;
  Material::Coordinates t_stress;
  Material::Tangent t_hessian;
  double energy;
  CHKERR Material::evaluateDeviator(parameters, t_zero, energy, &t_stress,
                                    &t_hessian);
  CHKERR checkNear(energy, 0., tolerance, "reference deviatoric energy");
  CHKERR checkCoordinates(t_stress, t_zero, tolerance, "reference stress");
  Material::InverseState inverse;
  CHKERR Material::evaluateInverse(parameters, t_zero, inverse);
  CHKERR checkCoordinates(inverse.tMaterialDeviator, t_zero, tolerance,
                          "reference inverse");
  CHKERR checkNear(inverse.conjugateEnergy, 0., tolerance,
                   "reference conjugate energy");
  constexpr auto t_delta = FTensor::Kronecker_Delta<int>();
  Material::Tangent t_error;
  t_error(A, B) = t_hessian(A, B) - 4. * parameters.c10 * t_delta(A, B);
  CHKERR checkNear(std::sqrt(t_error(A, B) * t_error(A, B)), 0., tolerance,
                   "reference Hessian");
  t_error(A, B) =
      inverse.tCompliance(A, B) - t_delta(A, B) / (4. * parameters.c10);
  CHKERR checkNear(std::sqrt(t_error(A, B) * t_error(A, B)), 0., tolerance,
                   "reference compliance");
 
  Material::VolumeState volume;
  CHKERR Material::evaluateVolume(parameters, 0., volume);
  CHKERR checkNear(volume.jacobian, 1., tolerance, "reference Jacobian");
  CHKERR checkNear(volume.energy, 0., tolerance, "reference volume energy");
  CHKERR checkNear(volume.firstDerivative, 0., tolerance,
                   "reference volume derivative");
  CHKERR checkNear(volume.secondDerivative, parameters.bulkModulus, tolerance,
                   "reference volume Hessian");
  CHKERR Material::evaluateVolume(parameters, std::log(0.25), volume);
  CHKERR checkNear(volume.secondDerivative, -parameters.bulkModulus / 8.,
                   tolerance, "unprojected compressed volume Hessian");
 
  // A very small K keeps the response finite even near the largest finite J.
  CHKERR Material::evaluateVolume(
      {parameters.c10, std::numeric_limits<double>::denorm_min()}, 709.7,
      volume);
  CHKERR checkNear(volume.energy / volume.firstDerivative, 0.5, tolerance,
                   "large finite volume energy");
  CHKERR checkNear(volume.secondDerivative / volume.firstDerivative, 2.,
                   tolerance, "large finite volume Hessian");
 
  for (int aa = 0; aa != 5; ++aa) {
    Material::Coordinates t_coordinate;
    t_coordinate(A) = 0.;
    t_coordinate(aa) = 1.;
    auto t_basis = Material::getTensor(t_coordinate);
    CHKERR checkNear(t_basis(i, i), 0., tolerance, "basis trace");
    CHKERR checkNear(t_basis(i, j) * t_basis(i, j), 1., tolerance,
                     "basis Frobenius norm");
    CHKERR checkCoordinates(Material::getCoordinates(t_basis), t_coordinate,
                            tolerance, "basis coordinate round trip");
    for (int bb = 0; bb != aa; ++bb) {
      t_coordinate(A) = 0.;
      t_coordinate(bb) = 1.;
      const auto t_other_basis = Material::getTensor(t_coordinate);
      CHKERR checkNear(t_basis(i, j) * t_other_basis(i, j), 0., tolerance,
                       "orthogonal tensor coordinates");
    }
  }
  const FTensor::Tensor2<double, 3, 3> t_full(2., 3., 1., -1., -1., 4., 5., 0.,
                                              2.);
  const auto t_projected =
      MoFEM::Tensor2SymmetricDeviatorBasis::getCoordinates(t_full);
  const double sqrt_two = std::sqrt(2.);
  const Material::Coordinates t_expected(3. / sqrt_two, -3. / std::sqrt(6.),
                                         sqrt_two, 3. * sqrt_two,
                                         2. * sqrt_two);
  CHKERR checkCoordinates(t_projected, t_expected, tolerance,
                          "full tensor symmetric deviatoric projection");
  CHKERR checkNear(t_projected(A) * t_projected(A), 34., tolerance,
                   "full tensor projected Frobenius norm");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkMaterialState(const Material::Parameters &parameters,
                                  const Material::Coordinates &t_deviator,
                                  const double fd_step,
                                  const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
  FTensor::Index<'C', 5> C;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  Material::Coordinates t_stress;
  Material::Tangent t_hessian;
  double energy;
  CHKERR Material::evaluateDeviator(parameters, t_deviator, energy, &t_stress,
                                    &t_hessian);
  Material::InverseState inverse;
  CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
  CHKERR checkCoordinates(inverse.tMaterialDeviator, t_deviator, tolerance,
                          "inverse round trip");
  CHKERR checkNear(inverse.deviatoricEnergy, energy, tolerance,
                   "inverse material energy");
  CHKERR checkNear(energy + inverse.conjugateEnergy,
                   t_stress(A) * t_deviator(A), tolerance, "Fenchel identity");
  double gap;
  CHKERR Material::evaluateFenchelGap(parameters, t_deviator, t_stress, inverse,
                                      gap);
  CHKERR checkNear(gap, 0., tolerance, "zero gap at the material copy");
  CHKERR checkPositiveSymmetric(t_hessian, tolerance, "deviatoric Hessian");
  CHKERR checkPositiveSymmetric(inverse.tCompliance, tolerance, "compliance");
 
  constexpr auto t_delta = FTensor::Kronecker_Delta<int>();
  Material::Tangent t_inverse_error;
  t_inverse_error(A, B) =
      t_hessian(A, C) * inverse.tCompliance(C, B) - t_delta(A, B);
  CHKERR checkNear(std::sqrt(t_inverse_error(A, B) * t_inverse_error(A, B)), 0.,
                   tolerance, "Hessian-compliance product");
 
  Material::Coordinates t_direction(0.31, -0.27, 0.41, 0.19, -0.37);
  const double direction_norm = std::sqrt(t_direction(B) * t_direction(B));
  t_direction(A) = t_direction(A) / direction_norm;
  const auto t_d = Material::getTensor(t_deviator);
  const auto t_v = Material::getTensor(t_direction);
  FTensor::Tensor2<double, 3, 3> t_commutator;
  t_commutator(i, j) = t_d(i, k) * t_v(k, j) - t_v(i, k) * t_d(k, j);
  if (std::sqrt(t_deviator(A) * t_deviator(A)) > 0.01 &&
      std::sqrt(t_commutator(i, j) * t_commutator(i, j)) < 1.e-3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Material derivative test direction must be noncoaxial");
 
  Material::Coordinates t_plus, t_minus, t_stress_plus, t_stress_minus;
  t_plus(A) = t_deviator(A) + fd_step * t_direction(A);
  t_minus(A) = t_deviator(A) - fd_step * t_direction(A);
  double energy_plus, energy_minus;
  CHKERR Material::evaluateDeviator(parameters, t_plus, energy_plus,
                                    &t_stress_plus);
  CHKERR Material::evaluateDeviator(parameters, t_minus, energy_minus,
                                    &t_stress_minus);
  CHKERR checkNear((energy_plus - energy_minus) / (2. * fd_step),
                   t_stress(A) * t_direction(A), tolerance,
                   "energy directional derivative");
  Material::Coordinates t_fd, t_exact;
  t_fd(A) = (t_stress_plus(A) - t_stress_minus(A)) / (2. * fd_step);
  t_exact(A) = t_hessian(A, B) * t_direction(B);
  CHKERR checkCoordinates(t_fd, t_exact, tolerance,
                          "forward stress directional derivative");
 
  t_plus(A) = t_stress(A) + fd_step * t_direction(A);
  t_minus(A) = t_stress(A) - fd_step * t_direction(A);
  Material::InverseState inverse_plus, inverse_minus;
  CHKERR Material::evaluateInverse(parameters, t_plus, inverse_plus);
  CHKERR Material::evaluateInverse(parameters, t_minus, inverse_minus);
  t_fd(A) =
      (inverse_plus.tMaterialDeviator(A) - inverse_minus.tMaterialDeviator(A)) /
      (2. * fd_step);
  t_exact(A) = inverse.tCompliance(A, B) * t_direction(B);
  CHKERR checkCoordinates(t_fd, t_exact, tolerance,
                          "inverse noncoaxial directional derivative");
  CHKERR checkNear(
      (inverse_plus.conjugateEnergy - inverse_minus.conjugateEnergy) /
          (2. * fd_step),
      t_deviator(A) * t_direction(A), tolerance,
      "conjugate energy directional derivative");
  t_plus(A) = t_deviator(A) + 0.2 * t_direction(A);
  CHKERR Material::evaluateFenchelGap(parameters, t_plus, t_stress, inverse,
                                      gap);
  if (gap <= 0.)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Distinct material copies must have a positive Fenchel gap");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkSpectralStates(const Material::Parameters &parameters,
                                   const double fd_step,
                                   const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Tensor2<double, 3, 3> t_rotation(0.36, -0.8, 0.48, 0.48, 0.6, 0.64,
                                            -0.8, 0., 0.6);
  for (const auto &diagonal : {std::array<double, 3>{0., 0., 0.},
                               {0.12, 0.12, -0.24},
                               {0.12, 0.120000001, -0.240000001},
                               {0.3, -0.12, -0.18}}) {
    Material::SymmetricTensor t_diagonal;
    t_diagonal(i, j) = 0.;
    for (int aa = 0; aa != 3; ++aa)
      t_diagonal(aa, aa) = diagonal[aa];
    auto t_coordinates = Material::getCoordinates(t_diagonal);
    CHKERR checkMaterialState(parameters, t_coordinates, fd_step, tolerance);
    FTensor::Tensor2<double, 3, 3> t_rotated_left;
    t_rotated_left(i, j) = t_rotation(i, k) * t_diagonal(k, j);
    Material::SymmetricTensor t_rotated;
    t_rotated(i, j) = (t_rotated_left(i, k) * t_rotation(j, k) ||
                       t_rotated_left(j, k) * t_rotation(i, k)) /
                      2.;
    auto t_rotated_coordinates = Material::getCoordinates(t_rotated);
    CHKERR checkMaterialState(parameters, t_rotated_coordinates, fd_step,
                              tolerance);
 
    Material::Coordinates t_stress, t_rotated_stress;
    double energy, rotated_energy;
    CHKERR Material::evaluateDeviator(parameters, t_coordinates, energy,
                                      &t_stress);
    CHKERR Material::evaluateDeviator(parameters, t_rotated_coordinates,
                                      rotated_energy, &t_rotated_stress);
    CHKERR checkNear(rotated_energy, energy, tolerance,
                     "rotation-invariant energy");
    auto t_s = Material::getTensor(t_stress);
    t_rotated_left(i, j) = t_rotation(i, k) * t_s(k, j);
    t_rotated(i, j) = (t_rotated_left(i, k) * t_rotation(j, k) ||
                       t_rotated_left(j, k) * t_rotation(i, k)) /
                      2.;
    CHKERR checkCoordinates(t_rotated_stress,
                            Material::getCoordinates(t_rotated), tolerance,
                            "rotated material stress");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkForwardReferences(const Material::Parameters &parameters,
                                      const double fd_step,
                                      const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  struct Reference {
    std::array<double, 3> diagonal;
    double energy;
    std::array<double, 3> piola;
  };
  const std::array<Reference, 3> references{
      {{{0.1, 0.1, 0.1},
        0.520205037263637,
        {3.63220735658965, 3.63220735658965, 3.63220735658965}},
       {{0.3, -0.12, 0.02},
        0.538123520405792,
        {3.03780092210880, 1.14937771717195, 1.74058832297094}},
       {{-0.3, 0.12, -0.02},
        0.441373550717383,
        {-3.68584071124042, 0.0674874416521112, -1.15285500880558}}}};
  for (const auto &reference : references) {
    Material::SymmetricTensor t_h;
    t_h(i, j) = 0.;
    for (int aa = 0; aa != 3; ++aa)
      t_h(aa, aa) = reference.diagonal[aa];
    const double theta = t_h(i, i);
    auto t_deviator = Material::getCoordinates(t_h);
    auto t_reconstructed = Material::getTensor(t_deviator, theta);
    Material::SymmetricTensor t_error;
    t_error(i, j) = t_reconstructed(i, j) - t_h(i, j);
    CHKERR checkNear(std::sqrt(t_error(i, j) * t_error(i, j)), 0., tolerance,
                     "reconstructed logarithmic stretch");
    Material::Coordinates t_stress;
    double energy;
    CHKERR Material::evaluateDeviator(parameters, t_deviator, energy,
                                      &t_stress);
    Material::VolumeState volume, volume_plus, volume_minus;
    CHKERR Material::evaluateVolume(parameters, theta, volume);
    CHKERR Material::evaluateVolume(parameters, theta + fd_step, volume_plus);
    CHKERR Material::evaluateVolume(parameters, theta - fd_step, volume_minus);
    CHKERR checkNear(energy + volume.energy, reference.energy, tolerance,
                     "independent forward energy");
    const auto t_s = Material::getTensor(t_stress);
    for (int aa = 0; aa != 3; ++aa)
      CHKERR checkNear((t_s(aa, aa) + volume.firstDerivative) /
                           std::exp(reference.diagonal[aa]),
                       reference.piola[aa], tolerance,
                       "independent principal Piola stress");
    CHKERR checkNear(
        (volume_plus.energy - volume_minus.energy) / (2. * fd_step),
        volume.firstDerivative, tolerance, "volumetric energy derivative");
    CHKERR checkNear(
        (volume_plus.firstDerivative - volume_minus.firstDerivative) /
            (2. * fd_step),
        volume.secondDerivative, tolerance, "volumetric tangent");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkLargeStress(const Material::Parameters &parameters,
                                const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  const Material::Coordinates t_direction(1., -0.5, 0.25, -0.75, 0.4);
  for (const double scale : {1.e4, 1.e8, 1.e100, 1.e200}) {
    Material::Coordinates t_stress;
    t_stress(A) = scale * t_direction(A);
    Material::InverseState inverse;
    CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
    double energy;
    Material::Coordinates t_recovered;
    CHKERR Material::evaluateDeviator(parameters, inverse.tMaterialDeviator,
                                      energy, &t_recovered);
    CHKERR checkCoordinates(t_recovered, t_stress, tolerance,
                            "large finite stress round trip");
    CHKERR checkNear(inverse.inverseResidual, 0., tolerance,
                     "large finite stress inverse residual");
    CHKERR checkNear(energy + inverse.conjugateEnergy,
                     t_stress(A) * inverse.tMaterialDeviator(A), tolerance,
                     "large stress Fenchel identity");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkSmallDeviator(const Material::Parameters &parameters,
                                  const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  const Material::Coordinates t_deviator(0.6e-10, -0.3e-10, 0.2e-10, -0.4e-10,
                                         0.5e-10);
  Material::Coordinates t_stress, t_error;
  double energy;
  CHKERR Material::evaluateDeviator(parameters, t_deviator, energy, &t_stress);
  const double norm_squared = t_deviator(A) * t_deviator(A);
  const double mu = 2. * parameters.c10;
  CHKERR checkNear(energy / (mu * norm_squared), 1., tolerance,
                   "small-deviator relative energy");
  t_error(A) = t_stress(A) / (2. * mu) - t_deviator(A);
  CHKERR checkNear(std::sqrt((t_error(A) * t_error(A)) / norm_squared), 0.,
                   tolerance, "small-deviator relative stress");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkEnergyScaling(const double tolerance) {
  MoFEMFunctionBegin;
  const Material::Parameters parameters{5.e307, 1.};
  const Material::SymmetricTensor t_d(0.65, 0., 0., -0.325, 0., -0.325);
  double energy;
  CHKERR Material::evaluateDeviator(parameters, Material::getCoordinates(t_d),
                                    energy);
  CHKERR checkNear(energy / 8.56694110570638e307, 1., tolerance,
                   "finite energy despite large exponential term");
  const Material::Coordinates t_small(1.e-200, 0., 0., 0., 0.);
  CHKERR Material::evaluateDeviator(parameters, t_small, energy);
  CHKERR checkNear(energy / 1.e-92, 1., tolerance,
                   "finite energy despite tiny squared strain");
  const Material::SymmetricTensor t_distorted(0.55, 0., 0., -0.275, 0., -0.275);
  Material::Coordinates t_stress;
  CHKERR Material::evaluateDeviator(
      parameters, Material::getCoordinates(t_distorted), energy, &t_stress);
  CHKERR checkNear(t_stress(0) / 1.7163010440184161e308, 1., tolerance,
                   "finite projected principal stress");
  CHKERR checkNear(t_stress(1) / 9.909068697744685e307, 1., tolerance,
                   "finite projected second stress coordinate");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkConditioning(const Material::Parameters &parameters,
                                 const double tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'A', 5> A;
  FTensor::Index<'B', 5> B;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  double previous_condition = 0.;
  for (const double amplitude : {0., 1., 2.}) {
    Material::SymmetricTensor t_diagonal;
    t_diagonal(i, j) = 0.;
    t_diagonal(0, 0) = -amplitude;
    t_diagonal(1, 1) = -amplitude;
    t_diagonal(2, 2) = 2. * amplitude;
    const auto t_deviator = Material::getCoordinates(t_diagonal);
    Material::Coordinates t_stress;
    double energy;
    CHKERR Material::evaluateDeviator(parameters, t_deviator, energy,
                                      &t_stress);
    Material::InverseState inverse;
    CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
    CHKERR checkCoordinates(inverse.tMaterialDeviator, t_deviator, tolerance,
                            "conditioned inverse round trip");
    CHKERR checkNear(inverse.tCompliance(2, 2),
                     std::exp(2. * amplitude) / (4. * parameters.c10),
                     tolerance, "repeated-plane shear compliance");
    Material::Tangent t_eigen_vectors;
    t_eigen_vectors(A, B) = inverse.tCompliance(A, B);
    Material::Coordinates t_eigen_values;
    CHKERR computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values);
    const double condition = t_eigen_values(4) / t_eigen_values(0);
    if (!std::isfinite(condition) || condition <= previous_condition)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Expected increasing compliance condition, got %.16e", condition);
    CHKERR PetscPrintf(PETSC_COMM_WORLD,
                       "Auxiliary inverse a=%.1f: compliance condition %.8e\n",
                       amplitude, condition);
    previous_condition = condition;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkInvalidInputs(const Material::Parameters &parameters) {
  MoFEMFunctionBegin;
  const double infinity = std::numeric_limits<double>::infinity();
  const double nan = std::numeric_limits<double>::quiet_NaN();
  for (const auto bad : {Material::Parameters{0., parameters.bulkModulus},
                         {-1., parameters.bulkModulus},
                         {parameters.c10, 0.},
                         {parameters.c10, -1.},
                         {nan, parameters.bulkModulus},
                         {parameters.c10, infinity}}) {
    CHKERR PetscPushErrorHandler(PetscReturnErrorHandler, nullptr);
    const auto error = Material::validateParameters(bad);
    CHKERR PetscPopErrorHandler();
    if (!error)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Invalid material parameters were accepted");
  }
  Material::Coordinates t_invalid(nan, 0., 0., 0., 0.);
  Material::InverseState inverse;
  Material::VolumeState volume;
  double energy;
  CHKERR PetscPushErrorHandler(PetscReturnErrorHandler, nullptr);
  const auto inverse_error =
      Material::evaluateInverse(parameters, t_invalid, inverse);
  const auto forward_error =
      Material::evaluateDeviator(parameters, t_invalid, energy);
  const auto volume_error =
      Material::evaluateVolume(parameters, infinity, volume);
  CHKERR PetscPopErrorHandler();
  if (!inverse_error || !forward_error || !volume_error)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Non-finite material input was accepted");
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] =
    "Neo-Hookean auxiliary logarithmic stress material atom.\n"
    "Use -json_config neohookean_auxiliary_log_stress_atom.json.\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
  try {
    moab::Core moab;
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
    auto *json_config = m_field.getInterface<JsonConfigManager>();
    auto *meshsets_manager = m_field.getInterface<MeshsetsManager>();
    // Only material parameters are needed; this atom has no element geometry.
    CHKERR meshsets_manager->addMeshset(BLOCKSET, 1001, "MAT_NEOHOOKEAN");
    CHKERR meshsets_manager->setMeshsetFromFile();
    const auto block =
        json_config->getParamsFromBlockset("MAT_NEOHOOKEAN", 1001);
    if (block.size() != 2 || !block.count("c10") || !block.count("k"))
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Atom requires MAT_NEOHOOKEAN JSON block 1001 with c10,k");
    const Material::Parameters parameters{block.at("c10"), block.at("k")};
    CHKERR Material::validateParameters(parameters);
    CHKERR checkNear(parameters.c10, 1.7, 1.e-14, "reference C10");
    CHKERR checkNear(parameters.bulkModulus, 8.5, 1.e-14, "reference K");
    PetscReal fd_step = 0., tolerance = 0.;
    CHKERR PetscOptionsGetReal(nullptr, nullptr, "-auxiliary_atom_fd_step",
                               &fd_step, nullptr);
    CHKERR PetscOptionsGetReal(nullptr, nullptr, "-auxiliary_atom_tolerance",
                               &tolerance, nullptr);
    if (!std::isfinite(fd_step) || fd_step <= 0. || !std::isfinite(tolerance) ||
        tolerance <= 0.)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Atom requires positive finite JSON finite-difference settings");
    CHKERR checkReference(parameters, tolerance);
    CHKERR checkReference({parameters.c10, 0.1}, tolerance);
    CHKERR checkSpectralStates(parameters, fd_step, tolerance);
    CHKERR checkForwardReferences(parameters, fd_step, tolerance);
    CHKERR checkLargeStress(parameters, tolerance);
    CHKERR checkSmallDeviator(parameters, tolerance);
    CHKERR checkEnergyScaling(tolerance);
    CHKERR checkConditioning(parameters, tolerance);
    CHKERR checkInvalidInputs(parameters);
    CHKERR PetscPrintf(
        PETSC_COMM_WORLD,
        "Neo-Hookean auxiliary logarithmic stress atom passed\n");
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}