MultiParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates for plastic models where the yield function and flow directions depend on multiple parameters (called "stress_params" in the documentation and sp in the code) that are themselves functions of stress. More...

#include <MultiParameterPlasticityStressUpdate.h>

Inheritance diagram for MultiParameterPlasticityStressUpdate:

[legend]

Classes
struct	yieldAndFlow
	Struct designed to hold info about a single yield function and its derivatives, as well as the flow directions. More...

Public Member Functions
	MultiParameterPlasticityStressUpdate (const InputParameters &parameters, unsigned num_sp, unsigned num_yf, unsigned num_intnl)

void	setQp (unsigned int qp)
	Sets the value of the global variable _qp for inheriting classes. More...

virtual bool	requiresIsotropicTensor ()=0
	Does the model require the elasticity tensor to be isotropic? More...

virtual bool	isIsotropic ()
	Is the implmented model isotropic? The safe default is 'false'. More...

virtual Real	computeTimeStepLimit ()


void	resetQpProperties () final
	Retained as empty methods to avoid a warning from Material.C in framework. These methods are unused in all inheriting classes and should not be overwritten. More...

void	resetProperties () final

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual void	initQpStatefulProperties () override

virtual void	updateState (RankTwoTensor &strain_increment, RankTwoTensor &inelastic_strain_increment, const RankTwoTensor &rotation_increment, RankTwoTensor &stress_new, const RankTwoTensor &stress_old, const RankFourTensor &elasticity_tensor, const RankTwoTensor &elastic_strain_old, bool compute_full_tangent_operator, RankFourTensor &tangent_operator) override
	Given a strain increment that results in a trial stress, perform some procedure (such as an iterative return-mapping process) to produce an admissible stress, an elastic strain increment and an inelastic strain increment. More...

virtual void	propagateQpStatefulProperties () override
	If updateState is not called during a timestep, this will be. More...

virtual TangentCalculationMethod	getTangentCalculationMethod () override

Real	yieldF (const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
	Computes the smoothed yield function. More...

Real	yieldF (const std::vector< Real > &yfs) const
	Computes the smoothed yield function. More...

Real	ismoother (Real f_diff) const
	Smooths yield functions. More...

Real	smoother (Real f_diff) const
	Derivative of ismoother. More...

Real	dsmoother (Real f_diff) const
	Derivative of smoother. More...

yieldAndFlow	smoothAllQuantities (const std::vector< Real > &stress_params, const std::vector< Real > &intnl) const
	Calculates all yield functions and derivatives, and then performs the smoothing scheme. More...

int	lineSearch (Real &res2, std::vector< Real > &stress_params, Real &gaE, const std::vector< Real > &trial_stress_params, yieldAndFlow &smoothed_q, const std::vector< Real > &intnl_ok, std::vector< Real > &intnl, std::vector< Real > &rhs, Real &linesearch_needed) const
	Performs a line-search to find stress_params and gaE Upon entry: More...

int	nrStep (const yieldAndFlow &smoothed_q, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, const std::vector< Real > &intnl, Real gaE, std::vector< Real > &rhs) const
	Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution. More...

Real	calculateRes2 (const std::vector< Real > &rhs) const
	Calculates the residual-squared for the Newton-Raphson + line-search. More...

void	calculateRHS (const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, std::vector< Real > &rhs) const
	Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j] ... More...

void	dnRHSdVar (const yieldAndFlow &smoothed_q, const std::vector< std::vector< Real >> &dintnl, const std::vector< Real > &stress_params, Real gaE, std::vector< double > &jac) const
	Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving the linear system using LAPACK gsev. More...

virtual void	errorHandler (const std::string &message) const
	Performs any necessary cleaning-up, then throw MooseException(message) More...

virtual void	yieldFunctionValuesV (const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const =0
	Computes the values of the yield functions, given stress_params and intnl parameters. More...

virtual void	computeAllQV (const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const =0
	Completely fills all_q with correct values. More...

virtual void	preReturnMapV (const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl)
	Derived classes may employ this function to record stuff or do other computations prior to the return-mapping algorithm. More...

virtual void	initializeVarsV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const
	Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm. More...

virtual void	setIntnlValuesV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const =0
	Sets the internal parameters based on the trial values of stress_params, their current values, and the old values of the internal parameters. More...

virtual void	setIntnlDerivativesV (const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const =0
	Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters. More...

virtual void	computeStressParams (const RankTwoTensor &stress, std::vector< Real > &stress_params) const =0
	Computes stress_params, given stress. More...

virtual void	initializeReturnProcess ()
	Derived classes may use this to perform calculations before any return-map process is performed, for instance, to initialize variables. More...

virtual void	finalizeReturnProcess (const RankTwoTensor &rotation_increment)
	Derived classes may use this to perform calculations after the return-map process has completed successfully in stress_param space but before the returned stress tensor has been calculcated. More...

virtual void	setStressAfterReturnV (const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const =0
	Sets stress from the admissible parameters. More...

virtual void	setInelasticStrainIncrementAfterReturn (const RankTwoTensor &stress_trial, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &elasticity_tensor, const RankTwoTensor &returned_stress, RankTwoTensor &inelastic_strain_increment) const
	Sets inelastic strain increment from the returned configuration This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called. More...

virtual void	consistentTangentOperatorV (const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto)
	Calculates the consistent tangent operator. More...

virtual std::vector< RankTwoTensor >	dstress_param_dstress (const RankTwoTensor &stress) const =0
	d(stress_param[i])/d(stress) at given stress More...

virtual std::vector< RankFourTensor >	d2stress_param_dstress (const RankTwoTensor &stress) const =0
	d2(stress_param[i])/d(stress)/d(stress) at given stress More...

virtual void	setEffectiveElasticity (const RankFourTensor &Eijkl)=0
	Sets _Eij and _En and _Cij. More...

void	dVardTrial (bool elastic_only, const std::vector< Real > &trial_stress_params, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, Real step_size, bool compute_full_tangent_operator, std::vector< std::vector< Real >> &dvar_dtrial) const
	Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_params for the (sub)strain increment. More...

bool	precisionLoss (const std::vector< Real > &solution, const std::vector< Real > &stress_params, Real gaE) const
	Check whether precision loss has occurred. More...

Protected Attributes
const unsigned	_num_sp
	Number of stress parameters. More...

const std::vector< Real >	_definitely_ok_sp
	An admissible value of stress_params at the initial time. More...

std::vector< std::vector< Real > >	_Eij
	E[i, j] in the system of equations to be solved. More...

Real	_En
	normalising factor More...

std::vector< std::vector< Real > >	_Cij
	_Cij[i, j] * _Eij[j, k] = 1 iff j == k More...

const unsigned	_num_yf
	Number of yield functions. More...

const unsigned	_num_intnl
	Number of internal parameters. More...

const unsigned	_max_nr_its
	Maximum number of Newton-Raphson iterations allowed in the return-map process. More...

const bool	_perform_finite_strain_rotations
	Whether to perform finite-strain rotations. More...

const Real	_smoothing_tol
	Smoothing tolerance: edges of the yield surface get smoothed by this amount. More...

const Real	_smoothing_tol2
	Square of the smoothing tolerance. More...

const Real	_f_tol
	The yield-function tolerance. More...

const Real	_f_tol2
	Square of the yield-function tolerance. More...

const Real	_min_step_size
	In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value. More...

bool	_step_one
	handles case of initial_stress that is inadmissible being supplied More...

const bool	_warn_about_precision_loss
	Output a warning message if precision loss is encountered during the return-map process. More...

MaterialProperty< RankTwoTensor > &	_plastic_strain
	plastic strain More...

const MaterialProperty< RankTwoTensor > &	_plastic_strain_old
	Old value of plastic strain. More...

MaterialProperty< std::vector< Real > > &	_intnl
	internal parameters More...

const MaterialProperty< std::vector< Real > > &	_intnl_old
	old values of internal parameters More...

MaterialProperty< std::vector< Real > > &	_yf
	yield functions More...

MaterialProperty< Real > &	_iter
	Number of Newton-Raphson iterations used in the return-map. More...

MaterialProperty< Real > &	_max_iter_used
	Maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation. More...

const MaterialProperty< Real > &	_max_iter_used_old
	Old value of maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation. More...

MaterialProperty< Real > &	_linesearch_needed
	Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) More...

const std::string	_base_name
	Name used as a prefix for all material properties related to the stress update model. More...

Static Protected Attributes
constexpr static unsigned	_tensor_dimensionality = 3
	Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) More...

Private Types
enum	SmootherFunctionType { SmootherFunctionType::cos, SmootherFunctionType::poly1, SmootherFunctionType::poly2, SmootherFunctionType::poly3 }
	The type of smoother function. More...

Private Attributes
std::vector< Real >	_trial_sp
	"Trial" value of stress_params that initializes the return-map process This is derived from stress = stress_old + Eijkl * strain_increment. More...

RankTwoTensor	_stress_trial
	"Trial" value of stress that is set at the beginning of the return-map process. More...

std::vector< Real >	_rhs
	0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i] ... More...

std::vector< std::vector< Real > >	_dvar_dtrial
	d({stress_param[i], gaE})/d(trial_stress_param[j]) More...

std::vector< Real >	_ok_sp
	The state (ok_sp, ok_intnl) is known to be admissible, so ok_sp are stress_params that are "OK". More...

std::vector< Real >	_ok_intnl
	The state (ok_sp, ok_intnl) is known to be admissible. More...

std::vector< Real >	_del_stress_params
	_del_stress_params = trial_stress_params - ok_sp This is fixed at the beginning of the return-map process, irrespective of substepping. More...

std::vector< Real >	_current_sp
	The current values of the stress params during the Newton-Raphson. More...

std::vector< Real >	_current_intnl
	The current values of the internal params during the Newton-Raphson. More...

enum MultiParameterPlasticityStressUpdate::SmootherFunctionType	_smoother_function_type

Detailed Description

MultiParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates for plastic models where the yield function and flow directions depend on multiple parameters (called "stress_params" in the documentation and sp in the code) that are themselves functions of stress.

Let the stress_params be S = {S[0], S[1], ... S[N-1]} and define _num_sp = N.

For instance, CappedDruckerPrager plasticity has _num_sp = 2 and defines S[0] = p = stress.trace() S[1] = q = sqrt(stress.secondInvariant())

The point of the stress_params is to reduce the number of degrees of freedom involved in the return-map algorithm (for computational efficiency as all the intensive computation occurs in "stress_param space") and to allow users to write their yield functions, etc, in forms that are clear to them to avoid errors.

Derived classes can describe plasticity via multiple yield functions. The number of yield function is _num_yf. The yield function(s) and flow potential(s) must be functions of the stress_params.

Derived classes can use any number of internal parameters. The number of internal parameters is _num_intnl. The derived classes must define the evolution of these internal parameters, which are typically functions of certain "distances" in the return-mapping procedure.

For instance, CappedDruckerPrager plasticity has three yield functions (a smoothed DruckerPrager cone, a tensile failure envelope, and a compressive failure envelope) and two internal parameters (plastic shear strain and plastic tensile strain). The strength parameters (cohesion, friction angle, dilation angle, tensile strength, compressive strength) are functions of these internal parameters.

A novel smoothing procedure is used to smooth the derived class's yield functions into a single, smooth yield function. The smoothing is also performed for the flow potential. All return-mapping, etc, processes are performed on this single yield function and flow potential.

The return-map algorithm implements the updateState method of StressUpdateBase. In particular, the following system of equations is solved: 0 = _rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = _rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j] ... 0 = _rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j] 0 = _rhs[N] = f(S, intnl) (here there is an implied sum over j) for the _num_sp stress_params S, and the single scalar ga (actually I solve for gaE = ga * _En, where _En is a normalising factor so that gaE is of similar magnitude to the S variables). In general E[a, b] = dS[a]/dstress(i, j) E(i, j, k, l) dS[b]/dstress(k, l), and in the code below it is assumed that the E[a, b] are independent of the stress_params, but adding a dependence would only require some Jacobian terms to be modified. There are N+1 rhs equations to solve. Here f is the smoothed yield function, so the last equation is the admissibility condition (the returned stress lies on the yield surface) and g is the flow potential so the other conditions implement the normality condition. It is up to the derived classes to defined E[i, j] so that the rhs really represents the normality condition.

