www.mooseframework.org
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
ComputeMultiPlasticityStress Class Reference

ComputeMultiPlasticityStress performs the return-map algorithm and associated stress updates for plastic models defined by a General User Objects. More...

#include <ComputeMultiPlasticityStress.h>

Inheritance diagram for ComputeMultiPlasticityStress:
[legend]

Public Member Functions

 ComputeMultiPlasticityStress (const InputParameters &parameters)
 
void outputAndCheckDebugParameters ()
 Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly. More...
 
void checkDerivatives ()
 Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations. More...
 
void checkJacobian (const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
 Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations. More...
 
void checkSolution (const RankFourTensor &E_inv)
 Checks that Ax does equal b in the NR procedure. More...
 

Protected Types

enum  TangentOperatorEnum { elastic, linear, nonlinear }
 The type of tangent operator to return. tangent operator = d(stress_rate)/d(strain_rate). More...
 
enum  DeactivationSchemeEnum {
  optimized, safe, dumb, optimized_to_safe,
  safe_to_dumb, optimized_to_safe_to_dumb, optimized_to_dumb
}
 
enum  quickStep_called_from_t { computeQpStress_function, returnMap_function }
 The functions from which quickStep can be called. More...
 

Protected Member Functions

virtual void computeQpStress ()
 
virtual void initQpStatefulProperties ()
 
virtual bool reinstateLinearDependentConstraints (std::vector< bool > &deactivated_due_to_ld)
 makes all deactivated_due_to_ld false, and if >0 of them were initially true, returns true More...
 
virtual unsigned int numberActive (const std::vector< bool > &active)
 counts the number of active constraints More...
 
virtual Real residual2 (const std::vector< Real > &pm, const std::vector< Real > &f, const RankTwoTensor &epp, const std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld)
 The residual-squared. More...
 
virtual bool returnMap (const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &f, unsigned int &iter, bool can_revert_to_dumb, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, bool final_step, RankFourTensor &consistent_tangent_operator, std::vector< Real > &cumulative_pm)
 Implements the return map. More...
 
virtual bool lineSearch (Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, const RankFourTensor &E_inv, RankTwoTensor &delta_dp, const RankTwoTensor &dstress, const std::vector< Real > &dpm, const std::vector< Real > &dintnl, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, bool &linesearch_needed)
 Performs a line search. More...
 
virtual bool singleStep (Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, RankTwoTensor &delta_dp, const RankFourTensor &E_inv, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, std::vector< bool > &active, DeactivationSchemeEnum deactivation_scheme, bool &linesearch_needed, bool &ld_encountered)
 Performs a single Newton-Raphson + linesearch step Constraints are deactivated and the step is re-done if deactivation_scheme is set appropriately. More...
 
virtual bool checkAdmissible (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &all_f)
 Checks whether the yield functions are in the admissible region. More...
 
void buildDumbOrder (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
 Builds the order which "dumb" activation will take. More...
 
virtual void incrementDumb (int &dumb_iteration, const std::vector< unsigned int > &dumb_order, std::vector< bool > &act)
 Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of dumb_iteration == 1) More...
 
virtual bool checkKuhnTucker (const std::vector< Real > &f, const std::vector< Real > &pm, const std::vector< bool > &active)
 Checks Kuhn-Tucker conditions, and alters "active" if appropriate. More...
 
virtual void applyKuhnTucker (const std::vector< Real > &f, const std::vector< Real > &pm, std::vector< bool > &active)
 Checks Kuhn-Tucker conditions, and alters "active" if appropriate. More...
 
virtual void preReturnMap ()
 
virtual void postReturnMap ()
 
virtual bool quickStep (const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, RankFourTensor &consistent_tangent_operator, const quickStep_called_from_t called_from, bool final_step)
 Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined through the TensorMechanicsPlasticXXXX.returnMap functions. More...
 
virtual bool plasticStep (const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, RankFourTensor &consistent_tangent_operator)
 performs a plastic step More...
 
bool canChangeScheme (DeactivationSchemeEnum current_deactivation_scheme, bool can_revert_to_dumb)
 
bool canIncrementDumb (int dumb_iteration)
 
void changeScheme (const std::vector< bool > &initial_act, bool can_revert_to_dumb, const RankTwoTensor &initial_stress, const std::vector< Real > &intnl_old, DeactivationSchemeEnum &current_deactivation_scheme, std::vector< bool > &act, int &dumb_iteration, std::vector< unsigned int > &dumb_order)
 
bool canAddConstraints (const std::vector< bool > &act, const std::vector< Real > &all_f)
 
unsigned int activeCombinationNumber (const std::vector< bool > &act)
 
RankFourTensor consistentTangentOperator (const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &pm_this_step, const std::vector< Real > &cumulative_pm)
 Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rate) More...
 
virtual void computeQpProperties () override
 
virtual void calculateConstraints (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
 The constraints. More...
 
virtual void calculateRHS (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
 Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic]) More...
 
