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

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

#include <ComputeMultiPlasticityStress.h>

Inheritance diagram for ComputeMultiPlasticityStress:
[legend]

Public Member Functions

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

Protected Types

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

Protected Member Functions

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

Protected Attributes

unsigned int _max_iter
 Maximum number of Newton-Raphson iterations allowed. More...
 
Real _min_stepsize
 Minimum fraction of applied strain that may be applied during adaptive stepsizing. More...
 
Real _max_stepsize_for_dumb
 "dumb" deactivation will only be used if the stepsize falls below this quantity More...
 
bool _ignore_failures
 Even if the returnMap fails, return the best values found for stress and internal parameters. More...
 
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
 
Real _epp_tol
 Tolerance on the plastic strain increment ("direction") constraint. More...
 
std::vector< Real > _dummy_pm
 dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStress More...
 
std::vector< Real > _cumulative_pm
 the sum of the plastic multipliers over all the sub-steps. More...
 
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
 
bool _n_supplied
 User supplied the transverse direction vector. More...
 
RealVectorValue _n_input
 the supplied transverse direction vector More...
 
RealTensorValue _rot
 rotation matrix that takes _n to (0, 0, 1) More...
 
bool _perform_finite_strain_rotations
 whether to perform the rotations necessary in finite-strain simulations More...
 
MaterialProperty< RankTwoTensor > & _plastic_strain
 plastic strain More...
 
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
 Old value of plastic strain. More...
 
MaterialProperty< std::vector< Real > > & _intnl
 internal parameters More...
 
const MaterialProperty< std::vector< Real > > & _intnl_old
 old values of internal parameters More...
 
MaterialProperty< std::vector< Real > > & _yf
 yield functions More...
 
MaterialProperty< Real > & _iter
 Number of Newton-Raphson iterations used in the return-map. More...
 
MaterialProperty< Real > & _linesearch_needed
 Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< Real > & _ld_encountered
 Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< Real > & _constraints_added
 Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise) More...
 
MaterialProperty< RealVectorValue > & _n
 current value of transverse direction More...
 
const MaterialProperty< RealVectorValue > & _n_old
 old value of transverse direction More...
 
