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:

Public Types
typedef DerivativeMaterialPropertyNameInterface::SymbolName	SymbolName

Public Member Functions
	ComputeMultiPlasticityStress (const InputParameters &parameters)
const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialProperty (const std::string &name)
const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialPropertyByName (const std::string &name)
void	validateDerivativeMaterialPropertyBase (const std::string &base)
const MaterialPropertyName	derivativePropertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const
const MaterialPropertyName	derivativePropertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const
const MaterialPropertyName	derivativePropertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const
const MaterialPropertyName	derivativePropertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)
GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)
const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")
void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)
void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")
void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)
void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")
void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")
void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")
const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const
const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const
const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const
const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const
const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const
const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const
const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const
const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const
void	outputAndCheckDebugParameters ()
	Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly.
void	checkDerivatives ()
	Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
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	checkSolution (const RankFourTensor &E_inv)
	Checks that Ax does equal b in the NR procedure.
UserObjectName	getUserObjectName (const std::string &param_name) const
const T &	getUserObject (const std::string &param_name, bool is_dependency=true) const
const T &	getUserObjectByName (const UserObjectName &object_name, bool is_dependency=true) const
const UserObjectBase &	getUserObjectBase (const std::string &param_name, bool is_dependency=true) const
const UserObjectBase &	getUserObjectBaseByName (const UserObjectName &object_name, bool is_dependency=true) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const

Static Public Member Functions
static InputParameters	validParams ()

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

Protected Member Functions
virtual void	computeQpStress ()
	Compute the stress and store it in the _stress material property for the current quadrature point.
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
virtual unsigned int	numberActive (const std::vector< bool > &active)
	counts the number of active constraints
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 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.
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 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.
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.
void	buildDumbOrder (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
	Builds the order which "dumb" activation will take.
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)
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	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	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 SolidMechanicsPlasticXXXX.returnMap functions.
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
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)
virtual void	computeQpProperties () override
virtual void	calculateConstraints (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
	The constraints.
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)
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.
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)
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	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)
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.
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	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].
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.
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.
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	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.
unsigned int	modelNumber (unsigned int surface)
	returns the model number, given the surface number
bool	anyActiveSurfaces (int model, const std::vector< bool > &active)
	returns true if any internal surfaces of the given model are active according to 'active'
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.
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.
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.
virtual void	addUserObjectDependencyHelper (const UserObjectBase &) const

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

Private Member Functions
RankTwoTensor	rot (const RankTwoTensor &tens)
bool	haveMaterialProperty (const std::string &prop_name)
std::vector< VariableName >	buildVariableVector (const VariableName &c1, const VariableName &c2, const VariableName &c3)
void	validateCouplingHelper (const MaterialPropertyName &base, const std::vector< VariableName > &c, const System &system, std::vector< VariableName > &missing)
bool	isNotObjectVariable (const VariableName &name)
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.
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)
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.
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	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.
bool	dof_included (unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
virtual int	singularValuesOfR (const std::vector< RankTwoTensor > &r, std::vector< Real > &s)
	Performs a singular-value decomposition of r and returns the singular values.
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 directions "r" If linear dependence is found, then deactivated_due_to_ld will contain 'true' entries where surfaces need to be deactivated_due_to_ld.
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.
const UserObjectBase &	getUserObjectFromFEProblem (const UserObjectName &object_name, const THREAD_ID tid=0) const
const T &	castUserObject (const UserObjectBase &uo_base, const std::string &param_name="") const
void	mooseObjectError (const std::string &param_name, std::stringstream &oss) const
const std::string &	userObjectType (const UserObjectBase &uo) const
const std::string &	userObjectName (const UserObjectBase &uo) const

Private Attributes
FEProblemBase &	_dmi_fe_problem
std::vector< unsigned int >	_model_given_surface
	given a surface number, this returns the model number
std::vector< unsigned int >	_model_surface_given_surface
	given a surface number, this returns the corresponding-model's internal surface number
const MooseObject &	_uoi_moose_object
const FEProblemBase &	_uoi_feproblem
const THREAD_ID	_uoi_tid

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

Member Enumeration Documentation

DeactivationSchemeEnum

enum ComputeMultiPlasticityStress::DeactivationSchemeEnum

protected

Enumerator
optimized
safe
dumb
optimized_to_safe
safe_to_dumb
optimized_to_safe_to_dumb
optimized_to_dumb

Definition at line 74 of file ComputeMultiPlasticityStress.h.

  {
    optimized,
    safe,
    dumb,
    optimized_to_safe,
    safe_to_dumb,
    optimized_to_safe_to_dumb,
    optimized_to_dumb
  } _deactivation_scheme;

quickStep_called_from_t

enum ComputeMultiPlasticityStress::quickStep_called_from_t

protected

The functions from which quickStep can be called.

Enumerator
computeQpStress_function
returnMap_function

Definition at line 440 of file ComputeMultiPlasticityStress.h.

  {
    computeQpStress_function,
    returnMap_function
  };

TangentOperatorEnum

enum ComputeMultiPlasticityStress::TangentOperatorEnum

protected

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

Enumerator
elastic
linear
nonlinear

Definition at line 49 of file ComputeMultiPlasticityStress.h.

  {
    elastic,
    linear,
    nonlinear
  } _tangent_operator_type;

Constructor & Destructor Documentation

ComputeMultiPlasticityStress()

ComputeMultiPlasticityStress::ComputeMultiPlasticityStress

(

const InputParameters &

parameters

)

Definition at line 97 of file ComputeMultiPlasticityStress.C.

  : ComputeStressBase(parameters),
    MultiPlasticityDebugger(this),
    _max_iter(getParam<unsigned int>("max_NR_iterations")),
    _min_stepsize(getParam<Real>("min_stepsize")),
    _max_stepsize_for_dumb(getParam<Real>("max_stepsize_for_dumb")),
    _ignore_failures(getParam<bool>("ignore_failures")),
 
    _tangent_operator_type((TangentOperatorEnum)(int)getParam<MooseEnum>("tangent_operator")),
 
    _epp_tol(getParam<Real>("ep_plastic_tolerance")),
 
    _dummy_pm(0),
 
    _cumulative_pm(0),
 
    _deactivation_scheme((DeactivationSchemeEnum)(int)getParam<MooseEnum>("deactivation_scheme")),
 
    _n_supplied(parameters.isParamValid("transverse_direction")),
    _n_input(_n_supplied ? getParam<RealVectorValue>("transverse_direction") : RealVectorValue()),
    _rot(RealTensorValue()),
 
    _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
 
    _elasticity_tensor_name(_base_name + "elasticity_tensor"),
    _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
    _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
    _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
    _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
    _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
    _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
    _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
                                                           // for visualisation i convert it to Real
    _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
                                                                            // boolean, but for
                                                                            // visualisation i
                                                                            // convert it to Real
    _ld_encountered(declareProperty<Real>(
        "plastic_linear_dependence_encountered")), // this is really a boolean, but for
                                                   // visualisation i convert it to Real
    _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
                                                                            // boolean, but for
                                                                            // visualisation i
                                                                            // convert it to Real
    _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
    _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
 
    _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
    _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
    _rotation_increment(
        getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
 
    _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
    _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
 
    // TODO: This design does NOT work. It makes these materials construction order dependent and it
    // disregards block restrictions.
    _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
              hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
    _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : nullptr),
    _elastic_flexural_rigidity_tensor(
        _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
                  : nullptr),
    _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : nullptr),
    _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress")
                                 : nullptr),
    _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
                                    : nullptr),
 
    _my_elasticity_tensor(RankFourTensor()),
    _my_strain_increment(RankTwoTensor()),
    _my_flexural_rigidity_tensor(RankFourTensor()),
    _my_curvature(RankTwoTensor())
{
  if (_epp_tol <= 0)
    mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
 
  if (_n_supplied)
  {
    // normalise the inputted transverse_direction
    if (_n_input.norm() == 0)
      mooseError(
          "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
    else
      _n_input /= _n_input.norm();
  }
 
  if (_num_surfaces == 1)
    _deactivation_scheme = safe;
}

Member Function Documentation

activeCombinationNumber()

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

protected

Definition at line 1572 of file ComputeMultiPlasticityStress.C.

{
  unsigned num = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (act[surface])
      num += (1 << surface); // (1 << x) = 2^x
 
  return num;
}

Referenced by returnMap().

activeModelSurfaces()

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

protectedinherited

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

Parameters

	model	the model number
	active	array with entries being 'true' if the surface is active
[out]	active_surfaces_of_model	the output

Definition at line 809 of file MultiPlasticityRawComponentAssembler.C.

{
  active_surfaces_of_model.resize(0);
  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
    if (active[_surfaces_given_model[model][model_surface]])
      active_surfaces_of_model.push_back(model_surface);
}

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().

activeSurfaces()

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

protectedinherited

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

Parameters

	model	the model number
	active	array with entries being 'true' if the surface is active
[out]	active_surfaces	the output

Definition at line 798 of file MultiPlasticityRawComponentAssembler.C.

{
  active_surfaces.resize(0);
  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
    if (active[_surfaces_given_model[model][model_surface]])
      active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
}

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

anyActiveSurfaces()

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

protectedinherited

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

Definition at line 789 of file MultiPlasticityRawComponentAssembler.C.

{
  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
    if (active[_surfaces_given_model[model][model_surface]])
      return true;
  return false;
}

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

applyKuhnTucker()

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

protectedvirtual

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

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

Parameters

f	values of the active yield functions
pm	values of all the plastic multipliers
active	the 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 1311 of file ComputeMultiPlasticityStress.C.

{
  bool turned_off = false;
  unsigned ind = 0;
 
  // turn off all active surfaces that have f<0 and pm!=0
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
        if (pm[surface] != 0)
        {
          turned_off = true;
          active[surface] = false;
        }
    }
    else if (pm[surface] != 0)
      mooseError("Crash due to plastic multiplier not being zero.  This occurred because of poor "
                 "coding!!");
  }
 
  // if didn't turn off anything yet, turn off surface with minimum pm
  if (!turned_off)
  {
    int surface_to_turn_off = -1;
    Real min_pm = 0;
    for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      if (pm[surface] < min_pm)
      {
        min_pm = pm[surface];
        surface_to_turn_off = surface;
      }
    if (surface_to_turn_off >= 0)
      active[surface_to_turn_off] = false;
  }
}

Referenced by returnMap().

buildActiveConstraints()

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

protectedvirtualinherited

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

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

Parameters

	f	yield functions (should be _num_surfaces of these)
	stress	stress tensor
	intnl	internal parameters
	Eijkl	elasticity tensor (stress = Eijkl*strain)
[out]	act	the set of active constraints (will be resized to _num_surfaces)

Definition at line 342 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(f.size() == _num_surfaces,
              "buildActiveConstraints called with f.size = " << f.size() << " while there are "
                                                             << _num_surfaces << " surfaces");
  mooseAssert(intnl.size() == _num_models,
              "buildActiveConstraints called with intnl.size = "
                  << intnl.size() << " while there are " << _num_models << " models");
 
  if (_specialIC == "rock")
    buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
  else if (_specialIC == "joint")
    buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
  else // no specialIC
  {
    act.resize(0);
    unsigned ind = 0;
    for (unsigned model = 0; model < _num_models; ++model)
    {
      std::vector<Real> model_f(0);
      for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
        model_f.push_back(f[ind++]);
      std::vector<bool> model_act;
      RankTwoTensor returned_stress;
      _f[model]->activeConstraints(
          model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
      for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
        act.push_back(model_act[model_surface]);
    }
  }
}

Referenced by returnMap().

buildActiveConstraintsJoint()

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

privateinherited

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

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

Parameters

	f	yield functions (should be _num_surfaces of these)
	stress	stress tensor
	intnl	internal parameters
	Eijkl	elasticity tensor (stress = Eijkl*strain)
[out]	act	the set of active constraints (will be resized to _num_surfaces)

Definition at line 379 of file MultiPlasticityRawComponentAssembler.C.

{
  act.assign(2, false);
 
  RankTwoTensor returned_stress;
  std::vector<bool> active_tensile;
  std::vector<bool> active_shear;
  std::vector<Real> f_single;
 
  // first try tensile alone
  f_single.assign(1, 0);
  f_single[0] = f[0];
  _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
  _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
  if (f_single[0] <= _f[1]->_f_tol)
  {
    act[0] = active_tensile[0];
    return;
  }
 
  // next try shear alone
  f_single.assign(1, 0);
  f_single[0] = f[1];
  _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_shear, returned_stress);
  _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
  if (f_single[0] <= _f[0]->_f_tol)
  {
    act[1] = active_shear[0];
    return;
  }
 
  // must be mixed
  act[0] = act[1] = true;
  return;
}

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints().

buildActiveConstraintsRock()

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

privateinherited

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

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

Parameters

	f	yield functions (should be _num_surfaces of these)
	stress	stress tensor
	intnl	internal parameters
	Eijkl	elasticity tensor (stress = Eijkl*strain)
[out]	act	the set of active constraints (will be resized to _num_surfaces)

Definition at line 420 of file MultiPlasticityRawComponentAssembler.C.

{
  act.assign(9, false);
 
  RankTwoTensor returned_stress;
  std::vector<bool> active_tensile;
  std::vector<bool> active_MC;
  std::vector<Real> f_single;
 
  // first try tensile alone
  f_single.assign(3, 0);
  f_single[0] = f[0];
  f_single[1] = f[1];
  f_single[2] = f[2];
  _f[0]->activeConstraints(f_single, stress, intnl[0], Eijkl, active_tensile, returned_stress);
  _f[1]->yieldFunctionV(returned_stress, intnl[1], f_single);
  if (f_single[0] <= _f[1]->_f_tol && f_single[1] <= _f[1]->_f_tol &&
      f_single[2] <= _f[1]->_f_tol && f_single[3] <= _f[1]->_f_tol &&
      f_single[4] <= _f[1]->_f_tol && f_single[5] <= _f[1]->_f_tol)
  {
    act[0] = active_tensile[0];
    act[1] = active_tensile[1];
    act[2] = active_tensile[2];
    return;
  }
 
  // next try MC alone
  f_single.assign(6, 0);
  f_single[0] = f[3];
  f_single[1] = f[4];
  f_single[2] = f[5];
  f_single[3] = f[6];
  f_single[4] = f[7];
  f_single[5] = f[8];
  _f[1]->activeConstraints(f_single, stress, intnl[1], Eijkl, active_MC, returned_stress);
  _f[0]->yieldFunctionV(returned_stress, intnl[0], f_single);
  if (f_single[0] <= _f[0]->_f_tol && f_single[1] <= _f[0]->_f_tol && f_single[2] <= _f[0]->_f_tol)
  {
    act[3] = active_MC[0];
    act[4] = active_MC[1];
    act[5] = active_MC[2];
    act[6] = active_MC[3];
    act[7] = active_MC[4];
    act[8] = active_MC[5];
    return;
  }
 
  // must be a mix.
  // The possibilities are enumerated below.
 
  // tensile=edge, MC=tip (two possibilities)
  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
      active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
  {
    act[1] = act[2] = act[6] = true;
    act[4] = true;
    return;
  }
  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
  {
    act[1] = act[2] = act[6] = true; // i don't think act[4] is necessary, is it?!
    return;
  }
 
  // tensile = edge, MC=edge (two possibilities)
  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
  {
    act[1] = act[2] = act[4] = act[6] = true;
    return;
  }
  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
  {
    act[1] = act[2] = act[4] = act[6] = true;
    return;
  }
 
  // tensile = edge, MC=face
  if (active_tensile[0] == false && active_tensile[1] == true && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
  {
    act[1] = act[2] = act[6] = true;
    return;
  }
 
  // tensile = face, MC=tip (two possibilities)
  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
      active_MC[0] == true && active_MC[1] == true && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
  {
    act[2] = act[6] = true;
    act[4] = true;
    act[8] = true;
    return;
  }
  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == true && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == true)
  {
    act[2] = act[6] = true;
    act[8] = true;
    return;
  }
 
  // tensile = face, MC=face
  if (active_tensile[0] == false && active_tensile[1] == false && active_tensile[2] == true &&
      active_MC[0] == false && active_MC[1] == false && active_MC[2] == false &&
      active_MC[3] == true && active_MC[4] == false && active_MC[5] == false)
  {
    act[1] = act[2] = act[6] = true;
    return;
  }
 
  // tensile = face, MC=edge (two possibilites).
  act[2] = true; // tensile face
  act[3] = active_MC[0];
  act[4] = active_MC[1];
  act[5] = active_MC[2];
  act[6] = active_MC[3];
  act[7] = active_MC[4];
  act[8] = active_MC[5];
  return;
}

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints().

buildDumbOrder()

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

protected

Builds the order which "dumb" activation will take.

Parameters

	stress	stress to evaluate yield functions and derivatives at
	intnl	internal parameters to evaluate yield functions and derivatives at
[out]	dumb_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).

Definition at line 1522 of file ComputeMultiPlasticityStress.C.

{
  if (dumb_order.size() != 0)
    return;
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  std::vector<Real> f;
  yieldFunction(stress, intnl, act, f);
  std::vector<RankTwoTensor> df_dstress;
  dyieldFunction_dstress(stress, intnl, act, df_dstress);
 
  typedef std::pair<Real, unsigned> pair_for_sorting;
  std::vector<pair_for_sorting> dist(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    dist[surface].first = f[surface] / df_dstress[surface].L2norm();
    dist[surface].second = surface;
  }
  std::sort(dist.begin(), dist.end()); // sorted in ascending order of f/df_dstress
 
  dumb_order.resize(_num_surfaces);
  for (unsigned i = 0; i < _num_surfaces; ++i)
    dumb_order[i] = dist[_num_surfaces - 1 - i].second;
  // now dumb_order[0] is the surface with the greatest f/df_dstress
}

