PlasticIncrementalOptimizationTests.cpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalOptimizationTests.cpp
 * @brief Atom tests for incremental 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 {
 
struct OpVerifyPlasticControl : public VolUserDataOperator {
  OpVerifyPlasticControl(boost::shared_ptr<MatrixDouble> values,
                         boost::shared_ptr<VectorDouble> kappa_increments,
                         const double expected_kappa_increment,
                         boost::shared_ptr<PetscInt> element_count)
      : VolUserDataOperator(NOSPACE, OPSPACE), valuesPtr(std::move(values)),
        kappaIncrementsPtr(std::move(kappa_increments)),
        expectedKappaIncrement(expected_kappa_increment),
        elementCount(std::move(element_count)) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    const int number_of_points = getGaussPts().size2();
    if (!valuesPtr || valuesPtr->size1() != number_of_points ||
        valuesPtr->size2() != 6)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Plastic-flow reconstruction has shape %zu x %zu, expected "
              "%d x 6",
              valuesPtr ? valuesPtr->size1() : 0,
              valuesPtr ? valuesPtr->size2() : 0, number_of_points);
    if (!kappaIncrementsPtr ||
        kappaIncrementsPtr->size() != number_of_points)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Plastic-kappa increment reconstruction has size %zu, "
              "expected %d",
              kappaIncrementsPtr ? kappaIncrementsPtr->size() : 0,
              number_of_points);
    if (number_of_points <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Plastic-flow reconstruction must exercise multiple "
              "quadrature points");
 
    auto t_flow = getFTensor2SymmetricFromMat<3>(*valuesPtr);
    const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
        t_coordinates(1., 2., 3., 4., 5.);
    const auto t_expected =
        plasticLogarithmicStretchTensorFromCoordinates(t_coordinates);
    for (int gg = 0; gg != number_of_points; ++gg) {
      for (int row = 0; row != 3; ++row)
        for (int column = row; column != 3; ++column)
          if (t_flow(row, column) != t_expected(row, column))
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Trace-free plastic-flow reconstruction failed at "
                    "point %d component (%d,%d): %.17g != %.17g",
                    gg, row, column, t_flow(row, column),
                    t_expected(row, column));
      if ((*kappaIncrementsPtr)[gg] != expectedKappaIncrement)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Vector-backed plastic-kappa increment reconstruction failed "
                "at point %d: %.17g != %.17g",
                gg, (*kappaIncrementsPtr)[gg], expectedKappaIncrement);
      ++t_flow;
    }
    ++*elementCount;
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> valuesPtr;
  boost::shared_ptr<VectorDouble> kappaIncrementsPtr;
  double expectedKappaIncrement;
  boost::shared_ptr<PetscInt> elementCount;
};
 
using TestControlBlockVisitor =
    std::function<MoFEMErrorCode(PetscInt, PetscScalar *)>;
using TestOpAssemblePlasticCellMeasure =
    FormsIntegrators<VolUserDataOperator>::Assembly<A>::LinearForm<
        GAUSS>::OpSource<1, 1>;
 
