PlasticIncrementalOptimizationOperators.cpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalOptimizationOperators.cpp
 * @brief Finite-element operators for incremental optimization
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <PlasticIncrementalOptimization.hpp>
#include <Lie.hpp>
#include <MatrixFunction.hpp>
#include <impl/PlasticIncrementalOptimizationInternal.hpp>
 
namespace EshelbianPlasticity {
 
using namespace PlasticIncrementalOptimizationInternal;
 
namespace PlasticIncrementalOptimizationInternal {
 
using MechanicalValueEle = VolumeElementForcesAndSourcesCore;
 
struct OpIntegratePlasticConjugateWork : public VolUserDataOperator {
  OpIntegratePlasticConjugateWork(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction_in,
      boost::shared_ptr<double> conjugate_work_ptr)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
        t_direction(std::move(t_direction_in)),
        conjugateWorkPtr(std::move(conjugate_work_ptr)) {
    if (!dataAtPts || !conjugateWorkPtr)
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "Plastic work operator has a null data pointer");
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override {
    MoFEMFunctionBegin;
    (void)side;
    (void)type;
    (void)data;
 
    constexpr int dim = SPACE_DIM;
    const int nb_integration_pts = getGaussPts().size2();
    auto t_p = dataAtPts->getFTensorApproxP(nb_integration_pts);
    auto t_f_e = dataAtPts->getFTensorSmallH(nb_integration_pts);
    auto t_h_p = dataAtPts->getFTensorPlasticH(nb_integration_pts);
    auto t_w = getFTensor0IntegrationWeight();
 
    FTENSOR_INDEXES(dim, i, j, k, l);
    const EigenMatrix::Fun<double> exp_fun = [](const double value) {
      return std::exp(value);
    };
 
    double local_work = 0;
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, dim> t_eigen_values;
      FTensor::Tensor2<double, dim, dim> t_eigen_vectors;
      t_eigen_vectors(i, j) = t_h_p(i, j);
      if (computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values) !=
          MB_SUCCESS)
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Failed to diagonalise the current plastic logarithmic "
                "stretch");
      const int nb_unique = getUniqNb<dim>(t_eigen_values);
      if (nb_unique < dim)
        CHKERR sortEigenVals<dim>(t_eigen_values, t_eigen_vectors);
 
      const auto t_diff_exp = EigenMatrix::getDiffMat(
          t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, nb_unique);
      FTensor::Tensor2<double, dim, dim> t_diff_f_p;
      t_diff_f_p(i, j) = t_diff_exp(i, j, k, l) * t_direction(k, l);
 
      // Eqs. (1.22) and (1.29)--(1.30), labels
      // eq:variation-Fp-from-Hp, eq:fixed-state-plastic-variation, and
      // eq:plastic-variation-conjugacy. This is the exact mixed-discrete
      // variation; the logarithmic Mandel rewrite needs continuum
      // compatibility and is retained separately as a diagnostic.
      local_work -= t_w * (t_p(i, j) * (t_f_e(i, k) * t_diff_f_p(k, j)));
 
      ++t_p;
      ++t_f_e;
      ++t_h_p;
      ++t_w;
    }
    *conjugateWorkPtr += getMeasure() * local_work;
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  FTensor::Tensor2_symmetric<double, SPACE_DIM> t_direction;
  boost::shared_ptr<double> conjugateWorkPtr;
};
 
using OpAssemblePlasticForce =
    FormsIntegrators<VolUserDataOperator>::Assembly<A>::LinearForm<
        GAUSS>::OpBaseTimesVector<1, plasticLogarithmicStretchCoordinateSize,
                                  plasticLogarithmicStretchCoordinateSize>;
 
using PlasticControlDofIndices =
    std::map<EntityHandle,
             std::array<PetscInt, plasticLogarithmicStretchCoordinateSize>>;
 
/**
 * Assemble -R_{H^p}^T p for one plastic P0 cell.
 *
 * The adjoint fields p are interpolated with the elastic DM.  Since the
 * plastic field is cellwise constant, this OPSPACE operator can integrate all
 * five coefficient directions and insert them directly into their control
 * DOFs without introducing a mixed auxiliary DM.
 */