Referenced by changeScheme(), and returnMap().

calculateConstraints()

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

protectedvirtualinherited

The constraints.

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

Parameters

	stress	The stress
	intnl_old	old values of the internal parameters
	intnl	internal parameters
	pm	Current value(s) of the plasticity multiplier(s) (consistency parameters)
	delta_dp	Change in plastic strain incurred so far during the return
[out]	f	Active yield function(s)
[out]	r	Active flow directions
[out]	epp	Plastic-strain increment constraint
[out]	ic	Active internal-parameter constraint
	active	The active constraints.

Definition at line 227 of file MultiPlasticityLinearSystem.C.

{
  // see comments at the start of .h file
 
  mooseAssert(intnl_old.size() == _num_models,
              "Size of intnl_old is " << intnl_old.size()
                                      << " which is incorrect in calculateConstraints");
  mooseAssert(intnl.size() == _num_models,
              "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
  mooseAssert(pm.size() == _num_surfaces,
              "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
  mooseAssert(active.size() == _num_surfaces,
              "Size of active is " << active.size()
                                   << " which is incorrect in calculateConstraints");
 
  // yield functions
  yieldFunction(stress, intnl, active, f);
 
  // flow directions and "epp"
  flowPotential(stress, intnl, active, r);
  epp = RankTwoTensor();
  unsigned ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface])
      epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
  epp -= delta_dp;
 
  // internal constraints
  std::vector<Real> h;
  hardPotential(stress, intnl, active, h);
  ic.resize(0);
  ind = 0;
  std::vector<unsigned int> active_surfaces;
  std::vector<unsigned int>::iterator active_surface;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeSurfaces(model, active, active_surfaces);
    if (active_surfaces.size() > 0)
    {
      // some surfaces are active in this model, so must form an internal constraint
      ic.push_back(intnl[model] - intnl_old[model]);
      for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
           ++active_surface)
        ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
                                                             // since it was constructed in the same
                                                             // manner
    }
  }
}

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

calculateJacobian()

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

protectedvirtualinherited

d(rhs)/d(dof)

Definition at line 366 of file MultiPlasticityLinearSystem.C.

{
  // see comments at the start of .h file
 
  mooseAssert(intnl.size() == _num_models,
              "Size of intnl is " << intnl.size() << " which is incorrect in calculateJacobian");
  mooseAssert(pm.size() == _num_surfaces,
              "Size of pm is " << pm.size() << " which is incorrect in calculateJacobian");
  mooseAssert(active.size() == _num_surfaces,
              "Size of active is " << active.size() << " which is incorrect in calculateJacobian");
  mooseAssert(deactivated_due_to_ld.size() == _num_surfaces,
              "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
                                                  << " which is incorrect in calculateJacobian");
 
  unsigned ind = 0;
  unsigned active_surface_ind = 0;
 
  std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
  unsigned num_active_surface = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active_surface[surface])
      num_active_surface++;
 
  std::vector<bool> active_model(
      _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
  for (unsigned model = 0; model < _num_models; ++model)
    active_model[model] = anyActiveSurfaces(model, active_surface);
 
  unsigned num_active_model = 0;
  for (unsigned model = 0; model < _num_models; ++model)
    if (active_model[model])
      num_active_model++;
 
  ind = 0;
  std::vector<unsigned int> active_model_index(_num_models);
  for (unsigned model = 0; model < _num_models; ++model)
    if (active_model[model])
      active_model_index[model] = ind++;
    else
      active_model_index[model] =
          _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
 
  std::vector<RankTwoTensor> df_dstress;
  dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
 
  std::vector<Real> df_dintnl;
  dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
 
  std::vector<RankTwoTensor> r;
  flowPotential(stress, intnl, active, r);
 
  std::vector<RankFourTensor> dr_dstress;
  dflowPotential_dstress(stress, intnl, active, dr_dstress);
 
  std::vector<RankTwoTensor> dr_dintnl;
  dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
 
  std::vector<Real> h;
  hardPotential(stress, intnl, active, h);
 
  std::vector<RankTwoTensor> dh_dstress;
  dhardPotential_dstress(stress, intnl, active, dh_dstress);
 
  std::vector<Real> dh_dintnl;
  dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
 
  // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
  RankFourTensor depp_dstress;
  ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface]) // includes deactivated_due_to_ld
      depp_dstress += pm[surface] * dr_dstress[ind++];
  depp_dstress += E_inv;
 
  // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
  std::vector<RankTwoTensor> depp_dpm;
  depp_dpm.resize(num_active_surface);
  ind = 0;
  active_surface_ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
                                   // dofs in the NR
        depp_dpm[active_surface_ind++] = r[ind];
      ind++;
    }
  }
 
  // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
  std::vector<RankTwoTensor> depp_dintnl;
  depp_dintnl.assign(num_active_model, RankTwoTensor());
  ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      unsigned int model_num = modelNumber(surface);
      if (active_model[model_num]) // only include models with surfaces which are still active after
                                   // deactivated_due_to_ld
        depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
      ind++;
    }
  }
 
  // df_dstress has been calculated above
  // df_dpm is always zero
  // df_dintnl has been calculated above, but only the active_surface+active_model stuff needs to be
  // included in Jacobian: see below
 
  std::vector<RankTwoTensor> dic_dstress;
  dic_dstress.assign(num_active_model, RankTwoTensor());
  ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      unsigned int model_num = modelNumber(surface);
      if (active_model[model_num]) // only include ic for models with active_surface (ie, if model
                                   // only contains deactivated_due_to_ld don't include it)
        dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
      ind++;
    }
  }
 
  std::vector<std::vector<Real>> dic_dpm;
  dic_dpm.resize(num_active_model);
  ind = 0;
  active_surface_ind = 0;
  for (unsigned model = 0; model < num_active_model; ++model)
    dic_dpm[model].assign(num_active_surface, 0);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
      {
        unsigned int model_num = modelNumber(surface);
        // if (active_model[model_num]) // do not need this check as if the surface has
        // active_surface, the model must be deemed active!
        dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
        active_surface_ind++;
      }
      ind++;
    }
  }
 
  std::vector<std::vector<Real>> dic_dintnl;
  dic_dintnl.resize(num_active_model);
  for (unsigned model = 0; model < num_active_model; ++model)
  {
    dic_dintnl[model].assign(num_active_model, 0);
    dic_dintnl[model][model] = 1; // deriv wrt internal parameter
  }
  ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      unsigned int model_num = modelNumber(surface);
      if (active_model[model_num]) // only the models that contain surfaces that are still active
                                   // after deactivation_due_to_ld
        dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
            pm[surface] * dh_dintnl[ind];
      ind++;
    }
  }
 
  unsigned int dim = 3;
  unsigned int system_size =
      6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
  jac.resize(system_size);
  for (unsigned i = 0; i < system_size; ++i)
    jac[i].assign(system_size, 0);
 
  unsigned int row_num = 0;
  unsigned int col_num = 0;
  for (unsigned i = 0; i < dim; ++i)
    for (unsigned j = 0; j <= i; ++j)
    {
      for (unsigned k = 0; k < dim; ++k)
        for (unsigned l = 0; l <= k; ++l)
          jac[col_num][row_num++] =
              depp_dstress(i, j, k, l) +
              (k != l ? depp_dstress(i, j, l, k)
                      : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
      for (unsigned surface = 0; surface < num_active_surface; ++surface)
        jac[col_num][row_num++] = depp_dpm[surface](i, j);
      for (unsigned a = 0; a < num_active_model; ++a)
        jac[col_num][row_num++] = depp_dintnl[a](i, j);
      row_num = 0;
      col_num++;
    }
 
  ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active_surface[surface])
    {
      for (unsigned k = 0; k < dim; ++k)
        for (unsigned l = 0; l <= k; ++l)
          jac[col_num][row_num++] =
              df_dstress[ind](k, l) +
              (k != l ? df_dstress[ind](l, k)
                      : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
      for (unsigned beta = 0; beta < num_active_surface; ++beta)
        jac[col_num][row_num++] = 0; // df_dpm
      for (unsigned model = 0; model < _num_models; ++model)
        if (active_model[model]) // only use df_dintnl for models in active_model
        {
          if (modelNumber(surface) == model)
            jac[col_num][row_num++] = df_dintnl[ind];
          else
            jac[col_num][row_num++] = 0;
        }
      ind++;
      row_num = 0;
      col_num++;
    }
 
  for (unsigned a = 0; a < num_active_model; ++a)
  {
    for (unsigned k = 0; k < dim; ++k)
      for (unsigned l = 0; l <= k; ++l)
        jac[col_num][row_num++] =
            dic_dstress[a](k, l) +
            (k != l ? dic_dstress[a](l, k)
                    : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
    for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
      jac[col_num][row_num++] = dic_dpm[a][alpha];
    for (unsigned b = 0; b < num_active_model; ++b)
      jac[col_num][row_num++] = dic_dintnl[a][b];
    row_num = 0;
    col_num++;
  }
 
  mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
}

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

calculateRHS()

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

protectedvirtualinherited

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

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

Parameters

	stress	The stress
	intnl_old	old values of the internal parameters
	intnl	internal parameters
	pm	Current value(s) of the plasticity multiplier(s) (consistency parameters)
	delta_dp	Change in plastic strain incurred so far during the return
[out]	rhs	the rhs
	active	The active constraints.
	eliminate_ld	Check 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_ld	constraints deactivated due to linear-dependence of flow directions

Definition at line 287 of file MultiPlasticityLinearSystem.C.

{
  // see comments at the start of .h file
 
  mooseAssert(intnl_old.size() == _num_models,
              "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
  mooseAssert(intnl.size() == _num_models,
              "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
  mooseAssert(pm.size() == _num_surfaces,
              "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
  mooseAssert(active.size() == _num_surfaces,
              "Size of active is " << active.size() << " which is incorrect in calculateRHS");
 
  std::vector<Real> f;  // the yield functions
  RankTwoTensor epp;    // the plastic-strain constraint ("direction constraint")
  std::vector<Real> ic; // the "internal constraints"
 
  std::vector<RankTwoTensor> r;
  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
 
  if (eliminate_ld)
    eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
  else
    deactivated_due_to_ld.assign(_num_surfaces, false);
 
  std::vector<bool> active_not_deact(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
 
  unsigned num_active_f = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active_not_deact[surface])
      num_active_f++;
 
  unsigned num_active_ic = 0;
  for (unsigned model = 0; model < _num_models; ++model)
    if (anyActiveSurfaces(model, active_not_deact))
      num_active_ic++;
 
  unsigned int dim = 3;
  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
                                                               // num_active_f comes from "f",
                                                               // num_active_f comes from "ic"
 
  rhs.resize(system_size);
 
  unsigned ind = 0;
  for (unsigned i = 0; i < dim; ++i)
    for (unsigned j = 0; j <= i; ++j)
      rhs[ind++] = -epp(i, j);
  unsigned active_surface = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface])
    {
      if (!deactivated_due_to_ld[surface])
        rhs[ind++] = -f[active_surface];
      active_surface++;
    }
  unsigned active_model = 0;
  for (unsigned model = 0; model < _num_models; ++model)
    if (anyActiveSurfaces(model, active))
    {
      if (anyActiveSurfaces(model, active_not_deact))
        rhs[ind++] = -ic[active_model];
      active_model++;
    }
 
  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
}

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

canAddConstraints()

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

protected

Definition at line 1013 of file ComputeMultiPlasticityStress.C.

{
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
      return true;
  return false;
}

Referenced by returnMap().

canChangeScheme()

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

protected

Definition at line 1023 of file ComputeMultiPlasticityStress.C.

{
  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
    return true;
 
  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
    return true;
 
  if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
      can_revert_to_dumb)
    return true;
 
  if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
      can_revert_to_dumb)
    return true;
 
  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
      can_revert_to_dumb)
    return true;
 
  return false;
}

Referenced by returnMap().

canIncrementDumb()

bool ComputeMultiPlasticityStress::canIncrementDumb ( int  dumb_iteration )

protected

Definition at line 1565 of file ComputeMultiPlasticityStress.C.

{
  // (1 << _num_surfaces) = 2^_num_surfaces
  return ((dumb_iteration + 1) < (1 << _num_surfaces));
}

Referenced by returnMap().

changeScheme()

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

protected

Definition at line 1048 of file ComputeMultiPlasticityStress.C.

{
  if (current_deactivation_scheme == optimized &&
      (_deactivation_scheme == optimized_to_safe ||
       _deactivation_scheme == optimized_to_safe_to_dumb))
  {
    current_deactivation_scheme = safe;
    for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      act[surface] = initial_act[surface];
  }
  else if ((current_deactivation_scheme == safe &&
            (_deactivation_scheme == optimized_to_safe_to_dumb ||
             _deactivation_scheme == safe_to_dumb) &&
            can_revert_to_dumb) ||
           (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
            can_revert_to_dumb))
  {
    current_deactivation_scheme = dumb;
    dumb_iteration = 0;
    buildDumbOrder(initial_stress, intnl_old, dumb_order);
    incrementDumb(dumb_iteration, dumb_order, act);
  }
}

Referenced by returnMap().

checkAdmissible()

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

protectedvirtual

Checks whether the yield functions are in the admissible region.

Parameters

	stress	stress
	intnl	internal parameters
[out]	all_f	the values of all the yield functions

Returns

return false if any yield functions exceed their tolerance

Definition at line 1269 of file ComputeMultiPlasticityStress.C.

{
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  yieldFunction(stress, intnl, act, all_f);
 
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
      return false;
 
  return true;
}

Referenced by returnMap().

checkDerivatives()

void MultiPlasticityDebugger::checkDerivatives ( )

inherited

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

Definition at line 85 of file MultiPlasticityDebugger.C.

{
  Moose::err
      << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
  outputAndCheckDebugParameters();
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  Moose::err << "\ndyieldFunction_dstress.  Relative L2 norms.\n";
  std::vector<RankTwoTensor> df_dstress;
  std::vector<RankTwoTensor> fddf_dstress;
  dyieldFunction_dstress(_fspb_debug_stress, _fspb_debug_intnl, act, df_dstress);
  fddyieldFunction_dstress(_fspb_debug_stress, _fspb_debug_intnl, fddf_dstress);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    Moose::err << "surface = " << surface << " Relative L2norm = "
               << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
                      (df_dstress[surface] + fddf_dstress[surface]).L2norm()
               << "\n";
    Moose::err << "Coded:\n";
    df_dstress[surface].print();
    Moose::err << "Finite difference:\n";
    fddf_dstress[surface].print();
  }
 
  Moose::err << "\ndyieldFunction_dintnl.\n";
  std::vector<Real> df_dintnl;
  dyieldFunction_dintnl(_fspb_debug_stress, _fspb_debug_intnl, act, df_dintnl);
  Moose::err << "Coded:\n";
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    Moose::err << df_dintnl[surface] << " ";
  Moose::err << "\n";
  std::vector<Real> fddf_dintnl;
  fddyieldFunction_dintnl(_fspb_debug_stress, _fspb_debug_intnl, fddf_dintnl);
  Moose::err << "Finite difference:\n";
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    Moose::err << fddf_dintnl[surface] << " ";
  Moose::err << "\n";
 
  Moose::err << "\ndflowPotential_dstress.  Relative L2 norms.\n";
  std::vector<RankFourTensor> dr_dstress;
  std::vector<RankFourTensor> fddr_dstress;
  dflowPotential_dstress(_fspb_debug_stress, _fspb_debug_intnl, act, dr_dstress);
  fddflowPotential_dstress(_fspb_debug_stress, _fspb_debug_intnl, fddr_dstress);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    Moose::err << "surface = " << surface << " Relative L2norm = "
               << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
                      (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
               << "\n";
    Moose::err << "Coded:\n";
    dr_dstress[surface].print();
    Moose::err << "Finite difference:\n";
    fddr_dstress[surface].print();
  }
 
  Moose::err << "\ndflowPotential_dintnl.  Relative L2 norms.\n";
  std::vector<RankTwoTensor> dr_dintnl;
  std::vector<RankTwoTensor> fddr_dintnl;
  dflowPotential_dintnl(_fspb_debug_stress, _fspb_debug_intnl, act, dr_dintnl);
  fddflowPotential_dintnl(_fspb_debug_stress, _fspb_debug_intnl, fddr_dintnl);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    Moose::err << "surface = " << surface << " Relative L2norm = "
               << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
                      (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
               << "\n";
    Moose::err << "Coded:\n";
    dr_dintnl[surface].print();
    Moose::err << "Finite difference:\n";
    fddr_dintnl[surface].print();
  }
 
  Moose::err << std::flush;
}

Referenced by initQpStatefulProperties(), and plasticStep().

checkJacobian()

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

inherited

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

Definition at line 163 of file MultiPlasticityDebugger.C.

