plastic_local_law_atom.cpp

Go to the documentation of this file.

/**
 * @file plastic_local_law_atom.cpp
 * @brief Pointwise checks of the ADOL-C two-field plastic local law.
 *
 * The law is the local density of the Chapter 1 positive-reconstruction
 * objective in doc/development/proximal_plasticity.tex (label
 * eq:tao-primal-objective). This atom needs no mesh: it records the tape,
 * then checks value, gradient and Hessian action against their closed forms
 * and against centred finite differences, using TestingInterface.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
#include <MatrixFunction.hpp>
#include <NeoHookeanLogarithmicMaterial.hpp>
#include <PlasticIncrementalOptimization.hpp>
#include <impl/PlasticIncrementalOptimizationInternal.hpp>
#include <impl/PlasticReducedHessianOperators.hpp>
using namespace EshelbianPlasticity;
using namespace EshelbianPlasticity::PlasticIncrementalOptimizationInternal;
 
static char help[] = "Pointwise checks of the two-field plastic local law.\n"
                     "  -test_nb 1,2,...  select the registered checks\n\n";
 
namespace {
 
constexpr int nbCoords = plasticLogarithmicStretchCoordinateSize;
constexpr int nbActive = plasticLocalLawActiveCount;
using Increment = PlasticLocalLaw::IncrementTensor;
using Gradient = PlasticLocalLaw::GradientTensor;
using Direction = PlasticLocalLaw::DirectionTensor;
using HessianAction = PlasticLocalLaw::HessianActionTensor;
using SymmetricTensor = FTensor::Tensor2_symmetric<double, SPACE_DIM>;
 
/** Exact first directional data for F_p = exp(H_p). */
struct PlasticExponentialDirectionalData {
  FTensor::Tensor2_symmetric<double, SPACE_DIM> plasticF;
  FTensor::Tensor2_symmetric<double, SPACE_DIM> inversePlasticF;
  FTensor::Tensor2_symmetric<double, SPACE_DIM> directionPlasticF;
  FTensor::Tensor2_symmetric<double, SPACE_DIM> directionInversePlasticF;
};
 
/** Evaluate exp(H), exp(-H), and their actions for derivative checks. */
MoFEMErrorCode evaluatePlasticExponentialDirectionalData(
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_plastic_h,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_direction,
    PlasticExponentialDirectionalData &data) {
  MoFEMFunctionBegin;
  static_assert(SPACE_DIM == 3,
                "Plastic exponential directional data is three-dimensional");
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
  FTensor::Tensor1<double, SPACE_DIM> t_eigen_values;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vectors;
  t_eigen_vectors(i, j) = t_plastic_h(i, j);
  if (computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values) !=
      MB_SUCCESS)
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Failed to diagonalise logarithmic plastic deformation for a "
            "directional exponential action");
  const int nb_unique = getUniqNb<SPACE_DIM>(t_eigen_values);
  if (nb_unique < SPACE_DIM)
    CHKERR sortEigenVals<SPACE_DIM>(t_eigen_values, t_eigen_vectors);
  const EigenMatrix::Fun<double> exp_fun = [](const double value) {
    return std::exp(value);
  };
  const EigenMatrix::Fun<double> inverse_exp_fun = [](const double value) {
    return std::exp(-value);
  };
  const EigenMatrix::Fun<double> derivative_inverse_exp_fun =
      [](const double value) { return -std::exp(-value); };
  const auto t_plastic_f =
      EigenMatrix::getMat(t_eigen_values, t_eigen_vectors, exp_fun);
  const auto t_inverse_plastic_f =
      EigenMatrix::getMat(t_eigen_values, t_eigen_vectors, inverse_exp_fun);
  const auto t_diff_exp = EigenMatrix::getDiffMat(
      t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, nb_unique);
  const auto t_diff_inverse_exp =
      EigenMatrix::getDiffMat(t_eigen_values, t_eigen_vectors, inverse_exp_fun,
                              derivative_inverse_exp_fun, nb_unique);
  data.plasticF(i, j) = t_plastic_f(i, j);
  data.inversePlasticF(i, j) = t_inverse_plastic_f(i, j);
  data.directionPlasticF(i, j) = t_diff_exp(i, j, k, l) * t_direction(k, l);
  data.directionInversePlasticF(i, j) =
      t_diff_inverse_exp(i, j, k, l) * t_direction(k, l);
  MoFEMFunctionReturn(0);
}
 
/** Contract D^2 exp(H^p) with two directions and the conjugate tensor
 * for the matrix-exponential second-derivative check. */