struct OpAssembleReducedHelmholtzGradient : public VolUserDataOperator {
  OpAssembleReducedHelmholtzGradient(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<MatrixDouble> adjoint_piola_ptr,
      boost::shared_ptr<MatrixDouble> adjoint_stretch_ptr,
      boost::shared_ptr<MatrixDouble> adjoint_rotation_ptr,
      boost::shared_ptr<MatrixDouble> adjoint_rotation_gradient_ptr,
      boost::shared_ptr<const PlasticControlDofIndices> control_indices_ptr,
      Vec gradient, const double alpha_omega)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
        adjointPiolaPtr(std::move(adjoint_piola_ptr)),
        adjointStretchPtr(std::move(adjoint_stretch_ptr)),
        adjointRotationPtr(std::move(adjoint_rotation_ptr)),
        adjointRotationGradientPtr(std::move(adjoint_rotation_gradient_ptr)),
        controlIndicesPtr(std::move(control_indices_ptr)), gradient(gradient),
        alphaOmega(alpha_omega) {
    if (!dataAtPts || !adjointPiolaPtr || !adjointStretchPtr ||
        !adjointRotationPtr || !adjointRotationGradientPtr ||
        !controlIndicesPtr || !gradient)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Reduced Helmholtz-gradient operator has a null data pointer");
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    constexpr int dim = SPACE_DIM;
    constexpr int coordinate_size =
        plasticLogarithmicStretchCoordinateSize;
    constexpr int stretch_size = dim * (dim + 1) / 2;
    const int nb_integration_pts = getGaussPts().size2();
    const auto indices_it = controlIndicesPtr->find(getFEEntityHandle());
    if (indices_it == controlIndicesPtr->end())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Plastic cell %llu has no control-DOF layout",
              static_cast<unsigned long long>(getFEEntityHandle()));
 
#ifndef NDEBUG
    const auto check_matrix =
        [nb_integration_pts](const boost::shared_ptr<MatrixDouble> &matrix_ptr,
                             const int columns, const char *name) {
          MoFEMFunctionBegin;
          if (matrix_ptr->size1() != nb_integration_pts ||
              matrix_ptr->size2() != columns)
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Inconsistent %s data in reduced Helmholtz gradient",
                    name);
          MoFEMFunctionReturn(0);
        };
    CHKERR check_matrix(adjointPiolaPtr, dim * dim, "Piola adjoint");
    CHKERR check_matrix(adjointStretchPtr, stretch_size, "stretch adjoint");
    CHKERR check_matrix(adjointRotationPtr, dim, "rotation adjoint");
    CHKERR check_matrix(adjointRotationGradientPtr, dim * dim,
                        "rotation-adjoint gradient");
