PlasticIncrementalOptimization.cpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalOptimization.cpp
 * @brief Incremental plasticity optimization
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <Lie.hpp>
#include <MatrixFunction.hpp>
#include <PlasticIncrementalOptimization.hpp>
#include <impl/PlasticIncrementalOptimizationInternal.hpp>
 
namespace EshelbianPlasticity {
 
using namespace PlasticIncrementalOptimizationInternal;
 
namespace PlasticIncrementalOptimizationInternal {
 
MoFEMErrorCode savePlasticHistory(EshelbianCore &ep, PlasticHistory &history) {
  MoFEMFunctionBegin;
  // Snapshot H_p^n from the quadrature history in Eq. (1.54), label
  // eq:quadrature-history-update.
  history.clear();
  auto save = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
    MoFEMFunctionBegin;
    auto field_data = field_entity_ptr->getEntFieldData();
    if (field_data.size() != plasticLogarithmicStretchCoordinateSize)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Equilibrated-value test requires five trace-free P0 plasticH "
              "coordinates per cell; got %zu",
              field_data.size());
    PackedPlasticCoordinates values;
    std::copy(field_data.begin(), field_data.end(), values.begin());
    history.emplace(field_entity_ptr->getEnt(), values);
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      save, ep.plasticHField, ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setPlasticHistory(
    EshelbianCore &ep, const PlasticHistory &history,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_direction,
    const double scale) {
  MoFEMFunctionBegin;
  // Apply the logarithmic update H_p^{n+1}=H_p^n+Delta H_p from Eq. (1.4),
  // label eq:logarithmic-update, to the saved P0 coefficients.
  const auto t_direction_coordinates =
      plasticLogarithmicStretchCoordinatesFromTensor(t_direction);
  auto set = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
    MoFEMFunctionBegin;
    const auto history_it = history.find(field_entity_ptr->getEnt());
    if (history_it == history.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Missing saved plasticH data for entity %llu",
              static_cast<unsigned long long>(field_entity_ptr->getEnt()));
    auto field_data = field_entity_ptr->getEntFieldData();
    if (field_data.size() != history_it->second.size())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Changed plasticH coefficient count while restoring history");
    std::copy(history_it->second.begin(), history_it->second.end(),
              field_data.begin());
    auto t_field = getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
        &field_data[0]);
    FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
    t_field(L) += scale * t_direction_coordinates(L);
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      set, ep.plasticHField, ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode savePlasticKappaHistory(EshelbianCore &ep,
                                       PlasticKappaHistory &history) {
  MoFEMFunctionBegin;
  history.clear();
  auto save = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
    MoFEMFunctionBegin;
    auto field_data = field_entity_ptr->getEntFieldData();
    if (field_data.size() != 1)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Committed plastic kappa requires one P0 coefficient per "
              "cell; got %zu",
              field_data.size());
    const double kappa = field_data[0];
    if (!std::isfinite(kappa) || kappa < 0)
      SETERRQ(ep.mField.get_comm(), MOFEM_INVALID_DATA,
              "Committed plastic kappa on entity %llu must be finite and "
              "non-negative; got %g",
              static_cast<unsigned long long>(field_entity_ptr->getEnt()),
              kappa);
    history.emplace(field_entity_ptr->getEnt(), kappa);
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      save, ep.plasticKappaField, ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setPlasticKappaHistory(EshelbianCore &ep,
                                      const PlasticKappaHistory &history) {
  MoFEMFunctionBegin;
  auto set = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
    MoFEMFunctionBegin;
    const auto history_it = history.find(field_entity_ptr->getEnt());
    if (history_it == history.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Missing saved plastic kappa for entity %llu",
              static_cast<unsigned long long>(field_entity_ptr->getEnt()));
    auto field_data = field_entity_ptr->getEntFieldData();
    if (field_data.size() != 1)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Changed plastic-kappa coefficient count while restoring "
              "history");
    field_data[0] = history_it->second;
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      set, ep.plasticKappaField, ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode clearPlasticIncrementFields(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnValues(
      [](const double) { return 0.; }, ep.plasticFlowField,
      ep.plasticVolumes.get());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setStateOnMesh(EshelbianCore &ep, Vec state) {
  MoFEMFunctionBegin;
  CHKERR VecGhostUpdateBegin(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, state, INSERT_VALUES,
                                 SCATTER_REVERSE, RowColData::COL);
  MoFEMFunctionReturn(0);
}
 
} // namespace PlasticIncrementalOptimizationInternal
 
PlasticIncrementalOptimizationProblem::PlasticIncrementalOptimizationProblem(
    EshelbianCore &ep, SmartPetscObj<TS> ts)
    : ep(ep), timeSolver(std::move(ts)) {}
 
std::string PlasticIncrementalOptimizationProblem::getName() const {
  return "plasticity";
}
 
DM PlasticIncrementalOptimizationProblem::getControlDM() const {
  return ep.dmIncrementalOptimization;
}
 
const std::string &
PlasticIncrementalOptimizationProblem::getPlasticFlowFieldName() const {
  return ep.plasticFlowField;
}
 
const std::string &
PlasticIncrementalOptimizationProblem::getPlasticKappaFieldName() const {
  return ep.plasticKappaField;
}
 
const std::string &
PlasticIncrementalOptimizationProblem::getVolumeElementName() const {
  return ep.elementVolumeName;
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::addConstitutiveGeometryOperators(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline) const {
  // Only transformed quadrature weights are needed. Keep the P0 bases
  // unchanged while evaluating vector-owned trial coefficients.
  return MoFEM::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pipeline, {NOSPACE}, ep.materialH1Positions);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::validate() {
  MoFEMFunctionBegin;
  if (!ep.plasticVolume || ep.interfaceCrack)
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Plastic incremental optimization requires -plastic_volume 1 "
            "and -cohesive_interface_on 0");
  if (!ep.dmIncrementalOptimization)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic incremental-optimization control DM is not available");
  if (!timeSolver)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic incremental optimization requires a state TS");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::createControlData(
    SmartPetscObj<Vec> &committed_control,
    SmartPetscObj<Vec> &candidate_control) {
  MoFEMFunctionBegin;
  committed_control =
      createDMVector(ep.dmIncrementalOptimization, RowColData::ROW);
  if (!ep.incrementalTrialControl)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic trial-control vector is not available");
  candidate_control = ep.incrementalTrialControl;
  // Mesh field data store committed kappa, whereas every vector with this DM
  // layout stores fresh increments (Delta H_p, Delta kappa).
  CHKERR VecZeroEntries(committed_control);
  CHKERR VecZeroEntries(candidate_control);
  fallbackJacobianSmoothGradient = vectorDuplicate(committed_control);
  CHKERR VecZeroEntries(fallbackJacobianSmoothGradient);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "",
                             "-incremental_optimization_yield_stress",
                             &initialYieldStress, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(
      PETSC_NULLPTR, "",
      "-incremental_optimization_isotropic_hardening_modulus",
      &isotropicHardeningModulus, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "",
                             "-incremental_optimization_dissipation_epsilon",
                             &dissipationRegularizationEpsilon, PETSC_NULLPTR);
  if (!std::isfinite(initialYieldStress) || !(initialYieldStress > 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_INVALID_DATA,
            "Initial plastic yield stress must be finite and positive");
  if (!std::isfinite(isotropicHardeningModulus) ||
      isotropicHardeningModulus < 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_INVALID_DATA,
            "Isotropic hardening modulus must be finite and non-negative");
  if (!std::isfinite(dissipationRegularizationEpsilon) ||
      dissipationRegularizationEpsilon < 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_INVALID_DATA,
            "Plastic dissipation epsilon must be finite and non-negative");
  PlasticKappaHistory kappa_history;
  CHKERR savePlasticKappaHistory(ep, kappa_history);
  MOFEM_LOG("EP", Sev::inform)
      << "Plastic model parameters: initial yield stress "
      << initialYieldStress << ", isotropic hardening modulus "
      << isotropicHardeningModulus
      << ", dissipation epsilon " << dissipationRegularizationEpsilon
      << ", constraint scaling sqrt(reference cell volume)";
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::configureConstraintDM(
    DM constraint_dm) {
  MoFEMFunctionBegin;
  if (!constraint_dm)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic constraint DM is null");
  CHKERR DMMoFEMAddSubFieldRow(constraint_dm, ep.plasticKappaField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddSubFieldCol(constraint_dm, ep.plasticFlowField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddSubFieldCol(constraint_dm, ep.plasticKappaField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddElement(constraint_dm, ep.elementVolumeName);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::initialiseInequalityMultipliers(
    DM constraint_dm, Vec multipliers) {
  return initialisePlasticInequalityMultipliers(*this, constraint_dm,
                                                multipliers);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::initialiseReferenceControl(
    Vec committed_control, Vec reference_control) {
  MoFEMFunctionBegin;
  (void)committed_control;
  // Both Delta H_p and Delta kappa are fresh increment variables. Committed
  // kappa remains in mesh field storage.
  CHKERR VecZeroEntries(reference_control);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::solveState(
    Vec state, EquilibriumSolveReport &report) {
  return solveEquilibriumStateTS(ep, timeSolver, state, report);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::restoreState(Vec state) {
  return setStateOnMesh(ep, state);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::evaluateConservativeValue(
    double &value) {
  MoFEMFunctionBegin;
  CHKERR evaluateEquilibratedMechanicalValue(ep, timeSolver, value);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::assembleConservativeGradient(
    Vec smooth_gradient) {
  return assembleReducedHelmholtzGradient(ep, timeSolver, smooth_gradient);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::evaluateDissipation(Vec control,
                                                           double &value) {
  // The generic "dissipation" slot is the complete local incremental
  // resistance: irreversible yield work plus stored hardening-energy change.
  return evaluatePlasticIncrementalResistanceValue(*this, control, value);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::assembleObjectiveGradient(
    Vec control, Vec smooth_gradient, Vec objective_gradient) {
  return assembleIncrementalObjectiveGradient(*this, control, smooth_gradient,
                                              objective_gradient);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::evaluateInequalityConstraints(
    Vec control, Vec constraints) {
  return evaluatePlasticInequalityConstraints(*this, control, constraints);
}
 
MoFEMErrorCode
PlasticIncrementalOptimizationProblem::assembleInequalityJacobian(
    Vec control, Vec smooth_gradient, Mat jacobian) {
  return assemblePlasticInequalityJacobian(*this, control, smooth_gradient,
                                           jacobian);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::
    assembleInequalityJacobianWithoutSmoothGradient(Vec control,
                                                     Mat jacobian) {
  MoFEMFunctionBegin;
  if (!fallbackJacobianSmoothGradient)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Fallback plastic Jacobian scratch vector is unavailable");
  CHKERR VecZeroEntries(fallbackJacobianSmoothGradient);
  // At the exact cone apex, zero is an admissible subgradient. Away from the
  // apex, and on the rounded branch, the row depends only on the control.
  CHKERR assemblePlasticInequalityJacobian(
      *this, control, fallbackJacobianSmoothGradient, jacobian);
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<IncrementalOptimizationContext>
createPlasticIncrementalOptimizationContext(EshelbianCore &ep,
                                            SmartPetscObj<TS> ts,
                                            SmartPetscObj<Vec> state) {
  auto problem = boost::make_shared<PlasticIncrementalOptimizationProblem>(
      ep, std::move(ts));
  return createIncrementalOptimizationContext(problem, std::move(state));
}
 
PlasticIncrementalOptimizationProblem &getPlasticIncrementalOptimizationProblem(
    const boost::shared_ptr<IncrementalOptimizationContext> &context) {
  if (!context || !context->problem)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Incremental-optimization context has no problem");
  auto plastic_problem =
      boost::dynamic_pointer_cast<PlasticIncrementalOptimizationProblem>(
          context->problem);
  if (!plastic_problem)
    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                      "Incremental-optimization problem is not plasticity");
  return *plastic_problem;
}
 
MoFEMErrorCode validateEquilibratedMechanicalValueScope(EshelbianCore &ep) {
  MoFEMFunctionBegin;
 
  if (!ep.physicalEquations ||
      !ep.physicalEquations->providesHelmholtzFreeEnergy())
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Selected material does not provide the Helmholtz free energy "
            "required by the equilibrated mechanical value");
  if ((ep.physicalEquations->getFeatures() & PhysicalEquations::noStretchMask)
          .any() ||
      ep.stretchSelector != StretchSelector::LOG ||
      ep.gradApproximator != RotSelector::NO_H1_CONFIGURATION ||
      ep.rotSelector != RotSelector::LARGE_ROT)
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Equilibrated mechanical value requires explicit logarithmic "
            "stretch, -grad no_h1, -rotations large, and -no_stretch 0");
 
  // Keep this threshold identical to setElasticElementToTs(), which selects
  // I2Function as soon as alphaRho exceeds machine epsilon.
  if (std::abs(ep.alphaRho) > std::numeric_limits<double>::epsilon())
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Incremental optimization currently supports first-order "
            "pseudo-time relaxation only; set -density_alpha_rho 0");
 
  if (ep.externalStrainVecPtr && !ep.externalStrainVecPtr->empty())
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "The material Helmholtz derivative with imposed external strain "
            "is not implemented");
  if (ep.crackingOn || (ep.contactFaces && !ep.contactFaces->empty()))
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Equilibrated mechanical value does not include cracking or "
            "contact contributions");
  if (!ep.internalStressTagName.empty())
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Equilibrated mechanical value does not include initial "
            "internal stress");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode solveEquilibriumStateTS(EshelbianCore &ep, TS ts, Vec state,
                                       EquilibriumSolveReport &report) {
  MoFEMFunctionBegin;
  // Fixed-control state map S(q) from Eq. (3.5), label gen:eq:state-map; its
  // plastic specialization is Eq. (1.56), label eq:discrete-state-solve.
  report = {};
  if (!ts || !state)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Fixed-control equilibrium requires a TS and a state vector");
 
  PetscReal initial_dt = 0;
  CHKERR TSGetTimeStep(ts, &initial_dt);
  // currentPhysicalTime is the fixed outer load coordinate. TS time is a
  // local pseudo-time, reset for every TAO trial, and all rates in the state
  // equations are derivatives with respect to this pseudo-time.
  CHKERR TSSetStepNumber(ts, 0);
  CHKERR TSSetTime(ts, 0);
  CHKERR TSSetSolution(ts, state);
  CHKERR TSRestartStep(ts);
  SNES snes = nullptr;
  PetscBool previous_force_iteration = PETSC_FALSE;
  CHKERR TSGetSNES(ts, &snes);
  CHKERR SNESGetForceIteration(snes, &previous_force_iteration);
  // Re-equilibrate every distinct TAO control. Taking at least one Newton step
  // keeps the reduced Helmholtz value and adjoint gradient consistent for
  // line search even when the previous state's residual is below the absolute
  // SNES tolerance.
  CHKERR SNESSetForceIteration(snes, PETSC_TRUE);
  const MoFEMErrorCode solve_error = TSSolve(ts, PETSC_NULLPTR);
  const MoFEMErrorCode restore_force_error =
      SNESSetForceIteration(snes, previous_force_iteration);
  const MoFEMErrorCode restore_dt_error = TSSetTimeStep(ts, initial_dt);
  CHKERR restore_force_error;
  CHKERR restore_dt_error;
  if (solve_error) {
    report.status = EquilibriumSolveStatus::Failed;
    report.errorCode = solve_error;
    MoFEMFunctionReturnHot(0);
  }
 
  TSConvergedReason reason;
  CHKERR TSGetConvergedReason(ts, &reason);
  if (reason < 0) {
    report.status = EquilibriumSolveStatus::Diverged;
    report.errorCode = MOFEM_OPERATION_UNSUCCESSFUL;
    MoFEMFunctionReturnHot(0);
  }
 
  PetscReal pseudo_time = 0;
  PetscReal final_pseudo_time = 0;
  CHKERR TSGetTime(ts, &pseudo_time);
  CHKERR TSGetMaxTime(ts, &final_pseudo_time);
  const double tolerance =
      100. * std::numeric_limits<double>::epsilon() *
      std::max(1., std::abs(static_cast<double>(final_pseudo_time)));
  if (pseudo_time + tolerance < final_pseudo_time) {
    if (reason != TS_CONVERGED_ITS) {
      report.status = EquilibriumSolveStatus::Diverged;
      report.errorCode = MOFEM_OPERATION_UNSUCCESSFUL;
      MoFEMFunctionReturnHot(0);
    }
    PetscInt steps = 0;
    CHKERR TSGetStepNumber(ts, &steps);
    MOFEM_LOG("EP", Sev::warning)
        << "Fixed-control equilibrium reached the TS step limit after " << steps
        << " steps at pseudo-time " << pseudo_time << " (target "
        << final_pseudo_time
        << "); continuing with the final state and checking its rate-free "
           "equilibrium residual";
  }
 
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, state, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(state, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state, INSERT_VALUES, SCATTER_FORWARD);
 
  // The last SNES norm belongs to the time-discrete, rate-dependent equation.
  // It can be small while damping still balances a nonzero mechanical
  // residual. The reduced Helmholtz gradient
  // is valid only at the rate-free stationary state, so assemble the same
  // IFunction with a zero pseudo-time rate and report that residual instead.
  auto zero_pseudo_rate = vectorDuplicate(state);
  auto rate_free_residual = vectorDuplicate(state);
  CHKERR VecZeroEntries(zero_pseudo_rate);
  CHKERR TSComputeIFunction(ts, pseudo_time, state, zero_pseudo_rate,
                           rate_free_residual, PETSC_FALSE);
  PetscReal pseudo_time_equation_norm = 0;
  CHKERR SNESGetFunctionNorm(snes, &pseudo_time_equation_norm);
  PetscReal rate_free_norm = 0;
  CHKERR VecNorm(rate_free_residual, NORM_2, &rate_free_norm);
  MOFEM_LOG("EP", Sev::verbose)
      << "Pseudo-time equation residual " << pseudo_time_equation_norm
      << ", rate-free equilibrium residual " << rate_free_norm
      << " at pseudo-time " << pseudo_time;
  report.equilibriumNorm = rate_free_norm;
  report.status = EquilibriumSolveStatus::Converged;
  MoFEMFunctionReturn(0);
}
 
namespace PlasticIncrementalOptimizationInternal {
 
namespace {
 
MoFEMErrorCode synchroniseCommitStageError(MPI_Comm comm,
                                            MoFEMErrorCode local_error,
                                            MoFEMErrorCode &global_error) {
  MoFEMFunctionBegin;
  const int local_code = static_cast<int>(local_error);
  int global_code = 0;
  CHKERR MPI_Allreduce(&local_code, &global_code, 1, MPI_INT, MPI_MAX, comm);
  global_error = static_cast<MoFEMErrorCode>(global_code);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
MoFEMErrorCode commitIncrementalOptimizationLoadStep(
    EshelbianCore &ep, Vec solution, Vec equilibrated_state,
    Vec committed_control, Vec candidate_control, Vec reference_control) {
  MoFEMFunctionBegin;
  // Commit Delta H_p and Delta kappa only after the complete solution passed
  // validation. Trial vectors never overwrite either mesh history.
  PlasticHistory plastic_before;
  PlasticKappaHistory kappa_before;
  MoFEMErrorCode local_stage_error =
      ep.plasticVolume ? savePlasticHistory(ep, plastic_before) : 0;
  MoFEMErrorCode stage_error = 0;
  CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                     stage_error);
  CHKERR stage_error;
  local_stage_error =
      ep.plasticVolume ? savePlasticKappaHistory(ep, kappa_before) : 0;
  CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                     stage_error);
  CHKERR stage_error;
  auto control_before = vectorDuplicate(committed_control);
  auto candidate_before = vectorDuplicate(candidate_control);
  auto reference_before = vectorDuplicate(reference_control);
  auto state_before = vectorDuplicate(equilibrated_state);
  auto mesh_control_before =
      createDMVector(ep.dmIncrementalOptimization, RowColData::ROW);
  const auto record_first_error = [](MoFEMErrorCode &first_error,
                                     MoFEMErrorCode next_error) {
    if (!first_error)
      first_error = next_error;
  };
  local_stage_error = VecCopy(committed_control, control_before);
  record_first_error(local_stage_error,
                     VecCopy(candidate_control, candidate_before));
  record_first_error(local_stage_error,
                     VecCopy(reference_control, reference_before));
  record_first_error(local_stage_error,
                     VecCopy(equilibrated_state, state_before));
  CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                     stage_error);
  CHKERR stage_error;
  local_stage_error = DMoFEMMeshToLocalVector(
      ep.dmIncrementalOptimization, mesh_control_before, INSERT_VALUES,
      SCATTER_FORWARD, RowColData::ROW);
  CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                     stage_error);
  CHKERR stage_error;
 
  const auto commit = [&]() {
    MoFEMFunctionBegin;
    local_stage_error = DMoFEMMeshToLocalVector(
        ep.dmIncrementalOptimization, solution, INSERT_VALUES,
        SCATTER_REVERSE, RowColData::ROW);
    CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                       stage_error);
    if (stage_error)
      MoFEMFunctionReturnHot(stage_error);
    std::map<EntityHandle, PackedPlasticCoordinates> plastic_increment;
    if (ep.plasticVolume) {
      auto save_increment = [&](boost::shared_ptr<FieldEntity> field_entity) {
        MoFEMFunctionBegin;
        auto data = field_entity->getEntFieldData();
        if (data.size() != plasticLogarithmicStretchCoordinateSize)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Plastic increment does not have five trace-free P0 "
                  "coordinates");
        PackedPlasticCoordinates values{};
        std::copy(data.begin(), data.end(), values.begin());
        plastic_increment.emplace(field_entity->getEnt(), values);
        MoFEMFunctionReturn(0);
      };
      local_stage_error =
          ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
              save_increment, ep.plasticFlowField, ep.plasticVolumes.get());
      CHKERR synchroniseCommitStageError(
          ep.mField.get_comm(), local_stage_error, stage_error);
      if (stage_error)
        MoFEMFunctionReturnHot(stage_error);
      auto add_increment = [&](boost::shared_ptr<FieldEntity> field_entity) {
        MoFEMFunctionBegin;
        auto data = field_entity->getEntFieldData();
        const auto increment = plastic_increment.find(field_entity->getEnt());
        if (increment == plastic_increment.end() ||
            data.size() != plasticLogarithmicStretchCoordinateSize)
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Plastic commit entity layout is inconsistent");
        for (int coordinate = 0;
             coordinate != plasticLogarithmicStretchCoordinateSize;
             ++coordinate)
          data[coordinate] += increment->second[coordinate];
        MoFEMFunctionReturn(0);
      };
      local_stage_error =
          ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
              add_increment, ep.plasticHField, ep.plasticVolumes.get());
      CHKERR synchroniseCommitStageError(
          ep.mField.get_comm(), local_stage_error, stage_error);
      if (stage_error)
        MoFEMFunctionReturnHot(stage_error);
 
      // The reverse scatter placed Delta kappa in the shared scalar field
      // layout. Replace it by the accumulated mesh history kappa_n+1.
      auto add_kappa_increment =
          [&](boost::shared_ptr<FieldEntity> field_entity) {
            MoFEMFunctionBegin;
            auto data = field_entity->getEntFieldData();
            const auto history = kappa_before.find(field_entity->getEnt());
            if (history == kappa_before.end() || data.size() != 1)
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                      "Plastic-kappa commit entity layout is inconsistent");
            const double delta_kappa = data[0];
            if (!std::isfinite(delta_kappa))
              SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                      "Accepted plastic-kappa increment on entity %llu must "
                      "be finite; got %g",
                      static_cast<unsigned long long>(field_entity->getEnt()),
                      delta_kappa);
            // Feasibility has already been checked with TAO's configured
            // tolerance. Clamp a tolerated negative residue so that the
            // committed history remains exactly irreversible.
            data[0] = history->second + std::max(0., delta_kappa);
            MoFEMFunctionReturn(0);
          };
      local_stage_error =
          ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
              add_kappa_increment, ep.plasticKappaField,
              ep.plasticVolumes.get());
      CHKERR synchroniseCommitStageError(
          ep.mField.get_comm(), local_stage_error, stage_error);
      if (stage_error)
        MoFEMFunctionReturnHot(stage_error);
    }
    // Both TAO blocks remain increment-only. Keep accepted Delta H_p on the
    // mesh until the step monitor writes it; mesh plasticKappa already stores
    // the accumulated value and is never cleared by the physical-step driver.
    local_stage_error = VecZeroEntries(committed_control);
    record_first_error(local_stage_error, VecZeroEntries(candidate_control));
    record_first_error(local_stage_error, VecZeroEntries(reference_control));
    CHKERR synchroniseCommitStageError(ep.mField.get_comm(), local_stage_error,
                                       stage_error);
    if (stage_error)
      MoFEMFunctionReturnHot(stage_error);
    MoFEMFunctionReturn(0);
  };
 
  const MoFEMErrorCode commit_error = commit();
  if (commit_error) {
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> t_zero(0., 0., 0., 0.,
                                                               0., 0.);
    const MoFEMErrorCode restore_plastic_error =
        ep.plasticVolume ? setPlasticHistory(ep, plastic_before, t_zero, 0.)
                         : 0;
    const MoFEMErrorCode restore_control_error =
        VecCopy(control_before, committed_control);
    const MoFEMErrorCode restore_candidate_error =
        VecCopy(candidate_before, candidate_control);
    const MoFEMErrorCode restore_reference_error =
        VecCopy(reference_before, reference_control);
    const MoFEMErrorCode restore_mesh_control_error = DMoFEMMeshToLocalVector(
        ep.dmIncrementalOptimization, mesh_control_before, INSERT_VALUES,
        SCATTER_REVERSE, RowColData::ROW);
    const MoFEMErrorCode restore_kappa_error =
        ep.plasticVolume ? setPlasticKappaHistory(ep, kappa_before) : 0;
    const MoFEMErrorCode restore_state_error =
        VecCopy(state_before, equilibrated_state);
    // Do not return between paired ghost-update calls during rollback. Every
    // rank completes the same communication sequence before errors are
    // reduced and reported below.
    const MoFEMErrorCode restore_state_ghost_begin_error =
        VecGhostUpdateBegin(equilibrated_state, INSERT_VALUES,
                            SCATTER_FORWARD);
    const MoFEMErrorCode restore_state_ghost_end_error =
        VecGhostUpdateEnd(equilibrated_state, INSERT_VALUES,
                          SCATTER_FORWARD);
    const MoFEMErrorCode restore_mesh_state_error = DMoFEMMeshToLocalVector(
        ep.dmElastic, equilibrated_state, INSERT_VALUES, SCATTER_REVERSE,
        RowColData::COL);
    MoFEMErrorCode local_restore_error = restore_plastic_error;
    if (!local_restore_error)
      local_restore_error = restore_control_error;
    if (!local_restore_error)
      local_restore_error = restore_candidate_error;
    if (!local_restore_error)
      local_restore_error = restore_reference_error;
    if (!local_restore_error)
      local_restore_error = restore_mesh_control_error;
    if (!local_restore_error)
      local_restore_error = restore_kappa_error;
    if (!local_restore_error)
      local_restore_error = restore_state_error;
    if (!local_restore_error)
      local_restore_error = restore_state_ghost_begin_error;
    if (!local_restore_error)
      local_restore_error = restore_state_ghost_end_error;
    if (!local_restore_error)
      local_restore_error = restore_mesh_state_error;
    MoFEMErrorCode restore_error = 0;
    CHKERR synchroniseCommitStageError(ep.mField.get_comm(),
                                       local_restore_error, restore_error);
    CHKERR restore_error;
    CHKERR commit_error;
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace PlasticIncrementalOptimizationInternal
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::validateSolution(
    DM constraint_dm, Vec control, Vec smooth_gradient,
    Vec inequality_multipliers, PetscReal gradient_tolerance,
    PetscReal constraint_tolerance, std::string &diagnostics) {
  return validatePlasticIncrementalSolution(
      *this, constraint_dm, control, smooth_gradient, inequality_multipliers,
      gradient_tolerance, constraint_tolerance, diagnostics);
}
 
MoFEMErrorCode PlasticIncrementalOptimizationProblem::commit(
    Vec solution, Vec equilibrated_state, Vec committed_control,
    Vec candidate_control, Vec reference_control) {
  return commitIncrementalOptimizationLoadStep(
      ep, solution, equilibrated_state, committed_control, candidate_control,
      reference_control);
}
 
MoFEMErrorCode EshelbianCore::runIncrementalOptimizationTAO(TS ts, Vec state) {
  MoFEMFunctionBegin;
  auto context = createPlasticIncrementalOptimizationContext(
      *this, SmartPetscObj<TS>(ts, true), SmartPetscObj<Vec>(state, true));
  CHKERR EshelbianPlasticity::solveIncrementalOptimizationTAO(context, state);
  MoFEMFunctionReturn(0);
}
 
} // namespace EshelbianPlasticity