For instance, CappedDruckerPrager has (with Eijkl being the elasticity tensor): E[0, 0] = _Epp = Eijkl.sum3x3() E[0, 1] = 0 E[1, 0] = 0 E[1, 1] = Eijkl(0, 1, 0, 1)

Definition at line 99 of file MultiParameterPlasticityStressUpdate.h.

Member Enumeration Documentation

◆ SmootherFunctionType

enum MultiParameterPlasticityStressUpdate::SmootherFunctionType

strongprivate

The type of smoother function.

Enumerator
cos
poly1
poly2
poly3

Definition at line 743 of file MultiParameterPlasticityStressUpdate.h.

   {
     cos,
     poly1,
     poly2,
     poly3
   } _smoother_function_type;

Constructor & Destructor Documentation

◆ MultiParameterPlasticityStressUpdate()

MultiParameterPlasticityStressUpdate::MultiParameterPlasticityStressUpdate	(	const InputParameters &	parameters,
		unsigned	num_sp,
		unsigned	num_yf,
		unsigned	num_intnl
	)

Definition at line 82 of file MultiParameterPlasticityStressUpdate.C.

   : StressUpdateBase(parameters),
     _num_sp(num_sp),
     _definitely_ok_sp(isParamValid("admissible_stress")
                           ? getParam<std::vector<Real>>("admissible_stress")
                           : std::vector<Real>(_num_sp, 0.0)),
     _Eij(num_sp, std::vector<Real>(num_sp)),
     _En(1.0),
     _Cij(num_sp, std::vector<Real>(num_sp)),
     _num_yf(num_yf),
     _num_intnl(num_intnl),
     _max_nr_its(getParam<unsigned>("max_NR_iterations")),
     _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
     _smoothing_tol(getParam<Real>("smoothing_tol")),
     _smoothing_tol2(Utility::pow<2>(getParam<Real>("smoothing_tol"))),
     _f_tol(getParam<Real>("yield_function_tol")),
     _f_tol2(Utility::pow<2>(getParam<Real>("yield_function_tol"))),
     _min_step_size(getParam<Real>("min_step_size")),
     _step_one(declareRestartableData<bool>("step_one", true)),
     _warn_about_precision_loss(getParam<bool>("warn_about_precision_loss")),
  
     _plastic_strain(declareProperty<RankTwoTensor>(_base_name + "plastic_strain")),
     _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "plastic_strain")),
     _intnl(declareProperty<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
     _intnl_old(
         getMaterialPropertyOld<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
     _yf(declareProperty<std::vector<Real>>(_base_name + "plastic_yield_function")),
     _iter(declareProperty<Real>(_base_name +
                                 "plastic_NR_iterations")), // this is really an unsigned int, but
                                                            // for visualisation i convert it to Real
     _max_iter_used(declareProperty<Real>(
         _base_name + "max_plastic_NR_iterations")), // this is really an unsigned int, but
                                                     // for visualisation i convert it to Real
     _max_iter_used_old(getMaterialPropertyOld<Real>(_base_name + "max_plastic_NR_iterations")),
     _linesearch_needed(
         declareProperty<Real>(_base_name + "plastic_linesearch_needed")), // this is really a
                                                                           // boolean, but for
                                                                           // visualisation i
                                                                           // convert it to Real
     _trial_sp(num_sp),
     _stress_trial(RankTwoTensor()),
     _rhs(num_sp + 1),
     _dvar_dtrial(num_sp + 1, std::vector<Real>(num_sp, 0.0)),
     _ok_sp(num_sp),
     _ok_intnl(num_intnl),
     _del_stress_params(num_sp),
     _current_sp(num_sp),
     _current_intnl(num_intnl),
     _smoother_function_type(
         parameters.get<MooseEnum>("smoother_function_type").getEnum<SmootherFunctionType>())
 {
   if (_definitely_ok_sp.size() != _num_sp)
     mooseError("MultiParameterPlasticityStressUpdate: admissible_stress parameter must consist of ",
                _num_sp,
                " numbers");
 }

Member Function Documentation

◆ calculateRes2()

Real MultiParameterPlasticityStressUpdate::calculateRes2 ( const std::vector< Real > & rhs ) const

protected

Calculates the residual-squared for the Newton-Raphson + line-search.

Parameters

rhs[in] The RHS vector

Returns: sum_i (rhs[i] * rhs[i])

Definition at line 802 of file MultiParameterPlasticityStressUpdate.C.

 {
   Real res2 = 0.0;
   for (const auto & r : rhs)
     res2 += r * r;
   return res2;
 }

Referenced by lineSearch(), and updateState().

◆ calculateRHS()

void MultiParameterPlasticityStressUpdate::calculateRHS	(	const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	stress_params,
		Real	gaE,
		const yieldAndFlow &	smoothed_q,
		std::vector< Real > &	rhs
	)		const

protected

Calculates the RHS in the following 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j] ...

0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j] 0 = rhs[N] = f(S, intnl) where N = _num_sp

Parameters

trial_stress_params[in]	The trial stress parameters for this (sub)strain increment, S[:]^trial
stress_params[in]	The current stress parameters, S[:]
gaE[in]	ga*_En (the normalisation with _En is so that gaE is of similar magnitude to S)
smoothed_q[in]	Holds the current value of yield function and derivatives evaluated at the current stress parameters and the current internal parameters
rhs[out]	The result

Definition at line 811 of file MultiParameterPlasticityStressUpdate.C.

 {
   const Real ga = gaE / _En;
   for (unsigned i = 0; i < _num_sp; ++i)
   {
     rhs[i] = stress_params[i] - trial_stress_params[i];
     for (unsigned j = 0; j < _num_sp; ++j)
       rhs[i] += ga * _Eij[i][j] * smoothed_q.dg[j];
   }
   rhs[_num_sp] = smoothed_q.f;
 }

Referenced by lineSearch(), and updateState().

◆ computeAllQV()

virtual void MultiParameterPlasticityStressUpdate::computeAllQV	(	const std::vector< Real > &	stress_params,
		const std::vector< Real > &	intnl,
		std::vector< yieldAndFlow > &	all_q
	)		const

protectedpure virtual

Completely fills all_q with correct values.

These values are: (1) the yield function values, yf[i] (2) d(yf[i])/d(stress_params[j]) (3) d(yf[i])/d(intnl[j]) (4) d(flowPotential[i])/d(stress_params[j]) (5) d2(flowPotential[i])/d(stress_params[j])/d(stress_params[k]) (6) d2(flowPotential[i])/d(stress_params[j])/d(intnl[k])

Parameters

	stress_params[in]	The stress parameters
	intnl[in]	The internal parameters
[out]	all_q	All the desired quantities

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by smoothAllQuantities().

◆ computeStressParams()

virtual void MultiParameterPlasticityStressUpdate::computeStressParams	(	const RankTwoTensor &	stress,
		std::vector< Real > &	stress_params
	)		const

protectedpure virtual

Computes stress_params, given stress.

Derived classes must override this

Parameters

stress[in]	Stress tensor
stress_params[out]	The compute stress_params

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by updateState().

◆ computeTimeStepLimit()

Real StressUpdateBase::computeTimeStepLimit ( )

virtualinherited

Reimplemented in RadialReturnStressUpdate.

Definition at line 56 of file StressUpdateBase.C.

 {
   return std::numeric_limits<Real>::max();
 }

◆ consistentTangentOperatorV()

void MultiParameterPlasticityStressUpdate::consistentTangentOperatorV	(	const RankTwoTensor &	stress_trial,
		const std::vector< Real > &	trial_stress_params,
		const RankTwoTensor &	stress,
		const std::vector< Real > &	stress_params,
		Real	gaE,
		const yieldAndFlow &	smoothed_q,
		const RankFourTensor &	Eijkl,
		bool	compute_full_tangent_operator,
		const std::vector< std::vector< Real >> &	dvar_dtrial,
		RankFourTensor &	cto
	)

protectedvirtual

Calculates the consistent tangent operator.

Derived classes may choose to override this for computational efficiency. The implementation in this class is quite expensive, even though it looks compact and clean, because of all the manipulations of RankFourTensors involved.

Parameters

	stress_trial[in]	the trial value of the stress tensor for this strain increment
	trial_stress_params[in]	the trial values of the stress_params for this strain increment
	stress[in]	the returned value of the stress tensor for this strain increment
	stress_params[in]	the returned value of the stress_params for this strain increment
	gaE[in]	the total value of that came from this strain increment
	smoothed_q[in]	contains the yield function and derivatives evaluated at (p, q)
	Eijkl[in]	The elasticity tensor
	compute_full_tangent_operator[in]	true if the full consistent tangent operator is needed, otherwise false
	dvar_dtrial[in]	dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j]) for this strain increment
[out]	cto	The consistent tangent operator

Reimplemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, TensileStressUpdate, and CappedMohrCoulombCosseratStressUpdate.