virtual void calculateJacobian (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
 d(rhs)/d(dof) More...
 
virtual void nrStep (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
 Performs one Newton-Raphson step. More...
 
virtual void yieldFunction (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
 The active yield function(s) More...
 
virtual void dyieldFunction_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
 The derivative of the active yield function(s) with respect to stress. More...
 
virtual void dyieldFunction_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
 The derivative of active yield function(s) with respect to their internal parameters (the user objects assume there is exactly one internal param per yield function) More...
 
virtual void flowPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
 The active flow potential(s) - one for each yield function. More...
 
virtual void dflowPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
 The derivative of the active flow potential(s) with respect to stress. More...
 
virtual void dflowPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
 The derivative of the active flow potentials with respect to the active internal parameters The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha]. More...
 
virtual void hardPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
 The active hardening potentials (one for each internal parameter and for each yield function) by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part of model a, so we only calculate those here. More...
 
virtual void dhardPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
 The derivative of the active hardening potentials with respect to stress By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only calculate those here. More...
 
virtual void dhardPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
 The derivative of the active hardening potentials with respect to the active internal parameters. More...
 
virtual void buildActiveConstraints (const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
 Constructs a set of active constraints, given the yield functions, f. More...
 
unsigned int modelNumber (unsigned int surface)
 returns the model number, given the surface number More...
 
bool anyActiveSurfaces (int model, const std::vector< bool > &active)
 returns true if any internal surfaces of the given model are active according to 'active' More...
 
void activeModelSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
 Returns the internal surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model. More...
 
void activeSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
 Returns the external surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model. More...
 
bool returnMapAll (const RankTwoTensor &trial_stress, const std::vector< Real > &intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &stress, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, RankTwoTensor &delta_dp, std::vector< Real > &yf, unsigned &num_successful_plastic_returns, unsigned &custom_model)
 Performs a returnMap for each plastic model using their inbuilt returnMap functions. More...
 

Protected Attributes

unsigned int _max_iter
 Maximum number of Newton-Raphson iterations allowed. More...
 
Real _min_stepsize
 Minimum fraction of applied strain that may be applied during adaptive stepsizing. More...
 
Real _max_stepsize_for_dumb
 "dumb" deactivation will only be used if the stepsize falls below this quantity More...
 
bool _ignore_failures
 Even if the returnMap fails, return the best values found for stress and internal parameters. More...
 
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
 
Real _epp_tol
 Tolerance on the plastic strain increment ("direction") constraint. More...
 
std::vector< Real > _dummy_pm
 dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStress More...
 
std::vector< Real > _cumulative_pm
 the sum of the plastic multipliers over all the sub-steps. More...
 
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
 
bool _n_supplied
 User supplied the transverse direction vector. More...
 
RealVectorValue _n_input
 the supplied transverse direction vector More...
 
RealTensorValue _rot
 rotation matrix that takes _n to (0, 0, 1) More...
 
bool _perform_finite_strain_rotations
 whether to perform the rotations necessary in finite-strain simulations More...
 
MaterialProperty< RankTwoTensor > & _plastic_strain
 plastic strain More...
 
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
 Old value of plastic strain. More...
 
MaterialProperty< std::vector< Real > > & _intnl
 internal parameters More...
 
const MaterialProperty< std::vector< Real > > & _intnl_old
 old values of internal parameters More...
 
MaterialProperty< std::vector< Real > > & _yf
 yield functions More...
 
MaterialProperty< Real > & _iter
 Number of Newton-Raphson iterations used in the return-map. More...
 
MaterialProperty< Real > & _linesearch_needed
 Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< Real > & _ld_encountered
 Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< Real > & _constraints_added
 Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< RealVectorValue > & _n
 current value of transverse direction More...
 
const MaterialProperty< RealVectorValue > & _n_old
 old value of transverse direction More...
 
const MaterialProperty< RankTwoTensor > & _strain_increment
 strain increment (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _total_strain_old
 Old value of total strain (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _rotation_increment
 Rotation increment (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _stress_old
 Old value of stress. More...
 
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
 Old value of elastic strain. More...
 
bool _cosserat
 whether Cosserat mechanics should be used More...
 
const MaterialProperty< RankTwoTensor > * _curvature
 The Cosserat curvature strain. More...
 
const MaterialProperty< RankFourTensor > * _elastic_flexural_rigidity_tensor
 The Cosserat elastic flexural rigidity tensor. More...
 
MaterialProperty< RankTwoTensor > * _couple_stress
 the Cosserat couple-stress More...
 
const MaterialProperty< RankTwoTensor > * _couple_stress_old
 the old value of Cosserat couple-stress More...
 
MaterialProperty< RankFourTensor > * _Jacobian_mult_couple
 derivative of couple-stress w.r.t. curvature More...
 
RankFourTensor _my_elasticity_tensor
 Elasticity tensor that can be rotated by this class (ie, its not const) More...
 
RankTwoTensor _my_strain_increment
 Strain increment that can be rotated by this class, and split into multiple increments (ie, its not const) More...
 
RankFourTensor _my_flexural_rigidity_tensor
 Flexual rigidity tensor that can be rotated by this class (ie, its not const) More...
 
RankTwoTensor _my_curvature
 Curvature that can be rotated by this class, and split into multiple increments (ie, its not const) More...
 
const std::string _base_name
 
const std::string _elasticity_tensor_name
 
const MaterialProperty< RankTwoTensor > & _mechanical_strain
 
MaterialProperty< RankTwoTensor > & _stress
 
MaterialProperty< RankTwoTensor > & _elastic_strain
 
const MaterialProperty< RankFourTensor > & _elasticity_tensor
 
const MaterialProperty< RankTwoTensor > & _extra_stress
 Extra stress tensor. More...
 
std::vector< Function * > _initial_stress_fcn
 initial stress components More...
 
MaterialProperty< RankFourTensor > & _Jacobian_mult
 derivative of stress w.r.t. strain (_dstress_dstrain) More...
 
MooseEnum _fspb_debug
 none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobian_and_linear_system - check entire jacobian and check that Ax=b More...
 
RankTwoTensor _fspb_debug_stress
 Debug the Jacobian entries at this stress. More...
 
std::vector< Real > _fspb_debug_pm
 Debug the Jacobian entires at these plastic multipliers. More...
 
std::vector< Real > _fspb_debug_intnl
 Debug the Jacobian entires at these internal parameters. More...
 
Real _fspb_debug_stress_change
 Debug finite-differencing parameter for the stress. More...
 
std::vector< Real > _fspb_debug_pm_change
 Debug finite-differencing parameters for the plastic multipliers. More...
 
std::vector< Real > _fspb_debug_intnl_change
 Debug finite-differencing parameters for the internal parameters. More...
 
Real _svd_tol
 Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent. More...
 
Real _min_f_tol
 Minimum value of the _f_tol parameters for the Yield Function User Objects. More...
 
const InputParameters & _params
 
unsigned int _num_models
 Number of plastic models for this material. More...
 
unsigned int _num_surfaces
 Number of surfaces within the plastic models. More...
 
std::vector< std::vector< unsigned int > > _surfaces_given_model
 _surfaces_given_model[model_number] = vector of surface numbers for this model More...
 
MooseEnum _specialIC
 Allows initial set of active constraints to be chosen optimally. More...
 
std::vector< const TensorMechanicsPlasticModel * > _f
 User objects that define the yield functions, flow potentials, etc. More...
 

Private Member Functions

RankTwoTensor rot (const RankTwoTensor &tens)
 

Detailed Description

ComputeMultiPlasticityStress performs the return-map algorithm and associated stress updates for plastic models defined by a General User Objects.

Note that if run in debug mode you might have to use the –no-trap-fpe flag because PETSc-LAPACK-BLAS explicitly compute 0/0 and 1/0, and this causes Libmesh to trap the floating-point exceptions

Definition at line 31 of file ComputeMultiPlasticityStress.h.

Member Enumeration Documentation

◆ DeactivationSchemeEnum

◆ quickStep_called_from_t

The functions from which quickStep can be called.

Enumerator
computeQpStress_function 
returnMap_function 

Definition at line 439 of file ComputeMultiPlasticityStress.h.

◆ TangentOperatorEnum

The type of tangent operator to return. tangent operator = d(stress_rate)/d(strain_rate).

Enumerator
elastic 
linear 
nonlinear 

Definition at line 53 of file ComputeMultiPlasticityStress.h.

Constructor & Destructor Documentation

◆ ComputeMultiPlasticityStress()

ComputeMultiPlasticityStress::ComputeMultiPlasticityStress ( const InputParameters &  parameters)

Definition at line 98 of file ComputeMultiPlasticityStress.C.

99  : ComputeStressBase(parameters),
101  _max_iter(getParam<unsigned int>("max_NR_iterations")),
102  _min_stepsize(getParam<Real>("min_stepsize")),
103  _max_stepsize_for_dumb(getParam<Real>("max_stepsize_for_dumb")),
104  _ignore_failures(getParam<bool>("ignore_failures")),
105 
106  _tangent_operator_type((TangentOperatorEnum)(int)getParam<MooseEnum>("tangent_operator")),
107 
108  _epp_tol(getParam<Real>("ep_plastic_tolerance")),
109 
110  _dummy_pm(0),
111 
112  _cumulative_pm(0),
113 
114  _deactivation_scheme((DeactivationSchemeEnum)(int)getParam<MooseEnum>("deactivation_scheme")),
115 
116  _n_supplied(parameters.isParamValid("transverse_direction")),
117  _n_input(_n_supplied ? getParam<RealVectorValue>("transverse_direction") : RealVectorValue()),
118  _rot(RealTensorValue()),
119 
120  _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
121 
122  _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
123  _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
124  _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
125  _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
126  _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
127  _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
128  // for visualisation i convert it to Real
129  _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
130  // boolean, but for
131  // visualisation i
132  // convert it to Real
133  _ld_encountered(declareProperty<Real>(
134  "plastic_linear_dependence_encountered")), // this is really a boolean, but for
135  // visualisation i convert it to Real
136  _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
137  // boolean, but for
138  // visualisation i
139  // convert it to Real
140  _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
141  _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
142 
143  _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
144  _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
146  getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
147 
148  _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
149  _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
150 
151  // TODO: This design does NOT work. It makes these materials construction order dependent and it
152  // disregards block restrictions.
153  _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
154  hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
155  _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : NULL),
157  _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
158  : NULL),
159  _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : NULL),
160  _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress") : NULL),
161  _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
162  : NULL),
163 
164  _my_elasticity_tensor(RankFourTensor()),
165  _my_strain_increment(RankTwoTensor()),
166  _my_flexural_rigidity_tensor(RankFourTensor()),
167  _my_curvature(RankTwoTensor())
168 {
169  if (_epp_tol <= 0)
170  mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
171 
172  if (_n_supplied)
173  {
174  // normalise the inputted transverse_direction
175  if (_n_input.norm() == 0)
176  mooseError(
177  "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
178  else
179  _n_input /= _n_input.norm();
180  }
181 
182  if (_num_surfaces == 1)
184 }
bool _perform_finite_strain_rotations
whether to perform the rotations necessary in finite-strain simulations
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment (coming from ComputeIncrementalSmallStrain, for example)
const MaterialProperty< RealVectorValue > & _n_old
old value of transverse direction
MaterialProperty< Real > & _constraints_added
Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
const MaterialProperty< RankTwoTensor > & _strain_increment
strain increment (coming from ComputeIncrementalSmallStrain, for example)
TangentOperatorEnum
The type of tangent operator to return. tangent operator = d(stress_rate)/d(strain_rate).
const MaterialProperty< RankTwoTensor > & _stress_old
Old value of stress.
unsigned int _max_iter
Maximum number of Newton-Raphson iterations allowed.
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
bool _ignore_failures
Even if the returnMap fails, return the best values found for stress and internal parameters...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Old value of elastic strain.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool _n_supplied
User supplied the transverse direction vector.
MaterialProperty< Real > & _ld_encountered
Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true...
MaterialProperty< std::vector< Real > > & _yf
yield functions
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
std::vector< Real > _dummy_pm
dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStre...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
const MaterialProperty< RankTwoTensor > * _couple_stress_old
the old value of Cosserat couple-stress
RealVectorValue _n_input
the supplied transverse direction vector
Real _max_stepsize_for_dumb
"dumb" deactivation will only be used if the stepsize falls below this quantity
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
const std::string _base_name
ComputeStressBase(const InputParameters &parameters)
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
Real _epp_tol
Tolerance on the plastic strain increment ("direction") constraint.
Real _min_stepsize
Minimum fraction of applied strain that may be applied during adaptive stepsizing.
const MaterialProperty< RankTwoTensor > * _curvature
The Cosserat curvature strain.
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
std::vector< Real > _cumulative_pm
the sum of the plastic multipliers over all the sub-steps.
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
MaterialProperty< RankTwoTensor > * _couple_stress
the Cosserat couple-stress
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
bool _cosserat
whether Cosserat mechanics should be used
MultiPlasticityDebugger(const MooseObject *moose_object)
const MaterialProperty< RankFourTensor > * _elastic_flexural_rigidity_tensor
The Cosserat elastic flexural rigidity tensor.
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
MaterialProperty< RankFourTensor > * _Jacobian_mult_couple
derivative of couple-stress w.r.t. curvature
const MaterialProperty< RankTwoTensor > & _total_strain_old
Old value of total strain (coming from ComputeIncrementalSmallStrain, for example) ...

Member Function Documentation

◆ activeCombinationNumber()

unsigned int ComputeMultiPlasticityStress::activeCombinationNumber ( const std::vector< bool > &  act)
protected

Definition at line 1570 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1571 {
1572  unsigned num = 0;
1573  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1574  if (act[surface])
1575  num += (1 << surface); // (1 << x) = 2^x
1576 
1577  return num;
1578 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ activeModelSurfaces()

void MultiPlasticityRawComponentAssembler::activeModelSurfaces ( int  model,
const std::vector< bool > &  active,
std::vector< unsigned int > &  active_surfaces_of_model 
)
protectedinherited

Returns the internal surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.

Parameters
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfaces_of_modelthe output

Definition at line 810 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), and MultiPlasticityRawComponentAssembler::yieldFunction().

814 {
815  active_surfaces_of_model.resize(0);
816  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
817  if (active[_surfaces_given_model[model][model_surface]])
818  active_surfaces_of_model.push_back(model_surface);
819 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ activeSurfaces()

void MultiPlasticityRawComponentAssembler::activeSurfaces ( int  model,
const std::vector< bool > &  active,
std::vector< unsigned int > &  active_surfaces 
)
protectedinherited

Returns the external surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.

Parameters
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfacesthe output

Definition at line 799 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

802 {
803  active_surfaces.resize(0);
804  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
805  if (active[_surfaces_given_model[model][model_surface]])
806  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
807 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ anyActiveSurfaces()

bool MultiPlasticityRawComponentAssembler::anyActiveSurfaces ( int  model,
const std::vector< bool > &  active 
)
protectedinherited

returns true if any internal surfaces of the given model are active according to 'active'

Definition at line 790 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), MultiPlasticityDebugger::checkSolution(), MultiPlasticityDebugger::dof_included(), MultiPlasticityLinearSystem::nrStep(), and residual2().

791 {
792  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
793  if (active[_surfaces_given_model[model][model_surface]])
794  return true;
795  return false;
796 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ applyKuhnTucker()

void ComputeMultiPlasticityStress::applyKuhnTucker ( const std::vector< Real > &  f,
const std::vector< Real > &  pm,
std::vector< bool > &  active 
)
protectedvirtual

Checks Kuhn-Tucker conditions, and alters "active" if appropriate.

Do not let the simplicity of this routine fool you! Explicitly: (1) checks that pm = 0 for all the f < 0. If not, then active is set to false for that constraint. This may be triggered if upon exit of the NR loops a constraint got deactivated due to linear dependence, and then f<0 and its pm>0. (2) checks that pm = 0 for all inactive constraints. This should always be true unless someone has screwed with the code. (3) if any pm < 0, then active is set to false for that constraint. This may be triggered if _deactivation_scheme!="optimized".

Parameters
fvalues of the active yield functions
pmvalues of all the plastic multipliers
activethe active constraints (true if active)
Returns
return false if any of the Kuhn-Tucker conditions were violated (and hence the set of active constraints was changed)

Definition at line 1309 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1312 {
1313  bool turned_off = false;
1314  unsigned ind = 0;
1315 
1316  // turn off all active surfaces that have f<0 and pm!=0
1317  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1318  {
1319  if (active[surface])
1320  {
1321  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1322  if (pm[surface] != 0)
1323  {
1324  turned_off = true;
1325  active[surface] = false;
1326  }
1327  }
1328  else if (pm[surface] != 0)
1329  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1330  "coding!!");
1331  }
1332 
1333  // if didn't turn off anything yet, turn off surface with minimum pm
1334  if (!turned_off)
1335  {
1336  int surface_to_turn_off = -1;
1337  Real min_pm = 0;
1338  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1339  if (pm[surface] < min_pm)
1340  {
1341  min_pm = pm[surface];
1342  surface_to_turn_off = surface;
1343  }
1344  if (surface_to_turn_off >= 0)
1345  active[surface_to_turn_off] = false;
1346  }
1347 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ buildActiveConstraints()

void MultiPlasticityRawComponentAssembler::buildActiveConstraints ( const std::vector< Real > &  f,
const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const RankFourTensor &  Eijkl,
std::vector< bool > &  act 
)
protectedvirtualinherited

Constructs a set of active constraints, given the yield functions, f.

This uses TensorMechanicsPlasticModel::activeConstraints to identify the active constraints for each model.

Parameters
fyield functions (should be _num_surfaces of these)
stressstress tensor
intnlinternal parameters
Eijklelasticity tensor (stress = Eijkl*strain)
[out]actthe set of active constraints (will be resized to _num_surfaces)

Definition at line 343 of file MultiPlasticityRawComponentAssembler.C.

Referenced by returnMap().

348 {
349  mooseAssert(f.size() == _num_surfaces,
350  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
351  << _num_surfaces << " surfaces");
352  mooseAssert(intnl.size() == _num_models,
353  "buildActiveConstraints called with intnl.size = "
354  << intnl.size() << " while there are " << _num_models << " models");
355 
356  if (_specialIC == "rock")
357  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
358  else if (_specialIC == "joint")
359  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
360  else // no specialIC
361  {
362  act.resize(0);
363  unsigned ind = 0;
364  for (unsigned model = 0; model < _num_models; ++model)
365  {
366  std::vector<Real> model_f(0);
367  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
368  model_f.push_back(f[ind++]);
369  std::vector<bool> model_act;
370  RankTwoTensor returned_stress;
371  _f[model]->activeConstraints(
372  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
373  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
374  act.push_back(model_act[model_surface]);
375  }
376  }
377 }
MooseEnum _specialIC
Allows initial set of active constraints to be chosen optimally.
void buildActiveConstraintsRock(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Rock" version Constructs a set of active constraints, given the yield functions, f...
void buildActiveConstraintsJoint(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Joint" version Constructs a set of active constraints, given the yield functions, f.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ buildDumbOrder()

void ComputeMultiPlasticityStress::buildDumbOrder ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
std::vector< unsigned int > &  dumb_order 
)
protected

Builds the order which "dumb" activation will take.

Parameters
stressstress to evaluate yield functions and derivatives at
intnlinternal parameters to evaluate yield functions and derivatives at
[out]dumb_orderdumb_order[0] will be the yield surface furthest away from (stress, intnl), dumb_order[1] will be the next yield surface, etc. The distance measure used is f/|df_dstress|. This array can then be fed into incrementDumb in order to first try the yield surfaces which are farthest away from the (stress, intnl).

Definition at line 1520 of file ComputeMultiPlasticityStress.C.

Referenced by changeScheme(), and returnMap().

1523 {
1524  if (dumb_order.size() != 0)
1525  return;
1526 
1527  std::vector<bool> act;
1528  act.assign(_num_surfaces, true);
1529 
1530  std::vector<Real> f;
1531  yieldFunction(stress, intnl, act, f);
1532  std::vector<RankTwoTensor> df_dstress;
1533  dyieldFunction_dstress(stress, intnl, act, df_dstress);
1534 
1535  typedef std::pair<Real, unsigned> pair_for_sorting;
1536  std::vector<pair_for_sorting> dist(_num_surfaces);
1537  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1538  {
1539  dist[surface].first = f[surface] / df_dstress[surface].L2norm();
1540  dist[surface].second = surface;
1541  }
1542  std::sort(dist.begin(), dist.end()); // sorted in ascending order of f/df_dstress
1543 
1544  dumb_order.resize(_num_surfaces);
1545  for (unsigned i = 0; i < _num_surfaces; ++i)
1546  dumb_order[i] = dist[_num_surfaces - 1 - i].second;
1547  // now dumb_order[0] is the surface with the greatest f/df_dstress
1548 }
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ calculateConstraints()

void MultiPlasticityLinearSystem::calculateConstraints ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankTwoTensor &  delta_dp,
std::vector< Real > &  f,
std::vector< RankTwoTensor > &  r,
RankTwoTensor &  epp,
std::vector< Real > &  ic,
const std::vector< bool > &  active 
)
protectedvirtualinherited

The constraints.

These are set to zero (or <=0 in the case of the yield functions) by the Newton-Raphson process, except in the case of linear-dependence which complicates things.

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]fActive yield function(s)
[out]rActive flow directions
[out]eppPlastic-strain increment constraint
[out]icActive internal-parameter constraint
activeThe active constraints.

Definition at line 227 of file MultiPlasticityLinearSystem.C.

Referenced by MultiPlasticityLinearSystem::calculateRHS(), and lineSearch().

237 {
238  // see comments at the start of .h file
239 
240  mooseAssert(intnl_old.size() == _num_models,
241  "Size of intnl_old is " << intnl_old.size()
242  << " which is incorrect in calculateConstraints");
243  mooseAssert(intnl.size() == _num_models,
244  "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
245  mooseAssert(pm.size() == _num_surfaces,
246  "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
247  mooseAssert(active.size() == _num_surfaces,
248  "Size of active is " << active.size()
249  << " which is incorrect in calculateConstraints");
250 
251  // yield functions
252  yieldFunction(stress, intnl, active, f);
253 
254  // flow directions and "epp"
255  flowPotential(stress, intnl, active, r);
256  epp = RankTwoTensor();
257  unsigned ind = 0;
258  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
259  if (active[surface])
260  epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
261  epp -= delta_dp;
262 
263  // internal constraints
264  std::vector<Real> h;
265  hardPotential(stress, intnl, active, h);
266  ic.resize(0);
267  ind = 0;
268  std::vector<unsigned int> active_surfaces;
269  std::vector<unsigned int>::iterator active_surface;
270  for (unsigned model = 0; model < _num_models; ++model)
271  {
272  activeSurfaces(model, active, active_surfaces);
273  if (active_surfaces.size() > 0)
274  {
275  // some surfaces are active in this model, so must form an internal constraint
276  ic.push_back(intnl[model] - intnl_old[model]);
277  for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
278  ++active_surface)
279  ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
280  // since it was constructed in the same
281  // manner
282  }
283  }
284 }
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
unsigned int _num_models
Number of plastic models for this material.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
void activeSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
Returns the external surface number(s) of the active surfaces of the given model This may be of size=...

◆ calculateJacobian()

void MultiPlasticityLinearSystem::calculateJacobian ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankFourTensor &  E_inv,
const std::vector< bool > &  active,
const std::vector< bool > &  deactivated_due_to_ld,
std::vector< std::vector< Real >> &  jac 
)
protectedvirtualinherited

d(rhs)/d(dof)

Definition at line 366 of file MultiPlasticityLinearSystem.C.

Referenced by MultiPlasticityDebugger::checkJacobian(), MultiPlasticityDebugger::checkSolution(), and MultiPlasticityLinearSystem::nrStep().

373 {
374  // see comments at the start of .h file
375 
376  mooseAssert(intnl.size() == _num_models,
377  "Size of intnl is " << intnl.size() << " which is incorrect in calculateJacobian");
378  mooseAssert(pm.size() == _num_surfaces,
379  "Size of pm is " << pm.size() << " which is incorrect in calculateJacobian");
380  mooseAssert(active.size() == _num_surfaces,
381  "Size of active is " << active.size() << " which is incorrect in calculateJacobian");
382  mooseAssert(deactivated_due_to_ld.size() == _num_surfaces,
383  "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
384  << " which is incorrect in calculateJacobian");
385 
386  unsigned ind = 0;
387  unsigned active_surface_ind = 0;
388 
389  std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
390  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
391  active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
392  unsigned num_active_surface = 0;
393  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
394  if (active_surface[surface])
395  num_active_surface++;
396 
397  std::vector<bool> active_model(
398  _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
399  for (unsigned model = 0; model < _num_models; ++model)
400  active_model[model] = anyActiveSurfaces(model, active_surface);
401 
402  unsigned num_active_model = 0;
403  for (unsigned model = 0; model < _num_models; ++model)
404  if (active_model[model])
405  num_active_model++;
406 
407  ind = 0;
408  std::vector<unsigned int> active_model_index(_num_models);
409  for (unsigned model = 0; model < _num_models; ++model)
410  if (active_model[model])
411  active_model_index[model] = ind++;
412  else
413  active_model_index[model] =
414  _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
415 
416  std::vector<RankTwoTensor> df_dstress;
417  dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
418 
419  std::vector<Real> df_dintnl;
420  dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
421 
422  std::vector<RankTwoTensor> r;
423  flowPotential(stress, intnl, active, r);
424 
425  std::vector<RankFourTensor> dr_dstress;
426  dflowPotential_dstress(stress, intnl, active, dr_dstress);
427 
428  std::vector<RankTwoTensor> dr_dintnl;
429  dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
430 
431  std::vector<Real> h;
432  hardPotential(stress, intnl, active, h);
433 
434  std::vector<RankTwoTensor> dh_dstress;
435  dhardPotential_dstress(stress, intnl, active, dh_dstress);
436 
437  std::vector<Real> dh_dintnl;
438  dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
439 
440  // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
441  RankFourTensor depp_dstress;
442  ind = 0;
443  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
444  if (active[surface]) // includes deactivated_due_to_ld
445  depp_dstress += pm[surface] * dr_dstress[ind++];
446  depp_dstress += E_inv;
447 
448  // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
449  std::vector<RankTwoTensor> depp_dpm;
450  depp_dpm.resize(num_active_surface);
451  ind = 0;
452  active_surface_ind = 0;
453  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
454  {
455  if (active[surface])
456  {
457  if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
458  // dofs in the NR
459  depp_dpm[active_surface_ind++] = r[ind];
460  ind++;
461  }
462  }
463 
464  // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
465  std::vector<RankTwoTensor> depp_dintnl;
466  depp_dintnl.assign(num_active_model, RankTwoTensor());
467  ind = 0;
468  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
469  {
470  if (active[surface])
471  {
472  unsigned int model_num = modelNumber(surface);
473  if (active_model[model_num]) // only include models with surfaces which are still active after
474  // deactivated_due_to_ld
475  depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
476  ind++;
477  }
478  }
479 
480  // df_dstress has been calculated above
481  // df_dpm is always zero
482  // df_dintnl has been calculated above, but only the active_surface+active_model stuff needs to be
483  // included in Jacobian: see below
484 
485  std::vector<RankTwoTensor> dic_dstress;
486  dic_dstress.assign(num_active_model, RankTwoTensor());
487  ind = 0;
488  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
489  {
490  if (active[surface])
491  {
492  unsigned int model_num = modelNumber(surface);
493  if (active_model[model_num]) // only include ic for models with active_surface (ie, if model
494  // only contains deactivated_due_to_ld don't include it)
495  dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
496  ind++;
497  }
498  }
499 
500  std::vector<std::vector<Real>> dic_dpm;
501  dic_dpm.resize(num_active_model);
502  ind = 0;
503  active_surface_ind = 0;
504  for (unsigned model = 0; model < num_active_model; ++model)
505  dic_dpm[model].assign(num_active_surface, 0);
506  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
507  {
508  if (active[surface])
509  {
510  if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
511  {
512  unsigned int model_num = modelNumber(surface);
513  // if (active_model[model_num]) // do not need this check as if the surface has
514  // active_surface, the model must be deemed active!
515  dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
516  active_surface_ind++;
517  }
518  ind++;
519  }
520  }
521 
522  std::vector<std::vector<Real>> dic_dintnl;
523  dic_dintnl.resize(num_active_model);
524  for (unsigned model = 0; model < num_active_model; ++model)
525  {
526  dic_dintnl[model].assign(num_active_model, 0);
527  dic_dintnl[model][model] = 1; // deriv wrt internal parameter
528  }
529  ind = 0;
530  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
531  {
532  if (active[surface])
533  {
534  unsigned int model_num = modelNumber(surface);
535  if (active_model[model_num]) // only the models that contain surfaces that are still active
536  // after deactivation_due_to_ld
537  dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
538  pm[surface] * dh_dintnl[ind];
539  ind++;
540  }
541  }
542 
543  unsigned int dim = 3;
544  unsigned int system_size =
545  6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
546  jac.resize(system_size);
547  for (unsigned i = 0; i < system_size; ++i)
548  jac[i].assign(system_size, 0);
549 
550  unsigned int row_num = 0;
551  unsigned int col_num = 0;
552  for (unsigned i = 0; i < dim; ++i)
553  for (unsigned j = 0; j <= i; ++j)
554  {
555  for (unsigned k = 0; k < dim; ++k)
556  for (unsigned l = 0; l <= k; ++l)
557  jac[col_num][row_num++] =
558  depp_dstress(i, j, k, l) +
559  (k != l ? depp_dstress(i, j, l, k)
560  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
561  for (unsigned surface = 0; surface < num_active_surface; ++surface)
562  jac[col_num][row_num++] = depp_dpm[surface](i, j);
563  for (unsigned a = 0; a < num_active_model; ++a)
564  jac[col_num][row_num++] = depp_dintnl[a](i, j);
565  row_num = 0;
566  col_num++;
567  }
568 
569  ind = 0;
570  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
571  if (active_surface[surface])
572  {
573  for (unsigned k = 0; k < dim; ++k)
574  for (unsigned l = 0; l <= k; ++l)
575  jac[col_num][row_num++] =
576  df_dstress[ind](k, l) +
577  (k != l ? df_dstress[ind](l, k)
578  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
579  for (unsigned beta = 0; beta < num_active_surface; ++beta)
580  jac[col_num][row_num++] = 0; // df_dpm
581  for (unsigned model = 0; model < _num_models; ++model)
582  if (active_model[model]) // only use df_dintnl for models in active_model
583  {
584  if (modelNumber(surface) == model)
585  jac[col_num][row_num++] = df_dintnl[ind];
586  else
587  jac[col_num][row_num++] = 0;
588  }
589  ind++;
590  row_num = 0;
591  col_num++;
592  }
593 
594  for (unsigned a = 0; a < num_active_model; ++a)
595  {
596  for (unsigned k = 0; k < dim; ++k)
597  for (unsigned l = 0; l <= k; ++l)
598  jac[col_num][row_num++] =
599  dic_dstress[a](k, l) +
600  (k != l ? dic_dstress[a](l, k)
601  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
602  for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
603  jac[col_num][row_num++] = dic_dpm[a][alpha];
604  for (unsigned b = 0; b < num_active_model; ++b)
605  jac[col_num][row_num++] = dic_dintnl[a][b];
606  row_num = 0;
607  col_num++;
608  }
609 
610  mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
611 }
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
virtual void dhardPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
The derivative of the active hardening potentials with respect to the active internal parameters...
virtual void dhardPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
The derivative of the active hardening potentials with respect to stress By assumption in the Userobj...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
unsigned int _num_models
Number of plastic models for this material.
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...

◆ calculateRHS()

void MultiPlasticityLinearSystem::calculateRHS ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankTwoTensor &  delta_dp,
std::vector< Real > &  rhs,
const std::vector< bool > &  active,
bool  eliminate_ld,
std::vector< bool > &  deactivated_due_to_ld 
)
protectedvirtualinherited

Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])

Note that the 'epp' components only contain the upper diagonal. These contain flow directions and plasticity-multipliers for all active surfaces, even the deactivated_due_to_ld surfaces. Note that the 'f' components only contain the active and not deactivated_due_to_ld surfaces Note that the 'ic' components only contain the internal constraints for models which contain active and not deactivated_due_to_ld surfaces. They contain hardening-potentials and plasticity-multipliers for the active surfaces, even the deactivated_due_to_ld surfaces

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]rhsthe rhs
activeThe active constraints.
eliminate_ldCheck for linear dependence of constraints and put the results into deactivated_due_to_ld. Usually this should be true, but for certain debug operations it should be false
[out]deactivated_due_to_ldconstraints deactivated due to linear-dependence of flow directions

Definition at line 287 of file MultiPlasticityLinearSystem.C.

Referenced by MultiPlasticityDebugger::checkSolution(), MultiPlasticityDebugger::fdJacobian(), and MultiPlasticityLinearSystem::nrStep().

296 {
297  // see comments at the start of .h file
298 
299  mooseAssert(intnl_old.size() == _num_models,
300  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
301  mooseAssert(intnl.size() == _num_models,
302  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
303  mooseAssert(pm.size() == _num_surfaces,
304  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
305  mooseAssert(active.size() == _num_surfaces,
306  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
307 
308  std::vector<Real> f; // the yield functions
309  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
310  std::vector<Real> ic; // the "internal constraints"
311 
312  std::vector<RankTwoTensor> r;
313  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
314 
315  if (eliminate_ld)
316  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
317  else
318  deactivated_due_to_ld.assign(_num_surfaces, false);
319 
320  std::vector<bool> active_not_deact(_num_surfaces);
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
323 
324  unsigned num_active_f = 0;
325  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
326  if (active_not_deact[surface])
327  num_active_f++;
328 
329  unsigned num_active_ic = 0;
330  for (unsigned model = 0; model < _num_models; ++model)
331  if (anyActiveSurfaces(model, active_not_deact))
332  num_active_ic++;
333 
334  unsigned int dim = 3;
335  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
336  // num_active_f comes from "f",
337  // num_active_f comes from "ic"
338 
339  rhs.resize(system_size);
340 
341  unsigned ind = 0;
342  for (unsigned i = 0; i < dim; ++i)
343  for (unsigned j = 0; j <= i; ++j)
344  rhs[ind++] = -epp(i, j);
345  unsigned active_surface = 0;
346  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
347  if (active[surface])
348  {
349  if (!deactivated_due_to_ld[surface])
350  rhs[ind++] = -f[active_surface];
351  active_surface++;
352  }
353  unsigned active_model = 0;
354  for (unsigned model = 0; model < _num_models; ++model)
355  if (anyActiveSurfaces(model, active))
356  {
357  if (anyActiveSurfaces(model, active_not_deact))
358  rhs[ind++] = -ic[active_model];
359  active_model++;
360  }
361 
362  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
363 }
virtual void eliminateLinearDependence(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs a number of singular-value decompositions to check for linear-dependence of the active direc...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.

◆ canAddConstraints()

bool ComputeMultiPlasticityStress::canAddConstraints ( const std::vector< bool > &  act,
const std::vector< Real > &  all_f 
)
protected

Definition at line 1011 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1013 {
1014  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1015  if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1016  return true;
1017  return false;
1018 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ canChangeScheme()

bool ComputeMultiPlasticityStress::canChangeScheme ( DeactivationSchemeEnum  current_deactivation_scheme,
bool  can_revert_to_dumb 
)
protected

Definition at line 1021 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1023 {
1024  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
1025  return true;
1026 
1027  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
1028  return true;
1029 
1030  if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
1031  can_revert_to_dumb)
1032  return true;
1033 
1034  if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
1035  can_revert_to_dumb)
1036  return true;
1037 
1038  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1039  can_revert_to_dumb)
1040  return true;
1041 
1042  return false;
1043 }
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme

◆ canIncrementDumb()

bool ComputeMultiPlasticityStress::canIncrementDumb ( int  dumb_iteration)
protected

Definition at line 1563 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1564 {
1565  // (1 << _num_surfaces) = 2^_num_surfaces
1566  return ((dumb_iteration + 1) < (1 << _num_surfaces));
1567 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ changeScheme()

void ComputeMultiPlasticityStress::changeScheme ( const std::vector< bool > &  initial_act,
bool  can_revert_to_dumb,
const RankTwoTensor &  initial_stress,
const std::vector< Real > &  intnl_old,
DeactivationSchemeEnum current_deactivation_scheme,
std::vector< bool > &  act,
int &  dumb_iteration,
std::vector< unsigned int > &  dumb_order 
)
protected

Definition at line 1046 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1054 {
1055  if (current_deactivation_scheme == optimized &&
1058  {
1059  current_deactivation_scheme = safe;
1060  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1061  act[surface] = initial_act[surface];
1062  }
1063  else if ((current_deactivation_scheme == safe &&
1066  can_revert_to_dumb) ||
1067  (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1068  can_revert_to_dumb))
1069  {
1070  current_deactivation_scheme = dumb;
1071  dumb_iteration = 0;
1072  buildDumbOrder(initial_stress, intnl_old, dumb_order);
1073  incrementDumb(dumb_iteration, dumb_order, act);
1074  }
1075 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void incrementDumb(int &dumb_iteration, const std::vector< unsigned int > &dumb_order, std::vector< bool > &act)
Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of...
void buildDumbOrder(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
Builds the order which "dumb" activation will take.
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme

◆ checkAdmissible()

bool ComputeMultiPlasticityStress::checkAdmissible ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
std::vector< Real > &  all_f 
)
protectedvirtual

Checks whether the yield functions are in the admissible region.

Parameters
stressstress
intnlinternal parameters
[out]all_fthe values of all the yield functions
Returns
return false if any yield functions exceed their tolerance

Definition at line 1267 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1270 {
1271  std::vector<bool> act;
1272  act.assign(_num_surfaces, true);
1273 
1274  yieldFunction(stress, intnl, act, all_f);
1275 
1276  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1277  if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1278  return false;
1279 
1280  return true;
1281 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ checkDerivatives()

void MultiPlasticityDebugger::checkDerivatives ( )
inherited

Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.

Definition at line 84 of file MultiPlasticityDebugger.C.

Referenced by initQpStatefulProperties(), and plasticStep().

85 {
86  Moose::err
87  << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
89 
90  std::vector<bool> act;
91  act.assign(_num_surfaces, true);
92 
93  Moose::err << "\ndyieldFunction_dstress. Relative L2 norms.\n";
94  std::vector<RankTwoTensor> df_dstress;
95  std::vector<RankTwoTensor> fddf_dstress;
98  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
99  {
100  Moose::err << "surface = " << surface << " Relative L2norm = "
101  << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
102  (df_dstress[surface] + fddf_dstress[surface]).L2norm()
103  << "\n";
104  Moose::err << "Coded:\n";
105  df_dstress[surface].print();
106  Moose::err << "Finite difference:\n";
107  fddf_dstress[surface].print();
108  }
109 
110  Moose::err << "\ndyieldFunction_dintnl.\n";
111  std::vector<Real> df_dintnl;
113  Moose::err << "Coded:\n";
114  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
115  Moose::err << df_dintnl[surface] << " ";
116  Moose::err << "\n";
117  std::vector<Real> fddf_dintnl;
119  Moose::err << "Finite difference:\n";
120  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
121  Moose::err << fddf_dintnl[surface] << " ";
122  Moose::err << "\n";
123 
124  Moose::err << "\ndflowPotential_dstress. Relative L2 norms.\n";
125  std::vector<RankFourTensor> dr_dstress;
126  std::vector<RankFourTensor> fddr_dstress;
129  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
130  {
131  Moose::err << "surface = " << surface << " Relative L2norm = "
132  << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
133  (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
134  << "\n";
135  Moose::err << "Coded:\n";
136  dr_dstress[surface].print();
137  Moose::err << "Finite difference:\n";
138  fddr_dstress[surface].print();
139  }
140 
141  Moose::err << "\ndflowPotential_dintnl. Relative L2 norms.\n";
142  std::vector<RankTwoTensor> dr_dintnl;
143  std::vector<RankTwoTensor> fddr_dintnl;
146  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
147  {
148  Moose::err << "surface = " << surface << " Relative L2norm = "
149  << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
150  (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
151  << "\n";
152  Moose::err << "Coded:\n";
153  dr_dintnl[surface].print();
154  Moose::err << "Finite difference:\n";
155  fddr_dintnl[surface].print();
156  }
157 }
void fddyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &df_dintnl)
The finite-difference derivative of yield function(s) with respect to internal parameter(s) ...
void fddyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &df_dstress)
The finite-difference derivative of yield function(s) with respect to stress.
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
virtual void fddflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &dr_dintnl)
The finite-difference derivative of the flow potentials with respect to internal parameters.
virtual void fddflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankFourTensor > &dr_dstress)
The finite-difference derivative of the flow potential(s) with respect to stress. ...
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
Real L2norm(const RankTwoTensor &r2tensor)
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...

◆ checkJacobian()

void MultiPlasticityDebugger::checkJacobian ( const RankFourTensor &  E_inv,
const std::vector< Real > &  intnl_old 
)
inherited

Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.

Definition at line 160 of file MultiPlasticityDebugger.C.

Referenced by computeQpStress(), and plasticStep().

162 {
163  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
165 
166  std::vector<bool> act;
167  act.assign(_num_surfaces, true);
168  std::vector<bool> deactivated_due_to_ld;
169  deactivated_due_to_ld.assign(_num_surfaces, false);
170 
171  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
172 
173  std::vector<std::vector<Real>> jac;
177  E_inv,
178  act,
179  deactivated_due_to_ld,
180  jac);
181 
182  std::vector<std::vector<Real>> fdjac;
184  intnl_old,
187  delta_dp,
188  E_inv,
189  false,
190  fdjac);
191 
192  Real L2_numer = 0;
193  Real L2_denom = 0;
194  for (unsigned row = 0; row < jac.size(); ++row)
195  for (unsigned col = 0; col < jac.size(); ++col)
196  {
197  L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
198  L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
199  }
200  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
201 
202  Moose::err << "\nHand-coded Jacobian:\n";
203  for (unsigned row = 0; row < jac.size(); ++row)
204  {
205  for (unsigned col = 0; col < jac.size(); ++col)
206  Moose::err << jac[row][col] << " ";
207  Moose::err << "\n";
208  }
209 
210  Moose::err << "Finite difference Jacobian:\n";
211  for (unsigned row = 0; row < fdjac.size(); ++row)
212  {
213  for (unsigned col = 0; col < fdjac.size(); ++col)
214  Moose::err << fdjac[row][col] << " ";
215  Moose::err << "\n";
216  }
217 }
virtual void fdJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
The Jacobian calculated using finite differences.
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

◆ checkKuhnTucker()

bool ComputeMultiPlasticityStress::checkKuhnTucker ( const std::vector< Real > &  f,
const std::vector< Real > &  pm,
const std::vector< bool > &  active 
)
protectedvirtual

Checks Kuhn-Tucker conditions, and alters "active" if appropriate.

Do not let the simplicity of this routine fool you! Explicitly: (1) checks that pm = 0 for all the f < 0. If not, then active is set to false for that constraint. This may be triggered if upon exit of the NR loops a constraint got deactivated due to linear dependence, and then f<0 and its pm>0. (2) checks that pm = 0 for all inactive constraints. This should always be true unless someone has screwed with the code. (3) if any pm < 0, then active is set to false for that constraint. This may be triggered if _deactivation_scheme!="optimized".

Parameters
fvalues of the active yield functions
pmvalues of all the plastic multipliers
activethe active constraints (true if active)
Returns
return false if any of the Kuhn-Tucker conditions were violated (and hence the set of active constraints was changed)

Definition at line 1284 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1287 {
1288  unsigned ind = 0;
1289  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1290  {
1291  if (active[surface])
1292  {
1293  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1294  if (pm[surface] != 0)
1295  return false;
1296  }
1297  else if (pm[surface] != 0)
1298  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1299  "coding!!");
1300  }
1301  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1302  if (pm[surface] < 0)
1303  return false;
1304 
1305  return true;
1306 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ checkSolution()

void MultiPlasticityDebugger::checkSolution ( const RankFourTensor &  E_inv)
inherited

Checks that Ax does equal b in the NR procedure.

Definition at line 366 of file MultiPlasticityDebugger.C.

Referenced by computeQpStress(), and plasticStep().

367 {
368  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
369  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
371 
372  std::vector<bool> act;
373  act.assign(_num_surfaces, true);
374  std::vector<bool> deactivated_due_to_ld;
375  deactivated_due_to_ld.assign(_num_surfaces, false);
376 
377  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
378 
379  std::vector<Real> orig_rhs;
384  delta_dp,
385  orig_rhs,
386  act,
387  true,
388  deactivated_due_to_ld);
389 
390  Moose::err << "\nb = ";
391  for (unsigned i = 0; i < orig_rhs.size(); ++i)
392  Moose::err << orig_rhs[i] << " ";
393  Moose::err << "\n\n";
394 
395  std::vector<std::vector<Real>> jac_coded;
399  E_inv,
400  act,
401  deactivated_due_to_ld,
402  jac_coded);
403 
404  Moose::err
405  << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
406  Moose::err
407  << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
408  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
409  "linear dependence.\n";
410  Moose::err << "Hand-coded Jacobian:\n";
411  for (unsigned row = 0; row < jac_coded.size(); ++row)
412  {
413  for (unsigned col = 0; col < jac_coded.size(); ++col)
414  Moose::err << jac_coded[row][col] << " ";
415  Moose::err << "\n";
416  }
417 
418  deactivated_due_to_ld.assign(_num_surfaces,
419  false); // this potentially gets changed by nrStep, below
420  RankTwoTensor dstress;
421  std::vector<Real> dpm;
422  std::vector<Real> dintnl;
427  E_inv,
428  delta_dp,
429  dstress,
430  dpm,
431  dintnl,
432  act,
433  deactivated_due_to_ld);
434 
435  std::vector<bool> active_not_deact(_num_surfaces);
436  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
437  active_not_deact[surface] = !deactivated_due_to_ld[surface];
438 
439  std::vector<Real> x;
440  x.assign(orig_rhs.size(), 0);
441  unsigned ind = 0;
442  for (unsigned i = 0; i < 3; ++i)
443  for (unsigned j = 0; j <= i; ++j)
444  x[ind++] = dstress(i, j);
445  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
446  if (active_not_deact[surface])
447  x[ind++] = dpm[surface];
448  for (unsigned model = 0; model < _num_models; ++model)
449  if (anyActiveSurfaces(model, active_not_deact))
450  x[ind++] = dintnl[model];
451 
452  mooseAssert(ind == orig_rhs.size(),
453  "Incorrect extracting of changes from NR solution in the "
454  "finite-difference checking of nrStep");
455 
456  Moose::err << "\nThis yields x =";
457  for (unsigned i = 0; i < orig_rhs.size(); ++i)
458  Moose::err << x[i] << " ";
459  Moose::err << "\n";
460 
461  std::vector<std::vector<Real>> jac_fd;
466  delta_dp,
467  E_inv,
468  true,
469  jac_fd);
470 
471  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
472  Moose::err << "in order to check that the solution is correct\n";
473  Moose::err << "Finite-difference Jacobian:\n";
474  for (unsigned row = 0; row < jac_fd.size(); ++row)
475  {
476  for (unsigned col = 0; col < jac_fd.size(); ++col)
477  Moose::err << jac_fd[row][col] << " ";
478  Moose::err << "\n";
479  }
480 
481  Real L2_numer = 0;
482  Real L2_denom = 0;
483  for (unsigned row = 0; row < jac_coded.size(); ++row)
484  for (unsigned col = 0; col < jac_coded.size(); ++col)
485  {
486  L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
487  L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
488  }
489  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
490  << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
491 
492  std::vector<Real> fd_times_x;
493  fd_times_x.assign(orig_rhs.size(), 0);
494  for (unsigned row = 0; row < orig_rhs.size(); ++row)
495  for (unsigned col = 0; col < orig_rhs.size(); ++col)
496  fd_times_x[row] += jac_fd[row][col] * x[col];
497 
498  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
499  for (unsigned i = 0; i < orig_rhs.size(); ++i)
500  Moose::err << fd_times_x[i] << " ";
501  Moose::err << "\n";
502  Moose::err << "Recall that b = \n";
503  for (unsigned i = 0; i < orig_rhs.size(); ++i)
504  Moose::err << orig_rhs[i] << " ";
505  Moose::err << "\n";
506 
507  L2_numer = 0;
508  L2_denom = 0;
509  for (unsigned i = 0; i < orig_rhs.size(); ++i)
510  {
511  L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
512  L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
513  }
514  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
515 }
virtual void fdJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
The Jacobian calculated using finite differences.
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.

◆ computeQpProperties()

void ComputeStressBase::computeQpProperties ( )
overrideprotectedvirtualinherited

Definition at line 51 of file ComputeStressBase.C.

52 {
54 
55  // Add in extra stress
56  _stress[_qp] += _extra_stress[_qp];
57 }
virtual void computeQpStress()=0
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankTwoTensor > & _extra_stress
Extra stress tensor.

◆ computeQpStress()

void ComputeMultiPlasticityStress::computeQpStress ( )
protectedvirtual

Implements ComputeStressBase.

Definition at line 217 of file ComputeMultiPlasticityStress.C.

218 {
219  // the following "_my" variables can get rotated by preReturnMap and postReturnMap
222  if (_cosserat)
223  {
224  _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
225  _my_curvature = (*_curvature)[_qp];
226  }
227 
228  if (_fspb_debug == "jacobian_and_linear_system")
229  {
230  // cannot do this at initQpStatefulProperties level since E_ijkl is not defined
231  checkJacobian(_elasticity_tensor[_qp].invSymm(), _intnl_old[_qp]);
232  checkSolution(_elasticity_tensor[_qp].invSymm());
233  mooseError("Finite-differencing completed. Exiting with no error");
234  }
235 
236  preReturnMap(); // do rotations to new frame if necessary
237 
238  unsigned int number_iterations;
239  bool linesearch_needed = false;
240  bool ld_encountered = false;
241  bool constraints_added = false;
242 
243  _cumulative_pm.assign(_num_surfaces, 0);
244  // try a "quick" return first - this can be purely elastic, or a customised plastic return defined
245  // by a TensorMechanicsPlasticXXXX UserObject
246  const bool found_solution = quickStep(rot(_stress_old[_qp]),
247  _stress[_qp],
248  _intnl_old[_qp],
249  _intnl[_qp],
250  _dummy_pm,
252  rot(_plastic_strain_old[_qp]),
253  _plastic_strain[_qp],
256  _yf[_qp],
257  number_iterations,
258  _Jacobian_mult[_qp],
260  true);
261 
262  // if not purely elastic or the customised stuff failed, do some plastic return
263  if (!found_solution)
265  _stress[_qp],
266  _intnl_old[_qp],
267  _intnl[_qp],
268  rot(_plastic_strain_old[_qp]),
269  _plastic_strain[_qp],
272  _yf[_qp],
273  number_iterations,
274  linesearch_needed,
275  ld_encountered,
276  constraints_added,
277  _Jacobian_mult[_qp]);
278 
279  if (_cosserat)
280  {
281  (*_couple_stress)[_qp] = (*_elastic_flexural_rigidity_tensor)[_qp] * _my_curvature;
282  (*_Jacobian_mult_couple)[_qp] = _my_flexural_rigidity_tensor;
283  }
284 
285  postReturnMap(); // rotate back from new frame if necessary
286 
287  _iter[_qp] = 1.0 * number_iterations;
288  _linesearch_needed[_qp] = linesearch_needed;
289  _ld_encountered[_qp] = ld_encountered;
290  _constraints_added[_qp] = constraints_added;
291 
292  // Update measures of strain
294  (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
295 
296  // Rotate the tensors to the current configuration
298  {
299  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
300  _elastic_strain[_qp] =
301  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
302  _plastic_strain[_qp] =
303  _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
304  }
305 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
bool _perform_finite_strain_rotations
whether to perform the rotations necessary in finite-strain simulations
MooseEnum _fspb_debug
none - don&#39;t do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment (coming from ComputeIncrementalSmallStrain, for example)
MaterialProperty< Real > & _constraints_added
Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankTwoTensor > & _strain_increment
strain increment (coming from ComputeIncrementalSmallStrain, for example)
const MaterialProperty< RankTwoTensor > & _stress_old
Old value of stress.
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Old value of elastic strain.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
MaterialProperty< Real > & _ld_encountered
Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true...
RankTwoTensor rot(const RankTwoTensor &tens)
MaterialProperty< std::vector< Real > > & _yf
yield functions
void checkSolution(const RankFourTensor &E_inv)
Checks that Ax does equal b in the NR procedure.
std::vector< Real > _dummy_pm
dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStre...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
const MaterialProperty< RankFourTensor > & _elasticity_tensor
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
std::vector< Real > _cumulative_pm
the sum of the plastic multipliers over all the sub-steps.
MaterialProperty< RankTwoTensor > & _elastic_strain
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
virtual bool quickStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, RankFourTensor &consistent_tangent_operator, const quickStep_called_from_t called_from, bool final_step)
Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined ...
bool _cosserat
whether Cosserat mechanics should be used
virtual bool plasticStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, RankFourTensor &consistent_tangent_operator)
performs a plastic step
void checkJacobian(const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.

◆ consistentTangentOperator()

RankFourTensor ComputeMultiPlasticityStress::consistentTangentOperator ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const RankFourTensor &  E_ijkl,
const std::vector< Real > &  pm_this_step,
const std::vector< Real > &  cumulative_pm 
)
protected

Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rate)

The computations performed depend upon _tangent_operator_type

Parameters
stressThe value of stress after the return map algorithm has converged
intnlThe internal parameters after the return map has converged
E_ijklThe elasticity tensor (in the case of no plasticity this is the jacobian)
pm_this_stepThe plastic multipliers coming from the final strain increment. In many cases these will be equal to cumulative_pm, but in the case where the returnMap algorithm had to be performed in multiple substeps of smaller applied strain increments, pm_this_step are just the plastic multipliers for the final application of the strain incrment
cumulative_pmThe plastic multipliers needed for this current Return (this is the sum of the plastic multipliers over all substeps if the strain increment was applied in small substeps)

Definition at line 1581 of file ComputeMultiPlasticityStress.C.

Referenced by quickStep(), and returnMap().

1586 {
1587 
1589  return E_ijkl;
1590 
1591  // Typically act_at_some_step = act, but it is possible
1592  // that when subdividing a strain increment, a surface
1593  // is only active for one sub-step
1594  std::vector<bool> act_at_some_step(_num_surfaces);
1595  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1596  act_at_some_step[surface] = (cumulative_pm[surface] > 0);
1597 
1598  // "act" might contain surfaces that are linearly dependent
1599  // with others. Only the plastic multipliers that are > 0
1600  // for this strain increment need to be varied to find
1601  // the consistent tangent operator
1602  std::vector<bool> act_vary(_num_surfaces);
1603  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1604  act_vary[surface] = (pm_this_step[surface] > 0);
1605 
1606  std::vector<RankTwoTensor> df_dstress;
1607  dyieldFunction_dstress(stress, intnl, act_vary, df_dstress);
1608  std::vector<Real> df_dintnl;
1609  dyieldFunction_dintnl(stress, intnl, act_vary, df_dintnl);
1610  std::vector<RankTwoTensor> r;
1611  flowPotential(stress, intnl, act_vary, r);
1612  std::vector<RankFourTensor> dr_dstress_at_some_step;
1613  dflowPotential_dstress(stress, intnl, act_at_some_step, dr_dstress_at_some_step);
1614  std::vector<RankTwoTensor> dr_dintnl_at_some_step;
1615  dflowPotential_dintnl(stress, intnl, act_at_some_step, dr_dintnl_at_some_step);
1616  std::vector<Real> h;
1617  hardPotential(stress, intnl, act_vary, h);
1618 
1619  unsigned ind1;
1620  unsigned ind2;
1621 
1622  // r_minus_stuff[alpha] = r[alpha] -
1623  // pm_cumulatve[gamma]*dr[gamma]_dintnl[a]_at_some_step*h[a][alpha], with alpha only being in
1624  // act_vary, but gamma being act_at_some_step
1625  std::vector<RankTwoTensor> r_minus_stuff;
1626  ind1 = 0;
1627  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1628  if (act_vary[surface1])
1629  {
1630  r_minus_stuff.push_back(r[ind1]);
1631  ind2 = 0;
1632  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1633  if (act_at_some_step[surface2])
1634  {
1635  if (modelNumber(surface1) == modelNumber(surface2))
1636  {
1637  r_minus_stuff.back() -=
1638  cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
1639  }
1640  ind2++;
1641  }
1642  ind1++;
1643  }
1644 
1645  unsigned int num_currently_active = 0;
1646  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1647  if (act_vary[surface])
1648  num_currently_active += 1;
1649 
1650  // zzz is a matrix in the form that can be easily
1651  // inverted by MatrixTools::inverse
1652  // Eg for num_currently_active = 3
1653  // (zzz[0] zzz[1] zzz[2])
1654  // (zzz[3] zzz[4] zzz[5])
1655  // (zzz[6] zzz[7] zzz[8])
1656  std::vector<PetscScalar> zzz;
1657  zzz.assign(num_currently_active * num_currently_active, 0.0);
1658 
1659  ind1 = 0;
1660  RankTwoTensor r2;
1661  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1662  if (act_vary[surface1])
1663  {
1664  ind2 = 0;
1665  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1666  if (act_vary[surface2])
1667  {
1668  r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
1669  zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
1670  if (modelNumber(surface1) == modelNumber(surface2))
1671  zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
1672  ind2++;
1673  }
1674  ind1++;
1675  }
1676 
1677  if (num_currently_active > 0)
1678  {
1679  // invert zzz, in place. if num_currently_active = 0 then zzz is not needed.
1680  try
1681  {
1682  MatrixTools::inverse(zzz, num_currently_active);
1683  }
1684  catch (const MooseException & e)
1685  {
1686  // in the very rare case of zzz being singular, just return the "elastic" tangent operator
1687  return E_ijkl;
1688  }
1689  }
1690 
1691  RankFourTensor strain_coeff = E_ijkl;
1692  ind1 = 0;
1693  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1694  if (act_vary[surface1])
1695  {
1696  RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1697  ind2 = 0;
1698  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1699  if (act_vary[surface2])
1700  {
1701  RankTwoTensor part2 = E_ijkl * df_dstress[ind2];
1702  for (unsigned i = 0; i < 3; i++)
1703  for (unsigned j = 0; j < 3; j++)
1704  for (unsigned k = 0; k < 3; k++)
1705  for (unsigned l = 0; l < 3; l++)
1706  strain_coeff(i, j, k, l) -=
1707  part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1708  ind2++;
1709  }
1710  ind1++;
1711  }
1712 
1714  return strain_coeff;
1715 
1716  RankFourTensor stress_coeff(RankFourTensor::initIdentitySymmetricFour);
1717 
1718  RankFourTensor part3;
1719  ind1 = 0;
1720  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1721  if (act_at_some_step[surface1])
1722  {
1723  part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
1724  ind1++;
1725  }
1726 
1727  stress_coeff += part3;
1728 
1729  part3 = part3.transposeMajor(); // this is because below i want df_dstress[ind2]*part3, and this
1730  // equals (part3.transposeMajor())*df_dstress[ind2]
1731 
1732  ind1 = 0;
1733  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1734  if (act_vary[surface1])
1735  {
1736  RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1737  ind2 = 0;
1738  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1739  if (act_vary[surface2])
1740  {
1741  RankTwoTensor part2 = part3 * df_dstress[ind2];
1742  for (unsigned i = 0; i < 3; i++)
1743  for (unsigned j = 0; j < 3; j++)
1744  for (unsigned k = 0; k < 3; k++)
1745  for (unsigned l = 0; l < 3; l++)
1746  stress_coeff(i, j, k, l) -=
1747  part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1748  ind2++;
1749  }
1750  ind1++;
1751  }
1752 
1753  // need to find the inverse of stress_coeff, but remember
1754  // stress_coeff does not have the symmetries commonly found
1755  // in tensor mechanics:
1756  // stress_coeff(i, j, k, l) = stress_coeff(j, i, k, l) = stress_coeff(i, j, l, k) !=
1757  // stress_coeff(k, l, i, j)
1758  // (note the final not-equals). We want s_inv, such that
1759  // s_inv(i, j, m, n)*stress_coeff(m, n, k, l) = (de_ik*de_jl + de_il*de_jk)/2
1760  // where de_ij = 1 if i=j, and 0 otherwise.
1761  RankFourTensor s_inv;
1762  try
1763  {
1764  s_inv = stress_coeff.invSymm();
1765  }
1766  catch (const MooseException & e)
1767  {
1768  return strain_coeff; // when stress_coeff is singular (perhaps for incompressible plasticity?)
1769  // return the "linear" tangent operator
1770  }
1771 
1772  return s_inv * strain_coeff;
1773 }
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...

◆ dflowPotential_dintnl()

void MultiPlasticityRawComponentAssembler::dflowPotential_dintnl ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  dr_dintnl 
)
protectedvirtualinherited

The derivative of the active flow potentials with respect to the active internal parameters The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha].

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dintnl"
[out]dr_dintnlthe derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 234 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityDebugger::checkDerivatives(), and consistentTangentOperator().