{
  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
  outputAndCheckDebugParameters();
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
  std::vector<bool> deactivated_due_to_ld;
  deactivated_due_to_ld.assign(_num_surfaces, false);
 
  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
 
  std::vector<std::vector<Real>> jac;
  calculateJacobian(_fspb_debug_stress,
                    _fspb_debug_intnl,
                    _fspb_debug_pm,
                    E_inv,
                    act,
                    deactivated_due_to_ld,
                    jac);
 
  std::vector<std::vector<Real>> fdjac;
  fdJacobian(_fspb_debug_stress,
             intnl_old,
             _fspb_debug_intnl,
             _fspb_debug_pm,
             delta_dp,
             E_inv,
             false,
             fdjac);
 
  Real L2_numer = 0;
  Real L2_denom = 0;
  for (unsigned row = 0; row < jac.size(); ++row)
    for (unsigned col = 0; col < jac.size(); ++col)
    {
      L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
      L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
    }
  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
 
  Moose::err << "\nHand-coded Jacobian:\n";
  for (unsigned row = 0; row < jac.size(); ++row)
  {
    for (unsigned col = 0; col < jac.size(); ++col)
      Moose::err << jac[row][col] << " ";
    Moose::err << "\n";
  }
 
  Moose::err << "Finite difference Jacobian:\n";
  for (unsigned row = 0; row < fdjac.size(); ++row)
  {
    for (unsigned col = 0; col < fdjac.size(); ++col)
      Moose::err << fdjac[row][col] << " ";
    Moose::err << "\n";
  }
 
  Moose::err << std::flush;
}

Referenced by computeQpStress(), and plasticStep().

checkKuhnTucker()

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

protectedvirtual

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

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

Parameters

f	values of the active yield functions
pm	values of all the plastic multipliers
active	the 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 1286 of file ComputeMultiPlasticityStress.C.

{
  unsigned ind = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (active[surface])
    {
      if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
        if (pm[surface] != 0)
          return false;
    }
    else if (pm[surface] != 0)
      mooseError("Crash due to plastic multiplier not being zero.  This occurred because of poor "
                 "coding!!");
  }
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (pm[surface] < 0)
      return false;
 
  return true;
}

Referenced by returnMap().

checkSolution()

void MultiPlasticityDebugger::checkSolution ( const RankFourTensor &  E_inv )

inherited

Checks that Ax does equal b in the NR procedure.

Definition at line 371 of file MultiPlasticityDebugger.C.

{
  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
  outputAndCheckDebugParameters();
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
  std::vector<bool> deactivated_due_to_ld;
  deactivated_due_to_ld.assign(_num_surfaces, false);
 
  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
 
  std::vector<Real> orig_rhs;
  calculateRHS(_fspb_debug_stress,
               _fspb_debug_intnl,
               _fspb_debug_intnl,
               _fspb_debug_pm,
               delta_dp,
               orig_rhs,
               act,
               true,
               deactivated_due_to_ld);
 
  Moose::err << "\nb = ";
  for (unsigned i = 0; i < orig_rhs.size(); ++i)
    Moose::err << orig_rhs[i] << " ";
  Moose::err << "\n\n";
 
  std::vector<std::vector<Real>> jac_coded;
  calculateJacobian(_fspb_debug_stress,
                    _fspb_debug_intnl,
                    _fspb_debug_pm,
                    E_inv,
                    act,
                    deactivated_due_to_ld,
                    jac_coded);
 
  Moose::err
      << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
  Moose::err
      << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
                "linear dependence.\n";
  Moose::err << "Hand-coded Jacobian:\n";
  for (unsigned row = 0; row < jac_coded.size(); ++row)
  {
    for (unsigned col = 0; col < jac_coded.size(); ++col)
      Moose::err << jac_coded[row][col] << " ";
    Moose::err << "\n";
  }
 
  deactivated_due_to_ld.assign(_num_surfaces,
                               false); // this potentially gets changed by nrStep, below
  RankTwoTensor dstress;
  std::vector<Real> dpm;
  std::vector<Real> dintnl;
  nrStep(_fspb_debug_stress,
         _fspb_debug_intnl,
         _fspb_debug_intnl,
         _fspb_debug_pm,
         E_inv,
         delta_dp,
         dstress,
         dpm,
         dintnl,
         act,
         deactivated_due_to_ld);
 
  std::vector<bool> active_not_deact(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_not_deact[surface] = !deactivated_due_to_ld[surface];
 
  std::vector<Real> x;
  x.assign(orig_rhs.size(), 0);
  unsigned ind = 0;
  for (unsigned i = 0; i < 3; ++i)
    for (unsigned j = 0; j <= i; ++j)
      x[ind++] = dstress(i, j);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active_not_deact[surface])
      x[ind++] = dpm[surface];
  for (unsigned model = 0; model < _num_models; ++model)
    if (anyActiveSurfaces(model, active_not_deact))
      x[ind++] = dintnl[model];
 
  mooseAssert(ind == orig_rhs.size(),
              "Incorrect extracting of changes from NR solution in the "
              "finite-difference checking of nrStep");
 
  Moose::err << "\nThis yields x =";
  for (unsigned i = 0; i < orig_rhs.size(); ++i)
    Moose::err << x[i] << " ";
  Moose::err << "\n";
 
  std::vector<std::vector<Real>> jac_fd;
  fdJacobian(_fspb_debug_stress,
             _fspb_debug_intnl,
             _fspb_debug_intnl,
             _fspb_debug_pm,
             delta_dp,
             E_inv,
             true,
             jac_fd);
 
  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
  Moose::err << "in order to check that the solution is correct\n";
  Moose::err << "Finite-difference Jacobian:\n";
  for (unsigned row = 0; row < jac_fd.size(); ++row)
  {
    for (unsigned col = 0; col < jac_fd.size(); ++col)
      Moose::err << jac_fd[row][col] << " ";
    Moose::err << "\n";
  }
 
  Real L2_numer = 0;
  Real L2_denom = 0;
  for (unsigned row = 0; row < jac_coded.size(); ++row)
    for (unsigned col = 0; col < jac_coded.size(); ++col)
    {
      L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
      L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
    }
  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
             << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
 
  std::vector<Real> fd_times_x;
  fd_times_x.assign(orig_rhs.size(), 0);
  for (unsigned row = 0; row < orig_rhs.size(); ++row)
    for (unsigned col = 0; col < orig_rhs.size(); ++col)
      fd_times_x[row] += jac_fd[row][col] * x[col];
 
  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
  for (unsigned i = 0; i < orig_rhs.size(); ++i)
    Moose::err << fd_times_x[i] << " ";
  Moose::err << "\n";
  Moose::err << "Recall that b = \n";
  for (unsigned i = 0; i < orig_rhs.size(); ++i)
    Moose::err << orig_rhs[i] << " ";
  Moose::err << "\n";
 
  L2_numer = 0;
  L2_denom = 0;
  for (unsigned i = 0; i < orig_rhs.size(); ++i)
  {
    L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
    L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
  }
  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5
             << std::endl;
}

Referenced by computeQpStress(), and plasticStep().

computeQpProperties()

void ComputeGeneralStressBase::computeQpProperties ( )

overrideprotectedvirtualinherited

Definition at line 44 of file ComputeGeneralStressBase.C.

{
  computeQpStress();
 
  // Add in extra stress
  _stress[_qp] += _extra_stress[_qp];
}

computeQpStress()

void ComputeMultiPlasticityStress::computeQpStress ( )

protectedvirtual

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

Implements ComputeGeneralStressBase.

Definition at line 219 of file ComputeMultiPlasticityStress.C.

{
  // the following "_my" variables can get rotated by preReturnMap and postReturnMap
  _my_elasticity_tensor = _elasticity_tensor[_qp];
  _my_strain_increment = _strain_increment[_qp];
  if (_cosserat)
  {
    _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
    _my_curvature = (*_curvature)[_qp];
  }
 
  if (_fspb_debug == "jacobian_and_linear_system")
  {
    // cannot do this at initQpStatefulProperties level since E_ijkl is not defined
    checkJacobian(_elasticity_tensor[_qp].invSymm(), _intnl_old[_qp]);
    checkSolution(_elasticity_tensor[_qp].invSymm());
    mooseError("Finite-differencing completed.  Exiting with no error");
  }
 
  preReturnMap(); // do rotations to new frame if necessary
 
  unsigned int number_iterations;
  bool linesearch_needed = false;
  bool ld_encountered = false;
  bool constraints_added = false;
 
  _cumulative_pm.assign(_num_surfaces, 0);
  // try a "quick" return first - this can be purely elastic, or a customised plastic return defined
  // by a SolidMechanicsPlasticXXXX UserObject
  const bool found_solution = quickStep(rot(_stress_old[_qp]),
                                        _stress[_qp],
                                        _intnl_old[_qp],
                                        _intnl[_qp],
                                        _dummy_pm,
                                        _cumulative_pm,
                                        rot(_plastic_strain_old[_qp]),
                                        _plastic_strain[_qp],
                                        _my_elasticity_tensor,
                                        _my_strain_increment,
                                        _yf[_qp],
                                        number_iterations,
                                        _Jacobian_mult[_qp],
                                        computeQpStress_function,
                                        true);
 
  // if not purely elastic or the customised stuff failed, do some plastic return
  if (!found_solution)
    plasticStep(rot(_stress_old[_qp]),
                _stress[_qp],
                _intnl_old[_qp],
                _intnl[_qp],
                rot(_plastic_strain_old[_qp]),
                _plastic_strain[_qp],
                _my_elasticity_tensor,
                _my_strain_increment,
                _yf[_qp],
                number_iterations,
                linesearch_needed,
                ld_encountered,
                constraints_added,
                _Jacobian_mult[_qp]);
 
  if (_cosserat)
  {
    (*_couple_stress)[_qp] = (*_elastic_flexural_rigidity_tensor)[_qp] * _my_curvature;
    (*_Jacobian_mult_couple)[_qp] = _my_flexural_rigidity_tensor;
  }
 
  postReturnMap(); // rotate back from new frame if necessary
 
  _iter[_qp] = 1.0 * number_iterations;
  _linesearch_needed[_qp] = linesearch_needed;
  _ld_encountered[_qp] = ld_encountered;
  _constraints_added[_qp] = constraints_added;
 
  // Update measures of strain
  _elastic_strain[_qp] = _elastic_strain_old[_qp] + _my_strain_increment -
                         (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
 
  // Rotate the tensors to the current configuration
  if (_perform_finite_strain_rotations)
  {
    _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
    _elastic_strain[_qp] =
        _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
    _plastic_strain[_qp] =
        _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
  }
}

consistentTangentOperator()

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

protected

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

The computations performed depend upon _tangent_operator_type

Parameters

stress	The value of stress after the return map algorithm has converged
intnl	The internal parameters after the return map has converged
E_ijkl	The elasticity tensor (in the case of no plasticity this is the jacobian)
pm_this_step	The 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_pm	The 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 1583 of file ComputeMultiPlasticityStress.C.

{
 
  if (_tangent_operator_type == elastic)
    return E_ijkl;
 
  // Typically act_at_some_step = act, but it is possible
  // that when subdividing a strain increment, a surface
  // is only active for one sub-step
  std::vector<bool> act_at_some_step(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    act_at_some_step[surface] = (cumulative_pm[surface] > 0);
 
  // "act" might contain surfaces that are linearly dependent
  // with others.  Only the plastic multipliers that are > 0
  // for this strain increment need to be varied to find
  // the consistent tangent operator
  std::vector<bool> act_vary(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    act_vary[surface] = (pm_this_step[surface] > 0);
 
  std::vector<RankTwoTensor> df_dstress;
  dyieldFunction_dstress(stress, intnl, act_vary, df_dstress);
  std::vector<Real> df_dintnl;
  dyieldFunction_dintnl(stress, intnl, act_vary, df_dintnl);
  std::vector<RankTwoTensor> r;
  flowPotential(stress, intnl, act_vary, r);
  std::vector<RankFourTensor> dr_dstress_at_some_step;
  dflowPotential_dstress(stress, intnl, act_at_some_step, dr_dstress_at_some_step);
  std::vector<RankTwoTensor> dr_dintnl_at_some_step;
  dflowPotential_dintnl(stress, intnl, act_at_some_step, dr_dintnl_at_some_step);
  std::vector<Real> h;
  hardPotential(stress, intnl, act_vary, h);
 
  unsigned ind1;
  unsigned ind2;
 
  // r_minus_stuff[alpha] = r[alpha] -
  // pm_cumulatve[gamma]*dr[gamma]_dintnl[a]_at_some_step*h[a][alpha], with alpha only being in
  // act_vary, but gamma being act_at_some_step
  std::vector<RankTwoTensor> r_minus_stuff;
  ind1 = 0;
  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
    if (act_vary[surface1])
    {
      r_minus_stuff.push_back(r[ind1]);
      ind2 = 0;
      for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
        if (act_at_some_step[surface2])
        {
          if (modelNumber(surface1) == modelNumber(surface2))
          {
            r_minus_stuff.back() -=
                cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
          }
          ind2++;
        }
      ind1++;
    }
 
  unsigned int num_currently_active = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (act_vary[surface])
      num_currently_active += 1;
 
  // zzz is a matrix in the form that can be easily
  // inverted by MatrixTools::inverse
  // Eg for num_currently_active = 3
  // (zzz[0] zzz[1] zzz[2])
  // (zzz[3] zzz[4] zzz[5])
  // (zzz[6] zzz[7] zzz[8])
  std::vector<PetscScalar> zzz;
  zzz.assign(num_currently_active * num_currently_active, 0.0);
 
  ind1 = 0;
  RankTwoTensor r2;
  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
    if (act_vary[surface1])
    {
      ind2 = 0;
      for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
        if (act_vary[surface2])
        {
          r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
          zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
          if (modelNumber(surface1) == modelNumber(surface2))
            zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
          ind2++;
        }
      ind1++;
    }
 
  if (num_currently_active > 0)
  {
    // invert zzz, in place.  if num_currently_active = 0 then zzz is not needed.
    try
    {
      MatrixTools::inverse(zzz, num_currently_active);
    }
    catch (const MooseException & e)
    {
      // in the very rare case of zzz being singular, just return the "elastic" tangent operator
      return E_ijkl;
    }
  }
 
  RankFourTensor strain_coeff = E_ijkl;
  ind1 = 0;
  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
    if (act_vary[surface1])
    {
      RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
      ind2 = 0;
      for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
        if (act_vary[surface2])
        {
          RankTwoTensor part2 = E_ijkl * df_dstress[ind2];
          for (unsigned i = 0; i < 3; i++)
            for (unsigned j = 0; j < 3; j++)
              for (unsigned k = 0; k < 3; k++)
                for (unsigned l = 0; l < 3; l++)
                  strain_coeff(i, j, k, l) -=
                      part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
          ind2++;
        }
      ind1++;
    }
 
  if (_tangent_operator_type == linear)
    return strain_coeff;
 
  RankFourTensor stress_coeff(RankFourTensor::initIdentitySymmetricFour);
 
  RankFourTensor part3;
  ind1 = 0;
  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
    if (act_at_some_step[surface1])
    {
      part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
      ind1++;
    }
 
  stress_coeff += part3;
 
  part3 = part3.transposeMajor(); // this is because below i want df_dstress[ind2]*part3, and this
                                  // equals (part3.transposeMajor())*df_dstress[ind2]
 
  ind1 = 0;
  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
    if (act_vary[surface1])
    {
      RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
      ind2 = 0;
      for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
        if (act_vary[surface2])
        {
          RankTwoTensor part2 = part3 * df_dstress[ind2];
          for (unsigned i = 0; i < 3; i++)
            for (unsigned j = 0; j < 3; j++)
              for (unsigned k = 0; k < 3; k++)
                for (unsigned l = 0; l < 3; l++)
                  stress_coeff(i, j, k, l) -=
                      part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
          ind2++;
        }
      ind1++;
    }
 
  // need to find the inverse of stress_coeff, but remember
  // stress_coeff does not have the symmetries commonly found
  // in tensor mechanics:
  // stress_coeff(i, j, k, l) = stress_coeff(j, i, k, l) = stress_coeff(i, j, l, k) !=
  // stress_coeff(k, l, i, j)
  // (note the final not-equals).  We want s_inv, such that
  // s_inv(i, j, m, n)*stress_coeff(m, n, k, l) = (de_ik*de_jl + de_il*de_jk)/2
  // where de_ij = 1 if i=j, and 0 otherwise.
  RankFourTensor s_inv;
  try
  {
    s_inv = stress_coeff.invSymm();
  }
  catch (const MooseException & e)
  {
    return strain_coeff; // when stress_coeff is singular (perhaps for incompressible plasticity?)
                         // return the "linear" tangent operator
  }
 
  return s_inv * strain_coeff;
}

Referenced by quickStep(), and returnMap().

dflowPotential_dintnl()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 233 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  dr_dintnl.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<RankTwoTensor> model_dr_dintnl;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
    }
  }
}

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

dflowPotential_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 205 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  dr_dstress.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<RankFourTensor> model_dr_dstress;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        dr_dstress.push_back(model_dr_dstress[*active_surface]);
    }
  }
}

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

dhardPotential_dintnl()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the hardening potentials
	intnl	vector of internal parameters
	active	set of active constraints - only the active derivatives are put into "dh_dintnl"
[out]	dh_dintnl	the 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 315 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  dh_dintnl.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<Real> model_dh_dintnl;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
    }
  }
}

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

dhardPotential_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 287 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  dh_dstress.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<RankTwoTensor> model_dh_dstress;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        dh_dstress.push_back(model_dh_dstress[*active_surface]);
    }
  }
}

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