Definition at line 725 of file MultiParameterPlasticityStressUpdate.C.

 {
   cto = elasticity_tensor;
   if (!compute_full_tangent_operator)
     return;
  
   const Real ga = gaE / _En;
  
   const std::vector<RankTwoTensor> dsp = dstress_param_dstress(stress);
   const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
  
   for (unsigned a = 0; a < _num_sp; ++a)
   {
     for (unsigned b = 0; b < _num_sp; ++b)
     {
       const RankTwoTensor t = elasticity_tensor * dsp_trial[a];
       RankTwoTensor s = _Cij[b][a] * dsp[b];
       for (unsigned c = 0; c < _num_sp; ++c)
         s -= dsp[b] * _Cij[b][c] * dvar_dtrial[c][a];
       s = elasticity_tensor * s;
       cto -= s.outerProduct(t);
     }
   }
  
   const std::vector<RankFourTensor> d2sp = d2stress_param_dstress(stress);
   RankFourTensor Tijab;
   for (unsigned i = 0; i < _num_sp; ++i)
     Tijab += smoothed_q.dg[i] * d2sp[i];
   Tijab = ga * elasticity_tensor * Tijab;
  
   RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
   try
   {
     inv = inv.transposeMajor().invSymm();
   }
   catch (const MooseException & e)
   {
     // Cannot form the inverse, so probably at some degenerate place in stress space.
     // Just return with the "best estimate" of the cto.
     mooseWarning("MultiParameterPlasticityStressUpdate: Cannot invert 1+T in consistent tangent "
                  "operator computation at quadpoint ",
                  _qp,
                  " of element ",
                  _current_elem->id());
     return;
   }
  
   cto = (cto.transposeMajor() * inv).transposeMajor();
 }

Referenced by updateState().

◆ d2stress_param_dstress()

virtual std::vector<RankFourTensor> MultiParameterPlasticityStressUpdate::d2stress_param_dstress ( const RankTwoTensor & stress ) const

protectedpure virtual

d2(stress_param[i])/d(stress)/d(stress) at given stress

Parameters

stress stress tensor

Returns: d2(stress_param[:])/d(stress)/d(stress)

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by consistentTangentOperatorV().

◆ dnRHSdVar()

void MultiParameterPlasticityStressUpdate::dnRHSdVar	(	const yieldAndFlow &	smoothed_q,
		const std::vector< std::vector< Real >> &	dintnl,
		const std::vector< Real > &	stress_params,
		Real	gaE,
		std::vector< double > &	jac
	)		const

protected

Derivative of -RHS with respect to the stress_params and gaE, placed into an array ready for solving the linear system using LAPACK gsev.

Parameters

smoothed_q[in]	Holds the current value of yield function and derivatives evaluated at the current values of the stress_params and the internal parameters
dintnl[in]	The derivatives of the internal parameters wrt the stress_params
stress_params[in]	The current value of the stress_params during the Newton-Raphson process
gaE[in]	The current value of gaE
jac[out]	The outputted derivatives

Definition at line 828 of file MultiParameterPlasticityStressUpdate.C.

 {
   for (auto & jac_entry : jac)
     jac_entry = 0.0;
  
   const Real ga = gaE / _En;
  
   unsigned ind = 0;
   for (unsigned var = 0; var < _num_sp; ++var)
   {
     for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
     {
       if (var == rhs)
         jac[ind] -= 1.0;
       for (unsigned j = 0; j < _num_sp; ++j)
       {
         jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g[j][var];
         for (unsigned k = 0; k < _num_intnl; ++k)
           jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g_di[j][k] * dintnl[k][var];
       }
       ind++;
     }
     // now rhs = _num_sp (that is, the yield function)
     jac[ind] -= smoothed_q.df[var];
     for (unsigned k = 0; k < _num_intnl; ++k)
       jac[ind] -= smoothed_q.df_di[k] * dintnl[k][var];
     ind++;
   }
  
   // now var = _num_sp (that is, gaE)
   for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
   {
     for (unsigned j = 0; j < _num_sp; ++j)
       jac[ind] -= (1.0 / _En) * _Eij[rhs][j] * smoothed_q.dg[j];
     ind++;
   }
   // now rhs = _num_sp (that is, the yield function)
   jac[ind] = 0.0;
 }

Referenced by dVardTrial(), and nrStep().

◆ dsmoother()

Real MultiParameterPlasticityStressUpdate::dsmoother ( Real f_diff ) const

protected

Derivative of smoother.

Definition at line 533 of file MultiParameterPlasticityStressUpdate.C.

 {
   if (std::abs(f_diff) >= _smoothing_tol)
     return 0.0;
   switch (_smoother_function_type)
   {
     case SmootherFunctionType::cos:
       return 0.25 * M_PI / _smoothing_tol * std::cos(f_diff * M_PI * 0.5 / _smoothing_tol);
     case SmootherFunctionType::poly1:
       return 0.75 / _smoothing_tol * (1.0 - Utility::pow<2>(f_diff / _smoothing_tol));
     case SmootherFunctionType::poly2:
       return 0.625 / _smoothing_tol * (1.0 - Utility::pow<4>(f_diff / _smoothing_tol));
     case SmootherFunctionType::poly3:
       return (7.0 / 12.0 / _smoothing_tol) * (1.0 - Utility::pow<6>(f_diff / _smoothing_tol));
     default:
       return 0.0;
   }
 }

Referenced by smoothAllQuantities().

◆ dstress_param_dstress()

virtual std::vector<RankTwoTensor> MultiParameterPlasticityStressUpdate::dstress_param_dstress ( const RankTwoTensor & stress ) const

protectedpure virtual

d(stress_param[i])/d(stress) at given stress

Parameters

stress stress tensor

Returns: d(stress_param[:])/d(stress)

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by consistentTangentOperatorV(), and setInelasticStrainIncrementAfterReturn().

◆ dVardTrial()

void MultiParameterPlasticityStressUpdate::dVardTrial	(	bool	elastic_only,
		const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	stress_params,
		Real	gaE,
		const std::vector< Real > &	intnl,
		const yieldAndFlow &	smoothed_q,
		Real	step_size,
		bool	compute_full_tangent_operator,
		std::vector< std::vector< Real >> &	dvar_dtrial
	)		const

protected

Calculates derivatives of the stress_params and gaE with repect to the trial values of the stress_params for the (sub)strain increment.

After the strain increment has been fully applied, dvar_dtrial will contain the result appropriate to the full strain increment. Before that time (if applying in sub-strain increments) it will contain the result appropriate to the amount of strain increment applied successfully.

Parameters

elastic_only[in]	whether this was an elastic step: if so then the updates to dvar_dtrial are fairly trivial
trial_stress_params[in]	Trial values of stress_params for this (sub)strain increment
stress_params[in]	Returned values of stress_params for this (sub)strain increment
gaE[in]	the value of gaE that came from this (sub)strain increment
intnl[in]	the value of the internal parameters at the returned position
smoothed_q[in]	contains the yield function and derivatives evaluated at (stress_params, intnl)
step_size[in]	size of this (sub)strain increment
compute_full_tangent_operator[in]	true if the full consistent tangent operator is needed, otherwise false
dvar_dtrial[out]	dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])

Definition at line 873 of file MultiParameterPlasticityStressUpdate.C.

 {
   if (!_fe_problem.currentlyComputingJacobian())
     return;
  
   if (!compute_full_tangent_operator)
     return;
  
   if (elastic_only)
   {
     // no change to gaE, and all off-diag stuff remains unchanged from previous step
     for (unsigned v = 0; v < _num_sp; ++v)
       dvar_dtrial[v][v] += step_size;
     return;
   }
  
   const Real ga = gaE / _En;
  
   std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
   setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
  
   // rhs is described elsewhere, the following are changes in rhs wrt the trial_stress_param
   // values
   // In the following we use d(intnl)/d(trial variable) = - d(intnl)/d(variable)
   std::vector<Real> rhs_cto((_num_sp + 1) * _num_sp);
  
   unsigned ind = 0;
   for (unsigned a = 0; a < _num_sp; ++a)
   {
     // change in RHS[b] wrt changes in stress_param_trial[a]
     for (unsigned b = 0; b < _num_sp; ++b)
     {
       if (a == b)
         rhs_cto[ind] -= 1.0;
       for (unsigned j = 0; j < _num_sp; ++j)
         for (unsigned k = 0; k < _num_intnl; ++k)
           rhs_cto[ind] -= ga * _Eij[b][j] * smoothed_q.d2g_di[j][k] * dintnl[k][a];
       ind++;
     }
     // now b = _num_sp (that is, the yield function)
     for (unsigned k = 0; k < _num_intnl; ++k)
       rhs_cto[ind] -= smoothed_q.df_di[k] * dintnl[k][a];
     ind++;
   }
  
   // jac = d(-rhs)/d(var)
   std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
   dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
  
   std::vector<int> ipiv(_num_sp + 1);
   int info;
   const int gesv_num_rhs = _num_sp + 1;
   const int gesv_num_pq = _num_sp;
   LAPACKgesv_(&gesv_num_rhs,
               &gesv_num_pq,
               &jac[0],
               &gesv_num_rhs,
               &ipiv[0],
               &rhs_cto[0],
               &gesv_num_rhs,
               &info);
   if (info != 0)
     errorHandler("MultiParameterPlasticityStressUpdate: PETSC LAPACK gsev routine returned with "
                  "error code " +
                  Moose::stringify(info));
  
   ind = 0;
   std::vector<std::vector<Real>> dvarn_dtrialn(_num_sp + 1, std::vector<Real>(_num_sp, 0.0));
   for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
   {
     for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
     {
       dvarn_dtrialn[v][spt] = rhs_cto[ind];
       ind++;
     }
     // the final NR variable is gaE
     dvarn_dtrialn[_num_sp][spt] = rhs_cto[ind];
     ind++;
   }
  
   const std::vector<std::vector<Real>> dvar_dtrial_old = dvar_dtrial;
  
   for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
   {
     for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
     {
       dvar_dtrial[v][spt] = step_size * dvarn_dtrialn[v][spt];
       for (unsigned a = 0; a < _num_sp; ++a)
         dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
     }
   }
   // for gaE the formulae are a little different
   const unsigned v = _num_sp;
   for (unsigned spt = 0; spt < _num_sp; ++spt)
   {
     dvar_dtrial[v][spt] += step_size * dvarn_dtrialn[v][spt]; // note +=
     for (unsigned a = 0; a < _num_sp; ++a)
       dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
   }
 }

Referenced by updateState().

◆ errorHandler()

void MultiParameterPlasticityStressUpdate::errorHandler ( const std::string & message ) const

protectedvirtual

Performs any necessary cleaning-up, then throw MooseException(message)

Parameters

message The message to using in MooseException

Definition at line 661 of file MultiParameterPlasticityStressUpdate.C.

 {
   throw MooseException(message);
 }

Referenced by dVardTrial(), and updateState().

◆ finalizeReturnProcess()

void MultiParameterPlasticityStressUpdate::finalizeReturnProcess ( const RankTwoTensor & rotation_increment )

protectedvirtual

Derived classes may use this to perform calculations after the return-map process has completed successfully in stress_param space but before the returned stress tensor has been calculcated.

Parameters

rotation_increment[in] The large-strain rotation increment

Reimplemented in CappedDruckerPragerStressUpdate, CappedWeakPlaneStressUpdate, and CappedWeakInclinedPlaneStressUpdate.

Definition at line 672 of file MultiParameterPlasticityStressUpdate.C.

674 {

675 }

Referenced by updateState().

◆ getTangentCalculationMethod()

virtual TangentCalculationMethod MultiParameterPlasticityStressUpdate::getTangentCalculationMethod ( )