238 {
239  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
240  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
241 
242  dr_dintnl.resize(0);
243  std::vector<unsigned int> active_surfaces_of_model;
244  std::vector<unsigned int>::iterator active_surface;
245  std::vector<RankTwoTensor> model_dr_dintnl;
246  for (unsigned model = 0; model < _num_models; ++model)
247  {
248  activeModelSurfaces(model, active, active_surfaces_of_model);
249  if (active_surfaces_of_model.size() > 0)
250  {
251  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
252  for (active_surface = active_surfaces_of_model.begin();
253  active_surface != active_surfaces_of_model.end();
254  ++active_surface)
255  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
256  }
257  }
258 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ dflowPotential_dstress()

void MultiPlasticityRawComponentAssembler::dflowPotential_dstress ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankFourTensor > &  dr_dstress 
)
protectedvirtualinherited

The derivative of the active flow potential(s) with respect to stress.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dstress"
[out]dr_dstressthe derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 206 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityDebugger::checkDerivatives(), and consistentTangentOperator().

211 {
212  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
213  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
214 
215  dr_dstress.resize(0);
216  std::vector<unsigned int> active_surfaces_of_model;
217  std::vector<unsigned int>::iterator active_surface;
218  std::vector<RankFourTensor> model_dr_dstress;
219  for (unsigned model = 0; model < _num_models; ++model)
220  {
221  activeModelSurfaces(model, active, active_surfaces_of_model);
222  if (active_surfaces_of_model.size() > 0)
223  {
224  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
225  for (active_surface = active_surfaces_of_model.begin();
226  active_surface != active_surfaces_of_model.end();
227  ++active_surface)
228  dr_dstress.push_back(model_dr_dstress[*active_surface]);
229  }
230  }
231 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ dhardPotential_dintnl()

void MultiPlasticityRawComponentAssembler::dhardPotential_dintnl ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  dh_dintnl 
)
protectedvirtualinherited