#endif
 
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto t_adjoint_piola =
        MatrixSizeHelper<GetFTensor2FromMatType<dim, dim, -1, DL>, DL>::get(
            *adjointPiolaPtr, nb_integration_pts)();
    auto t_adjoint_stretch =
        MatrixSizeHelper<GetFTensor1FromMatType<stretch_size, -1, DL>, DL>::get(
            *adjointStretchPtr, nb_integration_pts)();
    auto t_adjoint_rotation =
        MatrixSizeHelper<GetFTensor1FromMatType<dim, -1, DL>, DL>::get(
            *adjointRotationPtr, nb_integration_pts)();
    auto t_adjoint_rotation_gradient =
        MatrixSizeHelper<GetFTensor2FromMatType<dim, dim, -1, DL>, DL>::get(
            *adjointRotationGradientPtr, nb_integration_pts)();
    auto t_reference_p = dataAtPts->getFTensorApproxP(nb_integration_pts);
    auto t_elastic_f = dataAtPts->getFTensorSmallH(nb_integration_pts);
    auto t_plastic_h = dataAtPts->getFTensorPlasticH(nb_integration_pts);
    auto t_plastic_f = dataAtPts->getFTensorPlasticF(nb_integration_pts);
    auto t_inv_plastic_f =
        dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
    auto t_adjoint_stress_d_p =
        dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
    auto t_levi_d_p =
        dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
    auto t_rotation_gradient =
        dataAtPts->getFTensorRotAxisGrad(nb_integration_pts);
    auto t_w = getFTensor0IntegrationWeight();
 
    std::array<PetscScalar, coordinate_size> local_gradient{};
    FTENSOR_INDEXES(dim, i, j, k, l, m, n);
    FTensor::Index<'S', stretch_size> S;
    FTensor::Index<'C', coordinate_size> C;
    const EigenMatrix::Fun<double> exp_fun = [](const double value) {
      return std::exp(value);
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, dim> t_eigen_values;
      FTensor::Tensor2<double, dim, dim> t_eigen_vectors;
      t_eigen_vectors(i, j) = t_plastic_h(i, j);
      if (computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values) !=
          MB_SUCCESS)
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Failed to diagonalise plastic logarithmic stretch while "
                "assembling its reduced gradient");
      const int nb_unique = getUniqNb<dim>(t_eigen_values);
      if (nb_unique < dim)
        CHKERR sortEigenVals<dim>(t_eigen_values, t_eigen_vectors);
      const auto t_diff_exp = EigenMatrix::getDiffMat(
          t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, nb_unique);
      const double det_plastic_f = determinantTensor3by3(t_plastic_f);
 
      FTensor::Tensor2<double, dim, dim> t_adjoint_rotation_gradient_pull;
      FTensor::Tensor2<double, dim, dim> t_rotation_gradient_pull;
      t_adjoint_rotation_gradient_pull(k, j) =
          t_adjoint_rotation_gradient(k, i) * t_inv_plastic_f(i, j);
      t_rotation_gradient_pull(k, j) =
          t_rotation_gradient(k, i) * t_inv_plastic_f(i, j);
 
      for (int coordinate = 0; coordinate != coordinate_size; ++coordinate) {
        FTensor::Tensor1<double, coordinate_size> t_coordinate_direction;
        t_coordinate_direction(C) = 0.;
        t_coordinate_direction(coordinate) = 1.;
        const auto t_direction =
            plasticLogarithmicStretchTensorFromCoordinates(
                t_coordinate_direction);
        FTensor::Tensor2<double, dim, dim> t_d_plastic_f;
        t_d_plastic_f(i, j) =
            t_diff_exp(i, j, k, l) * t_direction(k, l);
        FTensor::Tensor2<double, dim, dim> t_d_inv_plastic_f;
        t_d_inv_plastic_f(i, j) =
            -t_inv_plastic_f(i, k) * t_d_plastic_f(k, l) *
            t_inv_plastic_f(l, j);
 
        // Piola/bubble compatibility block.  The divergence transform is
        // invariant because H^p and every coordinate direction are
        // trace-free, hence J_p=1 and dJ_p=0.
        double residual_derivative_contraction =
            -t_adjoint_piola(i, k) * t_elastic_f(i, j) *
            t_d_plastic_f(j, k);
 
        FTensor::Tensor2<double, dim, dim> t_d_intermediate_p;
        t_d_intermediate_p(i, j) =
            t_reference_p(i, k) * t_d_plastic_f(j, k) / det_plastic_f;
        FTensor::Tensor1<double, stretch_size> t_d_adjoint_stress;
        t_d_adjoint_stress(S) =
            t_adjoint_stress_d_p(i, j, S) * t_d_intermediate_p(i, j);
        // The explicit-stretch residual contains minus the stress conjugate
        // induced by the intermediate Piola transform.
        residual_derivative_contraction -=
            det_plastic_f *
            (t_adjoint_stretch(S) * t_d_adjoint_stress(S));
 
        FTensor::Tensor1<double, dim> t_d_levi;
        t_d_levi(k) = t_levi_d_p(k, i, j) * t_d_intermediate_p(i, j);
        residual_derivative_contraction -=
            det_plastic_f * t_adjoint_rotation(k) * t_d_levi(k);
 
        if (alphaOmega != 0.) {
          FTensor::Tensor2<double, dim, dim>
              t_d_adjoint_rotation_gradient_pull;
          FTensor::Tensor2<double, dim, dim> t_d_rotation_gradient_pull;
          t_d_adjoint_rotation_gradient_pull(k, j) =
              t_adjoint_rotation_gradient(k, i) *
              t_d_inv_plastic_f(i, j);
          t_d_rotation_gradient_pull(k, j) =
              t_rotation_gradient(k, i) * t_d_inv_plastic_f(i, j);
          residual_derivative_contraction +=
              det_plastic_f * alphaOmega *
              (t_d_adjoint_rotation_gradient_pull(k, j) *
                   t_rotation_gradient_pull(k, j) +
               t_adjoint_rotation_gradient_pull(k, j) *
                   t_d_rotation_gradient_pull(k, j));
        }
 
        // L=Psi-p.R has derivative -p.R_Hp because the bulk Helmholtz
        // energy has no explicit H^p dependence in this independent elastic-
        // stretch, trace-free discretization.
        local_gradient[coordinate] -=
            getMeasure() * t_w * residual_derivative_contraction;
      }
 
      ++t_adjoint_piola;
      ++t_adjoint_stretch;
      ++t_adjoint_rotation;
      ++t_adjoint_rotation_gradient;
      ++t_reference_p;
      ++t_elastic_f;
      ++t_plastic_h;
      ++t_plastic_f;
      ++t_inv_plastic_f;
      ++t_adjoint_stress_d_p;
      ++t_levi_d_p;
      ++t_rotation_gradient;
      ++t_w;
    }
 
    CHKERR VecSetValues(gradient, coordinate_size, indices_it->second.data(),
                        local_gradient.data(), ADD_VALUES);
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<MatrixDouble> adjointPiolaPtr;
  boost::shared_ptr<MatrixDouble> adjointStretchPtr;
  boost::shared_ptr<MatrixDouble> adjointRotationPtr;
  boost::shared_ptr<MatrixDouble> adjointRotationGradientPtr;
  boost::shared_ptr<const PlasticControlDofIndices> controlIndicesPtr;
  Vec gradient;
  const double alphaOmega;
};
 