inlineoverrideprotectedvirtual

Reimplemented from StressUpdateBase.

Definition at line 123 of file MultiParameterPlasticityStressUpdate.h.

   {
     return TangentCalculationMethod::FULL;
   }

◆ initializeReturnProcess()

void MultiParameterPlasticityStressUpdate::initializeReturnProcess ( )

protectedvirtual

Derived classes may use this to perform calculations before any return-map process is performed, for instance, to initialize variables.

This is called at the very start of updateState, even before any checking for admissible stresses, etc, is performed

Reimplemented in CappedDruckerPragerStressUpdate, CappedWeakPlaneStressUpdate, and CappedWeakInclinedPlaneStressUpdate.

Definition at line 667 of file MultiParameterPlasticityStressUpdate.C.

668 {

669 }

Referenced by updateState().

◆ initializeVarsV()

void MultiParameterPlasticityStressUpdate::initializeVarsV	(	const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	intnl_old,
		std::vector< Real > &	stress_params,
		Real &	gaE,
		std::vector< Real > &	intnl
	)		const

protectedvirtual

Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.

The purpose of these "good guesses" is to speed up the Newton-Raphson process by providing it with a good initial guess. Derived classes may choose to override this if their plastic models are easy enough to solve approximately. The default values, provided by this class, are simply gaE = 0, stress_params = trial_stress_params, that is, the "good guess" is just the trial point for this (sub)strain increment.

Parameters

trial_stress_params[in]	The stress_params at the trial point
intnl_old[in]	The internal parameters before applying the (sub)strain increment
stress_params[out]	The "good guess" value of the stress_params
gaE[out]	The "good guess" value of gaE
intnl[out]	The "good guess" value of the internal parameters

Reimplemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Definition at line 713 of file MultiParameterPlasticityStressUpdate.C.

 {
   gaE = 0.0;
   std::copy(trial_stress_params.begin(), trial_stress_params.end(), stress_params.begin());
   std::copy(intnl_old.begin(), intnl_old.end(), intnl.begin());
 }

Referenced by updateState().

◆ initQpStatefulProperties()

void MultiParameterPlasticityStressUpdate::initQpStatefulProperties ( )

overrideprotectedvirtual

Reimplemented in CappedWeakInclinedPlaneStressUpdate.

Definition at line 141 of file MultiParameterPlasticityStressUpdate.C.

 {
   _plastic_strain[_qp].zero();
   _intnl[_qp].assign(_num_intnl, 0);
   _yf[_qp].assign(_num_yf, 0);
   _iter[_qp] = 0.0;
   _max_iter_used[_qp] = 0.0;
   _linesearch_needed[_qp] = 0.0;
 }

Referenced by CappedWeakInclinedPlaneStressUpdate::initQpStatefulProperties(), and updateState().

◆ isIsotropic()

virtual bool StressUpdateBase::isIsotropic ( )

inlinevirtualinherited

Is the implmented model isotropic? The safe default is 'false'.

Reimplemented in RadialReturnStressUpdate, CappedDruckerPragerStressUpdate, and CappedMohrCoulombStressUpdate.

Definition at line 112 of file StressUpdateBase.h.

112 { return false; };

Referenced by ComputeMultipleInelasticStress::initialSetup().

◆ ismoother()

Real MultiParameterPlasticityStressUpdate::ismoother ( Real f_diff ) const

protected

Smooths yield functions.

The returned value must be zero if abs(f_diff) >= _smoothing_tol and otherwise must satisfy, over -_smoothing_tol <= f_diff <= _smoothing_tol: (1) C2 (2) zero at f_diff = +/- _smoothing_tol (3) derivative is +/-0.5 at f_diff = +/- _smoothing_tol (4) derivative must be in [-0.5, 0.5] (5) second derivative is zero at f_diff = +/- _smoothing_tol (6) second derivative must be non-negative in order to ensure C2 differentiability and convexity of the smoothed yield surface.

Definition at line 484 of file MultiParameterPlasticityStressUpdate.C.

 {
   if (std::abs(f_diff) >= _smoothing_tol)
     return 0.0;
   switch (_smoother_function_type)
   {
     case SmootherFunctionType::cos:
       return -_smoothing_tol / M_PI * std::cos(0.5 * M_PI * f_diff / _smoothing_tol);
     case SmootherFunctionType::poly1:
       return 0.75 / _smoothing_tol *
              (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
               (_smoothing_tol2 / 12.0) * (Utility::pow<4>(f_diff / _smoothing_tol) - 1.0));
     case SmootherFunctionType::poly2:
       return 0.625 / _smoothing_tol *
              (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
               (_smoothing_tol2 / 30.0) * (Utility::pow<6>(f_diff / _smoothing_tol) - 1.0));
     case SmootherFunctionType::poly3:
       return (7.0 / 12.0 / _smoothing_tol) *
              (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
               (_smoothing_tol2 / 56.0) * (Utility::pow<8>(f_diff / _smoothing_tol) - 1.0));
     default:
       return 0.0;
   }
 }

Referenced by smoothAllQuantities(), and yieldF().

◆ lineSearch()

int MultiParameterPlasticityStressUpdate::lineSearch	(	Real &	res2,
		std::vector< Real > &	stress_params,
		Real &	gaE,
		const std::vector< Real > &	trial_stress_params,
		yieldAndFlow &	smoothed_q,
		const std::vector< Real > &	intnl_ok,
		std::vector< Real > &	intnl,
		std::vector< Real > &	rhs,
		Real &	linesearch_needed
	)		const

protected

Performs a line-search to find stress_params and gaE Upon entry:

rhs contains the solution to the Newton-Raphson (ie nrStep should have been called). If a full Newton step is used then stress_params[:] += rhs[0:_num_sp-1] and gaE += rhs[_num_sp]
res2 contains the residual-squared before applying any of solution
stress_params contains the stress_params before applying any of the solution
gaE contains gaE before applying any of the solution (that is contained in rhs) Upon exit:
stress_params will be the stress_params after applying the solution
gaE will be the stress_params after applying the solution
rhs will contain the updated rhs values (after applying the solution) ready for the next Newton-Raphson step,
res2 will be the residual-squared after applying the solution
intnl will contain the internal variables corresponding to the return from trial_stress_params to stress_params (and starting from intnl_ok)
linesearch_needed will be 1.0 if a linesearch was needed

smoothed_q will contain the value of the yield function and its derivatives, etc, at (stress_params, intnl)

Parameters

res2[in,out]	the residual-squared, both as an input and output
stress_params[in,out]	Upon input the value of the stress_params before the current Newton-Raphson process was initiated. Upon exit this will hold the values coming from the line search.
trial_stress_params[in]	Trial value for the stress_params for this (sub)strain increment
gaE[in,out]	Upon input the value of gaE before the current Newton-Raphson iteration was initiated. Upon exit this will hold the value coming from the line-search
smoothed_q[in,out]	Upon input, the value of the smoothed yield function and derivatives at the prior-to-Newton configuration. Upon exit this is evaluated at the new (stress_params, intnl)
intnl_ok[in]	The value of the internal parameters from the start of this (sub)strain increment
intnl[in,out]	The value of the internal parameters after the line-search has converged
rhs[in,out]	Upon entry this contains the solution to the Newton-Raphson. Upon exit this contains the updated rhs values

Returns: 0 if successful, 1 otherwise

