PlasticIncrementalConstitutive.cpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalConstitutive.cpp
 * @brief Local constitutive maps for incremental optimization
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <PlasticIncrementalOptimization.hpp>
#include <impl/PlasticIncrementalOptimizationInternal.hpp>
#include <sstream>
 
namespace EshelbianPlasticity {
 
using namespace PlasticIncrementalOptimizationInternal;
 
namespace PlasticIncrementalOptimizationInternal {
 
namespace {
 
using PlasticIncrementalEle = VolumeElementForcesAndSourcesCore;
using PlasticIncrementalOp = PlasticIncrementalEle::UserDataOperator;
 
void addTetrahedronOperator(const boost::shared_ptr<PlasticIncrementalEle> &fe,
                            ForcesAndSourcesCore::UserDataOperator *op) {
  op->doEntities.fill(false);
  op->doEntities[MBTET] = true;
  fe->getOpPtrVector().push_back(op);
}
 
double getCellMeasure(PlasticIncrementalOp &op) {
  auto t_weight = op.getFTensor0IntegrationWeight();
  double weight_sum = 0;
  for (int gg = 0; gg != op.getGaussPts().size2(); ++gg) {
    weight_sum += t_weight;
    ++t_weight;
  }
  return op.getMeasure() * weight_sum;
}
 
boost::shared_ptr<PlasticIncrementalEle> createPlasticIncrementalElement(
    const PlasticIncrementalOptimizationProblem &problem) {
  auto fe = boost::make_shared<PlasticIncrementalEle>(
      *getInterfacePtr(problem.getControlDM()));
  // The parent volume FE also carries the USER_BASE bubble field. Although
  // these algebraic P0 operators do not use basis values directly, the FE
  // lifecycle prepares every registered approximation base.
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  // Keep the constitutive P0 measure identical to the mechanical volume
  // integration for both affine and high-order geometry.
  fe->getRuleHook = [](int, int, int order_data) {
    return 2 * (order_data + 1);
  };
  CHK_THROW_MESSAGE(
      problem.addConstitutiveGeometryOperators(fe->getOpPtrVector()),
      "Add plastic constitutive geometry operators");
  return fe;
}
 
struct PlasticValidationMetrics {
  double minWork = std::numeric_limits<double>::max();
  double maxYieldFunction = 0;
  double maxYieldExcess = 0;
  double minConstraint = std::numeric_limits<double>::max();
  double minMultiplier = std::numeric_limits<double>::max();
  double maxActiveConstraint = 0;
  double maxEtaStationarityResidual = 0;
  double maxComplementarityResidual = 0;
  double maxKktResidual = 0;
  double maxEquivalentIncrement = 0;
  int nonfinite = 0;
  int invalidMeasure = 0;
};
 
} // namespace
 
MoFEMErrorCode evaluatePlasticIncrementalResistanceValue(
    const PlasticIncrementalOptimizationProblem &problem, Vec control,
    double &value) {
  MoFEMFunctionBegin;
  value = 0;
  const auto &kappa_field = problem.getPlasticKappaFieldName();
 
  auto fe = createPlasticIncrementalElement(problem);
  DM control_dm = problem.getControlDM();
  const auto control_data_dm = SmartPetscObj<DM>(control_dm, true);
  const auto control_vec = SmartPetscObj<Vec>(control, true);
  auto kappa_increment_ptr = boost::make_shared<VectorDouble>();
  auto committed_kappa_ptr = boost::make_shared<VectorDouble>();
  auto local_value_ptr = boost::make_shared<double>(0);
 
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 control_data_dm, PlasticIncrementalOp::OPCOL,
                                 kappa_field, kappa_increment_ptr, control_vec,
                                 MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 kappa_field, committed_kappa_ptr, MBTET));
 