dof_included()

bool MultiPlasticityDebugger::dof_included ( unsigned int  dof, const std::vector< bool > &  deactivated_due_to_ld  )

privateinherited

Definition at line 353 of file MultiPlasticityDebugger.C.

{
  if (dof < unsigned(6))
    // these are the stress components
    return true;
  unsigned eff_dof = dof - 6;
  if (eff_dof < _num_surfaces)
    // these are the plastic multipliers, pm
    return !deactivated_due_to_ld[eff_dof];
  eff_dof -= _num_surfaces; // now we know the dof is an intnl parameter
  std::vector<bool> active_surface(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_surface[surface] = !deactivated_due_to_ld[surface];
  return anyActiveSurfaces(eff_dof, active_surface);
}

Referenced by MultiPlasticityDebugger::fdJacobian().

dyieldFunction_dintnl()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the yield function
	intnl	vector of internal parameters
	active	set of active constraints - only the active derivatives are put into "df_dintnl"
[out]	df_dintnl	the derivatives. df_dstress[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 151 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  df_dintnl.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<Real> model_df_dintnl;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        df_dintnl.push_back(model_df_dintnl[*active_surface]);
    }
  }
}

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

dyieldFunction_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 123 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  df_dstress.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<RankTwoTensor> model_df_dstress;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        df_dstress.push_back(model_df_dstress[*active_surface]);
    }
  }
}

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

eliminateLinearDependence()

void MultiPlasticityLinearSystem::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  )

privatevirtualinherited

Performs a number of singular-value decompositions to check for linear-dependence of the active directions "r" If linear dependence is found, then deactivated_due_to_ld will contain 'true' entries where surfaces need to be deactivated_due_to_ld.

Parameters

	stress	the current stress
	intnl	the current values of internal parameters
	f	Active yield function values
	r	the flow directions that for those yield functions that are active upon entry to this function
	active	true if active
[out]	deactivated_due_to_ld	Yield functions deactivated due to linearly-dependent flow directions

Definition at line 106 of file MultiPlasticityLinearSystem.C.

{
  deactivated_due_to_ld.resize(_num_surfaces, false);
 
  unsigned num_active = r.size();
 
  if (num_active <= 1)
    return;
 
  std::vector<double> s;
  int info = singularValuesOfR(r, s);
  if (info != 0)
    mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
               "returned with error code ",
               info);
 
  // num_lin_dep are the number of linearly dependent
  // "r vectors", if num_active <= 6
  unsigned int num_lin_dep = 0;
 
  unsigned i = s.size();
  while (i-- > 0)
    if (s[i] < _svd_tol * s[0])
      num_lin_dep++;
    else
      break;
 
  if (num_lin_dep == 0 && num_active <= 6)
    return;
 
  // From here on, some flow directions are linearly dependent
 
  // Find the signed "distance" of the current (stress, internal) configuration
  // from the yield surfaces.  This distance will not be precise, but
  // i want to preferentially deactivate yield surfaces that are close
  // to the current stress point.
  std::vector<RankTwoTensor> df_dstress;
  dyieldFunction_dstress(stress, intnl, active, df_dstress);
 
  typedef std::pair<Real, unsigned> pair_for_sorting;
  std::vector<pair_for_sorting> dist(num_active);
  for (unsigned i = 0; i < num_active; ++i)
  {
    dist[i].first = f[i] / df_dstress[i].L2norm();
    dist[i].second = i;
  }
  std::sort(dist.begin(), dist.end()); // sorted in ascending order
 
  // There is a potential problem when we have equal f[i], for it can give oscillations
  bool equals_detected = false;
  for (unsigned i = 0; i < num_active - 1; ++i)
    if (std::abs(dist[i].first - dist[i + 1].first) < _min_f_tol)
    {
      equals_detected = true;
      dist[i].first += _min_f_tol * (MooseRandom::rand() - 0.5);
    }
  if (equals_detected)
    std::sort(dist.begin(), dist.end()); // sorted in ascending order
 
  std::vector<bool> scheduled_for_deactivation;
  scheduled_for_deactivation.assign(num_active, false);
 
  // In the following loop we go through all the flow directions, from
  // the one with the largest dist, to the one with the smallest dist,
  // adding them one-by-one into r_tmp.  Upon each addition we check
  // for linear-dependence.  if LD is found, we schedule the most
  // recently added flow direction for deactivation, and pop it
  // back off r_tmp
  unsigned current_yf;
  current_yf = dist[num_active - 1].second;
  // the one with largest dist
  std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
 
  unsigned num_kept_active = 1;
  for (unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
  {
    current_yf = dist[num_active - yf_to_try].second;
    if (num_active == 2) // shortcut to we don't have to singularValuesOfR
      scheduled_for_deactivation[current_yf] = true;
    else if (num_kept_active >= 6) // shortcut to we don't have to singularValuesOfR: there can
                                   // never be > 6 linearly-independent r vectors
      scheduled_for_deactivation[current_yf] = true;
    else
    {
      r_tmp.push_back(r[current_yf]);
      info = singularValuesOfR(r_tmp, s);
      if (info != 0)
        mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
                   "returned with error code ",
                   info);
      if (s[s.size() - 1] < _svd_tol * s[0])
      {
        scheduled_for_deactivation[current_yf] = true;
        r_tmp.pop_back();
        num_lin_dep--;
      }
      else
        num_kept_active++;
      if (num_lin_dep == 0 && num_active <= 6)
        // have taken out all the vectors that were linearly dependent
        // so no point continuing
        break;
    }
  }
 
  unsigned int old_active_number = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface])
    {
      if (scheduled_for_deactivation[old_active_number])
        deactivated_due_to_ld[surface] = true;
      old_active_number++;
    }
}

Referenced by MultiPlasticityLinearSystem::calculateRHS().

fddflowPotential_dintnl()

void MultiPlasticityDebugger::fddflowPotential_dintnl ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl, std::vector< RankTwoTensor > &  dr_dintnl  )

privatevirtualinherited

The finite-difference derivative of the flow potentials with respect to internal parameters.

Parameters

	stress	the stress at which to calculate the flow potential
	intnl	vector of internal parameters
[out]	dr_dintnl	the derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 612 of file MultiPlasticityDebugger.C.

{
  dr_dintnl.resize(_num_surfaces);
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  std::vector<RankTwoTensor> origr;
  flowPotential(stress, intnl, act, origr);
 
  std::vector<Real> intnlep;
  intnlep.resize(_num_models);
  for (unsigned model = 0; model < _num_models; ++model)
    intnlep[model] = intnl[model];
  Real ep;
  std::vector<RankTwoTensor> rep;
  unsigned int model;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    model = modelNumber(surface);
    ep = _fspb_debug_intnl_change[model];
    intnlep[model] += ep;
    flowPotential(stress, intnlep, act, rep);
    dr_dintnl[surface] = (rep[surface] - origr[surface]) / ep;
    intnlep[model] -= ep;
  }
}

Referenced by MultiPlasticityDebugger::checkDerivatives().

fddflowPotential_dstress()

void MultiPlasticityDebugger::fddflowPotential_dstress ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl, std::vector< RankFourTensor > &  dr_dstress  )

privatevirtualinherited

The finite-difference derivative of the flow potential(s) with respect to stress.

Parameters

	stress	the stress at which to calculate the flow potential
	intnl	vector of internal parameters
[out]	dr_dstress	the derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 583 of file MultiPlasticityDebugger.C.

{
  dr_dstress.assign(_num_surfaces, RankFourTensor());
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  Real ep = _fspb_debug_stress_change;
  RankTwoTensor stressep;
  std::vector<RankTwoTensor> rep, rep_minus;
  for (unsigned i = 0; i < 3; ++i)
    for (unsigned j = 0; j < 3; ++j)
    {
      stressep = stress;
      // do a central difference
      stressep(i, j) += ep / 2.0;
      flowPotential(stressep, intnl, act, rep);
      stressep(i, j) -= ep;
      flowPotential(stressep, intnl, act, rep_minus);
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        for (unsigned k = 0; k < 3; ++k)
          for (unsigned l = 0; l < 3; ++l)
            dr_dstress[surface](k, l, i, j) = (rep[surface](k, l) - rep_minus[surface](k, l)) / ep;
    }
}

Referenced by MultiPlasticityDebugger::checkDerivatives().

fddyieldFunction_dintnl()

void MultiPlasticityDebugger::fddyieldFunction_dintnl ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl, std::vector< Real > &  df_dintnl  )

privateinherited

The finite-difference derivative of yield function(s) with respect to internal parameter(s)

Parameters

	stress	the stress at which to calculate the yield function
	intnl	vector of internal parameters
[out]	df_dintnl	the derivative (or derivatives in the case of multisurface plasticity). df_dintnl[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 552 of file MultiPlasticityDebugger.C.

{
  df_dintnl.resize(_num_surfaces);
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  std::vector<Real> origf;
  yieldFunction(stress, intnl, act, origf);
 
  std::vector<Real> intnlep;
  intnlep.resize(_num_models);
  for (unsigned model = 0; model < _num_models; ++model)
    intnlep[model] = intnl[model];
  Real ep;
  std::vector<Real> fep;
  unsigned int model;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    model = modelNumber(surface);
    ep = _fspb_debug_intnl_change[model];
    intnlep[model] += ep;
    yieldFunction(stress, intnlep, act, fep);
    df_dintnl[surface] = (fep[surface] - origf[surface]) / ep;
    intnlep[model] -= ep;
  }
}

Referenced by MultiPlasticityDebugger::checkDerivatives().

fddyieldFunction_dstress()

void MultiPlasticityDebugger::fddyieldFunction_dstress ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl, std::vector< RankTwoTensor > &  df_dstress  )

privateinherited

The finite-difference derivative of yield function(s) with respect to stress.

Parameters

	stress	the stress at which to calculate the yield function
	intnl	vector of internal parameters
[out]	df_dstress	the derivative (or derivatives in the case of multisurface plasticity). df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Definition at line 524 of file MultiPlasticityDebugger.C.

{
  df_dstress.assign(_num_surfaces, RankTwoTensor());
 
  std::vector<bool> act;
  act.assign(_num_surfaces, true);
 
  Real ep = _fspb_debug_stress_change;
  RankTwoTensor stressep;
  std::vector<Real> fep, fep_minus;
  for (unsigned i = 0; i < 3; ++i)
    for (unsigned j = 0; j < 3; ++j)
    {
      stressep = stress;
      // do a central difference to attempt to capture discontinuities
      // such as those encountered in tensile and Mohr-Coulomb
      stressep(i, j) += ep / 2.0;
      yieldFunction(stressep, intnl, act, fep);
      stressep(i, j) -= ep;
      yieldFunction(stressep, intnl, act, fep_minus);
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        df_dstress[surface](i, j) = (fep[surface] - fep_minus[surface]) / ep;
    }
}

Referenced by MultiPlasticityDebugger::checkDerivatives().

fdJacobian()

void MultiPlasticityDebugger::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  )

privatevirtualinherited

The Jacobian calculated using finite differences.

The output should be equal to calculateJacobian(...) if everything is coded correctly.

Parameters

	stress	the stress at which to calculate the Jacobian
	intnl_old	the old values of internal variables (jacobian is inependent of these, but they are needed to do the finite-differencing cleanly)
	intnl	the vector of internal parameters at which to calculate the Jacobian
	pm	the plasticity multipliers at which to calculate the Jacobian
	delta_dp	plastic_strain - plastic_strain_old (Jacobian is independent of this, but it is needed to do the finite-differencing cleanly)
	E_inv	inverse of the elasticity tensor
	eliminate_ld	only calculate the Jacobian for the linearly independent constraints
[out]	jac	the finite-difference Jacobian

Definition at line 225 of file MultiPlasticityDebugger.C.

{
  std::vector<bool> active;
  active.assign(_num_surfaces, true);
 
  std::vector<bool> deactivated_due_to_ld;
  std::vector<bool> deactivated_due_to_ld_ep;
 
  std::vector<Real> orig_rhs;
  calculateRHS(stress,
               intnl_old,
               intnl,
               pm,
               delta_dp,
               orig_rhs,
               active,
               eliminate_ld,
               deactivated_due_to_ld); // this calculates RHS, and also set deactivated_due_to_ld.
                                       // The latter stays fixed for the rest of this routine
 
  unsigned int whole_system_size = 6 + _num_surfaces + _num_models;
  unsigned int system_size =
      orig_rhs.size(); // will be = whole_system_size if eliminate_ld = false, since all active=true
  jac.resize(system_size);
  for (unsigned row = 0; row < system_size; ++row)
    jac[row].assign(system_size, 0);
 
  std::vector<Real> rhs_ep;
  unsigned col = 0;
 
  RankTwoTensor stressep;
  RankTwoTensor delta_dpep;
  Real ep = _fspb_debug_stress_change;
  for (unsigned i = 0; i < 3; ++i)
    for (unsigned j = 0; j <= i; ++j)
    {
      stressep = stress;
      stressep(i, j) += ep;
      if (i != j)
        stressep(j, i) += ep;
      delta_dpep = delta_dp;
      for (unsigned k = 0; k < 3; ++k)
        for (unsigned l = 0; l < 3; ++l)
        {
          delta_dpep(k, l) -= E_inv(k, l, i, j) * ep;
          if (i != j)
            delta_dpep(k, l) -= E_inv(k, l, j, i) * ep;
        }
      active.assign(_num_surfaces, true);
      calculateRHS(stressep,
                   intnl_old,
                   intnl,
                   pm,
                   delta_dpep,
                   rhs_ep,
                   active,
                   false,
                   deactivated_due_to_ld_ep);
      unsigned row = 0;
      for (unsigned dof = 0; dof < whole_system_size; ++dof)
        if (dof_included(dof, deactivated_due_to_ld))
        {
          jac[row][col] =
              -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
          row++;
        }
      col++; // all of the first 6 columns are dof_included since they're stresses
    }
 
  std::vector<Real> pmep;
  pmep.resize(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    pmep[surface] = pm[surface];
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
  {
    if (!dof_included(6 + surface, deactivated_due_to_ld))
      continue;
    ep = _fspb_debug_pm_change[surface];
    pmep[surface] += ep;
    active.assign(_num_surfaces, true);
    calculateRHS(
        stress, intnl_old, intnl, pmep, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
    unsigned row = 0;
    for (unsigned dof = 0; dof < whole_system_size; ++dof)
      if (dof_included(dof, deactivated_due_to_ld))
      {
        jac[row][col] =
            -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
        row++;
      }
    pmep[surface] -= ep;
    col++;
  }
 
  std::vector<Real> intnlep;
  intnlep.resize(_num_models);
  for (unsigned model = 0; model < _num_models; ++model)
    intnlep[model] = intnl[model];
  for (unsigned model = 0; model < _num_models; ++model)
  {
    if (!dof_included(6 + _num_surfaces + model, deactivated_due_to_ld))
      continue;
    ep = _fspb_debug_intnl_change[model];
    intnlep[model] += ep;
    active.assign(_num_surfaces, true);
    calculateRHS(
        stress, intnl_old, intnlep, pm, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
    unsigned row = 0;
    for (unsigned dof = 0; dof < whole_system_size; ++dof)
      if (dof_included(dof, deactivated_due_to_ld))
      {
        jac[row][col] =
            -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
        row++;
      }
    intnlep[model] -= ep;
    col++;
  }
}

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

flowPotential()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the flow potential
	intnl	vector of internal parameters
	active	set of active constraints - only the active flow potentials are put into "r"
[out]	r	the flow potential (flow potentials in the multi-surface case)

Definition at line 178 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  r.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<RankTwoTensor> model_r;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->flowPotentialV(stress, intnl[model], model_r);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        r.push_back(model_r[*active_surface]);
    }
  }
}

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

hardPotential()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the hardening potential
	intnl	vector of internal parameters
	active	set of active constraints - only the active hardening potentials are put into "h"
[out]	h	the 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 260 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  h.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<Real> model_h;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->hardPotentialV(stress, intnl[model], model_h);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        h.push_back(model_h[*active_surface]);
    }
  }
}

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

incrementDumb()

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

protectedvirtual

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

Parameters

[in,out]	dumb_iteration	Used to set act bitwise - the "dumb" scheme tries all possible combinations of act until a successful return
[in]	dumb_order	dumb_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]	act	active constraints

Definition at line 1553 of file ComputeMultiPlasticityStress.C.

{
  dumb_iteration += 1;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    act[dumb_order[surface]] =
        (dumb_iteration &
         (1 << surface)); // returns true if the surface_th bit of dumb_iteration == 1
}

Referenced by changeScheme(), and returnMap().

initQpStatefulProperties()

void ComputeMultiPlasticityStress::initQpStatefulProperties ( )

protectedvirtual

Reimplemented from ComputeGeneralStressBase.

Definition at line 189 of file ComputeMultiPlasticityStress.C.

{
  ComputeStressBase::initQpStatefulProperties();
 
  _plastic_strain[_qp].zero();
 
  _intnl[_qp].assign(_num_models, 0);
 
  _yf[_qp].assign(_num_surfaces, 0);
 
  _dummy_pm.assign(_num_surfaces, 0);
 
  _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
  _linesearch_needed[_qp] = 0;
  _ld_encountered[_qp] = 0;
  _constraints_added[_qp] = 0;
 
  _n[_qp] = _n_input;
 
  if (_cosserat)
    (*_couple_stress)[_qp].zero();
 
  if (_fspb_debug == "jacobian")
  {
    checkDerivatives();
    mooseError("Finite-differencing completed.  Exiting with no error");
  }
}