Definition at line 553 of file MultiParameterPlasticityStressUpdate.C.

 {
   const Real res2_old = res2;
   const std::vector<Real> sp_params_old = stress_params;
   const Real gaE_old = gaE;
   const std::vector<Real> delta_nr_params = rhs;
  
   Real lam = 1.0;                     // line-search parameter
   const Real lam_min = 1E-10;         // minimum value of lam allowed
   const Real slope = -2.0 * res2_old; // "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 = res2_old; // 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
  
   while (true)
   {
     // update variables using the current line-search parameter
     for (unsigned i = 0; i < _num_sp; ++i)
       stress_params[i] = sp_params_old[i] + lam * delta_nr_params[i];
     gaE = gaE_old + lam * delta_nr_params[_num_sp];
  
     // and internal parameters
     setIntnlValuesV(trial_stress_params, stress_params, intnl_ok, intnl);
  
     smoothed_q = smoothAllQuantities(stress_params, intnl);
  
     // update rhs for next-time through
     calculateRHS(trial_stress_params, stress_params, gaE, smoothed_q, rhs);
     res2 = calculateRes2(rhs);
  
     // do the line-search
     if (res2 < res2_old + 1E-4 * lam * slope)
       break;
     else if (lam < lam_min)
       return 1;
     else if (lam == 1.0)
     {
       // model as a quadratic
       tmp_lam = -0.5 * slope / (res2 - res2_old - slope);
     }
     else
     {
       // model as a cubic
       const Real rhs1 = res2 - res2_old - lam * slope;
       const Real rhs2 = f2 - res2_old - lam2 * slope;
       const Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
       const Real b =
           (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
       if (a == 0.0)
         tmp_lam = -slope / (2.0 * b);
       else
       {
         const Real disc = Utility::pow<2>(b) - 3.0 * 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 = res2;
     lam = std::max(tmp_lam, 0.1 * lam);
   }
  
   if (lam < 1.0)
     linesearch_needed = 1.0;
   return 0;
 }

Referenced by updateState().

◆ nrStep()

int MultiParameterPlasticityStressUpdate::nrStep	(	const yieldAndFlow &	smoothed_q,
		const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	stress_params,
		const std::vector< Real > &	intnl,
		Real	gaE,
		std::vector< Real > &	rhs
	)		const

protected

Performs a Newton-Raphson step to attempt to zero rhs Upon return, rhs will contain the solution.

Parameters

smoothed_q[in]	The value of the smoothed yield function and derivatives prior to this Newton-Raphson step
trial_stress_params[in]	Trial value for the stress_params for this (sub)strain increment
stress_params[in]	The current value of the stress_params
intnl[in]	The current value of the internal parameters
gaE[in]	The current value of gaE
rhs[in,out]	Upon entry, the rhs to zero using Newton-Raphson. Upon exit, the solution to the Newton-Raphson problem

Returns: 0 if successful, 1 otherwise

Definition at line 637 of file MultiParameterPlasticityStressUpdate.C.

 {
   std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
   setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
  
   std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
   dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
  
   // use LAPACK to solve the linear system
   const int nrhs = 1;
   std::vector<int> ipiv(_num_sp + 1);
   int info;
   const int gesv_num_rhs = _num_sp + 1;
   LAPACKgesv_(
       &gesv_num_rhs, &nrhs, &jac[0], &gesv_num_rhs, &ipiv[0], &rhs[0], &gesv_num_rhs, &info);
   return info;
 }

Referenced by updateState().

◆ precisionLoss()

bool MultiParameterPlasticityStressUpdate::precisionLoss	(	const std::vector< Real > &	solution,
		const std::vector< Real > &	stress_params,
		Real	gaE
	)		const

protected

Check whether precision loss has occurred.

Parameters

[in]	solution	The solution to the Newton-Raphson system
[in]	stress_params	The currect values of the stress_params for this (sub)strain increment
[in]	gaE	The currenct value of gaE for this (sub)strain increment

Returns: true if precision loss has occurred

Definition at line 983 of file MultiParameterPlasticityStressUpdate.C.

 {
   if (std::abs(solution[_num_sp]) > 1E-13 * std::abs(gaE))
     return false;
   for (unsigned i = 0; i < _num_sp; ++i)
     if (std::abs(solution[i]) > 1E-13 * std::abs(stress_params[i]))
       return false;
   return true;
 }

Referenced by updateState().

◆ preReturnMapV()

void MultiParameterPlasticityStressUpdate::preReturnMapV	(	const std::vector< Real > &	trial_stress_params,
		const RankTwoTensor &	stress_trial,
		const std::vector< Real > &	intnl_old,
		const std::vector< Real > &	yf,
		const RankFourTensor &	Eijkl
	)

protectedvirtual

Derived classes may employ this function to record stuff or do other computations prior to the return-mapping algorithm.

We know that (trial_stress_params, intnl_old) is inadmissible when this is called

Parameters

trial_stress_params[in]	The trial values of the stress parameters
stress_trial[in]	Trial stress tensor
intnl_old[in]	Old value of the internal parameters.
yf[in]	The yield functions at (p_trial, q_trial, intnl_old)
Eijkl[in]	The elasticity tensor

Reimplemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, TensileStressUpdate, and CappedMohrCoulombCosseratStressUpdate.

Definition at line 678 of file MultiParameterPlasticityStressUpdate.C.

684 {

685 }

Referenced by updateState().

◆ propagateQpStatefulProperties()

void MultiParameterPlasticityStressUpdate::propagateQpStatefulProperties ( )

overrideprotectedvirtual

If updateState is not called during a timestep, this will be.

This method allows derived classes to set internal parameters from their Old values, for instance

Reimplemented from StressUpdateBase.

Definition at line 152 of file MultiParameterPlasticityStressUpdate.C.

 {
   _plastic_strain[_qp] = _plastic_strain_old[_qp];
   std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
   _max_iter_used[_qp] = _max_iter_used_old[_qp];
 }

◆ requiresIsotropicTensor()

virtual bool StressUpdateBase::requiresIsotropicTensor ( )

pure virtualinherited

Does the model require the elasticity tensor to be isotropic?

Implemented in RadialReturnStressUpdate, CappedDruckerPragerStressUpdate, CappedDruckerPragerCosseratStressUpdate, LinearViscoelasticStressUpdate, CappedMohrCoulombCosseratStressUpdate, CappedWeakPlaneStressUpdate, CappedWeakInclinedPlaneStressUpdate, CappedWeakPlaneCosseratStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by ComputeMultipleInelasticStress::initialSetup().

◆ resetProperties()

void StressUpdateBase::resetProperties ( )

inlinefinalinherited

Definition at line 123 of file StressUpdateBase.h.

123 {}

◆ resetQpProperties()

void StressUpdateBase::resetQpProperties ( )

inlinefinalinherited

Retained as empty methods to avoid a warning from Material.C in framework. These methods are unused in all inheriting classes and should not be overwritten.

Definition at line 122 of file StressUpdateBase.h.

122 {}

◆ setEffectiveElasticity()

virtual void MultiParameterPlasticityStressUpdate::setEffectiveElasticity ( const RankFourTensor & Eijkl )

protectedpure virtual

Sets _Eij and _En and _Cij.

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, TensileStressUpdate, and CappedMohrCoulombCosseratStressUpdate.

Referenced by updateState().

◆ setInelasticStrainIncrementAfterReturn()

void MultiParameterPlasticityStressUpdate::setInelasticStrainIncrementAfterReturn	(	const RankTwoTensor &	stress_trial,
		Real	gaE,
		const yieldAndFlow &	smoothed_q,
		const RankFourTensor &	elasticity_tensor,
		const RankTwoTensor &	returned_stress,
		RankTwoTensor &	inelastic_strain_increment
	)		const

protectedvirtual

Sets inelastic strain increment from the returned configuration This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called.

Derived classes may override this function

Parameters

stress_trial[in]	The trial value of stress
gaE[in]	The value of gaE after the return-map process has completed successfully
smoothed_q[in]	Holds the current value of yield function and derivatives evaluated at the returned configuration
elasticity_tensor[in]	The elasticity tensor
returned_stress[in]	The stress after the return-map process
inelastic_strain_increment[out]	The inelastic strain increment resulting from this return-map

Reimplemented in TwoParameterPlasticityStressUpdate.

Definition at line 786 of file MultiParameterPlasticityStressUpdate.C.

 {
   const std::vector<RankTwoTensor> dsp_dstress = dstress_param_dstress(returned_stress);
   inelastic_strain_increment = RankTwoTensor();
   for (unsigned i = 0; i < _num_sp; ++i)
     inelastic_strain_increment += smoothed_q.dg[i] * dsp_dstress[i];
   inelastic_strain_increment *= gaE / _En;
 }

Referenced by updateState().

◆ setIntnlDerivativesV()

virtual void MultiParameterPlasticityStressUpdate::setIntnlDerivativesV	(	const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	current_stress_params,
		const std::vector< Real > &	intnl,
		std::vector< std::vector< Real >> &	dintnl
	)		const

protectedpure virtual

Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.

Derived classes must override this.

Parameters

trial_stress_params[in]	The trial stress parameters
current_stress_params[in]	The current stress parameters
intnl[in]	The current value of the internal parameters
dintnl[out]	The derivatives dintnl[i][j] = d(intnl[i])/d(stress_param j)

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by dVardTrial(), and nrStep().

◆ setIntnlValuesV()

virtual void MultiParameterPlasticityStressUpdate::setIntnlValuesV	(	const std::vector< Real > &	trial_stress_params,
		const std::vector< Real > &	current_stress_params,
		const std::vector< Real > &	intnl_old,
		std::vector< Real > &	intnl
	)		const

protectedpure virtual

Sets the internal parameters based on the trial values of stress_params, their current values, and the old values of the internal parameters.

Derived classes must override this.

Parameters

trial_stress_params[in]	The trial stress parameters (eg trial_p and trial_q)
current_stress_params[in]	The current stress parameters (eg p and q)
intnl_old[out]	Old value of internal parameters
intnl[out]	The value of internal parameters to be set

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, and TensileStressUpdate.

Referenced by lineSearch(), and updateState().

◆ setQp()

void StressUpdateBase::setQp ( unsigned int qp )

inherited

Sets the value of the global variable _qp for inheriting classes.

Definition at line 43 of file StressUpdateBase.C.

 {
   _qp = qp;
 }

◆ setStressAfterReturnV()

virtual void MultiParameterPlasticityStressUpdate::setStressAfterReturnV	(	const RankTwoTensor &	stress_trial,
		const std::vector< Real > &	stress_params,
		Real	gaE,
		const std::vector< Real > &	intnl,
		const yieldAndFlow &	smoothed_q,
		const RankFourTensor &	Eijkl,
		RankTwoTensor &	stress
	)		const

protectedpure virtual

Sets stress from the admissible parameters.

This is called after the return-map process has completed successfully in stress_param space, just after finalizeReturnProcess has been called. Derived classes must override this function

Parameters

stress_trial[in]	The trial value of stress
stress_params[in]	The value of the stress_params after the return-map process has completed successfully
gaE[in]	The value of gaE after the return-map process has completed successfully
intnl[in]	The value of the internal parameters after the return-map process has completed successfully
smoothed_q[in]	Holds the current value of yield function and derivatives evaluated at the returned state
Eijkl[in]	The elasticity tensor
stress[out]	The returned value of the stress tensor

Implemented in TwoParameterPlasticityStressUpdate, CappedMohrCoulombStressUpdate, TensileStressUpdate, and CappedMohrCoulombCosseratStressUpdate.

Referenced by updateState().

◆ smoothAllQuantities()

MultiParameterPlasticityStressUpdate::yieldAndFlow MultiParameterPlasticityStressUpdate::smoothAllQuantities	(	const std::vector< Real > &	stress_params,
		const std::vector< Real > &	intnl
	)		const

protected

Calculates all yield functions and derivatives, and then performs the smoothing scheme.

Parameters

stress_params[in]	The stress parameters (eg stress_params[0] = stress_zz and stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
intnl[in]	Internal parameters

Returns: The smoothed yield function and derivatives

Definition at line 400 of file MultiParameterPlasticityStressUpdate.C.

 {
   std::vector<yieldAndFlow> all_q(_num_yf, yieldAndFlow(_num_sp, _num_intnl));
   computeAllQV(stress_params, intnl, all_q);
  
   /* This routine holds the key to my smoothing strategy.  It
    * may be proved that this smoothing strategy produces a
    * yield surface that is both C2 differentiable and convex,
    * assuming the individual yield functions are C2 and
    * convex too.
    * Of course all the derivatives must also be smoothed.
    * Also, I assume that d(flow potential)/dstress gets smoothed
    * by the Yield Function (which produces a C2 flow potential).
    * See the line identified in the loop below.
    * Only time will tell whether this is a good strategy, but it
    * works well in all tests so far.  Convexity is irrelevant
    * for the non-associated case, but at least the return-map
    * problem should always have a unique solution.
    * For two yield functions+flows, labelled 1 and 2, we
    * should have
    * d(g1 - g2) . d(f1 - f2) >= 0
    * If not then the return-map problem for even the
    * multi-surface plasticity with no smoothing won't have a
    * unique solution.  If the multi-surface plasticity has
    * a unique solution then the smoothed version defined
    * below will too.
    */
  
   // res_f is the index that contains the smoothed yieldAndFlow
   std::size_t res_f = 0;
  
   for (std::size_t a = 1; a < all_q.size(); ++a)
   {
     if (all_q[res_f].f >= all_q[a].f + _smoothing_tol)
       // no smoothing is needed: res_f is already indexes the largest yield function
       continue;
     else if (all_q[a].f >= all_q[res_f].f + _smoothing_tol)
     {
       // no smoothing is needed, and res_f needs to index to all_q[a]
       res_f = a;
       continue;
     }
     else
     {
       // smoothing is required
       const Real f_diff = all_q[res_f].f - all_q[a].f;
       const Real ism = ismoother(f_diff);
       const Real sm = smoother(f_diff);
       const Real dsm = dsmoother(f_diff);
       // we want: all_q[res_f].f = 0.5 * all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism,
       // but we have to do the derivatives first
       for (unsigned i = 0; i < _num_sp; ++i)
       {
         for (unsigned j = 0; j < _num_sp; ++j)
           all_q[res_f].d2g[i][j] =
               0.5 * (all_q[res_f].d2g[i][j] + all_q[a].d2g[i][j]) +
               dsm * (all_q[res_f].df[j] - all_q[a].df[j]) * (all_q[res_f].dg[i] - all_q[a].dg[i]) +
               sm * (all_q[res_f].d2g[i][j] - all_q[a].d2g[i][j]);
         for (unsigned j = 0; j < _num_intnl; ++j)
           all_q[res_f].d2g_di[i][j] = 0.5 * (all_q[res_f].d2g_di[i][j] + all_q[a].d2g_di[i][j]) +
                                       dsm * (all_q[res_f].df_di[j] - all_q[a].df_di[j]) *
                                           (all_q[res_f].dg[i] - all_q[a].dg[i]) +
                                       sm * (all_q[res_f].d2g_di[i][j] - all_q[a].d2g_di[i][j]);
       }
       for (unsigned i = 0; i < _num_sp; ++i)
       {
         all_q[res_f].df[i] = 0.5 * (all_q[res_f].df[i] + all_q[a].df[i]) +
                              sm * (all_q[res_f].df[i] - all_q[a].df[i]);
         // whether the following (smoothing g with f's smoother) is a good strategy remains to be
         // seen...
         all_q[res_f].dg[i] = 0.5 * (all_q[res_f].dg[i] + all_q[a].dg[i]) +
                              sm * (all_q[res_f].dg[i] - all_q[a].dg[i]);
       }
       for (unsigned i = 0; i < _num_intnl; ++i)
         all_q[res_f].df_di[i] = 0.5 * (all_q[res_f].df_di[i] + all_q[a].df_di[i]) +
                                 sm * (all_q[res_f].df_di[i] - all_q[a].df_di[i]);
       all_q[res_f].f = 0.5 * (all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism;
     }
   }
   return all_q[res_f];
 }

Referenced by lineSearch(), and updateState().

◆ smoother()

Real MultiParameterPlasticityStressUpdate::smoother ( Real f_diff ) const

protected

Derivative of ismoother.

Definition at line 510 of file MultiParameterPlasticityStressUpdate.C.

 {
   if (std::abs(f_diff) >= _smoothing_tol)
     return 0.0;
   switch (_smoother_function_type)
   {
     case SmootherFunctionType::cos:
       return 0.5 * std::sin(f_diff * M_PI * 0.5 / _smoothing_tol);
     case SmootherFunctionType::poly1:
       return 0.75 / _smoothing_tol *
              (f_diff - (_smoothing_tol / 3.0) * Utility::pow<3>(f_diff / _smoothing_tol));
     case SmootherFunctionType::poly2:
       return 0.625 / _smoothing_tol *
              (f_diff - (_smoothing_tol / 5.0) * Utility::pow<5>(f_diff / _smoothing_tol));
     case SmootherFunctionType::poly3:
       return (7.0 / 12.0 / _smoothing_tol) *
              (f_diff - (_smoothing_tol / 7.0) * Utility::pow<7>(f_diff / _smoothing_tol));
     default:
       return 0.0;
   }
 }

Referenced by smoothAllQuantities().

◆ updateState()

void MultiParameterPlasticityStressUpdate::updateState	(	RankTwoTensor &	strain_increment,
		RankTwoTensor &	inelastic_strain_increment,
		const RankTwoTensor &	rotation_increment,
		RankTwoTensor &	stress_new,
		const RankTwoTensor &	stress_old,
		const RankFourTensor &	elasticity_tensor,
		const RankTwoTensor &	elastic_strain_old,
		bool	compute_full_tangent_operator,
		RankFourTensor &	tangent_operator
	)

overrideprotectedvirtual

Given a strain increment that results in a trial stress, perform some procedure (such as an iterative return-mapping process) to produce an admissible stress, an elastic strain increment and an inelastic strain increment.

If _fe_problem.currentlyComputingJacobian() = true, then updateState also computes d(stress)/d(strain) (or some approximation to it).

This method is called by ComputeMultipleInelasticStress. This method is pure virutal: all inheriting classes must overwrite this method.

Parameters

strain_increment	Upon input: the strain increment. Upon output: the elastic strain increment
inelastic_strain_increment	The inelastic_strain resulting from the interative procedure
rotation_increment	The finite-strain rotation increment
stress_new	Upon input: the trial stress that results from applying strain_increment as an elastic strain. Upon output: the admissible stress
stress_old	The old value of stress
elasticity_tensor	The elasticity tensor
compute_full_tangent_operator	The calling routine would like the full consistent tangent operator to be placed in tangent_operator, if possible. This is irrelevant if _fe_problem.currentlyComputingJacobian() = false
tangent_operator	d(stress)/d(strain), or some approximation to it If compute_full_tangent_operator=false, then tangent_operator=elasticity_tensor is an appropriate choice. tangent_operator is only computed if _fe_problem.currentlyComputingJacobian() = true

Implements StressUpdateBase.

Definition at line 160 of file MultiParameterPlasticityStressUpdate.C.

 {
   // Size _yf[_qp] appropriately
   _yf[_qp].assign(_num_yf, 0);
   // _plastic_strain and _intnl are usually sized appropriately because they are stateful, but this
   // Material may be used from a DiracKernel where stateful materials are not allowed.  The best we
   // can do is:
   if (_intnl[_qp].size() != _num_intnl)
     initQpStatefulProperties();
  
   initializeReturnProcess();
  
   if (_t_step >= 2)
     _step_one = false;
  
   // initially assume an elastic deformation
   std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
  
   _iter[_qp] = 0.0;
   _max_iter_used[_qp] = std::max(_max_iter_used[_qp], _max_iter_used_old[_qp]);
   _linesearch_needed[_qp] = 0.0;
  
   computeStressParams(stress_new, _trial_sp);
   yieldFunctionValuesV(_trial_sp, _intnl[_qp], _yf[_qp]);
  
   if (yieldF(_yf[_qp]) <= _f_tol)
   {
     _plastic_strain[_qp] = _plastic_strain_old[_qp];
     inelastic_strain_increment.zero();
     if (_fe_problem.currentlyComputingJacobian())
       tangent_operator = elasticity_tensor;
     return;
   }
  
   _stress_trial = stress_new;
   /* The trial stress must be inadmissible
    * so we need to return to the yield surface.  The following
    * equations must be satisfied.
    *
    * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
    * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
    * ...
    * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
    * 0 = rhs[N] = f(S, intnl)
    *
    * as well as equations defining intnl parameters as functions of
    * stress_params, trial_stress_params and intnl_old
    *
    * The unknowns are S[0], ..., S[N-1], gaE, and the intnl parameters.
    * Here gaE = ga * _En (the _En serves to make gaE similar magnitude to S)
    * I find it convenient to solve the first N+1 equations for p, q and gaE,
    * while substituting the "intnl parameters" equations into these during the solve process
    */
  
   for (auto & deriv : _dvar_dtrial)
     deriv.assign(_num_sp, 0.0);
  
   preReturnMapV(_trial_sp, stress_new, _intnl_old[_qp], _yf[_qp], elasticity_tensor);
  
   setEffectiveElasticity(elasticity_tensor);
  
   if (_step_one)
     std::copy(_definitely_ok_sp.begin(), _definitely_ok_sp.end(), _ok_sp.begin());
   else
     computeStressParams(stress_old, _ok_sp);
   std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _ok_intnl.begin());
  
   // Return-map problem: must apply the following changes in stress_params,
   // and find the returned stress_params and gaE
   for (unsigned i = 0; i < _num_sp; ++i)
     _del_stress_params[i] = _trial_sp[i] - _ok_sp[i];
  
   Real step_taken = 0.0; // amount of del_stress_params that we've applied and the return-map
                          // problem has succeeded
   Real step_size = 1.0;  // potentially can apply del_stress_params in substeps
   Real gaE_total = 0.0;
  
   // current values of the yield function, derivatives, etc
   yieldAndFlow smoothed_q;
  
   // In the following sub-stepping procedure it is possible that
   // the last step is an elastic step, and therefore smoothed_q won't
   // be computed on the last step, so we have to compute it.
   bool smoothed_q_calculated = false;
  
   while (step_taken < 1.0 && step_size >= _min_step_size)
   {
     if (1.0 - step_taken < step_size)
       // prevent over-shoots of substepping
       step_size = 1.0 - step_taken;
  
     // trial variables in terms of admissible variables
     for (unsigned i = 0; i < _num_sp; ++i)
       _trial_sp[i] = _ok_sp[i] + step_size * _del_stress_params[i];
  
     // initialize variables that are to be found via Newton-Raphson
     _current_sp = _trial_sp;
     Real gaE = 0.0;
  
     // flags indicating failure of Newton-Raphson and line-search
     int nr_failure = 0;
     int ls_failure = 0;
  
     // NR iterations taken in this substep
     unsigned step_iter = 0;
  
     // The residual-squared for the line-search
     Real res2 = 0.0;
  
     if (step_size < 1.0 && yieldF(_trial_sp, _ok_intnl) <= _f_tol)
       // This is an elastic step
       // The "step_size < 1.0" in above condition is for efficiency: we definitely
       // know that this is a plastic step if step_size = 1.0
       smoothed_q_calculated = false;
     else
     {
       // this is a plastic step
  
       // initialize current_sp, gaE and current_intnl based on the non-smoothed situation
       initializeVarsV(_trial_sp, _ok_intnl, _current_sp, gaE, _current_intnl);
       // and find the smoothed yield function, flow potential and derivatives
       smoothed_q = smoothAllQuantities(_current_sp, _current_intnl);
       smoothed_q_calculated = true;
       calculateRHS(_trial_sp, _current_sp, gaE, smoothed_q, _rhs);
       res2 = calculateRes2(_rhs);
  
       // Perform a Newton-Raphson with linesearch to get current_sp, gaE, and also smoothed_q
       while (res2 > _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0)
       {
         // solve the linear system and store the answer (the "updates") in rhs
         nr_failure = nrStep(smoothed_q, _trial_sp, _current_sp, _current_intnl, gaE, _rhs);
         if (nr_failure != 0)
           break;
  
         // handle precision loss
         if (precisionLoss(_rhs, _current_sp, gaE))
         {
           if (_warn_about_precision_loss)
           {
             Moose::err << "MultiParameterPlasticityStressUpdate: precision-loss in element "
                        << _current_elem->id() << " quadpoint=" << _qp << ".  Stress_params =";
             for (unsigned i = 0; i < _num_sp; ++i)
               Moose::err << " " << _current_sp[i];
             Moose::err << " gaE = " << gaE << "\n";
           }
           res2 = 0.0;
           break;
         }
  
         // apply (parts of) the updates, re-calculate smoothed_q, and res2
         ls_failure = lineSearch(res2,
                                 _current_sp,
                                 gaE,
                                 _trial_sp,
                                 smoothed_q,
                                 _ok_intnl,
                                 _current_intnl,
                                 _rhs,
                                 _linesearch_needed[_qp]);
         step_iter++;
       }
     }
     if (res2 <= _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0 &&
         gaE >= 0.0)
     {
       // this Newton-Raphson worked fine, or this was an elastic step
       std::copy(_current_sp.begin(), _current_sp.end(), _ok_sp.begin());
       gaE_total += gaE;
       step_taken += step_size;
       setIntnlValuesV(_trial_sp, _ok_sp, _ok_intnl, _intnl[_qp]);
       std::copy(_intnl[_qp].begin(), _intnl[_qp].end(), _ok_intnl.begin());
       // calculate dp/dp_trial, dp/dq_trial, etc, for Jacobian
       dVardTrial(!smoothed_q_calculated,
                  _trial_sp,
                  _ok_sp,
                  gaE,
                  _ok_intnl,
                  smoothed_q,
                  step_size,
                  compute_full_tangent_operator,
                  _dvar_dtrial);
       if (static_cast<Real>(step_iter) > _iter[_qp])
         _iter[_qp] = static_cast<Real>(step_iter);
       if (static_cast<Real>(step_iter) > _max_iter_used[_qp])
         _max_iter_used[_qp] = static_cast<Real>(step_iter);
       step_size *= 1.1;
     }
     else
     {
       // Newton-Raphson + line-search process failed
       step_size *= 0.5;
     }
   }
  
   if (step_size < _min_step_size)
     errorHandler("MultiParameterPlasticityStressUpdate: Minimum step-size violated");
  
   // success!
   finalizeReturnProcess(rotation_increment);
   yieldFunctionValuesV(_ok_sp, _intnl[_qp], _yf[_qp]);
  
   if (!smoothed_q_calculated)
     smoothed_q = smoothAllQuantities(_ok_sp, _intnl[_qp]);
  
   setStressAfterReturnV(
       _stress_trial, _ok_sp, gaE_total, _intnl[_qp], smoothed_q, elasticity_tensor, stress_new);
  
   setInelasticStrainIncrementAfterReturn(_stress_trial,
                                          gaE_total,
                                          smoothed_q,
                                          elasticity_tensor,
                                          stress_new,
                                          inelastic_strain_increment);
  
   strain_increment = strain_increment - inelastic_strain_increment;
   _plastic_strain[_qp] = _plastic_strain_old[_qp] + inelastic_strain_increment;
  
   if (_fe_problem.currentlyComputingJacobian())
     // for efficiency, do not compute the tangent operator if not currently computing Jacobian
     consistentTangentOperatorV(_stress_trial,
                                _trial_sp,
                                stress_new,
                                _ok_sp,
                                gaE_total,
                                smoothed_q,
                                elasticity_tensor,
                                compute_full_tangent_operator,
                                _dvar_dtrial,
                                tangent_operator);
 }

◆ validParams()

InputParameters MultiParameterPlasticityStressUpdate::validParams ( )

static

Definition at line 23 of file MultiParameterPlasticityStressUpdate.C.

 {
   InputParameters params = StressUpdateBase::validParams();
   params.addClassDescription("Return-map and Jacobian algorithms for plastic models where the "
                              "yield function and flow directions depend on multiple functions of "
                              "stress");
   params.addRangeCheckedParam<unsigned int>(
       "max_NR_iterations",
       20,
       "max_NR_iterations>0",
       "Maximum number of Newton-Raphson iterations allowed during the return-map algorithm");
   params.addParam<bool>("perform_finite_strain_rotations",
                         false,
                         "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.addRequiredParam<Real>(
       "smoothing_tol",
       "Intersections of the yield surfaces will be smoothed by this amount (this "
       "is measured in units of stress).  Often this is related to other physical "
       "parameters (eg, 0.1*cohesion) but it is important to set this small enough "
       "so that the individual yield surfaces do not mix together in the smoothing "
       "process to produce a result where no stress is admissible (for example, "
       "mixing together tensile and compressive failure envelopes).");
   params.addRequiredParam<Real>("yield_function_tol",
                                 "The return-map process will be deemed to have converged if all "
                                 "yield functions are within yield_function_tol of zero.  If this "
                                 "is set very low then precision-loss might be encountered: if the "
                                 "code detects precision loss then it also deems the return-map "
                                 "process has converged.");
   MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
   params.addParam<Real>("min_step_size",
                         1.0,
                         "In order to help the Newton-Raphson procedure, the applied strain "
                         "increment may be applied in sub-increments of size greater than this "
                         "value.  Usually it is better for Moose's nonlinear convergence to "
                         "increase max_NR_iterations rather than decrease this parameter.");
   params.addParam<bool>("warn_about_precision_loss",
                         false,
                         "Output a message to the console every "
                         "time precision-loss is encountered "
                         "during the Newton-Raphson process");
   params.addParam<std::vector<Real>>("admissible_stress",
                                      "A single admissible value of the value of the stress "
                                      "parameters for internal parameters = 0.  This is used "
                                      "to initialize the return-mapping algorithm during the first "
                                      "nonlinear iteration.  If not given then it is assumed that "
                                      "stress parameters = 0 is admissible.");
   MooseEnum smoother_fcn_enum("cos poly1 poly2 poly3", "cos");
   params.addParam<MooseEnum>("smoother_function_type",
                              smoother_fcn_enum,
                              "Type of smoother function to use.  'cos' means (-a/pi)cos(pi x/2/a), "
                              "'polyN' means a polynomial of degree 2N+2");
   params.addParamNamesToGroup("smoother_function_type", "Advanced");
   return params;
 }

Referenced by CappedMohrCoulombStressUpdate::validParams(), TensileStressUpdate::validParams(), and TwoParameterPlasticityStressUpdate::validParams().

◆ yieldF() [1/2]

Real MultiParameterPlasticityStressUpdate::yieldF	(	const std::vector< Real > &	stress_params,
		const std::vector< Real > &	intnl
	)		const

protected

Computes the smoothed yield function.

Parameters

stress_params	The stress parameters (eg stress_params[0] = stress_zz and stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
intnl	The internal parameters (eg intnl[0] is shear, intnl[1] is tensile)

Returns: The smoothed yield function value

Definition at line 688 of file MultiParameterPlasticityStressUpdate.C.

 {
   std::vector<Real> yfs(_num_yf);
   yieldFunctionValuesV(stress_params, intnl, yfs);
   return yieldF(yfs);
 }

Referenced by updateState().

◆ yieldF() [2/2]

Real MultiParameterPlasticityStressUpdate::yieldF ( const std::vector< Real > & yfs ) const

protected

Computes the smoothed yield function.

Parameters

yfs	The values of the individual yield functions

Returns: The smoothed yield function value

Definition at line 697 of file MultiParameterPlasticityStressUpdate.C.

 {
   Real yf = yfs[0];
   for (std::size_t i = 1; i < yfs.size(); ++i)
     if (yf >= yfs[i] + _smoothing_tol)
       // no smoothing is needed, and yf is the biggest yield function
       continue;
     else if (yfs[i] >= yf + _smoothing_tol)
       // no smoothing is needed, and yfs[i] is the biggest yield function
       yf = yfs[i];
     else
       yf = 0.5 * (yf + yfs[i] + _smoothing_tol) + ismoother(yf - yfs[i]);
   return yf;
 }

◆ yieldFunctionValuesV()

virtual void MultiParameterPlasticityStressUpdate::yieldFunctionValuesV	(	const std::vector< Real > &	stress_params,
		const std::vector< Real > &	intnl,
		std::vector< Real > &	yf
	)		const

protectedpure virtual

Computes the values of the yield functions, given stress_params and intnl parameters.

Derived classes must override this, to provide the values of the yield functions in yf.

Parameters

	stress_params[in]	The stress parameters
	intnl[in]	The internal parameters
[out]	yf	The yield function values

Implemented in CappedMohrCoulombStressUpdate, TwoParameterPlasticityStressUpdate, and TensileStressUpdate.

Referenced by updateState(), and yieldF().

Member Data Documentation

◆ _base_name

const std::string StressUpdateBase::_base_name

protectedinherited

Name used as a prefix for all material properties related to the stress update model.

Definition at line 128 of file StressUpdateBase.h.

◆ _Cij

std::vector<std::vector<Real> > MultiParameterPlasticityStressUpdate::_Cij

protected

_Cij[i, j] * _Eij[j, k] = 1 iff j == k

Definition at line 144 of file MultiParameterPlasticityStressUpdate.h.

Referenced by consistentTangentOperatorV(), CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(), TensileStressUpdate::setEffectiveElasticity(), CappedMohrCoulombStressUpdate::setEffectiveElasticity(), and TwoParameterPlasticityStressUpdate::setEffectiveElasticity().

◆ _current_intnl

std::vector<Real> MultiParameterPlasticityStressUpdate::_current_intnl

private

The current values of the internal params during the Newton-Raphson.

Definition at line 737 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _current_sp

std::vector<Real> MultiParameterPlasticityStressUpdate::_current_sp

private

The current values of the stress params during the Newton-Raphson.

Definition at line 732 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _definitely_ok_sp

const std::vector<Real> MultiParameterPlasticityStressUpdate::_definitely_ok_sp

protected

An admissible value of stress_params at the initial time.

Definition at line 135 of file MultiParameterPlasticityStressUpdate.h.

Referenced by MultiParameterPlasticityStressUpdate(), and updateState().

◆ _del_stress_params

std::vector<Real> MultiParameterPlasticityStressUpdate::_del_stress_params

private

_del_stress_params = trial_stress_params - ok_sp This is fixed at the beginning of the return-map process, irrespective of substepping.

The return-map problem is: apply del_stress_params to stress_prams, and then find an admissible (returned) stress_params and gaE

Definition at line 727 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _dvar_dtrial

std::vector<std::vector<Real> > MultiParameterPlasticityStressUpdate::_dvar_dtrial

private

d({stress_param[i], gaE})/d(trial_stress_param[j])

Definition at line 703 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _Eij

std::vector<std::vector<Real> > MultiParameterPlasticityStressUpdate::_Eij

protected

E[i, j] in the system of equations to be solved.

Definition at line 138 of file MultiParameterPlasticityStressUpdate.h.

Referenced by calculateRHS(), dnRHSdVar(), dVardTrial(), CappedMohrCoulombStressUpdate::initializeVarsV(), CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(), TensileStressUpdate::setEffectiveElasticity(), CappedMohrCoulombStressUpdate::setEffectiveElasticity(), TwoParameterPlasticityStressUpdate::setEffectiveElasticity(), TensileStressUpdate::setIntnlDerivativesV(), CappedMohrCoulombStressUpdate::setIntnlDerivativesV(), TensileStressUpdate::setIntnlValuesV(), and CappedMohrCoulombStressUpdate::setIntnlValuesV().

◆ _En

Real MultiParameterPlasticityStressUpdate::_En

protected

normalising factor

Definition at line 141 of file MultiParameterPlasticityStressUpdate.h.

Referenced by calculateRHS(), consistentTangentOperatorV(), dnRHSdVar(), dVardTrial(), CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(), TensileStressUpdate::setEffectiveElasticity(), CappedMohrCoulombStressUpdate::setEffectiveElasticity(), TwoParameterPlasticityStressUpdate::setEffectiveElasticity(), and setInelasticStrainIncrementAfterReturn().

◆ _f_tol

const Real MultiParameterPlasticityStressUpdate::_f_tol

protected

The yield-function tolerance.

Definition at line 165 of file MultiParameterPlasticityStressUpdate.h.

Referenced by CappedMohrCoulombStressUpdate::initializeVarsV(), and updateState().

◆ _f_tol2

const Real MultiParameterPlasticityStressUpdate::_f_tol2

protected

Square of the yield-function tolerance.

Definition at line 168 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _intnl

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

protected

internal parameters

Definition at line 190 of file MultiParameterPlasticityStressUpdate.h.

Referenced by CappedWeakPlaneCosseratStressUpdate::consistentTangentOperator(), CappedWeakPlaneStressUpdate::consistentTangentOperator(), initQpStatefulProperties(), propagateQpStatefulProperties(), and updateState().

◆ _intnl_old

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

protected

old values of internal parameters

Definition at line 193 of file MultiParameterPlasticityStressUpdate.h.

Referenced by propagateQpStatefulProperties(), and updateState().

◆ _iter

MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_iter

protected

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

Definition at line 199 of file MultiParameterPlasticityStressUpdate.h.

Referenced by initQpStatefulProperties(), and updateState().

◆ _linesearch_needed

MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_linesearch_needed

protected

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

Definition at line 208 of file MultiParameterPlasticityStressUpdate.h.

Referenced by initQpStatefulProperties(), and updateState().

◆ _max_iter_used

MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_max_iter_used

protected

Maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation.

Definition at line 202 of file MultiParameterPlasticityStressUpdate.h.

Referenced by initQpStatefulProperties(), propagateQpStatefulProperties(), and updateState().

◆ _max_iter_used_old

const MaterialProperty<Real>& MultiParameterPlasticityStressUpdate::_max_iter_used_old

protected

Old value of maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation.

Definition at line 205 of file MultiParameterPlasticityStressUpdate.h.

Referenced by propagateQpStatefulProperties(), and updateState().

◆ _max_nr_its

const unsigned MultiParameterPlasticityStressUpdate::_max_nr_its

protected

Maximum number of Newton-Raphson iterations allowed in the return-map process.

Definition at line 153 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _min_step_size

const Real MultiParameterPlasticityStressUpdate::_min_step_size

protected

In order to help the Newton-Raphson procedure, the applied strain increment may be applied in sub-increments of size greater than this value.

Definition at line 175 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _num_intnl

const unsigned MultiParameterPlasticityStressUpdate::_num_intnl

protected

Number of internal parameters.

Definition at line 150 of file MultiParameterPlasticityStressUpdate.h.

Referenced by CappedMohrCoulombStressUpdate::computeAllQV(), dnRHSdVar(), dVardTrial(), initQpStatefulProperties(), nrStep(), smoothAllQuantities(), and updateState().

◆ _num_sp

const unsigned MultiParameterPlasticityStressUpdate::_num_sp

protected

Number of stress parameters.

Definition at line 132 of file MultiParameterPlasticityStressUpdate.h.

Referenced by calculateRHS(), TensileStressUpdate::computeAllQV(), CappedMohrCoulombStressUpdate::computeAllQV(), TensileStressUpdate::consistentTangentOperatorV(), CappedMohrCoulombStressUpdate::consistentTangentOperatorV(), consistentTangentOperatorV(), dnRHSdVar(), dVardTrial(), TensileStressUpdate::initializeVarsV(), CappedMohrCoulombStressUpdate::initializeVarsV(), lineSearch(), MultiParameterPlasticityStressUpdate(), nrStep(), precisionLoss(), CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(), TensileStressUpdate::setEffectiveElasticity(), CappedMohrCoulombStressUpdate::setEffectiveElasticity(), setInelasticStrainIncrementAfterReturn(), CappedMohrCoulombStressUpdate::setIntnlDerivativesV(), smoothAllQuantities(), and updateState().

◆ _num_yf

const unsigned MultiParameterPlasticityStressUpdate::_num_yf

protected

Number of yield functions.

Definition at line 147 of file MultiParameterPlasticityStressUpdate.h.

Referenced by TensileStressUpdate::computeAllQV(), CappedMohrCoulombStressUpdate::computeAllQV(), initQpStatefulProperties(), smoothAllQuantities(), updateState(), and yieldF().

◆ _ok_intnl

std::vector<Real> MultiParameterPlasticityStressUpdate::_ok_intnl

private

The state (ok_sp, ok_intnl) is known to be admissible.

Definition at line 718 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _ok_sp

std::vector<Real> MultiParameterPlasticityStressUpdate::_ok_sp

private

The state (ok_sp, ok_intnl) is known to be admissible, so ok_sp are stress_params that are "OK".

If the strain_increment is applied in substeps then ok_sp is updated after each sub strain_increment is applied and the return-map is successful. At the end of the entire return-map process _ok_sp will contain the stress_params where (_ok_sp, _intnl) is admissible.

Definition at line 713 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _perform_finite_strain_rotations

const bool MultiParameterPlasticityStressUpdate::_perform_finite_strain_rotations

protected

Whether to perform finite-strain rotations.

Definition at line 156 of file MultiParameterPlasticityStressUpdate.h.

Referenced by CappedWeakInclinedPlaneStressUpdate::finalizeReturnProcess(), and CappedWeakInclinedPlaneStressUpdate::initializeReturnProcess().

◆ _plastic_strain

MaterialProperty<RankTwoTensor>& MultiParameterPlasticityStressUpdate::_plastic_strain

protected

plastic strain

Definition at line 184 of file MultiParameterPlasticityStressUpdate.h.

Referenced by initQpStatefulProperties(), propagateQpStatefulProperties(), and updateState().

◆ _plastic_strain_old

const MaterialProperty<RankTwoTensor>& MultiParameterPlasticityStressUpdate::_plastic_strain_old

protected

Old value of plastic strain.

Definition at line 187 of file MultiParameterPlasticityStressUpdate.h.

Referenced by propagateQpStatefulProperties(), and updateState().

◆ _rhs

std::vector<Real> MultiParameterPlasticityStressUpdate::_rhs

private

0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i] 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i] ...

0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i] 0 = rhs[N] = f(S, intnl) Here N = num_sp

Definition at line 698 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _smoother_function_type

enum MultiParameterPlasticityStressUpdate::SmootherFunctionType MultiParameterPlasticityStressUpdate::_smoother_function_type

private

Referenced by dsmoother(), ismoother(), and smoother().

◆ _smoothing_tol

const Real MultiParameterPlasticityStressUpdate::_smoothing_tol

protected

Smoothing tolerance: edges of the yield surface get smoothed by this amount.

Definition at line 159 of file MultiParameterPlasticityStressUpdate.h.

Referenced by CappedDruckerPragerStressUpdate::CappedDruckerPragerStressUpdate(), CappedWeakPlaneStressUpdate::CappedWeakPlaneStressUpdate(), dsmoother(), ismoother(), smoothAllQuantities(), smoother(), and yieldF().

◆ _smoothing_tol2

const Real MultiParameterPlasticityStressUpdate::_smoothing_tol2

protected

Square of the smoothing tolerance.

Definition at line 162 of file MultiParameterPlasticityStressUpdate.h.

Referenced by ismoother().

◆ _step_one

bool MultiParameterPlasticityStressUpdate::_step_one

protected

handles case of initial_stress that is inadmissible being supplied

Definition at line 178 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _stress_trial

RankTwoTensor MultiParameterPlasticityStressUpdate::_stress_trial

private

"Trial" value of stress that is set at the beginning of the return-map process.

It is fixed at stress_old + Eijkl * strain_increment irrespective of any sub-stepping

Definition at line 688 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _tensor_dimensionality

constexpr static unsigned MultiParameterPlasticityStressUpdate::_tensor_dimensionality = 3

staticconstexprprotected

Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics)

Definition at line 129 of file MultiParameterPlasticityStressUpdate.h.

◆ _trial_sp

std::vector<Real> MultiParameterPlasticityStressUpdate::_trial_sp

private

"Trial" value of stress_params that initializes the return-map process This is derived from stress = stress_old + Eijkl * strain_increment.

However, since the return-map process can fail and be restarted by applying strain_increment in multiple substeps, _trial_sp can vary from substep to substep.

Definition at line 681 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _warn_about_precision_loss

const bool MultiParameterPlasticityStressUpdate::_warn_about_precision_loss

protected

Output a warning message if precision loss is encountered during the return-map process.

Definition at line 181 of file MultiParameterPlasticityStressUpdate.h.

Referenced by updateState().

◆ _yf

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

protected

yield functions

Definition at line 196 of file MultiParameterPlasticityStressUpdate.h.

Referenced by initQpStatefulProperties(), and updateState().

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

tensor_mechanics/include/materials/MultiParameterPlasticityStressUpdate.h
tensor_mechanics/src/materials/MultiParameterPlasticityStressUpdate.C

Classes

Public Member Functions

Static Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Types

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ SmootherFunctionType

Constructor & Destructor Documentation

◆ MultiParameterPlasticityStressUpdate()

Member Function Documentation

◆ calculateRes2()

◆ calculateRHS()

◆ computeAllQV()

◆ computeStressParams()

◆ computeTimeStepLimit()

◆ consistentTangentOperatorV()

◆ d2stress_param_dstress()

◆ dnRHSdVar()

◆ dsmoother()

◆ dstress_param_dstress()

◆ dVardTrial()

◆ errorHandler()

◆ finalizeReturnProcess()

◆ getTangentCalculationMethod()

◆ initializeReturnProcess()

◆ initializeVarsV()

◆ initQpStatefulProperties()

◆ isIsotropic()

◆ ismoother()

◆ lineSearch()

◆ nrStep()

◆ precisionLoss()

◆ preReturnMapV()

◆ propagateQpStatefulProperties()

◆ requiresIsotropicTensor()

◆ resetProperties()

◆ resetQpProperties()

◆ setEffectiveElasticity()

◆ setInelasticStrainIncrementAfterReturn()

◆ setIntnlDerivativesV()

◆ setIntnlValuesV()

◆ setQp()

◆ setStressAfterReturnV()

◆ smoothAllQuantities()

◆ smoother()

◆ updateState()

◆ validParams()

◆ yieldF() [1/2]

◆ yieldF() [2/2]

◆ yieldFunctionValuesV()

Member Data Documentation

◆ _base_name

◆ _Cij

◆ _current_intnl

◆ _current_sp

◆ _definitely_ok_sp

◆ _del_stress_params

◆ _dvar_dtrial

◆ _Eij

◆ _En

◆ _f_tol

◆ _f_tol2

◆ _intnl

◆ _intnl_old

◆ _iter

◆ _linesearch_needed

◆ _max_iter_used

◆ _max_iter_used_old

◆ _max_nr_its

◆ _min_step_size

◆ _num_intnl

◆ _num_sp

◆ _num_yf

◆ _ok_intnl

◆ _ok_sp

◆ _perform_finite_strain_rotations