MultiParameterPlasticityStressUpdate.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "StressUpdateBase.h"
 
#include <array>
 
class MultiParameterPlasticityStressUpdate : public StressUpdateBase
{
public:
  static InputParameters validParams();
 
  MultiParameterPlasticityStressUpdate(const InputParameters & parameters,
                                       unsigned num_sp,
                                       unsigned num_yf,
                                       unsigned num_intnl);
 
protected:
  virtual void initQpStatefulProperties() override;
  using StressUpdateBase::updateState;
  virtual void updateState(RankTwoTensor & strain_increment,
                           RankTwoTensor & inelastic_strain_increment,
                           const RankTwoTensor & rotation_increment,
                           RankTwoTensor & stress_new,
                           const RankTwoTensor & stress_old,
                           const RankFourTensor & elasticity_tensor,
                           const RankTwoTensor & elastic_strain_old,
                           bool compute_full_tangent_operator,
                           RankFourTensor & tangent_operator) override;
 
  virtual void propagateQpStatefulProperties() override;
 
  virtual TangentCalculationMethod getTangentCalculationMethod() override
  {
    return TangentCalculationMethod::FULL;
  }
 
  constexpr static unsigned _tensor_dimensionality = 3;
 
  const unsigned _num_sp;
 
  const std::vector<Real> _definitely_ok_sp;
 
  std::vector<std::vector<Real>> _Eij;
 
  Real _En;
 
  std::vector<std::vector<Real>> _Cij;
 
  const unsigned _num_yf;
 
  const unsigned _num_intnl;
 
  const unsigned _max_nr_its;
 
  const bool _perform_finite_strain_rotations;
 
  const Real _smoothing_tol;
 
  const Real _smoothing_tol2;
 
  const Real _f_tol;
 
  const Real _f_tol2;
 
  const Real _min_step_size;
 
  bool _step_one;
 
  const bool _warn_about_precision_loss;
 
  MaterialProperty<RankTwoTensor> & _plastic_strain;
 
  const MaterialProperty<RankTwoTensor> & _plastic_strain_old;
 
  MaterialProperty<std::vector<Real>> & _intnl;
 
  const MaterialProperty<std::vector<Real>> & _intnl_old;
 
  MaterialProperty<std::vector<Real>> & _yf;
 
  MaterialProperty<Real> & _iter;
 
  MaterialProperty<Real> & _max_iter_used;
 
  const MaterialProperty<Real> & _max_iter_used_old;
 
  MaterialProperty<Real> & _linesearch_needed;
 
  struct yieldAndFlow
  {
    Real f;                                // yield function value
    std::vector<Real> df;                  // df/d(stress_param[i])
    std::vector<Real> df_di;               // df/d(intnl_variable[a])
    std::vector<Real> dg;                  // d(flow)/d(stress_param[i])
    std::vector<std::vector<Real>> d2g;    // d^2(flow)/d(sp[i])/d(sp[j])
    std::vector<std::vector<Real>> d2g_di; // d^2(flow)/d(sp[i])/d(intnl[a])
 
    yieldAndFlow() : yieldAndFlow(0, 0) {}
 
    yieldAndFlow(unsigned num_var, unsigned num_intnl)
      : f(0.0),
        df(num_var),
        df_di(num_intnl),
        dg(num_var),
        d2g(num_var, std::vector<Real>(num_var, 0.0)),
        d2g_di(num_var, std::vector<Real>(num_intnl, 0.0))
    {
    }
 
    // Zero everything without resizing/reallocating
    void reset()
    {
      f = 0.0;
      std::fill(df.begin(), df.end(), 0.0);
      std::fill(df_di.begin(), df_di.end(), 0.0);
      std::fill(dg.begin(), dg.end(), 0.0);
      for (auto & row : d2g)
        std::fill(row.begin(), row.end(), 0.0);
      for (auto & row : d2g_di)
        std::fill(row.begin(), row.end(), 0.0);
    }
 
    // this may be involved in the smoothing of a group of yield functions
    bool operator<(const yieldAndFlow & fd) const { return f < fd.f; }
  };
 