/**
 * Produce the five trace-free coefficient forces in Gauss-by-coefficients
 * layout.
 *
 * Implements Eqs. (1.29)--(1.30), labels
 * eq:fixed-state-plastic-variation and eq:plastic-variation-conjugacy, in the
 * coefficient form of Eq. (1.69), label eq:baseline-smooth-gradient.
 */
struct OpCalculatePlasticForce : public VolUserDataOperator {
  OpCalculatePlasticForce(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          boost::shared_ptr<MatrixDouble> coefficient_force_ptr)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(std::move(data_ptr)),
        coefficientForcePtr(std::move(coefficient_force_ptr)) {
    if (!dataAtPts || !coefficientForcePtr)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "Plastic-force operator has a null integration-point data or "
          "coefficient-force pointer");
  }
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    constexpr int dim = SPACE_DIM;
    constexpr int coordinate_size = plasticLogarithmicStretchCoordinateSize;
    const int nb_integration_pts = getGaussPts().size2();
    using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
    auto get_force =
        MatrixSizeHelper<GetFTensor1FromMatType<coordinate_size, coordinate_size,
                                               DL>,
                         DL>::size(*coefficientForcePtr, nb_integration_pts);
    auto t_force = get_force();
    auto t_p = dataAtPts->getFTensorApproxP(nb_integration_pts);
    auto t_f_e = dataAtPts->getFTensorSmallH(nb_integration_pts);
    auto t_h_p = dataAtPts->getFTensorPlasticH(nb_integration_pts);
 
    FTENSOR_INDEXES(dim, i, j, k, m, n);
    FTENSOR_INDEX(coordinate_size, L);
    const EigenMatrix::Fun<double> exp_fun = [](const double value) {
      return std::exp(value);
    };
 
    for (int gg = 0; gg != nb_integration_pts; ++gg) {
      FTensor::Tensor1<double, dim> t_eigen_values;
      FTensor::Tensor2<double, dim, dim> t_eigen_vectors;
      t_eigen_vectors(i, j) = t_h_p(i, j);
      if (computeEigenValuesSymmetric(t_eigen_vectors, t_eigen_values) !=
          MB_SUCCESS)
        SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                "Failed to diagonalise the current plastic logarithmic "
                "stretch");
      const int nb_unique = getUniqNb<dim>(t_eigen_values);
      if (nb_unique < dim)
        CHKERR sortEigenVals<dim>(t_eigen_values, t_eigen_vectors);
      const auto t_diff_exp = EigenMatrix::getDiffMat(
          t_eigen_values, t_eigen_vectors, exp_fun, exp_fun, nb_unique);
 
      FTensor::Tensor2<double, dim, dim> t_elastic_pull_back;
      FTensor::Tensor2_symmetric<double, dim> t_tensor_force;
      t_elastic_pull_back(k, j) = t_f_e(i, k) * t_p(i, j);
      t_tensor_force(m, n) =
          -t_diff_exp(k, j, m, n) * t_elastic_pull_back(k, j);
      const auto t_coordinate_force =
          plasticLogarithmicStretchCoordinatesFromTensor(t_tensor_force);
      t_force(L) = t_coordinate_force(L);
 
      ++t_force;
      ++t_p;
      ++t_f_e;
      ++t_h_p;
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<MatrixDouble> coefficientForcePtr;
};
 
} // namespace PlasticIncrementalOptimizationInternal
 
