PlasticIncrementalOptimizationInternal.hpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalOptimizationInternal.hpp
 * @brief Shared implementation details for plastic incremental optimization
 */
 
#ifndef MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_INTERNAL_HPP
#define MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_INTERNAL_HPP
 
 
namespace EshelbianPlasticity::PlasticIncrementalOptimizationInternal {
 
/** Restore a borrowed PETSc view while unwinding from a failed callback. */
template <typename Restore> class PetscRestoreGuard {
public:
  explicit PetscRestoreGuard(Restore restore)
      : restoreFunction(std::move(restore)) {}
  PetscRestoreGuard(const PetscRestoreGuard &) = delete;
  PetscRestoreGuard &operator=(const PetscRestoreGuard &) = delete;
 
  ~PetscRestoreGuard() {
    if (active)
      (void)restoreFunction();
  }
 
  MoFEMErrorCode restore() {
    active = false;
    return restoreFunction();
  }
 
private:
  Restore restoreFunction;
  bool active = true;
};
 
/** Reconstruct a trace-free tensor in the orthonormal basis of Eq. (1.52),
 * label eq:sym-basis. */
template <typename T>
FTensor::Tensor2_symmetric<double, SPACE_DIM>
plasticLogarithmicStretchTensorFromCoordinates(
    const FTensor::Tensor1<T, plasticLogarithmicStretchCoordinateSize>
        &t_coordinates) {
  static_assert(SPACE_DIM == 3,
                "The five-coordinate plastic basis is three-dimensional");
  return Tensor2SymmetricDeviatorBasis::getTensor(t_coordinates);
}
 
/** Contract a symmetric tensor with the orthonormal basis of Eq. (1.52),
 * label eq:sym-basis. */
template <typename T>
FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
plasticLogarithmicStretchCoordinatesFromTensor(
    const FTensor::Tensor2_symmetric<T, SPACE_DIM> &t_values) {
  static_assert(SPACE_DIM == 3,
                "The five-coordinate plastic basis is three-dimensional");
  return Tensor2SymmetricDeviatorBasis::getCoordinates(t_values);
}
 
/** Contract a full tensor with the symmetric trace-free basis of Eq. (1.52),
 * label eq:sym-basis. Only its symmetric part contributes. */
template <typename T>
FTensor::Tensor1<double, plasticLogarithmicStretchCoordinateSize>
plasticLogarithmicStretchCoordinatesFromTensor(
    const FTensor::Tensor2<T, SPACE_DIM, SPACE_DIM> &t_values) {
  static_assert(SPACE_DIM == 3,
                "The five-coordinate plastic basis is three-dimensional");
  return Tensor2SymmetricDeviatorBasis::getCoordinates(t_values);
}
 
/** Euclidean coordinate norm, equal to the Frobenius norm by Eq. (1.53),
 * label eq:symmetric-basis-gram. */
template <typename T>
double plasticCoordinateNorm(
    const FTensor::Tensor1<T, plasticLogarithmicStretchCoordinateSize>
        &t_values) {
  FTENSOR_INDEX(plasticLogarithmicStretchCoordinateSize, L);
  return std::sqrt(t_values(L) * t_values(L));
}
 
/** Frobenius norm of a symmetric tensor. */
template <typename T>
double
plasticFrobeniusNorm(const FTensor::Tensor2_symmetric<T, SPACE_DIM> &t_values) {
  FTENSOR_INDEXES(SPACE_DIM, i, j);
  return std::sqrt(t_values(i, j) * t_values(i, j));
}
 
constexpr double plasticEquivalentIncrementScale = 0.81649658092772603273;
constexpr double plasticEquivalentIncrementScaleSquared = 2. / 3.;
constexpr double plasticYieldFunctionScale = 1.22474487139158904910;
 
/** Evaluate alpha_epsilon from Eq. (1.81), including the exact epsilon=0
 * limit.  The quotient avoids cancellation close to the rounded apex. */
inline double plasticEquivalentIncrement(const double coordinate_norm,
                                         const double epsilon) {
  const double scaled_norm = plasticEquivalentIncrementScale * coordinate_norm;
  if (epsilon == 0)
    return scaled_norm;
  return scaled_norm * scaled_norm /
         (std::hypot(scaled_norm, epsilon) + epsilon);
}
 
using PackedPlasticCoordinates =
    std::array<double, plasticLogarithmicStretchCoordinateSize>;
using PlasticHistory = std::map<EntityHandle, PackedPlasticCoordinates>;
using PlasticKappaHistory = std::map<EntityHandle, double>;
 
MoFEMErrorCode savePlasticHistory(EshelbianCore &ep, PlasticHistory &history);
 
MoFEMErrorCode setPlasticHistory(
    EshelbianCore &ep, const PlasticHistory &history,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_direction,
    double scale);
 
MoFEMErrorCode savePlasticKappaHistory(EshelbianCore &ep,
                                       PlasticKappaHistory &history);
 
MoFEMErrorCode setPlasticKappaHistory(EshelbianCore &ep,
                                      const PlasticKappaHistory &history);
 
/** Clear only increment-valued mesh fields, preserving committed kappa. */
MoFEMErrorCode clearPlasticIncrementFields(EshelbianCore &ep);
 
MoFEMErrorCode setStateOnMesh(EshelbianCore &ep, Vec state);
 
MoFEMErrorCode evaluatePlasticIncrementalResistanceValue(
    const PlasticIncrementalOptimizationProblem &problem, Vec control,
    double &value);
 
MoFEMErrorCode assembleIncrementalObjectiveGradient(
    PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec smooth_gradient, Vec objective_gradient);
 
MoFEMErrorCode initialisePlasticInequalityMultipliers(
    const PlasticIncrementalOptimizationProblem &problem, DM constraint_dm,
    Vec multipliers);
 
MoFEMErrorCode evaluatePlasticInequalityConstraints(
    const PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec constraints);
 
MoFEMErrorCode assemblePlasticInequalityJacobian(
    PlasticIncrementalOptimizationProblem &problem, Vec control,
    Vec smooth_gradient, Mat jacobian);
 
MoFEMErrorCode validatePlasticIncrementalSolution(
    PlasticIncrementalOptimizationProblem &problem, DM constraint_dm,
    Vec control, Vec smooth_gradient, Vec inequality_multipliers,
    PetscReal gradient_tolerance, PetscReal constraint_tolerance,
    std::string &diagnostics);
 
/** Assemble the fixed-state mixed-kinematic plastic conjugate diagnostic. */
MoFEMErrorCode assemblePlasticForce(EshelbianCore &ep, Vec force);
 
/** Assemble the reduced bulk Helmholtz-energy plastic gradient by adjoint. */
MoFEMErrorCode assembleReducedHelmholtzGradient(EshelbianCore &ep, TS ts,
                                               Vec gradient);
 
} // namespace EshelbianPlasticity::PlasticIncrementalOptimizationInternal
 
#endif // MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_INTERNAL_HPP