auxiliary_logarithmic_stress_space_atom.cpp

Go to the documentation of this file.

/**
 * @file auxiliary_logarithmic_stress_space_atom.cpp
 * @brief Weak mesh-stress closure, observability and coefficient elimination.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <EshelbianPlasticity.hpp>
#include <NeoHookeanLogarithmicMaterial.hpp>
#include <IntegrationRules.hpp>
 
namespace {
 
using Material = EshelbianPlasticity::NeoHookeanLogarithmicMaterial;
using Coordinates = Material::Coordinates;
using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
 
MoFEMErrorCode checkNear(const double value, const double expected,
                         const double tolerance, const char *description) {
  MoFEMFunctionBegin;
  if (!std::isfinite(value) ||
      std::abs(value - expected) > tolerance * std::max(1., std::abs(expected)))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "%s: expected %.16e, obtained %.16e", description, expected, value);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode makeQuadrature(const int order, MatrixDouble &points,
                             VectorDouble &weights) {
  MoFEMFunctionBegin;
  const auto *rule = IntRules::XiaoGimbutas::getTetrahedronRule(order);
  if (!rule || rule->numBarycentricCoordinates != 4)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Space atom requires a tetrahedral Xiao-Gimbutas rule");
  points.resize(3, rule->numPoints, false);
  weights.resize(rule->numPoints, false);
  double sum = 0.;
  for (unsigned int gg = 0; gg != rule->numPoints; ++gg) {
    points(0, gg) = rule->points[4 * gg + 1];
    points(1, gg) = rule->points[4 * gg + 2];
    points(2, gg) = rule->points[4 * gg + 3];
    weights[gg] = rule->weights[gg] / 6.;
    if (!(weights[gg] > 0.) || !std::isfinite(weights[gg]))
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Space comparisons require strictly positive quadrature weights");
    sum += weights[gg];
  }
  CHKERR checkNear(sum, 1. / 6., 1.e-13, "reference tetrahedron volume");
  MoFEMFunctionReturn(0);
}
 
/** Evaluate the production USER_BASE directly on the reference tetrahedron. */
MoFEMErrorCode evaluateBasis(MatrixDouble &points, const int order,
                             MatrixDouble &basis) {
  MoFEMFunctionBegin;
  EntitiesFieldData data(MBTET);
  data.dataOnEntities[MBTET][0].getOrder() = order;
  auto context = boost::make_shared<EntPolynomialBaseCtx>(
      data, L2, DISCONTINUOUS, USER_BASE);
  EshelbianPlasticity::CGGUserPolynomialBase user_base(nullptr, true);
  CHKERR user_base.getValue(points, context);
  basis = data.dataOnEntities[MBTET][0].getN(USER_BASE);
  if (basis.size2() != (order + 1) * (order + 2) * (order + 3) / 6)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Unexpected tetrahedral USER_BASE dimension");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode makeKinematics(const MatrixDouble &basis, MatrixDouble &values,
                             const bool constant = false) {
  MoFEMFunctionBegin;
  FTensor::Index<'a', 5> a;
  const Coordinates t_constant(.12, -.07, .04, .03, -.05);
  const Coordinates t_first(.22, .09, -.13, .07, .11);
  const Coordinates t_second(-.08, .18, .07, -.10, .05);
  const Coordinates t_third(.12, -.15, .09, .03, -.08);
  auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::size(
      values, basis.size1())();
  for (unsigned int gg = 0; gg != basis.size1(); ++gg) {
    t_d(a) = t_constant(a);
    if (!constant)
      t_d(a) += basis(gg, 1) * t_first(a) + basis(gg, 2) * t_second(a) +
                basis(gg, 3) * t_third(a);
    ++t_d;
  }
  MoFEMFunctionReturn(0);
}
 
struct CoefficientState {
  VectorDouble residual;
  MatrixDouble compliance;
  double gap = 0.;
  double minimumGap = 0.;
  double maximumMismatch = 0.;
};
 