MoFEMErrorCode addCalculatePlasticLogarithmicStretchFieldValues(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    const std::string &field_name,
    boost::shared_ptr<MatrixDouble> tensor_values,
    const EntityType zero_type, SmartPetscObj<DM> data_dm,
    SmartPetscObj<Vec> data_vector) {
  MoFEMFunctionBegin;
  if (!tensor_values)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Plastic logarithmic-stretch output is null");
  if (data_dm && !data_vector)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Cross-DM plastic logarithmic-stretch evaluation requires a "
            "data vector");
 
  if (data_dm)
    pipeline.push_back(new OpCalculateTensor2SymmetricDeviatorFieldValues<SPACE_DIM>(
        data_dm, ForcesAndSourcesCore::UserDataOperator::OPROW, field_name,
        tensor_values, data_vector, zero_type));
  else
    pipeline.push_back(new OpCalculateTensor2SymmetricDeviatorFieldValues<SPACE_DIM>(
        field_name, tensor_values, data_vector, zero_type));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
evaluateEquilibratedMechanicalValue(EshelbianCore &ep, TS ts,
                                    double &value) {
  MoFEMFunctionBegin;
  CHKERR validateEquilibratedMechanicalValueScope(ep);
  if (!ep.physicalEquations || !ep.dataAtPts)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Material equations or integration-point data are not set");
 
  auto local_helmholtz_free_energy = boost::make_shared<double>(0);
  auto fe = boost::make_shared<MechanicalValueEle>(ep.mField);
  fe->data_ctx |= PetscData::CTX_SET_TIME;
  if (ts) {
    CHKERR TSGetTime(ts, &fe->ts_t);
    CHKERR TSGetTimeStep(ts, &fe->ts_dt);
  }
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, fe);
  auto *helmholtz_free_energy_op =
      ep.physicalEquations->returnOpCalculateHelmholtzFreeEnergy(
          ep.dataAtPts, local_helmholtz_free_energy);
  if (!helmholtz_free_energy_op)
    SETERRQ(ep.mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
            "Selected material does not provide a Helmholtz free energy");
  fe->getOpPtrVector().push_back(helmholtz_free_energy_op);
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({fe.get()});
  const MoFEMErrorCode loop_error =
      DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({fe.get()});
  CHKERR loop_error;
 
  double global_helmholtz_free_energy = 0;
  const int mpi_error =
      MPI_Allreduce(local_helmholtz_free_energy.get(),
                    &global_helmholtz_free_energy, 1, MPI_DOUBLE, MPI_SUM,
                    ep.mField.get_comm());
  if (mpi_error != MPI_SUCCESS)
    SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
            "MPI reduction of the equilibrated mechanical value failed");
 
  // This is intentionally only the bulk material Helmholtz free energy.
  // Springs, loads, and discretization-stabilization terms belong to the
  // equilibrium constraint; they affect the reduced derivative through the
  // static adjoint but are never added to this objective value.
  value = global_helmholtz_free_energy;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode evaluatePlasticConjugateWork(
    EshelbianCore &ep,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_direction,
    double &conjugate_work) {
  MoFEMFunctionBegin;
  CHKERR validateEquilibratedMechanicalValueScope(ep);
  if (!ep.plasticVolume || !ep.plasticVolumes || ep.plasticVolumes->empty())
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic conjugate work requires a non-empty plastic volume");
 
  auto local_work = boost::make_shared<double>(0);
  auto fe = boost::make_shared<MechanicalValueEle>(ep.mField);
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, fe);
  fe->getOpPtrVector().push_back(new OpIntegratePlasticConjugateWork(
      ep.dataAtPts, t_direction, local_work));
  fe->exeTestHook = [&ep](FEMethod *fe_method_ptr) {
    return ep.plasticVolumes->find(fe_method_ptr->getFEEntityHandle()) !=
           ep.plasticVolumes->end();
  };
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({fe.get()});
  const MoFEMErrorCode loop_error =
      DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({fe.get()});
  CHKERR loop_error;
 
  const int mpi_error =
      MPI_Allreduce(local_work.get(), &conjugate_work, 1, MPI_DOUBLE, MPI_SUM,
                    ep.mField.get_comm());
  if (mpi_error != MPI_SUCCESS)
    SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
            "MPI reduction of the plastic conjugate work failed");
 
  MoFEMFunctionReturn(0);
}
 