  Real yieldF(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;
 
  Real yieldF(const std::vector<Real> & yfs) const;
 
  Real ismoother(Real f_diff) const;
 
  Real smoother(Real f_diff) const;
 
  Real dsmoother(Real f_diff) const;
 
  yieldAndFlow smoothAllQuantities(const std::vector<Real> & stress_params,
                                   const std::vector<Real> & intnl) const;
 
  int lineSearch(Real & res2,
                 std::vector<Real> & stress_params,
                 Real & gaE,
                 const std::vector<Real> & trial_stress_params,
                 yieldAndFlow & smoothed_q,
                 const std::vector<Real> & intnl_ok,
                 std::vector<Real> & intnl,
                 std::vector<Real> & rhs,
                 Real & linesearch_needed) const;
 
  int nrStep(const yieldAndFlow & smoothed_q,
             const std::vector<Real> & trial_stress_params,
             const std::vector<Real> & stress_params,
             const std::vector<Real> & intnl,
             Real gaE,
             std::vector<Real> & rhs) const;
 
  Real calculateRes2(const std::vector<Real> & rhs) const;
 
  void calculateRHS(const std::vector<Real> & trial_stress_params,
                    const std::vector<Real> & stress_params,
                    Real gaE,
                    const yieldAndFlow & smoothed_q,
                    std::vector<Real> & rhs) const;
 
  void dnRHSdVar(const yieldAndFlow & smoothed_q,
                 const std::vector<std::vector<Real>> & dintnl,
                 const std::vector<Real> & stress_params,
                 Real gaE,
                 std::vector<double> & jac) const;
 
  virtual void errorHandler(const std::string & message) const;
 
  virtual void yieldFunctionValuesV(const std::vector<Real> & stress_params,
                                    const std::vector<Real> & intnl,
                                    std::vector<Real> & yf) const = 0;
 
  virtual void computeAllQV(const std::vector<Real> & stress_params,
                            const std::vector<Real> & intnl,
                            std::vector<yieldAndFlow> & all_q) const = 0;
 
  virtual void preReturnMapV(const std::vector<Real> & trial_stress_params,
                             const RankTwoTensor & stress_trial,
                             const std::vector<Real> & intnl_old,
                             const std::vector<Real> & yf,
                             const RankFourTensor & Eijkl);
 
  virtual void initializeVarsV(const std::vector<Real> & trial_stress_params,
                               const std::vector<Real> & intnl_old,
                               std::vector<Real> & stress_params,
                               Real & gaE,
                               std::vector<Real> & intnl) const;
 
  virtual void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
                               const std::vector<Real> & current_stress_params,
                               const std::vector<Real> & intnl_old,
                               std::vector<Real> & intnl) const = 0;
 
  virtual void setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
                                    const std::vector<Real> & current_stress_params,
                                    const std::vector<Real> & intnl,
                                    std::vector<std::vector<Real>> & dintnl) const = 0;
 
  virtual void computeStressParams(const RankTwoTensor & stress,
                                   std::vector<Real> & stress_params) const = 0;
 
  virtual void initializeReturnProcess();
 
  virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment);
 
  virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial,
                                     const std::vector<Real> & stress_params,
                                     Real gaE,
                                     const std::vector<Real> & intnl,
                                     const yieldAndFlow & smoothed_q,
                                     const RankFourTensor & Eijkl,
                                     RankTwoTensor & stress) const = 0;
 
  virtual void setInelasticStrainIncrementAfterReturn(const RankTwoTensor & stress_trial,
                                                      Real gaE,
                                                      const yieldAndFlow & smoothed_q,
                                                      const RankFourTensor & elasticity_tensor,
                                                      const RankTwoTensor & returned_stress,
                                                      RankTwoTensor & inelastic_strain_increment);
 
  virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial,
                                          const std::vector<Real> & trial_stress_params,
                                          const RankTwoTensor & stress,
                                          const std::vector<Real> & stress_params,
                                          Real gaE,
                                          const yieldAndFlow & smoothed_q,
                                          const RankFourTensor & Eijkl,
                                          bool compute_full_tangent_operator,
                                          const std::vector<std::vector<Real>> & dvar_dtrial,
                                          RankFourTensor & cto);
 
  virtual void dstressparam_dstress(const RankTwoTensor & stress,
                                    std::vector<RankTwoTensor> & dsp) const = 0;
 
  virtual void d2stressparam_dstress(const RankTwoTensor & stress,
                                     std::vector<RankFourTensor> & d2sp) const = 0;
 
  virtual void setEffectiveElasticity(const RankFourTensor & Eijkl) = 0;
 
  void dVardTrial(bool elastic_only,
                  const std::vector<Real> & trial_stress_params,
                  const std::vector<Real> & stress_params,
                  Real gaE,
                  const std::vector<Real> & intnl,
                  const yieldAndFlow & smoothed_q,
                  Real step_size,
                  bool compute_full_tangent_operator,
                  std::vector<std::vector<Real>> & dvar_dtrial) const;
 
  bool precisionLoss(const std::vector<Real> & solution,
                     const std::vector<Real> & stress_params,
                     Real gaE) const;
 
private:
  std::vector<Real> _trial_sp;
 
  RankTwoTensor _stress_trial;
 
  std::vector<Real> _rhs;
 
  std::vector<std::vector<Real>> _dvar_dtrial;
 
  std::vector<Real> _ok_sp;
 
  std::vector<Real> _ok_intnl;
 
  std::vector<Real> _del_stress_params;
 
  std::vector<Real> _current_sp;
 
  std::vector<Real> _current_intnl;
 
  yieldAndFlow _smoothed_q;
 
  mutable std::vector<RankTwoTensor> _dsp_scratch;
 
  mutable std::vector<RankTwoTensor> _dsp_trial_scratch;
 
  mutable std::vector<RankFourTensor> _d2sp_scratch;
 
  mutable std::vector<Real> _yfs_scratch;
 
  mutable std::vector<Real> _sp_params_old_scratch;
 
  mutable std::vector<Real> _delta_nr_params_scratch;
 
  mutable std::vector<std::vector<Real>> _dintnl_scratch;
 
  mutable std::vector<double> _jac_scratch;
 
  mutable std::vector<PetscBLASInt> _ipiv_scratch;
 
  mutable std::vector<yieldAndFlow> _all_q_scratch;
 
private:
  enum class SmootherFunctionType
  {
    cos,
    poly1,
    poly2,
    poly3
  } _smoother_function_type;
};