/** Dense coefficient oracle, with tensor components contracted by FTensor.
 * These arrays represent the small reference-cell system, not global FE
 * assembly. Basis indices remain explicit and stress is interpolated before
 * every material inverse; it is never replaced by the pointwise direct law.
 */
MoFEMErrorCode evaluateCoefficients(
    const Material::Parameters &parameters, const MatrixDouble &basis,
    const VectorDouble &weights, MatrixDouble &kinematics,
    VectorDouble &coefficients, CoefficientState &state, const bool tangent) {
  MoFEMFunctionBegin;
  const int nb_basis = basis.size2();
  const int nb_gauss = basis.size1();
  if (coefficients.size() != 5 * nb_basis || weights.size() != nb_gauss)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Inconsistent coefficient-space test data");
  state.residual.resize(5 * nb_basis, false);
  state.residual.clear();
  if (tangent) {
    state.compliance.resize(5 * nb_basis, 5 * nb_basis, false);
    state.compliance.clear();
  }
  state.gap = 0.;
  state.minimumGap = std::numeric_limits<double>::infinity();
  state.maximumMismatch = 0.;
  FTensor::Index<'a', 5> a;
  FTensor::Index<'b', 5> b;
  auto t_kinematic =
      MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
          kinematics, nb_gauss)();
  for (int gg = 0; gg != nb_gauss; ++gg) {
    Coordinates t_stress, t_d, t_mismatch;
    t_stress(a) = 0.;
    auto t_coefficient = getFTensor1FromArray<5, 5>(coefficients);
    for (int bb = 0; bb != nb_basis; ++bb) {
      t_stress(a) += basis(gg, bb) * t_coefficient(a);
      ++t_coefficient;
    }
    t_d(a) = t_kinematic(a);
    Material::InverseState inverse;
    CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
    double gap;
    CHKERR Material::evaluateFenchelGap(parameters, t_d, t_stress, inverse, gap);
    t_mismatch(a) = t_d(a) - inverse.tMaterialDeviator(a);
    state.gap += weights[gg] * gap;
    state.minimumGap = std::min(state.minimumGap, gap);
    state.maximumMismatch = std::max(
        state.maximumMismatch, std::sqrt(t_mismatch(a) * t_mismatch(a)));
    auto t_residual = getFTensor1FromArray<5, 5>(state.residual);
    for (int rr = 0; rr != nb_basis; ++rr) {
      const double alpha = weights[gg] * basis(gg, rr);
      t_residual(a) += alpha * t_mismatch(a);
      if (tangent && alpha != 0.) {
        auto t_c = getFTensor2FromArray<5, 5, 5>(state.compliance, 5 * rr);
        for (int cc = 0; cc != nb_basis; ++cc) {
          t_c(a, b) += (alpha * basis(gg, cc)) * inverse.tCompliance(a, b);
          ++t_c;
        }
      }
      ++t_residual;
    }
    ++t_kinematic;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode solveCoefficients(
    const Material::Parameters &parameters, const MatrixDouble &basis,
    const VectorDouble &weights, MatrixDouble &kinematics,
    VectorDouble &coefficients, const double tolerance,
    CoefficientState &state) {
  MoFEMFunctionBegin;
  for (int iteration = 0; iteration != 30; ++iteration) {
    CHKERR evaluateCoefficients(parameters, basis, weights, kinematics,
                                coefficients, state, true);
    const double residual_norm = norm_2(state.residual);
    if (residual_norm < tolerance)
      MoFEMFunctionReturnHot(0);
    VectorDouble step = state.residual;
    CHKERR solveLinearSystem(static_cast<const MatrixDouble &>(state.compliance),
                             step);
    bool accepted = false;
    for (double alpha = 1.; alpha >= 1.e-7; alpha *= .5) {
      VectorDouble trial = coefficients + alpha * step;
      CoefficientState trial_state;
      CHKERR evaluateCoefficients(parameters, basis, weights, kinematics, trial,
                                  trial_state, false);
      if (norm_2(trial_state.residual) <=
          (1. - 1.e-4 * alpha) * residual_norm) {
        coefficients.swap(trial);
        accepted = true;
        break;
      }
    }
    if (!accepted)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Stress coefficient Newton line search failed");
  }
  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
          "Stress coefficient moments did not converge");
  MoFEMFunctionReturn(0);
}
 