The derivative of the active hardening potentials with respect to the active internal parameters.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dintnl"
[out]dh_dintnlthe derivatives. dh_dintnl[a][alpha][b] = dh[a][alpha]/dintnl[b]. Note that the userobjects assume that there is exactly one internal parameter per yield function, so the derivative is only nonzero for a=alpha=b, so that is all we calculate

Definition at line 316 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

320 {
321  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
322  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
323 
324  dh_dintnl.resize(0);
325  std::vector<unsigned int> active_surfaces_of_model;
326  std::vector<unsigned int>::iterator active_surface;
327  std::vector<Real> model_dh_dintnl;
328  for (unsigned model = 0; model < _num_models; ++model)
329  {
330  activeModelSurfaces(model, active, active_surfaces_of_model);
331  if (active_surfaces_of_model.size() > 0)
332  {
333  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
334  for (active_surface = active_surfaces_of_model.begin();
335  active_surface != active_surfaces_of_model.end();
336  ++active_surface)
337  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
338  }
339  }
340 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ dhardPotential_dstress()

void MultiPlasticityRawComponentAssembler::dhardPotential_dstress ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  dh_dstress 
)
protectedvirtualinherited

The derivative of the active hardening potentials with respect to stress By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only calculate those here.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dstress"
[out]dh_dstressthe derivative. dh_dstress[a](i, j) = dh[a]/dstress(k, l)

Definition at line 288 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

293 {
294  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
295  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
296 
297  dh_dstress.resize(0);
298  std::vector<unsigned int> active_surfaces_of_model;
299  std::vector<unsigned int>::iterator active_surface;
300  std::vector<RankTwoTensor> model_dh_dstress;
301  for (unsigned model = 0; model < _num_models; ++model)
302  {
303  activeModelSurfaces(model, active, active_surfaces_of_model);
304  if (active_surfaces_of_model.size() > 0)
305  {
306  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
307  for (active_surface = active_surfaces_of_model.begin();
308  active_surface != active_surfaces_of_model.end();
309  ++active_surface)
310  dh_dstress.push_back(model_dh_dstress[*active_surface]);
311  }
312  }
313 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ dyieldFunction_dintnl()

void MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  df_dintnl 
)
protectedvirtualinherited

The derivative of active yield function(s) with respect to their internal parameters (the user objects assume there is exactly one internal param per yield function)

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dintnl"
[out]df_dintnlthe derivatives. df_dstress[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 152 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityDebugger::checkDerivatives(), and consistentTangentOperator().

156 {
157  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
158  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
159 
160  df_dintnl.resize(0);
161  std::vector<unsigned int> active_surfaces_of_model;
162  std::vector<unsigned int>::iterator active_surface;
163  std::vector<Real> model_df_dintnl;
164  for (unsigned model = 0; model < _num_models; ++model)
165  {
166  activeModelSurfaces(model, active, active_surfaces_of_model);
167  if (active_surfaces_of_model.size() > 0)
168  {
169  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
170  for (active_surface = active_surfaces_of_model.begin();
171  active_surface != active_surfaces_of_model.end();
172  ++active_surface)
173  df_dintnl.push_back(model_df_dintnl[*active_surface]);
174  }
175  }
176 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ dyieldFunction_dstress()

void MultiPlasticityRawComponentAssembler::dyieldFunction_dstress ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  df_dstress 
)
protectedvirtualinherited

The derivative of the active yield function(s) with respect to stress.

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dstress"
[out]df_dstressthe derivative (or derivatives in the case of multisurface plasticity). df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Definition at line 124 of file MultiPlasticityRawComponentAssembler.C.

Referenced by buildDumbOrder(), MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityDebugger::checkDerivatives(), consistentTangentOperator(), and MultiPlasticityLinearSystem::eliminateLinearDependence().

129 {
130  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
131  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
132 
133  df_dstress.resize(0);
134  std::vector<unsigned int> active_surfaces_of_model;
135  std::vector<unsigned int>::iterator active_surface;
136  std::vector<RankTwoTensor> model_df_dstress;
137  for (unsigned model = 0; model < _num_models; ++model)
138  {
139  activeModelSurfaces(model, active, active_surfaces_of_model);
140  if (active_surfaces_of_model.size() > 0)
141  {
142  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
143  for (active_surface = active_surfaces_of_model.begin();
144  active_surface != active_surfaces_of_model.end();
145  ++active_surface)
146  df_dstress.push_back(model_df_dstress[*active_surface]);
147  }
148  }
149 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ flowPotential()

void MultiPlasticityRawComponentAssembler::flowPotential ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  r 
)
protectedvirtualinherited

The active flow potential(s) - one for each yield function.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active flow potentials are put into "r"
[out]rthe flow potential (flow potentials in the multi-surface case)

Definition at line 179 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), consistentTangentOperator(), MultiPlasticityDebugger::fddflowPotential_dintnl(), and MultiPlasticityDebugger::fddflowPotential_dstress().

183 {
184  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
185  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
186 
187  r.resize(0);
188  std::vector<unsigned int> active_surfaces_of_model;
189  std::vector<unsigned int>::iterator active_surface;
190  std::vector<RankTwoTensor> model_r;
191  for (unsigned model = 0; model < _num_models; ++model)
192  {
193  activeModelSurfaces(model, active, active_surfaces_of_model);
194  if (active_surfaces_of_model.size() > 0)
195  {
196  _f[model]->flowPotentialV(stress, intnl[model], model_r);
197  for (active_surface = active_surfaces_of_model.begin();
198  active_surface != active_surfaces_of_model.end();
199  ++active_surface)
200  r.push_back(model_r[*active_surface]);
201  }
202  }
203 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ hardPotential()

void MultiPlasticityRawComponentAssembler::hardPotential ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  h 
)
protectedvirtualinherited

The active hardening potentials (one for each internal parameter and for each yield function) by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part of model a, so we only calculate those here.

Parameters
stressthe stress at which to calculate the hardening potential
intnlvector of internal parameters
activeset of active constraints - only the active hardening potentials are put into "h"
[out]hthe hardening potentials. h[alpha] = hardening potential for yield fcn alpha (and, by the above assumption we know which hardening parameter, a, this belongs to)

Definition at line 261 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), and consistentTangentOperator().

265 {
266  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
267  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
268 
269  h.resize(0);
270  std::vector<unsigned int> active_surfaces_of_model;
271  std::vector<unsigned int>::iterator active_surface;
272  std::vector<Real> model_h;
273  for (unsigned model = 0; model < _num_models; ++model)
274  {
275  activeModelSurfaces(model, active, active_surfaces_of_model);
276  if (active_surfaces_of_model.size() > 0)
277  {
278  _f[model]->hardPotentialV(stress, intnl[model], model_h);
279  for (active_surface = active_surfaces_of_model.begin();
280  active_surface != active_surfaces_of_model.end();
281  ++active_surface)
282  h.push_back(model_h[*active_surface]);
283  }
284  }
285 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ incrementDumb()

void ComputeMultiPlasticityStress::incrementDumb ( int &  dumb_iteration,
const std::vector< unsigned int > &  dumb_order,
std::vector< bool > &  act 
)
protectedvirtual

Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of dumb_iteration == 1)

Parameters
[in,out]dumb_iterationUsed to set act bitwise - the "dumb" scheme tries all possible combinations of act until a successful return
[in]dumb_orderdumb_order dumb_order[0] will be the yield surface furthest away from (stress, intnl), dumb_order[1] will be the next yield surface, etc. The distance measure used is f/|df_dstress|. This array can then be fed into incrementDumb in order to first try the yield surfaces which are farthest away from the (stress, intnl).
[out]actactive constraints