MoFEMErrorCode evaluatePlasticExponentialSecondContraction(
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_plastic_h,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_first_direction,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_second_direction,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_conjugate,
    double &value) {
  MoFEMFunctionBegin;
  static_assert(SPACE_DIM == 3,
                "Plastic exponential Hessian action is three-dimensional");
  FTENSOR_INDEXES(SPACE_DIM, i, j, k, l);
  FTensor::Tensor1<double, SPACE_DIM> t_eigen_values;
  FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_eigen_vectors;
  t_eigen_vectors(i, j) = t_plastic_h(i, j);
  if (computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values) !=
      MB_SUCCESS)
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "Failed to diagonalise logarithmic plastic deformation for a "
            "second directional exponential action");
  const int nb_unique = getUniqNb<SPACE_DIM>(t_eigen_values);
  if (nb_unique < SPACE_DIM)
    CHKERR sortEigenVals<SPACE_DIM>(t_eigen_values, t_eigen_vectors);
  const EigenMatrix::Fun<double> exp_fun = [](const double argument) {
    return std::exp(argument);
  };
  FTensor::Tensor2_symmetric<double, SPACE_DIM> t_mutable_conjugate;
  t_mutable_conjugate(i, j) = t_conjugate(i, j);
  const auto t_contracted_hessian = EigenMatrix::getDiffDiffMat(
      t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, exp_fun,
      t_mutable_conjugate, nb_unique);
  value = t_contracted_hessian(i, j, k, l) * t_first_direction(i, j) *
          t_second_direction(k, l);
  if (!std::isfinite(value))
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Contracted second directional plastic exponential is not finite");
  MoFEMFunctionReturn(0);
}
 
/** Closed-form gradient: g_H = f_r (2/3) h / ell_eps,
 * g_psi = e^psi [f_r + alpha (psi - psi^k)]. */
Gradient closedFormGradient(const Increment &t_h, double psi,
                            const PlasticLocalLawParameters &p) {
  const auto terms = evaluatePlasticLocalLawTerms(t_h, psi, p);
  Gradient t_g;
  for (int l = 0; l != nbCoords; ++l)
    t_g(l) = terms.resistanceForce * plasticEquivalentIncrementScaleSquared *
             t_h(l) / terms.regularisedMagnitude;
  t_g(nbCoords) =
      terms.slack *
      (terms.resistanceForce + p.proximalCoefficient * (psi - p.latentCentre));
  return t_g;
}
 
HessianAction closedFormHessianAction(const Increment &t_h, double psi,
                                      const Direction &t_direction,
                                      const PlasticLocalLawParameters &p) {
  const auto terms = evaluatePlasticLocalLawTerms(t_h, psi, p);
  std::array<double, nbCoords> normal;
  double delta_ell = 0;
  for (int l = 0; l != nbCoords; ++l) {
    normal[l] = plasticEquivalentIncrementScaleSquared * t_h(l) /
                terms.regularisedMagnitude;
    delta_ell += normal[l] * t_direction(l);
  }
  const double direction_latent = t_direction(nbCoords);
  const double delta_eta = delta_ell + terms.slack * direction_latent;
  const double delta_resistance = p.hardeningModulus * delta_eta;
  HessianAction t_action;
  for (int l = 0; l != nbCoords; ++l) {
    const double delta_normal =
        plasticEquivalentIncrementScaleSquared * t_direction(l) /
            terms.regularisedMagnitude -
        normal[l] * delta_ell / terms.regularisedMagnitude;
    t_action(l) =
        delta_resistance * normal[l] + terms.resistanceForce * delta_normal;
  }
  t_action(nbCoords) =
      terms.slack * delta_resistance +
      terms.slack *
          (terms.resistanceForce +
           p.proximalCoefficient * (psi - p.latentCentre + 1.)) *
          direction_latent;
  return t_action;
}
 
double maxAbs(const Gradient &t_a) {
  double m = 0;
  for (int a = 0; a != nbActive; ++a)
    m = std::max(m, std::abs(t_a(a)));
  return m;
}
 
double maxAbsDifference(const Gradient &t_a, const Gradient &t_b) {
  double m = 0;
  for (int a = 0; a != nbActive; ++a)
    m = std::max(m, std::abs(t_a(a) - t_b(a)));
  return m;
}
 
struct Point {
  Increment t_h;
  double psi;
};
 
void perturb(const Point &base, const Gradient &t_direction, double step,
             Point &out) {
  for (int l = 0; l != nbCoords; ++l)
    out.t_h(l) = base.t_h(l) + step * t_direction(l);
  out.psi = base.psi + step * t_direction(nbCoords);
}
 