namespace PlasticIncrementalOptimizationInternal {
 
MoFEMErrorCode assemblePlasticForce(EshelbianCore &ep, Vec force) {
  MoFEMFunctionBegin;
  // Exact coefficient gradient: Eq. (1.69), label
  // eq:baseline-smooth-gradient, from Eqs. (1.29)--(1.30).
  if (!force)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Plastic-force vector is null");
  CHKERR VecZeroEntries(force);
  if (!ep.plasticVolume)
    MoFEMFunctionReturnHot(0);
 
  auto force_at_points = boost::make_shared<MatrixDouble>();
  auto fe = boost::make_shared<MechanicalValueEle>(ep.mField);
  fe->ksp_f = force;
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, fe);
  fe->getOpPtrVector().push_back(
      new OpCalculatePlasticForce(ep.dataAtPts, force_at_points));
  fe->getOpPtrVector().push_back(
      new OpAssemblePlasticForce(ep.plasticFlowField, force_at_points));
  fe->exeTestHook = [&ep](FEMethod *method) {
    return ep.plasticVolumes->find(method->getFEEntityHandle()) !=
           ep.plasticVolumes->end();
  };
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({fe.get()});
  const MoFEMErrorCode loop_error = DMoFEMLoopFiniteElements(
      ep.dmIncrementalOptimization, ep.elementVolumeName, fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({fe.get()});
  CHKERR loop_error;
  CHKERR VecAssemblyBegin(force);
  CHKERR VecAssemblyEnd(force);
  MoFEMFunctionReturn(0);
}
 
namespace {
 
MoFEMErrorCode createPlasticControlDofIndices(
    EshelbianCore &ep,
    boost::shared_ptr<PlasticControlDofIndices> &control_indices_ptr) {
  MoFEMFunctionBegin;
  control_indices_ptr = boost::make_shared<PlasticControlDofIndices>();
  const Problem *problem_ptr = nullptr;
  CHKERR DMMoFEMGetProblemPtr(ep.dmIncrementalOptimization, &problem_ptr);
  if (!problem_ptr || !problem_ptr->getNumeredRowDofsPtr())
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Incremental-optimization problem has no numbered row DOFs");
 
  const auto field_bit =
      ep.mField.get_field_bit_number(ep.plasticFlowField);
  const auto &dofs =
      problem_ptr->getNumeredRowDofsPtr()->get<Unique_mi_tag>();
  auto dof_it =
      dofs.lower_bound(FieldEntity::getLoBitNumberUId(field_bit));
  const auto dof_end =
      dofs.upper_bound(FieldEntity::getHiBitNumberUId(field_bit));
  for (; dof_it != dof_end; ++dof_it) {
    const auto &dof = *dof_it;
    if (dof->getDofOrder() != 0 ||
        dof->getDofCoeffIdx() >= plasticLogarithmicStretchCoordinateSize)
      continue;
    auto [indices_it, inserted] =
        control_indices_ptr->try_emplace(dof->getEnt());
    if (inserted)
      indices_it->second.fill(-1);
    auto &index = indices_it->second[dof->getDofCoeffIdx()];
    if (index >= 0)
      SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Duplicate plastic control coefficient %d on entity %llu",
              dof->getDofCoeffIdx(),
              static_cast<unsigned long long>(dof->getEnt()));
    index = dof->getPetscGlobalDofIdx();
  }
 
