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") : nullptr),
157  _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
158  : nullptr),
159  _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : nullptr),
160  _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress")
161  : nullptr),
162  _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
163  : nullptr),
164 
169 {
170  if (_epp_tol <= 0)
171  mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
172 
173  if (_n_supplied)
174  {
175  // normalise the inputted transverse_direction
176  if (_n_input.norm() == 0)
177  mooseError(
178  "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
179  else
180  _n_input /= _n_input.norm();
181  }
182 
183  if (_num_surfaces == 1)
185 }
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 1571 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1572 {
1573  unsigned num = 0;
1574  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1575  if (act[surface])
1576  num += (1 << surface); // (1 << x) = 2^x
1577 
1578  return num;
1579 }
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 1310 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1313 {
1314  bool turned_off = false;
1315  unsigned ind = 0;
1316 
1317  // turn off all active surfaces that have f<0 and pm!=0
1318  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1319  {
1320  if (active[surface])
1321  {
1322  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1323  if (pm[surface] != 0)
1324  {
1325  turned_off = true;
1326  active[surface] = false;
1327  }
1328  }
1329  else if (pm[surface] != 0)
1330  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1331  "coding!!");
1332  }
1333 
1334  // if didn't turn off anything yet, turn off surface with minimum pm
1335  if (!turned_off)
1336  {
1337  int surface_to_turn_off = -1;
1338  Real min_pm = 0;
1339  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1340  if (pm[surface] < min_pm)
1341  {
1342  min_pm = pm[surface];
1343  surface_to_turn_off = surface;
1344  }
1345  if (surface_to_turn_off >= 0)
1346  active[surface_to_turn_off] = false;
1347  }
1348 }
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 1521 of file ComputeMultiPlasticityStress.C.

Referenced by changeScheme(), and returnMap().

1524 {
1525  if (dumb_order.size() != 0)
1526  return;
1527 
1528  std::vector<bool> act;
1529  act.assign(_num_surfaces, true);
1530 
1531  std::vector<Real> f;
1532  yieldFunction(stress, intnl, act, f);
1533  std::vector<RankTwoTensor> df_dstress;
1534  dyieldFunction_dstress(stress, intnl, act, df_dstress);
1535 
1536  typedef std::pair<Real, unsigned> pair_for_sorting;
1537  std::vector<pair_for_sorting> dist(_num_surfaces);
1538  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1539  {
1540  dist[surface].first = f[surface] / df_dstress[surface].L2norm();
1541  dist[surface].second = surface;
1542  }
1543  std::sort(dist.begin(), dist.end()); // sorted in ascending order of f/df_dstress
1544 
1545  dumb_order.resize(_num_surfaces);
1546  for (unsigned i = 0; i < _num_surfaces; ++i)
1547  dumb_order[i] = dist[_num_surfaces - 1 - i].second;
1548  // now dumb_order[0] is the surface with the greatest f/df_dstress
1549 }
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.
RankTwoTensorTempl< Real > RankTwoTensor
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.
RankTwoTensorTempl< Real > RankTwoTensor
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 1012 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1014 {
1015  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1016  if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1017  return true;
1018  return false;
1019 }
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 1022 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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

◆ canIncrementDumb()

bool ComputeMultiPlasticityStress::canIncrementDumb ( int  dumb_iteration)
protected

Definition at line 1564 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1565 {
1566  // (1 << _num_surfaces) = 2^_num_surfaces
1567  return ((dumb_iteration + 1) < (1 << _num_surfaces));
1568 }
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 1047 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1055 {
1056  if (current_deactivation_scheme == optimized &&
1059  {
1060  current_deactivation_scheme = safe;
1061  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1062  act[surface] = initial_act[surface];
1063  }
1064  else if ((current_deactivation_scheme == safe &&
1067  can_revert_to_dumb) ||
1068  (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1069  can_revert_to_dumb))
1070  {
1071  current_deactivation_scheme = dumb;
1072  dumb_iteration = 0;
1073  buildDumbOrder(initial_stress, intnl_old, dumb_order);
1074  incrementDumb(dumb_iteration, dumb_order, act);
1075  }
1076 }
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 1268 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1271 {
1272  std::vector<bool> act;
1273  act.assign(_num_surfaces, true);
1274 
1275  yieldFunction(stress, intnl, act, all_f);
1276 
1277  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1278  if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1279  return false;
1280 
1281  return true;
1282 }
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 1285 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1288 {
1289  unsigned ind = 0;
1290  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1291  {
1292  if (active[surface])
1293  {
1294  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1295  if (pm[surface] != 0)
1296  return false;
1297  }
1298  else if (pm[surface] != 0)
1299  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1300  "coding!!");
1301  }
1302  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1303  if (pm[surface] < 0)
1304  return false;
1305 
1306  return true;
1307 }
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 218 of file ComputeMultiPlasticityStress.C.