double maxAbsDifference(const SymmetricTensor &t_a,
                        const SymmetricTensor &t_b) {
  double result = 0;
  for (int i = 0; i != SPACE_DIM; ++i)
    for (int j = 0; j != SPACE_DIM; ++j)
      result = std::max(result, std::abs(t_a(i, j) - t_b(i, j)));
  return result;
}
 
double maxAbs(const SymmetricTensor &t_a) {
  double result = 0;
  for (int i = 0; i != SPACE_DIM; ++i)
    for (int j = 0; j != SPACE_DIM; ++j)
      result = std::max(result, std::abs(t_a(i, j)));
  return result;
}
 
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
  try {
    moab::Core moab_core;
    moab::Interface &moab = moab_core;
    MoFEM::Core mofem_core(moab);
    MoFEM::Interface &m_field = mofem_core;
    auto testing = m_field.getInterface<TestingInterface>();
 
    auto law = createPlasticLocalLaw();
    CHKERR law->recordTape();
 
    // Reference parameters: hardening on, plastic slack open.
    const PlasticLocalLawParameters hardening{/*yieldStress*/ 0.1,
                                              /*hardeningModulus*/ 0.02,
                                              /*committedKappa*/ 0.03,
                                              /*latentCentre*/ std::log(0.01),
                                              /*proximalCoefficient*/ 0.1,
                                              /*regularisationEpsilon*/ 1e-4};
    PlasticLocalLawParameters perfect = hardening;
    perfect.hardeningModulus = 0;
    const std::array<const PlasticLocalLawParameters *, 2> parameter_sets{
        &hardening, &perfect};
 
    const Point active{Increment{0.012, -0.004, 0.007, 0.002, -0.009},
                       std::log(0.003)};
    const Point elastic{Increment{0., 0., 0., 0., 0.}, std::log(0.0005)};
 
    const std::array<Gradient, 3> directions{
        Gradient{0.3, -0.7, 0.2, 0.5, -0.4, 0.},
        Gradient{0., 0., 0., 0., 0., 1.},
        Gradient{-0.2, 0.4, 0.6, -0.3, 0.1, 0.8}};
 
    // 1. Tape value equals the shared formula, and its terms are consistent.
    CHKERR testing->addTestCommandLine("", 1, "tape value", [&]() {
      MoFEMFunctionBegin;
      for (const auto *p : parameter_sets) {
        CHKERR law->setParameters(*p);
        for (const auto *pt : {&active, &elastic}) {
          PlasticLocalLawTerms terms;
          CHKERR law->evaluateValue(pt->t_h, pt->psi, terms);
          const auto direct =
              evaluatePlasticLocalLawTerms(pt->t_h, pt->psi, *p);
          CHKERR testing->checkValue<TestType::SCALAR>(
              "density tape vs direct", terms.density, direct.density,
              {1e-15, 1e-13});
          CHKERR testing->checkValue<TestType::SCALAR>(
              "density = resistance + proximal", terms.density,
              terms.resistance + terms.proximal, {1e-15, 1e-13});
          CHKERR testing->checkError<TestType::SCALAR>(
              "proximal Bregman term non-negative",
              std::max(0., -terms.proximal), 1., {0., 0.});
        }
      }
      MoFEMFunctionReturn(0);
    });
 
    // 2. Reconstruction admissibility: eta - ell_eps = e^psi > 0, ell_eps >=
    //    ell, and ell_eps(0) = eps.
    CHKERR testing->addTestCommandLine("", 2, "reconstruction", [&]() {
      MoFEMFunctionBegin;
      CHKERR law->setParameters(hardening);
      PlasticLocalLawTerms terms;
      CHKERR law->evaluateValue(active.t_h, active.psi, terms);
      CHKERR testing->checkValue<TestType::SCALAR>(
          "slack equals exp(psi)",
          terms.equivalentIncrement - terms.regularisedMagnitude,
          std::exp(active.psi), {1e-15, 1e-13});
      CHKERR testing->checkError<TestType::SCALAR>(
          "slack positive", terms.slack > 0 ? 0. : 1., 1., {0., 0.});
      const double exact_norm =
          plasticEquivalentIncrementScale * plasticCoordinateNorm(active.t_h);
      CHKERR testing->checkError<TestType::SCALAR>(
          "ell_eps >= ell",
          std::max(0., exact_norm - terms.regularisedMagnitude), exact_norm,
          {0., 0.});
      CHKERR testing->checkValue<TestType::SCALAR>(
          "ell_eps^2 = ell^2 + eps^2",
          terms.regularisedMagnitude * terms.regularisedMagnitude,
          exact_norm * exact_norm +
              hardening.regularisationEpsilon * hardening.regularisationEpsilon,
          {1e-15, 1e-13});
      CHKERR law->evaluateValue(elastic.t_h, elastic.psi, terms);
      CHKERR testing->checkValue<TestType::SCALAR>(
          "ell_eps at origin equals eps", terms.regularisedMagnitude,
          hardening.regularisationEpsilon, {1e-15, 1e-13});
      CHKERR testing->checkValue<TestType::SCALAR>(
          "f_r at origin", terms.resistanceForce,
          hardening.yieldStress +
              hardening.hardeningModulus *
                  (hardening.committedKappa + hardening.regularisationEpsilon +
                   std::exp(elastic.psi)),
          {1e-15, 1e-13});
      MoFEMFunctionReturn(0);
    });
 
    // 3. Tape gradient equals the closed form of Eq. (tao-primal-gradient)
    //    and centred finite differences of the tape value.
    CHKERR testing->addTestCommandLine("", 3, "gradient", [&]() {
      MoFEMFunctionBegin;
      for (const auto *p : parameter_sets) {
        CHKERR law->setParameters(*p);
        for (const auto *pt : {&active, &elastic}) {
          Gradient t_g;
          CHKERR law->evaluateGradient(pt->t_h, pt->psi, t_g);
          const auto t_closed = closedFormGradient(pt->t_h, pt->psi, *p);
          CHKERR testing->checkError<TestType::DERIVATIVE>(
              "gradient tape vs closed form", maxAbsDifference(t_g, t_closed),
              std::max(maxAbs(t_g), maxAbs(t_closed)), {1e-16, 1e-12});
          for (const auto &t_v : directions) {
            double assembled = 0;
            for (int a = 0; a != nbActive; ++a)
              assembled += t_g(a) * t_v(a);
            for (const double step : {1e-5, 1e-6}) {
              Point plus, minus;
              perturb(*pt, t_v, step, plus);
              perturb(*pt, t_v, -step, minus);
              PlasticLocalLawTerms terms_plus, terms_minus;
              CHKERR law->evaluateValue(plus.t_h, plus.psi, terms_plus);
              CHKERR law->evaluateValue(minus.t_h, minus.psi, terms_minus);
              const double fd =
                  (terms_plus.density - terms_minus.density) / (2 * step);
              // Truncation is O(step^2 |phi'''|); phi''' ~ f_r / ell^2.
              const double tolerance = step == 1e-5 ? 1e-5 : 1e-7;
              CHKERR testing->checkError<TestType::DERIVATIVE>(
                  "gradient centred difference", std::abs(fd - assembled),
                  std::max(std::abs(fd), std::abs(assembled)),
                  {1e-12, tolerance});
            }
          }
        }
      }
      MoFEMFunctionReturn(0);
    });
 
    // 4. Perfect plasticity: the latent stationarity point is explicit,
    //    psi = psi^k - sigma_y / alpha, and the tensor gradient does not
    //    depend on psi (Eq. perfect-plasticity-latent-solution).
    CHKERR testing->addTestCommandLine("", 4, "perfect plasticity", [&]() {
      MoFEMFunctionBegin;
      CHKERR law->setParameters(perfect);
      const double psi_star = perfect.latentCentre -
                              perfect.yieldStress / perfect.proximalCoefficient;
      Gradient t_g_star, t_g_other;
      CHKERR law->evaluateGradient(active.t_h, psi_star, t_g_star);
      CHKERR law->evaluateGradient(active.t_h, active.psi, t_g_other);
      CHKERR testing->checkError<TestType::DERIVATIVE>(
          "latent gradient vanishes at the explicit root",
          std::abs(t_g_star(nbCoords)),
          std::exp(psi_star) * perfect.yieldStress, {1e-16, 1e-12});
      double tensor_difference = 0, tensor_scale = 0;
      for (int l = 0; l != nbCoords; ++l) {
        tensor_difference =
            std::max(tensor_difference, std::abs(t_g_star(l) - t_g_other(l)));
        tensor_scale = std::max(tensor_scale, std::abs(t_g_other(l)));
      }
      CHKERR testing->checkError<TestType::DERIVATIVE>(
          "tensor gradient independent of psi", tensor_difference, tensor_scale,
          {1e-16, 1e-14});
      // The latent gradient has the sign of sigma_y + alpha (psi - psi^k).
      Gradient t_g_above;
      CHKERR law->evaluateGradient(active.t_h, psi_star + 1., t_g_above);
      CHKERR testing->checkError<TestType::SCALAR>(
          "latent gradient positive above the root",
          t_g_above(nbCoords) > 0 ? 0. : 1., 1., {0., 0.});
      MoFEMFunctionReturn(0);
    });
 
    // 5. Hardening couples the fields only through A_h e^psi: closing the
    //    slack changes the tensor gradient by A_h (e^psi_1 - e^psi_2) G.
    CHKERR testing->addTestCommandLine("", 5, "hardening coupling", [&]() {
      MoFEMFunctionBegin;
      CHKERR law->setParameters(hardening);
      const double psi_open = active.psi;
      const double psi_closed = std::log(1e-7);
      Gradient t_g_open, t_g_closed;
      CHKERR law->evaluateGradient(active.t_h, psi_open, t_g_open);
      CHKERR law->evaluateGradient(active.t_h, psi_closed, t_g_closed);
      const auto terms =
          evaluatePlasticLocalLawTerms(active.t_h, psi_open, hardening);
      const double factor = hardening.hardeningModulus *
                            (std::exp(psi_open) - std::exp(psi_closed));
      double error = 0, scale = 0;
      for (int l = 0; l != nbCoords; ++l) {
        const double predicted = factor *
                                 plasticEquivalentIncrementScaleSquared *
                                 active.t_h(l) / terms.regularisedMagnitude;
        error = std::max(error,
                         std::abs((t_g_open(l) - t_g_closed(l)) - predicted));
        scale = std::max(scale, std::abs(t_g_open(l)));
      }
      CHKERR testing->checkError<TestType::DERIVATIVE>(
          "tensor gradient shift equals A_h delta(e^psi) G", error, scale,
          {1e-16, 1e-12});
      MoFEMFunctionReturn(0);
    });
 
    // 6. Invalid parameters are rejected: eps = 0 and non-positive alpha.
    CHKERR testing->addTestCommandLine("", 6, "parameter validation", [&]() {
      MoFEMFunctionBegin;
      auto rejects = [&](PlasticLocalLawParameters bad) {
        CHKERR PetscPushErrorHandler(PetscIgnoreErrorHandler, nullptr);
        const auto ierr = law->setParameters(bad);
        CHKERR PetscPopErrorHandler();
        return ierr != 0 ? 0. : 1.;
      };
      auto bad_eps = hardening;
      bad_eps.regularisationEpsilon = 0;
      CHKERR testing->checkError<TestType::SCALAR>(
          "epsilon zero rejected", rejects(bad_eps), 1., {0., 0.});
      auto bad_alpha = hardening;
      bad_alpha.proximalCoefficient = 0;
      CHKERR testing->checkError<TestType::SCALAR>(
          "zero proximal coefficient rejected", rejects(bad_alpha), 1.,
          {0., 0.});
      auto bad_hardening = hardening;
      bad_hardening.hardeningModulus = -1;
      CHKERR testing->checkError<TestType::SCALAR>(
          "negative hardening rejected", rejects(bad_hardening), 1., {0., 0.});
      // The valid parameters still apply after a rejected update.
      CHKERR law->setParameters(hardening);
      CHKERR testing->checkValue<TestType::SCALAR>(
          "valid parameters restored",
          law->getParameters().regularisationEpsilon,
          hardening.regularisationEpsilon, {0., 0.});
      MoFEMFunctionReturn(0);
    });
 
    // 7. The ADOL-C Hessian-vector product equals the analytical directional
    //    derivative and a centred difference of the pointwise gradient.
    CHKERR testing->addTestCommandLine("", 7, "hessian action", [&]() {
      MoFEMFunctionBegin;
      for (const auto *p : parameter_sets) {
        CHKERR law->setParameters(*p);
        for (const auto *pt : {&active, &elastic}) {
          for (const auto &t_v : directions) {
            HessianAction t_action;
            CHKERR law->evaluateHessianAction(pt->t_h, pt->psi, t_v, t_action);
            const auto t_closed =
                closedFormHessianAction(pt->t_h, pt->psi, t_v, *p);
            CHKERR testing->checkError<TestType::DERIVATIVE>(
                "Hessian action tape vs closed form",
                maxAbsDifference(t_action, t_closed),
                std::max(maxAbs(t_action), maxAbs(t_closed)), {1e-12, 1e-11});
 
            // At h = 0 the regularised norm has curvature O(1 / eps), so use
            // a smaller perturbation to keep the centred-difference truncation
            // error below the derivative tolerance.
            const double step = pt == &elastic ? 1e-8 : 1e-7;
            Point plus, minus;
            perturb(*pt, t_v, step, plus);
            perturb(*pt, t_v, -step, minus);
            Gradient t_gradient_plus, t_gradient_minus;
            CHKERR law->evaluateGradient(plus.t_h, plus.psi, t_gradient_plus);
            CHKERR law->evaluateGradient(minus.t_h, minus.psi,
                                         t_gradient_minus);
            HessianAction t_difference;
            for (int a = 0; a != nbActive; ++a)
              t_difference(a) =
                  (t_gradient_plus(a) - t_gradient_minus(a)) / (2. * step);
            CHKERR testing->checkError<TestType::DERIVATIVE>(
                "Hessian action centred gradient difference",
                maxAbsDifference(t_action, t_difference),
                std::max(maxAbs(t_action), maxAbs(t_difference)), {1e-8, 2e-7});
          }
        }
      }
 
      CHKERR law->setParameters(perfect);
      HessianAction t_perfect_cross;
      CHKERR law->evaluateHessianAction(active.t_h, active.psi, directions[1],
                                        t_perfect_cross);
      double tensor_cross_norm = 0;
      for (int l = 0; l != nbCoords; ++l)
        tensor_cross_norm =
            std::max(tensor_cross_norm, std::abs(t_perfect_cross(l)));
      CHKERR testing->checkError<TestType::DERIVATIVE>(
          "perfect-plastic tensor-latent Hessian coupling vanishes",
          tensor_cross_norm, std::abs(t_perfect_cross(nbCoords)),
          {1e-16, 1e-13});
 
      CHKERR law->setParameters(hardening);
      HessianAction t_origin_action;
      CHKERR law->evaluateHessianAction(elastic.t_h, elastic.psi, directions[0],
                                        t_origin_action);
      const auto origin_terms =
          evaluatePlasticLocalLawTerms(elastic.t_h, elastic.psi, hardening);
      double origin_error = 0;
      double origin_scale = 0;
      for (int l = 0; l != nbCoords; ++l) {
        const double expected = origin_terms.resistanceForce *
                                plasticEquivalentIncrementScaleSquared *
                                directions[0](l) /
                                hardening.regularisationEpsilon;
        origin_error =
            std::max(origin_error, std::abs(t_origin_action(l) - expected));
        origin_scale = std::max(origin_scale, std::abs(expected));
      }
      CHKERR testing->checkError<TestType::DERIVATIVE>(
          "regularised tensor curvature at the origin", origin_error,
          origin_scale, {1e-10, 1e-12});
      MoFEMFunctionReturn(0);
    });
 
    // 8. Check first and contracted second matrix-exponential derivatives
    //    against centred finite differences.
    CHKERR testing->addTestCommandLine(
        "", 8, "plastic exponential directional actions", [&]() {
          MoFEMFunctionBegin;
          Increment t_u_coordinates;
          Increment t_v_coordinates;
          for (int l = 0; l != nbCoords; ++l) {
            t_u_coordinates(l) = directions[0](l);
            t_v_coordinates(l) = directions[2](l);
          }
          const auto t_h =
              plasticLogarithmicStretchTensorFromCoordinates(active.t_h);
          const auto t_u =
              plasticLogarithmicStretchTensorFromCoordinates(t_u_coordinates);
          const auto t_v =
              plasticLogarithmicStretchTensorFromCoordinates(t_v_coordinates);
          SymmetricTensor t_s{0.7, -0.2, 0.1, 0.5, 0.3, -0.4};
          PlasticExponentialDirectionalData data;
          CHKERR evaluatePlasticExponentialDirectionalData(t_h, t_u, data);
 
          constexpr double step = 2e-6;
          SymmetricTensor t_plus_h;
          SymmetricTensor t_minus_h;
          for (int i = 0; i != SPACE_DIM; ++i)
            for (int j = 0; j <= i; ++j) {
              t_plus_h(i, j) = t_h(i, j) + step * t_u(i, j);
              t_minus_h(i, j) = t_h(i, j) - step * t_u(i, j);
            }
          PlasticExponentialDirectionalData plus_data;
          PlasticExponentialDirectionalData minus_data;
          CHKERR evaluatePlasticExponentialDirectionalData(t_plus_h, t_v,
                                                           plus_data);
          CHKERR evaluatePlasticExponentialDirectionalData(t_minus_h, t_v,
                                                           minus_data);
 
          SymmetricTensor t_first_difference_f;
          SymmetricTensor t_first_difference_inverse;
          for (int i = 0; i != SPACE_DIM; ++i)
            for (int j = 0; j <= i; ++j) {
              t_first_difference_f(i, j) =
                  (plus_data.plasticF(i, j) - minus_data.plasticF(i, j)) /
                  (2. * step);
              t_first_difference_inverse(i, j) =
                  (plus_data.inversePlasticF(i, j) -
                   minus_data.inversePlasticF(i, j)) /
                  (2. * step);
            }
          CHKERR testing->checkError<TestType::DERIVATIVE>(
              "D exp(H)[U]",
              maxAbsDifference(data.directionPlasticF, t_first_difference_f),
              std::max(maxAbs(data.directionPlasticF),
                       maxAbs(t_first_difference_f)),
              {1e-10, 2e-8});
          CHKERR testing->checkError<TestType::DERIVATIVE>(
              "D exp(-H)[U]",
              maxAbsDifference(data.directionInversePlasticF,
                               t_first_difference_inverse),
              std::max(maxAbs(data.directionInversePlasticF),
                       maxAbs(t_first_difference_inverse)),
              {1e-10, 2e-8});
 
          double exact_second_contraction = 0;
          CHKERR evaluatePlasticExponentialSecondContraction(
              t_h, t_u, t_v, t_s, exact_second_contraction);
          double difference_second_contraction = 0;
          for (int i = 0; i != SPACE_DIM; ++i)
            for (int j = 0; j != SPACE_DIM; ++j)
              difference_second_contraction +=
                  t_s(i, j) *
                  (plus_data.directionPlasticF(i, j) -
                   minus_data.directionPlasticF(i, j)) /
                  (2. * step);
          CHKERR testing->checkError<TestType::DERIVATIVE>(
              "S:D2 exp(H)[U,V]",
              std::abs(exact_second_contraction -
                       difference_second_contraction),
              std::max(std::abs(exact_second_contraction),
                       std::abs(difference_second_contraction)),
              {1e-10, 3e-8});
          MoFEMFunctionReturn(0);
        });
 
    CHKERR testing->addTestCommandLine(
        "", 9, "local elastic logarithmic sensitivity", [&]() {
          MoFEMFunctionBegin;
          FTENSOR_INDEXES(3, i, j, k);
          FTensor::Index<'a', 5> a;
          FTensor::Index<'b', 5> b;
          const auto function = [](const SymmetricTensor &tensor,
                                   EigenMatrix::Fun<double> f) {
            FTENSOR_INDEXES(3, m, n);
            FTensor::Tensor2<double, 3, 3> vectors;
            FTensor::Tensor1<double, 3> values;
            vectors(m, n) = tensor(m, n);
            CHK_THROW_MESSAGE(computeEigenValuesSymmetric(vectors, values),
                              "Diagonalise local sensitivity test tensor");
            return EigenMatrix::getMat(values, vectors, f);
          };
          const EigenMatrix::Fun<double> exp_f = [](double x) {
            return std::exp(x);
          };
          const EigenMatrix::Fun<double> inv_f = [](double x) {
            return std::exp(-x);
          };
          const EigenMatrix::Fun<double> log_f = [](double x) {
            return .5 * std::log(x);
          };
          const std::array<SymmetricTensor, 4> elastic_logs{
              {{0., 0., 0., 0., 0., 0.},
               {.12, 0., 0., .12, 0., -.08},
               {.12, 0., 0., .12 + 1e-12, 0., -.08},
               {.18, .07, -.04, -.09, .03, .02}}};
          const std::array<SymmetricTensor, 3> plastic_logs{
              {{0., 0., 0., 0., 0., 0.},
               {.1, 0., 0., .1, 0., -.2},
               {.21, -.09, .04, -.13, .06, -.08}}};
          double worst_error = 0.;
          for (const auto &he : elastic_logs) {
            const auto ue = function(he, exp_f);
            for (const auto &hp : plastic_logs) {
              FTensor::Tensor2<double, 5, 5> map;
              CHKERR evaluateLocalElasticLogarithmicSensitivity(hp, ue, map);
              const auto fp = function(hp, exp_f);
              FTensor::Tensor2<double, 3, 3> anchor;
              anchor(i, j) = ue(i, k) * fp(k, j);
              for (const auto &source : {directions[0], directions[2]}) {
                Increment v, exact;
                for (int n = 0; n != 5; ++n)
                  v(n) = source(n);
                const auto tv =
                    plasticLogarithmicStretchTensorFromCoordinates(v);
                exact(a) = map(a, b) * v(b);
                for (const double step : {2e-5, 2e-6}) {
                  const auto evaluate = [&](double shift) {
                    SymmetricTensor trial_hp, c;
                    trial_hp(i, j) = hp(i, j) + shift * tv(i, j);
                    const auto inverse = function(trial_hp, inv_f);
                    FTensor::Tensor2<double, 3, 3> fe, metric;
                    fe(i, j) = anchor(i, k) * inverse(k, j);
                    metric(i, j) = fe(k, i) * fe(k, j);
                    c(i, j) = (metric(i, j) || metric(j, i)) / 2.;
                    return function(c, log_f);
                  };
                  const auto plus = evaluate(step), minus = evaluate(-step);
                  SymmetricTensor derivative;
                  derivative(i, j) = (plus(i, j) - minus(i, j)) / (2. * step);
                  const auto reference =
                      plasticLogarithmicStretchCoordinatesFromTensor(
                          derivative);
                  double error = 0.;
                  for (int n = 0; n != 5; ++n)
                    error = std::max(error, std::abs(exact(n) - reference(n)));
                  worst_error = std::max(worst_error, error);
                  CHKERR testing->checkError<TestType::DERIVATIVE>(
                      "Direct Ue sensitivity vs fixed-F differences", error,
                      plasticCoordinateNorm(reference), {1e-9, 2e-7});
                  CHKERR testing->checkError<TestType::DERIVATIVE>(
                      "Isochoric elastic logarithmic sensitivity",
                      std::abs(derivative(i, i)), 1., {2e-8, 2e-8});
                }
              }
            }
            using Material = NeoHookeanLogarithmicMaterial;
            const Material::Parameters parameters{1.7, 8.5};
            const auto coordinates = Material::getCoordinates(he);
            Material::Tangent tangent;
            double energy = 0.;
            CHKERR Material::evaluateDeviator(parameters, coordinates, energy,
                                              nullptr, &tangent);
            for (int n = 0; n != 5; ++n) {
              constexpr double step = 2e-6;
              Increment unit, plus, minus, gp, gm, action;
              unit(a) = 0.;
              unit(n) = 1.;
              plus(a) = coordinates(a) + step * unit(a);
              minus(a) = coordinates(a) - step * unit(a);
              CHKERR Material::evaluateDeviator(parameters, plus, energy, &gp);
              CHKERR Material::evaluateDeviator(parameters, minus, energy, &gm);
              action(a) = tangent(a, b) * unit(b);
              double error = 0.;
              for (int col = 0; col != 5; ++col)
                error =
                    std::max(error, std::abs(action(col) - (gp(col) - gm(col)) /
                                                               (2. * step)));
              CHKERR testing->checkError<TestType::DERIVATIVE>(
                  "Pure logarithmic material tangent", error,
                  plasticCoordinateNorm(action), {2e-8, 2e-7});
            }
          }
          SymmetricTensor zero{0., 0., 0., 0., 0., 0.};
          SymmetricTensor identity{1., 0., 0., 1., 0., 1.};
          FTensor::Tensor2<double, 5, 5> map;
          CHKERR evaluateLocalElasticLogarithmicSensitivity(zero, identity,
                                                            map);
          for (int row = 0; row != 5; ++row)
            for (int col = 0; col != 5; ++col)
              if (std::abs(map(row, col) + (row == col ? 1. : 0.)) > 1e-12)
                SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                        "Identity elastic sensitivity is not minus identity");
          NeoHookeanLogarithmicMaterial::Tangent tangent;
          double energy = 0.;
          CHKERR NeoHookeanLogarithmicMaterial::evaluateDeviator(
              {1.7, 8.5}, NeoHookeanLogarithmicMaterial::getCoordinates(zero),
              energy, nullptr, &tangent);
          FTensor::Index<'d', 5> d;
          FTensor::Tensor2<double, 5, 5> tangent_map, estimate;
          tangent_map(a, b) = tangent(a, d) * map(d, b);
          estimate(a, b) = map(d, a) * tangent_map(d, b);
          for (int row = 0; row != 5; ++row)
            for (int col = 0; col != 5; ++col)
              if (std::abs(estimate(row, col) - (row == col ? 6.8 : 0.)) >
                  1e-11)
                SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                        "Identity local material estimate is not 2*mu*I");
          MOFEM_LOG("WORLD", Sev::inform)
              << "Local stretch sensitivity maximum FD error " << worst_error;
          MoFEMFunctionReturn(0);
        });
 
    CHKERR testing->runTests();
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}