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 ()
 Compute the stress and store it in the _stress material property for the current quadrature point. More...
 
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...
 
const std::string _elasticity_tensor_name
 Name of the elasticity tensor material property. More...
 
const MaterialProperty< RankFourTensor > & _elasticity_tensor
 Elasticity tensor material property. 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
 Base name prepended to all material property names to allow for multi-material systems. More...
 
const MaterialProperty< RankTwoTensor > & _mechanical_strain
 Mechanical strain material property. More...
 
MaterialProperty< RankTwoTensor > & _stress
 Stress material property. More...
 
MaterialProperty< RankTwoTensor > & _elastic_strain
 Elastic strain material property. More...
 
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 444 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  _elasticity_tensor_name(_base_name + "elasticity_tensor"),
123  _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
124  _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
125  _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
126  _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
127  _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
128  _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
129  _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
130  // for visualisation i convert it to Real
131  _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
132  // boolean, but for
133  // visualisation i
134  // convert it to Real
135  _ld_encountered(declareProperty<Real>(
136  "plastic_linear_dependence_encountered")), // this is really a boolean, but for
137  // visualisation i convert it to Real
138  _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
139  // boolean, but for
140  // visualisation i
141  // convert it to Real
142  _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
143  _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
144 
145  _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
146  _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
148  getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
149 
150  _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
151  _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
152 
153  // TODO: This design does NOT work. It makes these materials construction order dependent and it
154  // disregards block restrictions.
155  _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
156  hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
157  _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : nullptr),
159  _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
160  : nullptr),
161  _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : nullptr),
162  _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress")
163  : nullptr),
164  _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
165  : nullptr),
166 
171 {
172  if (_epp_tol <= 0)
173  mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
174 
175  if (_n_supplied)
176  {
177  // normalise the inputted transverse_direction
178  if (_n_input.norm() == 0)
179  mooseError(
180  "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
181  else
182  _n_input /= _n_input.norm();
183  }
184 
185  if (_num_surfaces == 1)
187 }
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)
const std::string _elasticity_tensor_name
Name of the elasticity tensor material property.
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
Base name prepended to all material property names to allow for multi-material systems.
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
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
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 1573 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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

Referenced by returnMap().

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

Referenced by changeScheme(), and returnMap().

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

Referenced by returnMap().

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

Referenced by returnMap().

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

◆ canIncrementDumb()

bool ComputeMultiPlasticityStress::canIncrementDumb ( int  dumb_iteration)
protected

Definition at line 1566 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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

Referenced by returnMap().

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

Referenced by returnMap().

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

Referenced by returnMap().

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

50 {
52 
53  // Add in extra stress
54  _stress[_qp] += _extra_stress[_qp];
55 }
virtual void computeQpStress()=0
Compute the stress and store it in the _stress material property for the current quadrature point...
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
const MaterialProperty< RankTwoTensor > & _extra_stress
Extra stress tensor.

◆ computeQpStress()

void ComputeMultiPlasticityStress::computeQpStress ( )
protectedvirtual

Compute the stress and store it in the _stress material property for the current quadrature point.

Implements ComputeStressBase.

Definition at line 220 of file ComputeMultiPlasticityStress.C.

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

Referenced by quickStep(), and returnMap().

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

Referenced by changeScheme(), and returnMap().

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

◆ initQpStatefulProperties()

void ComputeMultiPlasticityStress::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 190 of file ComputeMultiPlasticityStress.C.

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

Referenced by singleStep().

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

Referenced by returnMap(), and singleStep().

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

Referenced by computeQpStress().

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

Referenced by computeQpStress().

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

Referenced by computeQpStress().

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

Referenced by computeQpStress(), and returnMap().

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

Referenced by singleStep().

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

Referenced by lineSearch(), and singleStep().

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

Referenced by plasticStep().

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

Referenced by computeQpStress().

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

Referenced by returnMap().

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

Base name prepended to all material property names to allow for multi-material systems.

Definition at line 44 of file ComputeStressBase.h.

Referenced by ComputeLinearElasticStress::initialSetup(), and ComputeCosseratLinearElasticStress::initialSetup().

◆ _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 131 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _cosserat

bool ComputeMultiPlasticityStress::_cosserat
protected

whether Cosserat mechanics should be used

Definition at line 155 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 164 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 167 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 158 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 161 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 152 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _elasticity_tensor

const MaterialProperty<RankFourTensor>& ComputeMultiPlasticityStress::_elasticity_tensor
protected

Elasticity tensor material property.

Definition at line 104 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _elasticity_tensor_name

const std::string ComputeMultiPlasticityStress::_elasticity_tensor_name
protected

Name of the elasticity tensor material property.

Definition at line 102 of file ComputeMultiPlasticityStress.h.

◆ _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 54 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 57 of file ComputeStressBase.h.

◆ _intnl

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

internal parameters

Definition at line 113 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 116 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 122 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 170 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 128 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 125 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 182 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 173 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 179 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 176 of file ComputeMultiPlasticityStress.h.

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

◆ _n

MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n
protected

current value of transverse direction

Definition at line 134 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 137 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 107 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 110 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 146 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 140 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

◆ _stress

MaterialProperty<RankTwoTensor>& ComputeStressBase::_stress
protectedinherited

Stress material property.

Definition at line 49 of file ComputeStressBase.h.

Referenced by ComputeMultipleInelasticCosseratStress::computeAdmissibleState(), ComputeMultipleInelasticStress::computeAdmissibleState(), ComputeStressBase::computeQpProperties(), ComputeStrainIncrementBasedStress::computeQpStress(), ComputeLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeDamageStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), FiniteStrainPlasticMaterial::computeQpStress(), ComputeLinearElasticPFFractureStress::computeQpStress(), computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration(), ComputeLinearElasticPFFractureStress::computeStrainSpectral(), ComputeLinearElasticPFFractureStress::computeStrainVolDev(), ComputeLinearElasticPFFractureStress::computeStressSpectral(), ComputeMultipleInelasticStress::finiteStrainRotation(), ComputeStressBase::initQpStatefulProperties(), FiniteStrainCrystalPlasticity::initQpStatefulProperties(), FiniteStrainUObasedCP::initQpStatefulProperties(), FiniteStrainHyperElasticViscoPlastic::initQpStatefulProperties(), postReturnMap(), FiniteStrainUObasedCP::postSolveQp(), FiniteStrainHyperElasticViscoPlastic::postSolveQp(), FiniteStrainCrystalPlasticity::postSolveQp(), ComputeSmearedCrackingStress::updateCrackingStateAndStress(), ComputeMultipleInelasticStress::updateQpState(), and ComputeMultipleInelasticStress::updateQpStateSingleModel().

◆ _stress_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_stress_old
protected

Old value of stress.

Definition at line 149 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 143 of file ComputeMultiPlasticityStress.h.

◆ _yf

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

yield functions

Definition at line 119 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().


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