const MaterialProperty< RankTwoTensor > & _strain_increment
 strain increment (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _total_strain_old
 Old value of total strain (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _rotation_increment
 Rotation increment (coming from ComputeIncrementalSmallStrain, for example) More...
 
const MaterialProperty< RankTwoTensor > & _stress_old
 Old value of stress. More...
 
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
 Old value of elastic strain. More...
 
bool _cosserat
 whether Cosserat mechanics should be used More...
 
const MaterialProperty< RankTwoTensor > * _curvature
 The Cosserat curvature strain. More...
 
const MaterialProperty< RankFourTensor > * _elastic_flexural_rigidity_tensor
 The Cosserat elastic flexural rigidity tensor. More...
 
MaterialProperty< RankTwoTensor > * _couple_stress
 the Cosserat couple-stress More...
 
const MaterialProperty< RankTwoTensor > * _couple_stress_old
 the old value of Cosserat couple-stress More...
 
MaterialProperty< RankFourTensor > * _Jacobian_mult_couple
 derivative of couple-stress w.r.t. curvature More...
 
RankFourTensor _my_elasticity_tensor
 Elasticity tensor that can be rotated by this class (ie, its not const) More...
 
RankTwoTensor _my_strain_increment
 Strain increment that can be rotated by this class, and split into multiple increments (ie, its not const) More...
 
RankFourTensor _my_flexural_rigidity_tensor
 Flexual rigidity tensor that can be rotated by this class (ie, its not const) More...
 
RankTwoTensor _my_curvature
 Curvature that can be rotated by this class, and split into multiple increments (ie, its not const) More...
 
const std::string _base_name
 
const std::string _elasticity_tensor_name
 
const MaterialProperty< RankTwoTensor > & _mechanical_strain
 
MaterialProperty< RankTwoTensor > & _stress
 
MaterialProperty< RankTwoTensor > & _elastic_strain
 
const MaterialProperty< RankFourTensor > & _elasticity_tensor
 
const MaterialProperty< RankTwoTensor > & _extra_stress
 Extra stress tensor. More...
 
std::vector< Function * > _initial_stress_fcn
 initial stress components More...
 
MaterialProperty< RankFourTensor > & _Jacobian_mult
 derivative of stress w.r.t. strain (_dstress_dstrain) More...
 
const bool _store_stress_old
 Parameter which decides whether to store old stress. This is required for HHT time integration and Rayleigh damping. More...
 
const bool _initial_stress_provided
 Whether initial stress was provided. InitialStress Deprecation: remove this. More...
 
MaterialProperty< RankTwoTensor > * _initial_stress
 Initial stress, if provided. InitialStress Deprecation: remove this. More...
 
const MaterialProperty< RankTwoTensor > * _initial_stress_old
 Old value of initial stress, which is needed to correctly implement finite-strain rotations. InitialStress Deprecation: remove this. 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)
 

Private Attributes

bool & _step_one
 True if this is the first timestep (timestep < 2). More...
 

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 28 of file ComputeMultiPlasticityStress.h.

Member Enumeration Documentation

The functions from which quickStep can be called.

Enumerator
computeQpStress_function 
returnMap_function 

Definition at line 436 of file ComputeMultiPlasticityStress.h.

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

Enumerator
elastic 
linear 
nonlinear 

Definition at line 50 of file ComputeMultiPlasticityStress.h.

Constructor & Destructor Documentation

ComputeMultiPlasticityStress::ComputeMultiPlasticityStress ( const InputParameters &  parameters)

Definition at line 93 of file ComputeMultiPlasticityStress.C.

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

Member Function Documentation

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

Definition at line 1578 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1579 {
1580  unsigned num = 0;
1581  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1582  if (act[surface])
1583  num += (1 << surface); // (1 << x) = 2^x
1584 
1585  return num;
1586 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
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 809 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().

813 {
814  active_surfaces_of_model.resize(0);
815  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816  if (active[_surfaces_given_model[model][model_surface]])
817  active_surfaces_of_model.push_back(model_surface);
818 }
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.
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 798 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

801 {
802  active_surfaces.resize(0);
803  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804  if (active[_surfaces_given_model[model][model_surface]])
805  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 }
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.
void ComputeStressBase::addQpInitialStress ( )
protectedinherited

InitialStress Deprecation: remove this method.

Adds initial stress, if it is provided, to _stress[_qp]. This function should NOT be used if you calculate stress using

stress = stress_old + elasticity * strain_increment

because stress_old will already include initial stress. However this function SHOULD be used if your code uses

stress = elasticity * (elastic_strain_old + strain_increment) or stress = elasticity * elastic_strain

since in these cases the elastic_strain and elastic_strain_old will not include any contribution from initial stress.

Definition at line 112 of file ComputeStressBase.C.

Referenced by ComputeLinearElasticStress::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration(), ComputeStressBase::initQpStatefulProperties(), ComputeMultipleInelasticStress::updateQpState(), and ComputeMultipleInelasticStress::updateQpStateSingleModel().

113 {
115  _stress[_qp] += (*_initial_stress)[_qp];
116 }
MaterialProperty< RankTwoTensor > & _stress
const bool _initial_stress_provided
Whether initial stress was provided. InitialStress Deprecation: remove this.
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 789 of file MultiPlasticityRawComponentAssembler.C.

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

790 {
791  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792  if (active[_surfaces_given_model[model][model_surface]])
793  return true;
794  return false;
795 }
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.
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 1317 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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

Referenced by returnMap().

344 {
345  mooseAssert(f.size() == _num_surfaces,
346  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
347  << _num_surfaces
348  << " surfaces");
349  mooseAssert(intnl.size() == _num_models,
350  "buildActiveConstraints called with intnl.size = " << intnl.size()
351  << " while there are "
352  << _num_models
353  << " models");
354 
355  if (_specialIC == "rock")
356  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357  else if (_specialIC == "joint")
358  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359  else // no specialIC
360  {
361  act.resize(0);
362  unsigned ind = 0;
363  for (unsigned model = 0; model < _num_models; ++model)
364  {
365  std::vector<Real> model_f(0);
366  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367  model_f.push_back(f[ind++]);
368  std::vector<bool> model_act;
369  RankTwoTensor returned_stress;
370  _f[model]->activeConstraints(
371  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373  act.push_back(model_act[model_surface]);
374  }
375  }
376 }
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.
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 1528 of file ComputeMultiPlasticityStress.C.

Referenced by changeScheme(), and returnMap().

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

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

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

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

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

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

292 {
293  // see comments at the start of .h file
294 
295  mooseAssert(intnl_old.size() == _num_models,
296  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
297  mooseAssert(intnl.size() == _num_models,
298  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
299  mooseAssert(pm.size() == _num_surfaces,
300  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
301  mooseAssert(active.size() == _num_surfaces,
302  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
303 
304  std::vector<Real> f; // the yield functions
305  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
306  std::vector<Real> ic; // the "internal constraints"
307 
308  std::vector<RankTwoTensor> r;
309  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
310 
311  if (eliminate_ld)
312  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
313  else
314  deactivated_due_to_ld.assign(_num_surfaces, false);
315 
316  std::vector<bool> active_not_deact(_num_surfaces);
317  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
318  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
319 
320  unsigned num_active_f = 0;
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  if (active_not_deact[surface])
323  num_active_f++;
324 
325  unsigned num_active_ic = 0;
326  for (unsigned model = 0; model < _num_models; ++model)
327  if (anyActiveSurfaces(model, active_not_deact))
328  num_active_ic++;
329 
330  unsigned int dim = 3;
331  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
332  // num_active_f comes from "f",
333  // num_active_f comes from "ic"
334 
335  rhs.resize(system_size);
336 
337  unsigned ind = 0;
338  for (unsigned i = 0; i < dim; ++i)
339  for (unsigned j = 0; j <= i; ++j)
340  rhs[ind++] = -epp(i, j);
341  unsigned active_surface = 0;
342  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
343  if (active[surface])
344  {
345  if (!deactivated_due_to_ld[surface])
346  rhs[ind++] = -f[active_surface];
347  active_surface++;
348  }
349  unsigned active_model = 0;
350  for (unsigned model = 0; model < _num_models; ++model)
351  if (anyActiveSurfaces(model, active))
352  {
353  if (anyActiveSurfaces(model, active_not_deact))
354  rhs[ind++] = -ic[active_model];
355  active_model++;
356  }
357 
358  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
359 }
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.
bool ComputeMultiPlasticityStress::canAddConstraints ( const std::vector< bool > &  act,
const std::vector< Real > &  all_f 
)
protected

Definition at line 1019 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1021 {
1022  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1023  if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1024  return true;
1025  return false;
1026 }
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 ComputeMultiPlasticityStress::canChangeScheme ( DeactivationSchemeEnum  current_deactivation_scheme,
bool  can_revert_to_dumb 
)
protected

Definition at line 1029 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1031 {
1032  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
1033  return true;
1034 
1035  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
1036  return true;
1037 
1038  if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
1039  can_revert_to_dumb)
1040  return true;
1041 
1042  if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
1043  can_revert_to_dumb)
1044  return true;
1045 
1046  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1047  can_revert_to_dumb)
1048  return true;
1049 
1050  return false;
1051 }
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
bool ComputeMultiPlasticityStress::canIncrementDumb ( int  dumb_iteration)
protected

Definition at line 1571 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1572 {
1573  // (1 << _num_surfaces) = 2^_num_surfaces
1574  return ((dumb_iteration + 1) < (1 << _num_surfaces));
1575 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
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 1054 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1062 {
1063  if (current_deactivation_scheme == optimized &&
1066  {
1067  current_deactivation_scheme = safe;
1068  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1069  act[surface] = initial_act[surface];
1070  }
1071  else if ((current_deactivation_scheme == safe &&
1074  can_revert_to_dumb) ||
1075  (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1076  can_revert_to_dumb))
1077  {
1078  current_deactivation_scheme = dumb;
1079  dumb_iteration = 0;
1080  buildDumbOrder(initial_stress, intnl_old, dumb_order);
1081  incrementDumb(dumb_iteration, dumb_order, act);
1082  }
1083 }
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
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 1275 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1278 {
1279  std::vector<bool> act;
1280  act.assign(_num_surfaces, true);
1281 
1282  yieldFunction(stress, intnl, act, all_f);
1283 
1284  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1285  if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1286  return false;
1287 
1288  return true;
1289 }
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.
void MultiPlasticityDebugger::checkDerivatives ( )
inherited

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

Definition at line 80 of file MultiPlasticityDebugger.C.

Referenced by initQpStatefulProperties(), and plasticStep().

81 {
82  Moose::err
83  << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
85 
86  std::vector<bool> act;
87  act.assign(_num_surfaces, true);
88 
89  Moose::err << "\ndyieldFunction_dstress. Relative L2 norms.\n";
90  std::vector<RankTwoTensor> df_dstress;
91  std::vector<RankTwoTensor> fddf_dstress;
94  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
95  {
96  Moose::err << "surface = " << surface << " Relative L2norm = "
97  << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
98  (df_dstress[surface] + fddf_dstress[surface]).L2norm()
99  << "\n";
100  Moose::err << "Coded:\n";
101  df_dstress[surface].print();
102  Moose::err << "Finite difference:\n";
103  fddf_dstress[surface].print();
104  }
105 
106  Moose::err << "\ndyieldFunction_dintnl.\n";
107  std::vector<Real> df_dintnl;
109  Moose::err << "Coded:\n";
110  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
111  Moose::err << df_dintnl[surface] << " ";
112  Moose::err << "\n";
113  std::vector<Real> fddf_dintnl;
115  Moose::err << "Finite difference:\n";
116  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
117  Moose::err << fddf_dintnl[surface] << " ";
118  Moose::err << "\n";
119 
120  Moose::err << "\ndflowPotential_dstress. Relative L2 norms.\n";
121  std::vector<RankFourTensor> dr_dstress;
122  std::vector<RankFourTensor> fddr_dstress;
125  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
126  {
127  Moose::err << "surface = " << surface << " Relative L2norm = "
128  << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
129  (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
130  << "\n";
131  Moose::err << "Coded:\n";
132  dr_dstress[surface].print();
133  Moose::err << "Finite difference:\n";
134  fddr_dstress[surface].print();
135  }
136 
137  Moose::err << "\ndflowPotential_dintnl. Relative L2 norms.\n";
138  std::vector<RankTwoTensor> dr_dintnl;
139  std::vector<RankTwoTensor> fddr_dintnl;
142  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
143  {
144  Moose::err << "surface = " << surface << " Relative L2norm = "
145  << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
146  (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
147  << "\n";
148  Moose::err << "Coded:\n";
149  dr_dintnl[surface].print();
150  Moose::err << "Finite difference:\n";
151  fddr_dintnl[surface].print();
152  }
153 }
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...
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 156 of file MultiPlasticityDebugger.C.

Referenced by computeQpStress(), and plasticStep().

158 {
159  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
161 
162  std::vector<bool> act;
163  act.assign(_num_surfaces, true);
164  std::vector<bool> deactivated_due_to_ld;
165  deactivated_due_to_ld.assign(_num_surfaces, false);
166 
167  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
168 
169  std::vector<std::vector<Real>> jac;
170  calculateJacobian(_fspb_debug_stress,
173  E_inv,
174  act,
175  deactivated_due_to_ld,
176  jac);
177 
178  std::vector<std::vector<Real>> fdjac;
179  fdJacobian(_fspb_debug_stress,
180  intnl_old,
183  delta_dp,
184  E_inv,
185  false,
186  fdjac);
187 
188  Real L2_numer = 0;
189  Real L2_denom = 0;
190  for (unsigned row = 0; row < jac.size(); ++row)
191  for (unsigned col = 0; col < jac.size(); ++col)
192  {
193  L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
194  L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
195  }
196  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
197 
198  Moose::err << "\nHand-coded Jacobian:\n";
199  for (unsigned row = 0; row < jac.size(); ++row)
200  {
201  for (unsigned col = 0; col < jac.size(); ++col)
202  Moose::err << jac[row][col] << " ";
203  Moose::err << "\n";
204  }
205 
206  Moose::err << "Finite difference Jacobian:\n";
207  for (unsigned row = 0; row < fdjac.size(); ++row)
208  {
209  for (unsigned col = 0; col < fdjac.size(); ++col)
210  Moose::err << fdjac[row][col] << " ";
211  Moose::err << "\n";
212  }
213 }
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.
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 1292 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

1295 {
1296  unsigned ind = 0;
1297  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1298  {
1299  if (active[surface])
1300  {
1301  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1302  if (pm[surface] != 0)
1303  return false;
1304  }
1305  else if (pm[surface] != 0)
1306  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1307  "coding!!");
1308  }
1309  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1310  if (pm[surface] < 0)
1311  return false;
1312 
1313  return true;
1314 }
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.
void MultiPlasticityDebugger::checkSolution ( const RankFourTensor &  E_inv)
inherited

Checks that Ax does equal b in the NR procedure.

Definition at line 362 of file MultiPlasticityDebugger.C.

Referenced by computeQpStress(), and plasticStep().

363 {
364  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
365  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
367 
368  std::vector<bool> act;
369  act.assign(_num_surfaces, true);
370  std::vector<bool> deactivated_due_to_ld;
371  deactivated_due_to_ld.assign(_num_surfaces, false);
372 
373  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
374 
375  std::vector<Real> orig_rhs;
376  calculateRHS(_fspb_debug_stress,
380  delta_dp,
381  orig_rhs,
382  act,
383  true,
384  deactivated_due_to_ld);
385 
386  Moose::err << "\nb = ";
387  for (unsigned i = 0; i < orig_rhs.size(); ++i)
388  Moose::err << orig_rhs[i] << " ";
389  Moose::err << "\n\n";
390 
391  std::vector<std::vector<Real>> jac_coded;
392  calculateJacobian(_fspb_debug_stress,
395  E_inv,
396  act,
397  deactivated_due_to_ld,
398  jac_coded);
399 
400  Moose::err
401  << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
402  Moose::err
403  << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
404  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
405  "linear dependence.\n";
406  Moose::err << "Hand-coded Jacobian:\n";
407  for (unsigned row = 0; row < jac_coded.size(); ++row)
408  {
409  for (unsigned col = 0; col < jac_coded.size(); ++col)
410  Moose::err << jac_coded[row][col] << " ";
411  Moose::err << "\n";
412  }
413 
414  deactivated_due_to_ld.assign(_num_surfaces,
415  false); // this potentially gets changed by nrStep, below
416  RankTwoTensor dstress;
417  std::vector<Real> dpm;
418  std::vector<Real> dintnl;
419  nrStep(_fspb_debug_stress,
423  E_inv,
424  delta_dp,
425  dstress,
426  dpm,
427  dintnl,
428  act,
429  deactivated_due_to_ld);
430 
431  std::vector<bool> active_not_deact(_num_surfaces);
432  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
433  active_not_deact[surface] = !deactivated_due_to_ld[surface];
434 
435  std::vector<Real> x;
436  x.assign(orig_rhs.size(), 0);
437  unsigned ind = 0;
438  for (unsigned i = 0; i < 3; ++i)
439  for (unsigned j = 0; j <= i; ++j)
440  x[ind++] = dstress(i, j);
441  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
442  if (active_not_deact[surface])
443  x[ind++] = dpm[surface];
444  for (unsigned model = 0; model < _num_models; ++model)
445  if (anyActiveSurfaces(model, active_not_deact))
446  x[ind++] = dintnl[model];
447 
448  mooseAssert(ind == orig_rhs.size(),
449  "Incorrect extracting of changes from NR solution in the "
450  "finite-difference checking of nrStep");
451 
452  Moose::err << "\nThis yields x =";
453  for (unsigned i = 0; i < orig_rhs.size(); ++i)
454  Moose::err << x[i] << " ";
455  Moose::err << "\n";
456 
457  std::vector<std::vector<Real>> jac_fd;
458  fdJacobian(_fspb_debug_stress,
462  delta_dp,
463  E_inv,
464  true,
465  jac_fd);
466 
467  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
468  Moose::err << "in order to check that the solution is correct\n";
469  Moose::err << "Finite-difference Jacobian:\n";
470  for (unsigned row = 0; row < jac_fd.size(); ++row)
471  {
472  for (unsigned col = 0; col < jac_fd.size(); ++col)
473  Moose::err << jac_fd[row][col] << " ";
474  Moose::err << "\n";
475  }
476 
477  Real L2_numer = 0;
478  Real L2_denom = 0;
479  for (unsigned row = 0; row < jac_coded.size(); ++row)
480  for (unsigned col = 0; col < jac_coded.size(); ++col)
481  {
482  L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
483  L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
484  }
485  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
486  << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
487 
488  std::vector<Real> fd_times_x;
489  fd_times_x.assign(orig_rhs.size(), 0);
490  for (unsigned row = 0; row < orig_rhs.size(); ++row)
491  for (unsigned col = 0; col < orig_rhs.size(); ++col)
492  fd_times_x[row] += jac_fd[row][col] * x[col];
493 
494  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
495  for (unsigned i = 0; i < orig_rhs.size(); ++i)
496  Moose::err << fd_times_x[i] << " ";
497  Moose::err << "\n";
498  Moose::err << "Recall that b = \n";
499  for (unsigned i = 0; i < orig_rhs.size(); ++i)
500  Moose::err << orig_rhs[i] << " ";
501  Moose::err << "\n";
502 
503  L2_numer = 0;
504  L2_denom = 0;
505  for (unsigned i = 0; i < orig_rhs.size(); ++i)
506  {
507  L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
508  L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
509  }
510  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
511 }
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.
void ComputeStressBase::computeQpProperties ( )
overrideprotectedvirtualinherited

Reimplemented in ComputeBirchMurnaghanEquationOfStress, and ComputeStressEosBase.

Definition at line 99 of file ComputeStressBase.C.

100 {
101  // InitialStress Deprecation: remove the following 2 lines
103  (*_initial_stress)[_qp] = (*_initial_stress_old)[_qp];
104 
105  computeQpStress();
106 
107  // Add in extra stress
108  _stress[_qp] += _extra_stress[_qp];
109 }
virtual void computeQpStress()=0
MaterialProperty< RankTwoTensor > & _stress
const bool _initial_stress_provided
Whether initial stress was provided. InitialStress Deprecation: remove this.
const MaterialProperty< RankTwoTensor > & _extra_stress
Extra stress tensor.
void ComputeMultiPlasticityStress::computeQpStress ( )
protectedvirtual

Implements ComputeStressBase.

Definition at line 214 of file ComputeMultiPlasticityStress.C.

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

Referenced by quickStep(), and returnMap().

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

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

234 {
235  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
236  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
237 
238  dr_dintnl.resize(0);
239  std::vector<unsigned int> active_surfaces_of_model;
240  std::vector<unsigned int>::iterator active_surface;
241  std::vector<RankTwoTensor> model_dr_dintnl;
242  for (unsigned model = 0; model < _num_models; ++model)
243  {
244  activeModelSurfaces(model, active, active_surfaces_of_model);
245  if (active_surfaces_of_model.size() > 0)
246  {
247  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
248  for (active_surface = active_surfaces_of_model.begin();
249  active_surface != active_surfaces_of_model.end();
250  ++active_surface)
251  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
252  }
253  }
254 }
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.
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 202 of file MultiPlasticityRawComponentAssembler.C.

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

207 {
208  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
209  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
210 
211  dr_dstress.resize(0);
212  std::vector<unsigned int> active_surfaces_of_model;
213  std::vector<unsigned int>::iterator active_surface;
214  std::vector<RankFourTensor> model_dr_dstress;
215  for (unsigned model = 0; model < _num_models; ++model)
216  {
217  activeModelSurfaces(model, active, active_surfaces_of_model);
218  if (active_surfaces_of_model.size() > 0)
219  {
220  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
221  for (active_surface = active_surfaces_of_model.begin();
222  active_surface != active_surfaces_of_model.end();
223  ++active_surface)
224  dr_dstress.push_back(model_dr_dstress[*active_surface]);
225  }
226  }
227 }
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.
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 312 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

316 {
317  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
318  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
319 
320  dh_dintnl.resize(0);
321  std::vector<unsigned int> active_surfaces_of_model;
322  std::vector<unsigned int>::iterator active_surface;
323  std::vector<Real> model_dh_dintnl;
324  for (unsigned model = 0; model < _num_models; ++model)
325  {
326  activeModelSurfaces(model, active, active_surfaces_of_model);
327  if (active_surfaces_of_model.size() > 0)
328  {
329  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
330  for (active_surface = active_surfaces_of_model.begin();
331  active_surface != active_surfaces_of_model.end();
332  ++active_surface)
333  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
334  }
335  }
336 }
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.
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 284 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

289 {
290  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
291  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
292 
293  dh_dstress.resize(0);
294  std::vector<unsigned int> active_surfaces_of_model;
295  std::vector<unsigned int>::iterator active_surface;
296  std::vector<RankTwoTensor> model_dh_dstress;
297  for (unsigned model = 0; model < _num_models; ++model)
298  {
299  activeModelSurfaces(model, active, active_surfaces_of_model);
300  if (active_surfaces_of_model.size() > 0)
301  {
302  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
303  for (active_surface = active_surfaces_of_model.begin();
304  active_surface != active_surfaces_of_model.end();
305  ++active_surface)
306  dh_dstress.push_back(model_dh_dstress[*active_surface]);
307  }
308  }
309 }
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.
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 148 of file MultiPlasticityRawComponentAssembler.C.

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

152 {
153  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
154  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
155 
156  df_dintnl.resize(0);
157  std::vector<unsigned int> active_surfaces_of_model;
158  std::vector<unsigned int>::iterator active_surface;
159  std::vector<Real> model_df_dintnl;
160  for (unsigned model = 0; model < _num_models; ++model)
161  {
162  activeModelSurfaces(model, active, active_surfaces_of_model);
163  if (active_surfaces_of_model.size() > 0)
164  {
165  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
166  for (active_surface = active_surfaces_of_model.begin();
167  active_surface != active_surfaces_of_model.end();
168  ++active_surface)
169  df_dintnl.push_back(model_df_dintnl[*active_surface]);
170  }
171  }
172 }
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.
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 120 of file MultiPlasticityRawComponentAssembler.C.

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

125 {
126  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
127  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
128 
129  df_dstress.resize(0);
130  std::vector<unsigned int> active_surfaces_of_model;
131  std::vector<unsigned int>::iterator active_surface;
132  std::vector<RankTwoTensor> model_df_dstress;
133  for (unsigned model = 0; model < _num_models; ++model)
134  {
135  activeModelSurfaces(model, active, active_surfaces_of_model);
136  if (active_surfaces_of_model.size() > 0)
137  {
138  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
139  for (active_surface = active_surfaces_of_model.begin();
140  active_surface != active_surfaces_of_model.end();
141  ++active_surface)
142  df_dstress.push_back(model_df_dstress[*active_surface]);
143  }
144  }
145 }
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.
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 175 of file MultiPlasticityRawComponentAssembler.C.

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

179 {
180  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
181  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
182 
183  r.resize(0);
184  std::vector<unsigned int> active_surfaces_of_model;
185  std::vector<unsigned int>::iterator active_surface;
186  std::vector<RankTwoTensor> model_r;
187  for (unsigned model = 0; model < _num_models; ++model)
188  {
189  activeModelSurfaces(model, active, active_surfaces_of_model);
190  if (active_surfaces_of_model.size() > 0)
191  {
192  _f[model]->flowPotentialV(stress, intnl[model], model_r);
193  for (active_surface = active_surfaces_of_model.begin();
194  active_surface != active_surfaces_of_model.end();
195  ++active_surface)
196  r.push_back(model_r[*active_surface]);
197  }
198  }
199 }
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.
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 257 of file MultiPlasticityRawComponentAssembler.C.

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

261 {
262  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
263  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
264 
265  h.resize(0);
266  std::vector<unsigned int> active_surfaces_of_model;
267  std::vector<unsigned int>::iterator active_surface;
268  std::vector<Real> model_h;
269  for (unsigned model = 0; model < _num_models; ++model)
270  {
271  activeModelSurfaces(model, active, active_surfaces_of_model);
272  if (active_surfaces_of_model.size() > 0)
273  {
274  _f[model]->hardPotentialV(stress, intnl[model], model_h);
275  for (active_surface = active_surfaces_of_model.begin();
276  active_surface != active_surfaces_of_model.end();
277  ++active_surface)
278  h.push_back(model_h[*active_surface]);
279  }
280  }
281 }
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.
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 1559 of file ComputeMultiPlasticityStress.C.

Referenced by changeScheme(), and returnMap().

1562 {
1563  dumb_iteration += 1;
1564  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1565  act[dumb_order[surface]] =
1566  (dumb_iteration &
1567  (1 << surface)); // returns true if the surface_th bit of dumb_iteration == 1
1568 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void ComputeMultiPlasticityStress::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 184 of file ComputeMultiPlasticityStress.C.

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

Referenced by singleStep().

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

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

621 {
622  // Calculate RHS and Jacobian
623  std::vector<Real> rhs;
624  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
625 
626  std::vector<std::vector<Real>> jac;
627  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
628 
629  // prepare for LAPACKgesv_ routine provided by PETSc
630  int system_size = rhs.size();
631 
632  std::vector<double> a(system_size * system_size);
633  // Fill in the a "matrix" by going down columns
634  unsigned ind = 0;
635  for (int col = 0; col < system_size; ++col)
636  for (int row = 0; row < system_size; ++row)
637  a[ind++] = jac[row][col];
638 
639  int nrhs = 1;
640  std::vector<int> ipiv(system_size);
641  int info;
642  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
643 
644  if (info != 0)
645  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
646  "routine returned with error code ",
647  info);
648 
649  // Extract the results back to dstress, dpm and dintnl
650  std::vector<bool> active_not_deact(_num_surfaces);
651  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
652  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
653 
654  unsigned int dim = 3;
655  ind = 0;
656 
657  for (unsigned i = 0; i < dim; ++i)
658  for (unsigned j = 0; j <= i; ++j)
659  dstress(i, j) = dstress(j, i) = rhs[ind++];
660  dpm.assign(_num_surfaces, 0);
661  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
662  if (active_not_deact[surface])
663  dpm[surface] = rhs[ind++];
664  dintnl.assign(_num_models, 0);
665  for (unsigned model = 0; model < _num_models; ++model)
666  if (anyActiveSurfaces(model, active_not_deact))
667  dintnl[model] = rhs[ind++];
668 
669  mooseAssert(static_cast<int>(ind) == system_size,
670  "Incorrect extracting of changes from NR solution in nrStep");
671 }
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.
unsigned ComputeMultiPlasticityStress::numberActive ( const std::vector< bool > &  active)
protectedvirtual

counts the number of active constraints

Definition at line 1265 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap(), and singleStep().

1266 {
1267  unsigned num_active = 0;
1268  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1269  if (active[surface])
1270  num_active++;
1271  return num_active;
1272 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
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 50 of file MultiPlasticityDebugger.C.

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

51 {
52  Moose::err << "Debug Parameters are as follows\n";
53  Moose::err << "stress = \n";
54  _fspb_debug_stress.print();
55 
56  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
59  mooseError("The debug parameters have the wrong size\n");
60 
61  Moose::err << "plastic multipliers =\n";
62  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
63  Moose::err << _fspb_debug_pm[surface] << "\n";
64 
65  Moose::err << "internal parameters =\n";
66  for (unsigned model = 0; model < _num_models; ++model)
67  Moose::err << _fspb_debug_intnl[model] << "\n";
68 
69  Moose::err << "finite-differencing parameter for stress-changes:\n"
70  << _fspb_debug_stress_change << "\n";
71  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
72  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
73  Moose::err << _fspb_debug_pm_change[surface] << "\n";
74  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
75  for (unsigned model = 0; model < _num_models; ++model)
76  Moose::err << _fspb_debug_intnl_change[model] << "\n";
77 }
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.
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 461 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

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

Definition at line 332 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

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

Definition at line 313 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

314 {
315  if (_n_supplied)
316  {
317  // this is a rotation matrix that will rotate _n to the "z" axis
318  _rot = RotationMatrix::rotVecToZ(_n[_qp]);
319 
320  // rotate the tensors to this frame
321  _my_elasticity_tensor.rotate(_rot);
322  _my_strain_increment.rotate(_rot);
323  if (_cosserat)
324  {
326  _my_curvature.rotate(_rot);
327  }
328  }
329 }
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
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

Definition at line 364 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress(), and returnMap().

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

Referenced by singleStep().

1254 {
1255  bool reinstated_actives = false;
1256  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1257  if (deactivated_due_to_ld[surface])
1258  reinstated_actives = true;
1259 
1260  deactivated_due_to_ld.assign(_num_surfaces, false);
1261  return reinstated_actives;
1262 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
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 1358 of file ComputeMultiPlasticityStress.C.

Referenced by lineSearch(), and singleStep().

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

Referenced by plasticStep().

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

Definition at line 597 of file MultiPlasticityRawComponentAssembler.C.

Referenced by quickStep().

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

Definition at line 305 of file ComputeMultiPlasticityStress.C.

Referenced by computeQpStress().

306 {
307  if (!_n_supplied)
308  return tens;
309  return tens.rotated(_rot);
310 }
bool _n_supplied
User supplied the transverse direction vector.
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
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 1086 of file ComputeMultiPlasticityStress.C.

Referenced by returnMap().

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

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

97 {
98  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
99  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
100 
101  f.resize(0);
102  std::vector<unsigned int> active_surfaces_of_model;
103  std::vector<unsigned int>::iterator active_surface;
104  std::vector<Real> model_f;
105  for (unsigned model = 0; model < _num_models; ++model)
106  {
107  activeModelSurfaces(model, active, active_surfaces_of_model);
108  if (active_surfaces_of_model.size() > 0)
109  {
110  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
111  for (active_surface = active_surfaces_of_model.begin();
112  active_surface != active_surfaces_of_model.end();
113  ++active_surface)
114  f.push_back(model_f[*active_surface]);
115  }
116  }
117 }
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

const std::string ComputeStressBase::_base_name
protectedinherited
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 123 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

bool ComputeMultiPlasticityStress::_cosserat
protected

whether Cosserat mechanics should be used

Definition at line 147 of file ComputeMultiPlasticityStress.h.

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

MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_couple_stress
protected

the Cosserat couple-stress

Definition at line 156 of file ComputeMultiPlasticityStress.h.

const MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_couple_stress_old
protected

the old value of Cosserat couple-stress

Definition at line 159 of file ComputeMultiPlasticityStress.h.

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 67 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

const MaterialProperty<RankTwoTensor>* ComputeMultiPlasticityStress::_curvature
protected

The Cosserat curvature strain.

Definition at line 150 of file ComputeMultiPlasticityStress.h.

enum ComputeMultiPlasticityStress::DeactivationSchemeEnum ComputeMultiPlasticityStress::_deactivation_scheme
protected
std::vector<Real> ComputeMultiPlasticityStress::_dummy_pm
protected

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

Definition at line 61 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

const MaterialProperty<RankFourTensor>* ComputeMultiPlasticityStress::_elastic_flexural_rigidity_tensor
protected

The Cosserat elastic flexural rigidity tensor.

Definition at line 153 of file ComputeMultiPlasticityStress.h.

MaterialProperty<RankTwoTensor>& ComputeStressBase::_elastic_strain
protectedinherited
const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_elastic_strain_old
protected

Old value of elastic strain.

Definition at line 144 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

const MaterialProperty<RankFourTensor>& ComputeStressBase::_elasticity_tensor
protectedinherited
const std::string ComputeStressBase::_elasticity_tensor_name
protectedinherited
Real ComputeMultiPlasticityStress::_epp_tol
protected

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

Definition at line 58 of file ComputeMultiPlasticityStress.h.

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

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_extra_stress
protectedinherited
std::vector<const TensorMechanicsPlasticModel *> MultiPlasticityRawComponentAssembler::_f
protectedinherited
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 58 of file MultiPlasticityDebugger.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl
protectedinherited
std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change
protectedinherited
std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm
protectedinherited
std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change
protectedinherited

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 73 of file MultiPlasticityDebugger.h.

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

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress
protectedinherited
Real MultiPlasticityDebugger::_fspb_debug_stress_change
protectedinherited
bool ComputeMultiPlasticityStress::_ignore_failures
protected

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

Definition at line 47 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

MaterialProperty<RankTwoTensor>* ComputeStressBase::_initial_stress
protectedinherited

Initial stress, if provided. InitialStress Deprecation: remove this.

Definition at line 74 of file ComputeStressBase.h.

Referenced by ComputeStressBase::initQpStatefulProperties().

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

initial stress components

Definition at line 62 of file ComputeStressBase.h.

Referenced by ComputeStressBase::ComputeStressBase(), and ComputeStressBase::initQpStatefulProperties().

const MaterialProperty<RankTwoTensor>* ComputeStressBase::_initial_stress_old
protectedinherited

Old value of initial stress, which is needed to correctly implement finite-strain rotations. InitialStress Deprecation: remove this.

Definition at line 77 of file ComputeStressBase.h.

const bool ComputeStressBase::_initial_stress_provided
protectedinherited
MaterialProperty<std::vector<Real> >& ComputeMultiPlasticityStress::_intnl
protected

internal parameters

Definition at line 105 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

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

old values of internal parameters

Definition at line 108 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

MaterialProperty<Real>& ComputeMultiPlasticityStress::_iter
protected

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

Definition at line 114 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

MaterialProperty<RankFourTensor>& ComputeStressBase::_Jacobian_mult
protectedinherited
MaterialProperty<RankFourTensor>* ComputeMultiPlasticityStress::_Jacobian_mult_couple
protected

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

Definition at line 162 of file ComputeMultiPlasticityStress.h.

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 120 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

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 117 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

unsigned int ComputeMultiPlasticityStress::_max_iter
protected

Maximum number of Newton-Raphson iterations allowed.

Definition at line 38 of file ComputeMultiPlasticityStress.h.

Referenced by returnMap().

Real ComputeMultiPlasticityStress::_max_stepsize_for_dumb
protected

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

Definition at line 44 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_mechanical_strain
protectedinherited
Real MultiPlasticityLinearSystem::_min_f_tol
protectedinherited

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

Definition at line 132 of file MultiPlasticityLinearSystem.h.

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

Real ComputeMultiPlasticityStress::_min_stepsize
protected

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

Definition at line 41 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

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 174 of file ComputeMultiPlasticityStress.h.

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

RankFourTensor ComputeMultiPlasticityStress::_my_elasticity_tensor
protected

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

Definition at line 165 of file ComputeMultiPlasticityStress.h.

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

RankFourTensor ComputeMultiPlasticityStress::_my_flexural_rigidity_tensor
protected

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

Definition at line 171 of file ComputeMultiPlasticityStress.h.

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

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 168 of file ComputeMultiPlasticityStress.h.

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

MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n
protected

current value of transverse direction

Definition at line 126 of file ComputeMultiPlasticityStress.h.

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

RealVectorValue ComputeMultiPlasticityStress::_n_input
protected

the supplied transverse direction vector

Definition at line 90 of file ComputeMultiPlasticityStress.h.

Referenced by ComputeMultiPlasticityStress(), and initQpStatefulProperties().

const MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n_old
protected

old value of transverse direction

Definition at line 129 of file ComputeMultiPlasticityStress.h.

Referenced by postReturnMap().

bool ComputeMultiPlasticityStress::_n_supplied
protected

User supplied the transverse direction vector.

Definition at line 87 of file ComputeMultiPlasticityStress.h.

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

unsigned int MultiPlasticityRawComponentAssembler::_num_models
protectedinherited
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 56 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().

const InputParameters& MultiPlasticityRawComponentAssembler::_params
protectedinherited
bool ComputeMultiPlasticityStress::_perform_finite_strain_rotations
protected

whether to perform the rotations necessary in finite-strain simulations

Definition at line 96 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain
protected

plastic strain

Definition at line 99 of file ComputeMultiPlasticityStress.h.

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

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain_old
protected

Old value of plastic strain.

Definition at line 102 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

RealTensorValue ComputeMultiPlasticityStress::_rot
protected

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

Definition at line 93 of file ComputeMultiPlasticityStress.h.

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

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_rotation_increment
protected

Rotation increment (coming from ComputeIncrementalSmallStrain, for example)

Definition at line 138 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and postReturnMap().

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC
protectedinherited
bool& ComputeMultiPlasticityStress::_step_one
private

True if this is the first timestep (timestep < 2).

In the first timestep, an initial stress is needed to subdivide. This boolean variable eliminates the use of the _app.isRestarting() in this class. This boolean is delcared as a reference so that the variable is restartable data: if we restart, the code will not think it is the first timestep again.

Definition at line 591 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

const bool ComputeStressBase::_store_stress_old
protectedinherited

Parameter which decides whether to store old stress. This is required for HHT time integration and Rayleigh damping.

Definition at line 68 of file ComputeStressBase.h.

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_strain_increment
protected

strain increment (coming from ComputeIncrementalSmallStrain, for example)

Definition at line 132 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

MaterialProperty<RankTwoTensor>& ComputeStressBase::_stress
protectedinherited

Definition at line 53 of file ComputeStressBase.h.

Referenced by ComputeStressBase::addQpInitialStress(), ComputeMultipleInelasticCosseratStress::computeAdmissibleState(), ComputeMultipleInelasticStress::computeAdmissibleState(), ComputeStressBase::computeQpProperties(), ComputeStressEosBase::computeQpProperties(), ComputeBirchMurnaghanEquationOfStress::computeQpProperties(), ComputeLinearElasticStress::computeQpStress(), ComputeStrainIncrementBasedStress::computeQpStress(), ComputeLinearElasticPFFractureStress::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeIsotropicLinearElasticPFFractureStress::computeQpStress(), ComputeFiniteStrainElasticStressBirchMurnaghan::computeQpStress(), FiniteStrainPlasticMaterial::computeQpStress(), computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration(), ComputeSmearedCrackingStress::crackingStressRotation(), ComputeMultipleInelasticStress::finiteStrainRotation(), ComputeStressBase::initQpStatefulProperties(), FiniteStrainCrystalPlasticity::initQpStatefulProperties(), FiniteStrainUObasedCP::initQpStatefulProperties(), FiniteStrainHyperElasticViscoPlastic::initQpStatefulProperties(), postReturnMap(), FiniteStrainUObasedCP::postSolveQp(), FiniteStrainHyperElasticViscoPlastic::postSolveQp(), FiniteStrainCrystalPlasticity::postSolveQp(), ComputeMultipleInelasticStress::updateQpState(), ComputeMultipleInelasticStress::updateQpStateSingleModel(), and LinearIsoElasticPFDamage::updateVar().

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_stress_old
protected

Old value of stress.

Definition at line 141 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

std::vector<std::vector<unsigned int> > MultiPlasticityRawComponentAssembler::_surfaces_given_model
protectedinherited
Real MultiPlasticityLinearSystem::_svd_tol
protectedinherited

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

Definition at line 129 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

enum ComputeMultiPlasticityStress::TangentOperatorEnum ComputeMultiPlasticityStress::_tangent_operator_type
protected
const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_total_strain_old
protected

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

Definition at line 135 of file ComputeMultiPlasticityStress.h.

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

yield functions

Definition at line 111 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().


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