lineSearch()

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

protectedvirtual

Performs a line search.

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

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

Parameters

[in,out]	nr_res2	The residual-squared
[out]	stress	The stress after returning to the yield surface
	intnl_old	The internal variables at the previous "time" step
[in,out]	intnl	The internal variables
[in,out]	pm	The plasticity multiplier(s) (consistency parameter(s))
	E_inv	inverse of the elasticity tensor
[in,out]	delta_dp	Change in plastic strain from start of "time" step to current configuration (plastic_strain - plastic_strain_old)
	dstress	Change in stress for a full Newton step
	dpm	Change in plasticity multiplier for a full Newton step
	dintnl	change in internal parameter(s) for a full Newton step
[in,out]	f	Yield function(s). In this routine, only the active constraints that are not deactivated_due_to_ld are contained in f.
[in,out]	epp	Plastic strain increment constraint
[in,out]	ic	Internal constraint. In this routine, only the active constraints that are not deactivated_due_to_ld are contained in ic.
	active	The active constraints.
	deactivated_due_to_ld	True if a constraint has temporarily been made deactive due to linear dependence.
[out]	linesearch_needed	True 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 1389 of file ComputeMultiPlasticityStress.C.

{
  // Line search algorithm straight out of "Numerical Recipes"
 
  bool success =
      true; // return value: will be false if linesearch couldn't reduce the residual-squared
 
  // Aim is to decrease residual2
 
  Real lam = 1.0; // the line-search parameter: 1.0 is a full Newton step
  Real lam_min =
      1E-10; // minimum value of lam allowed - perhaps this should be dynamically calculated?
  Real f0 = nr_res2;         // initial value of residual2
  Real slope = -2 * nr_res2; // "Numerical Recipes" uses -b*A*x, in order to check for roundoff, but
                             // i hope the nrStep would warn if there were problems.
  Real tmp_lam;              // cached value of lam used in quadratic & cubic line search
  Real f2 = nr_res2;         // cached value of f = residual2 used in the cubic in the line search
  Real lam2 = lam;           // cached value of lam used in the cubic in the line search
 
  // pm during the line-search
  std::vector<Real> ls_pm;
  ls_pm.resize(pm.size());
 
  // delta_dp during the line-search
  RankTwoTensor ls_delta_dp;
 
  // internal parameter during the line-search
  std::vector<Real> ls_intnl;
  ls_intnl.resize(intnl.size());
 
  // stress during the line-search
  RankTwoTensor ls_stress;
 
  // flow directions (not used in line search, but calculateConstraints returns this parameter)
  std::vector<RankTwoTensor> r;
 
  while (true)
  {
    // update the variables using this line-search parameter
    for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
      ls_pm[alpha] = pm[alpha] + dpm[alpha] * lam;
    ls_delta_dp = delta_dp - E_inv * dstress * lam;
    for (unsigned a = 0; a < intnl.size(); ++a)
      ls_intnl[a] = intnl[a] + dintnl[a] * lam;
    ls_stress = stress + dstress * lam;
 
    // calculate the new active yield functions, epp and active internal constraints
    calculateConstraints(ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
 
    // calculate the new residual-squared
    nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
 
    if (nr_res2 < f0 + 1E-4 * lam * slope)
      break;
    else if (lam < lam_min)
    {
      success = false;
      // restore plastic multipliers, yield functions, etc to original values
      for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
        ls_pm[alpha] = pm[alpha];
      ls_delta_dp = delta_dp;
      for (unsigned a = 0; a < intnl.size(); ++a)
        ls_intnl[a] = intnl[a];
      ls_stress = stress;
      calculateConstraints(
          ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
      nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
      break;
    }
    else if (lam == 1.0)
    {
      // model as a quadratic
      tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
    }
    else
    {
      // model as a cubic
      Real rhs1 = nr_res2 - f0 - lam * slope;
      Real rhs2 = f2 - f0 - lam2 * slope;
      Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
      Real b =
          (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
      if (a == 0)
        tmp_lam = -slope / (2.0 * b);
      else
      {
        Real disc = Utility::pow<2>(b) - 3 * a * slope;
        if (disc < 0)
          tmp_lam = 0.5 * lam;
        else if (b <= 0)
          tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
        else
          tmp_lam = -slope / (b + std::sqrt(disc));
      }
      if (tmp_lam > 0.5 * lam)
        tmp_lam = 0.5 * lam;
    }
    lam2 = lam;
    f2 = nr_res2;
    lam = std::max(tmp_lam, 0.1 * lam);
  }
 
  if (lam < 1.0)
    linesearch_needed = true;
 
  // assign the quantities found in the line-search
  // back to the originals
  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
    pm[alpha] = ls_pm[alpha];
  delta_dp = ls_delta_dp;
  for (unsigned a = 0; a < intnl.size(); ++a)
    intnl[a] = ls_intnl[a];
  stress = ls_stress;
 
  return success;
}

Referenced by singleStep().

modelNumber()

unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface )

protectedinherited

returns the model number, given the surface number

Definition at line 783 of file MultiPlasticityRawComponentAssembler.C.

{
  return _model_given_surface[surface];
}

Referenced by applyKuhnTucker(), MultiPlasticityLinearSystem::calculateJacobian(), canAddConstraints(), checkAdmissible(), checkKuhnTucker(), consistentTangentOperator(), MultiPlasticityDebugger::fddflowPotential_dintnl(), MultiPlasticityDebugger::fddyieldFunction_dintnl(), residual2(), returnMap(), and singleStep().

nrStep()

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

protectedvirtualinherited

Performs one Newton-Raphson step.

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

Parameters

	stress	Current value of stress
	intnl_old	The internal variables at the previous "time" step
	intnl	Current value of the internal variables
	pm	Current value of the plasticity multipliers (consistency parameters)
	E_inv	inverse of the elasticity tensor
	delta_dp	Current value of the plastic-strain increment (ie plastic_strain - plastic_strain_old)
[out]	dstress	The change in stress for a full Newton step
[out]	dpm	The change in all plasticity multipliers for a full Newton step
[out]	dintnl	The change in all internal variables for a full Newton step
	active	The active constraints
[out]	deactivated_due_to_ld	The constraints deactivated due to linear-dependence of the flow directions

Definition at line 614 of file MultiPlasticityLinearSystem.C.

{
  // Calculate RHS and Jacobian
  std::vector<Real> rhs;
  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
 
  std::vector<std::vector<Real>> jac;
  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
 
  // prepare for LAPACKgesv_ routine provided by PETSc
  PetscBLASInt system_size = rhs.size();
 
  std::vector<double> a(system_size * system_size);
  // Fill in the a "matrix" by going down columns
  unsigned ind = 0;
  for (int col = 0; col < system_size; ++col)
    for (int row = 0; row < system_size; ++row)
      a[ind++] = jac[row][col];
 
  PetscBLASInt nrhs = 1;
  std::vector<PetscBLASInt> ipiv(system_size);
  PetscBLASInt info;
  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
 
  if (info != 0)
    mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
               "routine returned with error code ",
               info);
 
  // Extract the results back to dstress, dpm and dintnl
  std::vector<bool> active_not_deact(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
 
  unsigned int dim = 3;
  ind = 0;
 
  for (unsigned i = 0; i < dim; ++i)
    for (unsigned j = 0; j <= i; ++j)
      dstress(i, j) = dstress(j, i) = rhs[ind++];
  dpm.assign(_num_surfaces, 0);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active_not_deact[surface])
      dpm[surface] = rhs[ind++];
  dintnl.assign(_num_models, 0);
  for (unsigned model = 0; model < _num_models; ++model)
    if (anyActiveSurfaces(model, active_not_deact))
      dintnl[model] = rhs[ind++];
 
  mooseAssert(static_cast<int>(ind) == system_size,
              "Incorrect extracting of changes from NR solution in nrStep");
}

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

numberActive()

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

protectedvirtual

counts the number of active constraints

Definition at line 1259 of file ComputeMultiPlasticityStress.C.

{
  unsigned num_active = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface])
      num_active++;
  return num_active;
}

Referenced by returnMap(), and singleStep().

outputAndCheckDebugParameters()

void MultiPlasticityDebugger::outputAndCheckDebugParameters ( )

inherited

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

Definition at line 53 of file MultiPlasticityDebugger.C.

{
  Moose::err << "Debug Parameters are as follows\n";
  Moose::err << "stress = \n";
  _fspb_debug_stress.print();
 
  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
      _fspb_debug_pm_change.size() != _num_surfaces ||
      _fspb_debug_intnl_change.size() != _num_models)
    mooseError("The debug parameters have the wrong size\n");
 
  Moose::err << "plastic multipliers =\n";
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    Moose::err << _fspb_debug_pm[surface] << "\n";
 
  Moose::err << "internal parameters =\n";
  for (unsigned model = 0; model < _num_models; ++model)
    Moose::err << _fspb_debug_intnl[model] << "\n";
 
  Moose::err << "finite-differencing parameter for stress-changes:\n"
             << _fspb_debug_stress_change << "\n";
  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    Moose::err << _fspb_debug_pm_change[surface] << "\n";
  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
  for (unsigned model = 0; model < _num_models; ++model)
    Moose::err << _fspb_debug_intnl_change[model] << "\n";
 
  Moose::err << std::flush;
}

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

plasticStep()

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

protectedvirtual

performs a plastic step

Parameters

	stress_old	The value of stress at the previous "time" step
[out]	stress	stress after returning to the yield surface
	intnl_old	The internal variables at the previous "time" step
[out]	intnl	internal variables after returning to the yield surface
	plastic_strain_old	The value of plastic strain at the previous "time" step
[out]	plastic_strain	plastic_strain after returning to the yield surface
	E_ijkl	The elasticity tensor.
	strain_increment	The applied strain increment
[out]	yf	All the yield functions at (stress, intnl)
[out]	iterations	The total number of Newton-Raphson iterations used
[out]	linesearch_needed	True if a linesearch was needed at any stage during the Newton-Raphson proceedure
[out]	ld_encountered	True if a linear-dependence of the flow directions was encountered at any stage during the Newton-Raphson proceedure
[out]	constraints_added	True if constraints were added into the active set at any stage during the Newton-Raphson proceedure
[out]	consistent_tangent_operator	The 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 466 of file ComputeMultiPlasticityStress.C.

{
  bool return_successful = false;
 
  // total number of Newton-Raphson iterations used
  unsigned int iter = 0;
 
  Real step_size = 1.0;
  Real time_simulated = 0.0;
 
  // the "good" variables hold the latest admissible stress
  // and internal parameters.
  RankTwoTensor stress_good = stress_old;
  RankTwoTensor plastic_strain_good = plastic_strain_old;
  std::vector<Real> intnl_good(_num_models);
  for (unsigned model = 0; model < _num_models; ++model)
    intnl_good[model] = intnl_old[model];
  std::vector<Real> yf_good(_num_surfaces);
 
  // Following is necessary because I want strain_increment to be "const"
  // but I also want to be able to subdivide an initial_stress
  RankTwoTensor this_strain_increment = strain_increment;
 
  RankTwoTensor dep = step_size * this_strain_increment;
 
  _cumulative_pm.assign(_num_surfaces, 0);
 
  unsigned int num_consecutive_successes = 0;
  while (time_simulated < 1.0 && step_size >= _min_stepsize)
  {
    iter = 0;
    return_successful = returnMap(stress_good,
                                  stress,
                                  intnl_good,
                                  intnl,
                                  plastic_strain_good,
                                  plastic_strain,
                                  E_ijkl,
                                  dep,
                                  yf,
                                  iter,
                                  step_size <= _max_stepsize_for_dumb,
                                  linesearch_needed,
                                  ld_encountered,
                                  constraints_added,
                                  time_simulated + step_size >= 1,
                                  consistent_tangent_operator,
                                  _cumulative_pm);
    iterations += iter;
 
    if (return_successful)
    {
      num_consecutive_successes += 1;
      time_simulated += step_size;
 
      if (time_simulated < 1.0) // this condition is just for optimization: if time_simulated=1 then
                                // the "good" quantities are no longer needed
      {
        stress_good = stress;
        plastic_strain_good = plastic_strain;
        for (unsigned model = 0; model < _num_models; ++model)
          intnl_good[model] = intnl[model];
        for (unsigned surface = 0; surface < _num_surfaces; ++surface)
          yf_good[surface] = yf[surface];
        if (num_consecutive_successes >= 2)
          step_size *= 1.2;
      }
      step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
    }
    else
    {
      step_size *= 0.5;
      num_consecutive_successes = 0;
      stress = stress_good;
      plastic_strain = plastic_strain_good;
      for (unsigned model = 0; model < _num_models; ++model)
        intnl[model] = intnl_good[model];
      yf.resize(_num_surfaces); // might have excited with junk
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        yf[surface] = yf_good[surface];
      dep = step_size * this_strain_increment;
    }
  }
 
  if (!return_successful)
  {
    if (_ignore_failures)
    {
      stress = stress_good;
      plastic_strain = plastic_strain_good;
      for (unsigned model = 0; model < _num_models; ++model)
        intnl[model] = intnl_good[model];
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        yf[surface] = yf_good[surface];
    }
    else
    {
      Moose::out << "After reducing the stepsize to " << step_size
                 << " with original strain increment with L2norm " << this_strain_increment.L2norm()
                 << " the returnMap algorithm failed" << std::endl;
 
      _fspb_debug_stress = stress_good + E_ijkl * dep;
      _fspb_debug_pm.assign(
          _num_surfaces,
          1); // this is chosen arbitrarily - please change if a more suitable value occurs to you!
      _fspb_debug_intnl.resize(_num_models);
      for (unsigned model = 0; model < _num_models; ++model)
        _fspb_debug_intnl[model] = intnl_good[model];
      checkDerivatives();
      checkJacobian(_my_elasticity_tensor.invSymm(), _intnl_old[_qp]);
      checkSolution(_my_elasticity_tensor.invSymm());
      mooseError("Exiting\n");
    }
  }
 
  return return_successful;
}

Referenced by computeQpStress().

postReturnMap()

void ComputeMultiPlasticityStress::postReturnMap ( )

protectedvirtual

Definition at line 337 of file ComputeMultiPlasticityStress.C.

{
  if (_n_supplied)
  {
    // this is a rotation matrix that will rotate "z" axis back to _n
    _rot = _rot.transpose();
 
    // rotate the tensors back to original frame where _n is correctly oriented
    _my_elasticity_tensor.rotate(_rot);
    _Jacobian_mult[_qp].rotate(_rot);
    _my_strain_increment.rotate(_rot);
    _stress[_qp].rotate(_rot);
    _plastic_strain[_qp].rotate(_rot);
    if (_cosserat)
    {
      _my_flexural_rigidity_tensor.rotate(_rot);
      (*_Jacobian_mult_couple)[_qp].rotate(_rot);
      _my_curvature.rotate(_rot);
      (*_couple_stress)[_qp].rotate(_rot);
    }
 
    // Rotate n by _rotation_increment
    for (const auto i : make_range(Moose::dim))
    {
      _n[_qp](i) = 0;
      for (const auto j : make_range(Moose::dim))
        _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
    }
  }
}

Referenced by computeQpStress().

preReturnMap()

void ComputeMultiPlasticityStress::preReturnMap ( )

protectedvirtual

Definition at line 318 of file ComputeMultiPlasticityStress.C.

{
  if (_n_supplied)
  {
    // this is a rotation matrix that will rotate _n to the "z" axis
    _rot = RotationMatrix::rotVecToZ(_n[_qp]);
 
    // rotate the tensors to this frame
    _my_elasticity_tensor.rotate(_rot);
    _my_strain_increment.rotate(_rot);
    if (_cosserat)
    {
      _my_flexural_rigidity_tensor.rotate(_rot);
      _my_curvature.rotate(_rot);
    }
  }
}

Referenced by computeQpStress().

quickStep()

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

protectedvirtual

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

Parameters

	stress_old	The value of stress at the previous "time" step
[out]	stress	If returnvalue=true then this is the returned value of stress. Otherwise, this is undefined
	intnl_old	The internal variables at the previous "time" step
[out]	intnl	If returnvalue=true then this is the value of the internal parameters after returning. Otherwise, this is undefined
[out]	pm	If returnvalue=true, this is the plastic multipliers needed to bring aout the return. Otherwise, this is undefined
	[in/out]	cumulative_pm If returnvalue=true, this is cumulative plastic multipliers, updated with pm. Otherwise, this is untouched by the algorithm
	plastic_strain_old	The value of plastic strain at the previous "time" step
[out]	plastic_strain	If returnvalue=true, this is the new plastic strain. Otherwise it is set to plastic_strain_old
	E_ijkl	The elasticity tensor.
	strain_increment	The applied strain increment
[out]	yf	If returnvalue=true, then all the yield functions at (stress, intnl). Otherwise, all the yield functions at (stress_old, intnl_old)
[out]	iterations	Number of NR iterations used, which is always zero in the current implementation.
	called_from	This can be called from computeQpStress, in which case it can actually provde an answer to the returnmap algorithm, or from returnMap in which case it is probably only providing an answer to a particular subdivision of the returnmap algorithm. The consistent tangent operator is calculated idfferently in each case
	final_step	The consistent tangent operator is calculated if this is true
