PlasticIncrementalOptimization.hpp

Go to the documentation of this file.

/**
 * @file PlasticIncrementalOptimization.hpp
 * @brief Plasticity implementation of incremental optimization
 */
 
#ifndef MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_HPP
#define MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_HPP
 
#include <IncrementalOptimization.hpp>
 
namespace EshelbianPlasticity {
 
struct EshelbianCore;
 
/** Number of coefficients in the trace-free approximation of H^p. */
constexpr int plasticLogarithmicStretchCoordinateSize =
    SPACE_DIM * (SPACE_DIM + 1) / 2 - 1;
 
/** Plasticity adapter; no plastic data is present in the generic context. */
struct PlasticIncrementalOptimizationProblem
    : public IncrementalOptimizationProblem {
  PlasticIncrementalOptimizationProblem(EshelbianCore &ep,
                                        SmartPetscObj<TS> ts);
 
  std::string getName() const override;
  MoFEMErrorCode validate() override;
  DM getControlDM() const override;
  MoFEMErrorCode
  createControlData(SmartPetscObj<Vec> &committed_control,
                    SmartPetscObj<Vec> &candidate_control) override;
  MoFEMErrorCode configureConstraintDM(DM constraint_dm) override;
  /** Initialise multipliers for sqrt(cell volume)-scaled inequality rows. */
  MoFEMErrorCode initialiseInequalityMultipliers(DM constraint_dm,
                                                 Vec multipliers) override;
  MoFEMErrorCode initialiseReferenceControl(Vec committed_control,
                                            Vec reference_control) override;
  MoFEMErrorCode solveState(Vec state, EquilibriumSolveReport &report) override;
  MoFEMErrorCode restoreState(Vec state) override;
  MoFEMErrorCode evaluateConservativeValue(double &value) override;
  MoFEMErrorCode assembleConservativeGradient(Vec smooth_gradient) override;
  MoFEMErrorCode evaluateDissipation(Vec control, double &value) override;
  MoFEMErrorCode assembleObjectiveGradient(Vec control, Vec smooth_gradient,
                                           Vec objective_gradient) override;
  /** Return sqrt(cell volume) times the exact or rounded epigraph gap. */
  MoFEMErrorCode evaluateInequalityConstraints(Vec control,
                                               Vec constraints) override;
  MoFEMErrorCode assembleInequalityJacobian(Vec control, Vec smooth_gradient,
                                            Mat jacobian) override;
  MoFEMErrorCode assembleInequalityJacobianWithoutSmoothGradient(
      Vec control, Mat jacobian) override;
  /** Recover physical multipliers and check unscaled plastic feasibility. */
  MoFEMErrorCode validateSolution(
      DM constraint_dm, Vec control, Vec smooth_gradient,
      Vec inequality_multipliers, PetscReal gradient_tolerance,
      PetscReal constraint_tolerance, std::string &diagnostics) override;
  MoFEMErrorCode commit(Vec solution, Vec equilibrated_state,
                        Vec committed_control, Vec candidate_control,
                        Vec reference_control) override;
 
  /** Plastic-flow control field; its vector values are Delta H^p. */
  const std::string &getPlasticFlowFieldName() const;
  /** Scalar control/constraint field; its vector values are Delta kappa. */
  const std::string &getPlasticKappaFieldName() const;
  double getInitialYieldStress() const { return initialYieldStress; }
  double getIsotropicHardeningModulus() const {
    return isotropicHardeningModulus;
  }
  double getDissipationRegularizationEpsilon() const {
    return dissipationRegularizationEpsilon;
  }
  const std::string &getVolumeElementName() const;
  /** Add geometry operators for local constitutive integration. */
  MoFEMErrorCode addConstitutiveGeometryOperators(
      boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline) const;
 
private:
  EshelbianCore &ep;
  SmartPetscObj<TS> timeSolver;
  /** Zero force selects a finite exact-cone row when no smooth force exists. */
  SmartPetscObj<Vec> fallbackJacobianSmoothGradient;
  double initialYieldStress = 1;
  double isotropicHardeningModulus = 0;
  double dissipationRegularizationEpsilon = 0;
};
 
boost::shared_ptr<IncrementalOptimizationContext>
createPlasticIncrementalOptimizationContext(EshelbianCore &ep,
                                            SmartPetscObj<TS> ts,
                                            SmartPetscObj<Vec> state);
 
PlasticIncrementalOptimizationProblem &getPlasticIncrementalOptimizationProblem(
    const boost::shared_ptr<IncrementalOptimizationContext> &context);
 
/** Reject combinations for which the material-energy adjoint is incomplete. */
MoFEMErrorCode validateEquilibratedMechanicalValueScope(EshelbianCore &ep);
 
/**
 * Evaluate only the bulk material Helmholtz free energy at equilibrium.
 * Boundary-device and stabilization energies are excluded; their influence
 * enters the reduced derivative through the equilibrium adjoint.
 */
MoFEMErrorCode evaluateEquilibratedMechanicalValue(EshelbianCore &ep, TS ts,
                                                   double &value);
 
/**
 * Fixed-control state solve corresponding to Eq. (3.5), label
 * gen:eq:state-map, which reports failure without aborting a trial.
 *
 * The outer physical/load time is EshelbianCore::currentPhysicalTime.  The TS
 * time used here is reset for every control trial and is only pseudo-time for
 * first-order dynamic relaxation.  Its rate terms are excluded from both the
 * Helmholtz value and the static sensitivity operator, and the returned
 * equilibrium norm is recomputed with all pseudo-time rates set to zero.
 * Reaching the TS step limit before the final pseudo-time emits a warning;
 * the caller still checks the final state's rate-free equilibrium norm.
 */
MoFEMErrorCode solveEquilibriumStateTS(EshelbianCore &ep, TS ts, Vec state,
                                       EquilibriumSolveReport &report);
 
/**
 * Evaluate five trace-free logarithmic-plastic-stretch coordinates as a full
 * symmetric tensor at volume quadrature points. Implements Eqs.
 * (1.50)--(1.52), labels eq:increment-expansion through
 * eq:symmetric-basis-gram.
 */
MoFEMErrorCode addCalculatePlasticLogarithmicStretchFieldValues(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
    const std::string &field_name,
    boost::shared_ptr<MatrixDouble> tensor_values,
    const EntityType zero_type = MBTET, SmartPetscObj<DM> data_dm = {},
    SmartPetscObj<Vec> data_vector = {});
 
/**
 * Evaluate the exact fixed-state discrete plastic variation
 * -integral P : Fe Dexp(Hp)[W] for a constant symmetric direction. See Eqs.
 * (1.29)--(1.30), labels eq:fixed-state-plastic-variation and
 * eq:plastic-variation-conjugacy.
 */
MoFEMErrorCode evaluatePlasticConjugateWork(
    EshelbianCore &ep,
    const FTensor::Tensor2_symmetric<double, SPACE_DIM> &t_direction,
    double &conjugate_work);
 
/**
 * Centred equilibrated-value/conjugate-force test for Eq. (1.87), label
 * eq:energy-force-check.
 */
MoFEMErrorCode testEquilibratedMechanicalValue(EshelbianCore &ep, TS ts,
                                               Vec state);
 
/**
 * Control layout, reconstruction, commit, and rollback test for Eqs.
 * (1.50)--(1.55), labels eq:history-space through eq:quadrature-weight.
 */
MoFEMErrorCode testIncrementalOptimizationLayout(EshelbianCore &ep, TS ts,
                                                 Vec state);
 
/**
 * State/value/force/cache/rollback transaction test for the invariants in
 * Eq. (3.25), label gen:eq:common-tests.
 */
MoFEMErrorCode testIncrementalOptimizationTransaction(EshelbianCore &ep, TS ts,
                                                      Vec state);
 
/**
 * Branch-controlled directional finite differences of the complete
 * incremental objective and the exact TAO-facing gradient.
 */
MoFEMErrorCode testIncrementalOptimizationObjectiveDerivative(EshelbianCore &ep,
                                                              TS ts, Vec state);
 
MoFEMErrorCode
testIncrementalOptimizationConstraintDerivative(EshelbianCore &ep, TS ts,
                                                Vec state);
 
} // namespace EshelbianPlasticity
 
#endif // MOFEM_PLASTIC_INCREMENTAL_OPTIMIZATION_HPP