219 {
220  // the following "_my" variables can get rotated by preReturnMap and postReturnMap
223  if (_cosserat)
224  {
225  _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
226  _my_curvature = (*_curvature)[_qp];
227  }
228 
229  if (_fspb_debug == "jacobian_and_linear_system")
230  {
231  // cannot do this at initQpStatefulProperties level since E_ijkl is not defined
232  checkJacobian(_elasticity_tensor[_qp].invSymm(), _intnl_old[_qp]);
233  checkSolution(_elasticity_tensor[_qp].invSymm());
234  mooseError("Finite-differencing completed. Exiting with no error");
235  }
236 
237  preReturnMap(); // do rotations to new frame if necessary
238 
239  unsigned int number_iterations;
240  bool linesearch_needed = false;
241  bool ld_encountered = false;
242  bool constraints_added = false;
243 
244  _cumulative_pm.assign(_num_surfaces, 0);
245  // try a "quick" return first - this can be purely elastic, or a customised plastic return defined
246  // by a TensorMechanicsPlasticXXXX UserObject
247  const bool found_solution = quickStep(rot(_stress_old[_qp]),
248  _stress[_qp],
249  _intnl_old[_qp],
250  _intnl[_qp],
251  _dummy_pm,
253  rot(_plastic_strain_old[_qp]),
254  _plastic_strain[_qp],
257  _yf[_qp],
258  number_iterations,
259  _Jacobian_mult[_qp],
261  true);
262 
263  // if not purely elastic or the customised stuff failed, do some plastic return
264  if (!found_solution)
266  _stress[_qp],
267  _intnl_old[_qp],
268  _intnl[_qp],
269  rot(_plastic_strain_old[_qp]),
270  _plastic_strain[_qp],
273  _yf[_qp],
274  number_iterations,
275  linesearch_needed,
276  ld_encountered,
277  constraints_added,
278  _Jacobian_mult[_qp]);
279 
280  if (_cosserat)
281  {
282  (*_couple_stress)[_qp] = (*_elastic_flexural_rigidity_tensor)[_qp] * _my_curvature;
283  (*_Jacobian_mult_couple)[_qp] = _my_flexural_rigidity_tensor;
284  }
285 
286  postReturnMap(); // rotate back from new frame if necessary
287 
288  _iter[_qp] = 1.0 * number_iterations;
289  _linesearch_needed[_qp] = linesearch_needed;
290  _ld_encountered[_qp] = ld_encountered;
291  _constraints_added[_qp] = constraints_added;
292 
293  // Update measures of strain
295  (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
296 
297  // Rotate the tensors to the current configuration
299  {
300  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
301  _elastic_strain[_qp] =
302  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
303  _plastic_strain[_qp] =
304  _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
305  }
306 }
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 1582 of file ComputeMultiPlasticityStress.C.

Referenced by quickStep(), and returnMap().

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

Referenced by changeScheme(), and returnMap().

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

◆ initQpStatefulProperties()

void ComputeMultiPlasticityStress::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 188 of file ComputeMultiPlasticityStress.C.

189 {
191 
192  _plastic_strain[_qp].zero();
193 
194  _intnl[_qp].assign(_num_models, 0);
195 
196  _yf[_qp].assign(_num_surfaces, 0);
197 
198  _dummy_pm.assign(_num_surfaces, 0);
199 
200  _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
201  _linesearch_needed[_qp] = 0;
202  _ld_encountered[_qp] = 0;
203  _constraints_added[_qp] = 0;
204 
205  _n[_qp] = _n_input;
206 
207  if (_cosserat)
208  (*_couple_stress)[_qp].zero();
209 
210  if (_fspb_debug == "jacobian")
211  {
213  mooseError("Finite-differencing completed. Exiting with no error");
214  }
215 }
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 1388 of file ComputeMultiPlasticityStress.C.

Referenced by singleStep().

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

Referenced by returnMap(), and singleStep().

1259 {
1260  unsigned num_active = 0;
1261  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1262  if (active[surface])
1263  num_active++;
1264  return num_active;
1265 }
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 465 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