[out]	consistent_tangent_operator	If final_step==true and returnvalue=true, then this is the consistent tangent operator d(stress_rate)/d(strain_rate). Otherwise it is undefined.

Returns

true if the (stress, intnl) are admissible, in which case the (stress_old, intnl_old) could have been admissible, or exactly one of the plastic models successfully used its custom returnMap function to provide the returned (stress, intnl) values and all other plastic models are admissible at that configuration. Or, false, then (stress_old, intnl_old) is not admissible according to >=1 plastic model and the custom returnMap functions failed in some way.

Definition at line 369 of file ComputeMultiPlasticityStress.C.

{
  iterations = 0;
 
  unsigned num_plastic_returns;
  RankTwoTensor delta_dp;
 
  // the following does the customized returnMap algorithm
  // for all the plastic models.
  unsigned custom_model = 0;
  bool successful_return = returnMapAll(stress_old + E_ijkl * strain_increment,
                                        intnl_old,
                                        E_ijkl,
                                        _epp_tol,
                                        stress,
                                        intnl,
                                        pm,
                                        cumulative_pm,
                                        delta_dp,
                                        yf,
                                        num_plastic_returns,
                                        custom_model);
 
  // the following updates the plastic_strain, when necessary
  // and calculates the consistent_tangent_operator, when necessary
  if (num_plastic_returns == 0)
  {
    // if successful_return = true, then a purely elastic step
    // if successful_return = false, then >=1 plastic model is in
    //    inadmissible zone and failed to return using its customized
    //    returnMap function.
    // In either case:
    plastic_strain = plastic_strain_old;
    if (successful_return && final_step)
    {
      if (called_from == computeQpStress_function)
        consistent_tangent_operator = E_ijkl;
      else // cannot necessarily use E_ijkl since different plastic models may have been active
           // during other substeps
        consistent_tangent_operator =
            consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
    }
    return successful_return;
  }
  else if (num_plastic_returns == 1 && successful_return)
  {
    // one model has successfully completed its custom returnMap algorithm
    // and the other models have signalled they are elastic at
    // the trial stress
    plastic_strain = plastic_strain_old + delta_dp;
    if (final_step)
    {
      if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
      {
        if (_tangent_operator_type == elastic)
          consistent_tangent_operator = E_ijkl;
        else
        {
          std::vector<Real> custom_model_pm;
          for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
            custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
          consistent_tangent_operator =
              _f[custom_model]->consistentTangentOperator(stress_old + E_ijkl * strain_increment,
                                                          intnl_old[custom_model],
                                                          stress,
                                                          intnl[custom_model],
                                                          E_ijkl,
                                                          custom_model_pm);
        }
      }
      else // cannot necessarily use the custom consistentTangentOperator since different plastic
           // models may have been active during other substeps or the custom model says not to use
           // its custom CTO algorithm
        consistent_tangent_operator =
            consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
    }
    return true;
  }
  else // presumably returnMapAll is incorrectly coded!
    mooseError("ComputeMultiPlasticityStress::quickStep   should not get here!");
}

Referenced by computeQpStress(), and returnMap().

reinstateLinearDependentConstraints()

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

protectedvirtual

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

Definition at line 1246 of file ComputeMultiPlasticityStress.C.

{
  bool reinstated_actives = false;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (deactivated_due_to_ld[surface])
      reinstated_actives = true;
 
  deactivated_due_to_ld.assign(_num_surfaces, false);
  return reinstated_actives;
}

Referenced by singleStep().

residual2()

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

protectedvirtual

The residual-squared.

Parameters

pm	the plastic multipliers for all constraints
f	the active yield function(s) (not including the ones that are deactivated_due_to_ld)
epp	the plastic strain increment constraint
ic	the active internal constraint(s) (not including the ones that are deactivated_due_to_ld)
active	true if constraint is active
deactivated_due_to_ld	true if constraint has been temporarily deactivated due to linear dependence of flow directions

Definition at line 1352 of file ComputeMultiPlasticityStress.C.

{
  Real nr_res2 = 0;
  unsigned ind = 0;
 
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (active[surface])
    {
      if (!deactivated_due_to_ld[surface])
      {
        if (!(pm[surface] == 0 && f[ind] <= 0))
          nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
      }
      else if (deactivated_due_to_ld[surface] && f[ind] > 0)
        nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
      ind++;
    }
 
  nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
 
  std::vector<bool> active_not_deact(_num_surfaces);
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
  ind = 0;
  for (unsigned model = 0; model < _num_models; ++model)
    if (anyActiveSurfaces(model, active_not_deact))
      nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
 
  return nr_res2;
}

Referenced by lineSearch(), and singleStep().

returnMap()

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

protectedvirtual

Implements the return map.

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

Parameters

	stress_old	The value of stress at the previous "time" step
[out]	stress	The stress after returning to the yield surface
	intnl_old	The internal variables at the previous "time" step
[out]	intnl	All the internal variables after returning to the yield surface
	plastic_strain_old	The value of plastic strain at the previous "time" step
[out]	plastic_strain	The value of plastic strain after returning to the yield surface
	E_ijkl	The elasticity tensor. If no plasticity then stress = stress_old + E_ijkl*strain_increment
	strain_increment	The applied strain increment
[out]	f	All the yield functions after returning to the yield surface
[out]	iter	The number of Newton-Raphson iterations used
	can_revert_to_dumb	If the _deactivation_scheme is set to revert to dumb, it will only be allowed to do so if this parameter is true
[out]	linesearch_needed	True if a linesearch was needed at any stage during the Newton-Raphson proceedure
[out]	ld_encountered	True if a linear-dependence of the flow directions was encountered at any stage during the Newton-Raphson proceedure
[out]	constraints_added	True if constraints were added into the active set at any stage during the Newton-Raphson proceedure
	final_step	Each 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_operator	The 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_pm	Upon 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 613 of file ComputeMultiPlasticityStress.C.

{
 
  // The "consistency parameters" (plastic multipliers)
  // Change in plastic strain in this timestep = pm*flowPotential
  // Each pm must be non-negative
  std::vector<Real> pm;
  pm.assign(_num_surfaces, 0.0);
 
  bool successful_return = quickStep(stress_old,
                                     stress,
                                     intnl_old,
                                     intnl,
                                     pm,
                                     cumulative_pm,
                                     plastic_strain_old,
                                     plastic_strain,
                                     E_ijkl,
                                     strain_increment,
                                     f,
                                     iter,
                                     consistent_tangent_operator,
                                     returnMap_function,
                                     final_step);
 
  if (successful_return)
    return successful_return;
 
  // Here we know that the trial stress and intnl_old
  // is inadmissible, and we have to return from those values
  // value to the yield surface.  There are three
  // types of constraints we have to satisfy, listed
  // below, and calculated in calculateConstraints(...)
  // f<=0, epp=0, ic=0 (up to tolerances), and these are
  //     guaranteed to hold if nr_res2<=0.5
  //
  // Kuhn-Tucker conditions must also be satisfied
  // These are:
  // pm>=0, which may not hold upon exit of the NR loops
  //     due to _deactivation_scheme!=optimized;
  // pm*f=0 (up to tolerances), which may not hold upon exit
  //     of the NR loops if a constraint got deactivated
  //     due to linear dependence, and then f<0, and its pm>0
 
  // Plastic strain constraint, L2 norm must be zero (up to a tolerance)
  RankTwoTensor epp;
 
  // Yield function constraint passed to this function as
  // std::vector<Real> & f
  // Each yield function must be <= 0 (up to tolerance)
  // Note that only the constraints that are active will be
  // contained in f until the final few lines of returnMap
 
  // Internal constraint(s), must be zero (up to a tolerance)
  // Note that only the constraints that are active will be
  // contained in ic.
  std::vector<Real> ic;
 
  // Record the stress before Newton-Raphson in case of failure-and-restart
  RankTwoTensor initial_stress = stress;
 
  iter = 0;
 
  // Initialize the set of active constraints
  // At this stage, the active constraints are
  // those that exceed their _f_tol
  // active constraints.
  std::vector<bool> act;
  buildActiveConstraints(f, stress, intnl, E_ijkl, act);
 
  // Inverse of E_ijkl (assuming symmetric)
  RankFourTensor E_inv = E_ijkl.invSymm();
 
  // convenience variable that holds the change in plastic strain incurred during the return
  // delta_dp = plastic_strain - plastic_strain_old
  // delta_dp = E^{-1}*(initial_stress - stress), where initial_stress = E*(strain -
  // plastic_strain_old)
  RankTwoTensor delta_dp = RankTwoTensor();
 
  // whether single step was successful (whether line search was successful, and whether turning off
  // constraints was successful)
  bool single_step_success = true;
 
  // deactivation scheme
  DeactivationSchemeEnum deact_scheme = _deactivation_scheme;
 
  // For complicated deactivation schemes we have to record the initial active set
  std::vector<bool> initial_act;
  initial_act.resize(_num_surfaces);
  if (_deactivation_scheme == optimized_to_safe ||
      _deactivation_scheme == optimized_to_safe_to_dumb ||
      _deactivation_scheme == optimized_to_dumb)
  {
    // if "optimized" fails we can change the deactivation scheme to "safe", etc
    deact_scheme = optimized;
    for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      initial_act[surface] = act[surface];
  }
 
  if (_deactivation_scheme == safe_to_dumb)
    deact_scheme = safe;
 
  // For "dumb" deactivation, the active set takes all combinations until a solution is found
  int dumb_iteration = 0;
  std::vector<unsigned int> dumb_order;
 
  if (_deactivation_scheme == dumb ||
      (_deactivation_scheme == optimized_to_safe_to_dumb && can_revert_to_dumb) ||
      (_deactivation_scheme == safe_to_dumb && can_revert_to_dumb) ||
      (_deactivation_scheme == optimized_to_dumb && can_revert_to_dumb))
    buildDumbOrder(stress, intnl, dumb_order);
 
  if (_deactivation_scheme == dumb)
  {
    incrementDumb(dumb_iteration, dumb_order, act);
    yieldFunction(stress, intnl, act, f);
  }
 
  // To avoid any re-trials of "act" combinations that
  // we've already tried and rejected, i record the
  // combinations in actives_tried
  std::set<unsigned int> actives_tried;
  actives_tried.insert(activeCombinationNumber(act));
 
  // The residual-squared that the line-search will reduce
  // Later it will get contributions from epp and ic, but
  // at present these are zero
  Real nr_res2 = 0;
  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
    if (act[surface])
      nr_res2 += 0.5 * Utility::pow<2>(f[surface] / _f[modelNumber(surface)]->_f_tol);
 
  successful_return = false;
 
  bool still_finding_solution = true;
  while (still_finding_solution)
  {
    single_step_success = true;
    unsigned int local_iter = 0;
 
    // The Newton-Raphson loops
    while (nr_res2 > 0.5 && local_iter++ < _max_iter && single_step_success)
      single_step_success = singleStep(nr_res2,
                                       stress,
                                       intnl_old,
                                       intnl,
                                       pm,
                                       delta_dp,
                                       E_inv,
                                       f,
                                       epp,
                                       ic,
                                       act,
                                       deact_scheme,
                                       linesearch_needed,
                                       ld_encountered);
 
    bool nr_good = (nr_res2 <= 0.5 && local_iter <= _max_iter && single_step_success);
 
    iter += local_iter;
 
    // 'act' might have changed due to using deact_scheme = optimized, so
    actives_tried.insert(activeCombinationNumber(act));
 
    if (!nr_good)
    {
      // failure of NR routine.
      // We might be able to change the deactivation_scheme and
      // then re-try the NR starting from the initial values
      // Or, if deact_scheme == "dumb", we just increarse the
      // dumb_iteration number and re-try
      bool change_scheme = false;
      bool increment_dumb = false;
      change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
      if (!change_scheme && deact_scheme == dumb)
        increment_dumb = canIncrementDumb(dumb_iteration);
 
      still_finding_solution = (change_scheme || increment_dumb);
 
      if (change_scheme)
        changeScheme(initial_act,
                     can_revert_to_dumb,
                     initial_stress,
                     intnl_old,
                     deact_scheme,
                     act,
                     dumb_iteration,
                     dumb_order);
 
      if (increment_dumb)
        incrementDumb(dumb_iteration, dumb_order, act);
 
      if (!still_finding_solution)
      {
        // we cannot change the scheme, or have run out of "dumb" options
        successful_return = false;
        break;
      }
    }
 
    bool kt_good = false;
    if (nr_good)
    {
      // check Kuhn-Tucker
      kt_good = checkKuhnTucker(f, pm, act);
      if (!kt_good)
      {
        if (deact_scheme != dumb)
        {
          applyKuhnTucker(f, pm, act);
 
          // true if we haven't tried this active set before
          still_finding_solution =
              (actives_tried.find(activeCombinationNumber(act)) == actives_tried.end());
          if (!still_finding_solution)
          {
            // must have tried turning off the constraints already.
            // so try changing the scheme
            if (canChangeScheme(deact_scheme, can_revert_to_dumb))
            {
              still_finding_solution = true;
              changeScheme(initial_act,
                           can_revert_to_dumb,
                           initial_stress,
                           intnl_old,
                           deact_scheme,
                           act,
                           dumb_iteration,
                           dumb_order);
            }
          }
        }
        else
        {
          bool increment_dumb = false;
          increment_dumb = canIncrementDumb(dumb_iteration);
          still_finding_solution = increment_dumb;
          if (increment_dumb)
            incrementDumb(dumb_iteration, dumb_order, act);
        }
 
        if (!still_finding_solution)
        {
          // have tried turning off the constraints already,
          // or have run out of "dumb" options
          successful_return = false;
          break;
        }
      }
    }
 
    bool admissible = false;
    if (nr_good && kt_good)
    {
      // check admissible
      std::vector<Real> all_f;
      if (_num_surfaces == 1)
        admissible = true; // for a single surface if NR has exited successfully then (stress,
                           // intnl) must be admissible
      else
        admissible = checkAdmissible(stress, intnl, all_f);
 
      if (!admissible)
      {
        // Not admissible.
        // We can try adding constraints back in
        // We can try changing the deactivation scheme
        // Or, if deact_scheme == dumb, just increase dumb_iteration
        bool add_constraints = canAddConstraints(act, all_f);
        if (add_constraints)
        {
          constraints_added = true;
          std::vector<bool> act_plus(_num_surfaces,
                                     false); // "act" with the positive constraints added in
          for (unsigned surface = 0; surface < _num_surfaces; ++surface)
            if (act[surface] ||
                (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol)))
              act_plus[surface] = true;
          if (actives_tried.find(activeCombinationNumber(act_plus)) == actives_tried.end())
          {
            // haven't tried this combination of actives yet
            constraints_added = true;
            for (unsigned surface = 0; surface < _num_surfaces; ++surface)
              act[surface] = act_plus[surface];
          }
          else
            add_constraints = false; // haven't managed to add a new combination
        }
 
        bool change_scheme = false;
        bool increment_dumb = false;
 
        if (!add_constraints)
          change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
 
        if (!add_constraints && !change_scheme && deact_scheme == dumb)
          increment_dumb = canIncrementDumb(dumb_iteration);
 
        still_finding_solution = (add_constraints || change_scheme || increment_dumb);
 
        if (change_scheme)
          changeScheme(initial_act,
                       can_revert_to_dumb,
                       initial_stress,
                       intnl_old,
                       deact_scheme,
                       act,
                       dumb_iteration,
                       dumb_order);
 
        if (increment_dumb)
          incrementDumb(dumb_iteration, dumb_order, act);
 
        if (!still_finding_solution)
        {
          // we cannot change the scheme, or have run out of "dumb" options
          successful_return = false;
          break;
        }
      }
    }
 
    successful_return = (nr_good && admissible && kt_good);
    if (successful_return)
      break;
 
    if (still_finding_solution)
    {
      stress = initial_stress;
      delta_dp = RankTwoTensor(); // back to zero change in plastic strain
      for (unsigned model = 0; model < _num_models; ++model)
        intnl[model] = intnl_old[model]; // back to old internal params
      pm.assign(_num_surfaces, 0.0);     // back to zero plastic multipliers
 
      unsigned num_active = numberActive(act);
      if (num_active == 0)
      {
        successful_return = false;
        break; // failure
      }
 
      actives_tried.insert(activeCombinationNumber(act));
 
      // Since "act" set has changed, either by changing deact_scheme, or by KT failing, so need to
      // re-calculate nr_res2
      yieldFunction(stress, intnl, act, f);
 
      nr_res2 = 0;
      unsigned ind = 0;
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        if (act[surface])
        {
          if (f[ind] > _f[modelNumber(surface)]->_f_tol)
            nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
          ind++;
        }
    }
  }
 
  // returned, with either success or failure
  if (successful_return)
  {
    plastic_strain += delta_dp;
 
    for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      cumulative_pm[surface] += pm[surface];
 
    if (final_step)
      consistent_tangent_operator =
          consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
 
    if (f.size() != _num_surfaces)
    {
      // at this stage f.size() = num_active, but we need to return with all the yield functions
      // evaluated, so:
      act.assign(_num_surfaces, true);
      yieldFunction(stress, intnl, act, f);
    }
  }
 
  return successful_return;
}

Referenced by plasticStep().

returnMapAll()

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

protectedinherited

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

Performs a returnMap for each plastic model.

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

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

Parameters

	trial_stress	the trial stress
	intnl_old	the old values of the internal parameters
	E_ijkl	the elasticity tensor
	ep_plastic_tolerance	the tolerance on the plastic strain