  auto *resistance_op =
      new PlasticIncrementalOp(kappa_field, PlasticIncrementalOp::OPROW);
  resistance_op->doWorkRhsHook =
      [&, local_value_ptr](
          DataOperator *base_op_ptr, int, EntityType,
          EntitiesFieldData::EntData &row_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    const auto t_kappa_increment = getFTensor0FromVec(*kappa_increment_ptr);
    const auto t_committed_kappa = getFTensor0FromVec(*committed_kappa_ptr);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    // Irreversible work plus the linear-hardening Helmholtz-energy increment.
    // Both P0 values are element-constant, so the first Gauss-point values
    // define this cell contribution.
    const double delta_kappa = t_kappa_increment;
    const double block_value =
        getCellMeasure(*op_ptr) *
        ((problem.getInitialYieldStress() +
          problem.getIsotropicHardeningModulus() * t_committed_kappa) *
             delta_kappa +
         0.5 * problem.getIsotropicHardeningModulus() * delta_kappa *
             delta_kappa);
    if (!std::isfinite(block_value))
      *local_value_ptr = std::numeric_limits<double>::infinity();
    else
      *local_value_ptr += block_value;
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, resistance_op);
 
  CHKERR DMoFEMLoopFiniteElements(control_dm, problem.getVolumeElementName(),
                                  fe);
  CHKERR MPI_Allreduce(local_value_ptr.get(), &value, 1, MPI_DOUBLE, MPI_SUM,
                       PetscObjectComm(reinterpret_cast<PetscObject>(control)));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode assembleIncrementalObjectiveGradient(
    PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec smooth_gradient, Vec objective_gradient) {
  MoFEMFunctionBegin;
  CHKERR VecCopy(smooth_gradient, objective_gradient);
  const auto &kappa_field = problem.getPlasticKappaFieldName();
 
  auto fe = createPlasticIncrementalElement(problem);
  fe->f = objective_gradient;
 
  DM control_dm = problem.getControlDM();
  const auto control_data_dm = SmartPetscObj<DM>(control_dm, true);
  const auto control_vec = SmartPetscObj<Vec>(control, true);
  auto kappa_increment_ptr = boost::make_shared<VectorDouble>();
  auto committed_kappa_ptr = boost::make_shared<VectorDouble>();
 
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 control_data_dm, PlasticIncrementalOp::OPCOL,
                                 kappa_field, kappa_increment_ptr, control_vec,
                                 MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 kappa_field, committed_kappa_ptr, MBTET));
 
  auto *kappa_op =
      new PlasticIncrementalOp(kappa_field, PlasticIncrementalOp::OPROW);
  kappa_op->doWorkRhsHook =
      [&](DataOperator *base_op_ptr, int, EntityType,
          EntitiesFieldData::EntData &row_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    const auto t_kappa_increment = getFTensor0FromVec(*kappa_increment_ptr);
    const auto t_committed_kappa = getFTensor0FromVec(*committed_kappa_ptr);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    // Second block of Eq. (1.70), label
    // eq:constrained-objective-gradient. These P0 producer values are
    // element-constant, so the first Gauss-point values define the row.
    const double local_gradient =
        getCellMeasure(*op_ptr) * (problem.getInitialYieldStress() +
                                   problem.getIsotropicHardeningModulus() *
                                       (t_committed_kappa + t_kappa_increment));
    CHKERR VecSetValues<AssemblyTypeSelector<PETSC>>(
        op_ptr->getKSPf(), row_data, &local_gradient, ADD_VALUES);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, kappa_op);
 
  CHKERR DMoFEMLoopFiniteElements(control_dm, problem.getVolumeElementName(),
                                  fe);
  CHKERR VecAssemblyBegin(objective_gradient);
  CHKERR VecAssemblyEnd(objective_gradient);
  CHKERR VecGhostUpdateBegin(objective_gradient, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(objective_gradient, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode initialisePlasticInequalityMultipliers(
    const PlasticIncrementalOptimizationProblem &problem, DM constraint_dm,
    Vec multipliers) {
  MoFEMFunctionBegin;
  if (!constraint_dm || !multipliers)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Plastic multiplier initialization requires a constraint DM and "
            "multiplier vector");
 
  const auto &kappa_field = problem.getPlasticKappaFieldName();
  auto fe = createPlasticIncrementalElement(problem);
  fe->f = multipliers;
  auto committed_kappa_ptr = boost::make_shared<VectorDouble>();
 
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 kappa_field, committed_kappa_ptr, MBTET));
  auto *multiplier_op =
      new PlasticIncrementalOp(kappa_field, PlasticIncrementalOp::OPROW);
  multiplier_op->doWorkRhsHook =
      [&](DataOperator *base_op_ptr, int, EntityType,
          EntitiesFieldData::EntData &row_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    const double committed_kappa =
        getFTensor0FromVec(*committed_kappa_ptr);
    if (!std::isfinite(committed_kappa) || committed_kappa < 0)
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Committed plastic kappa must be finite and non-negative; got "
              "%g",
              committed_kappa);
    const double current_yield =
        problem.getInitialYieldStress() +
        problem.getIsotropicHardeningModulus() * committed_kappa;
    if (!std::isfinite(current_yield) || !(current_yield > 0))
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "Current plastic yield stress must be finite and positive; got "
              "%g",
              current_yield);
    const double constraint_scale = std::sqrt(getCellMeasure(*op_ptr));
    // J_eta=sqrt(V), so eta stationarity requires lambda_TAO=sqrt(V)*sigma_y.
    const double local_multiplier = constraint_scale * current_yield;
    CHKERR VecSetValues<AssemblyTypeSelector<PETSC>>(
        op_ptr->getKSPf(), row_data, &local_multiplier, INSERT_VALUES);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, multiplier_op);
 
  CHKERR DMoFEMLoopFiniteElements(constraint_dm,
                                  problem.getVolumeElementName(), fe);
  CHKERR VecAssemblyBegin(multipliers);
  CHKERR VecAssemblyEnd(multipliers);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode evaluatePlasticInequalityConstraints(
    const PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec constraints) {
  MoFEMFunctionBegin;
  const auto &plastic_field = problem.getPlasticFlowFieldName();
  const auto &kappa_field = problem.getPlasticKappaFieldName();
  DM constraint_dm = nullptr;
  CHKERR VecGetDM(constraints, &constraint_dm);
  if (!constraint_dm)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Plastic inequality vector has no constraint DM");
 
  auto fe = createPlasticIncrementalElement(problem);
  fe->f = constraints;
  const auto control_dm = SmartPetscObj<DM>(problem.getControlDM(), true);
  const auto control_vec = SmartPetscObj<Vec>(control, true);
  auto plastic_increment_ptr = boost::make_shared<MatrixDouble>();
  auto kappa_increment_ptr = boost::make_shared<VectorDouble>();
 
  addTetrahedronOperator(fe, new OpCalculateVectorFieldValues<
                                 plasticLogarithmicStretchCoordinateSize>(
                                 control_dm, PlasticIncrementalOp::OPCOL,
                                 plastic_field, plastic_increment_ptr,
                                 control_vec, MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 control_dm, PlasticIncrementalOp::OPCOL,
                                 kappa_field, kappa_increment_ptr, control_vec,
                                 MBTET));
 
  const double regularization_epsilon =
      problem.getDissipationRegularizationEpsilon();
 
  auto *constraint_op =
      new PlasticIncrementalOp(kappa_field, PlasticIncrementalOp::OPROW);
  constraint_op->doWorkRhsHook =
      [&](DataOperator *base_op_ptr, int, EntityType,
          EntitiesFieldData::EntData &row_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    const int nb_integration_pts = op_ptr->getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_increment = MatrixSizeHelper<
        GetFTensor1FromMatType<plasticLogarithmicStretchCoordinateSize, -1, DL>,
        DL>::get(*plastic_increment_ptr, nb_integration_pts);
    const auto t_increment = get_increment();
    const auto t_kappa_increment = getFTensor0FromVec(*kappa_increment_ptr);
    const double constraint_scale = std::sqrt(getCellMeasure(*op_ptr));
    // Scale the exact or rounded P0 cone row by the reference-volume root.
    const double local_constraint =
        constraint_scale *
        (t_kappa_increment -
         plasticEquivalentIncrement(plasticCoordinateNorm(t_increment),
                                    regularization_epsilon));
    CHKERR VecSetValues<AssemblyTypeSelector<PETSC>>(
        op_ptr->getKSPf(), row_data, &local_constraint, INSERT_VALUES);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, constraint_op);
 
  CHKERR DMoFEMLoopFiniteElements(constraint_dm, problem.getVolumeElementName(),
                                  fe);
  CHKERR VecAssemblyBegin(constraints);
  CHKERR VecAssemblyEnd(constraints);
  CHKERR VecGhostUpdateBegin(constraints, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(constraints, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode assemblePlasticInequalityJacobian(
    PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec smooth_gradient, Mat jacobian) {
  MoFEMFunctionBegin;
  const auto &plastic_field = problem.getPlasticFlowFieldName();
  const auto &kappa_field = problem.getPlasticKappaFieldName();
  DM constraint_dm = nullptr;
  CHKERR MatGetDM(jacobian, &constraint_dm);
  if (!constraint_dm)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Plastic inequality Jacobian has no constraint DM");
 
  auto fe = createPlasticIncrementalElement(problem);
  fe->ksp_B = jacobian;
  const double regularization_epsilon =
      problem.getDissipationRegularizationEpsilon();
 
  const auto control_dm = SmartPetscObj<DM>(problem.getControlDM(), true);
  const auto control_vec = SmartPetscObj<Vec>(control, true);
  const auto gradient_vec = SmartPetscObj<Vec>(smooth_gradient, true);
  auto increment_ptr = boost::make_shared<MatrixDouble>();
  auto gradient_ptr = boost::make_shared<MatrixDouble>();
  auto committed_kappa_ptr = boost::make_shared<VectorDouble>();
 
  addTetrahedronOperator(fe,
                         new OpCalculateVectorFieldValues<
                             plasticLogarithmicStretchCoordinateSize>(
                             control_dm, PlasticIncrementalOp::OPCOL,
                             plastic_field, increment_ptr, control_vec, MBTET));
  addTetrahedronOperator(fe,
                         new OpCalculateVectorFieldValues<
                             plasticLogarithmicStretchCoordinateSize>(
                             control_dm, PlasticIncrementalOp::OPCOL,
                             plastic_field, gradient_ptr, gradient_vec, MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 kappa_field, committed_kappa_ptr, MBTET));
 
  auto *flow_op = new PlasticIncrementalOp(
      kappa_field, plastic_field, PlasticIncrementalOp::OPROWCOL, false);
  flow_op->doWorkLhsHook =
      [&](DataOperator *base_op_ptr, int, int, EntityType, EntityType,
          EntitiesFieldData::EntData &row_data,
          EntitiesFieldData::EntData &col_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty() || col_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    const int nb_integration_pts = op_ptr->getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_increment = MatrixSizeHelper<
        GetFTensor1FromMatType<plasticLogarithmicStretchCoordinateSize, -1, DL>,
        DL>::get(*increment_ptr, nb_integration_pts);
    auto get_gradient = MatrixSizeHelper<
        GetFTensor1FromMatType<plasticLogarithmicStretchCoordinateSize, -1, DL>,
        DL>::get(*gradient_ptr, nb_integration_pts);
    const auto t_increment = get_increment();
    const auto t_gradient = get_gradient();
    const auto t_committed_kappa = getFTensor0FromVec(*committed_kappa_ptr);
    // These P0 producer values are element-constant, so the first Gauss-point
    // view supplies the complete local constitutive block.
    const double increment_norm = plasticCoordinateNorm(t_increment);
    FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
        t_constraint_slope;
    FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
    const double gradient_denominator =
        std::hypot(plasticEquivalentIncrementScale * increment_norm,
                   regularization_epsilon);
    // Eq. (1.83), label eq:regularised-dissipation-gradient, when epsilon>0;
    // otherwise the exact selection in Eq. (1.74), label
    // eq:exact-constraint-branch-selection.
    if (gradient_denominator > 0) {
      t_constraint_slope(L) = plasticEquivalentIncrementScaleSquared *
                              t_increment(L) / gradient_denominator;
    } else {
      const double cell_measure = getCellMeasure(*op_ptr);
      if (!(cell_measure > 0))
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Plastic P0 block has no positive cell measure");
      FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
          t_cell_force;
      t_cell_force(L) = -t_gradient(L) / cell_measure;
      const double q =
          plasticYieldFunctionScale * plasticCoordinateNorm(t_cell_force);
      const double current_yield =
          problem.getInitialYieldStress() +
          problem.getIsotropicHardeningModulus() * t_committed_kappa;
      t_constraint_slope(L) = t_cell_force(L) / std::max(current_yield, q);
    }
 
    const double constraint_scale = std::sqrt(getCellMeasure(*op_ptr));
    MatrixDouble local_matrix(1, plasticLogarithmicStretchCoordinateSize);
    auto t_local_matrix =
        getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
            local_matrix.data().data());
    t_local_matrix(L) = -constraint_scale * t_constraint_slope(L);
    CHKERR MatSetValues<AssemblyTypeSelector<PETSC>>(
        op_ptr->getKSPB(), row_data, col_data, local_matrix.data().data(),
        ADD_VALUES);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, flow_op);
 
  auto *kappa_op = new PlasticIncrementalOp(
      kappa_field, kappa_field, PlasticIncrementalOp::OPROWCOL, false);
  kappa_op->doWorkLhsHook =
      [&](DataOperator *base_op_ptr, int, int, EntityType, EntityType,
          EntitiesFieldData::EntData &row_data,
          EntitiesFieldData::EntData &col_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty() || col_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    const double constraint_scale = std::sqrt(getCellMeasure(*op_ptr));
    MatrixDouble local_matrix(1, 1);
    local_matrix(0, 0) = constraint_scale;
    CHKERR MatSetValues<AssemblyTypeSelector<PETSC>>(
        op_ptr->getKSPB(), row_data, col_data, local_matrix.data().data(),
        ADD_VALUES);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, kappa_op);
 
  CHKERR DMoFEMLoopFiniteElements(constraint_dm, problem.getVolumeElementName(),
                                  fe);
  CHKERR MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode validatePlasticIncrementalSolution(
    PlasticIncrementalOptimizationProblem &problem, DM constraint_dm,
    Vec control, Vec smooth_gradient, Vec inequality_multipliers,
    PetscReal gradient_tolerance, PetscReal constraint_tolerance,
    std::string &diagnostics) {
  MoFEMFunctionBegin;
  if (!constraint_dm || !control || !smooth_gradient ||
      !inequality_multipliers)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Plastic solution validation requires control and constraint DMs "
            "with control, gradient, and multiplier vectors");
 
  const auto &plastic_field = problem.getPlasticFlowFieldName();
  const auto &kappa_field = problem.getPlasticKappaFieldName();
  const auto control_dm = SmartPetscObj<DM>(problem.getControlDM(), true);
  const auto constraint_data_dm = SmartPetscObj<DM>(constraint_dm, true);
  const auto control_vec = SmartPetscObj<Vec>(control, true);
  const auto gradient_vec = SmartPetscObj<Vec>(smooth_gradient, true);
  const auto multiplier_vec = SmartPetscObj<Vec>(inequality_multipliers, true);
 
  auto plastic_increment_ptr = boost::make_shared<MatrixDouble>();
  auto kappa_increment_ptr = boost::make_shared<VectorDouble>();
  auto gradient_ptr = boost::make_shared<MatrixDouble>();
  auto multiplier_ptr = boost::make_shared<VectorDouble>();
  auto committed_kappa_ptr = boost::make_shared<VectorDouble>();
  auto metrics_ptr = boost::make_shared<PlasticValidationMetrics>();
 
  auto fe = createPlasticIncrementalElement(problem);
  addTetrahedronOperator(fe, new OpCalculateVectorFieldValues<
                                 plasticLogarithmicStretchCoordinateSize>(
                                 control_dm, PlasticIncrementalOp::OPCOL,
                                 plastic_field, plastic_increment_ptr,
                                 control_vec, MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 control_dm, PlasticIncrementalOp::OPCOL,
                                 kappa_field, kappa_increment_ptr, control_vec,
                                 MBTET));
  addTetrahedronOperator(fe,
                         new OpCalculateVectorFieldValues<
                             plasticLogarithmicStretchCoordinateSize>(
                             control_dm, PlasticIncrementalOp::OPCOL,
                             plastic_field, gradient_ptr, gradient_vec, MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 constraint_data_dm,
                                 PlasticIncrementalOp::OPROW, kappa_field,
                                 multiplier_ptr, multiplier_vec, MBTET));
  addTetrahedronOperator(fe, new OpCalculateScalarFieldValues(
                                 kappa_field, committed_kappa_ptr, MBTET));
 
  auto *validation_op =
      new PlasticIncrementalOp(kappa_field, PlasticIncrementalOp::OPROW);
  validation_op->doWorkRhsHook =
      [&, metrics_ptr](DataOperator *base_op_ptr, int, EntityType,
                       EntitiesFieldData::EntData &row_data) -> MoFEMErrorCode {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    auto *op_ptr = static_cast<PlasticIncrementalOp *>(base_op_ptr);
    const int nb_integration_pts = op_ptr->getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_increment = MatrixSizeHelper<
        GetFTensor1FromMatType<plasticLogarithmicStretchCoordinateSize, -1, DL>,
        DL>::get(*plastic_increment_ptr, nb_integration_pts);
    auto get_gradient = MatrixSizeHelper<
        GetFTensor1FromMatType<plasticLogarithmicStretchCoordinateSize, -1, DL>,
        DL>::get(*gradient_ptr, nb_integration_pts);
    const auto t_increment = get_increment();
    const auto t_gradient = get_gradient();
    const double delta_kappa = getFTensor0FromVec(*kappa_increment_ptr);
    const double kappa_n = getFTensor0FromVec(*committed_kappa_ptr);
    const double cell_measure = getCellMeasure(*op_ptr);
    // Recover the multiplier of the unscaled cone for physical KKT checks.
    const double multiplier =
        std::sqrt(cell_measure) * getFTensor0FromVec(*multiplier_ptr);
    const double current_yield =
        problem.getInitialYieldStress() +
        problem.getIsotropicHardeningModulus() * (kappa_n + delta_kappa);
 
    FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
    bool finite_block =
        std::isfinite(cell_measure) && std::isfinite(delta_kappa) &&
        std::isfinite(kappa_n) && kappa_n >= 0 &&
        std::isfinite(current_yield) &&
        current_yield > 0 && std::isfinite(multiplier) &&
        std::isfinite(plasticCoordinateNorm(t_increment)) &&
        std::isfinite(plasticCoordinateNorm(t_gradient));
    if (!finite_block) {
      metrics_ptr->nonfinite = 1;
      MoFEMFunctionReturnHot(0);
    }
    if (!(cell_measure > std::numeric_limits<double>::epsilon())) {
      metrics_ptr->invalidMeasure = 1;
      MoFEMFunctionReturnHot(0);
    }
 
    FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
        t_force_coordinates;
    t_force_coordinates(L) = -t_gradient(L) / cell_measure;
    const auto t_flow =
        plasticLogarithmicStretchTensorFromCoordinates(t_increment);
    const auto t_force =
        plasticLogarithmicStretchTensorFromCoordinates(t_force_coordinates);
    FTENSOR_INDEXES(SPACE_DIM, i, j);
    const double work = cell_measure * t_force(i, j) * t_flow(i, j);
    const double q =
        plasticYieldFunctionScale * plasticFrobeniusNorm(t_force);
    const double flow_norm = plasticCoordinateNorm(t_increment);
    const double equivalent_increment = plasticEquivalentIncrement(
        flow_norm, problem.getDissipationRegularizationEpsilon());
    const double constraint = delta_kappa - equivalent_increment;
    FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize> t_slope;
    const double gradient_denominator = std::hypot(
        plasticEquivalentIncrementScale * flow_norm,
        problem.getDissipationRegularizationEpsilon());
    if (gradient_denominator > 0)
      t_slope(L) = plasticEquivalentIncrementScaleSquared * t_increment(L) /
                   gradient_denominator;
    else
      t_slope(L) =
          t_force_coordinates(L) / std::max(current_yield, q);
    FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
        t_kkt_residual;
    // Unscaled c and its physical multiplier retain the original tolerances.
    t_kkt_residual(L) = t_gradient(L) + multiplier * t_slope(L);
    const double eta_stationarity_residual =
        std::abs(cell_measure * current_yield - multiplier);
    const double complementarity_residual =
        std::abs(multiplier * constraint);
 
    metrics_ptr->minWork = std::min(metrics_ptr->minWork, work);
    metrics_ptr->maxYieldFunction =
        std::max(metrics_ptr->maxYieldFunction, q);
    metrics_ptr->maxYieldExcess =
        std::max(metrics_ptr->maxYieldExcess, q - current_yield);
    metrics_ptr->minConstraint =
        std::min(metrics_ptr->minConstraint, constraint);
    metrics_ptr->minMultiplier =
        std::min(metrics_ptr->minMultiplier, multiplier);
    metrics_ptr->maxActiveConstraint =
        std::max(metrics_ptr->maxActiveConstraint, std::abs(constraint));
    metrics_ptr->maxEtaStationarityResidual =
        std::max(metrics_ptr->maxEtaStationarityResidual,
                 eta_stationarity_residual);
    metrics_ptr->maxComplementarityResidual =
        std::max(metrics_ptr->maxComplementarityResidual,
                 complementarity_residual);
    metrics_ptr->maxKktResidual =
        std::max(metrics_ptr->maxKktResidual,
                 plasticCoordinateNorm(t_kkt_residual));
    metrics_ptr->maxEquivalentIncrement =
        std::max(metrics_ptr->maxEquivalentIncrement, equivalent_increment);
    MoFEMFunctionReturn(0);
  };
  addTetrahedronOperator(fe, validation_op);
 
  CHKERR DMoFEMLoopFiniteElements(constraint_dm,
                                  problem.getVolumeElementName(), fe);
 
  PlasticValidationMetrics global_metrics;
  const MPI_Comm comm =
      PetscObjectComm(reinterpret_cast<PetscObject>(constraint_dm));
  const auto reduce_max = [&](const double local, double &global) {
    return MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_MAX, comm);
  };
  const auto reduce_min = [&](const double local, double &global) {
    return MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_MIN, comm);
  };
  CHKERR reduce_min(metrics_ptr->minWork, global_metrics.minWork);
  CHKERR reduce_max(metrics_ptr->maxYieldFunction,
                    global_metrics.maxYieldFunction);
  CHKERR reduce_max(metrics_ptr->maxYieldExcess,
                    global_metrics.maxYieldExcess);
  CHKERR reduce_min(metrics_ptr->minConstraint, global_metrics.minConstraint);
  CHKERR reduce_min(metrics_ptr->minMultiplier, global_metrics.minMultiplier);
  CHKERR reduce_max(metrics_ptr->maxActiveConstraint,
                    global_metrics.maxActiveConstraint);
  CHKERR reduce_max(metrics_ptr->maxEtaStationarityResidual,
                    global_metrics.maxEtaStationarityResidual);
  CHKERR reduce_max(metrics_ptr->maxComplementarityResidual,
                    global_metrics.maxComplementarityResidual);
  CHKERR reduce_max(metrics_ptr->maxKktResidual,
                    global_metrics.maxKktResidual);
  CHKERR reduce_max(metrics_ptr->maxEquivalentIncrement,
                    global_metrics.maxEquivalentIncrement);
  CHKERR MPI_Allreduce(&metrics_ptr->nonfinite, &global_metrics.nonfinite, 1,
                       MPI_INT, MPI_MAX, comm);
  CHKERR MPI_Allreduce(&metrics_ptr->invalidMeasure,
                       &global_metrics.invalidMeasure, 1, MPI_INT, MPI_MAX,
                       comm);
 
  const double kkt_tolerance =
      10 * std::max<double>(gradient_tolerance, 1e-10);
  const double feasibility_tolerance =
      10 * std::max<double>(constraint_tolerance, 1e-10);
  const bool plastic_flow_active =
      global_metrics.maxEquivalentIncrement > feasibility_tolerance;
  if (global_metrics.nonfinite || global_metrics.invalidMeasure ||
      global_metrics.minWork < -kkt_tolerance ||
      global_metrics.minConstraint < -feasibility_tolerance ||
      global_metrics.minMultiplier < -kkt_tolerance ||
      global_metrics.maxActiveConstraint > feasibility_tolerance ||
      global_metrics.maxEtaStationarityResidual > kkt_tolerance ||
      global_metrics.maxComplementarityResidual > kkt_tolerance ||
      global_metrics.maxKktResidual > kkt_tolerance)
    SETERRQ(comm, MOFEM_OPERATION_UNSUCCESSFUL,
            "Incremental-optimization plastic acceptance failed: non-finite "
            "%d, invalid measure %d, minimum work %g, max "
            "yield excess %g, minimum constraint %g, minimum multiplier %g, "
            "active residual %g, eta stationarity %g, complementarity %g, "
            "plastic KKT residual %g",
            global_metrics.nonfinite, global_metrics.invalidMeasure,
            global_metrics.minWork,
            global_metrics.maxYieldExcess, global_metrics.minConstraint,
            global_metrics.minMultiplier, global_metrics.maxActiveConstraint,
            global_metrics.maxEtaStationarityResidual,
            global_metrics.maxComplementarityResidual,
            global_metrics.maxKktResidual);
 
  std::ostringstream stream;
  stream << "max yield function " << global_metrics.maxYieldFunction
         << ", max yield excess " << global_metrics.maxYieldExcess
         << ", min plastic work " << global_metrics.minWork
         << ", min constraint " << global_metrics.minConstraint
         << ", min multiplier " << global_metrics.minMultiplier
         << ", active residual " << global_metrics.maxActiveConstraint
         << ", eta stationarity "
         << global_metrics.maxEtaStationarityResidual
         << ", complementarity "
         << global_metrics.maxComplementarityResidual
         << ", plastic KKT residual " << global_metrics.maxKktResidual
         << ", max equivalent plastic increment "
         << global_metrics.maxEquivalentIncrement << ", plastic flow active "
         << (plastic_flow_active ? "true" : "false")
         << ", dissipation epsilon "
         << problem.getDissipationRegularizationEpsilon();
  diagnostics = stream.str();
  MoFEMFunctionReturn(0);
}
 
} // namespace PlasticIncrementalOptimizationInternal
 
} // namespace EshelbianPlasticity