MatrixDouble crossMass(const MatrixDouble &row_basis,
                       const MatrixDouble &col_basis,
                       const VectorDouble &weights) {
  MatrixDouble matrix(5 * row_basis.size2(), 5 * col_basis.size2());
  matrix.clear();
  FTensor::Index<'a', 5> a;
  FTensor::Index<'b', 5> b;
  constexpr auto t_identity = FTensor::Kronecker_Delta<int>();
  for (unsigned int rr = 0; rr != row_basis.size2(); ++rr) {
    auto t_m = getFTensor2FromArray<5, 5, 5>(matrix, 5 * rr);
    for (unsigned int cc = 0; cc != col_basis.size2(); ++cc) {
      double pairing = 0.;
      for (unsigned int gg = 0; gg != weights.size(); ++gg)
        pairing += weights[gg] * row_basis(gg, rr) * col_basis(gg, cc);
      t_m(a, b) = pairing * t_identity(a, b);
      ++t_m;
    }
  }
  return matrix;
}
 
MoFEMErrorCode checkPositiveDefinite(const MatrixDouble &matrix,
                                     const char *description) {
  MoFEMFunctionBegin;
  VectorDouble eigenvalues;
  MatrixDouble eigenvectors;
  CHKERR computeEigenValuesSymmetric(matrix, eigenvalues, eigenvectors);
  if (!(eigenvalues[0] > 1.e-12 * eigenvalues[eigenvalues.size() - 1]))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "%s is not observable with the selected positive quadrature",
            description);
  const double condition = eigenvalues[eigenvalues.size() - 1] / eigenvalues[0];
  MatrixDouble scaled(matrix);
  for (unsigned int row = 0; row != scaled.size1(); ++row)
    for (unsigned int col = 0; col != scaled.size2(); ++col)
      scaled(row, col) /= std::sqrt(matrix(row, row) * matrix(col, col));
  CHKERR computeEigenValuesSymmetric(scaled, eigenvalues, eigenvectors);
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "%s (%d coefficients): condition %.6e, "
                     "diagonally scaled condition %.6e\n",
                     description, static_cast<int>(matrix.size1()), condition,
                     eigenvalues[eigenvalues.size() - 1] / eigenvalues[0]);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkExactCases(
    const Material::Parameters &parameters, const MatrixDouble &constant_basis,
    const VectorDouble &weights, MatrixDouble &kinematics,
    const double tolerance, const double newton_tolerance) {
  MoFEMFunctionBegin;
  FTensor::Index<'a', 5> a;
  MatrixDouble constant_kinematics;
  CHKERR makeKinematics(constant_basis, constant_kinematics, true);
  VectorDouble coefficients(5);
  coefficients.clear();
  CoefficientState state;
  CHKERR solveCoefficients(parameters, constant_basis, weights,
                           constant_kinematics, coefficients, newton_tolerance,
                           state);
  CHKERR checkNear(state.maximumMismatch, 0., tolerance, "P0 constant copy");
  CHKERR checkNear(state.gap, 0., tolerance, "P0 constant Fenchel gap");
  const Coordinates t_constant(.12, -.07, .04, .03, -.05);
  Coordinates t_exact, t_error;
  double energy;
  CHKERR Material::evaluateDeviator(parameters, t_constant, energy, &t_exact);
  auto t_coefficient = getFTensor1FromArray<5, 5>(coefficients);
  t_error(a) = t_coefficient(a) - t_exact(a);
  CHKERR checkNear(std::sqrt(t_error(a) * t_error(a)), 0., tolerance,
                   "P0 coefficient direct-law recovery");
 
  // Only this separate reference space has independent quadrature copies.
  const int nb_gauss = weights.size();
  MatrixDouble independent_basis(nb_gauss, nb_gauss);
  independent_basis.clear();
  for (int gg = 0; gg != nb_gauss; ++gg)
    independent_basis(gg, gg) = 1.;
  coefficients.resize(5 * nb_gauss, false);
  coefficients.clear();
  CHKERR solveCoefficients(parameters, independent_basis, weights, kinematics,
                           coefficients, newton_tolerance, state);
  CHKERR checkNear(state.maximumMismatch, 0., tolerance,
                   "independent quadrature material copies");
  CHKERR checkNear(state.gap, 0., tolerance, "independent-copy Fenchel gap");
  auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
      kinematics, nb_gauss)();
  auto t_independent = getFTensor1FromArray<5, 5>(coefficients);
  for (int gg = 0; gg != nb_gauss; ++gg) {
    Coordinates t_kinematic;
    t_kinematic(a) = t_d(a);
    CHKERR Material::evaluateDeviator(parameters, t_kinematic, energy, &t_exact);
    t_error(a) = t_independent(a) - t_exact(a);
    CHKERR checkNear(std::sqrt(t_error(a) * t_error(a)), 0., tolerance,
                     "independent-copy direct-law recovery");
    ++t_d;
    ++t_independent;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkLostMode(
    const Material::Parameters &parameters, const MatrixDouble &constant_basis,
    const MatrixDouble &linear_basis, const VectorDouble &weights,
    const double tolerance) {
  MoFEMFunctionBegin;
  const auto cross = crossMass(constant_basis, linear_basis, weights);
  CHKERR checkPositiveDefinite(crossMass(linear_basis, linear_basis, weights),
                               "P1 kinematic mass");
  // P0 stress sees five averages, whereas P1 kinematics has twenty DOFs.
  VectorDouble lost_mode(20);
  lost_mode.clear();
  lost_mode[5 + 2] = 1.;
  lost_mode[2] = -cross(2, 5 + 2) / cross(2, 2);
  const VectorDouble weak = prod(cross, lost_mode);
  CHKERR checkNear(norm_2(weak), 0., tolerance, "underresolved cross-mass mode");
  const auto mass = crossMass(linear_basis, linear_basis, weights);
  const VectorDouble mass_mode = prod(mass, lost_mode);
  if (!(inner_prod(lost_mode, mass_mode) > 1.e-3))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "The allegedly lost kinematic mode has zero physical norm");
  for (unsigned int gg = 0; gg != weights.size(); ++gg) {
    const double value = linear_basis(gg, 1) +
                         lost_mode[2] * linear_basis(gg, 0);
    Coordinates t_d(0., 0., value, 0., 0.), t_stress, t_error;
    double energy;
    CHKERR Material::evaluateDeviator(parameters, t_d, energy, &t_stress);
    Material::InverseState inverse;
    CHKERR Material::evaluateInverse(parameters, t_stress, inverse);
    FTensor::Index<'a', 5> a;
    t_error(a) = inverse.tMaterialDeviator(a) - t_d(a);
    CHKERR checkNear(std::sqrt(t_error(a) * t_error(a)), 0., tolerance,
                     "material inverse on the unobservable kinematic mode");
  }
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "Underresolved P0/P1 pairing: rank 5 of 20, "
                     "nonzero deviatoric mode invisible; all inverses valid\n");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkCoefficientElimination(
    const Material::Parameters &parameters, const MatrixDouble &basis,
    const VectorDouble &weights, MatrixDouble &kinematics,
    const VectorDouble &equilibrated_stress, const double tolerance) {
  MoFEMFunctionBegin;
  VectorDouble stress = .6 * equilibrated_stress;
  CoefficientState state;
  CHKERR evaluateCoefficients(parameters, basis, weights, kinematics, stress,
                              state, true);
  if (!(norm_2(state.residual) > 1.e-4))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Elimination test requires a nonzero copy residual");
  const MatrixDouble coupling = crossMass(basis, basis, weights);
  // A positive mass block supplies a fixed mechanical test problem. The
  // constitutive row, its compliance and both cross blocks use the actual
  // mesh-stress discretisation at the deliberately unequilibrated state.
  const MatrixDouble mechanical = 1.3 * coupling;
  const int size = stress.size();
  VectorDouble mechanical_residual(size);
  for (int ii = 0; ii != size; ++ii)
    mechanical_residual[ii] = .01 * std::sin(ii + 1.);
 
  MatrixDouble enlarged(2 * size, 2 * size);
  VectorDouble explicit_step(2 * size);
  for (int rr = 0; rr != size; ++rr) {
    explicit_step[rr] = -mechanical_residual[rr];
    explicit_step[size + rr] = -state.residual[rr];
    for (int cc = 0; cc != size; ++cc) {
      enlarged(rr, cc) = mechanical(rr, cc);
      enlarged(rr, size + cc) = coupling(rr, cc);
      enlarged(size + rr, cc) = coupling(cc, rr);
      enlarged(size + rr, size + cc) = -state.compliance(rr, cc);
    }
  }
  CHKERR solveLinearSystem(enlarged, explicit_step);
 
  MatrixDouble inverse_coupling(size, size);
  for (int cc = 0; cc != size; ++cc) {
    VectorDouble column(size);
    for (int rr = 0; rr != size; ++rr)
      column[rr] = coupling(cc, rr);
    CHKERR solveLinearSystem(static_cast<const MatrixDouble &>(state.compliance),
                             column);
    for (int rr = 0; rr != size; ++rr)
      inverse_coupling(rr, cc) = column[rr];
  }
  VectorDouble inverse_residual = state.residual;
  CHKERR solveLinearSystem(static_cast<const MatrixDouble &>(state.compliance),
                           inverse_residual);
  MatrixDouble schur = mechanical + prod(coupling, inverse_coupling);
  VectorDouble rhs_correction = prod(coupling, inverse_residual);
  VectorDouble condensed_d = -mechanical_residual - rhs_correction;
  CHKERR solveLinearSystem(static_cast<const MatrixDouble &>(schur), condensed_d);
  const VectorDouble condensed_t =
      inverse_residual + prod(inverse_coupling, condensed_d);
  VectorDouble difference = explicit_step;
  for (int ii = 0; ii != size; ++ii) {
    difference[ii] -= condensed_d[ii];
    difference[size + ii] -= condensed_t[ii];
  }
  CHKERR checkNear(norm_2(difference) / std::max(1., norm_2(explicit_step)), 0.,
                   tolerance, "explicit/coefficient-Schur Newton step");
  VectorDouble omitted_correction = -mechanical_residual;
  CHKERR solveLinearSystem(static_cast<const MatrixDouble &>(schur),
                           omitted_correction);
  omitted_correction -= condensed_d;
  if (!(norm_2(omitted_correction) > 1.e-4))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Omitting the nonzero copy-residual correction did not change "
            "the condensed Newton step");
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "Coefficient Schur step agrees; copy residual %.6e, "
                     "omitted-correction error %.6e\n",
                     norm_2(state.residual), norm_2(omitted_correction));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runSpaceAtom(const Material::Parameters &parameters,
                           const int quadrature_order, const double tolerance,
                           const double newton_tolerance) {
  MoFEMFunctionBegin;
  MatrixDouble points;
  VectorDouble weights;
  CHKERR makeQuadrature(quadrature_order, points, weights);
  std::array<MatrixDouble, 3> bases;
  for (int order = 0; order != 3; ++order) {
    CHKERR evaluateBasis(points, order, bases[order]);
    CHKERR checkPositiveDefinite(crossMass(bases[order], bases[order], weights),
                                 "nested stress-space mass");
    if (order)
      for (unsigned int gg = 0; gg != weights.size(); ++gg)
        for (unsigned int bb = 0; bb != bases[order - 1].size2(); ++bb)
          CHKERR checkNear(bases[order](gg, bb), bases[order - 1](gg, bb),
                           tolerance, "actual USER_BASE nesting");
  }
  MatrixDouble kinematics;
  CHKERR makeKinematics(bases[1], kinematics);
  CHKERR checkExactCases(parameters, bases[0], weights, kinematics, tolerance,
                         newton_tolerance);
  VectorDouble coefficients, previous, linear_solution;
  double previous_gap = std::numeric_limits<double>::infinity();
  for (int order = 0; order != 3; ++order) {
    coefficients.resize(5 * bases[order].size2(), false);
    coefficients.clear();
    std::copy(previous.begin(), previous.end(), coefficients.begin());
    CoefficientState state;
    CHKERR solveCoefficients(parameters, bases[order], weights, kinematics,
                             coefficients, newton_tolerance, state);
    CHKERR checkPositiveDefinite(state.compliance,
                                 "assembled stress coefficient compliance");
    if (state.minimumGap < -tolerance || state.gap < -tolerance ||
        (order && !(state.gap < previous_gap - 1.e-3 * tolerance)))
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Nested mesh-stress Fenchel gaps are negative or fail to decrease");
    if (order == 1) {
      linear_solution = coefficients;
      if (!(state.maximumMismatch > 1.e-5))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Weak P1 closure unexpectedly imposed pointwise material copies");
    }
    CHKERR PetscPrintf(PETSC_COMM_WORLD,
                       "Mesh stress P%d: %d coefficients, moment %.6e, "
                       "gap %.12e, maximum copy mismatch %.6e\n",
                       order, static_cast<int>(coefficients.size()),
                       norm_2(state.residual), state.gap, state.maximumMismatch);
    previous = coefficients;
    previous_gap = state.gap;
  }
  CHKERR checkLostMode(parameters, bases[0], bases[1], weights, tolerance);
  CHKERR checkCoefficientElimination(parameters, bases[1], weights, kinematics,
                                     linear_solution, tolerance);
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "Auxiliary logarithmic stress space atom passed "
                     "with %d fixed positive quadrature points\n",
                     static_cast<int>(weights.size()));
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] =
    "Check weak auxiliary stress spaces and coefficient Schur elimination.\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 = m_field.getInterface<JsonConfigManager>();
    auto *meshsets = m_field.getInterface<MeshsetsManager>();
    CHKERR meshsets->addMeshset(BLOCKSET, 1001, "MAT_NEOHOOKEAN");
    CHKERR meshsets->setMeshsetFromFile();
    const auto block = json->getParamsFromBlockset("MAT_NEOHOOKEAN", 1001);
    if (block.size() != 2 || !block.count("c10") || !block.count("k"))
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Space atom requires MAT_NEOHOOKEAN JSON block 1001 with c10,k");
    const Material::Parameters parameters{block.at("c10"), block.at("k")};
    CHKERR Material::validateParameters(parameters);
    PetscInt quadrature_order = 0;
    PetscReal tolerance = 0., newton_tolerance = 0.;
    CHKERR PetscOptionsGetInt(nullptr, nullptr, "-auxiliary_space_quadrature_order",
                              &quadrature_order, nullptr);
    CHKERR PetscOptionsGetReal(nullptr, nullptr, "-auxiliary_space_tolerance",
                               &tolerance, nullptr);
    CHKERR PetscOptionsGetReal(nullptr, nullptr,
                               "-auxiliary_space_newton_tolerance",
                               &newton_tolerance, nullptr);
    if (quadrature_order < 4 || !std::isfinite(tolerance) || tolerance <= 0. ||
        !std::isfinite(newton_tolerance) || newton_tolerance <= 0. ||
        newton_tolerance >= tolerance)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Space atom requires positive finite JSON tolerances and "
              "quadrature order at least four");
    CHKERR runSpaceAtom(parameters, quadrature_order, tolerance, newton_tolerance);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}