[out]	stress	is set to trial_stress in case (A) or (C), and the returned value of stress in case (B).
[out]	intnl	is set to intnl_old in case (A) or (C), and the returned value of intnl in case (B)
[out]	pm	Zero in case (A) or (C), otherwise the plastic multipliers needed to bring about the returnMap in case (B)
	[in/out]	cumulative_pm cumulative plastic multipliers, updated in case (B), otherwise left untouched
[out]	delta_dp	is unchanged in case (A) or (C), and is set to the change in plastic strain in case(B)
[out]	yf	will contain the yield function values at (stress, intnl)
[out]	num_successful_plastic_returns	will be 0 for (A) and (C), and 1 for (B)

Returns

true in case (A) and (B), and false in case (C)

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

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

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

Definition at line 597 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl_old.size() == _num_models,
              "returnMapAll: Incorrect size of internal parameters");
  mooseAssert(intnl.size() == _num_models, "returnMapAll: Incorrect size of internal parameters");
  mooseAssert(pm.size() == _num_surfaces, "returnMapAll: Incorrect size of pm");
 
  num_successful_plastic_returns = 0;
  yf.resize(0);
  pm.assign(_num_surfaces, 0.0);
 
  RankTwoTensor returned_stress; // each model will give a returned_stress.  if only one model is
                                 // plastically active, i set stress=returned_stress, so as to
                                 // record this returned value
  std::vector<Real> model_f;
  RankTwoTensor model_delta_dp;
  std::vector<Real> model_pm;
  bool trial_stress_inadmissible;
  bool successful_return = true;
  unsigned the_single_plastic_model = 0;
  bool using_custom_return_map = true;
 
  // run through all the plastic models, performing their
  // returnMap algorithms.
  // If one finds (trial_stress, intnl) inadmissible and
  // successfully returns, break from the loop to evaluate
  // all the others at that returned stress
  for (unsigned model = 0; model < _num_models; ++model)
  {
    if (using_custom_return_map)
    {
      model_pm.assign(_f[model]->numberSurfaces(), 0.0);
      bool model_returned = _f[model]->returnMap(trial_stress,
                                                 intnl_old[model],
                                                 E_ijkl,
                                                 ep_plastic_tolerance,
                                                 returned_stress,
                                                 intnl[model],
                                                 model_pm,
                                                 model_delta_dp,
                                                 model_f,
                                                 trial_stress_inadmissible);
      if (!trial_stress_inadmissible)
      {
        // in the elastic zone: record the yield-function values (returned_stress, intnl, model_pm
        // and model_delta_dp are undefined)
        for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
             ++model_surface)
          yf.push_back(model_f[model_surface]);
      }
      else if (trial_stress_inadmissible && !model_returned)
      {
        // in the plastic zone, and the customized returnMap failed
        // for some reason (or wasn't implemented).  The coder
        // should have correctly returned model_f(trial_stress, intnl_old)
        // so record them
        // (returned_stress, intnl, model_pm and model_delta_dp are undefined)
        for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
             ++model_surface)
          yf.push_back(model_f[model_surface]);
        // now there's almost zero point in using the custom
        // returnMap functions
        using_custom_return_map = false;
        successful_return = false;
      }
      else
      {
        // in the plastic zone, and the customized returnMap
        // succeeded.
        // record the first returned_stress and delta_dp if everything is going OK
        // as they could be the actual answer
        if (trial_stress_inadmissible)
          num_successful_plastic_returns++;
        the_single_plastic_model = model;
        stress = returned_stress;
        // note that i break here, and don't push_back
        // model_f to yf.  So now yf contains only the values of
        // model_f from previous models to the_single_plastic_model
        // also i don't set delta_dp = model_delta_dp yet, because
        // i might find problems later on
        // also, don't increment cumulative_pm for the same reason
 
        break;
      }
    }
    else
    {
      // not using custom returnMap functions because one
      // has already failed and that one said trial_stress
      // was inadmissible.  So now calculate the yield functions
      // at the trial stress
      _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
      for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
        yf.push_back(model_f[model_surface]);
    }
  }
 
  if (num_successful_plastic_returns == 0)
  {
    // here either all the models were elastic (successful_return=true),
    // or some were plastic and either the customized returnMap failed
    // or wasn't implemented (successful_return=false).
    // In either case, have to set the following:
    stress = trial_stress;
    for (unsigned model = 0; model < _num_models; ++model)
      intnl[model] = intnl_old[model];
    return successful_return;
  }
 
  // Now we know that num_successful_plastic_returns == 1 and all the other
  // models (with model number < the_single_plastic_model) must have been
  // admissible at (trial_stress, intnl).  However, all models might
  // not be admissible at (trial_stress, intnl), so must check that
  std::vector<Real> yf_at_returned_stress(0);
  bool all_admissible = true;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    if (model == the_single_plastic_model)
    {
      // no need to spend time calculating the yield function: we know its admissible
      for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
        yf_at_returned_stress.push_back(model_f[model_surface]);
      continue;
    }
    _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
    for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
    {
      if (model_f[model_surface] > _f[model]->_f_tol)
        // bummer, this model is not admissible at the returned_stress
        all_admissible = false;
      yf_at_returned_stress.push_back(model_f[model_surface]);
    }
    if (!all_admissible)
      // no point in continuing computing yield functions
      break;
  }
 
  if (!all_admissible)
  {
    // we tried using the returned value of stress predicted by
    // the_single_plastic_model, but it wasn't admissible according
    // to other plastic models.  We need to set:
    stress = trial_stress;
    for (unsigned model = 0; model < _num_models; ++model)
      intnl[model] = intnl_old[model];
    // and calculate the remainder of the yield functions at trial_stress
    for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
    {
      _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
      for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
        yf.push_back(model_f[model_surface]);
    }
    num_successful_plastic_returns = 0;
    return false;
  }
 
  // So the customized returnMap algorithm can provide a returned
  // (stress, intnl) configuration, and that is admissible according
  // to all plastic models
  yf.resize(0);
  for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
    yf.push_back(yf_at_returned_stress[surface]);
  delta_dp = model_delta_dp;
  for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
       ++model_surface)
  {
    cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
        model_pm[model_surface];
    pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
  }
  custom_model = the_single_plastic_model;
  return true;
}

Referenced by quickStep().

rot()

RankTwoTensor ComputeMultiPlasticityStress::rot ( const RankTwoTensor &  tens )

private

Definition at line 310 of file ComputeMultiPlasticityStress.C.

{
  if (!_n_supplied)
    return tens;
  return tens.rotated(_rot);
}

Referenced by computeQpStress().

singleStep()

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

protectedvirtual

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

Parameters

[in,out]	nr_res2	Residual-squared that the line-search will reduce
[in,out]	stress	stress
[in]	intnl_old	old values of the internal parameters
[in,out]	intnl	internal parameters
[in,out]	pm	plastic multipliers
[in,out]	delta_dp	Change in plastic strain from start of "time" step to current configuration (plastic_strain - plastic_strain_old)
[in]	E_inv	Inverse of the elasticity tensor
[in,out]	f	Yield function(s). Upon successful exit only the active constraints are contained in f
[in,out]	epp	Plastic strain increment constraint
[in,out]	ic	Internal constraint. Upon successful exit only the active constraints are contained in ic
	active	The active constraints. This is may be modified, depending upon deactivation_scheme
	deactivation_scheme	The scheme used for deactivating constraints
[out]	linesearch_needed	True if a linesearch was employed during this Newton-Raphson step
[out]	ld_encountered	True 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 1080 of file ComputeMultiPlasticityStress.C.

{
  bool successful_step; // return value
 
  Real nr_res2_before_step = nr_res2;
  RankTwoTensor stress_before_step;
  std::vector<Real> intnl_before_step;
  std::vector<Real> pm_before_step;
  RankTwoTensor delta_dp_before_step;
 
  if (deactivation_scheme == optimized)
  {
    // we potentially use the "before_step" quantities, so record them here
    stress_before_step = stress;
    intnl_before_step.resize(_num_models);
    for (unsigned model = 0; model < _num_models; ++model)
      intnl_before_step[model] = intnl[model];
    pm_before_step.resize(_num_surfaces);
    for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      pm_before_step[surface] = pm[surface];
    delta_dp_before_step = delta_dp;
  }
 
  // During the Newton-Raphson procedure, we'll be
  // changing the following parameters in order to
  // (attempt to) satisfy the constraints.
  RankTwoTensor dstress; // change in stress
  std::vector<Real> dpm; // change in plasticity multipliers ("consistency parameters").  For ALL
                         // contraints (active and deactive)
  std::vector<Real>
      dintnl; // change in internal parameters.  For ALL internal params (active and deactive)
 
  // The constraints that have been deactivated for this NR step
  // due to the flow directions being linearly dependent
  std::vector<bool> deact_ld;
  deact_ld.assign(_num_surfaces, false);
 
  /* After NR and linesearch, if _deactivation_scheme == "optimized", the
   * active plasticity multipliers are checked for non-negativity.  If some
   * are negative then they are deactivated forever, and the NR step is
   * re-done starting from the _before_step quantities recorded above
   */
  bool constraints_changing = true;
  bool reinstated_actives;
  while (constraints_changing)
  {
    // calculate dstress, dpm and dintnl for one full Newton-Raphson step
    nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
 
    for (unsigned surface = 0; surface < deact_ld.size(); ++surface)
      if (deact_ld[surface])
      {
        ld_encountered = true;
        break;
      }
 
    // perform a line search
    // The line-search will exit with updated values
    successful_step = lineSearch(nr_res2,
                                 stress,
                                 intnl_old,
                                 intnl,
                                 pm,
                                 E_inv,
                                 delta_dp,
                                 dstress,
                                 dpm,
                                 dintnl,
                                 f,
                                 epp,
                                 ic,
                                 active,
                                 deact_ld,
                                 linesearch_needed);
 
    if (!successful_step)
      // completely bomb out
      return successful_step;
 
    // See if any active constraints need to be removed, and the step re-done
    constraints_changing = false;
    if (deactivation_scheme == optimized)
    {
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        if (active[surface] && pm[surface] < 0.0)
          constraints_changing = true;
    }
 
    if (constraints_changing)
    {
      stress = stress_before_step;
      delta_dp = delta_dp_before_step;
      nr_res2 = nr_res2_before_step;
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
      {
        if (active[surface] && pm[surface] < 0.0)
        {
          // turn off the constraint forever
          active[surface] = false;
          pm_before_step[surface] = 0.0;
          intnl_before_step[modelNumber(surface)] =
              intnl_old[modelNumber(surface)]; // don't want to muck-up hardening!
        }
        intnl[modelNumber(surface)] = intnl_before_step[modelNumber(surface)];
        pm[surface] = pm_before_step[surface];
      }
      if (numberActive(active) == 0)
      {
        // completely bomb out
        successful_step = false;
        return successful_step;
      }
    }
 
    // reinstate any active values that have been turned off due to linear-dependence
    reinstated_actives = reinstateLinearDependentConstraints(deact_ld);
  } // ends "constraints_changing" loop
 
  // if active constraints were reinstated then nr_res2 needs to be re-calculated so it is correct
  // upson returning
  if (reinstated_actives)
  {
    bool completely_converged = true;
    if (successful_step && nr_res2 < 0.5)
    {
      // Here we have converged to the correct solution if
      // all the yield functions are < 0.  Excellent!
      //
      // This is quite tricky - perhaps i can refactor to make it more obvious.
      //
      // Because actives are now reinstated, the residual2
      // calculation below will give nr_res2 > 0.5, because it won't
      // realise that we only had to set the active-but-not-deactivated f = 0,
      // and not the entire active set.  If we pass that nr_res2 back from
      // this function then the calling function will not realise we've converged!
      // Therefore, check for this case
      unsigned ind = 0;
      for (unsigned surface = 0; surface < _num_surfaces; ++surface)
        if (active[surface])
          if (f[ind++] > _f[modelNumber(surface)]->_f_tol)
            completely_converged = false;
    }
    else
      completely_converged = false;
 
    if (!completely_converged)
      nr_res2 = residual2(pm, f, epp, ic, active, deact_ld);
  }
 
  return successful_step;
}

Referenced by returnMap().

singularValuesOfR()

int MultiPlasticityLinearSystem::singularValuesOfR ( const std::vector< RankTwoTensor > &  r, std::vector< Real > &  s  )

privatevirtualinherited

Performs a singular-value decomposition of r and returns the singular values.

Example: If r has size 5 then the singular values of the following matrix are returned: ( r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2) ) ( r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2) ) a = ( r[2](0,0) r[2](0,1) r[2](0,2) r[2](1,1) r[2](1,2) r[2](2,2) ) ( r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2) ) ( r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2) )

Parameters

	r	The flow directions
[out]	s	The singular values

Returns

The return value from the PETSc LAPACK gesvd reoutine

Definition at line 45 of file MultiPlasticityLinearSystem.C.

{
  PetscBLASInt bm = r.size();
  PetscBLASInt bn = 6;
 
  s.resize(std::min(bm, bn));
 
  // prepare for gesvd or gesdd routine provided by PETSc
  // Want to find the singular values of matrix
  //     (  r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2)  )
  //     (  r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2)  )
  // a = (  r[2](0,0) r[2](0,1) r[2](0,2) r[2](1,1) r[2](1,2) r[2](2,2)  )
  //     (  r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2)  )
  //     (  r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2)  )
  // bm = 5
 
  std::vector<double> a(bm * 6);
  // Fill in the a "matrix" by going down columns
  unsigned ind = 0;
  for (int col = 0; col < 3; ++col)
    for (int row = 0; row < bm; ++row)
      a[ind++] = r[row](0, col);
  for (int col = 3; col < 5; ++col)
    for (int row = 0; row < bm; ++row)
      a[ind++] = r[row](1, col - 2);
  for (int row = 0; row < bm; ++row)
    a[ind++] = r[row](2, 2);
 
  // u and vt are dummy variables because they won't
  // get referenced due to the "N" and "N" choices
  PetscBLASInt sizeu = 1;
  std::vector<double> u(sizeu);
  PetscBLASInt sizevt = 1;
  std::vector<double> vt(sizevt);
 
  PetscBLASInt sizework =
      16 * (bm + 6); // this is above the lowerbound specified in the LAPACK doco
  std::vector<double> work(sizework);
 
  PetscBLASInt info;
 
  LAPACKgesvd_("N",
               "N",
               &bm,
               &bn,
               &a[0],
               &bm,
               &s[0],
               &u[0],
               &sizeu,
               &vt[0],
               &sizevt,
               &work[0],
               &sizework,
               &info);
 
  return info;
}

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

validParams()

InputParameters ComputeMultiPlasticityStress::validParams ( )

static

Definition at line 22 of file ComputeMultiPlasticityStress.C.

{
  InputParameters params = ComputeStressBase::validParams();
  params += MultiPlasticityDebugger::validParams();
  params.addClassDescription("Base class for multi-surface finite-strain plasticity");
  params.addRangeCheckedParam<unsigned int>("max_NR_iterations",
                                            20,
                                            "max_NR_iterations>0",
                                            "Maximum number of Newton-Raphson iterations allowed");
  params.addRequiredParam<Real>("ep_plastic_tolerance",
                                "The Newton-Raphson process is only deemed "
                                "converged if the plastic strain increment "
                                "constraints have L2 norm less than this.");
  params.addRangeCheckedParam<Real>(
      "min_stepsize",
      0.01,
      "min_stepsize>0 & min_stepsize<=1",
      "If ordinary Newton-Raphson + line-search fails, then the applied strain increment is "
      "subdivided, and the return-map is tried again.  This parameter is the minimum fraction of "
      "applied strain increment that may be applied before the algorithm gives up entirely");
  params.addRangeCheckedParam<Real>("max_stepsize_for_dumb",
                                    0.01,
                                    "max_stepsize_for_dumb>0 & max_stepsize_for_dumb<=1",
                                    "If your deactivation_scheme is 'something_to_dumb', then "
                                    "'dumb' will only be used if the stepsize falls below this "
                                    "value.  This parameter is useful because the 'dumb' scheme is "
                                    "computationally expensive");
  MooseEnum deactivation_scheme("optimized safe dumb optimized_to_safe safe_to_dumb "
                                "optimized_to_safe_to_dumb optimized_to_dumb",
                                "optimized");
  params.addParam<MooseEnum>(
      "deactivation_scheme",
      deactivation_scheme,
      "Scheme by which constraints are deactivated.  (NOTE: This is irrelevant if there is only "
      "one yield surface.)  safe: return to the yield surface and then deactivate constraints with "
      "negative plasticity multipliers.  optimized: deactivate a constraint as soon as its "
      "plasticity multiplier becomes negative.  dumb: iteratively try all combinations of active "
      "constraints until the solution is found.  You may specify fall-back options.  Eg "
      "optimized_to_safe: first use 'optimized', and if that fails, try the return with 'safe'.");
  params.addParam<RealVectorValue>(
      "transverse_direction",
      "If this parameter is provided, before the return-map algorithm is "
      "called a rotation is performed so that the 'z' axis in the new "
      "frame lies along the transverse_direction in the original frame.  "
      "After returning, the inverse rotation is performed.  The "
      "transverse_direction will itself rotate with large strains.  This "
      "is so that transversely-isotropic plasticity models may be easily "
      "defined in the frame where the isotropy holds in the x-y plane.");
  params.addParam<bool>("ignore_failures",
                        false,
                        "The return-map algorithm will return with the best admissible "
                        "stresses and internal parameters that it can, even if they don't "
                        "fully correspond to the applied strain increment.  To speed "
                        "computations, this flag can be set to true, the max_NR_iterations "
                        "set small, and the min_stepsize large.");
  MooseEnum tangent_operator("elastic linear nonlinear", "nonlinear");
  params.addParam<MooseEnum>("tangent_operator",
                             tangent_operator,
                             "Type of tangent operator to return.  'elastic': return the "
                             "elasticity tensor.  'linear': return the consistent tangent operator "
                             "that is correct for plasticity with yield functions linear in "
                             "stress.  'nonlinear': return the full, general consistent tangent "
                             "operator.  The calculations assume the hardening potentials are "
                             "independent of stress and hardening parameters.");
  params.addParam<bool>("perform_finite_strain_rotations",
                        true,
                        "Tensors are correctly rotated in "
                        "finite-strain simulations.  For "
                        "optimal performance you can set "
                        "this to 'false' if you are only "
                        "ever using small strains");
  params.addClassDescription("Material for multi-surface finite-strain plasticity");
  return params;
}