479 {
495  bool return_successful = false;
496 
497  // total number of Newton-Raphson iterations used
498  unsigned int iter = 0;
499 
500  Real step_size = 1.0;
501  Real time_simulated = 0.0;
502 
503  // the "good" variables hold the latest admissible stress
504  // and internal parameters.
505  RankTwoTensor stress_good = stress_old;
506  RankTwoTensor plastic_strain_good = plastic_strain_old;
507  std::vector<Real> intnl_good(_num_models);
508  for (unsigned model = 0; model < _num_models; ++model)
509  intnl_good[model] = intnl_old[model];
510  std::vector<Real> yf_good(_num_surfaces);
511 
512  // Following is necessary because I want strain_increment to be "const"
513  // but I also want to be able to subdivide an initial_stress
514  RankTwoTensor this_strain_increment = strain_increment;
515 
516  RankTwoTensor dep = step_size * this_strain_increment;
517 
518  _cumulative_pm.assign(_num_surfaces, 0);
519 
520  unsigned int num_consecutive_successes = 0;
521  while (time_simulated < 1.0 && step_size >= _min_stepsize)
522  {
523  iter = 0;
524  return_successful = returnMap(stress_good,
525  stress,
526  intnl_good,
527  intnl,
528  plastic_strain_good,
529  plastic_strain,
530  E_ijkl,
531  dep,
532  yf,
533  iter,
534  step_size <= _max_stepsize_for_dumb,
535  linesearch_needed,
536  ld_encountered,
537  constraints_added,
538  time_simulated + step_size >= 1,
539  consistent_tangent_operator,
541  iterations += iter;
542 
543  if (return_successful)
544  {
545  num_consecutive_successes += 1;
546  time_simulated += step_size;
547 
548  if (time_simulated < 1.0) // this condition is just for optimization: if time_simulated=1 then
549  // the "good" quantities are no longer needed
550  {
551  stress_good = stress;
552  plastic_strain_good = plastic_strain;
553  for (unsigned model = 0; model < _num_models; ++model)
554  intnl_good[model] = intnl[model];
555  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
556  yf_good[surface] = yf[surface];
557  if (num_consecutive_successes >= 2)
558  step_size *= 1.2;
559  }
560  step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
561  }
562  else
563  {
564  step_size *= 0.5;
565  num_consecutive_successes = 0;
566  stress = stress_good;
567  plastic_strain = plastic_strain_good;
568  for (unsigned model = 0; model < _num_models; ++model)
569  intnl[model] = intnl_good[model];
570  yf.resize(_num_surfaces); // might have excited with junk
571  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
572  yf[surface] = yf_good[surface];
573  dep = step_size * this_strain_increment;
574  }
575  }
576 
577  if (!return_successful)
578  {
579  if (_ignore_failures)
580  {
581  stress = stress_good;
582  plastic_strain = plastic_strain_good;
583  for (unsigned model = 0; model < _num_models; ++model)
584  intnl[model] = intnl_good[model];
585  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
586  yf[surface] = yf_good[surface];
587  }
588  else
589  {
590  Moose::out << "After reducing the stepsize to " << step_size
591  << " with original strain increment with L2norm " << this_strain_increment.L2norm()
592  << " the returnMap algorithm failed\n";
593 
594  _fspb_debug_stress = stress_good + E_ijkl * dep;
595  _fspb_debug_pm.assign(
597  1); // this is chosen arbitrarily - please change if a more suitable value occurs to you!
599  for (unsigned model = 0; model < _num_models; ++model)
600  _fspb_debug_intnl[model] = intnl_good[model];
604  mooseError("Exiting\n");
605  }
606  }
607 
608  return return_successful;
609 }
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 336 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

337 {
338  if (_n_supplied)
339  {
340  // this is a rotation matrix that will rotate "z" axis back to _n
341  _rot = _rot.transpose();
342 
343  // rotate the tensors back to original frame where _n is correctly oriented
344  _my_elasticity_tensor.rotate(_rot);
345  _Jacobian_mult[_qp].rotate(_rot);
346  _my_strain_increment.rotate(_rot);
347  _stress[_qp].rotate(_rot);
348  _plastic_strain[_qp].rotate(_rot);
349  if (_cosserat)
350  {
352  (*_Jacobian_mult_couple)[_qp].rotate(_rot);
353  _my_curvature.rotate(_rot);
354  (*_couple_stress)[_qp].rotate(_rot);
355  }
356 
357  // Rotate n by _rotation_increment
358  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
359  {
360  _n[_qp](i) = 0;
361  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
362  _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
363  }
364  }
365 }
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 317 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

318 {
319  if (_n_supplied)
320  {
321  // this is a rotation matrix that will rotate _n to the "z" axis
322  _rot = RotationMatrix::rotVecToZ(_n[_qp]);
323 
324  // rotate the tensors to this frame
325  _my_elasticity_tensor.rotate(_rot);
326  _my_strain_increment.rotate(_rot);
327  if (_cosserat)
328  {
330  _my_curvature.rotate(_rot);
331  }
332  }
333 }
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 368 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress(), and returnMap().

383 {
384  iterations = 0;
385 
386  unsigned num_plastic_returns;
387  RankTwoTensor delta_dp;
388 
389  // the following does the customized returnMap algorithm
390  // for all the plastic models.
391  unsigned custom_model = 0;
392  bool successful_return = returnMapAll(stress_old + E_ijkl * strain_increment,
393  intnl_old,
394  E_ijkl,
395  _epp_tol,
396  stress,
397  intnl,
398  pm,
399  cumulative_pm,
400  delta_dp,
401  yf,
402  num_plastic_returns,
403  custom_model);
404 
405  // the following updates the plastic_strain, when necessary
406  // and calculates the consistent_tangent_operator, when necessary
407  if (num_plastic_returns == 0)
408  {
409  // if successful_return = true, then a purely elastic step
410  // if successful_return = false, then >=1 plastic model is in
411  // inadmissible zone and failed to return using its customized
412  // returnMap function.
413  // In either case:
414  plastic_strain = plastic_strain_old;
415  if (successful_return && final_step)
416  {
417  if (called_from == computeQpStress_function)
418  consistent_tangent_operator = E_ijkl;
419  else // cannot necessarily use E_ijkl since different plastic models may have been active
420  // during other substeps
421  consistent_tangent_operator =
422  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
423  }
424  return successful_return;
425  }
426  else if (num_plastic_returns == 1 && successful_return)
427  {
428  // one model has successfully completed its custom returnMap algorithm
429  // and the other models have signalled they are elastic at
430  // the trial stress
431  plastic_strain = plastic_strain_old + delta_dp;
432  if (final_step)
433  {
434  if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
435  {
437  consistent_tangent_operator = E_ijkl;
438  else
439  {
440  std::vector<Real> custom_model_pm;
441  for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
442  custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
443  consistent_tangent_operator =
444  _f[custom_model]->consistentTangentOperator(stress_old + E_ijkl * strain_increment,
445  intnl_old[custom_model],
446  stress,
447  intnl[custom_model],
448  E_ijkl,
449  custom_model_pm);
450  }
451  }
452  else // cannot necessarily use the custom consistentTangentOperator since different plastic
453  // models may have been active during other substeps or the custom model says not to use
454  // its custom CTO algorithm
455  consistent_tangent_operator =
456  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
457  }
458  return true;
459  }
460  else // presumably returnMapAll is incorrectly coded!
461  mooseError("ComputeMultiPlasticityStress::quickStep should not get here!");
462 }
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 1245 of file ComputeMultiPlasticityStress.C.

Referenced by singleStep().

1247 {
1248  bool reinstated_actives = false;
1249  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1250  if (deactivated_due_to_ld[surface])
1251  reinstated_actives = true;
1252 
1253  deactivated_due_to_ld.assign(_num_surfaces, false);
1254  return reinstated_actives;
1255 }
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 1351 of file ComputeMultiPlasticityStress.C.

Referenced by lineSearch(), and singleStep().

1357 {
1358  Real nr_res2 = 0;
1359  unsigned ind = 0;
1360 
1361  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1362  if (active[surface])
1363  {
1364  if (!deactivated_due_to_ld[surface])
1365  {
1366  if (!(pm[surface] == 0 && f[ind] <= 0))
1367  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1368  }
1369  else if (deactivated_due_to_ld[surface] && f[ind] > 0)
1370  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1371  ind++;
1372  }
1373 
1374  nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
1375 
1376  std::vector<bool> active_not_deact(_num_surfaces);
1377  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1378  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1379  ind = 0;
1380  for (unsigned model = 0; model < _num_models; ++model)
1381  if (anyActiveSurfaces(model, active_not_deact))
1382  nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
1383 
1384  return nr_res2;
1385 }
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 612 of file ComputeMultiPlasticityStress.C.

Referenced by plasticStep().

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

Referenced by computeQpStress().

310 {
311  if (!_n_supplied)
312  return tens;
313  return tens.rotated(_rot);
314 }
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 1079 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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