Definition at line 1551 of file ComputeMultiPlasticityStress.C.

Referenced by changeScheme(), and returnMap().

1554 {
1555  dumb_iteration += 1;
1556  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1557  act[dumb_order[surface]] =
1558  (dumb_iteration &
1559  (1 << surface)); // returns true if the surface_th bit of dumb_iteration == 1
1560 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ initQpStatefulProperties()

void ComputeMultiPlasticityStress::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 187 of file ComputeMultiPlasticityStress.C.

188 {
190 
191  _plastic_strain[_qp].zero();
192 
193  _intnl[_qp].assign(_num_models, 0);
194 
195  _yf[_qp].assign(_num_surfaces, 0);
196 
197  _dummy_pm.assign(_num_surfaces, 0);
198 
199  _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
200  _linesearch_needed[_qp] = 0;
201  _ld_encountered[_qp] = 0;
202  _constraints_added[_qp] = 0;
203 
204  _n[_qp] = _n_input;
205 
206  if (_cosserat)
207  (*_couple_stress)[_qp].zero();
208 
209  if (_fspb_debug == "jacobian")
210  {
212  mooseError("Finite-differencing completed. Exiting with no error");
213  }
214 }
MooseEnum _fspb_debug
none - don&#39;t do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
MaterialProperty< Real > & _constraints_added
Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
MaterialProperty< Real > & _ld_encountered
Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true...
MaterialProperty< std::vector< Real > > & _yf
yield functions
void checkDerivatives()
Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
std::vector< Real > _dummy_pm
dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStre...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
RealVectorValue _n_input
the supplied transverse direction vector
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
unsigned int _num_models
Number of plastic models for this material.
virtual void initQpStatefulProperties() override
bool _cosserat
whether Cosserat mechanics should be used

◆ lineSearch()

bool ComputeMultiPlasticityStress::lineSearch ( Real &  nr_res2,
RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl,
std::vector< Real > &  pm,
const RankFourTensor &  E_inv,
RankTwoTensor &  delta_dp,
const RankTwoTensor &  dstress,
const std::vector< Real > &  dpm,
const std::vector< Real > &  dintnl,
std::vector< Real > &  f,
RankTwoTensor &  epp,
std::vector< Real > &  ic,
const std::vector< bool > &  active,
const std::vector< bool > &  deactivated_due_to_ld,
bool &  linesearch_needed 
)
protectedvirtual

Performs a line search.

Algorithm is taken straight from "Numerical Recipes". Given the changes dstress, dpm and dintnl provided by the nrStep routine, a line-search looks for an appropriate under-relaxation that reduces the residual-squared (nr_res2).

Most variables are input/output variables: they enter the function with their values at the start of the Newton step, and they exit the function with values attained after applying the under-relaxation

Parameters
[in,out]nr_res2The residual-squared
[out]stressThe stress after returning to the yield surface
intnl_oldThe internal variables at the previous "time" step
[in,out]intnlThe internal variables
[in,out]pmThe plasticity multiplier(s) (consistency parameter(s))
E_invinverse of the elasticity tensor
[in,out]delta_dpChange in plastic strain from start of "time" step to current configuration (plastic_strain - plastic_strain_old)
dstressChange in stress for a full Newton step
dpmChange in plasticity multiplier for a full Newton step
dintnlchange in internal parameter(s) for a full Newton step
[in,out]fYield function(s). In this routine, only the active constraints that are not deactivated_due_to_ld are contained in f.
[in,out]eppPlastic strain increment constraint
[in,out]icInternal constraint. In this routine, only the active constraints that are not deactivated_due_to_ld are contained in ic.
activeThe active constraints.
deactivated_due_to_ldTrue if a constraint has temporarily been made deactive due to linear dependence.
[out]linesearch_neededTrue if the full Newton-Raphson step was cut by the linesearch
Returns
true if successfully found a step that reduces the residual-squared

Definition at line 1387 of file ComputeMultiPlasticityStress.C.

Referenced by singleStep().

1403 {
1404  // Line search algorithm straight out of "Numerical Recipes"
1405 
1406  bool success =
1407  true; // return value: will be false if linesearch couldn't reduce the residual-squared
1408 
1409  // Aim is to decrease residual2
1410 
1411  Real lam = 1.0; // the line-search parameter: 1.0 is a full Newton step
1412  Real lam_min =
1413  1E-10; // minimum value of lam allowed - perhaps this should be dynamically calculated?
1414  Real f0 = nr_res2; // initial value of residual2
1415  Real slope = -2 * nr_res2; // "Numerical Recipes" uses -b*A*x, in order to check for roundoff, but
1416  // i hope the nrStep would warn if there were problems.
1417  Real tmp_lam; // cached value of lam used in quadratic & cubic line search
1418  Real f2 = nr_res2; // cached value of f = residual2 used in the cubic in the line search
1419  Real lam2 = lam; // cached value of lam used in the cubic in the line search
1420 
1421  // pm during the line-search
1422  std::vector<Real> ls_pm;
1423  ls_pm.resize(pm.size());
1424 
1425  // delta_dp during the line-search
1426  RankTwoTensor ls_delta_dp;
1427 
1428  // internal parameter during the line-search
1429  std::vector<Real> ls_intnl;
1430  ls_intnl.resize(intnl.size());
1431 
1432  // stress during the line-search
1433  RankTwoTensor ls_stress;
1434 
1435  // flow directions (not used in line search, but calculateConstraints returns this parameter)
1436  std::vector<RankTwoTensor> r;
1437 
1438  while (true)
1439  {
1440  // update the variables using this line-search parameter
1441  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1442  ls_pm[alpha] = pm[alpha] + dpm[alpha] * lam;
1443  ls_delta_dp = delta_dp - E_inv * dstress * lam;
1444  for (unsigned a = 0; a < intnl.size(); ++a)
1445  ls_intnl[a] = intnl[a] + dintnl[a] * lam;
1446  ls_stress = stress + dstress * lam;
1447 
1448  // calculate the new active yield functions, epp and active internal constraints
1449  calculateConstraints(ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1450 
1451  // calculate the new residual-squared
1452  nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1453 
1454  if (nr_res2 < f0 + 1E-4 * lam * slope)
1455  break;
1456  else if (lam < lam_min)
1457  {
1458  success = false;
1459  // restore plastic multipliers, yield functions, etc to original values
1460  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1461  ls_pm[alpha] = pm[alpha];
1462  ls_delta_dp = delta_dp;
1463  for (unsigned a = 0; a < intnl.size(); ++a)
1464  ls_intnl[a] = intnl[a];
1465  ls_stress = stress;
1467  ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1468  nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1469  break;
1470  }
1471  else if (lam == 1.0)
1472  {
1473  // model as a quadratic
1474  tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
1475  }
1476  else
1477  {
1478  // model as a cubic
1479  Real rhs1 = nr_res2 - f0 - lam * slope;
1480  Real rhs2 = f2 - f0 - lam2 * slope;
1481  Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1482  Real b =
1483  (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1484  if (a == 0)
1485  tmp_lam = -slope / (2.0 * b);
1486  else
1487  {
1488  Real disc = Utility::pow<2>(b) - 3 * a * slope;
1489  if (disc < 0)
1490  tmp_lam = 0.5 * lam;
1491  else if (b <= 0)
1492  tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
1493  else
1494  tmp_lam = -slope / (b + std::sqrt(disc));
1495  }
1496  if (tmp_lam > 0.5 * lam)
1497  tmp_lam = 0.5 * lam;
1498  }
1499  lam2 = lam;
1500  f2 = nr_res2;
1501  lam = std::max(tmp_lam, 0.1 * lam);
1502  }
1503 
1504  if (lam < 1.0)
1505  linesearch_needed = true;
1506 
1507  // assign the quantities found in the line-search
1508  // back to the originals
1509  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1510  pm[alpha] = ls_pm[alpha];
1511  delta_dp = ls_delta_dp;
1512  for (unsigned a = 0; a < intnl.size(); ++a)
1513  intnl[a] = ls_intnl[a];
1514  stress = ls_stress;
1515 
1516  return success;
1517 }
virtual Real residual2(const std::vector< Real > &pm, const std::vector< Real > &f, const RankTwoTensor &epp, const std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld)
The residual-squared.
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.

◆ modelNumber()

unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface)
protectedinherited

◆ nrStep()

void MultiPlasticityLinearSystem::nrStep ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankFourTensor &  E_inv,
const RankTwoTensor &  delta_dp,
RankTwoTensor &  dstress,
std::vector< Real > &  dpm,
std::vector< Real > &  dintnl,
const std::vector< bool > &  active,
std::vector< bool > &  deactivated_due_to_ld 
)
protectedvirtualinherited

Performs one Newton-Raphson step.

The purpose here is to find the changes, dstress, dpm and dintnl according to the Newton-Raphson procedure

Parameters
stressCurrent value of stress
intnl_oldThe internal variables at the previous "time" step
intnlCurrent value of the internal variables
pmCurrent value of the plasticity multipliers (consistency parameters)
E_invinverse of the elasticity tensor
delta_dpCurrent value of the plastic-strain increment (ie plastic_strain - plastic_strain_old)
[out]dstressThe change in stress for a full Newton step
[out]dpmThe change in all plasticity multipliers for a full Newton step
[out]dintnlThe change in all internal variables for a full Newton step
activeThe active constraints
[out]deactivated_due_to_ldThe constraints deactivated due to linear-dependence of the flow directions

Definition at line 614 of file MultiPlasticityLinearSystem.C.

Referenced by MultiPlasticityDebugger::checkSolution(), and singleStep().

625 {
626  // Calculate RHS and Jacobian
627  std::vector<Real> rhs;
628  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
629 
630  std::vector<std::vector<Real>> jac;
631  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
632 
633  // prepare for LAPACKgesv_ routine provided by PETSc
634  int system_size = rhs.size();
635 
636  std::vector<double> a(system_size * system_size);
637  // Fill in the a "matrix" by going down columns
638  unsigned ind = 0;
639  for (int col = 0; col < system_size; ++col)
640  for (int row = 0; row < system_size; ++row)
641  a[ind++] = jac[row][col];
642 
643  int nrhs = 1;
644  std::vector<int> ipiv(system_size);
645  int info;
646  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
647 
648  if (info != 0)
649  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
650  "routine returned with error code ",
651  info);
652 
653  // Extract the results back to dstress, dpm and dintnl
654  std::vector<bool> active_not_deact(_num_surfaces);
655  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
656  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
657 
658  unsigned int dim = 3;
659  ind = 0;
660 
661  for (unsigned i = 0; i < dim; ++i)
662  for (unsigned j = 0; j <= i; ++j)
663  dstress(i, j) = dstress(j, i) = rhs[ind++];
664  dpm.assign(_num_surfaces, 0);
665  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
666  if (active_not_deact[surface])
667  dpm[surface] = rhs[ind++];
668  dintnl.assign(_num_models, 0);
669  for (unsigned model = 0; model < _num_models; ++model)
670  if (anyActiveSurfaces(model, active_not_deact))
671  dintnl[model] = rhs[ind++];
672 
673  mooseAssert(static_cast<int>(ind) == system_size,
674  "Incorrect extracting of changes from NR solution in nrStep");
675 }
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.

◆ numberActive()

unsigned ComputeMultiPlasticityStress::numberActive ( const std::vector< bool > &  active)
protectedvirtual

counts the number of active constraints

Definition at line 1257 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap(), and singleStep().