yieldFunction()

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

protectedvirtualinherited

The active yield function(s)

Parameters

	stress	the stress at which to calculate the yield function
	intnl	vector of internal parameters
	active	set of active constraints - only the active yield functions are put into "f"
[out]	f	the yield function (or functions in the case of multisurface plasticity)

Definition at line 96 of file MultiPlasticityRawComponentAssembler.C.

{
  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
 
  f.resize(0);
  std::vector<unsigned int> active_surfaces_of_model;
  std::vector<unsigned int>::iterator active_surface;
  std::vector<Real> model_f;
  for (unsigned model = 0; model < _num_models; ++model)
  {
    activeModelSurfaces(model, active, active_surfaces_of_model);
    if (active_surfaces_of_model.size() > 0)
    {
      _f[model]->yieldFunctionV(stress, intnl[model], model_f);
      for (active_surface = active_surfaces_of_model.begin();
           active_surface != active_surfaces_of_model.end();
           ++active_surface)
        f.push_back(model_f[*active_surface]);
    }
  }
}

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

Member Data Documentation

_base_name

const std::string ComputeGeneralStressBase::_base_name

protectedinherited

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

Definition at line 43 of file ComputeGeneralStressBase.h.

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

_constraints_added

MaterialProperty<Real>& ComputeMultiPlasticityStress::_constraints_added

protected

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

Definition at line 127 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_cosserat

bool ComputeMultiPlasticityStress::_cosserat

protected

whether Cosserat mechanics should be used

Definition at line 151 of file ComputeMultiPlasticityStress.h.

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

_couple_stress

MaterialProperty<RankTwoTensor>* const ComputeMultiPlasticityStress::_couple_stress

protected

the Cosserat couple-stress

Definition at line 160 of file ComputeMultiPlasticityStress.h.

_couple_stress_old

const MaterialProperty<RankTwoTensor>* const ComputeMultiPlasticityStress::_couple_stress_old

protected

the old value of Cosserat couple-stress

Definition at line 163 of file ComputeMultiPlasticityStress.h.

_cumulative_pm

std::vector<Real> ComputeMultiPlasticityStress::_cumulative_pm

protected

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

This is used for calculating the consistent tangent operator

Definition at line 66 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

_curvature

const MaterialProperty<RankTwoTensor>* const ComputeMultiPlasticityStress::_curvature

protected

The Cosserat curvature strain.

Definition at line 154 of file ComputeMultiPlasticityStress.h.

_deactivation_scheme

enum ComputeMultiPlasticityStress::DeactivationSchemeEnum ComputeMultiPlasticityStress::_deactivation_scheme

protected

Referenced by canChangeScheme(), changeScheme(), ComputeMultiPlasticityStress(), and returnMap().

_dummy_pm

std::vector<Real> ComputeMultiPlasticityStress::_dummy_pm

protected

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

Definition at line 60 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_elastic_flexural_rigidity_tensor

const MaterialProperty<RankFourTensor>* const ComputeMultiPlasticityStress::_elastic_flexural_rigidity_tensor

protected

The Cosserat elastic flexural rigidity tensor.

Definition at line 157 of file ComputeMultiPlasticityStress.h.

_elastic_strain

MaterialProperty<RankTwoTensor>& ComputeGeneralStressBase::_elastic_strain

protectedinherited

Elastic strain material property.

Definition at line 50 of file ComputeGeneralStressBase.h.

Referenced by ComputeSmearedCrackingStress::computeCrackStrainAndOrientation(), computeQpStress(), FiniteStrainPlasticMaterial::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeLinearElasticStress::computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeMultipleInelasticStressBase::computeQpStressIntermediateConfiguration(), ComputeMultipleInelasticStressBase::finiteStrainRotation(), and ComputeGeneralStressBase::initQpStatefulProperties().

_elastic_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_elastic_strain_old

protected

Old value of elastic strain.

Definition at line 148 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_elasticity_tensor

const MaterialProperty<RankFourTensor>& ComputeMultiPlasticityStress::_elasticity_tensor

protected

Elasticity tensor material property.

Definition at line 100 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_elasticity_tensor_name

const std::string ComputeMultiPlasticityStress::_elasticity_tensor_name

protected

Name of the elasticity tensor material property.

Definition at line 98 of file ComputeMultiPlasticityStress.h.

_epp_tol

Real ComputeMultiPlasticityStress::_epp_tol

protected

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

Definition at line 57 of file ComputeMultiPlasticityStress.h.

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

_extra_stress

const MaterialProperty<RankTwoTensor>& ComputeGeneralStressBase::_extra_stress

protectedinherited

Extra stress tensor.

Definition at line 53 of file ComputeGeneralStressBase.h.

Referenced by ComputeGeneralStressBase::computeQpProperties().

_f

std::vector<const SolidMechanicsPlasticModel *> MultiPlasticityRawComponentAssembler::_f

protectedinherited

User objects that define the yield functions, flow potentials, etc.

Definition at line 71 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::activeModelSurfaces(), MultiPlasticityRawComponentAssembler::activeSurfaces(), MultiPlasticityRawComponentAssembler::anyActiveSurfaces(), applyKuhnTucker(), MultiPlasticityRawComponentAssembler::buildActiveConstraints(), MultiPlasticityRawComponentAssembler::buildActiveConstraintsJoint(), MultiPlasticityRawComponentAssembler::buildActiveConstraintsRock(), canAddConstraints(), checkAdmissible(), checkKuhnTucker(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), quickStep(), residual2(), returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

_fspb_debug

MooseEnum MultiPlasticityDebugger::_fspb_debug

protectedinherited

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

Definition at line 57 of file MultiPlasticityDebugger.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_fspb_debug_intnl

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl

protectedinherited

Debug the Jacobian entires at these internal parameters.

Definition at line 66 of file MultiPlasticityDebugger.h.

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

_fspb_debug_intnl_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change

protectedinherited

Debug finite-differencing parameters for the internal parameters.

Definition at line 75 of file MultiPlasticityDebugger.h.

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

_fspb_debug_pm

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm

protectedinherited

Debug the Jacobian entires at these plastic multipliers.

Definition at line 63 of file MultiPlasticityDebugger.h.

Referenced by MultiPlasticityDebugger::checkJacobian(), MultiPlasticityDebugger::checkSolution(), MultiPlasticityDebugger::outputAndCheckDebugParameters(), and plasticStep().

_fspb_debug_pm_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change

protectedinherited

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 72 of file MultiPlasticityDebugger.h.

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

_fspb_debug_stress

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress

protectedinherited

Debug the Jacobian entries at this stress.

Definition at line 60 of file MultiPlasticityDebugger.h.

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

_fspb_debug_stress_change

Real MultiPlasticityDebugger::_fspb_debug_stress_change

protectedinherited

Debug finite-differencing parameter for the stress.

Definition at line 69 of file MultiPlasticityDebugger.h.

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

_ignore_failures

bool ComputeMultiPlasticityStress::_ignore_failures

protected

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

Definition at line 46 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

_initial_stress_fcn

std::vector<const Function *> ComputeGeneralStressBase::_initial_stress_fcn

protectedinherited

initial stress components

Definition at line 56 of file ComputeGeneralStressBase.h.

_intnl

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

protected

internal parameters

Definition at line 109 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_intnl_old

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

protected

old values of internal parameters

Definition at line 112 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and plasticStep().

_iter

MaterialProperty<Real>& ComputeMultiPlasticityStress::_iter

protected

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

Definition at line 118 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_Jacobian_mult

MaterialProperty<RankFourTensor>& ComputeGeneralStressBase::_Jacobian_mult

protectedinherited

derivative of stress w.r.t. strain (_dstress_dstrain)

Definition at line 59 of file ComputeGeneralStressBase.h.

Referenced by ComputeStrainIncrementBasedStress::computeQpJacobian(), FiniteStrainHyperElasticViscoPlastic::computeQpJacobian(), ComputeMultipleInelasticStressBase::computeQpJacobianMult(), ComputeMultipleInelasticCosseratStress::computeQpJacobianMult(), computeQpStress(), FiniteStrainPlasticMaterial::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeLinearElasticPFFractureStress::computeQpStress(), ComputeLinearElasticStress::computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeMultipleInelasticStressBase::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeMultipleCrystalPlasticityStress::computeQpStress(), ComputeMultipleInelasticStressBase::computeQpStressIntermediateConfiguration(), ComputeLinearElasticPFFractureStress::computeStrainSpectral(), ComputeLinearElasticPFFractureStress::computeStrainVolDev(), ComputeLinearElasticPFFractureStress::computeStressSpectral(), FiniteStrainUObasedCP::elasticTangentModuli(), FiniteStrainUObasedCP::elastoPlasticTangentModuli(), ComputeMultipleInelasticStressBase::finiteStrainRotation(), postReturnMap(), FiniteStrainCrystalPlasticity::postSolveQp(), FiniteStrainCrystalPlasticity::preSolveQp(), and ComputeMultipleInelasticStressBase::updateQpStateSingleModel().

_Jacobian_mult_couple

MaterialProperty<RankFourTensor>* const ComputeMultiPlasticityStress::_Jacobian_mult_couple

protected

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

Definition at line 166 of file ComputeMultiPlasticityStress.h.

_ld_encountered

MaterialProperty<Real>& ComputeMultiPlasticityStress::_ld_encountered

protected

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

Definition at line 124 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_linesearch_needed

MaterialProperty<Real>& ComputeMultiPlasticityStress::_linesearch_needed

protected

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

Definition at line 121 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

_max_iter

unsigned int ComputeMultiPlasticityStress::_max_iter

protected

Maximum number of Newton-Raphson iterations allowed.

Definition at line 37 of file ComputeMultiPlasticityStress.h.

Referenced by returnMap().

_max_stepsize_for_dumb

Real ComputeMultiPlasticityStress::_max_stepsize_for_dumb

protected

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

Definition at line 43 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

_mechanical_strain

const MaterialProperty<RankTwoTensor>& ComputeGeneralStressBase::_mechanical_strain

protectedinherited

Mechanical strain material property.

Definition at line 46 of file ComputeGeneralStressBase.h.

Referenced by ComputeStrainIncrementBasedStress::computeQpStress(), FiniteStrainPlasticMaterial::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeLinearElasticPFFractureStress::computeQpStress(), ComputeLinearElasticStress::computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeLinearElasticPFFractureStress::computeStrainSpectral(), ComputeLinearElasticPFFractureStress::computeStrainVolDev(), and ComputeLinearElasticPFFractureStress::computeStressSpectral().

_min_f_tol

Real MultiPlasticityLinearSystem::_min_f_tol

protectedinherited

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

Definition at line 131 of file MultiPlasticityLinearSystem.h.

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

_min_stepsize

Real ComputeMultiPlasticityStress::_min_stepsize

protected

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

Definition at line 40 of file ComputeMultiPlasticityStress.h.

Referenced by plasticStep().

_model_given_surface

std::vector<unsigned int> MultiPlasticityRawComponentAssembler::_model_given_surface

privateinherited

given a surface number, this returns the model number

Definition at line 291 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::modelNumber(), and MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_model_surface_given_surface

std::vector<unsigned int> MultiPlasticityRawComponentAssembler::_model_surface_given_surface

privateinherited

given a surface number, this returns the corresponding-model's internal surface number

Definition at line 294 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_my_curvature

RankTwoTensor ComputeMultiPlasticityStress::_my_curvature

protected

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

Definition at line 178 of file ComputeMultiPlasticityStress.h.

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

_my_elasticity_tensor

RankFourTensor ComputeMultiPlasticityStress::_my_elasticity_tensor

protected

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

Definition at line 169 of file ComputeMultiPlasticityStress.h.

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

_my_flexural_rigidity_tensor

RankFourTensor ComputeMultiPlasticityStress::_my_flexural_rigidity_tensor

protected

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

Definition at line 175 of file ComputeMultiPlasticityStress.h.

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

_my_strain_increment

RankTwoTensor ComputeMultiPlasticityStress::_my_strain_increment

protected

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

Definition at line 172 of file ComputeMultiPlasticityStress.h.

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

_n

MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n

protected

current value of transverse direction

Definition at line 130 of file ComputeMultiPlasticityStress.h.

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

_n_input

RealVectorValue ComputeMultiPlasticityStress::_n_input

protected

the supplied transverse direction vector

Definition at line 89 of file ComputeMultiPlasticityStress.h.

Referenced by ComputeMultiPlasticityStress(), and initQpStatefulProperties().

_n_old

const MaterialProperty<RealVectorValue>& ComputeMultiPlasticityStress::_n_old

protected

old value of transverse direction

Definition at line 133 of file ComputeMultiPlasticityStress.h.

Referenced by postReturnMap().

_n_supplied

bool ComputeMultiPlasticityStress::_n_supplied

protected

User supplied the transverse direction vector.

Definition at line 86 of file ComputeMultiPlasticityStress.h.

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

_num_models

unsigned int MultiPlasticityRawComponentAssembler::_num_models

protectedinherited

Number of plastic models for this material.

Definition at line 53 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints(), MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), MultiPlasticityDebugger::checkSolution(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityDebugger::fddflowPotential_dintnl(), MultiPlasticityDebugger::fddyieldFunction_dintnl(), MultiPlasticityDebugger::fdJacobian(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), initQpStatefulProperties(), MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), MultiPlasticityLinearSystem::nrStep(), MultiPlasticityDebugger::outputAndCheckDebugParameters(), plasticStep(), residual2(), returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

_num_surfaces

unsigned int MultiPlasticityRawComponentAssembler::_num_surfaces

protectedinherited

Number of surfaces within the plastic models.

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

Definition at line 62 of file MultiPlasticityRawComponentAssembler.h.

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

_params

const InputParameters& MultiPlasticityRawComponentAssembler::_params

protectedinherited

Definition at line 50 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_perform_finite_strain_rotations

bool ComputeMultiPlasticityStress::_perform_finite_strain_rotations

protected

whether to perform the rotations necessary in finite-strain simulations

Definition at line 95 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_plastic_strain

MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain

protected

plastic strain

Definition at line 103 of file ComputeMultiPlasticityStress.h.

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

_plastic_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_plastic_strain_old

protected

Old value of plastic strain.

Definition at line 106 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_rot

RealTensorValue ComputeMultiPlasticityStress::_rot

protected

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

Definition at line 92 of file ComputeMultiPlasticityStress.h.

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

_rotation_increment

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_rotation_increment

protected

Rotation increment (coming from ComputeIncrementalStrain, for example)

Definition at line 142 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and postReturnMap().

_specialIC

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC

protectedinherited

Allows initial set of active constraints to be chosen optimally.

Definition at line 68 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints(), and MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_strain_increment

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_strain_increment

protected

strain increment (coming from ComputeIncrementalStrain, for example)

Definition at line 136 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_stress

MaterialProperty<RankTwoTensor>& ComputeGeneralStressBase::_stress

protectedinherited

Stress material property.

Definition at line 48 of file ComputeGeneralStressBase.h.

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

_stress_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_stress_old

protected

Old value of stress.

Definition at line 145 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress().

_surfaces_given_model

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

protectedinherited

_surfaces_given_model[model_number] = vector of surface numbers for this model

Definition at line 65 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::activeModelSurfaces(), MultiPlasticityRawComponentAssembler::activeSurfaces(), MultiPlasticityRawComponentAssembler::anyActiveSurfaces(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), quickStep(), and MultiPlasticityRawComponentAssembler::returnMapAll().

_svd_tol

Real MultiPlasticityLinearSystem::_svd_tol

protectedinherited

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

Definition at line 128 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

_tangent_operator_type

enum ComputeMultiPlasticityStress::TangentOperatorEnum ComputeMultiPlasticityStress::_tangent_operator_type

protected

Referenced by consistentTangentOperator(), and quickStep().

_total_strain_old

const MaterialProperty<RankTwoTensor>& ComputeMultiPlasticityStress::_total_strain_old

protected

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

Definition at line 139 of file ComputeMultiPlasticityStress.h.

_yf

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

protected

yield functions

Definition at line 115 of file ComputeMultiPlasticityStress.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

---

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

• solid_mechanics/include/materials/ComputeMultiPlasticityStress.h
• solid_mechanics/src/materials/ComputeMultiPlasticityStress.C