  for (const auto &[entity, indices] : *control_indices_ptr)
    for (int coordinate = 0;
         coordinate != plasticLogarithmicStretchCoordinateSize;
         ++coordinate)
      if (indices[coordinate] < 0)
        SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Plastic control entity %llu has no global index for "
                "coefficient %d",
                static_cast<unsigned long long>(entity), coordinate);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode assembleHelmholtzStateGradient(EshelbianCore &ep,
                                              Vec state_gradient) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(state_gradient);
  auto fe = boost::make_shared<MechanicalValueEle>(ep.mField);
  fe->ksp_f = state_gradient;
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, fe);
  CHKERR ep.physicalEquations->pushHelmholtzStateGradient(
      ep, fe->getOpPtrVector(), ep.dataAtPts);
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({fe.get()});
  const MoFEMErrorCode loop_error =
      DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({fe.get()});
  CHKERR loop_error;
  CHKERR VecAssemblyBegin(state_gradient);
  CHKERR VecAssemblyEnd(state_gradient);
  CHKERR VecGhostUpdateBegin(state_gradient, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(state_gradient, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateBegin(state_gradient, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(state_gradient, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode solveStaticHelmholtzAdjoint(EshelbianCore &ep, TS ts,
                                          Vec state_gradient, Vec adjoint) {
  MoFEMFunctionBegin;
  Vec state = nullptr;
  CHKERR TSGetSolution(ts, &state);
  if (!state)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "State TS has no solution for the Helmholtz adjoint");
  auto zero_rate = vectorDuplicate(state);
  CHKERR VecZeroEntries(zero_rate);
  PetscReal pseudo_time = 0;
  CHKERR TSGetTime(ts, &pseudo_time);
 
  Mat jacobian = nullptr;
  Mat preconditioner = nullptr;
  CHKERR TSGetIJacobian(ts, &jacobian, &preconditioner, PETSC_NULLPTR,
                       PETSC_NULLPTR);
  if (!jacobian || !preconditioner)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "State TS has no IJacobian matrices for the static adjoint");
  // The outer optimization differentiates the rate-free equilibrium map.
  // shift=0 removes all pseudo-time damping from the sensitivity operator.
  CHKERR TSComputeIJacobian(ts, pseudo_time, state, zero_rate, 0., jacobian,
                           preconditioner, PETSC_FALSE);
 
  SNES snes = nullptr;
  KSP ksp = nullptr;
  CHKERR TSGetSNES(ts, &snes);
  CHKERR SNESGetKSP(snes, &ksp);
  CHKERR KSPSetOperators(ksp, jacobian, preconditioner);
  CHKERR KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
  CHKERR VecZeroEntries(adjoint);
  CHKERR KSPSolveTranspose(ksp, state_gradient, adjoint);
  KSPConvergedReason reason;
  CHKERR KSPGetConvergedReason(ksp, &reason);
  if (reason <= 0)
    SETERRQ(ep.mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
            "Static Helmholtz adjoint diverged with KSP reason %d",
            static_cast<int>(reason));
  CHKERR VecGhostUpdateBegin(adjoint, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(adjoint, INSERT_VALUES, SCATTER_FORWARD);
 
#ifndef NDEBUG
  auto adjoint_residual = vectorDuplicate(state_gradient);
  CHKERR MatMultTranspose(jacobian, adjoint, adjoint_residual);
  CHKERR VecAXPY(adjoint_residual, -1., state_gradient);
  PetscReal state_gradient_norm = 0;
  PetscReal residual_norm = 0;
  CHKERR VecNorm(state_gradient, NORM_2, &state_gradient_norm);
  CHKERR VecNorm(adjoint_residual, NORM_2, &residual_norm);
  MOFEM_LOG("EP", Sev::verbose)
      << "Static Helmholtz adjoint residual " << residual_norm
      << ", material state-gradient norm " << state_gradient_norm;
#endif
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode assemblePlasticAdjointContribution(EshelbianCore &ep,
                                                  Vec adjoint,
                                                  Vec gradient) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(gradient);
  auto control_indices = boost::shared_ptr<PlasticControlDofIndices>();
  CHKERR createPlasticControlDofIndices(ep, control_indices);
  auto adjoint_piola = boost::make_shared<MatrixDouble>();
  auto adjoint_stretch = ep.dataAtPts->getVarLogStreachPts();
  auto adjoint_rotation = boost::make_shared<MatrixDouble>();
  auto adjoint_rotation_gradient = boost::make_shared<MatrixDouble>();
  auto adjoint_vector = SmartPetscObj<Vec>(adjoint, true);
 
  auto fe = boost::make_shared<MechanicalValueEle>(ep.mField);
  // OpCalculateRotationAndSpatialGradient fills the P derivatives consumed
  // below only in Jacobian context.  Set it explicitly instead of relying on
  // values retained from the preceding static-Jacobian assembly.
  fe->ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, fe);
  fe->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
          ep.piolaStress, adjoint_piola, nullptr, adjoint_vector));
  fe->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
          ep.bubbleField, adjoint_piola, nullptr, adjoint_vector, MBMAXTYPE));
  CHKERR ep.physicalEquations->pushMaterialVariation(
      ep, fe->getOpPtrVector(), ep.dataAtPts, adjoint_vector);
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
      ep.rotAxis, adjoint_rotation, adjoint_vector, MBTET));
  fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
          ep.rotAxis, adjoint_rotation_gradient, adjoint_vector, MBTET));
  fe->getOpPtrVector().push_back(new OpAssembleReducedHelmholtzGradient(
      ep.dataAtPts, adjoint_piola, adjoint_stretch, adjoint_rotation,
      adjoint_rotation_gradient, control_indices, gradient, ep.alphaOmega));
  fe->exeTestHook = [&ep](FEMethod *method) {
    return ep.plasticVolumes->find(method->getFEEntityHandle()) !=
           ep.plasticVolumes->end();
  };
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({fe.get()});
  const MoFEMErrorCode loop_error =
      DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({fe.get()});
  CHKERR loop_error;
  CHKERR VecAssemblyBegin(gradient);
  CHKERR VecAssemblyEnd(gradient);
  CHKERR VecGhostUpdateBegin(gradient, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(gradient, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
MoFEMErrorCode assembleReducedHelmholtzGradient(EshelbianCore &ep, TS ts,
                                               Vec gradient) {
  MoFEMFunctionBegin;
  if (!ts || !gradient)
    SETERRQ(ep.mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Reduced Helmholtz gradient requires a TS and output vector");
  auto state_gradient = createDMVector(ep.dmElastic, RowColData::ROW);
  auto adjoint = vectorDuplicate(state_gradient);
  CHKERR assembleHelmholtzStateGradient(ep, state_gradient);
  CHKERR solveStaticHelmholtzAdjoint(ep, ts, state_gradient, adjoint);
  CHKERR assemblePlasticAdjointContribution(ep, adjoint, gradient);
  MoFEMFunctionReturn(0);
}
 
} // namespace PlasticIncrementalOptimizationInternal
 
} // namespace EshelbianPlasticity