1258 {
1259  unsigned num_active = 0;
1260  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1261  if (active[surface])
1262  num_active++;
1263  return num_active;
1264 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ outputAndCheckDebugParameters()

void MultiPlasticityDebugger::outputAndCheckDebugParameters ( )
inherited

Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly.

Definition at line 54 of file MultiPlasticityDebugger.C.

Referenced by MultiPlasticityDebugger::checkDerivatives(), MultiPlasticityDebugger::checkJacobian(), and MultiPlasticityDebugger::checkSolution().

55 {
56  Moose::err << "Debug Parameters are as follows\n";
57  Moose::err << "stress = \n";
58  _fspb_debug_stress.print();
59 
60  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
63  mooseError("The debug parameters have the wrong size\n");
64 
65  Moose::err << "plastic multipliers =\n";
66  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
67  Moose::err << _fspb_debug_pm[surface] << "\n";
68 
69  Moose::err << "internal parameters =\n";
70  for (unsigned model = 0; model < _num_models; ++model)
71  Moose::err << _fspb_debug_intnl[model] << "\n";
72 
73  Moose::err << "finite-differencing parameter for stress-changes:\n"
74  << _fspb_debug_stress_change << "\n";
75  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
76  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
77  Moose::err << _fspb_debug_pm_change[surface] << "\n";
78  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
79  for (unsigned model = 0; model < _num_models; ++model)
80  Moose::err << _fspb_debug_intnl_change[model] << "\n";
81 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

◆ plasticStep()

bool ComputeMultiPlasticityStress::plasticStep ( const RankTwoTensor &  stress_old,
RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl,
const RankTwoTensor &  plastic_strain_old,
RankTwoTensor &  plastic_strain,
const RankFourTensor &  E_ijkl,
const RankTwoTensor &  strain_increment,
std::vector< Real > &  yf,
unsigned int &  iterations,
bool &  linesearch_needed,
bool &  ld_encountered,
bool &  constraints_added,
RankFourTensor &  consistent_tangent_operator 
)
protectedvirtual

performs a plastic step

Parameters
stress_oldThe value of stress at the previous "time" step
[out]stressstress after returning to the yield surface
intnl_oldThe internal variables at the previous "time" step
[out]intnlinternal variables after returning to the yield surface
plastic_strain_oldThe value of plastic strain at the previous "time" step
[out]plastic_strainplastic_strain after returning to the yield surface
E_ijklThe elasticity tensor.
strain_incrementThe applied strain increment
[out]yfAll the yield functions at (stress, intnl)
[out]iterationsThe total number of Newton-Raphson iterations used
[out]linesearch_neededTrue if a linesearch was needed at any stage during the Newton-Raphson proceedure
[out]ld_encounteredTrue if a linear-dependence of the flow directions was encountered at any stage during the Newton-Raphson proceedure
[out]constraints_addedTrue if constraints were added into the active set at any stage during the Newton-Raphson proceedure
[out]consistent_tangent_operatorThe consistent tangent operator d(stress_rate)/d(strain_rate)
Returns
true if the (stress, intnl) are admissible. Otherwise, if _ignore_failures==true, the output variables will be the best admissible ones found during the return-map. Otherwise, if _ignore_failures==false, this routine will perform some finite-diference checks and call mooseError

the idea in the following is to potentially subdivide the strain increment into smaller fractions, of size "step_size". First step_size = 1 is tried, and if that fails then 0.5 is tried, then 0.25, etc. As soon as a step is successful, its results are put into the "good" variables, which are used if a subsequent step fails. If >= 2 consecutive steps are successful, the step_size is increased by 1.2

The point of all this is that I hope that the time-step for the entire mesh need not be cut if there are only a few "bad" quadpoints where the return-map is difficult

Definition at line 464 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

478 {
494  bool return_successful = false;
495 
496  // total number of Newton-Raphson iterations used
497  unsigned int iter = 0;
498 
499  Real step_size = 1.0;
500  Real time_simulated = 0.0;
501 
502  // the "good" variables hold the latest admissible stress
503  // and internal parameters.
504  RankTwoTensor stress_good = stress_old;
505  RankTwoTensor plastic_strain_good = plastic_strain_old;
506  std::vector<Real> intnl_good(_num_models);
507  for (unsigned model = 0; model < _num_models; ++model)
508  intnl_good[model] = intnl_old[model];
509  std::vector<Real> yf_good(_num_surfaces);
510 
511  // Following is necessary because I want strain_increment to be "const"
512  // but I also want to be able to subdivide an initial_stress
513  RankTwoTensor this_strain_increment = strain_increment;
514 
515  RankTwoTensor dep = step_size * this_strain_increment;
516 
517  _cumulative_pm.assign(_num_surfaces, 0);
518 
519  unsigned int num_consecutive_successes = 0;
520  while (time_simulated < 1.0 && step_size >= _min_stepsize)
521  {
522  iter = 0;
523  return_successful = returnMap(stress_good,
524  stress,
525  intnl_good,
526  intnl,
527  plastic_strain_good,
528  plastic_strain,
529  E_ijkl,
530  dep,
531  yf,
532  iter,
533  step_size <= _max_stepsize_for_dumb,
534  linesearch_needed,
535  ld_encountered,
536  constraints_added,
537  time_simulated + step_size >= 1,
538  consistent_tangent_operator,
540  iterations += iter;
541 
542  if (return_successful)
543  {
544  num_consecutive_successes += 1;
545  time_simulated += step_size;
546 
547  if (time_simulated < 1.0) // this condition is just for optimization: if time_simulated=1 then
548  // the "good" quantities are no longer needed
549  {
550  stress_good = stress;
551  plastic_strain_good = plastic_strain;
552  for (unsigned model = 0; model < _num_models; ++model)
553  intnl_good[model] = intnl[model];
554  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
555  yf_good[surface] = yf[surface];
556  if (num_consecutive_successes >= 2)
557  step_size *= 1.2;
558  }
559  step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
560  }
561  else
562  {
563  step_size *= 0.5;
564  num_consecutive_successes = 0;
565  stress = stress_good;
566  plastic_strain = plastic_strain_good;
567  for (unsigned model = 0; model < _num_models; ++model)
568  intnl[model] = intnl_good[model];
569  yf.resize(_num_surfaces); // might have excited with junk
570  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
571  yf[surface] = yf_good[surface];
572  dep = step_size * this_strain_increment;
573  }
574  }
575 
576  if (!return_successful)
577  {
578  if (_ignore_failures)
579  {
580  stress = stress_good;
581  plastic_strain = plastic_strain_good;
582  for (unsigned model = 0; model < _num_models; ++model)
583  intnl[model] = intnl_good[model];
584  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
585  yf[surface] = yf_good[surface];
586  }
587  else
588  {
589  Moose::out << "After reducing the stepsize to " << step_size
590  << " with original strain increment with L2norm " << this_strain_increment.L2norm()
591  << " the returnMap algorithm failed\n";
592 
593  _fspb_debug_stress = stress_good + E_ijkl * dep;
594  _fspb_debug_pm.assign(
596  1); // this is chosen arbitrarily - please change if a more suitable value occurs to you!
598  for (unsigned model = 0; model < _num_models; ++model)
599  _fspb_debug_intnl[model] = intnl_good[model];
603  mooseError("Exiting\n");
604  }
605  }
606 
607  return return_successful;
608 }
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
bool _ignore_failures
Even if the returnMap fails, return the best values found for stress and internal parameters...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void checkSolution(const RankFourTensor &E_inv)
Checks that Ax does equal b in the NR procedure.
void checkDerivatives()
Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
Real _max_stepsize_for_dumb
"dumb" deactivation will only be used if the stepsize falls below this quantity
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
virtual bool returnMap(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &f, unsigned int &iter, bool can_revert_to_dumb, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, bool final_step, RankFourTensor &consistent_tangent_operator, std::vector< Real > &cumulative_pm)
Implements the return map.
unsigned int _num_models
Number of plastic models for this material.
Real _min_stepsize
Minimum fraction of applied strain that may be applied during adaptive stepsizing.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
std::vector< Real > _cumulative_pm
the sum of the plastic multipliers over all the sub-steps.
void checkJacobian(const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.

◆ postReturnMap()

void ComputeMultiPlasticityStress::postReturnMap ( )
protectedvirtual

Definition at line 335 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

336 {
337  if (_n_supplied)
338  {
339  // this is a rotation matrix that will rotate "z" axis back to _n
340  _rot = _rot.transpose();
341 
342  // rotate the tensors back to original frame where _n is correctly oriented
343  _my_elasticity_tensor.rotate(_rot);
344  _Jacobian_mult[_qp].rotate(_rot);
345  _my_strain_increment.rotate(_rot);
346  _stress[_qp].rotate(_rot);
347  _plastic_strain[_qp].rotate(_rot);
348  if (_cosserat)
349  {
351  (*_Jacobian_mult_couple)[_qp].rotate(_rot);
352  _my_curvature.rotate(_rot);
353  (*_couple_stress)[_qp].rotate(_rot);
354  }
355 
356  // Rotate n by _rotation_increment
357  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
358  {
359  _n[_qp](i) = 0;
360  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
361  _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
362  }
363  }
364 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment (coming from ComputeIncrementalSmallStrain, for example)
const MaterialProperty< RealVectorValue > & _n_old
old value of transverse direction
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< RankTwoTensor > & _stress
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
bool _n_supplied
User supplied the transverse direction vector.
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
bool _cosserat
whether Cosserat mechanics should be used

◆ preReturnMap()

void ComputeMultiPlasticityStress::preReturnMap ( )
protectedvirtual

Definition at line 316 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

317 {
318  if (_n_supplied)
319  {
320  // this is a rotation matrix that will rotate _n to the "z" axis
321  _rot = RotationMatrix::rotVecToZ(_n[_qp]);
322 
323  // rotate the tensors to this frame
324  _my_elasticity_tensor.rotate(_rot);
325  _my_strain_increment.rotate(_rot);
326  if (_cosserat)
327  {
329  _my_curvature.rotate(_rot);
330  }
331  }
332 }
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
bool _n_supplied
User supplied the transverse direction vector.
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
bool _cosserat
whether Cosserat mechanics should be used

◆ quickStep()

bool ComputeMultiPlasticityStress::quickStep ( const RankTwoTensor &  stress_old,
RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl,
std::vector< Real > &  pm,
std::vector< Real > &  cumulative_pm,
const RankTwoTensor &  plastic_strain_old,
RankTwoTensor &  plastic_strain,
const RankFourTensor &  E_ijkl,
const RankTwoTensor &  strain_increment,
std::vector< Real > &  yf,
unsigned int &  iterations,
RankFourTensor &  consistent_tangent_operator,
const quickStep_called_from_t  called_from,
bool  final_step 
)
protectedvirtual

Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined through the TensorMechanicsPlasticXXXX.returnMap functions.

Parameters
stress_oldThe value of stress at the previous "time" step
[out]stressIf returnvalue=true then this is the returned value of stress. Otherwise, this is undefined
intnl_oldThe internal variables at the previous "time" step
[out]intnlIf returnvalue=true then this is the value of the internal parameters after returning. Otherwise, this is undefined
[out]pmIf returnvalue=true, this is the plastic multipliers needed to bring aout the return. Otherwise, this is undefined
[in/out]cumulative_pm If returnvalue=true, this is cumulative plastic multipliers, updated with pm. Otherwise, this is untouched by the algorithm
plastic_strain_oldThe value of plastic strain at the previous "time" step
[out]plastic_strainIf returnvalue=true, this is the new plastic strain. Otherwise it is set to plastic_strain_old
E_ijklThe elasticity tensor.
strain_incrementThe applied strain increment
[out]yfIf returnvalue=true, then all the yield functions at (stress, intnl). Otherwise, all the yield functions at (stress_old, intnl_old)
[out]iterationsNumber of NR iterations used, which is always zero in the current implementation.
called_fromThis can be called from computeQpStress, in which case it can actually provde an answer to the returnmap algorithm, or from returnMap in which case it is probably only providing an answer to a particular subdivision of the returnmap algorithm. The consistent tangent operator is calculated idfferently in each case
final_stepThe consistent tangent operator is calculated if this is true
[out]consistent_tangent_operatorIf final_step==true and returnvalue=true, then this is the consistent tangent operator d(stress_rate)/d(strain_rate). Otherwise it is undefined.
Returns
true if the (stress, intnl) are admissible, in which case the (stress_old, intnl_old) could have been admissible, or exactly one of the plastic models successfully used its custom returnMap function to provide the returned (stress, intnl) values and all other plastic models are admissible at that configuration. Or, false, then (stress_old, intnl_old) is not admissible according to >=1 plastic model and the custom returnMap functions failed in some way.

Definition at line 367 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress(), and returnMap().

382 {
383  iterations = 0;
384 
385  unsigned num_plastic_returns;
386  RankTwoTensor delta_dp;
387 
388  // the following does the customized returnMap algorithm
389  // for all the plastic models.
390  unsigned custom_model = 0;
391  bool successful_return = returnMapAll(stress_old + E_ijkl * strain_increment,
392  intnl_old,
393  E_ijkl,
394  _epp_tol,
395  stress,
396  intnl,
397  pm,
398  cumulative_pm,
399  delta_dp,
400  yf,
401  num_plastic_returns,
402  custom_model);
403 
404  // the following updates the plastic_strain, when necessary
405  // and calculates the consistent_tangent_operator, when necessary
406  if (num_plastic_returns == 0)
407  {
408  // if successful_return = true, then a purely elastic step
409  // if successful_return = false, then >=1 plastic model is in
410  // inadmissible zone and failed to return using its customized
411  // returnMap function.
412  // In either case:
413  plastic_strain = plastic_strain_old;
414  if (successful_return && final_step)
415  {
416  if (called_from == computeQpStress_function)
417  consistent_tangent_operator = E_ijkl;
418  else // cannot necessarily use E_ijkl since different plastic models may have been active
419  // during other substeps
420  consistent_tangent_operator =
421  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
422  }
423  return successful_return;
424  }
425  else if (num_plastic_returns == 1 && successful_return)
426  {
427  // one model has successfully completed its custom returnMap algorithm
428  // and the other models have signalled they are elastic at
429  // the trial stress
430  plastic_strain = plastic_strain_old + delta_dp;
431  if (final_step)
432  {
433  if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
434  {
436  consistent_tangent_operator = E_ijkl;
437  else
438  {
439  std::vector<Real> custom_model_pm;
440  for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
441  custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
442  consistent_tangent_operator =
443  _f[custom_model]->consistentTangentOperator(stress_old + E_ijkl * strain_increment,
444  intnl_old[custom_model],
445  stress,
446  intnl[custom_model],
447  E_ijkl,
448  custom_model_pm);
449  }
450  }
451  else // cannot necessarily use the custom consistentTangentOperator since different plastic
452  // models may have been active during other substeps or the custom model says not to use
453  // its custom CTO algorithm
454  consistent_tangent_operator =
455  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
456  }
457  return true;
458  }
459  else // presumably returnMapAll is incorrectly coded!
460  mooseError("ComputeMultiPlasticityStress::quickStep should not get here!");
461 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
bool returnMapAll(const RankTwoTensor &trial_stress, const std::vector< Real > &intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &stress, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, RankTwoTensor &delta_dp, std::vector< Real > &yf, unsigned &num_successful_plastic_returns, unsigned &custom_model)
Performs a returnMap for each plastic model using their inbuilt returnMap functions.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
Real _epp_tol
Tolerance on the plastic strain increment ("direction") constraint.
RankFourTensor consistentTangentOperator(const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &pm_this_step, const std::vector< Real > &cumulative_pm)
Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rat...

◆ reinstateLinearDependentConstraints()

bool ComputeMultiPlasticityStress::reinstateLinearDependentConstraints ( std::vector< bool > &  deactivated_due_to_ld)
protectedvirtual

makes all deactivated_due_to_ld false, and if >0 of them were initially true, returns true

Definition at line 1244 of file ComputeMultiPlasticityStress.C.

Referenced by singleStep().

1246 {
1247  bool reinstated_actives = false;
1248  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1249  if (deactivated_due_to_ld[surface])
1250  reinstated_actives = true;
1251 
1252  deactivated_due_to_ld.assign(_num_surfaces, false);
1253  return reinstated_actives;
1254 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.

◆ residual2()

Real ComputeMultiPlasticityStress::residual2 ( const std::vector< Real > &  pm,
const std::vector< Real > &  f,
const RankTwoTensor &  epp,
const std::vector< Real > &  ic,
const std::vector< bool > &  active,
const std::vector< bool > &  deactivated_due_to_ld 
)
protectedvirtual

The residual-squared.

Parameters
pmthe plastic multipliers for all constraints
fthe active yield function(s) (not including the ones that are deactivated_due_to_ld)
eppthe plastic strain increment constraint
icthe active internal constraint(s) (not including the ones that are deactivated_due_to_ld)
activetrue if constraint is active
deactivated_due_to_ldtrue if constraint has been temporarily deactivated due to linear dependence of flow directions

Definition at line 1350 of file ComputeMultiPlasticityStress.C.

Referenced by lineSearch(), and singleStep().

1356 {
1357  Real nr_res2 = 0;
1358  unsigned ind = 0;
1359 
1360  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1361  if (active[surface])
1362  {
1363  if (!deactivated_due_to_ld[surface])
1364  {
1365  if (!(pm[surface] == 0 && f[ind] <= 0))
1366  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1367  }
1368  else if (deactivated_due_to_ld[surface] && f[ind] > 0)
1369  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1370  ind++;
1371  }
1372 
1373  nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
1374 
1375  std::vector<bool> active_not_deact(_num_surfaces);
1376  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1377  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1378  ind = 0;
1379  for (unsigned model = 0; model < _num_models; ++model)
1380  if (anyActiveSurfaces(model, active_not_deact))
1381  nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
1382 
1383  return nr_res2;
1384 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.
Real _epp_tol
Tolerance on the plastic strain increment ("direction") constraint.

◆ returnMap()

bool ComputeMultiPlasticityStress::returnMap ( const RankTwoTensor &  stress_old,
RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl,
const RankTwoTensor &  plastic_strain_old,
RankTwoTensor &  plastic_strain,
const RankFourTensor &  E_ijkl,
const RankTwoTensor &  strain_increment,
std::vector< Real > &  f,
unsigned int &  iter,
bool  can_revert_to_dumb,
bool &  linesearch_needed,
bool &  ld_encountered,
bool &  constraints_added,
bool  final_step,
RankFourTensor &  consistent_tangent_operator,
std::vector< Real > &  cumulative_pm 
)
protectedvirtual

Implements the return map.

Note that this algorithm doesn't do any rotations. In order to find the final stress and plastic_strain must be rotated using _rotation_increment. This is usually done in computeQpStress

Parameters
stress_oldThe value of stress at the previous "time" step
[out]stressThe stress after returning to the yield surface
intnl_oldThe internal variables at the previous "time" step
[out]intnlAll the internal variables after returning to the yield surface
plastic_strain_oldThe value of plastic strain at the previous "time" step
[out]plastic_strainThe value of plastic strain after returning to the yield surface
E_ijklThe elasticity tensor. If no plasiticity then stress = stress_old + E_ijkl*strain_increment
strain_incrementThe applied strain increment
[out]fAll the yield functions after returning to the yield surface
[out]iterThe number of Newton-Raphson iterations used
can_revert_to_dumbIf the _deactivation_scheme is set to revert to dumb, it will only be allowed to do so if this parameter is true
[out]linesearch_neededTrue if a linesearch was needed at any stage during the Newton-Raphson proceedure
[out]ld_encounteredTrue if a linear-dependence of the flow directions was encountered at any stage during the Newton-Raphson proceedure
[out]constraints_addedTrue if constraints were added into the active set at any stage during the Newton-Raphson proceedure
final_stepEach strain increment may be decomposed into a sum of smaller increments if the return-map algorithm fails. This flag indicates whether this is the last application of incremental strain
[out]consistent_tangent_operatorThe consistent tangent operator d(stress_rate)/d(strain_rate). This is only output if final_step=true, and the return value of returnMap is also true.
[in,out]cumulative_pmUpon input: the plastic multipliers before the return map. Upon output: the plastic multipliers after this return map, if the return map was successful
Returns
true if the stress was successfully returned to the yield surface

Definition at line 611 of file ComputeMultiPlasticityStress.C.

Referenced by plasticStep().

628 {
629 
630  // The "consistency parameters" (plastic multipliers)
631  // Change in plastic strain in this timestep = pm*flowPotential
632  // Each pm must be non-negative
633  std::vector<Real> pm;
634  pm.assign(_num_surfaces, 0.0);
635 
636  bool successful_return = quickStep(stress_old,
637  stress,
638  intnl_old,
639  intnl,
640  pm,
641  cumulative_pm,
642  plastic_strain_old,
643  plastic_strain,
644  E_ijkl,
645  strain_increment,
646  f,
647  iter,
648  consistent_tangent_operator,
650  final_step);
651 
652  if (successful_return)
653  return successful_return;
654 
655  // Here we know that the trial stress and intnl_old
656  // is inadmissible, and we have to return from those values
657  // value to the yield surface. There are three
658  // types of constraints we have to satisfy, listed
659  // below, and calculated in calculateConstraints(...)
660  // f<=0, epp=0, ic=0 (up to tolerances), and these are
661  // guaranteed to hold if nr_res2<=0.5
662  //
663  // Kuhn-Tucker conditions must also be satisfied
664  // These are:
665  // pm>=0, which may not hold upon exit of the NR loops
666  // due to _deactivation_scheme!=optimized;
667  // pm*f=0 (up to tolerances), which may not hold upon exit
668  // of the NR loops if a constraint got deactivated
669  // due to linear dependence, and then f<0, and its pm>0
670 
671  // Plastic strain constraint, L2 norm must be zero (up to a tolerance)
672  RankTwoTensor epp;
673 
674  // Yield function constraint passed to this function as
675  // std::vector<Real> & f
676  // Each yield function must be <= 0 (up to tolerance)
677  // Note that only the constraints that are active will be
678  // contained in f until the final few lines of returnMap
679 
680  // Internal constraint(s), must be zero (up to a tolerance)
681  // Note that only the constraints that are active will be
682  // contained in ic.
683  std::vector<Real> ic;
684 
685  // Record the stress before Newton-Raphson in case of failure-and-restart
686  RankTwoTensor initial_stress = stress;
687 
688  iter = 0;
689 
690  // Initialize the set of active constraints
691  // At this stage, the active constraints are
692  // those that exceed their _f_tol
693  // active constraints.
694  std::vector<bool> act;
695  buildActiveConstraints(f, stress, intnl, E_ijkl, act);
696 
697  // Inverse of E_ijkl (assuming symmetric)
698  RankFourTensor E_inv = E_ijkl.invSymm();
699 
700  // convenience variable that holds the change in plastic strain incurred during the return
701  // delta_dp = plastic_strain - plastic_strain_old
702  // delta_dp = E^{-1}*(initial_stress - stress), where initial_stress = E*(strain -
703  // plastic_strain_old)
704  RankTwoTensor delta_dp = RankTwoTensor();
705 
706  // whether single step was successful (whether line search was successful, and whether turning off
707  // constraints was successful)
708  bool single_step_success = true;
709 
710  // deactivation scheme
712 
713  // For complicated deactivation schemes we have to record the initial active set
714  std::vector<bool> initial_act;
715  initial_act.resize(_num_surfaces);
719  {
720  // if "optimized" fails we can change the deactivation scheme to "safe", etc
721  deact_scheme = optimized;
722  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
723  initial_act[surface] = act[surface];
724  }
725 
727  deact_scheme = safe;
728 
729  // For "dumb" deactivation, the active set takes all combinations until a solution is found
730  int dumb_iteration = 0;
731  std::vector<unsigned int> dumb_order;
732 
733  if (_deactivation_scheme == dumb ||
734  (_deactivation_scheme == optimized_to_safe_to_dumb && can_revert_to_dumb) ||
735  (_deactivation_scheme == safe_to_dumb && can_revert_to_dumb) ||
736  (_deactivation_scheme == optimized_to_dumb && can_revert_to_dumb))
737  buildDumbOrder(stress, intnl, dumb_order);
738 
739  if (_deactivation_scheme == dumb)
740  {
741  incrementDumb(dumb_iteration, dumb_order, act);
742  yieldFunction(stress, intnl, act, f);
743  }
744 
745  // To avoid any re-trials of "act" combinations that
746  // we've already tried and rejected, i record the
747  // combinations in actives_tried
748  std::set<unsigned int> actives_tried;
749  actives_tried.insert(activeCombinationNumber(act));
750 
751  // The residual-squared that the line-search will reduce
752  // Later it will get contributions from epp and ic, but
753  // at present these are zero
754  Real nr_res2 = 0;
755  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
756  if (act[surface])
757  nr_res2 += 0.5 * Utility::pow<2>(f[surface] / _f[modelNumber(surface)]->_f_tol);
758 
759  successful_return = false;
760 
761  bool still_finding_solution = true;
762  while (still_finding_solution)
763  {
764  single_step_success = true;
765  unsigned int local_iter = 0;
766 
767  // The Newton-Raphson loops
768  while (nr_res2 > 0.5 && local_iter++ < _max_iter && single_step_success)
769  single_step_success = singleStep(nr_res2,
770  stress,
771  intnl_old,
772  intnl,
773  pm,
774  delta_dp,
775  E_inv,
776  f,
777  epp,
778  ic,
779  act,
780  deact_scheme,
781  linesearch_needed,
782  ld_encountered);
783 
784  bool nr_good = (nr_res2 <= 0.5 && local_iter <= _max_iter && single_step_success);
785 
786  iter += local_iter;
787 
788  // 'act' might have changed due to using deact_scheme = optimized, so
789  actives_tried.insert(activeCombinationNumber(act));
790 
791  if (!nr_good)
792  {
793  // failure of NR routine.
794  // We might be able to change the deactivation_scheme and
795  // then re-try the NR starting from the initial values
796  // Or, if deact_scheme == "dumb", we just increarse the
797  // dumb_iteration number and re-try
798  bool change_scheme = false;
799  bool increment_dumb = false;
800  change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
801  if (!change_scheme && deact_scheme == dumb)
802  increment_dumb = canIncrementDumb(dumb_iteration);
803 
804  still_finding_solution = (change_scheme || increment_dumb);
805 
806  if (change_scheme)
807  changeScheme(initial_act,
808  can_revert_to_dumb,
809  initial_stress,
810  intnl_old,
811  deact_scheme,
812  act,
813  dumb_iteration,
814  dumb_order);
815 
816  if (increment_dumb)
817  incrementDumb(dumb_iteration, dumb_order, act);
818 
819  if (!still_finding_solution)
820  {
821  // we cannot change the scheme, or have run out of "dumb" options
822  successful_return = false;
823  break;
824  }
825  }
826 
827  bool kt_good = false;
828  if (nr_good)
829  {
830  // check Kuhn-Tucker
831  kt_good = checkKuhnTucker(f, pm, act);
832  if (!kt_good)
833  {
834  if (deact_scheme != dumb)
835  {
836  applyKuhnTucker(f, pm, act);
837 
838  // true if we haven't tried this active set before
839  still_finding_solution =
840  (actives_tried.find(activeCombinationNumber(act)) == actives_tried.end());
841  if (!still_finding_solution)
842  {
843  // must have tried turning off the constraints already.
844  // so try changing the scheme
845  if (canChangeScheme(deact_scheme, can_revert_to_dumb))
846  {
847  still_finding_solution = true;
848  changeScheme(initial_act,
849  can_revert_to_dumb,
850  initial_stress,
851  intnl_old,
852  deact_scheme,
853  act,
854  dumb_iteration,
855  dumb_order);
856  }
857  }
858  }
859  else
860  {
861  bool increment_dumb = false;
862  increment_dumb = canIncrementDumb(dumb_iteration);
863  still_finding_solution = increment_dumb;
864  if (increment_dumb)
865  incrementDumb(dumb_iteration, dumb_order, act);
866  }
867 
868  if (!still_finding_solution)
869  {
870  // have tried turning off the constraints already,
871  // or have run out of "dumb" options
872  successful_return = false;
873  break;
874  }
875  }
876  }
877 
878  bool admissible = false;
879  if (nr_good && kt_good)
880  {
881  // check admissible
882  std::vector<Real> all_f;
883  if (_num_surfaces == 1)
884  admissible = true; // for a single surface if NR has exited successfully then (stress,
885  // intnl) must be admissible
886  else
887  admissible = checkAdmissible(stress, intnl, all_f);
888 
889  if (!admissible)
890  {
891  // Not admissible.
892  // We can try adding constraints back in
893  // We can try changing the deactivation scheme
894  // Or, if deact_scheme == dumb, just increase dumb_iteration
895  bool add_constraints = canAddConstraints(act, all_f);
896  if (add_constraints)
897  {
898  constraints_added = true;
899  std::vector<bool> act_plus(_num_surfaces,
900  false); // "act" with the positive constraints added in
901  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
902  if (act[surface] ||
903  (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol)))
904  act_plus[surface] = true;
905  if (actives_tried.find(activeCombinationNumber(act_plus)) == actives_tried.end())
906  {
907  // haven't tried this combination of actives yet
908  constraints_added = true;
909  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
910  act[surface] = act_plus[surface];
911  }
912  else
913  add_constraints = false; // haven't managed to add a new combination
914  }
915 
916  bool change_scheme = false;
917  bool increment_dumb = false;
918 
919  if (!add_constraints)
920  change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
921 
922  if (!add_constraints && !change_scheme && deact_scheme == dumb)
923  increment_dumb = canIncrementDumb(dumb_iteration);
924 
925  still_finding_solution = (add_constraints || change_scheme || increment_dumb);
926 
927  if (change_scheme)
928  changeScheme(initial_act,
929  can_revert_to_dumb,
930  initial_stress,
931  intnl_old,
932  deact_scheme,
933  act,
934  dumb_iteration,
935  dumb_order);
936 
937  if (increment_dumb)
938  incrementDumb(dumb_iteration, dumb_order, act);
939 
940  if (!still_finding_solution)
941  {
942  // we cannot change the scheme, or have run out of "dumb" options
943  successful_return = false;
944  break;
945  }
946  }
947  }
948 
949  successful_return = (nr_good && admissible && kt_good);
950  if (successful_return)
951  break;
952 
953  if (still_finding_solution)
954  {
955  stress = initial_stress;
956  delta_dp = RankTwoTensor(); // back to zero change in plastic strain
957  for (unsigned model = 0; model < _num_models; ++model)
958  intnl[model] = intnl_old[model]; // back to old internal params
959  pm.assign(_num_surfaces, 0.0); // back to zero plastic multipliers
960 
961  unsigned num_active = numberActive(act);
962  if (num_active == 0)
963  {
964  successful_return = false;
965  break; // failure
966  }
967 
968  actives_tried.insert(activeCombinationNumber(act));
969 
970  // Since "act" set has changed, either by changing deact_scheme, or by KT failing, so need to
971  // re-calculate nr_res2
972  yieldFunction(stress, intnl, act, f);
973 
974  nr_res2 = 0;
975  unsigned ind = 0;
976  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
977  if (act[surface])
978  {
979  if (f[ind] > _f[modelNumber(surface)]->_f_tol)
980  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
981  ind++;
982  }
983  }
984  }
985 
986  // returned, with either success or failure
987  if (successful_return)
988  {
989  plastic_strain += delta_dp;
990 
991  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
992  cumulative_pm[surface] += pm[surface];
993 
994  if (final_step)
995  consistent_tangent_operator =
996  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
997 
998  if (f.size() != _num_surfaces)
999  {
1000  // at this stage f.size() = num_active, but we need to return with all the yield functions
1001  // evaluated, so:
1002  act.assign(_num_surfaces, true);
1003  yieldFunction(stress, intnl, act, f);
1004  }
1005  }
1006 
1007  return successful_return;
1008 }
bool canAddConstraints(const std::vector< bool > &act, const std::vector< Real > &all_f)
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
virtual bool singleStep(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, RankTwoTensor &delta_dp, const RankFourTensor &E_inv, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, std::vector< bool > &active, DeactivationSchemeEnum deactivation_scheme, bool &linesearch_needed, bool &ld_encountered)
Performs a single Newton-Raphson + linesearch step Constraints are deactivated and the step is re-don...
virtual void buildActiveConstraints(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
Constructs a set of active constraints, given the yield functions, f.
bool canChangeScheme(DeactivationSchemeEnum current_deactivation_scheme, bool can_revert_to_dumb)
unsigned int _max_iter
Maximum number of Newton-Raphson iterations allowed.
unsigned int activeCombinationNumber(const std::vector< bool > &act)
virtual void applyKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
void changeScheme(const std::vector< bool > &initial_act, bool can_revert_to_dumb, const RankTwoTensor &initial_stress, const std::vector< Real > &intnl_old, DeactivationSchemeEnum &current_deactivation_scheme, std::vector< bool > &act, int &dumb_iteration, std::vector< unsigned int > &dumb_order)
virtual bool checkKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, const std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
virtual void incrementDumb(int &dumb_iteration, const std::vector< unsigned int > &dumb_order, std::vector< bool > &act)
Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of...
virtual bool checkAdmissible(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &all_f)
Checks whether the yield functions are in the admissible region.
unsigned int _num_models
Number of plastic models for this material.
RankFourTensor consistentTangentOperator(const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &pm_this_step, const std::vector< Real > &cumulative_pm)
Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rat...
virtual unsigned int numberActive(const std::vector< bool > &active)
counts the number of active constraints
virtual bool quickStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, RankFourTensor &consistent_tangent_operator, const quickStep_called_from_t called_from, bool final_step)
Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined ...
void buildDumbOrder(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
Builds the order which "dumb" activation will take.
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme

◆ returnMapAll()

bool MultiPlasticityRawComponentAssembler::returnMapAll ( const RankTwoTensor &  trial_stress,
const std::vector< Real > &  intnl_old,
const RankFourTensor &  E_ijkl,
Real  ep_plastic_tolerance,
RankTwoTensor &  stress,
std::vector< Real > &  intnl,
std::vector< Real > &  pm,
std::vector< Real > &  cumulative_pm,
RankTwoTensor &  delta_dp,
std::vector< Real > &  yf,
unsigned &  num_successful_plastic_returns,
unsigned &  custom_model 
)
protectedinherited

Performs a returnMap for each plastic model using their inbuilt returnMap functions.

Performs a returnMap for each plastic model.

This may be used to quickly ascertain whether a (trial_stress, intnl_old) configuration is admissible, or whether a single model's customized returnMap function can provide a solution to the return-map problem, or whether a full Newton-Raphson approach such as implemented in ComputeMultiPlasticityStress is needed.

There are three cases mentioned below: (A) The (trial_stress, intnl_old) configuration is admissible according to all plastic models (B) The (trial_stress, intnl_old) configuration is inadmissible to exactly one plastic model, and that model can successfully use its customized returnMap function to provide a returned (stress, intnl) configuration, and that configuration is admissible according to all plastic models (C) All other cases. This includes customized returnMap functions failing, or more than one plastic_model being inadmissible, etc

Parameters
trial_stressthe trial stress
intnl_oldthe old values of the internal parameters
E_ijklthe elasticity tensor
ep_plastic_tolerancethe tolerance on the plastic strain
[out]stressis set to trial_stress in case (A) or (C), and the returned value of stress in case (B).
[out]intnlis set to intnl_old in case (A) or (C), and the returned value of intnl in case (B)
[out]pmZero in case (A) or (C), otherwise the plastic multipliers needed to bring about the returnMap in case (B)
[in/out]cumulative_pm cumulative plastic multipliers, updated in case (B), otherwise left untouched
[out]delta_dpis unchanged in case (A) or (C), and is set to the change in plastic strain in case(B)
[out]yfwill contain the yield function values at (stress, intnl)
[out]num_successful_plastic_returnswill be 0 for (A) and (C), and 1 for (B)
Returns
true in case (A) and (B), and false in case (C)

If all models actually signal "elastic" by returning true from their returnMap, and by returning model_plastically_active=0, then yf will contain the yield function values num_successful_plastic_returns will be zero intnl = intnl_old delta_dp will be unchanged from its input value stress will be set to trial_stress pm will be zero cumulative_pm will be unchanged return value will be true num_successful_plastic_returns = 0

If only one model signals "plastically active" by returning true from its returnMap, and by returning model_plastically_active=1, then yf will contain the yield function values num_successful_plastic_returns will be one intnl will be set by the returnMap algorithm delta_dp will be set by the returnMap algorithm stress will be set by the returnMap algorithm pm will be nonzero for the single model, and zero for other models cumulative_pm will be updated return value will be true num_successful_plastic_returns = 1

If >1 model signals "plastically active" or if >=1 model's returnMap fails, then yf will contain the yield function values num_successful_plastic_returns will be set appropriately intnl = intnl_old delta_dp will be unchanged from its input value stress will be set to trial_stress pm will be zero cumulative_pm will be unchanged return value will be true if all returnMap functions returned true, otherwise it will be false num_successful_plastic_returns is set appropriately

Definition at line 598 of file MultiPlasticityRawComponentAssembler.C.

Referenced by quickStep().

610 {
611  mooseAssert(intnl_old.size() == _num_models,
612  "returnMapAll: Incorrect size of internal parameters");
613  mooseAssert(intnl.size() == _num_models, "returnMapAll: Incorrect size of internal parameters");
614  mooseAssert(pm.size() == _num_surfaces, "returnMapAll: Incorrect size of pm");
615 
616  num_successful_plastic_returns = 0;
617  yf.resize(0);
618  pm.assign(_num_surfaces, 0.0);
619 
620  RankTwoTensor returned_stress; // each model will give a returned_stress. if only one model is
621  // plastically active, i set stress=returned_stress, so as to
622  // record this returned value
623  std::vector<Real> model_f;
624  RankTwoTensor model_delta_dp;
625  std::vector<Real> model_pm;
626  bool trial_stress_inadmissible;
627  bool successful_return = true;
628  unsigned the_single_plastic_model = 0;
629  bool using_custom_return_map = true;
630 
631  // run through all the plastic models, performing their
632  // returnMap algorithms.
633  // If one finds (trial_stress, intnl) inadmissible and
634  // successfully returns, break from the loop to evaluate
635  // all the others at that returned stress
636  for (unsigned model = 0; model < _num_models; ++model)
637  {
638  if (using_custom_return_map)
639  {
640  model_pm.assign(_f[model]->numberSurfaces(), 0.0);
641  bool model_returned = _f[model]->returnMap(trial_stress,
642  intnl_old[model],
643  E_ijkl,
644  ep_plastic_tolerance,
645  returned_stress,
646  intnl[model],
647  model_pm,
648  model_delta_dp,
649  model_f,
650  trial_stress_inadmissible);
651  if (!trial_stress_inadmissible)
652  {
653  // in the elastic zone: record the yield-function values (returned_stress, intnl, model_pm
654  // and model_delta_dp are undefined)
655  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
656  ++model_surface)
657  yf.push_back(model_f[model_surface]);
658  }
659  else if (trial_stress_inadmissible && !model_returned)
660  {
661  // in the plastic zone, and the customized returnMap failed
662  // for some reason (or wasn't implemented). The coder
663  // should have correctly returned model_f(trial_stress, intnl_old)
664  // so record them
665  // (returned_stress, intnl, model_pm and model_delta_dp are undefined)
666  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
667  ++model_surface)
668  yf.push_back(model_f[model_surface]);
669  // now there's almost zero point in using the custom
670  // returnMap functions
671  using_custom_return_map = false;
672  successful_return = false;
673  }
674  else
675  {
676  // in the plastic zone, and the customized returnMap
677  // succeeded.
678  // record the first returned_stress and delta_dp if everything is going OK
679  // as they could be the actual answer
680  if (trial_stress_inadmissible)
681  num_successful_plastic_returns++;
682  the_single_plastic_model = model;
683  stress = returned_stress;
684  // note that i break here, and don't push_back
685  // model_f to yf. So now yf contains only the values of
686  // model_f from previous models to the_single_plastic_model
687  // also i don't set delta_dp = model_delta_dp yet, because
688  // i might find problems later on
689  // also, don't increment cumulative_pm for the same reason
690 
691  break;
692  }
693  }
694  else
695  {
696  // not using custom returnMap functions because one
697  // has already failed and that one said trial_stress
698  // was inadmissible. So now calculate the yield functions
699  // at the trial stress
700  _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
701  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
702  yf.push_back(model_f[model_surface]);
703  }
704  }
705 
706  if (num_successful_plastic_returns == 0)
707  {
708  // here either all the models were elastic (successful_return=true),
709  // or some were plastic and either the customized returnMap failed
710  // or wasn't implemented (successful_return=false).
711  // In either case, have to set the following:
712  stress = trial_stress;
713  for (unsigned model = 0; model < _num_models; ++model)
714  intnl[model] = intnl_old[model];
715  return successful_return;
716  }
717 
718  // Now we know that num_successful_plastic_returns == 1 and all the other
719  // models (with model number < the_single_plastic_model) must have been
720  // admissible at (trial_stress, intnl). However, all models might
721  // not be admissible at (trial_stress, intnl), so must check that
722  std::vector<Real> yf_at_returned_stress(0);
723  bool all_admissible = true;
724  for (unsigned model = 0; model < _num_models; ++model)
725  {
726  if (model == the_single_plastic_model)
727  {
728  // no need to spend time calculating the yield function: we know its admissible
729  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
730  yf_at_returned_stress.push_back(model_f[model_surface]);
731  continue;
732  }
733  _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
734  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
735  {
736  if (model_f[model_surface] > _f[model]->_f_tol)
737  // bummer, this model is not admissible at the returned_stress
738  all_admissible = false;
739  yf_at_returned_stress.push_back(model_f[model_surface]);
740  }
741  if (!all_admissible)
742  // no point in continuing computing yield functions
743  break;
744  }
745 
746  if (!all_admissible)
747  {
748  // we tried using the returned value of stress predicted by
749  // the_single_plastic_model, but it wasn't admissible according
750  // to other plastic models. We need to set:
751  stress = trial_stress;
752  for (unsigned model = 0; model < _num_models; ++model)
753  intnl[model] = intnl_old[model];
754  // and calculate the remainder of the yield functions at trial_stress
755  for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
756  {
757  _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
758  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
759  yf.push_back(model_f[model_surface]);
760  }
761  num_successful_plastic_returns = 0;
762  return false;
763  }
764 
765  // So the customized returnMap algorithm can provide a returned
766  // (stress, intnl) configuration, and that is admissible according
767  // to all plastic models
768  yf.resize(0);
769  for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
770  yf.push_back(yf_at_returned_stress[surface]);
771  delta_dp = model_delta_dp;
772  for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
773  ++model_surface)
774  {
775  cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
776  model_pm[model_surface];
777  pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
778  }
779  custom_model = the_single_plastic_model;
780  return true;
781 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ rot()

RankTwoTensor ComputeMultiPlasticityStress::rot ( const RankTwoTensor &  tens)
private

Definition at line 308 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

309 {
310  if (!_n_supplied)
311  return tens;
312  return tens.rotated(_rot);
313 }
bool _n_supplied
User supplied the transverse direction vector.
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)

◆ singleStep()

bool ComputeMultiPlasticityStress::singleStep ( Real &  nr_res2,
RankTwoTensor &  stress,
const std::vector< Real > &  intnl_old,
std::vector< Real > &  intnl,
std::vector< Real > &  pm,
RankTwoTensor &  delta_dp,
const RankFourTensor &  E_inv,
std::vector< Real > &  f,
RankTwoTensor &  epp,
std::vector< Real > &  ic,
std::vector< bool > &  active,
DeactivationSchemeEnum  deactivation_scheme,
bool &  linesearch_needed,
bool &  ld_encountered 
)
protectedvirtual

Performs a single Newton-Raphson + linesearch step Constraints are deactivated and the step is re-done if deactivation_scheme is set appropriately.

Parameters
[in,out]nr_res2Residual-squared that the line-search will reduce
[in,out]stressstress
[in]intnl_oldold values of the internal parameters
[in,out]intnlinternal parameters
[in,out]pmplastic multipliers
[in,out]delta_dpChange in plastic strain from start of "time" step to current configuration (plastic_strain - plastic_strain_old)
[in]E_invInverse of the elasticity tensor
[in,out]fYield function(s). Upon successful exit only the active constraints are contained in f
[in,out]eppPlastic strain increment constraint
[in,out]icInternal constraint. Upon successful exit only the active constraints are contained in ic
activeThe active constraints. This is may be modified, depending upon deactivation_scheme
deactivation_schemeThe scheme used for deactivating constraints
[out]linesearch_neededTrue if a linesearch was employed during this Newton-Raphson step
[out]ld_encounteredTrue if a linear-dependence of the flow directions was encountered at any stage during the Newton-Raphson proceedure
Returns
true if the step was successful, ie, if the linesearch was successful and the number of constraints wasn't reduced to zero via deactivation

Definition at line 1078 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1092 {
1093  bool successful_step; // return value
1094 
1095  Real nr_res2_before_step = nr_res2;
1096  RankTwoTensor stress_before_step;
1097  std::vector<Real> intnl_before_step;
1098  std::vector<Real> pm_before_step;
1099  RankTwoTensor delta_dp_before_step;
1100 
1101  if (deactivation_scheme == optimized)
1102  {
1103  // we potentially use the "before_step" quantities, so record them here
1104  stress_before_step = stress;
1105  intnl_before_step.resize(_num_models);
1106  for (unsigned model = 0; model < _num_models; ++model)
1107  intnl_before_step[model] = intnl[model];
1108  pm_before_step.resize(_num_surfaces);
1109  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1110  pm_before_step[surface] = pm[surface];
1111  delta_dp_before_step = delta_dp;
1112  }
1113 
1114  // During the Newton-Raphson procedure, we'll be
1115  // changing the following parameters in order to
1116  // (attempt to) satisfy the constraints.
1117  RankTwoTensor dstress; // change in stress
1118  std::vector<Real> dpm; // change in plasticity multipliers ("consistency parameters"). For ALL
1119  // contraints (active and deactive)
1120  std::vector<Real>
1121  dintnl; // change in internal parameters. For ALL internal params (active and deactive)
1122 
1123  // The constraints that have been deactivated for this NR step
1124  // due to the flow directions being linearly dependent
1125  std::vector<bool> deact_ld;
1126  deact_ld.assign(_num_surfaces, false);
1127 
1128  /* After NR and linesearch, if _deactivation_scheme == "optimized", the
1129  * active plasticity multipliers are checked for non-negativity. If some
1130  * are negative then they are deactivated forever, and the NR step is
1131  * re-done starting from the _before_step quantities recorded above
1132  */
1133  bool constraints_changing = true;
1134  bool reinstated_actives;
1135  while (constraints_changing)
1136  {
1137  // calculate dstress, dpm and dintnl for one full Newton-Raphson step
1138  nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
1139 
1140  for (unsigned surface = 0; surface < deact_ld.size(); ++surface)
1141  if (deact_ld[surface])
1142  {
1143  ld_encountered = true;
1144  break;
1145  }
1146 
1147  // perform a line search
1148  // The line-search will exit with updated values
1149  successful_step = lineSearch(nr_res2,
1150  stress,
1151  intnl_old,
1152  intnl,
1153  pm,
1154  E_inv,
1155  delta_dp,
1156  dstress,
1157  dpm,
1158  dintnl,
1159  f,
1160  epp,
1161  ic,
1162  active,
1163  deact_ld,
1164  linesearch_needed);
1165 
1166  if (!successful_step)
1167  // completely bomb out
1168  return successful_step;
1169 
1170  // See if any active constraints need to be removed, and the step re-done
1171  constraints_changing = false;
1172  if (deactivation_scheme == optimized)
1173  {
1174  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1175  if (active[surface] && pm[surface] < 0.0)
1176  constraints_changing = true;
1177  }
1178 
1179  if (constraints_changing)
1180  {
1181  stress = stress_before_step;
1182  delta_dp = delta_dp_before_step;
1183  nr_res2 = nr_res2_before_step;
1184  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1185  {
1186  if (active[surface] && pm[surface] < 0.0)
1187  {
1188  // turn off the constraint forever
1189  active[surface] = false;
1190  pm_before_step[surface] = 0.0;
1191  intnl_before_step[modelNumber(surface)] =
1192  intnl_old[modelNumber(surface)]; // don't want to muck-up hardening!
1193  }
1194  intnl[modelNumber(surface)] = intnl_before_step[modelNumber(surface)];
1195  pm[surface] = pm_before_step[surface];
1196  }
1197  if (numberActive(active) == 0)
1198  {
1199  // completely bomb out
1200  successful_step = false;
1201  return successful_step;
1202  }
1203  }
1204 
1205  // reinstate any active values that have been turned off due to linear-dependence
1206  reinstated_actives = reinstateLinearDependentConstraints(deact_ld);
1207  } // ends "constraints_changing" loop
1208 
1209  // if active constraints were reinstated then nr_res2 needs to be re-calculated so it is correct
1210  // upson returning
1211  if (reinstated_actives)
1212  {
1213  bool completely_converged = true;
1214  if (successful_step && nr_res2 < 0.5)
1215  {
1216  // Here we have converged to the correct solution if
1217  // all the yield functions are < 0. Excellent!
1218  //
1219  // This is quite tricky - perhaps i can refactor to make it more obvious.
1220  //
1221  // Because actives are now reinstated, the residual2
1222  // calculation below will give nr_res2 > 0.5, because it won't
1223  // realise that we only had to set the active-but-not-deactivated f = 0,
1224  // and not the entire active set. If we pass that nr_res2 back from
1225  // this function then the calling function will not realise we've converged!
1226  // Therefore, check for this case
1227  unsigned ind = 0;
1228  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1229  if (active[surface])
1230  if (f[ind++] > _f[modelNumber(surface)]->_f_tol)
1231  completely_converged = false;
1232  }
1233  else
1234  completely_converged = false;
1235 
1236  if (!completely_converged)
1237  nr_res2 = residual2(pm, f, epp, ic, active, deact_ld);
1238  }
1239 
1240  return successful_step;
1241 }
virtual bool reinstateLinearDependentConstraints(std::vector< bool > &deactivated_due_to_ld)
makes all deactivated_due_to_ld false, and if >0 of them were initially true, returns true ...
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
virtual bool lineSearch(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, const RankFourTensor &E_inv, RankTwoTensor &delta_dp, const RankTwoTensor &dstress, const std::vector< Real > &dpm, const std::vector< Real > &dintnl, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, bool &linesearch_needed)
Performs a line search.
virtual Real residual2(const std::vector< Real > &pm, const std::vector< Real > &f, const RankTwoTensor &epp, const std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld)
The residual-squared.
unsigned int _num_models
Number of plastic models for this material.
virtual unsigned int numberActive(const std::vector< bool > &active)
counts the number of active constraints
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.

◆ yieldFunction()

void MultiPlasticityRawComponentAssembler::yieldFunction ( const RankTwoTensor &  stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  f 
)
protectedvirtualinherited

The active yield function(s)

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active yield functions are put into "f"
[out]fthe yield function (or functions in the case of multisurface plasticity)

Definition at line 97 of file MultiPlasticityRawComponentAssembler.C.

Referenced by buildDumbOrder(), MultiPlasticityLinearSystem::calculateConstraints(), checkAdmissible(), MultiPlasticityDebugger::fddyieldFunction_dintnl(), MultiPlasticityDebugger::fddyieldFunction_dstress(), and returnMap().

101 {
102  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
103  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
104 
105  f.resize(0);
106  std::vector<unsigned int> active_surfaces_of_model;
107  std::vector<unsigned int>::iterator active_surface;
108  std::vector<Real> model_f;
109  for (unsigned model = 0; model < _num_models; ++model)
110  {
111  activeModelSurfaces(model, active, active_surfaces_of_model);
112  if (active_surfaces_of_model.size() > 0)
113  {
114  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
115  for (active_surface = active_surfaces_of_model.begin();
116  active_surface != active_surfaces_of_model.end();
117  ++active_surface)
118  f.push_back(model_f[*active_surface]);
119  }
120  }
121 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

Member Data Documentation

◆ _base_name

const std::string ComputeStressBase::_base_name
protectedinherited

◆ _constraints_added

MaterialProperty<Real>& ComputeMultiPlasticityStress::_constraints_added
protected

Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)

Definition at line 126 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _cosserat

bool ComputeMultiPlasticityStress::_cosserat
protected

whether Cosserat mechanics should be used

Definition at line 150 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), initQpStatefulProperties(), postReturnMap(), and preReturnMap().

◆ _couple_stress

MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_couple_stress
protected

the Cosserat couple-stress

Definition at line 159 of file ComputeMultiPlasticityStress.h.

◆ _couple_stress_old

const MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_couple_stress_old
protected

the old value of Cosserat couple-stress

Definition at line 162 of file ComputeMultiPlasticityStress.h.

◆ _cumulative_pm

std::vector<Real> ComputeMultiPlasticityStress::_cumulative_pm
protected

the sum of the plastic multipliers over all the sub-steps.

This is used for calculating the consistent tangent operator

Definition at line 70 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

◆ _curvature

const MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_curvature
protected

The Cosserat curvature strain.

Definition at line 153 of file ComputeMultiPlasticityStress.h.

◆ _deactivation_scheme

enum ComputeMultiPlasticityStress::DeactivationSchemeEnum ComputeMultiPlasticityStress::_deactivation_scheme
protected

◆ _dummy_pm

std::vector<Real> ComputeMultiPlasticityStress::_dummy_pm
protected

dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStress

Definition at line 64 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _elastic_flexural_rigidity_tensor

const MaterialProperty<RankFourTensor>* ComputeMultiPlasticityStress::_elastic_flexural_rigidity_tensor
protected

The Cosserat elastic flexural rigidity tensor.

Definition at line 156 of file ComputeMultiPlasticityStress.h.

◆ _elastic_strain

MaterialProperty<RankTwoTensor>& ComputeStressBase::_elastic_strain
protectedinherited

◆ _elastic_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_elastic_strain_old
protected

Old value of elastic strain.

Definition at line 147 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _elasticity_tensor

const MaterialProperty<RankFourTensor>& ComputeStressBase::_elasticity_tensor
protectedinherited

◆ _elasticity_tensor_name

const std::string ComputeStressBase::_elasticity_tensor_name
protectedinherited

◆ _epp_tol

Real ComputeMultiPlasticityStress::_epp_tol
protected

Tolerance on the plastic strain increment ("direction") constraint.

Definition at line 61 of file ComputeMultiPlasticityStress.h.

Referenced by ComputeMultiPlasticityStress(), quickStep(), and residual2().

◆ _extra_stress

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_extra_stress
protectedinherited

Extra stress tensor.

Definition at line 47 of file ComputeStressBase.h.

Referenced by ComputeStressBase::computeQpProperties().

◆ _f

std::vector<const TensorMechanicsPlasticModel *> MultiPlasticityRawComponentAssembler::_f
protectedinherited

◆ _fspb_debug

MooseEnum MultiPlasticityDebugger::_fspb_debug
protectedinherited

none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobian_and_linear_system - check entire jacobian and check that Ax=b

Definition at line 61 of file MultiPlasticityDebugger.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _fspb_debug_intnl

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl
protectedinherited

◆ _fspb_debug_intnl_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change
protectedinherited

◆ _fspb_debug_pm

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm
protectedinherited

◆ _fspb_debug_pm_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change
protectedinherited

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 76 of file MultiPlasticityDebugger.h.

Referenced by MultiPlasticityDebugger::fdJacobian(), and MultiPlasticityDebugger::outputAndCheckDebugParameters().

◆ _fspb_debug_stress

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress
protectedinherited

◆ _fspb_debug_stress_change

Real MultiPlasticityDebugger::_fspb_debug_stress_change
protectedinherited

◆ _ignore_failures

bool ComputeMultiPlasticityStress::_ignore_failures
protected

Even if the returnMap fails, return the best values found for stress and internal parameters.

Definition at line 50 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

◆ _initial_stress_fcn

std::vector<Function *> ComputeStressBase::_initial_stress_fcn
protectedinherited

initial stress components

Definition at line 50 of file ComputeStressBase.h.

◆ _intnl

MaterialProperty<std::vector<Real> >& ComputeMultiPlasticityStress::_intnl
protected

internal parameters

Definition at line 108 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _intnl_old

const MaterialProperty<std::vector<Real> >& ComputeMultiPlasticityStress::_intnl_old
protected

old values of internal parameters

Definition at line 111 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

◆ _iter

MaterialProperty<Real>& ComputeMultiPlasticityStress::_iter
protected

Number of Newton-Raphson iterations used in the return-map.

Definition at line 117 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _Jacobian_mult

MaterialProperty<RankFourTensor>& ComputeStressBase::_Jacobian_mult
protectedinherited

◆ _Jacobian_mult_couple

MaterialProperty<RankFourTensor>* ComputeMultiPlasticityStress::_Jacobian_mult_couple
protected

derivative of couple-stress w.r.t. curvature

Definition at line 165 of file ComputeMultiPlasticityStress.h.

◆ _ld_encountered

MaterialProperty<Real>& ComputeMultiPlasticityStress::_ld_encountered
protected

Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true, 0 otherwise)

Definition at line 123 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _linesearch_needed

MaterialProperty<Real>& ComputeMultiPlasticityStress::_linesearch_needed
protected

Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise)

Definition at line 120 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _max_iter

unsigned int ComputeMultiPlasticityStress::_max_iter
protected

Maximum number of Newton-Raphson iterations allowed.

Definition at line 41 of file ComputeMultiPlasticityStress.h.

Referenced by returnMap().

◆ _max_stepsize_for_dumb

Real ComputeMultiPlasticityStress::_max_stepsize_for_dumb
protected

"dumb" deactivation will only be used if the stepsize falls below this quantity

Definition at line 47 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

◆ _mechanical_strain

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_mechanical_strain
protectedinherited

◆ _min_f_tol

Real MultiPlasticityLinearSystem::_min_f_tol
protectedinherited

Minimum value of the _f_tol parameters for the Yield Function User Objects.

Definition at line 135 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence(), and MultiPlasticityLinearSystem::MultiPlasticityLinearSystem().

◆ _min_stepsize

Real ComputeMultiPlasticityStress::_min_stepsize
protected

Minimum fraction of applied strain that may be applied during adaptive stepsizing.

Definition at line 44 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

◆ _my_curvature

RankTwoTensor ComputeMultiPlasticityStress::_my_curvature
protected

Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)

Definition at line 177 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), postReturnMap(), and preReturnMap().