MoFEMErrorCode getTestFieldIS(DM dm, const std::string &field_name,
                              SmartPetscObj<IS> &field_is) {
  MoFEMFunctionBegin;
  if (!dm)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Cannot create a field IS from a null DM");
  IS raw_is = nullptr;
  CHKERR DMMoFEMGetFieldIS(dm, RowColData::ROW, field_name.c_str(), &raw_is);
  field_is = SmartPetscObj<IS>(raw_is);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode getTestFieldBlockCounts(DM dm,
                                       const std::string &field_name,
                                       const PetscInt block_size,
                                       PetscInt &local_count,
                                       PetscInt &global_count) {
  MoFEMFunctionBegin;
  SmartPetscObj<IS> field_is;
  CHKERR getTestFieldIS(dm, field_name, field_is);
  PetscInt local_size = 0;
  PetscInt global_size = 0;
  CHKERR ISGetLocalSize(field_is, &local_size);
  CHKERR ISGetSize(field_is, &global_size);
  if (block_size <= 0 || local_size % block_size || global_size % block_size)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Invalid %s field layout: local/global sizes %d/%d, block %d",
            field_name.c_str(), static_cast<int>(local_size),
            static_cast<int>(global_size), static_cast<int>(block_size));
  local_count = local_size / block_size;
  global_count = global_size / block_size;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode forEachOwnedTestFieldBlock(
    DM dm, const std::string &field_name, const PetscInt block_size, Vec vec,
    const TestControlBlockVisitor &visitor) {
  MoFEMFunctionBegin;
  if (!vec || block_size <= 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Invalid vector or block size for field %s", field_name.c_str());
  SmartPetscObj<IS> field_is;
  CHKERR getTestFieldIS(dm, field_name, field_is);
  Vec field_vec = nullptr;
  CHKERR VecGetSubVector(vec, field_is, &field_vec);
  PetscRestoreGuard field_vec_guard(
      [&] { return VecRestoreSubVector(vec, field_is, &field_vec); });
  PetscInt local_size = 0;
  CHKERR VecGetLocalSize(field_vec, &local_size);
  if (local_size % block_size) {
    CHKERR field_vec_guard.restore();
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Field %s local size %d is not divisible by block size %d",
            field_name.c_str(), static_cast<int>(local_size),
            static_cast<int>(block_size));
  }
  PetscScalar *array = nullptr;
  CHKERR VecGetArray(field_vec, &array);
  PetscRestoreGuard array_guard(
      [&] { return VecRestoreArray(field_vec, &array); });
  MoFEMErrorCode visitor_error = 0;
  for (PetscInt block = 0; block != local_size / block_size; ++block) {
    visitor_error = visitor(block, array + block * block_size);
    if (visitor_error)
      break;
  }
  const MoFEMErrorCode restore_array_error = array_guard.restore();
  const MoFEMErrorCode restore_vector_error = field_vec_guard.restore();
  CHKERR visitor_error;
  CHKERR restore_array_error;
  CHKERR restore_vector_error;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode getOwnedTestFieldEntities(
    EshelbianCore &ep, DM dm, const std::string &field_name,
    const PetscInt block_size, std::vector<EntityHandle> &entities) {
  MoFEMFunctionBegin;
  entities.clear();
  SmartPetscObj<IS> field_is;
  CHKERR getTestFieldIS(dm, field_name, field_is);
  PetscInt local_size = 0;
  CHKERR ISGetLocalSize(field_is, &local_size);
  if (block_size <= 0 || local_size % block_size)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Invalid owned %s field layout", field_name.c_str());
  const Problem *problem_ptr = nullptr;
  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
  if (!problem_ptr || !problem_ptr->getNumeredRowDofsPtr())
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Field %s DM has no numbered row DOFs", field_name.c_str());
  const auto &dofs =
      problem_ptr->getNumeredRowDofsPtr()->get<PetscGlobalIdx_mi_tag>();
  const PetscInt *indices = nullptr;
  CHKERR ISGetIndices(field_is, &indices);
  PetscRestoreGuard indices_guard(
      [&] { return ISRestoreIndices(field_is, &indices); });
  entities.reserve(local_size / block_size);
  for (PetscInt base = 0; base != local_size; base += block_size) {
    const auto first_dof = dofs.find(indices[base]);
    if (first_dof == dofs.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Missing field %s DOF %d", field_name.c_str(),
              static_cast<int>(indices[base]));
    const EntityHandle entity = (*first_dof)->getEnt();
    entities.emplace_back(entity);
    for (PetscInt coefficient = 0; coefficient != block_size; ++coefficient) {
      const auto dof = dofs.find(indices[base + coefficient]);
      if (dof == dofs.end() || (*dof)->getEnt() != entity ||
          (*dof)->getDofOrder() != 0 ||
          (*dof)->getDofCoeffIdx() != coefficient)
        SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Field %s is not an entity-ordered P0 block layout",
                field_name.c_str());
    }
  }
  CHKERR indices_guard.restore();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode getOwnedTestKappaValues(EshelbianCore &ep, DM control_dm,
                                       std::vector<double> &values) {
  MoFEMFunctionBegin;
  std::vector<EntityHandle> entities;
  CHKERR getOwnedTestFieldEntities(ep, control_dm, ep.plasticKappaField, 1,
                                   entities);
  PlasticKappaHistory history;
  CHKERR savePlasticKappaHistory(ep, history);
  values.resize(entities.size());
  for (std::size_t block = 0; block != entities.size(); ++block) {
    const auto history_it = history.find(entities[block]);
    if (history_it == history.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Missing committed plastic kappa for entity %llu",
              static_cast<unsigned long long>(entities[block]));
    values[block] = history_it->second;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode createTestPlasticCellMeasure(EshelbianCore &ep, DM control_dm,
                                            SmartPetscObj<Vec> &measure) {
  MoFEMFunctionBegin;
  measure = createDMVector(control_dm, RowColData::ROW);
  CHKERR VecZeroEntries(measure);
  if (!ep.plasticVolume)
    MoFEMFunctionReturnHot(0);
  auto fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  fe->ksp_f = measure;
  fe->getRuleHook = [](int, int, int order_data) {
    return 2 * (order_data + 1);
  };
  MoFEM::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {NOSPACE}, ep.materialH1Positions);
  fe->getOpPtrVector().push_back(new TestOpAssemblePlasticCellMeasure(
      ep.plasticKappaField, [](double, double, double) { return 1.; }));
  fe->exeTestHook = [&ep](FEMethod *method) {
    return ep.plasticVolumes->find(method->getFEEntityHandle()) !=
           ep.plasticVolumes->end();
  };
  CHKERR DMoFEMLoopFiniteElements(control_dm, ep.elementVolumeName, fe);
  CHKERR VecAssemblyBegin(measure);
  CHKERR VecAssemblyEnd(measure);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode getTestPlasticConstraintScales(
    Vec cell_measure, IS kappa_is, std::vector<double> &scales) {
  MoFEMFunctionBegin;
  Vec measures = nullptr;
  CHKERR VecGetSubVector(cell_measure, kappa_is, &measures);
  PetscInt size = 0;
  CHKERR VecGetLocalSize(measures, &size);
  const PetscScalar *values = nullptr;
  CHKERR VecGetArrayRead(measures, &values);
  scales.resize(size);
  for (PetscInt block = 0; block != size; ++block)
    scales[block] = std::sqrt(PetscRealPart(values[block]));
  CHKERR VecRestoreArrayRead(measures, &values);
  CHKERR VecRestoreSubVector(cell_measure, kappa_is, &measures);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode setControlField(EshelbianCore &ep, const std::string &field,
                               Range entities, const std::size_t block_size) {
  MoFEMFunctionBegin;
  auto set = [&](boost::shared_ptr<FieldEntity> field_entity) {
    MoFEMFunctionBegin;
    auto data = field_entity->getEntFieldData();
    if (data.size() != block_size)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Field %s has %zu coefficients on entity %llu; expected %zu",
              field.c_str(), data.size(),
              static_cast<unsigned long long>(field_entity->getEnt()),
              block_size);
    for (std::size_t coefficient = 0; coefficient != block_size; ++coefficient)
      data[coefficient] = coefficient + 1.;
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(set, field,
                                                                    &entities);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
checkControlLayout(EshelbianCore &ep, DM control_dm,
                   const std::string &field_name, Range entities,
                   const PetscInt expected_block_size) {
  MoFEMFunctionBegin;
  ParallelComm *parallel =
      ParallelComm::get_pcomm(&ep.mField.get_moab(), MYPCOMM_INDEX);
  if (!parallel)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "ParallelComm is unavailable while checking control layout");
  Range owned_entities;
  CHKERR parallel->filter_pstatus(entities, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
                                  &owned_entities);
  PetscInt expected_local = owned_entities.size();
  PetscInt expected_global = 0;
  const int mpi_error = MPI_Allreduce(&expected_local, &expected_global, 1,
                                      MPIU_INT, MPI_SUM, ep.mField.get_comm());
  if (mpi_error != MPI_SUCCESS)
    SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
            "MPI reduction of control entity counts failed");
 
  PetscInt actual_local = 0;
  PetscInt actual_global = 0;
  CHKERR getTestFieldBlockCounts(control_dm, field_name, expected_block_size,
                                 actual_local, actual_global);
  if (actual_local != expected_local || actual_global != expected_global)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Field %s layout mismatch: block/local/global = %d/%d/%d, "
            "expected %d/%d/%d",
            field_name.c_str(), static_cast<int>(expected_block_size),
            static_cast<int>(actual_local), static_cast<int>(actual_global),
            static_cast<int>(expected_block_size),
            static_cast<int>(expected_local),
            static_cast<int>(expected_global));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkPlasticControlEntityPairing(
    EshelbianCore &ep, const PlasticIncrementalOptimizationProblem &problem) {
  MoFEMFunctionBegin;
  const Problem *dm_problem_ptr = nullptr;
  CHKERR DMMoFEMGetProblemPtr(ep.dmIncrementalOptimization, &dm_problem_ptr);
  if (!dm_problem_ptr)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Incremental-optimization DM has no problem");
 
  SmartPetscObj<IS> plastic_is;
  SmartPetscObj<IS> kappa_is;
  CHKERR getTestFieldIS(problem.getControlDM(),
                        problem.getPlasticFlowFieldName(), plastic_is);
  CHKERR getTestFieldIS(problem.getControlDM(),
                        problem.getPlasticKappaFieldName(), kappa_is);
  PetscInt plastic_size = 0;
  PetscInt epigraph_size = 0;
  CHKERR ISGetLocalSize(plastic_is, &plastic_size);
  CHKERR ISGetLocalSize(kappa_is, &epigraph_size);
  if (plastic_size != plasticLogarithmicStretchCoordinateSize * epigraph_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic and epigraph local layouts cannot be paired");
 
  const PetscInt *plastic_indices = nullptr;
  const PetscInt *epigraph_indices = nullptr;
  CHKERR ISGetIndices(plastic_is, &plastic_indices);
  CHKERR ISGetIndices(kappa_is, &epigraph_indices);
  const auto &dofs_by_global =
      dm_problem_ptr->getNumeredRowDofsPtr()->get<PetscGlobalIdx_mi_tag>();
  std::vector<double> entity_measures(epigraph_size);
  for (PetscInt block = 0; block != epigraph_size; ++block) {
    const auto epigraph_dof = dofs_by_global.find(epigraph_indices[block]);
    if (epigraph_dof == dofs_by_global.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Missing epigraph DOF %d in the control DM",
              static_cast<int>(epigraph_indices[block]));
    const EntityHandle entity = (*epigraph_dof)->getEnt();
    if ((*epigraph_dof)->getDofCoeffIdx() != 0)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Epigraph block %d does not start at coefficient zero",
              static_cast<int>(block));
    for (PetscInt coordinate = 0;
         coordinate != plasticLogarithmicStretchCoordinateSize; ++coordinate) {
      const PetscInt offset =
          block * plasticLogarithmicStretchCoordinateSize + coordinate;
      const auto plastic_dof = dofs_by_global.find(plastic_indices[offset]);
      if (plastic_dof == dofs_by_global.end())
        SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Missing plastic DOF %d in the control DM",
                static_cast<int>(plastic_indices[offset]));
      if ((*plastic_dof)->getEnt() != entity ||
          (*plastic_dof)->getDofCoeffIdx() != coordinate)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Plastic block %d and epigraph block %d refer to different "
                "entities or coefficient orderings",
                static_cast<int>(block), static_cast<int>(block));
    }
 
    const EntityHandle *connectivity = nullptr;
    int number_of_nodes = 0;
    std::vector<EntityHandle> connectivity_storage;
    CHKERR ep.mField.get_moab().get_connectivity(
        entity, connectivity, number_of_nodes, true, &connectivity_storage);
    if (number_of_nodes != 4)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Plastic P0 entity %llu is not a four-node tetrahedron",
              static_cast<unsigned long long>(entity));
    double coordinates[12];
    CHKERR ep.mField.get_moab().get_coords(connectivity, number_of_nodes,
                                           coordinates);
    const double a_x = coordinates[3] - coordinates[0];
    const double a_y = coordinates[4] - coordinates[1];
    const double a_z = coordinates[5] - coordinates[2];
    const double b_x = coordinates[6] - coordinates[0];
    const double b_y = coordinates[7] - coordinates[1];
    const double b_z = coordinates[8] - coordinates[2];
    const double c_x = coordinates[9] - coordinates[0];
    const double c_y = coordinates[10] - coordinates[1];
    const double c_z = coordinates[11] - coordinates[2];
    entity_measures[block] =
        std::abs(a_x * (b_y * c_z - b_z * c_y) - a_y * (b_x * c_z - b_z * c_x) +
                 a_z * (b_x * c_y - b_y * c_x)) /
        6.;
  }
  CHKERR ISRestoreIndices(kappa_is, &epigraph_indices);
  CHKERR ISRestoreIndices(plastic_is, &plastic_indices);
 
  SmartPetscObj<Vec> plastic_cell_measure;
  CHKERR createTestPlasticCellMeasure(ep, problem.getControlDM(),
                                      plastic_cell_measure);
  Vec measures = nullptr;
  CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &measures);
  PetscInt measure_size = 0;
  CHKERR VecGetLocalSize(measures, &measure_size);
  if (measure_size != epigraph_size) {
    CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is, &measures);
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Cell-measure and epigraph entity layouts differ");
  }
  const PetscScalar *measure_array = nullptr;
  CHKERR VecGetArrayRead(measures, &measure_array);
  double maximum_measure_error_local = 0;
  double minimum_entity_measure_local = std::numeric_limits<double>::max();
  double maximum_entity_measure_local = 0;
  for (PetscInt block = 0; block != measure_size; ++block) {
    const double assembled_measure = PetscRealPart(measure_array[block]);
    const double entity_measure = entity_measures[block];
    if (!std::isfinite(assembled_measure) || !(assembled_measure > 0) ||
        !std::isfinite(entity_measure) || !(entity_measure > 0))
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Epigraph entity block %d has invalid assembled/entity cell "
              "measure %g/%g",
              static_cast<int>(block), assembled_measure, entity_measure);
    maximum_measure_error_local =
        std::max(maximum_measure_error_local,
                 std::abs(assembled_measure - entity_measure));
    minimum_entity_measure_local =
        std::min(minimum_entity_measure_local, entity_measure);
    maximum_entity_measure_local =
        std::max(maximum_entity_measure_local, entity_measure);
  }
  CHKERR VecRestoreArrayRead(measures, &measure_array);
  CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is, &measures);
  double maximum_measure_error = 0;
  double minimum_entity_measure = 0;
  double maximum_entity_measure = 0;
  CHKERR MPI_Allreduce(&maximum_measure_error_local, &maximum_measure_error, 1,
                       MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
  CHKERR MPI_Allreduce(&minimum_entity_measure_local, &minimum_entity_measure,
                       1, MPI_DOUBLE, MPI_MIN, ep.mField.get_comm());
  CHKERR MPI_Allreduce(&maximum_entity_measure_local, &maximum_entity_measure,
                       1, MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
  const double measure_tolerance = 1e-11 * std::max(1., maximum_entity_measure);
  if (maximum_measure_error > measure_tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Epigraph measure is not paired with its tetrahedral entity: "
            "error %g > %g",
            maximum_measure_error, measure_tolerance);
  MOFEM_LOG("EP", Sev::inform)
      << "Entity-paired P0 measure range [" << minimum_entity_measure << ", "
      << maximum_entity_measure << "], maximum assembly error "
      << maximum_measure_error;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkSyntheticActiveFlowValidation(
    EshelbianCore &ep,
    const boost::shared_ptr<IncrementalOptimizationContext> &context) {
  MoFEMFunctionBegin;
  auto &problem = getPlasticIncrementalOptimizationProblem(context);
  auto control = vectorDuplicate(context->referenceControl);
  auto smooth_gradient = vectorDuplicate(context->gradient);
  auto multipliers = vectorDuplicate(context->inequalityConstraints);
  auto initial_multipliers = vectorDuplicate(context->inequalityConstraints);
  CHKERR VecZeroEntries(control);
  CHKERR VecZeroEntries(smooth_gradient);
  CHKERR VecZeroEntries(multipliers);
  CHKERR VecZeroEntries(initial_multipliers);
  CHKERR problem.initialiseInequalityMultipliers(context->constraintDM,
                                               initial_multipliers);
  SmartPetscObj<IS> plastic_is;
  SmartPetscObj<IS> kappa_is;
  CHKERR getTestFieldIS(problem.getControlDM(),
                        problem.getPlasticFlowFieldName(), plastic_is);
  CHKERR getTestFieldIS(problem.getControlDM(),
                        problem.getPlasticKappaFieldName(), kappa_is);
  SmartPetscObj<Vec> plastic_cell_measure;
  CHKERR createTestPlasticCellMeasure(ep, problem.getControlDM(),
                                      plastic_cell_measure);
  std::vector<double> committed_kappa;
  CHKERR getOwnedTestKappaValues(ep, problem.getControlDM(), committed_kappa);
 
  Vec plastic = nullptr;
  Vec epigraph = nullptr;
  Vec plastic_gradient = nullptr;
  Vec measures = nullptr;
  CHKERR VecGetSubVector(control, plastic_is, &plastic);
  CHKERR VecGetSubVector(control, kappa_is, &epigraph);
  CHKERR VecGetSubVector(smooth_gradient, plastic_is, &plastic_gradient);
  CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &measures);
  PetscInt plastic_size = 0;
  PetscInt epigraph_size = 0;
  PetscInt gradient_size = 0;
  PetscInt measure_size = 0;
  PetscInt multiplier_size = 0;
  CHKERR VecGetLocalSize(plastic, &plastic_size);
  CHKERR VecGetLocalSize(epigraph, &epigraph_size);
  CHKERR VecGetLocalSize(plastic_gradient, &gradient_size);
  CHKERR VecGetLocalSize(measures, &measure_size);
  CHKERR VecGetLocalSize(multipliers, &multiplier_size);
  if (plastic_size != plasticLogarithmicStretchCoordinateSize * epigraph_size ||
      gradient_size != plastic_size || measure_size != epigraph_size ||
      multiplier_size != epigraph_size ||
      static_cast<PetscInt>(committed_kappa.size()) != epigraph_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Synthetic active-flow vector layouts are inconsistent");
 
  PetscScalar *plastic_array = nullptr;
  PetscScalar *epigraph_array = nullptr;
  PetscScalar *gradient_array = nullptr;
  PetscScalar *multiplier_array = nullptr;
  const PetscScalar *initial_multiplier_array = nullptr;
  const PetscScalar *measure_array = nullptr;
  CHKERR VecGetArray(plastic, &plastic_array);
  CHKERR VecGetArray(epigraph, &epigraph_array);
  CHKERR VecGetArray(plastic_gradient, &gradient_array);
  CHKERR VecGetArrayRead(measures, &measure_array);
  CHKERR VecGetArray(multipliers, &multiplier_array);
  CHKERR VecGetArrayRead(initial_multipliers, &initial_multiplier_array);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_direction(.31, -.17, .23, -.29, .19);
  const double direction_norm = plasticCoordinateNorm(t_direction);
  FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
  double initial_multiplier_error_local = 0;
  for (PetscInt block = 0; block != epigraph_size; ++block) {
    auto t_flow = getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
        plastic_array + block * plasticLogarithmicStretchCoordinateSize);
    auto t_gradient =
        getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
            gradient_array + block * plasticLogarithmicStretchCoordinateSize);
    const double flow_norm = .02 * (1. + .05 * block);
    t_flow(L) = flow_norm * t_direction(L) / direction_norm;
    const double denominator =
        std::hypot(plasticEquivalentIncrementScale * flow_norm,
                   problem.getDissipationRegularizationEpsilon());
    const double alpha =
        denominator - problem.getDissipationRegularizationEpsilon();
    const double current_yield =
        problem.getInitialYieldStress() +
        problem.getIsotropicHardeningModulus() *
            (committed_kappa[block] + alpha);
    const double lambda = PetscRealPart(measure_array[block]) * current_yield;
    epigraph_array[block] = alpha;
    multiplier_array[block] =
        lambda / std::sqrt(PetscRealPart(measure_array[block]));
    const double initial_yield =
        problem.getInitialYieldStress() +
        problem.getIsotropicHardeningModulus() * committed_kappa[block];
    initial_multiplier_error_local = std::max(
        initial_multiplier_error_local,
        std::abs(std::sqrt(PetscRealPart(measure_array[block])) *
                     PetscRealPart(initial_multiplier_array[block]) -
                 PetscRealPart(measure_array[block]) * initial_yield));
    t_gradient(L) = -lambda * plasticEquivalentIncrementScaleSquared *
                    t_flow(L) / denominator;
  }
  CHKERR VecRestoreArrayRead(initial_multipliers, &initial_multiplier_array);
  CHKERR VecRestoreArray(multipliers, &multiplier_array);
  CHKERR VecRestoreArrayRead(measures, &measure_array);
  CHKERR VecRestoreArray(plastic_gradient, &gradient_array);
  CHKERR VecRestoreArray(epigraph, &epigraph_array);
  CHKERR VecRestoreArray(plastic, &plastic_array);
  CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is, &measures);
  CHKERR VecRestoreSubVector(smooth_gradient, plastic_is, &plastic_gradient);
  CHKERR VecRestoreSubVector(control, kappa_is, &epigraph);
  CHKERR VecRestoreSubVector(control, plastic_is, &plastic);
  CHKERR VecGhostUpdateBegin(control, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(control, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(smooth_gradient, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(smooth_gradient, INSERT_VALUES, SCATTER_FORWARD);
 
  constexpr PetscReal gradient_tolerance = 1e-10;
  constexpr PetscReal constraint_tolerance = 1e-10;
  double initial_multiplier_error = 0;
  CHKERR MPI_Allreduce(&initial_multiplier_error_local,
                       &initial_multiplier_error, 1, MPI_DOUBLE, MPI_MAX,
                       ep.mField.get_comm());
  if (!std::isfinite(initial_multiplier_error) ||
      initial_multiplier_error > gradient_tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Volume-scaled multiplier initialization violates eta "
            "stationarity: error %g",
            initial_multiplier_error);
 
  auto objective_gradient = vectorDuplicate(smooth_gradient);
  auto stationarity = vectorDuplicate(smooth_gradient);
  CHKERR problem.assembleObjectiveGradient(control, smooth_gradient,
                                           objective_gradient);
  CHKERR MatZeroEntries(context->inequalityJacobian);
  CHKERR problem.assembleInequalityJacobian(control, smooth_gradient,
                                            context->inequalityJacobian);
  CHKERR MatMultTranspose(context->inequalityJacobian, multipliers,
                          stationarity);
  CHKERR VecAXPY(stationarity, -1., objective_gradient);
  PetscReal stationarity_norm = 0;
  CHKERR VecNorm(stationarity, NORM_INFINITY, &stationarity_norm);
  if (!std::isfinite(stationarity_norm) ||
      stationarity_norm > gradient_tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Volume-scaled Jacobian and multipliers violate physical "
            "stationarity: residual %g",
            stationarity_norm);
 
  std::string diagnostics;
  CHKERR problem.validateSolution(context->constraintDM, control,
                                  smooth_gradient, multipliers,
                                  gradient_tolerance, constraint_tolerance,
                                  diagnostics);
  if (diagnostics.find("plastic flow active true") == std::string::npos)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Synthetic KKT point passed validation without reporting active "
            "plastic flow: %s",
            diagnostics.c_str());
 
  MOFEM_LOG("EP", Sev::inform)
      << "Synthetic active-flow validation atom test passed";
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
MoFEMErrorCode testIncrementalConstitutiveKernels() {
  MoFEMFunctionBegin;
  // Orthonormal trace-free coordinate basis in Eq. (1.52), label
  // eq:sym-basis.
  FTENSOR_INDEX(SPACE_DIM, i);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_coordinates(1., 2., 3., 4., 5.);
  const auto t_trace_free =
      plasticLogarithmicStretchTensorFromCoordinates(t_coordinates);
  const auto t_round_trip =
      plasticLogarithmicStretchCoordinatesFromTensor(t_trace_free);
  FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
  FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_round_trip_error;
  t_round_trip_error(L) = t_round_trip(L) - t_coordinates(L);
  if (std::abs(t_trace_free(i, i)) > 1e-14 ||
      std::abs(plasticFrobeniusNorm(t_trace_free) -
               plasticCoordinateNorm(t_coordinates)) > 1e-14 ||
      plasticCoordinateNorm(t_round_trip_error) > 1e-14)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Trace-free plastic-coordinate basis is inconsistent");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testIncrementalOptimizationLayout(EshelbianCore &ep, TS ts,
                                                 Vec state) {
  MoFEMFunctionBegin;
  // Verify labels eq:history-space through eq:quadrature-weight (Eqs.
  // (1.49)--(1.54)).
  CHKERR testIncrementalConstitutiveKernels();
  if (!ep.plasticVolume || ep.interfaceCrack)
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "The incremental-optimization layout test currently requires "
            "-plastic_volume 1 and -cohesive_interface_on 0");
  CHKERR setControlField(ep, ep.plasticFlowField, *ep.plasticVolumes,
                         plasticLogarithmicStretchCoordinateSize);
  CHKERR setControlField(ep, ep.plasticKappaField, *ep.plasticVolumes, 1);
 
  auto context = createPlasticIncrementalOptimizationContext(
      ep, SmartPetscObj<TS>(ts, true), SmartPetscObj<Vec>(state, true));
  auto &problem = getPlasticIncrementalOptimizationProblem(context);
  DM control_dm = problem.getControlDM();
  SmartPetscObj<IS> kappa_is;
  CHKERR getTestFieldIS(control_dm, problem.getPlasticKappaFieldName(),
                        kappa_is);
  if (!context->constraintDM)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Incremental optimization has no constraint DM");
  PetscBool square_constraint_problem = PETSC_TRUE;
  CHKERR DMMoFEMGetSquareProblem(context->constraintDM,
                                 &square_constraint_problem);
  if (square_constraint_problem)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic constraint DM must be rectangular");
 
  DM constraint_vector_dm = nullptr;
  DM constraint_matrix_dm = nullptr;
  CHKERR VecGetDM(context->inequalityConstraints, &constraint_vector_dm);
  CHKERR MatGetDM(context->inequalityJacobian, &constraint_matrix_dm);
  if (constraint_vector_dm != context->constraintDM.get() ||
      constraint_matrix_dm != context->constraintDM.get())
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic constraint vector or Jacobian lost its constraint DM");
 
  PetscInt local_constraint_blocks = 0;
  PetscInt global_constraint_blocks = 0;
  CHKERR getTestFieldBlockCounts(control_dm,
                                 problem.getPlasticKappaFieldName(), 1,
                                 local_constraint_blocks,
                                 global_constraint_blocks);
  PetscInt local_constraint_size = 0;
  PetscInt global_constraint_size = 0;
  PetscInt local_control_size = 0;
  PetscInt global_control_size = 0;
  PetscInt local_matrix_rows = 0;
  PetscInt local_matrix_columns = 0;
  PetscInt global_matrix_rows = 0;
  PetscInt global_matrix_columns = 0;
  CHKERR VecGetLocalSize(context->inequalityConstraints,
                         &local_constraint_size);
  CHKERR VecGetSize(context->inequalityConstraints, &global_constraint_size);
  CHKERR VecGetLocalSize(context->committedControl, &local_control_size);
  CHKERR VecGetSize(context->committedControl, &global_control_size);
  CHKERR MatGetLocalSize(context->inequalityJacobian, &local_matrix_rows,
                         &local_matrix_columns);
  CHKERR MatGetSize(context->inequalityJacobian, &global_matrix_rows,
                    &global_matrix_columns);
  if (local_constraint_size != local_constraint_blocks ||
      global_constraint_size != global_constraint_blocks ||
      local_matrix_rows != local_constraint_size ||
      global_matrix_rows != global_constraint_size ||
      local_matrix_columns != local_control_size ||
      global_matrix_columns != global_control_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic constraint DM vector/matrix dimensions are inconsistent");
 
  IS raw_sub_row_is = nullptr;
  CHKERR DMMoFEMGetSubRowIS(context->constraintDM, &raw_sub_row_is);
  SmartPetscObj<IS> sub_row_is(raw_sub_row_is);
  PetscBool rows_match_kappa = PETSC_FALSE;
  CHKERR ISEqual(sub_row_is.get(), kappa_is, &rows_match_kappa);
  if (!rows_match_kappa)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Constraint-DM rows do not match the plastic kappa field");
 
  IS raw_sub_col_is = nullptr;
  CHKERR DMMoFEMGetSubColIS(context->constraintDM, &raw_sub_col_is);
  SmartPetscObj<IS> sub_col_is(raw_sub_col_is);
  PetscInt local_sub_columns = 0;
  PetscInt global_sub_columns = 0;
  CHKERR ISGetLocalSize(sub_col_is, &local_sub_columns);
  CHKERR ISGetSize(sub_col_is, &global_sub_columns);
  if (local_sub_columns != local_control_size ||
      global_sub_columns != global_control_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Constraint-DM columns do not span the plastic control space");
  PetscInt control_start = 0;
  PetscInt control_end = 0;
  CHKERR VecGetOwnershipRange(context->committedControl, &control_start,
                              &control_end);
  if (control_end - control_start != local_control_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "TAO control ownership is inconsistent with its local size");
  const PetscInt *sub_column_indices = nullptr;
  CHKERR ISGetIndices(sub_col_is, &sub_column_indices);
  for (PetscInt column = 0; column != local_sub_columns; ++column)
    if (sub_column_indices[column] != control_start + column)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Constraint-DM column numbering differs from the TAO control");
  CHKERR ISRestoreIndices(sub_col_is, &sub_column_indices);
  if (context->candidateControl.get() == context->committedControl.get())
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "TAO candidate aliases the committed control vector");
 
  CHKERR checkControlLayout(ep, control_dm,
                            problem.getPlasticFlowFieldName(),
                            *ep.plasticVolumes,
                            plasticLogarithmicStretchCoordinateSize);
  CHKERR checkControlLayout(ep, control_dm,
                            problem.getPlasticKappaFieldName(),
                            *ep.plasticVolumes, 1);
  CHKERR checkPlasticControlEntityPairing(ep, problem);
 
  auto check_zero_block = [&](const PetscInt block_size, PetscInt,
                              PetscScalar *block) {
    MoFEMFunctionBegin;
    for (PetscInt component = 0; component != block_size; ++component)
      if (PetscRealPart(block[component]) != 0)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Incremental TAO control was initialized from mesh history");
    MoFEMFunctionReturn(0);
  };
  const auto check_zero_control = [&](Vec control) {
    MoFEMFunctionBegin;
    CHKERR forEachOwnedTestFieldBlock(
        control_dm, problem.getPlasticFlowFieldName(),
        plasticLogarithmicStretchCoordinateSize, control,
        [&](PetscInt block, PetscScalar *values) {
          return check_zero_block(plasticLogarithmicStretchCoordinateSize,
                                  block, values);
        });
    CHKERR forEachOwnedTestFieldBlock(
        control_dm, problem.getPlasticKappaFieldName(), 1, control,
        [&](PetscInt block, PetscScalar *values) {
          return check_zero_block(1, block, values);
        });
    MoFEMFunctionReturn(0);
  };
  CHKERR check_zero_control(context->committedControl);
  CHKERR check_zero_control(context->candidateControl);
  CHKERR check_zero_control(context->referenceControl);
  PlasticKappaHistory initial_kappa_history;
  CHKERR savePlasticKappaHistory(ep, initial_kappa_history);
  for (const auto &[entity, kappa] : initial_kappa_history)
    if (kappa != 1.)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Committed mesh kappa on entity %llu was not preserved "
              "independently of the TAO increment vector",
              static_cast<unsigned long long>(entity));
 
  constexpr double trial_delta_kappa = 0.25;
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, context->candidateControl,
      [](PetscInt, PetscScalar *values) {
        for (PetscInt coordinate = 0;
             coordinate != plasticLogarithmicStretchCoordinateSize;
             ++coordinate)
          values[coordinate] = coordinate + 1.;
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1,
      context->candidateControl,
      [trial_delta_kappa](PetscInt, PetscScalar *value) {
        value[0] = trial_delta_kappa;
        return MoFEMErrorCode(0);
      });
  CHKERR VecGhostUpdateBegin(context->candidateControl, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(context->candidateControl, INSERT_VALUES,
                           SCATTER_FORWARD);
 
  PlasticKappaHistory mesh_kappa_before_trial;
  CHKERR savePlasticKappaHistory(ep, mesh_kappa_before_trial);
  for (const auto &[entity, kappa] : mesh_kappa_before_trial)
    if (kappa != 1.)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "TAO-vector initialization changed mesh kappa on entity %llu",
              static_cast<unsigned long long>(entity));
 
  const int rule = 2 * ep.spaceOrder;
  if (ep.plasticVolume) {
    auto values = boost::make_shared<MatrixDouble>();
    auto kappa_increments = boost::make_shared<VectorDouble>();
    auto count = boost::make_shared<PetscInt>(0);
    auto fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
    auto bubble_cache = boost::make_shared<CGGUserPolynomialBase::CachePhi>(
        0, 0, MatrixDouble());
    fe->getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
        new CGGUserPolynomialBase(bubble_cache));
    fe->getRuleHook = [rule](int, int, int) { return rule; };
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        fe->getOpPtrVector(), {L2}, ep.materialH1Positions, ep.frontAdjEdges);
    CHKERR addCalculatePlasticLogarithmicStretchFieldValues(
        fe->getOpPtrVector(), ep.plasticFlowField, values, MBTET,
        ep.dmIncrementalOptimization, context->candidateControl);
    fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
        ep.dmIncrementalOptimization,
        ForcesAndSourcesCore::UserDataOperator::OPROW, ep.plasticKappaField,
        kappa_increments, context->candidateControl, MBTET));
    fe->getOpPtrVector().push_back(new OpVerifyPlasticControl(
        values, kappa_increments, trial_delta_kappa, count));
    fe->exeTestHook = [&ep](FEMethod *method) {
      return ep.plasticVolumes->find(method->getFEEntityHandle()) !=
             ep.plasticVolumes->end();
    };
    CHKERR DMoFEMLoopFiniteElements(ep.dmIncrementalOptimization,
                                    ep.elementVolumeName, fe);
    PetscInt verified_global = 0;
    const PetscInt verified_local = *count;
    if (MPI_Allreduce(&verified_local, &verified_global, 1, MPIU_INT, MPI_SUM,
                      ep.mField.get_comm()) != MPI_SUCCESS)
      SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
              "MPI reduction of reconstructed plastic cells failed");
    PetscInt unused_local = 0;
    PetscInt global_count = 0;
    CHKERR getTestFieldBlockCounts(control_dm,
                                   problem.getPlasticFlowFieldName(),
                                   plasticLogarithmicStretchCoordinateSize,
                                   unused_local, global_count);
    if (verified_global != global_count)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Plastic reconstruction visited %d blocks; expected %d",
              static_cast<int>(verified_global),
              static_cast<int>(global_count));
  }
 
  auto committed_before = vectorDuplicate(context->committedControl);
  CHKERR VecCopy(context->committedControl, committed_before);
  CHKERR context->rollbackTrial();
  PetscBool candidate_equal = PETSC_FALSE;
  PetscBool committed_equal = PETSC_FALSE;
  CHKERR VecEqual(context->candidateControl, context->committedControl,
                  &candidate_equal);
  CHKERR VecEqual(context->committedControl, committed_before,
                  &committed_equal);
  if (!candidate_equal || !committed_equal)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Trial rollback changed committed fields or failed to restore the "
            "candidate bitwise");
  CHKERR context->rollbackTrial();
  CHKERR VecEqual(context->candidateControl, context->committedControl,
                  &candidate_equal);
  if (!candidate_equal)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Post-commit rollback failed to restore the candidate bitwise");
 
  PlasticHistory plastic_history_before_commit;
  PlasticKappaHistory kappa_history_before_commit;
  CHKERR savePlasticHistory(ep, plastic_history_before_commit);
  CHKERR savePlasticKappaHistory(ep, kappa_history_before_commit);
  auto accepted_solution = vectorDuplicate(context->referenceControl);
  auto accepted_solution_before = vectorDuplicate(context->referenceControl);
  auto state_before_commit = vectorDuplicate(state);
  CHKERR VecZeroEntries(accepted_solution);
  constexpr double accepted_delta_kappa = 0.02;
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, accepted_solution,
      [accepted_delta_kappa](PetscInt, PetscScalar *values) {
        for (PetscInt coordinate = 0;
             coordinate != plasticLogarithmicStretchCoordinateSize;
             ++coordinate)
          values[coordinate] = 0.;
        values[0] = accepted_delta_kappa /
                    plasticEquivalentIncrementScale;
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1, accepted_solution,
      [accepted_delta_kappa](PetscInt, PetscScalar *value) {
        value[0] = accepted_delta_kappa;
        return MoFEMErrorCode(0);
      });
  CHKERR VecGhostUpdateBegin(accepted_solution, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(accepted_solution, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR VecCopy(accepted_solution, accepted_solution_before);
  CHKERR VecCopy(state, state_before_commit);
  CHKERR problem.commit(accepted_solution, state, context->committedControl,
                        context->candidateControl, context->referenceControl);
 
  PetscBool accepted_solution_unchanged = PETSC_FALSE;
  PetscBool state_unchanged = PETSC_FALSE;
  CHKERR VecEqual(accepted_solution, accepted_solution_before,
                  &accepted_solution_unchanged);
  CHKERR VecEqual(state, state_before_commit, &state_unchanged);
  if (!accepted_solution_unchanged || !state_unchanged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic commit modified its read-only solution or state input");
  CHKERR check_zero_control(context->committedControl);
  CHKERR check_zero_control(context->candidateControl);
  CHKERR check_zero_control(context->referenceControl);
 
  PlasticHistory plastic_history_after_commit;
  PlasticKappaHistory kappa_history_after_commit;
  CHKERR savePlasticHistory(ep, plastic_history_after_commit);
  CHKERR savePlasticKappaHistory(ep, kappa_history_after_commit);
  for (const auto &[entity, before] : plastic_history_before_commit) {
    const auto after = plastic_history_after_commit.find(entity);
    if (after == plastic_history_after_commit.end())
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Committed plastic-H entity %llu disappeared",
              static_cast<unsigned long long>(entity));
    for (PetscInt coordinate = 0;
         coordinate != plasticLogarithmicStretchCoordinateSize;
         ++coordinate) {
      const double increment =
          coordinate == 0
              ? accepted_delta_kappa / plasticEquivalentIncrementScale
              : 0.;
      if (std::abs(after->second[coordinate] -
                   (before[coordinate] + increment)) > 1e-14)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Plastic-H additive commit failed on entity %llu",
                static_cast<unsigned long long>(entity));
    }
  }
  for (const auto &[entity, before] : kappa_history_before_commit) {
    const auto after = kappa_history_after_commit.find(entity);
    if (after == kappa_history_after_commit.end() ||
        std::abs(after->second - (before + accepted_delta_kappa)) > 1e-14)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Plastic-kappa additive commit failed on entity %llu",
              static_cast<unsigned long long>(entity));
  }
  CHKERR clearPlasticIncrementFields(ep);
  auto check_zero_flow = [&](boost::shared_ptr<FieldEntity> field_entity) {
    MoFEMFunctionBegin;
    for (const double value : field_entity->getEntFieldData())
      if (value != 0.)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Transient plastic-flow mesh data were not cleared");
    MoFEMFunctionReturn(0);
  };
  CHKERR ep.mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
      check_zero_flow, ep.plasticFlowField, ep.plasticVolumes.get());
  PlasticKappaHistory kappa_history_after_clear;
  CHKERR savePlasticKappaHistory(ep, kappa_history_after_clear);
  if (kappa_history_after_clear != kappa_history_after_commit)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Clearing transient plastic flow erased committed kappa");
 
  PetscInt global_size = 0;
  CHKERR VecGetSize(context->committedControl, &global_size);
  PetscInt local_blocks = 0;
  PetscInt global_blocks = 0;
  CHKERR getTestFieldBlockCounts(control_dm,
                                 problem.getPlasticFlowFieldName(),
                                 plasticLogarithmicStretchCoordinateSize,
                                 local_blocks, global_blocks);
  if (global_size !=
      (plasticLogarithmicStretchCoordinateSize + 1) * global_blocks)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Linear hardening added an unnecessary TAO control block");
  MOFEM_LOG("EP", Sev::inform)
      << "Incremental-optimization layout verified: global control size "
      << global_size << ", plastic blocks " << ep.plasticVolumes->size();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testEquilibratedMechanicalValue(EshelbianCore &ep, TS ts,
                                               Vec state) {
  MoFEMFunctionBegin;
  // Centred finite-difference verification of the reduced bulk Helmholtz
  // derivative.  The state is re-equilibrated at both perturbed controls, so
  // the assembled comparison includes the static equilibrium adjoint.
  CHKERR validateEquilibratedMechanicalValueScope(ep);
  if (!ep.plasticVolume || ep.interfaceCrack)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "The equilibrated-value identity test is plastic-only");
 
  double fd_epsilon = 1e-6;
  double relative_tolerance = 1e-4;
  double absolute_tolerance = 1e-8;
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-equilibrated_value_fd_epsilon", &fd_epsilon,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-equilibrated_value_gradient_rtol",
                             &relative_tolerance, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-equilibrated_value_gradient_atol",
                             &absolute_tolerance, PETSC_NULLPTR);
  if (!(fd_epsilon > 0) || !(relative_tolerance >= 0) ||
      !(absolute_tolerance >= 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_INVALID_DATA,
            "Invalid equilibrated-value derivative test tolerances");
 
  constexpr double inv_sqrt_two = 0.70710678118654752440;
  const FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction(
      inv_sqrt_two, 0., 0., -inv_sqrt_two, 0., 0.);
  PlasticHistory history;
  CHKERR savePlasticHistory(ep, history);
  auto base_state = vectorDuplicate(state);
 
  auto restore = [&]() {
    MoFEMFunctionBegin;
    CHKERR setPlasticHistory(ep, history, t_direction, 0);
    CHKERR VecCopy(base_state, state);
    CHKERR setStateOnMesh(ep, state);
    MoFEMFunctionReturn(0);
  };
  auto solve_equilibrium = [&]() {
    MoFEMFunctionBegin;
    EquilibriumSolveReport report;
    CHKERR solveEquilibriumStateTS(ep, ts, state, report);
    if (report.status != EquilibriumSolveStatus::Converged)
      SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
              "Equilibrated-value fixed-control solve did not converge "
              "(status %d, error %d)",
              static_cast<int>(report.status), report.errorCode);
    MoFEMFunctionReturn(0);
  };
  auto solve_at_scale = [&](const double scale, double &value) {
    MoFEMFunctionBegin;
    CHKERR setPlasticHistory(ep, history, t_direction, scale);
    CHKERR VecCopy(base_state, state);
    CHKERR setStateOnMesh(ep, state);
    CHKERR solve_equilibrium();
    CHKERR evaluateEquilibratedMechanicalValue(ep, ts, value);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_equilibrium();
  CHKERR VecCopy(state, base_state);
 
  auto reduced_gradient =
      createDMVector(ep.dmIncrementalOptimization, RowColData::ROW);
  CHKERR assembleReducedHelmholtzGradient(ep, ts, reduced_gradient);
  auto direction_vector = vectorDuplicate(reduced_gradient);
  CHKERR VecZeroEntries(direction_vector);
  const auto t_direction_coordinates =
      plasticLogarithmicStretchCoordinatesFromTensor(t_direction);
  CHKERR forEachOwnedTestFieldBlock(
      ep.dmIncrementalOptimization, ep.plasticFlowField,
      plasticLogarithmicStretchCoordinateSize, direction_vector,
      [&t_direction_coordinates](PetscInt, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
                values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        t_values(L) = t_direction_coordinates(L);
        return MoFEMErrorCode(0);
      });
  PetscScalar assembled_derivative_scalar = 0;
  CHKERR VecDot(reduced_gradient, direction_vector,
                &assembled_derivative_scalar);
  const double assembled_derivative =
      PetscRealPart(assembled_derivative_scalar);
  double value_plus = 0;
  double value_minus = 0;
  const MoFEMErrorCode plus_error = solve_at_scale(fd_epsilon, value_plus);
  const MoFEMErrorCode minus_error =
      plus_error ? 0 : solve_at_scale(-fd_epsilon, value_minus);
  const MoFEMErrorCode restore_error = restore();
  CHKERR plus_error;
  CHKERR minus_error;
  CHKERR restore_error;
 
  const double finite_difference =
      (value_plus - value_minus) / (2 * fd_epsilon);
  const double absolute_error =
      std::abs(finite_difference - assembled_derivative);
  const double tolerance =
      absolute_tolerance +
      relative_tolerance *
          std::max(std::abs(finite_difference), std::abs(assembled_derivative));
  MOFEM_LOG("EP", Sev::inform)
      << "Reduced Helmholtz plastic derivative: finite difference "
         "= "
      << finite_difference << ", assembled = " << assembled_derivative
      << ", absolute error = " << absolute_error
      << ", tolerance = " << tolerance;
  if (!std::isfinite(finite_difference) ||
      !std::isfinite(assembled_derivative) || absolute_error > tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Reduced Helmholtz value/adjoint-gradient identity is "
            "inconsistent: error = %g, tolerance = %g",
            absolute_error, tolerance);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testIncrementalOptimizationTransaction(EshelbianCore &ep, TS ts,
                                                      Vec state) {
  MoFEMFunctionBegin;
  // Transaction-level checks for Eq. (3.25), label gen:eq:common-tests,
  // including Eq. (1.88), label eq:smooth-gradient-check.
  auto context = createPlasticIncrementalOptimizationContext(
      ep, SmartPetscObj<TS>(ts, true), SmartPetscObj<Vec>(state, true));
  const auto &problem = getPlasticIncrementalOptimizationProblem(context);
  DM control_dm = problem.getControlDM();
 
  IncrementalObjectiveEvaluation value_only;
  CHKERR evaluateIncrementalObjective(context, context->referenceControl,
                                       value_only);
  if (value_only.stateStatus != EquilibriumSolveStatus::Converged ||
      !context->baselineValid || value_only.cacheHit ||
      value_only.stateSolveCount != 1 || value_only.smoothGradient ||
      value_only.objectiveGradient || context->cachedEvaluation.smoothGradient ||
      context->gradientAssemblySeconds != 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Objective-only transaction must equilibrate without assembling "
            "a gradient");
  auto value_state = vectorDuplicate(context->equilibratedState);
  CHKERR VecCopy(context->equilibratedState, value_state);
  PetscReal reference_objective = 0;
  IncrementalObjectiveEvaluation cached_value;
  CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
      context, context->referenceControl, reference_objective, nullptr,
      cached_value);
  if (!cached_value.cacheHit || cached_value.stateSolveCount != 1 ||
      reference_objective != value_only.dissipativeValue ||
      cached_value.smoothGradient || cached_value.objectiveGradient ||
      context->gradientAssemblySeconds != 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "TAO objective-only callback did not reuse the value cache");
 
  IncrementalObjectiveEvaluation first;
  CHKERR evaluateIncrementalObjectiveAndGradient(
      context, context->referenceControl, first);
  PetscBool state_reused = PETSC_FALSE;
  CHKERR VecEqual(value_state, context->equilibratedState, &state_reused);
  if (first.stateStatus != EquilibriumSolveStatus::Converged ||
      !first.cacheHit || first.stateSolveCount != 1 || !state_reused ||
      first.conservativeValue != value_only.conservativeValue ||
      first.dissipativeValue != value_only.dissipativeValue ||
      !first.smoothGradient)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Gradient request did not reuse its objective-only equilibrium");
  if (!std::isfinite(first.conservativeValue) ||
      !std::isfinite(first.equilibriumNorm) ||
      !std::isfinite(first.smoothGradientNorm) ||
      first.objectiveGradientNorm != std::numeric_limits<double>::max())
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Transaction returned invalid diagnostics or reported "
            "outer stationarity before the complete objective was active");
 
  auto first_gradient = vectorDuplicate(context->gradient);
  auto first_state = vectorDuplicate(context->equilibratedState);
  CHKERR VecCopy(context->gradient, first_gradient);
  CHKERR VecCopy(context->equilibratedState, first_state);
  const auto first_gradient_seconds = context->gradientAssemblySeconds;
  if (ep.plasticVolume) {
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction(
        0.31, -0.17, 0.23, -0.24, 0.19, -0.07);
    const auto t_direction_coordinates =
        plasticLogarithmicStretchCoordinatesFromTensor(t_direction);
    auto direction_vector = vectorDuplicate(context->committedControl);
    CHKERR VecZeroEntries(direction_vector);
    CHKERR forEachOwnedTestFieldBlock(
        control_dm, problem.getPlasticFlowFieldName(),
        plasticLogarithmicStretchCoordinateSize, direction_vector,
        [&t_direction_coordinates](PetscInt, PetscScalar *values) {
          auto t_values =
              getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
                  values);
          FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
          t_values(L) = t_direction_coordinates(L);
          return MoFEMErrorCode(0);
        });
    auto fixed_state_force = vectorDuplicate(context->gradient);
    CHKERR assemblePlasticForce(ep, fixed_state_force);
    PetscScalar vector_derivative = 0;
    CHKERR VecDot(fixed_state_force, direction_vector, &vector_derivative);
    double scalar_derivative = 0;
    CHKERR evaluatePlasticConjugateWork(ep, t_direction, scalar_derivative);
    const double force_error =
        std::abs(PetscRealPart(vector_derivative) - scalar_derivative);
    const double force_tolerance =
        1e-11 * std::max({1., std::abs(PetscRealPart(vector_derivative)),
                          std::abs(scalar_derivative)});
    if (force_error > force_tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Fixed-state mixed plastic conjugate does not reproduce the "
              "exact discrete directional variation: error %g > %g",
              force_error, force_tolerance);
  }
 
  IncrementalObjectiveEvaluation cached;
  CHKERR evaluateIncrementalObjectiveAndGradient(
      context, context->referenceControl, cached);
  PetscBool cached_gradient_equal = PETSC_FALSE;
  PetscBool cached_state_equal = PETSC_FALSE;
  CHKERR VecEqual(context->gradient, first_gradient, &cached_gradient_equal);
  CHKERR VecEqual(context->equilibratedState, first_state, &cached_state_equal);
  if (!cached.cacheHit || cached.stateSolveCount != 1 ||
      cached.conservativeValue != first.conservativeValue ||
      !cached_gradient_equal || !cached_state_equal ||
      context->gradientAssemblySeconds != first_gradient_seconds)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Repeated transaction evaluation was not served deterministically "
            "from the state/value/gradient cache");
 
  auto first_objective_gradient = vectorDuplicate(context->objectiveGradient);
  PetscReal first_objective = 0;
  IncrementalObjectiveEvaluation complete;
  CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
      context, context->referenceControl, first_objective,
      first_objective_gradient, complete);
  // ALMM can replace this workspace with its augmented-Lagrangian gradient.
  CHKERR VecShift(context->objectiveGradient, 1);
  PetscReal repeated_objective = 0;
  CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
      context, context->referenceControl, repeated_objective,
      context->objectiveGradient, complete);
  PetscBool objective_gradient_restored = PETSC_FALSE;
  CHKERR VecEqual(context->objectiveGradient, first_objective_gradient,
                  &objective_gradient_restored);
  if (!complete.cacheHit || complete.stateSolveCount != 1 ||
      repeated_objective != first_objective || !objective_gradient_restored ||
      context->gradientAssemblySeconds != first_gradient_seconds)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Complete objective gradient was not restored from the cached "
            "conservative gradient after TAO modified its workspace");
 
  if (ep.plasticVolume) {
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction(
        0.31, -0.17, 0.23, -0.24, 0.19, -0.07);
    const auto t_direction_coordinates =
        plasticLogarithmicStretchCoordinatesFromTensor(t_direction);
    auto direction_vector = vectorDuplicate(context->referenceControl);
    CHKERR VecZeroEntries(direction_vector);
    CHKERR forEachOwnedTestFieldBlock(
        control_dm, problem.getPlasticFlowFieldName(),
        plasticLogarithmicStretchCoordinateSize, direction_vector,
        [&t_direction_coordinates](PetscInt, PetscScalar *values) {
          auto t_values =
              getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
                  values);
          FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
          t_values(L) = t_direction_coordinates(L);
          return MoFEMErrorCode(0);
        });
    constexpr double epsilon = 1e-5;
    auto plus_control = vectorDuplicate(context->referenceControl);
    auto minus_control = vectorDuplicate(context->referenceControl);
    CHKERR VecWAXPY(plus_control, epsilon, direction_vector,
                    context->referenceControl);
    CHKERR VecWAXPY(minus_control, -epsilon, direction_vector,
                    context->referenceControl);
    const auto state_solve = context->stateSolve;
    const auto state_count = context->stateSolveCount;
    context->stateSolve = [](Vec trial_state, EquilibriumSolveReport &report) {
      MoFEMFunctionBegin;
      CHKERR VecSet(trial_state, 17);
      report.status = EquilibriumSolveStatus::Diverged;
      report.errorCode = MOFEM_OPERATION_UNSUCCESSFUL;
      report.equilibriumNorm = 1;
      MoFEMFunctionReturn(0);
    };
    PetscReal rejected_objective = 0;
    IncrementalObjectiveEvaluation rejected;
    const auto rejection_error =
        evaluateIncrementalOptimizationTAOObjectiveAndGradient(
            context, plus_control, rejected_objective, nullptr, rejected);
    context->stateSolve = state_solve;
    CHKERR rejection_error;
    PetscBool control_rolled_back = PETSC_FALSE;
    PetscBool state_rolled_back = PETSC_FALSE;
    CHKERR VecEqual(context->candidateControl, context->committedControl,
                    &control_rolled_back);
    CHKERR VecEqual(context->equilibratedState, first_state,
                    &state_rolled_back);
    if (rejected.stateStatus != EquilibriumSolveStatus::Diverged ||
        !std::isinf(rejected_objective) || rejected_objective < 0 ||
        rejected.stateSolveCount != state_count + 1 || context->cacheValid ||
        !control_rolled_back || !state_rolled_back ||
        context->gradientAssemblySeconds != first_gradient_seconds)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Rejected objective-only trial did not restore state and "
              "controls without gradient assembly");
 
    IncrementalObjectiveEvaluation plus;
    IncrementalObjectiveEvaluation minus;
    CHKERR evaluateIncrementalObjective(context, plus_control, plus);
    CHKERR evaluateIncrementalObjective(context, minus_control, minus);
    if (plus.stateStatus != EquilibriumSolveStatus::Converged ||
        minus.stateStatus != EquilibriumSolveStatus::Converged ||
        plus.cacheHit || plus.stateSolveCount != state_count + 2 ||
        plus.smoothGradient || minus.smoothGradient ||
        context->gradientAssemblySeconds != first_gradient_seconds)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Objective-only plastic-flow trials did not converge without "
              "gradient assembly");
    PetscScalar assembled = 0;
    CHKERR VecDot(first_gradient, direction_vector, &assembled);
    const double finite_difference =
        (plus.conservativeValue - minus.conservativeValue) / (2 * epsilon);
    const double tolerance =
        2e-6 + 5e-4 * std::max(std::abs(finite_difference),
                               std::abs(PetscRealPart(assembled)));
    if (std::abs(finite_difference - PetscRealPart(assembled)) > tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Direct trial DeltaHp incremental-objective gradient finite "
              "difference "
              "failed: %g versus %g",
              finite_difference, PetscRealPart(assembled));
    IncrementalObjectiveEvaluation restored;
    CHKERR evaluateIncrementalObjectiveAndGradient(
        context, context->referenceControl, restored);
    if (restored.stateStatus != EquilibriumSolveStatus::Converged)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Could not restore reference state after plastic trials");
  }
 
  CHKERR testEquilibratedMechanicalValue(ep, ts, state);
 
  MOFEM_LOG("EP", Sev::inform)
      << "Incremental-optimization transaction/cache atom test passed";
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode testIncrementalOptimizationObjectiveDerivative(EshelbianCore &ep,
                                                              TS ts,
                                                              Vec state) {
  MoFEMFunctionBegin;
  PlasticHistory committed_history;
  PlasticKappaHistory committed_kappa_history;
  CHKERR savePlasticHistory(ep, committed_history);
  CHKERR savePlasticKappaHistory(ep, committed_kappa_history);
  auto original_state = vectorDuplicate(state);
  CHKERR VecCopy(state, original_state);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_history_coordinates(.13, -.29, .17, .31, -.11);
  const auto t_history =
      plasticLogarithmicStretchTensorFromCoordinates(t_history_coordinates);
  CHKERR setPlasticHistory(ep, committed_history, t_history, 2e-3);
  CHKERR setControlField(ep, ep.plasticKappaField, *ep.plasticVolumes, 1);
  auto context = createPlasticIncrementalOptimizationContext(
      ep, SmartPetscObj<TS>(ts, true), SmartPetscObj<Vec>(state, true));
  const auto &problem = getPlasticIncrementalOptimizationProblem(context);
  DM control_dm = problem.getControlDM();
  SmartPetscObj<IS> kappa_is;
  CHKERR getTestFieldIS(control_dm, problem.getPlasticKappaFieldName(),
                        kappa_is);
  SmartPetscObj<Vec> plastic_cell_measure;
  CHKERR createTestPlasticCellMeasure(ep, control_dm, plastic_cell_measure);
  std::vector<double> committed_kappa;
  CHKERR getOwnedTestKappaValues(ep, control_dm, committed_kappa);
 
  double epsilon = 1e-5;
  double relative_tolerance = 5e-4;
  double absolute_tolerance = 2e-6;
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_epsilon",
                             &epsilon, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_rtol",
                             &relative_tolerance, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_atol",
                             &absolute_tolerance, PETSC_NULLPTR);
 
  IncrementalObjectiveEvaluation baseline;
  CHKERR evaluateIncrementalObjectiveAndGradient(
      context, context->referenceControl, baseline);
  if (baseline.stateStatus != EquilibriumSolveStatus::Converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Objective derivative test could not establish equilibrium");
 
  auto operators_tester = ep.mField.getInterface<OperatorsTester>();
  auto centre = vectorDuplicate(context->referenceControl);
  CHKERR VecCopy(context->referenceControl, centre);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_centre(.31, -.17, .23, -.24, .19);
  const double centre_norm = plasticCoordinateNorm(t_centre);
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, centre,
      [&](PetscInt, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        t_values(L) = 1e-3 * t_centre(L) / centre_norm;
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1, centre,
      [](PetscInt, PetscScalar *value) {
        value[0] = 2e-3;
        return MoFEMErrorCode(0);
      });
  CHKERR VecGhostUpdateBegin(centre, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(centre, INSERT_VALUES, SCATTER_FORWARD);
 
  auto gradient = vectorDuplicate(context->objectiveGradient);
  PetscReal objective = 0;
  IncrementalObjectiveEvaluation centre_evaluation;
  CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
      context, centre, objective, gradient, centre_evaluation);
  if (centre_evaluation.stateStatus != EquilibriumSolveStatus::Converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Objective derivative centre did not equilibrate");
 
  Vec epigraph = nullptr;
  Vec epigraph_gradient = nullptr;
  Vec measures = nullptr;
  CHKERR VecGetSubVector(centre, kappa_is, &epigraph);
  CHKERR VecGetSubVector(gradient, kappa_is, &epigraph_gradient);
  CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &measures);
  PetscInt epigraph_size = 0;
  PetscInt gradient_size = 0;
  PetscInt measure_size = 0;
  CHKERR VecGetLocalSize(epigraph, &epigraph_size);
  CHKERR VecGetLocalSize(epigraph_gradient, &gradient_size);
  CHKERR VecGetLocalSize(measures, &measure_size);
  if (gradient_size != epigraph_size || measure_size != epigraph_size ||
      static_cast<PetscInt>(committed_kappa.size()) != epigraph_size) {
    CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is, &measures);
    CHKERR VecRestoreSubVector(gradient, kappa_is, &epigraph_gradient);
    CHKERR VecRestoreSubVector(centre, kappa_is, &epigraph);
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Objective test kappa layouts are inconsistent");
  }
  const PetscScalar *epigraph_array = nullptr;
  const PetscScalar *gradient_array = nullptr;
  const PetscScalar *measure_array = nullptr;
  CHKERR VecGetArrayRead(epigraph, &epigraph_array);
  CHKERR VecGetArrayRead(epigraph_gradient, &gradient_array);
  CHKERR VecGetArrayRead(measures, &measure_array);
  double expected_dissipation_local = 0;
  double maximum_epigraph_gradient_error_local = 0;
  for (PetscInt block = 0; block != epigraph_size; ++block) {
    const double delta_kappa = PetscRealPart(epigraph_array[block]);
    const double kappa_n = committed_kappa[block];
    const double hardening_modulus =
        problem.getIsotropicHardeningModulus();
    const double expected_gradient =
        PetscRealPart(measure_array[block]) *
        (problem.getInitialYieldStress() +
         hardening_modulus * (kappa_n + delta_kappa));
    expected_dissipation_local += PetscRealPart(measure_array[block]) *
                                  ((problem.getInitialYieldStress() +
                                    hardening_modulus * kappa_n) *
                                       delta_kappa +
                                   0.5 * hardening_modulus * delta_kappa *
                                       delta_kappa);
    maximum_epigraph_gradient_error_local = std::max(
        maximum_epigraph_gradient_error_local,
        std::abs(PetscRealPart(gradient_array[block]) - expected_gradient));
  }
  CHKERR VecRestoreArrayRead(measures, &measure_array);
  CHKERR VecRestoreArrayRead(epigraph_gradient, &gradient_array);
  CHKERR VecRestoreArrayRead(epigraph, &epigraph_array);
  CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is, &measures);
  CHKERR VecRestoreSubVector(gradient, kappa_is, &epigraph_gradient);
  CHKERR VecRestoreSubVector(centre, kappa_is, &epigraph);
  double expected_dissipation = 0;
  double maximum_epigraph_gradient_error = 0;
  CHKERR MPI_Allreduce(&expected_dissipation_local, &expected_dissipation, 1,
                       MPI_DOUBLE, MPI_SUM, ep.mField.get_comm());
  CHKERR MPI_Allreduce(&maximum_epigraph_gradient_error_local,
                       &maximum_epigraph_gradient_error, 1, MPI_DOUBLE, MPI_MAX,
                       ep.mField.get_comm());
  double assembled_dissipation = 0;
  CHKERR evaluatePlasticIncrementalResistanceValue(problem, centre,
                                                   assembled_dissipation);
  const double dissipation_tolerance =
      1e-13 * std::max(1., std::abs(expected_dissipation));
  if (std::abs(assembled_dissipation - expected_dissipation) >
          dissipation_tolerance ||
      maximum_epigraph_gradient_error > dissipation_tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Kappa-increment objective value or gradient is inconsistent "
            "with linear isotropic hardening");
 
  auto plastic_direction = vectorDuplicate(context->referenceControl);
  auto plastic_direction_two = vectorDuplicate(context->referenceControl);
  auto plastic_direction_three = vectorDuplicate(context->referenceControl);
  auto epigraph_direction = vectorDuplicate(context->referenceControl);
  auto mixed_direction = vectorDuplicate(context->referenceControl);
  CHKERR VecZeroEntries(plastic_direction);
  CHKERR VecZeroEntries(plastic_direction_two);
  CHKERR VecZeroEntries(plastic_direction_three);
  CHKERR VecZeroEntries(epigraph_direction);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_plastic_direction(.37, -.19, .11, .29, -.23);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_plastic_direction_two(-.07, .43, -.31, .13, .29);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_plastic_direction_three(.19, .11, -.41, .37, -.17);
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, plastic_direction,
      [&t_plastic_direction](PetscInt block, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        t_values(L) = (1. + .125 * block) * t_plastic_direction(L);
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, plastic_direction_two,
      [&t_plastic_direction_two](PetscInt block, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        const double sign = block % 2 ? -1. : 1.;
        t_values(L) = sign * (1. + .075 * block) * t_plastic_direction_two(L);
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, plastic_direction_three,
      [&t_plastic_direction_three](PetscInt block, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        const double scale[] = {-1.25, .5, 1.75};
        t_values(L) = scale[block % 3] * t_plastic_direction_three(L);
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1, epigraph_direction,
      [](PetscInt block, PetscScalar *value) {
        value[0] = block % 2 ? -.75 : 1.;
        return MoFEMErrorCode(0);
      });
  PetscReal plastic_direction_norm = 0;
  PetscReal plastic_direction_two_norm = 0;
  PetscReal plastic_direction_three_norm = 0;
  PetscReal epigraph_direction_norm = 0;
  CHKERR VecNormalize(plastic_direction, &plastic_direction_norm);
  CHKERR VecNormalize(plastic_direction_two, &plastic_direction_two_norm);
  CHKERR VecNormalize(plastic_direction_three, &plastic_direction_three_norm);
  CHKERR VecNormalize(epigraph_direction, &epigraph_direction_norm);
  if (!(plastic_direction_norm > 0) || !(plastic_direction_two_norm > 0) ||
      !(plastic_direction_three_norm > 0) || !(epigraph_direction_norm > 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Deterministic objective direction is zero");
  PetscScalar dot_one_two = 0;
  PetscScalar dot_one_three = 0;
  PetscScalar dot_two_three = 0;
  CHKERR VecDot(plastic_direction, plastic_direction_two, &dot_one_two);
  CHKERR VecDot(plastic_direction, plastic_direction_three, &dot_one_three);
  CHKERR VecDot(plastic_direction_two, plastic_direction_three, &dot_two_three);
  const double a = PetscRealPart(dot_one_two);
  const double b = PetscRealPart(dot_one_three);
  const double c = PetscRealPart(dot_two_three);
  const double gram_determinant = 1. + 2. * a * b * c - a * a - b * b - c * c;
  if (!(gram_determinant > 1e-8))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic objective directions are not linearly independent: "
            "Gram determinant %g",
            gram_determinant);
  CHKERR VecWAXPY(mixed_direction, 1., epigraph_direction, plastic_direction);
  PetscReal mixed_direction_norm = 0;
  CHKERR VecNormalize(mixed_direction, &mixed_direction_norm);
  if (!(mixed_direction_norm > 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Mixed objective direction is zero");
 
  auto reset_state = [&]() {
    MoFEMFunctionBegin;
    context->cacheValid = false;
    CHKERR VecCopy(context->baselineState, context->lastValidState);
    CHKERR VecCopy(context->baselineState, context->equilibratedState);
    CHKERR setStateOnMesh(ep, context->equilibratedState);
    MoFEMFunctionReturn(0);
  };
  const OperatorsTester::ScalarFunction evaluate = [&](Vec control,
                                                       PetscReal &value) {
    MoFEMFunctionBegin;
    CHKERR reset_state();
    const auto gradient_seconds = context->gradientAssemblySeconds;
    IncrementalObjectiveEvaluation evaluation;
    CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
        context, control, value, nullptr, evaluation);
    if (evaluation.stateStatus != EquilibriumSolveStatus::Converged ||
        evaluation.smoothGradient || evaluation.objectiveGradient ||
        context->gradientAssemblySeconds != gradient_seconds)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Perturbed objective did not equilibrate without gradient "
              "assembly");
    MoFEMFunctionReturn(0);
  };
  double maximum_error = 0;
  auto check_direction = [&](const char *name, Vec direction) {
    MoFEMFunctionBegin;
    PetscReal finite_difference = 0;
    PetscReal assembled = 0;
    CHKERR operators_tester->checkScalarCentralFiniteDifference(
        centre, direction, gradient, epsilon, evaluate, finite_difference,
        assembled);
    const double error = std::abs(finite_difference - assembled);
    const double tolerance =
        absolute_tolerance +
        relative_tolerance *
            std::max(std::abs(finite_difference), std::abs(assembled));
    if (!std::isfinite(finite_difference) || !std::isfinite(assembled) ||
        error > tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "%s objective derivative mismatch: error %g > %g", name, error,
              tolerance);
    maximum_error = std::max(maximum_error, error);
    MOFEM_LOG("EP", Sev::inform)
        << "OperatorsTester " << name
        << " objective derivative: finite difference " << finite_difference
        << ", assembled " << assembled << ", error " << error;
    MoFEMFunctionReturn(0);
  };
  CHKERR check_direction("plastic-only", plastic_direction);
  CHKERR check_direction("plastic-only-two", plastic_direction_two);
  CHKERR check_direction("plastic-only-three", plastic_direction_three);
  CHKERR check_direction("epigraph-only", epigraph_direction);
  CHKERR check_direction("mixed", mixed_direction);
 
  const FTensor::Tensor2_symmetric<double, SPACE_DIM> t_zero(0., 0., 0., 0., 0.,
                                                             0.);
  CHKERR setPlasticHistory(ep, committed_history, t_zero, 0.);
  CHKERR setPlasticKappaHistory(ep, committed_kappa_history);
  CHKERR VecCopy(original_state, state);
  CHKERR setStateOnMesh(ep, state);
  MOFEM_LOG("EP", Sev::inform)
      << "Maximum objective directional-derivative error " << maximum_error;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkPlasticConstraintSubsetAssembly(
    EshelbianCore &ep, PlasticIncrementalOptimizationProblem &problem,
    Vec control, Vec smooth_gradient) {
  MoFEMFunctionBegin;
  auto coupled_dm = createDM(ep.mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(coupled_dm, ep.dmIncrementalOptimization,
                            "INCREMENTAL_OPTIMIZATION_COUPLED_TEST");
  CHKERR DMMoFEMSetDestroyProblem(coupled_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(coupled_dm, PETSC_FALSE);
  // The flow rows stand in for another dissipative process.  Plasticity must
  // populate only its kappa-row subset of this larger constraint system.
  CHKERR DMMoFEMAddSubFieldRow(coupled_dm, ep.plasticFlowField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddSubFieldRow(coupled_dm, ep.plasticKappaField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddSubFieldCol(coupled_dm, ep.plasticFlowField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddSubFieldCol(coupled_dm, ep.plasticKappaField,
                               ep.plasticVolumes);
  CHKERR DMMoFEMAddElement(coupled_dm, ep.elementVolumeName);
  CHKERR DMSetUp(coupled_dm);
 
  auto coupled_constraints = createDMVector(coupled_dm, RowColData::ROW);
  auto coupled_jacobian = createDMMatrix(coupled_dm);
  CHKERR VecZeroEntries(coupled_constraints);
  CHKERR evaluatePlasticInequalityConstraints(problem, control,
                                               coupled_constraints);
  CHKERR MatZeroEntries(coupled_jacobian);
  CHKERR assemblePlasticInequalityJacobian(
      problem, control, smooth_gradient, coupled_jacobian);
 
  IS raw_extra_row_is = nullptr;
  CHKERR DMMoFEMGetFieldIS(coupled_dm, RowColData::ROW,
                           ep.plasticFlowField.c_str(), &raw_extra_row_is);
  SmartPetscObj<IS> extra_row_is(raw_extra_row_is);
  Vec extra_constraints = nullptr;
  CHKERR VecGetSubVector(coupled_constraints, extra_row_is,
                         &extra_constraints);
  PetscReal extra_constraint_norm = 0;
  CHKERR VecNorm(extra_constraints, NORM_INFINITY, &extra_constraint_norm);
  CHKERR VecRestoreSubVector(coupled_constraints, extra_row_is,
                             &extra_constraints);
 
  auto unit_control = vectorDuplicate(control);
  CHKERR VecSet(unit_control, 1);
  CHKERR VecGhostUpdateBegin(unit_control, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(unit_control, INSERT_VALUES, SCATTER_FORWARD);
  auto coupled_action = createDMVector(coupled_dm, RowColData::ROW);
  CHKERR MatMult(coupled_jacobian, unit_control, coupled_action);
  Vec extra_action = nullptr;
  CHKERR VecGetSubVector(coupled_action, extra_row_is, &extra_action);
  PetscReal extra_action_norm = 0;
  CHKERR VecNorm(extra_action, NORM_INFINITY, &extra_action_norm);
  CHKERR VecRestoreSubVector(coupled_action, extra_row_is, &extra_action);
  if (extra_constraint_norm > 1e-14 || extra_action_norm > 1e-14)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic contribution modified non-plastic constraint rows");
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
testIncrementalOptimizationConstraintDerivative(EshelbianCore &ep, TS ts,
                                                Vec state) {
  MoFEMFunctionBegin;
  auto context = createPlasticIncrementalOptimizationContext(
      ep, SmartPetscObj<TS>(ts, true), SmartPetscObj<Vec>(state, true));
  auto &problem = getPlasticIncrementalOptimizationProblem(context);
  PetscBool active_validation_only = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_active_validation_only",
                             &active_validation_only, PETSC_NULLPTR);
  if (active_validation_only) {
    CHKERR checkSyntheticActiveFlowValidation(ep, context);
    MoFEMFunctionReturnHot(0);
  }
  DM control_dm = problem.getControlDM();
  SmartPetscObj<IS> plastic_is;
  SmartPetscObj<IS> kappa_is;
  CHKERR getTestFieldIS(control_dm, problem.getPlasticFlowFieldName(),
                        plastic_is);
  CHKERR getTestFieldIS(control_dm, problem.getPlasticKappaFieldName(),
                        kappa_is);
  SmartPetscObj<Vec> plastic_cell_measure;
  CHKERR createTestPlasticCellMeasure(ep, control_dm, plastic_cell_measure);
  const double regularization_epsilon =
      problem.getDissipationRegularizationEpsilon();
  std::vector<double> constraint_scales;
  CHKERR getTestPlasticConstraintScales(plastic_cell_measure, kappa_is,
                                       constraint_scales);
  double epsilon = 1e-7;
  double relative_tolerance = 1e-7;
  double absolute_tolerance = 1e-10;
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_epsilon",
                             &epsilon, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_rtol",
                             &relative_tolerance, PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-incremental_optimization_derivative_atol",
                             &absolute_tolerance, PETSC_NULLPTR);
 
  CHKERR evaluateIncrementalOptimizationInequalityConstraints(
      context, context->referenceControl, context->inequalityConstraints);
  CHKERR assembleIncrementalOptimizationInequalityJacobian(
      context, context->referenceControl, context->inequalityJacobian);
  if (context->stateSolveCount || context->gradientAssemblySeconds != 0 ||
      context->cachedEvaluation.smoothGradient)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Independent plastic constraints requested equilibrium or "
            "adjoint evaluation");
 
  auto baseline_objective_gradient =
      vectorDuplicate(context->objectiveGradient);
  PetscReal baseline_objective = 0;
  IncrementalObjectiveEvaluation baseline;
  CHKERR evaluateIncrementalOptimizationTAOObjectiveAndGradient(
      context, context->referenceControl, baseline_objective,
      baseline_objective_gradient, baseline);
  if (baseline.stateStatus != EquilibriumSolveStatus::Converged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Constraint derivative test could not establish equilibrium");
  auto preserved_baseline_gradient = vectorDuplicate(baseline.smoothGradient);
  CHKERR VecCopy(baseline.smoothGradient, preserved_baseline_gradient);
  baseline.smoothGradient = preserved_baseline_gradient;
 
  Vec baseline_gradient = nullptr;
  Vec baseline_measures = nullptr;
  CHKERR VecGetSubVector(baseline.smoothGradient, plastic_is,
                         &baseline_gradient);
  CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &baseline_measures);
  const PetscScalar *baseline_gradient_array = nullptr;
  const PetscScalar *baseline_measure_array = nullptr;
  PetscInt baseline_measure_size = 0;
  CHKERR VecGetLocalSize(baseline_measures, &baseline_measure_size);
  CHKERR VecGetArrayRead(baseline_gradient, &baseline_gradient_array);
  CHKERR VecGetArrayRead(baseline_measures, &baseline_measure_array);
  double maximum_baseline_q_local = 0;
  for (PetscInt block = 0; block != baseline_measure_size; ++block) {
    const double measure = PetscRealPart(baseline_measure_array[block]);
    const auto t_coefficient_force =
        getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
            baseline_gradient_array +
            block * plasticLogarithmicStretchCoordinateSize);
    FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize> t_force;
    FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
    t_force(L) = -t_coefficient_force(L) / measure;
    maximum_baseline_q_local =
        std::max(maximum_baseline_q_local,
                 plasticYieldFunctionScale * plasticCoordinateNorm(t_force));
  }
  CHKERR VecRestoreArrayRead(baseline_measures, &baseline_measure_array);
  CHKERR VecRestoreArrayRead(baseline_gradient, &baseline_gradient_array);
  CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is,
                             &baseline_measures);
  CHKERR VecRestoreSubVector(baseline.smoothGradient, plastic_is,
                             &baseline_gradient);
  double maximum_baseline_q = 0;
  CHKERR MPI_Allreduce(&maximum_baseline_q_local, &maximum_baseline_q, 1,
                       MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
  MOFEM_LOG("EP", Sev::inform)
      << "Constraint derivative baseline maximum q " << maximum_baseline_q;
 
  if (regularization_epsilon == 0) {
    // The evaluation owns aliases into the context. Preserve its matching
    // smooth gradient, overwrite the context workspace, and prove that the
    // adapter still assembles the force-selected exact-apex row from the
    // explicitly supplied objective transaction.
    CHKERR VecZeroEntries(context->gradient);
    CHKERR assembleIncrementalOptimizationInequalityJacobianFromEvaluation(
        context, context->referenceControl, baseline,
        context->inequalityJacobian);
    Vec apex_gradient = nullptr;
    Vec apex_measures = nullptr;
    CHKERR VecGetSubVector(baseline.smoothGradient, plastic_is,
                           &apex_gradient);
    CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &apex_measures);
    PetscInt apex_measure_size = 0;
    CHKERR VecGetLocalSize(apex_measures, &apex_measure_size);
    const PetscScalar *apex_gradient_array = nullptr;
    const PetscScalar *apex_measure_array = nullptr;
    const PetscInt *apex_plastic_indices = nullptr;
    const PetscInt *apex_epigraph_indices = nullptr;
    CHKERR VecGetArrayRead(apex_gradient, &apex_gradient_array);
    CHKERR VecGetArrayRead(apex_measures, &apex_measure_array);
    CHKERR ISGetIndices(plastic_is, &apex_plastic_indices);
    CHKERR ISGetIndices(kappa_is, &apex_epigraph_indices);
    PetscInt apex_row_start = 0;
    PetscInt apex_row_end = 0;
    CHKERR MatGetOwnershipRange(context->inequalityJacobian, &apex_row_start,
                                &apex_row_end);
    if (apex_row_end - apex_row_start != apex_measure_size)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Objective-evaluation apex rows do not match cell measures");
    double adapter_error_local = 0;
    for (PetscInt block = 0; block != apex_measure_size; ++block) {
      const double measure = PetscRealPart(apex_measure_array[block]);
      const auto t_coefficient_force =
          getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
              apex_gradient_array +
              block * plasticLogarithmicStretchCoordinateSize);
      FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize> t_force;
      FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
      t_force(L) = -t_coefficient_force(L) / measure;
      const double q =
          plasticYieldFunctionScale * plasticCoordinateNorm(t_force);
      const double denominator = std::max(problem.getInitialYieldStress(), q);
      PetscScalar row_values[plasticLogarithmicStretchCoordinateSize];
      const PetscInt row = apex_row_start + block;
      CHKERR MatGetValues(context->inequalityJacobian, 1, &row,
                          plasticLogarithmicStretchCoordinateSize,
                          apex_plastic_indices +
                              block * plasticLogarithmicStretchCoordinateSize,
                          row_values);
      for (PetscInt coordinate = 0;
           coordinate != plasticLogarithmicStretchCoordinateSize; ++coordinate)
        adapter_error_local =
            std::max(adapter_error_local,
                     std::abs(PetscRealPart(row_values[coordinate]) +
                              constraint_scales[block] *
                                  t_force(coordinate) / denominator));
      PetscScalar epigraph_value = 0;
      CHKERR MatGetValues(context->inequalityJacobian, 1, &row, 1,
                          apex_epigraph_indices + block, &epigraph_value);
      adapter_error_local = std::max(
          adapter_error_local,
          std::abs(PetscRealPart(epigraph_value) - constraint_scales[block]));
    }
    CHKERR ISRestoreIndices(kappa_is, &apex_epigraph_indices);
    CHKERR ISRestoreIndices(plastic_is, &apex_plastic_indices);
    CHKERR VecRestoreArrayRead(apex_measures, &apex_measure_array);
    CHKERR VecRestoreArrayRead(apex_gradient, &apex_gradient_array);
    CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is,
                               &apex_measures);
    CHKERR VecRestoreSubVector(baseline.smoothGradient, plastic_is,
                               &apex_gradient);
    double adapter_error = 0;
    CHKERR MPI_Allreduce(&adapter_error_local, &adapter_error, 1, MPI_DOUBLE,
                         MPI_MAX, ep.mField.get_comm());
    if (adapter_error > 1e-12)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Same-control objective-evaluation apex adapter row error %g",
              adapter_error);
    CHKERR VecCopy(preserved_baseline_gradient, context->gradient);
  }
 
  auto operators_tester = ep.mField.getInterface<OperatorsTester>();
  auto centre = vectorDuplicate(context->referenceControl);
  CHKERR VecCopy(context->referenceControl, centre);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_centre(.31, -.17, .23, -.24, .19);
  const double centre_norm = plasticCoordinateNorm(t_centre);
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, centre,
      [&](PetscInt block, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        t_values(L) = 1e-2 * (1. + .125 * block) * t_centre(L) / centre_norm;
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1, centre,
      [](PetscInt block, PetscScalar *value) {
        value[0] = 1.5e-2 + 5e-4 * block;
        return MoFEMErrorCode(0);
      });
  CHKERR VecGhostUpdateBegin(centre, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(centre, INSERT_VALUES, SCATTER_FORWARD);
 
  const auto baseline_state_count = context->stateSolveCount;
  const auto baseline_gradient_seconds = context->gradientAssemblySeconds;
  CHKERR assembleIncrementalOptimizationInequalityJacobian(
      context, centre, context->inequalityJacobian);
  CHKERR evaluateIncrementalOptimizationInequalityConstraints(
      context, centre, context->inequalityConstraints);
  PetscBool smooth_gradient_unchanged = PETSC_FALSE;
  PetscBool objective_gradient_unchanged = PETSC_FALSE;
  CHKERR VecEqual(context->gradient, preserved_baseline_gradient,
                  &smooth_gradient_unchanged);
  CHKERR VecEqual(context->objectiveGradient, baseline_objective_gradient,
                  &objective_gradient_unchanged);
  if (context->stateSolveCount != baseline_state_count ||
      context->gradientAssemblySeconds != baseline_gradient_seconds ||
      !smooth_gradient_unchanged || !objective_gradient_unchanged)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Plastic constraint trials recomputed or modified objective "
            "gradient data");
  CHKERR checkPlasticConstraintSubsetAssembly(ep, problem, centre,
                                               context->gradient);
 
  Vec plastic_centre = nullptr;
  Vec epigraph_centre = nullptr;
  CHKERR VecGetSubVector(centre, plastic_is, &plastic_centre);
  CHKERR VecGetSubVector(centre, kappa_is, &epigraph_centre);
  PetscInt plastic_size = 0;
  PetscInt epigraph_size = 0;
  PetscInt constraint_size = 0;
  CHKERR VecGetLocalSize(plastic_centre, &plastic_size);
  CHKERR VecGetLocalSize(epigraph_centre, &epigraph_size);
  CHKERR VecGetLocalSize(context->inequalityConstraints, &constraint_size);
  if (plastic_size != plasticLogarithmicStretchCoordinateSize * epigraph_size ||
      constraint_size != epigraph_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Constraint test layouts are inconsistent");
  const PetscScalar *plastic_array = nullptr;
  const PetscScalar *epigraph_array = nullptr;
  const PetscScalar *constraint_array = nullptr;
  const PetscInt *plastic_indices = nullptr;
  const PetscInt *epigraph_indices = nullptr;
  CHKERR VecGetArrayRead(plastic_centre, &plastic_array);
  CHKERR VecGetArrayRead(epigraph_centre, &epigraph_array);
  CHKERR VecGetArrayRead(context->inequalityConstraints, &constraint_array);
  CHKERR ISGetIndices(plastic_is, &plastic_indices);
  CHKERR ISGetIndices(kappa_is, &epigraph_indices);
  PetscInt centre_row_start = 0;
  PetscInt centre_row_end = 0;
  CHKERR MatGetOwnershipRange(context->inequalityJacobian, &centre_row_start,
                              &centre_row_end);
  if (centre_row_end - centre_row_start != epigraph_size)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Constraint matrix rows do not pair with epigraph blocks");
  double constitutive_error_local = 0;
  for (PetscInt block = 0; block != epigraph_size; ++block) {
    const PetscInt base = block * plasticLogarithmicStretchCoordinateSize;
    const auto t_increment =
        getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
            plastic_array + base);
    const double increment_norm = plasticCoordinateNorm(t_increment);
    const double scaled_norm = plasticEquivalentIncrementScale * increment_norm;
    const double expected_equivalent_increment =
        regularization_epsilon > 0
            ? std::hypot(scaled_norm, regularization_epsilon) -
                  regularization_epsilon
            : scaled_norm;
    const double expected_constraint =
        constraint_scales[block] *
        (PetscRealPart(epigraph_array[block]) - expected_equivalent_increment);
    constitutive_error_local = std::max(
        constitutive_error_local,
        std::abs(PetscRealPart(constraint_array[block]) - expected_constraint));
 
    PetscScalar plastic_values[plasticLogarithmicStretchCoordinateSize];
    const PetscInt row = centre_row_start + block;
    CHKERR MatGetValues(context->inequalityJacobian, 1, &row,
                        plasticLogarithmicStretchCoordinateSize,
                        plastic_indices + base, plastic_values);
    const double denominator = std::hypot(scaled_norm, regularization_epsilon);
    for (PetscInt coordinate = 0;
         coordinate != plasticLogarithmicStretchCoordinateSize; ++coordinate) {
      const double expected =
          -constraint_scales[block] * plasticEquivalentIncrementScaleSquared *
          t_increment(coordinate) / denominator;
      constitutive_error_local = std::max(
          constitutive_error_local,
          std::abs(PetscRealPart(plastic_values[coordinate]) - expected));
    }
    PetscScalar epigraph_value = 0;
    CHKERR MatGetValues(context->inequalityJacobian, 1, &row, 1,
                        epigraph_indices + block, &epigraph_value);
    constitutive_error_local = std::max(
        constitutive_error_local,
        std::abs(PetscRealPart(epigraph_value) - constraint_scales[block]));
  }
  CHKERR ISRestoreIndices(kappa_is, &epigraph_indices);
  CHKERR ISRestoreIndices(plastic_is, &plastic_indices);
  CHKERR VecRestoreArrayRead(context->inequalityConstraints, &constraint_array);
  CHKERR VecRestoreArrayRead(epigraph_centre, &epigraph_array);
  CHKERR VecRestoreArrayRead(plastic_centre, &plastic_array);
  CHKERR VecRestoreSubVector(centre, kappa_is, &epigraph_centre);
  CHKERR VecRestoreSubVector(centre, plastic_is, &plastic_centre);
  double constitutive_error = 0;
  CHKERR MPI_Allreduce(&constitutive_error_local, &constitutive_error, 1,
                       MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
  if (!std::isfinite(constitutive_error) || constitutive_error > 1e-12)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Independent constraint/Jacobian formula check failed: error %g",
            constitutive_error);
 
  auto plastic_direction = vectorDuplicate(context->referenceControl);
  auto epigraph_direction = vectorDuplicate(context->referenceControl);
  auto mixed_direction = vectorDuplicate(context->referenceControl);
  CHKERR VecZeroEntries(plastic_direction);
  CHKERR VecZeroEntries(epigraph_direction);
  const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
      t_direction(.41, -.13, .23, -.31, .17);
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticFlowFieldName(),
      plasticLogarithmicStretchCoordinateSize, plastic_direction,
      [&t_direction](PetscInt block, PetscScalar *values) {
        auto t_values =
            getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(values);
        FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
        t_values(L) = (1. + .1 * block) * t_direction(L);
        return MoFEMErrorCode(0);
      });
  CHKERR forEachOwnedTestFieldBlock(
      control_dm, problem.getPlasticKappaFieldName(), 1, epigraph_direction,
      [](PetscInt block, PetscScalar *value) {
        value[0] = block % 2 ? -.5 : 1.;
        return MoFEMErrorCode(0);
      });
  PetscReal plastic_direction_norm = 0;
  PetscReal epigraph_direction_norm = 0;
  CHKERR VecNormalize(plastic_direction, &plastic_direction_norm);
  CHKERR VecNormalize(epigraph_direction, &epigraph_direction_norm);
  if (!(plastic_direction_norm > 0) || !(epigraph_direction_norm > 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Deterministic constraint direction is zero");
  CHKERR VecWAXPY(mixed_direction, 1., epigraph_direction, plastic_direction);
  PetscReal mixed_direction_norm = 0;
  CHKERR VecNormalize(mixed_direction, &mixed_direction_norm);
  if (!(mixed_direction_norm > 0))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Mixed constraint direction is zero");
 
  auto error = vectorDuplicate(context->inequalityConstraints);
  const OperatorsTester::VectorFunction evaluate = [&](Vec control,
                                                       Vec constraints) {
    return evaluateIncrementalOptimizationInequalityConstraints(
        context, control, constraints);
  };
  PetscReal error_norm = 0;
  double tolerance = 0;
  auto check_direction = [&](const char *name, Vec direction) {
    MoFEMFunctionBegin;
    CHKERR operators_tester->checkVectorCentralFiniteDifference(
        centre, direction, context->inequalityConstraints,
        context->inequalityJacobian, epsilon, evaluate, error);
    PetscReal direction_error_norm = 0;
    CHKERR VecNorm(error, NORM_INFINITY, &direction_error_norm);
    auto assembled = vectorDuplicate(context->inequalityConstraints);
    CHKERR MatMult(context->inequalityJacobian, direction, assembled);
    PetscReal assembled_norm = 0;
    CHKERR VecNorm(assembled, NORM_INFINITY, &assembled_norm);
    const double direction_tolerance =
        absolute_tolerance + relative_tolerance * assembled_norm;
    if (!std::isfinite(direction_error_norm) ||
        direction_error_norm > direction_tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "%s constraint Jacobian mismatch: error %g > %g", name,
              direction_error_norm, direction_tolerance);
    error_norm = std::max(error_norm, direction_error_norm);
    tolerance = std::max(tolerance, direction_tolerance);
    MoFEMFunctionReturn(0);
  };
  CHKERR check_direction("plastic-only", plastic_direction);
  CHKERR check_direction("epigraph-only", epigraph_direction);
  CHKERR check_direction("mixed", mixed_direction);
 
  // A rounded cone has the classical zero derivative at its apex.  The exact
  // cone instead selects a subgradient from the cell force.  Exercise both
  // the subyield (KKT-cancelling) and overstress (saturated) selections.
  if (regularization_epsilon > 0) {
    CHKERR assembleIncrementalOptimizationInequalityJacobian(
        context, context->referenceControl, context->inequalityJacobian);
    const PetscInt *apex_plastic_indices = nullptr;
    const PetscInt *apex_epigraph_indices = nullptr;
    CHKERR ISGetIndices(plastic_is, &apex_plastic_indices);
    CHKERR ISGetIndices(kappa_is, &apex_epigraph_indices);
    PetscInt row_start = 0;
    PetscInt row_end = 0;
    CHKERR MatGetOwnershipRange(context->inequalityJacobian, &row_start,
                                &row_end);
    double apex_row_error_local = 0;
    for (PetscInt block = 0; block != row_end - row_start; ++block) {
      const PetscInt row = row_start + block;
      PetscScalar values[plasticLogarithmicStretchCoordinateSize];
      CHKERR MatGetValues(context->inequalityJacobian, 1, &row,
                          plasticLogarithmicStretchCoordinateSize,
                          apex_plastic_indices +
                              block * plasticLogarithmicStretchCoordinateSize,
                          values);
      for (const auto value : values)
        apex_row_error_local =
            std::max(apex_row_error_local, std::abs(PetscRealPart(value)));
      PetscScalar eta_value = 0;
      CHKERR MatGetValues(context->inequalityJacobian, 1, &row, 1,
                          apex_epigraph_indices + block, &eta_value);
      apex_row_error_local = std::max(
          apex_row_error_local,
          std::abs(PetscRealPart(eta_value) - constraint_scales[block]));
    }
    CHKERR ISRestoreIndices(kappa_is, &apex_epigraph_indices);
    CHKERR ISRestoreIndices(plastic_is, &apex_plastic_indices);
    double apex_row_error = 0;
    CHKERR MPI_Allreduce(&apex_row_error_local, &apex_row_error, 1, MPI_DOUBLE,
                         MPI_MAX, ep.mField.get_comm());
    if (apex_row_error > 1e-12)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Regularized apex row is not sqrt(V)*[0,...,0,1]: error %g",
              apex_row_error);
 
    auto apex_error = vectorDuplicate(context->inequalityConstraints);
    CHKERR operators_tester->checkVectorCentralFiniteDifference(
        context->referenceControl, mixed_direction,
        context->inequalityConstraints, context->inequalityJacobian, epsilon,
        evaluate, apex_error);
    PetscReal apex_error_norm = 0;
    CHKERR VecNorm(apex_error, NORM_INFINITY, &apex_error_norm);
    auto apex_assembled = vectorDuplicate(context->inequalityConstraints);
    CHKERR MatMult(context->inequalityJacobian, mixed_direction,
                   apex_assembled);
    PetscReal apex_assembled_norm = 0;
    CHKERR VecNorm(apex_assembled, NORM_INFINITY, &apex_assembled_norm);
    const double apex_tolerance =
        absolute_tolerance + relative_tolerance * apex_assembled_norm;
    if (!std::isfinite(apex_error_norm) || apex_error_norm > apex_tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Regularized apex constraint Jacobian mismatch: error %g > %g",
              apex_error_norm, apex_tolerance);
    error_norm = std::max(error_norm, apex_error_norm);
  } else {
    const FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
        t_force_direction(.31, -.17, .23, -.29, .19);
    const double force_direction_norm =
        plasticCoordinateNorm(t_force_direction);
    auto check_exact_apex = [&](const double yield_ratio) {
      MoFEMFunctionBegin;
      FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize> t_force;
      FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
      const double force_norm =
          yield_ratio * problem.getInitialYieldStress() /
          plasticYieldFunctionScale;
      t_force(L) = force_norm * t_force_direction(L) / force_direction_norm;
 
      auto synthetic_gradient = vectorDuplicate(context->gradient);
      CHKERR VecZeroEntries(synthetic_gradient);
      Vec apex_measures = nullptr;
      CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &apex_measures);
      const PetscScalar *apex_measure_array = nullptr;
      CHKERR VecGetArrayRead(apex_measures, &apex_measure_array);
      CHKERR forEachOwnedTestFieldBlock(
          control_dm, problem.getPlasticFlowFieldName(),
          plasticLogarithmicStretchCoordinateSize, synthetic_gradient,
          [&](PetscInt block, PetscScalar *gradient_values) {
            auto t_gradient =
                getFTensor1FromPtr<plasticLogarithmicStretchCoordinateSize>(
                    gradient_values);
            t_gradient(L) =
                -PetscRealPart(apex_measure_array[block]) * t_force(L);
            return MoFEMErrorCode(0);
          });
      CHKERR VecRestoreArrayRead(apex_measures, &apex_measure_array);
      CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is,
                                 &apex_measures);
 
      CHKERR MatZeroEntries(context->inequalityJacobian);
      CHKERR assemblePlasticInequalityJacobian(
          problem, context->referenceControl, synthetic_gradient,
          context->inequalityJacobian);
      const PetscInt *apex_plastic_indices = nullptr;
      const PetscInt *apex_epigraph_indices = nullptr;
      CHKERR ISGetIndices(plastic_is, &apex_plastic_indices);
      CHKERR ISGetIndices(kappa_is, &apex_epigraph_indices);
      CHKERR VecGetSubVector(plastic_cell_measure, kappa_is, &apex_measures);
      CHKERR VecGetArrayRead(apex_measures, &apex_measure_array);
      PetscInt row_start = 0;
      PetscInt row_end = 0;
      CHKERR MatGetOwnershipRange(context->inequalityJacobian, &row_start,
                                  &row_end);
      double selection_error_local = 0;
      double kkt_error_local = 0;
      const double q = plasticYieldFunctionScale * force_norm;
      const double slope_denominator =
          std::max(problem.getInitialYieldStress(), q);
      for (PetscInt block = 0; block != row_end - row_start; ++block) {
        const PetscInt row = row_start + block;
        PetscScalar values[plasticLogarithmicStretchCoordinateSize];
        CHKERR MatGetValues(context->inequalityJacobian, 1, &row,
                            plasticLogarithmicStretchCoordinateSize,
                            apex_plastic_indices +
                                block * plasticLogarithmicStretchCoordinateSize,
                            values);
        double kkt_norm_squared = 0;
        const double measure = PetscRealPart(apex_measure_array[block]);
        for (PetscInt coordinate = 0;
             coordinate != plasticLogarithmicStretchCoordinateSize;
             ++coordinate) {
          const double expected_row =
              -constraint_scales[block] * t_force(coordinate) /
              slope_denominator;
          selection_error_local = std::max(
              selection_error_local,
              std::abs(PetscRealPart(values[coordinate]) - expected_row));
          const double kkt_component = -measure * t_force(coordinate) -
                                       measure / constraint_scales[block] *
                                           problem.getInitialYieldStress() *
                                           PetscRealPart(values[coordinate]);
          kkt_norm_squared += kkt_component * kkt_component;
        }
        const double expected_kkt_norm =
            measure * std::max(0., yield_ratio - 1.) *
            problem.getInitialYieldStress() / plasticYieldFunctionScale;
        kkt_error_local =
            std::max(kkt_error_local,
                     std::abs(std::sqrt(kkt_norm_squared) - expected_kkt_norm));
        PetscScalar eta_value = 0;
        CHKERR MatGetValues(context->inequalityJacobian, 1, &row, 1,
                            apex_epigraph_indices + block, &eta_value);
        selection_error_local = std::max(
            selection_error_local,
            std::abs(PetscRealPart(eta_value) - constraint_scales[block]));
      }
      CHKERR VecRestoreArrayRead(apex_measures, &apex_measure_array);
      CHKERR VecRestoreSubVector(plastic_cell_measure, kappa_is,
                                 &apex_measures);
      CHKERR ISRestoreIndices(kappa_is, &apex_epigraph_indices);
      CHKERR ISRestoreIndices(plastic_is, &apex_plastic_indices);
      double selection_error = 0;
      double kkt_error = 0;
      CHKERR MPI_Allreduce(&selection_error_local, &selection_error, 1,
                           MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
      CHKERR MPI_Allreduce(&kkt_error_local, &kkt_error, 1, MPI_DOUBLE, MPI_MAX,
                           ep.mField.get_comm());
      const double apex_tolerance =
          1e-12 * std::max(1., problem.getInitialYieldStress());
      if (selection_error > apex_tolerance || kkt_error > apex_tolerance)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Exact apex force selection at q/sigma_y=%g failed: row "
                "error %g, KKT error %g",
                yield_ratio, selection_error, kkt_error);
      MoFEMFunctionReturn(0);
    };
    CHKERR check_exact_apex(.5);
    CHKERR check_exact_apex(2.);
  }
  MOFEM_LOG("EP", Sev::inform)
      << "OperatorsTester constraint derivative error " << error_norm
      << ", tolerance " << tolerance << ", dissipation epsilon "
      << regularization_epsilon;
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode EshelbianCore::solveTestIncrementalOptimizationLayout(TS ts,
                                                                     Vec state,
                                                                     int,
                                                                     double) {
  return testIncrementalOptimizationLayout(*this, ts, state);
}
 
} // namespace EshelbianPlasticity