◆ _my_elasticity_tensor

RankFourTensor ComputeMultiPlasticityStress::_my_elasticity_tensor
protected

Elasticity tensor that can be rotated by this class (ie, its not const)

Definition at line 168 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), plasticStep(), postReturnMap(), and preReturnMap().

◆ _my_flexural_rigidity_tensor

RankFourTensor ComputeMultiPlasticityStress::_my_flexural_rigidity_tensor
protected

Flexual rigidity tensor that can be rotated by this class (ie, its not const)

Definition at line 174 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), postReturnMap(), and preReturnMap().

◆ _my_strain_increment

RankTwoTensor ComputeMultiPlasticityStress::_my_strain_increment
protected

Strain increment that can be rotated by this class, and split into multiple increments (ie, its not const)

Definition at line 171 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), postReturnMap(), and preReturnMap().

◆ _n

MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n
protected

current value of transverse direction

Definition at line 129 of file ComputeMultiPlasticityStress.h.

Referenced by initQpStatefulProperties(), postReturnMap(), and preReturnMap().

◆ _n_input

RealVectorValue ComputeMultiPlasticityStress::_n_input
protected

the supplied transverse direction vector

Definition at line 93 of file ComputeMultiPlasticityStress.h.

Referenced by ComputeMultiPlasticityStress(), and initQpStatefulProperties().

◆ _n_old

const MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n_old
protected

old value of transverse direction

Definition at line 132 of file ComputeMultiPlasticityStress.h.

Referenced by postReturnMap().

◆ _n_supplied

bool ComputeMultiPlasticityStress::_n_supplied
protected

User supplied the transverse direction vector.

Definition at line 90 of file ComputeMultiPlasticityStress.h.

Referenced by ComputeMultiPlasticityStress(), postReturnMap(), preReturnMap(), and rot().

◆ _num_models

unsigned int MultiPlasticityRawComponentAssembler::_num_models
protectedinherited

◆ _num_surfaces

unsigned int MultiPlasticityRawComponentAssembler::_num_surfaces
protectedinherited

Number of surfaces within the plastic models.

For many situations this will be = _num_models since each model will contain just one surface. More generally it is >= _num_models. For instance, Mohr-Coulomb is a single model with 6 surfaces

Definition at line 64 of file MultiPlasticityRawComponentAssembler.h.

Referenced by activeCombinationNumber(), applyKuhnTucker(), MultiPlasticityRawComponentAssembler::buildActiveConstraints(), buildDumbOrder(), MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), canAddConstraints(), canIncrementDumb(), changeScheme(), checkAdmissible(), MultiPlasticityDebugger::checkDerivatives(), MultiPlasticityDebugger::checkJacobian(), checkKuhnTucker(), MultiPlasticityDebugger::checkSolution(), ComputeMultiPlasticityStress(), computeQpStress(), consistentTangentOperator(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityDebugger::dof_included(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityLinearSystem::eliminateLinearDependence(), MultiPlasticityDebugger::fddflowPotential_dintnl(), MultiPlasticityDebugger::fddflowPotential_dstress(), MultiPlasticityDebugger::fddyieldFunction_dintnl(), MultiPlasticityDebugger::fddyieldFunction_dstress(), MultiPlasticityDebugger::fdJacobian(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), incrementDumb(), initQpStatefulProperties(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), MultiPlasticityLinearSystem::nrStep(), numberActive(), MultiPlasticityDebugger::outputAndCheckDebugParameters(), plasticStep(), reinstateLinearDependentConstraints(), residual2(), returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

◆ _params

const InputParameters& MultiPlasticityRawComponentAssembler::_params
protectedinherited

◆ _perform_finite_strain_rotations

bool ComputeMultiPlasticityStress::_perform_finite_strain_rotations
protected

whether to perform the rotations necessary in finite-strain simulations

Definition at line 99 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _plastic_strain

MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain
protected

plastic strain

Definition at line 102 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), initQpStatefulProperties(), and postReturnMap().

◆ _plastic_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain_old
protected

Old value of plastic strain.

Definition at line 105 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _rot

RealTensorValue ComputeMultiPlasticityStress::_rot
protected

rotation matrix that takes _n to (0, 0, 1)

Definition at line 96 of file ComputeMultiPlasticityStress.h.

Referenced by postReturnMap(), preReturnMap(), and rot().

◆ _rotation_increment

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_rotation_increment
protected

Rotation increment (coming from ComputeIncrementalSmallStrain, for example)

Definition at line 141 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and postReturnMap().

◆ _specialIC

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC
protectedinherited

◆ _strain_increment

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_strain_increment
protected

strain increment (coming from ComputeIncrementalSmallStrain, for example)

Definition at line 135 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _stress

MaterialProperty<RankTwoTensor>& ComputeStressBase::_stress
protectedinherited

◆ _stress_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_stress_old
protected

Old value of stress.

Definition at line 144 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _surfaces_given_model

std::vector<std::vector<unsigned int> > MultiPlasticityRawComponentAssembler::_surfaces_given_model
protectedinherited

◆ _svd_tol

Real MultiPlasticityLinearSystem::_svd_tol
protectedinherited

Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent.

Definition at line 132 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

◆ _tangent_operator_type

enum ComputeMultiPlasticityStress::TangentOperatorEnum ComputeMultiPlasticityStress::_tangent_operator_type
protected

◆ _total_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_total_strain_old
protected

Old value of total strain (coming from ComputeIncrementalSmallStrain, for example)

Definition at line 138 of file ComputeMultiPlasticityStress.h.

◆ _yf

MaterialProperty<std::vector<Real> >& ComputeMultiPlasticityStress::_yf
protected

yield functions

Definition at line 114 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().


The documentation for this class was generated from the following files: