FiniteStrainMohrCoulombMulti implements rate-independent non-associative mohr-coulomb with hardening/softening in the finite-strain framework, using planar (non-smoothed) surfaces. More...

#include <SolidMechanicsPlasticMohrCoulombMulti.h>

Inheritance diagram for SolidMechanicsPlasticMohrCoulombMulti:

[legend]

Public Types
typedef DataFileName	DataFileParameterType

Public Member Functions
	SolidMechanicsPlasticMohrCoulombMulti (const InputParameters &parameters)

virtual unsigned int	numberSurfaces () const override
	The number of yield surfaces for this plasticity model. More...

virtual void	yieldFunctionV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &f) const override
	Calculates the yield functions. More...

virtual void	dyieldFunction_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &df_dstress) const override
	The derivative of yield functions with respect to stress. More...

virtual void	dyieldFunction_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &df_dintnl) const override
	The derivative of yield functions with respect to the internal parameter. More...

virtual void	flowPotentialV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &r) const override
	The flow potentials. More...

virtual void	dflowPotential_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankFourTensor > &dr_dstress) const override
	The derivative of the flow potential with respect to stress. More...

virtual void	dflowPotential_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dr_dintnl) const override
	The derivative of the flow potential with respect to the internal parameter. More...

virtual void	activeConstraints (const std::vector< Real > &f, const RankTwoTensor &stress, Real intnl, const RankFourTensor &Eijkl, std::vector< bool > &act, RankTwoTensor &returned_stress) const override
	The active yield surfaces, given a vector of yield functions. More...

virtual std::string	modelName () const override

virtual bool	useCustomReturnMap () const override
	Returns false. You will want to override this in your derived class if you write a custom returnMap function. More...

virtual bool	returnMap (const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const override
	Performs a custom return-map. More...

void	initialize ()

void	execute ()

void	finalize ()

virtual void	hardPotentialV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &h) const
	The hardening potential. More...

virtual void	dhardPotential_dstressV (const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dh_dstress) const
	The derivative of the hardening potential with respect to stress. More...

virtual void	dhardPotential_dintnlV (const RankTwoTensor &stress, Real intnl, std::vector< Real > &dh_dintnl) const
	The derivative of the hardening potential with respect to the internal parameter. More...

virtual bool	useCustomCTO () const
	Returns false. You will want to override this in your derived class if you write a custom consistent tangent operator function. More...

virtual RankFourTensor	consistentTangentOperator (const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const
	Calculates a custom consistent tangent operator. More...

bool	KuhnTuckerSingleSurface (Real yf, Real dpm, Real dpm_tol) const
	Returns true if the Kuhn-Tucker conditions for the single surface are satisfied. More...

SubProblem &	getSubProblem () const

bool	shouldDuplicateInitialExecution () const

virtual Real	spatialValue (const Point &) const

virtual const std::vector< Point >	spatialPoints () const

void	gatherSum (T &value)

void	gatherMax (T &value)

void	gatherMin (T &value)

void	gatherProxyValueMax (T1 &proxy, T2 &value)

void	gatherProxyValueMin (T1 &proxy, T2 &value)

void	setPrimaryThreadCopy (UserObject *primary)

UserObject *	primaryThreadCopy ()

std::set< UserObjectName >	getDependObjects () const

virtual bool	needThreadedCopy () const

const std::set< std::string > &	getRequestedItems () override

const std::set< std::string > &	getSuppliedItems () override

unsigned int	systemNumber () const

virtual bool	enabled () const

std::shared_ptr< MooseObject >	getSharedPtr ()

std::shared_ptr< const MooseObject >	getSharedPtr () const

MooseApp &	getMooseApp () const

const std::string &	type () const

virtual const std::string &	name () const

std::string	typeAndName () const

std::string	errorPrefix (const std::string &error_type) const

void	callMooseError (std::string msg, const bool with_prefix) const

MooseObjectParameterName	uniqueParameterName (const std::string &parameter_name) const

const InputParameters &	parameters () const

MooseObjectName	uniqueName () const

const T &	getParam (const std::string &name) const

std::vector< std::pair< T1, T2 > >	getParam (const std::string &param1, const std::string &param2) const

const T &	getRenamedParam (const std::string &old_name, const std::string &new_name) const

T	getCheckedPointerParam (const std::string &name, const std::string &error_string="") const

bool	isParamValid (const std::string &name) const

bool	isParamSetByUser (const std::string &nm) const

void	paramError (const std::string &param, Args... args) const

void	paramWarning (const std::string &param, Args... args) const

void	paramInfo (const std::string &param, Args... args) const

void	connectControllableParams (const std::string &parameter, const std::string &object_type, const std::string &object_name, const std::string &object_parameter) const

void	mooseError (Args &&... args) const

void	mooseErrorNonPrefixed (Args &&... args) const

void	mooseDocumentedError (const std::string &repo_name, const unsigned int issue_num, Args &&... args) const

void	mooseWarning (Args &&... args) const

void	mooseWarningNonPrefixed (Args &&... args) const

void	mooseDeprecated (Args &&... args) const

void	mooseInfo (Args &&... args) const

std::string	getDataFileName (const std::string &param) const

std::string	getDataFileNameByName (const std::string &relative_path) const

std::string	getDataFilePath (const std::string &relative_path) const

virtual void	initialSetup ()

virtual void	timestepSetup ()

virtual void	jacobianSetup ()

virtual void	residualSetup ()

virtual void	customSetup (const ExecFlagType &)

const ExecFlagEnum &	getExecuteOnEnum () const

UserObjectName	getUserObjectName (const std::string &param_name) const

const T &	getUserObject (const std::string &param_name, bool is_dependency=true) const

const T &	getUserObjectByName (const UserObjectName &object_name, bool is_dependency=true) const

const UserObject &	getUserObjectBase (const std::string &param_name, bool is_dependency=true) const

const UserObject &	getUserObjectBaseByName (const UserObjectName &object_name, bool is_dependency=true) const

const std::vector< MooseVariableScalar *> &	getCoupledMooseScalarVars ()

const std::set< TagID > &	getScalarVariableCoupleableVectorTags () const

const std::set< TagID > &	getScalarVariableCoupleableMatrixTags () const

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, MaterialData &material_data, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, MaterialData &material_data, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, const unsigned int state=0)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name, MaterialData &material_data)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data, const unsigned int state)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name)

std::pair< const MaterialProperty< T > *, std::set< SubdomainID > >	getBlockMaterialProperty (const MaterialPropertyName &name)

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialProperty (const std::string &name)

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialProperty ()

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialPropertyByName (const std::string &prop_name)

const MaterialProperty< T > &	getZeroMaterialProperty (Ts... args)

std::set< SubdomainID >	getMaterialPropertyBlocks (const std::string &name)

std::vector< SubdomainName >	getMaterialPropertyBlockNames (const std::string &name)

std::set< BoundaryID >	getMaterialPropertyBoundaryIDs (const std::string &name)

std::vector< BoundaryName >	getMaterialPropertyBoundaryNames (const std::string &name)

void	checkBlockAndBoundaryCompatibility (std::shared_ptr< MaterialBase > discrete)

std::unordered_map< SubdomainID, std::vector< MaterialBase *> >	buildRequiredMaterials (bool allow_stateful=true)

void	statefulPropertiesAllowed (bool)

bool	getMaterialPropertyCalled () const

const std::unordered_set< unsigned int > &	getMatPropDependencies () const

virtual void	resolveOptionalProperties ()

const GenericMaterialProperty< T, is_ad > &	getPossiblyConstantGenericMaterialPropertyByName (const MaterialPropertyName &prop_name, MaterialData &material_data, const unsigned int state)

bool	isImplicit ()

Moose::StateArg	determineState () const

virtual void	threadJoin (const UserObject &) override

virtual void	threadJoin (const UserObject &) override

virtual void	subdomainSetup () override

virtual void	subdomainSetup () override

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

const GenericOptionalMaterialProperty< T, is_ad > &	getGenericOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const GenericOptionalMaterialProperty< T, is_ad > &	getGenericOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalMaterialProperty< T > &	getOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalMaterialProperty< T > &	getOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalADMaterialProperty< T > &	getOptionalADMaterialProperty (const std::string &name)

const OptionalADMaterialProperty< T > &	getOptionalADMaterialProperty (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOld (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOld (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOlder (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOlder (const std::string &name)

MaterialBase &	getMaterial (const std::string &name)

MaterialBase &	getMaterial (const std::string &name)

MaterialBase &	getMaterialByName (const std::string &name, bool no_warn=false)

MaterialBase &	getMaterialByName (const std::string &name, bool no_warn=false)

bool	hasMaterialProperty (const std::string &name)

bool	hasMaterialProperty (const std::string &name)

bool	hasMaterialPropertyByName (const std::string &name)

bool	hasMaterialPropertyByName (const std::string &name)

bool	hasADMaterialProperty (const std::string &name)

bool	hasADMaterialProperty (const std::string &name)

bool	hasADMaterialPropertyByName (const std::string &name)

bool	hasADMaterialPropertyByName (const std::string &name)

bool	hasGenericMaterialProperty (const std::string &name)

bool	hasGenericMaterialProperty (const std::string &name)

bool	hasGenericMaterialPropertyByName (const std::string &name)

bool	hasGenericMaterialPropertyByName (const std::string &name)

const Function &	getFunction (const std::string &name) const

const Function &	getFunctionByName (const FunctionName &name) const

bool	hasFunction (const std::string &param_name) const

bool	hasFunctionByName (const FunctionName &name) const

bool	isDefaultPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

bool	hasPostprocessor (const std::string &param_name, const unsigned int index=0) const

bool	hasPostprocessorByName (const PostprocessorName &name) const

std::size_t	coupledPostprocessors (const std::string &param_name) const

const PostprocessorName &	getPostprocessorName (const std::string &param_name, const unsigned int index=0) const

const VectorPostprocessorValue &	getVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name, bool needs_broadcast) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name) const

bool	hasVectorPostprocessor (const std::string &param_name, const std::string &vector_name) const

bool	hasVectorPostprocessor (const std::string &param_name) const

bool	hasVectorPostprocessorByName (const VectorPostprocessorName &name, const std::string &vector_name) const

bool	hasVectorPostprocessorByName (const VectorPostprocessorName &name) const

const VectorPostprocessorName &	getVectorPostprocessorName (const std::string &param_name) const

T &	getSampler (const std::string &name)

Sampler &	getSampler (const std::string &name)

T &	getSamplerByName (const SamplerName &name)

Sampler &	getSamplerByName (const SamplerName &name)

virtual void	meshChanged ()

PerfGraph &	perfGraph ()

const PostprocessorValue &	getPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOld (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOld (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOlder (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOlder (const std::string &param_name, const unsigned int index=0) const

virtual const PostprocessorValue &	getPostprocessorValueByName (const PostprocessorName &name) const

virtual const PostprocessorValue &	getPostprocessorValueByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOldByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOldByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOlderByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOlderByName (const PostprocessorName &name) const

bool	isVectorPostprocessorDistributed (const std::string &param_name) const

bool	isVectorPostprocessorDistributed (const std::string &param_name) const

bool	isVectorPostprocessorDistributedByName (const VectorPostprocessorName &name) const

bool	isVectorPostprocessorDistributedByName (const VectorPostprocessorName &name) const

const Distribution &	getDistribution (const std::string &name) const

const T &	getDistribution (const std::string &name) const

const Distribution &	getDistribution (const std::string &name) const

const T &	getDistribution (const std::string &name) const

const Distribution &	getDistributionByName (const DistributionName &name) const

const T &	getDistributionByName (const std::string &name) const

const Distribution &	getDistributionByName (const DistributionName &name) const

const T &	getDistributionByName (const std::string &name) const

const Parallel::Communicator &	comm () const

processor_id_type	n_processors () const

processor_id_type	processor_id () const

Static Public Member Functions
static InputParameters	validParams ()

static void	sort (typename std::vector< T > &vector)

static void	sortDFS (typename std::vector< T > &vector)

static void	cyclicDependencyError (CyclicDependencyException< T2 > &e, const std::string &header)

Public Attributes
const Real	_f_tol
	Tolerance on yield function. More...

const Real	_ic_tol
	Tolerance on internal constraint. More...

const ConsoleStream	_console

Static Public Attributes
static constexpr PropertyValue::id_type	default_property_id

static constexpr PropertyValue::id_type	zero_property_id

static constexpr auto	SYSTEM

static constexpr auto	NAME

Protected Member Functions
virtual Real	cohesion (const Real internal_param) const
	cohesion as a function of residual value, rate, and internal_param More...

virtual Real	dcohesion (const Real internal_param) const
	d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param More...

virtual Real	phi (const Real internal_param) const
	phi as a function of residual value, rate, and internal_param More...

virtual Real	dphi (const Real internal_param) const
	d(phi)/d(internal_param) as a function of residual value, rate, and internal_param More...

virtual Real	psi (const Real internal_param) const
	psi as a function of residual value, rate, and internal_param More...

virtual Real	dpsi (const Real internal_param) const
	d(psi)/d(internal_param) as a function of residual value, rate, and internal_param More...

virtual Real	yieldFunction (const RankTwoTensor &stress, Real intnl) const
	The following functions are what you should override when building single-plasticity models. More...

virtual RankTwoTensor	dyieldFunction_dstress (const RankTwoTensor &stress, Real intnl) const
	The derivative of yield function with respect to stress. More...

virtual Real	dyieldFunction_dintnl (const RankTwoTensor &stress, Real intnl) const
	The derivative of yield function with respect to the internal parameter. More...

virtual RankTwoTensor	flowPotential (const RankTwoTensor &stress, Real intnl) const
	The flow potential. More...

virtual RankFourTensor	dflowPotential_dstress (const RankTwoTensor &stress, Real intnl) const
	The derivative of the flow potential with respect to stress. More...

virtual RankTwoTensor	dflowPotential_dintnl (const RankTwoTensor &stress, Real intnl) const
	The derivative of the flow potential with respect to the internal parameter. More...

virtual Real	hardPotential (const RankTwoTensor &stress, Real intnl) const
	The hardening potential. More...

virtual RankTwoTensor	dhardPotential_dstress (const RankTwoTensor &stress, Real intnl) const
	The derivative of the hardening potential with respect to stress. More...

virtual Real	dhardPotential_dintnl (const RankTwoTensor &stress, Real intnl) const
	The derivative of the hardening potential with respect to the internal parameter. More...

virtual void	addPostprocessorDependencyHelper (const PostprocessorName &name) const override

virtual void	addVectorPostprocessorDependencyHelper (const VectorPostprocessorName &name) const override

virtual void	addUserObjectDependencyHelper (const UserObject &uo) const override

void	addReporterDependencyHelper (const ReporterName &reporter_name) override

const ReporterName &	getReporterName (const std::string &param_name) const

T &	declareRestartableData (const std::string &data_name, Args &&... args)

ManagedValue< T >	declareManagedRestartableDataWithContext (const std::string &data_name, void *context, Args &&... args)

const T &	getRestartableData (const std::string &data_name) const

T &	declareRestartableDataWithContext (const std::string &data_name, void *context, Args &&... args)

T &	declareRecoverableData (const std::string &data_name, Args &&... args)

T &	declareRestartableDataWithObjectName (const std::string &data_name, const std::string &object_name, Args &&... args)

T &	declareRestartableDataWithObjectNameWithContext (const std::string &data_name, const std::string &object_name, void *context, Args &&... args)

std::string	restartableName (const std::string &data_name) const

const T &	getMeshProperty (const std::string &data_name, const std::string &prefix)

const T &	getMeshProperty (const std::string &data_name)

bool	hasMeshProperty (const std::string &data_name, const std::string &prefix) const

bool	hasMeshProperty (const std::string &data_name, const std::string &prefix) const

bool	hasMeshProperty (const std::string &data_name) const

bool	hasMeshProperty (const std::string &data_name) const

std::string	meshPropertyName (const std::string &data_name) const

PerfID	registerTimedSection (const std::string &section_name, const unsigned int level) const

PerfID	registerTimedSection (const std::string &section_name, const unsigned int level, const std::string &live_message, const bool print_dots=true) const

std::string	timedSectionName (const std::string &section_name) const

bool	isCoupledScalar (const std::string &var_name, unsigned int i=0) const

unsigned int	coupledScalarComponents (const std::string &var_name) const

unsigned int	coupledScalar (const std::string &var_name, unsigned int comp=0) const

libMesh::Order	coupledScalarOrder (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarValue (const std::string &var_name, unsigned int comp=0) const

const ADVariableValue &	adCoupledScalarValue (const std::string &var_name, unsigned int comp=0) const

const GenericVariableValue< is_ad > &	coupledGenericScalarValue (const std::string &var_name, unsigned int comp=0) const

const GenericVariableValue< false > &	coupledGenericScalarValue (const std::string &var_name, const unsigned int comp) const

const GenericVariableValue< true > &	coupledGenericScalarValue (const std::string &var_name, const unsigned int comp) const

const VariableValue &	coupledVectorTagScalarValue (const std::string &var_name, TagID tag, unsigned int comp=0) const

const VariableValue &	coupledMatrixTagScalarValue (const std::string &var_name, TagID tag, unsigned int comp=0) const

const VariableValue &	coupledScalarValueOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarValueOlder (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDot (const std::string &var_name, unsigned int comp=0) const

const ADVariableValue &	adCoupledScalarDot (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDot (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDotOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDu (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDotDu (const std::string &var_name, unsigned int comp=0) const

const MooseVariableScalar *	getScalarVar (const std::string &var_name, unsigned int comp) const

virtual void	checkMaterialProperty (const std::string &name, const unsigned int state)

void	markMatPropRequested (const std::string &)

MaterialPropertyName	getMaterialPropertyName (const std::string &name) const

void	checkExecutionStage ()

const T &	getReporterValue (const std::string &param_name, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, ReporterMode mode, const std::size_t time_index=0)

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

const GenericMaterialProperty< T, is_ad > *	defaultGenericMaterialProperty (const std::string &name)

const GenericMaterialProperty< T, is_ad > *	defaultGenericMaterialProperty (const std::string &name)

const MaterialProperty< T > *	defaultMaterialProperty (const std::string &name)

const MaterialProperty< T > *	defaultMaterialProperty (const std::string &name)

const ADMaterialProperty< T > *	defaultADMaterialProperty (const std::string &name)

const ADMaterialProperty< T > *	defaultADMaterialProperty (const std::string &name)

Static Protected Member Functions
static std::string	meshPropertyName (const std::string &data_name, const std::string &prefix)

Protected Attributes
SubProblem &	_subproblem

FEProblemBase &	_fe_problem

SystemBase &	_sys

const THREAD_ID	_tid

Assembly &	_assembly

const Moose::CoordinateSystemType &	_coord_sys

const bool	_duplicate_initial_execution

std::set< std::string >	_depend_uo

const bool &	_enabled

MooseApp &	_app

const std::string	_type

const std::string	_name

const InputParameters &	_pars

Factory &	_factory

ActionFactory &	_action_factory

const ExecFlagEnum &	_execute_enum

const ExecFlagType &	_current_execute_flag

MooseApp &	_restartable_app

const std::string	_restartable_system_name

const THREAD_ID	_restartable_tid

const bool	_restartable_read_only

FEProblemBase &	_mci_feproblem

MooseApp &	_pg_moose_app

const std::string	_prefix

FEProblemBase &	_sc_fe_problem

const THREAD_ID	_sc_tid

const Real &	_real_zero

const VariableValue &	_scalar_zero

const Point &	_point_zero

const InputParameters &	_mi_params

const std::string	_mi_name

const MooseObjectName	_mi_moose_object_name

FEProblemBase &	_mi_feproblem

SubProblem &	_mi_subproblem

const THREAD_ID	_mi_tid

const Moose::MaterialDataType	_material_data_type

MaterialData &	_material_data

bool	_stateful_allowed

bool	_get_material_property_called

std::vector< std::unique_ptr< PropertyValue > >	_default_properties

std::unordered_set< unsigned int >	_material_property_dependencies

const MaterialPropertyName	_get_suffix

const bool	_use_interpolated_state

const InputParameters &	_ti_params

FEProblemBase &	_ti_feproblem

bool	_is_implicit

Real &	_t

const Real &	_t_old

int &	_t_step

Real &	_dt

Real &	_dt_old

bool	_is_transient

const Parallel::Communicator &	_communicator

Static Protected Attributes
static const std::string	_interpolated_old

static const std::string	_interpolated_older

Private Types
enum	return_type { tip110100 = 0, tip010101 = 1, edge010100 = 2, edge000101 = 3, plane000100 = 4 }

Private Member Functions
void	yieldFunctionEigvals (Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
	Calculates the yield functions given the eigenvalues of stress. More...

void	df_dsig (const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
	this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi), or sin_angle = sin(psi), respectively More...

void	perturbStress (const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
	perturbs the stress tensor in the case of almost-equal eigenvalues. More...

bool	KuhnTuckerOK (const std::vector< Real > &yf, const std::vector< Real > &dpm, Real ep_plastic_tolerance) const
	Returns true if the Kuhn-Tucker conditions are satisfied. More...

bool	doReturnMap (const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
	See doco for returnMap function. More...

bool	returnTip (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
	Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are returned. More...

bool	returnPlane (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
	Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...

bool	returnEdge000101 (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
	Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...

bool	returnEdge010100 (const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
	Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk. More...

Private Attributes
const SolidMechanicsHardeningModel &	_cohesion
	Hardening model for cohesion. More...

const SolidMechanicsHardeningModel &	_phi
	Hardening model for phi. More...

const SolidMechanicsHardeningModel &	_psi
	Hardening model for psi. More...

const unsigned int	_max_iters
	Maximum Newton-Raphison iterations in the custom returnMap algorithm. More...

const Real	_shift
	yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalues More...

const bool	_use_custom_returnMap
	Whether to use the custom return-map algorithm. More...

Detailed Description

FiniteStrainMohrCoulombMulti implements rate-independent non-associative mohr-coulomb with hardening/softening in the finite-strain framework, using planar (non-smoothed) surfaces.

Definition at line 19 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Member Enumeration Documentation

◆ return_type

enum SolidMechanicsPlasticMohrCoulombMulti::return_type

private

Enumerator
tip110100
tip010101
edge010100
edge000101
plane000100

Definition at line 315 of file SolidMechanicsPlasticMohrCoulombMulti.h.

   {
     tip110100 = 0,
     tip010101 = 1,
     edge010100 = 2,
     edge000101 = 3,
     plane000100 = 4
   };

Constructor & Destructor Documentation

◆ SolidMechanicsPlasticMohrCoulombMulti()

SolidMechanicsPlasticMohrCoulombMulti::SolidMechanicsPlasticMohrCoulombMulti ( const InputParameters & parameters )

Definition at line 61 of file SolidMechanicsPlasticMohrCoulombMulti.C.

   : SolidMechanicsPlasticModel(parameters),
     _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
     _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
     _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
     _max_iters(getParam<unsigned int>("max_iterations")),
     _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
     _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
 {
   if (_shift < 0)
     mooseError("Value of 'shift' in SolidMechanicsPlasticMohrCoulombMulti must not be negative\n");
   if (_shift > _f_tol)
     _console << "WARNING: value of 'shift' in SolidMechanicsPlasticMohrCoulombMulti is probably "
                 "set too high"
              << std::endl;
   if (LIBMESH_DIM != 3)
     mooseError("SolidMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
   MooseRandom::seed(0);
 }

Member Function Documentation

◆ activeConstraints()

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

overridevirtual

The active yield surfaces, given a vector of yield functions.

This is used by FiniteStrainMultiPlasticity to determine the initial set of active constraints at the trial (stress, intnl) configuration. It is up to you (the coder) to determine how accurate you want the returned_stress to be. Currently it is only used by FiniteStrainMultiPlasticity to estimate a good starting value for the Newton-Rahson procedure, so currently it may not need to be super perfect.

Parameters

	f	values of the yield functions
	stress	stress tensor
	intnl	internal parameter
	Eijkl	elasticity tensor (stress = Eijkl*strain)
[out]	act	act[i] = true if the i_th yield function is active
[out]	returned_stress	Approximate value of the returned stress

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 248 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   act.assign(6, false);
 
   if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
       f[5] <= _f_tol)
   {
     returned_stress = stress;
     return;
   }
 
   Real returned_intnl;
   std::vector<Real> dpm(6);
   RankTwoTensor delta_dp;
   std::vector<Real> yf(6);
   bool trial_stress_inadmissible;
   doReturnMap(stress,
               intnl,
               Eijkl,
               0.0,
               returned_stress,
               returned_intnl,
               dpm,
               delta_dp,
               yf,
               trial_stress_inadmissible);
 
   for (unsigned i = 0; i < 6; ++i)
     act[i] = (dpm[i] > 0);
 }

◆ cohesion()

Real SolidMechanicsPlasticMohrCoulombMulti::cohesion ( const Real internal_param ) const

protectedvirtual

cohesion as a function of residual value, rate, and internal_param

Definition at line 285 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), returnTip(), and yieldFunctionV().

 {
   return _cohesion.value(internal_param);
 }

◆ consistentTangentOperator()

RankFourTensor SolidMechanicsPlasticModel::consistentTangentOperator	(	const RankTwoTensor &	trial_stress,
		Real	intnl_old,
		const RankTwoTensor &	stress,
		Real	intnl,
		const RankFourTensor &	E_ijkl,
		const std::vector< Real > &	cumulative_pm
	)		const

virtualinherited

Calculates a custom consistent tangent operator.

You may choose to over-ride this in your derived SolidMechanicsPlasticXXXX class.

(Note, if you over-ride returnMap, you will probably want to override consistentTangentOpertor too, otherwise it will default to E_ijkl.)

Parameters

stress_old	trial stress before returning
intnl_old	internal parameter before returning
stress	current returned stress state
intnl	internal parameter
E_ijkl	elasticity tensor
cumulative_pm	the cumulative plastic multipliers

Returns: the consistent tangent operator: E_ijkl if not over-ridden

Reimplemented in SolidMechanicsPlasticTensileMulti, SolidMechanicsPlasticDruckerPragerHyperbolic, SolidMechanicsPlasticMeanCapTC, and SolidMechanicsPlasticJ2.

Definition at line 252 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticJ2::consistentTangentOperator(), SolidMechanicsPlasticDruckerPragerHyperbolic::consistentTangentOperator(), SolidMechanicsPlasticMeanCapTC::consistentTangentOperator(), and SolidMechanicsPlasticTensileMulti::consistentTangentOperator().

 {
   return E_ijkl;
 }

◆ dcohesion()

Real SolidMechanicsPlasticMohrCoulombMulti::dcohesion ( const Real internal_param ) const

protectedvirtual

d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 291 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

 {
   return _cohesion.derivative(internal_param);
 }

◆ df_dsig()

void SolidMechanicsPlasticMohrCoulombMulti::df_dsig	(	const RankTwoTensor &	stress,
		Real	sin_angle,
		std::vector< RankTwoTensor > &	df
	)		const

private

this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi), or sin_angle = sin(psi), respectively

Definition at line 144 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dstressV(), and flowPotentialV().

 {
   std::vector<Real> eigvals;
   std::vector<RankTwoTensor> deigvals;
   stress.dsymmetricEigenvalues(eigvals, deigvals);
 
   if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
     perturbStress(stress, eigvals, deigvals);
 
   df.resize(6);
   df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
   df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
   df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
   df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
   df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
   df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
 }

◆ dflowPotential_dintnl()

RankTwoTensor SolidMechanicsPlasticModel::dflowPotential_dintnl	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of the flow potential with respect to the internal parameter.

Parameters

stress	the stress at which to calculate the flow potential
intnl	internal parameter

Returns: dr_dintnl(i, j) = dr(i, j)/dintnl

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticJ2, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 131 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dflowPotential_dintnlV().

 {
   return RankTwoTensor();
 }

◆ dflowPotential_dintnlV()

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

overridevirtual

The derivative of the flow potential with respect to the internal parameter.

Parameters

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

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 228 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   const Real cos_angle = std::cos(psi(intnl));
   const Real dsin_angle = cos_angle * dpsi(intnl);
 
   std::vector<Real> eigvals;
   std::vector<RankTwoTensor> deigvals;
   stress.dsymmetricEigenvalues(eigvals, deigvals);
 
   if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
     perturbStress(stress, eigvals, deigvals);
 
   dr_dintnl.resize(6);
   dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
   dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
   dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
 }

◆ dflowPotential_dstress()

RankFourTensor SolidMechanicsPlasticModel::dflowPotential_dstress	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of the flow potential with respect to stress.

Parameters

stress	the stress at which to calculate the flow potential
intnl	internal parameter

Returns: dr_dstress(i, j, k, l) = dr(i, j)/dstress(k, l)

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticIsotropicSD, SolidMechanicsPlasticJ2, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticOrthotropic, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, SolidMechanicsPlasticDruckerPragerHyperbolic, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 117 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dflowPotential_dstressV().

 {
   return RankFourTensor();
 }

◆ dflowPotential_dstressV()

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

overridevirtual

The derivative of the flow potential with respect to stress.

Parameters

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

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 204 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   std::vector<RankFourTensor> d2eigvals;
   stress.d2symmetricEigenvalues(d2eigvals);
 
   const Real sinpsi = std::sin(psi(intnl));
 
   dr_dstress.resize(6);
   dr_dstress[0] =
       0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
   dr_dstress[1] =
       0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
   dr_dstress[2] =
       0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
   dr_dstress[3] =
       0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
   dr_dstress[4] =
       0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
   dr_dstress[5] =
       0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
 }

◆ dhardPotential_dintnl()

Real SolidMechanicsPlasticModel::dhardPotential_dintnl	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of the hardening potential with respect to the internal parameter.

Parameters

stress	the stress at which to calculate the hardening potentials
intnl	internal parameter

Returns: the derivative

Reimplemented in SolidMechanicsPlasticMeanCapTC.

Definition at line 172 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dhardPotential_dintnlV().

 {
   return 0.0;
 }

◆ dhardPotential_dintnlV()

void SolidMechanicsPlasticModel::dhardPotential_dintnlV	(	const RankTwoTensor &	stress,
		Real	intnl,
		std::vector< Real > &	dh_dintnl
	)		const

virtualinherited

The derivative of the hardening potential with respect to the internal parameter.

Parameters

	stress	the stress at which to calculate the hardening potentials
	intnl	internal parameter
[out]	dh_dintnl	dh_dintnl[alpha] = dh[alpha]/dintnl

Definition at line 178 of file SolidMechanicsPlasticModel.C.

 {
   dh_dintnl.resize(numberSurfaces(), dhardPotential_dintnl(stress, intnl));
 }

◆ dhardPotential_dstress()

RankTwoTensor SolidMechanicsPlasticModel::dhardPotential_dstress	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of the hardening potential with respect to stress.

Parameters

stress	the stress at which to calculate the hardening potentials
intnl	internal parameter

Returns: dh_dstress(i, j) = dh/dstress(i, j)

Reimplemented in SolidMechanicsPlasticMeanCapTC.

Definition at line 158 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dhardPotential_dstressV().

 {
   return RankTwoTensor();
 }

◆ dhardPotential_dstressV()

void SolidMechanicsPlasticModel::dhardPotential_dstressV	(	const RankTwoTensor &	stress,
		Real	intnl,
		std::vector< RankTwoTensor > &	dh_dstress
	)		const

virtualinherited

The derivative of the hardening potential with respect to stress.

Parameters

	stress	the stress at which to calculate the hardening potentials
	intnl	internal parameter
[out]	dh_dstress	dh_dstress[alpha](i, j) = dh[alpha]/dstress(i, j)

Definition at line 164 of file SolidMechanicsPlasticModel.C.

 {
   dh_dstress.assign(numberSurfaces(), dhardPotential_dstress(stress, intnl));
 }

◆ doReturnMap()

bool SolidMechanicsPlasticMohrCoulombMulti::doReturnMap	(	const RankTwoTensor &	trial_stress,
		Real	intnl_old,
		const RankFourTensor &	E_ijkl,
		Real	ep_plastic_tolerance,
		RankTwoTensor &	returned_stress,
		Real &	returned_intnl,
		std::vector< Real > &	dpm,
		RankTwoTensor &	delta_dp,
		std::vector< Real > &	yf,
		bool &	trial_stress_inadmissible
	)		const

private

See doco for returnMap function.

The interface is identical to this one. This one can be called internally regardless of the value of _use_custom_returnMap

Definition at line 363 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by activeConstraints(), and returnMap().

 {
   mooseAssert(dpm.size() == 6,
               "SolidMechanicsPlasticMohrCoulombMulti size of dpm should be 6 but it is "
                   << dpm.size());
 
   std::vector<Real> eigvals;
   RankTwoTensor eigvecs;
   trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
   eigvals[0] += _shift;
   eigvals[2] -= _shift;
 
   Real sinphi = std::sin(phi(intnl_old));
   Real cosphi = std::cos(phi(intnl_old));
   Real coh = cohesion(intnl_old);
   Real cohcos = coh * cosphi;
 
   yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
 
   if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol && yf[3] <= _f_tol && yf[4] <= _f_tol &&
       yf[5] <= _f_tol)
   {
     // purely elastic (trial_stress, intnl_old)
     trial_stress_inadmissible = false;
     return true;
   }
 
   trial_stress_inadmissible = true;
   delta_dp.zero();
   returned_stress = RankTwoTensor();
 
   // these are the normals to the 6 yield surfaces, which are const because of the assumption of no
   // psi hardening
   std::vector<RealVectorValue> norm(6);
   const Real sinpsi = std::sin(psi(intnl_old));
   const Real oneminus = 0.5 * (1 - sinpsi);
   const Real oneplus = 0.5 * (1 + sinpsi);
   norm[0](0) = oneplus;
   norm[0](1) = -oneminus;
   norm[0](2) = 0;
   norm[1](0) = -oneminus;
   norm[1](1) = oneplus;
   norm[1](2) = 0;
   norm[2](0) = oneplus;
   norm[2](1) = 0;
   norm[2](2) = -oneminus;
   norm[3](0) = -oneminus;
   norm[3](1) = 0;
   norm[3](2) = oneplus;
   norm[4](0) = 0;
   norm[4](1) = oneplus;
   norm[4](2) = -oneminus;
   norm[5](0) = 0;
   norm[5](1) = -oneminus;
   norm[5](2) = oneplus;
 
   // the flow directions are these norm multiplied by Eijkl.
   // I call the flow directions "n".
   // In the following I assume that the Eijkl is
   // for an isotropic situation.  Then I don't have to
   // rotate to the principal-stress frame, and i don't
   // have to worry about strange off-diagonal things
   std::vector<RealVectorValue> n(6);
   for (unsigned ys = 0; ys < 6; ++ys)
     for (unsigned i = 0; i < 3; ++i)
       for (unsigned j = 0; j < 3; ++j)
         n[ys](i) += E_ijkl(i, i, j, j) * norm[ys](j);
   const Real mag_E = E_ijkl(0, 0, 0, 0);
 
   // With non-zero Poisson's ratio and hardening
   // it is not computationally cheap to know whether
   // the trial stress will return to the tip, edge,
   // or plane.  The following at least
   // gives a not-completely-stupid guess
   // trial_order[0] = type of return to try first
   // trial_order[1] = type of return to try second
   // trial_order[2] = type of return to try third
   // trial_order[3] = type of return to try fourth
   // trial_order[4] = type of return to try fifth
   // In the following the "binary" stuff indicates the
   // deactive (0) and active (1) surfaces, eg
   // 110100 means that surfaces 0, 1 and 3 are active
   // and 2, 4 and 5 are deactive
   const unsigned int number_of_return_paths = 5;
   std::vector<int> trial_order(number_of_return_paths);
   if (yf[1] > _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
   {
     trial_order[0] = tip110100;
     trial_order[1] = edge010100;
     trial_order[2] = plane000100;
     trial_order[3] = edge000101;
     trial_order[4] = tip010101;
   }
   else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
   {
     trial_order[0] = edge000101;
     trial_order[1] = plane000100;
     trial_order[2] = tip110100;
     trial_order[3] = tip010101;
     trial_order[4] = edge010100;
   }
   else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] <= _f_tol)
   {
     trial_order[0] = plane000100;
     trial_order[1] = edge000101;
     trial_order[2] = edge010100;
     trial_order[3] = tip110100;
     trial_order[4] = tip010101;
   }
   else
   {
     trial_order[0] = edge010100;
     trial_order[1] = plane000100;
     trial_order[2] = edge000101;
     trial_order[3] = tip110100;
     trial_order[4] = tip010101;
   }
 
   unsigned trial;
   bool nr_converged = false;
   bool kt_success = false;
   std::vector<RealVectorValue> ntip(3);
   std::vector<Real> dpmtip(3);
 
   for (trial = 0; trial < number_of_return_paths; ++trial)
   {
     switch (trial_order[trial])
     {
       case tip110100:
         for (unsigned int i = 0; i < 3; ++i)
         {
           ntip[0](i) = n[0](i);
           ntip[1](i) = n[1](i);
           ntip[2](i) = n[3](i);
         }
         kt_success = returnTip(eigvals,
                                ntip,
                                dpmtip,
                                returned_stress,
                                intnl_old,
                                sinphi,
                                cohcos,
                                0,
                                nr_converged,
                                ep_plastic_tolerance,
                                yf);
         if (nr_converged && kt_success)
         {
           dpm[0] = dpmtip[0];
           dpm[1] = dpmtip[1];
           dpm[3] = dpmtip[2];
           dpm[2] = dpm[4] = dpm[5] = 0;
         }
         break;
 
       case tip010101:
         for (unsigned int i = 0; i < 3; ++i)
         {
           ntip[0](i) = n[1](i);
           ntip[1](i) = n[3](i);
           ntip[2](i) = n[5](i);
         }
         kt_success = returnTip(eigvals,
                                ntip,
                                dpmtip,
                                returned_stress,
                                intnl_old,
                                sinphi,
                                cohcos,
                                0,
                                nr_converged,
                                ep_plastic_tolerance,
                                yf);
         if (nr_converged && kt_success)
         {
           dpm[1] = dpmtip[0];
           dpm[3] = dpmtip[1];
           dpm[5] = dpmtip[2];
           dpm[0] = dpm[2] = dpm[4] = 0;
         }
         break;
 
       case edge000101:
         kt_success = returnEdge000101(eigvals,
                                       n,
                                       dpm,
                                       returned_stress,
                                       intnl_old,
                                       sinphi,
                                       cohcos,
                                       0,
                                       mag_E,
                                       nr_converged,
                                       ep_plastic_tolerance,
                                       yf);
         break;
 
       case edge010100:
         kt_success = returnEdge010100(eigvals,
                                       n,
                                       dpm,
                                       returned_stress,
                                       intnl_old,
                                       sinphi,
                                       cohcos,
                                       0,
                                       mag_E,
                                       nr_converged,
                                       ep_plastic_tolerance,
                                       yf);
         break;
 
       case plane000100:
         kt_success = returnPlane(eigvals,
                                  n,
                                  dpm,
                                  returned_stress,
                                  intnl_old,
                                  sinphi,
                                  cohcos,
                                  0,
                                  nr_converged,
                                  ep_plastic_tolerance,
                                  yf);
         break;
     }
 
     if (nr_converged && kt_success)
       break;
   }
 
   if (trial == number_of_return_paths)
   {
     sinphi = std::sin(phi(intnl_old));
     cosphi = std::cos(phi(intnl_old));
     coh = cohesion(intnl_old);
     cohcos = coh * cosphi;
     yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
     Moose::err << "Trial stress = \n";
     trial_stress.print(Moose::err);
     Moose::err << "which has eigenvalues = " << eigvals[0] << " " << eigvals[1] << " " << eigvals[2]
                << "\n";
     Moose::err << "and yield functions = " << yf[0] << " " << yf[1] << " " << yf[2] << " " << yf[3]
                << " " << yf[4] << " " << yf[5] << "\n";
     Moose::err << "Internal parameter = " << intnl_old << std::endl;
     mooseError("SolidMechanicsPlasticMohrCoulombMulti: FAILURE!  You probably need to implement a "
                "line search if your hardening is too severe, or you need to tune your tolerances "
                "(eg, yield_function_tolerance should be a little smaller than (young "
                "modulus)*ep_plastic_tolerance).\n");
     return false;
   }
 
   // success
 
   returned_intnl = intnl_old;
   for (unsigned i = 0; i < 6; ++i)
     returned_intnl += dpm[i];
   for (unsigned i = 0; i < 6; ++i)
     for (unsigned j = 0; j < 3; ++j)
       delta_dp(j, j) += dpm[i] * norm[i](j);
   returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
   delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
   return true;
 }

◆ dphi()

Real SolidMechanicsPlasticMohrCoulombMulti::dphi ( const Real internal_param ) const

protectedvirtual

d(phi)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 303 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by dyieldFunction_dintnlV(), returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

 {
   return _phi.derivative(internal_param);
 }

◆ dpsi()

Real SolidMechanicsPlasticMohrCoulombMulti::dpsi ( const Real internal_param ) const

protectedvirtual

d(psi)/d(internal_param) as a function of residual value, rate, and internal_param

Definition at line 315 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by dflowPotential_dintnlV().

 {
   return _psi.derivative(internal_param);
 }

◆ dyieldFunction_dintnl()

Real SolidMechanicsPlasticModel::dyieldFunction_dintnl	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of yield function with respect to the internal parameter.

Parameters

stress	the stress at which to calculate the yield function
intnl	internal parameter

Returns: the derivative

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticJ2, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 90 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dyieldFunction_dintnlV().

 {
   return 0.0;
 }

◆ dyieldFunction_dintnlV()

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

overridevirtual

The derivative of yield functions with respect to the internal parameter.

Parameters

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	df_dintnl	df_dintnl[alpha] = df[alpha]/dintnl

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 173 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   std::vector<Real> eigvals;
   stress.symmetricEigenvalues(eigvals);
   eigvals[0] += _shift;
   eigvals[2] -= _shift;
 
   const Real sin_angle = std::sin(phi(intnl));
   const Real cos_angle = std::cos(phi(intnl));
   const Real dsin_angle = cos_angle * dphi(intnl);
   const Real dcos_angle = -sin_angle * dphi(intnl);
   const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
 
   df_dintnl.resize(6);
   df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
   df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
   df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
 }

◆ dyieldFunction_dstress()

RankTwoTensor SolidMechanicsPlasticModel::dyieldFunction_dstress	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The derivative of yield function with respect to stress.

Parameters

stress	the stress at which to calculate the yield function
intnl	internal parameter

Returns: df_dstress(i, j) = dyieldFunction/dstress(i, j)

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticIsotropicSD, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticJ2, SolidMechanicsPlasticOrthotropic, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 75 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::dyieldFunction_dstressV().

 {
   return RankTwoTensor();
 }

◆ dyieldFunction_dstressV()

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

overridevirtual

The derivative of yield functions with respect to stress.

Parameters

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	df_dstress	df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 165 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   const Real sinphi = std::sin(phi(intnl));
   df_dsig(stress, sinphi, df_dstress);
 }

◆ execute()

void SolidMechanicsPlasticModel::execute ( )

virtualinherited

Implements GeneralUserObject.

Definition at line 45 of file SolidMechanicsPlasticModel.C.

46 {

47 }

◆ finalize()

void SolidMechanicsPlasticModel::finalize ( )

virtualinherited

Implements GeneralUserObject.

Definition at line 50 of file SolidMechanicsPlasticModel.C.

51 {

52 }

◆ flowPotential()

RankTwoTensor SolidMechanicsPlasticModel::flowPotential	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The flow potential.

Parameters

stress	the stress at which to calculate the flow potential
intnl	internal parameter

Returns: the flow potential

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticIsotropicSD, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticJ2, SolidMechanicsPlasticOrthotropic, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 104 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::flowPotentialV().

 {
   return RankTwoTensor();
 }

◆ flowPotentialV()

void SolidMechanicsPlasticMohrCoulombMulti::flowPotentialV	(	const RankTwoTensor &	stress,
		Real	intnl,
		std::vector< RankTwoTensor > &	r
	)		const

overridevirtual

The flow potentials.

Parameters

	stress	the stress at which to calculate the flow potential
	intnl	internal parameter
[out]	r	r[alpha] is the flow potential for the "alpha" yield function

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 195 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   const Real sinpsi = std::sin(psi(intnl));
   df_dsig(stress, sinpsi, r);
 }

◆ hardPotential()

Real SolidMechanicsPlasticModel::hardPotential	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The hardening potential.

Parameters

stress	the stress at which to calculate the hardening potential
intnl	internal parameter

Returns: the hardening potential

Reimplemented in SolidMechanicsPlasticMeanCapTC.

Definition at line 145 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::hardPotentialV().

 {
   return -1.0;
 }

◆ hardPotentialV()

void SolidMechanicsPlasticModel::hardPotentialV	(	const RankTwoTensor &	stress,
		Real	intnl,
		std::vector< Real > &	h
	)		const

virtualinherited

The hardening potential.

Parameters

	stress	the stress at which to calculate the hardening potential
	intnl	internal parameter
[out]	h	h[alpha] is the hardening potential for the "alpha" yield function

Definition at line 150 of file SolidMechanicsPlasticModel.C.

 {
   h.assign(numberSurfaces(), hardPotential(stress, intnl));
 }

◆ initialize()

void SolidMechanicsPlasticModel::initialize ( )

virtualinherited

Implements GeneralUserObject.

Definition at line 40 of file SolidMechanicsPlasticModel.C.

41 {

42 }

◆ KuhnTuckerOK()

bool SolidMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK	(	const std::vector< Real > &	yf,
		const std::vector< Real > &	dpm,
		Real	ep_plastic_tolerance
	)		const

private

Returns true if the Kuhn-Tucker conditions are satisfied.

Parameters

yf	The six yield function values
dpm	The six plastic multipliers
ep_plastic_tolerance	The tolerance on the plastic strain (if dpm>-ep_plastic_tolerance then it is grouped as "non-negative" in the Kuhn-Tucker conditions).

Definition at line 1145 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   mooseAssert(
       yf.size() == 6,
       "SolidMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
           << yf.size());
   mooseAssert(
       dpm.size() == 6,
       "SolidMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
           << dpm.size());
   for (unsigned i = 0; i < 6; ++i)
     if (!SolidMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
       return false;
   return true;
 }

◆ KuhnTuckerSingleSurface()

bool SolidMechanicsPlasticModel::KuhnTuckerSingleSurface	(	Real	yf,
		Real	dpm,
		Real	dpm_tol
	)		const

inherited

Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.

Parameters

yf	Yield function value
dpm	plastic multiplier
dpm_tol	tolerance on plastic multiplier: viz dpm>-dpm_tol means "dpm is non-negative"

Definition at line 246 of file SolidMechanicsPlasticModel.C.

Referenced by KuhnTuckerOK(), SolidMechanicsPlasticTensileMulti::KuhnTuckerOK(), and SolidMechanicsPlasticModel::returnMap().

 {
   return (dpm == 0 && yf <= _f_tol) || (dpm > -dpm_tol && yf <= _f_tol && yf >= -_f_tol);
 }

◆ modelName()

std::string SolidMechanicsPlasticMohrCoulombMulti::modelName ( ) const

overridevirtual

Implements SolidMechanicsPlasticModel.

Definition at line 321 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   return "MohrCoulombMulti";
 }

◆ numberSurfaces()

unsigned int SolidMechanicsPlasticMohrCoulombMulti::numberSurfaces ( ) const

overridevirtual

The number of yield surfaces for this plasticity model.

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 83 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   return 6;
 }

◆ perturbStress()

void SolidMechanicsPlasticMohrCoulombMulti::perturbStress	(	const RankTwoTensor &	stress,
		std::vector< Real > &	eigvals,
		std::vector< RankTwoTensor > &	deigvals
	)		const

private

perturbs the stress tensor in the case of almost-equal eigenvalues.

Note that, upon entry, this algorithm assumes that eigvals are the eigvenvalues of stress

Parameters

	stress	input stress
[in,out]	eigvals	eigenvalues after perturbing.
[in,out]	deigvals	d(eigenvalues)/dstress_ij after perturbing.

Definition at line 124 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by df_dsig(), and dflowPotential_dintnlV().

 {
   Real small_perturbation;
   RankTwoTensor shifted_stress = stress;
   while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
   {
     for (unsigned i = 0; i < 3; ++i)
       for (unsigned j = 0; j <= i; ++j)
       {
         small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
         shifted_stress(i, j) += small_perturbation;
         shifted_stress(j, i) += small_perturbation;
       }
     shifted_stress.dsymmetricEigenvalues(eigvals, deigvals);
   }
 }

◆ phi()

Real SolidMechanicsPlasticMohrCoulombMulti::phi ( const Real internal_param ) const

protectedvirtual

phi as a function of residual value, rate, and internal_param

Definition at line 297 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), dyieldFunction_dintnlV(), dyieldFunction_dstressV(), returnEdge000101(), returnEdge010100(), returnPlane(), returnTip(), and yieldFunctionV().

 {
   return _phi.value(internal_param);
 }

◆ psi()

Real SolidMechanicsPlasticMohrCoulombMulti::psi ( const Real internal_param ) const

protectedvirtual

psi as a function of residual value, rate, and internal_param

Definition at line 309 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by dflowPotential_dintnlV(), dflowPotential_dstressV(), doReturnMap(), and flowPotentialV().

 {
   return _psi.value(internal_param);
 }

◆ returnEdge000101()

bool SolidMechanicsPlasticMohrCoulombMulti::returnEdge000101	(	const std::vector< Real > &	eigvals,
		const std::vector< RealVectorValue > &	n,
		std::vector< Real > &	dpm,
		RankTwoTensor &	returned_stress,
		Real	intnl_old,
		Real &	sinphi,
		Real &	cohcos,
		Real	initial_guess,
		Real	mag_E,
		bool &	nr_converged,
		Real	ep_plastic_tolerance,
		std::vector< Real > &	yf
	)		const

private

Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters

eigvals	The three stress eigenvalues, sorted in ascending order
n	The six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]	The six plastic multipliers resulting from the return-map to the edge
returned_stress[out]	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]	The new value of sin(friction angle) after returning.
cohcos[out]	The new value of cohesion*cos(friction angle) after returning.
mag_E	An approximate value for the magnitude of the Young's modulus. This is used to set appropriate tolerances in the Newton-Raphson procedure
initial_guess	A guess for sum(dpm)
nr_converged[out]	Whether the internal Newton-Raphson process converged
ep_plastic_tolerance	The user-set tolerance on the plastic strain.
yf[out]	The yield functions after return

Definition at line 911 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

 {
   // This returns to the Mohr-Coulomb edge
   // with 000101 being active.  This means that
   // f3=f5=0.  Denoting the returned stress by "a"
   // (in principal stress space), this means that
   // a0=a1 = (2Ccosphi + a2(1+sinphi))/(sinphi-1)
   //
   // Also, a = eigvals - dpm3*n3 - dpm5*n5,
   // where the dpm are the plastic multipliers
   // and the n are the flow directions.
   //
   // Hence we have 5 equations and 5 unknowns (a,
   // dpm3 and dpm5).  Eliminating all unknowns
   // except for x = dpm3+dpm5 gives the residual below
   // (I used mathematica)
 
   Real x = initial_guess;
   sinphi = std::sin(phi(intnl_old + x));
   Real cosphi = std::cos(phi(intnl_old + x));
   Real coh = cohesion(intnl_old + x);
   cohcos = coh * cosphi;
 
   const Real numer_const =
       -n[3](2) * eigvals[0] - n[5](1) * eigvals[0] + n[5](2) * eigvals[0] + n[3](2) * eigvals[1] +
       n[5](0) * eigvals[1] - n[5](2) * eigvals[1] - n[5](0) * eigvals[2] + n[5](1) * eigvals[2] +
       n[3](0) * (-eigvals[1] + eigvals[2]) - n[3](1) * (-eigvals[0] + eigvals[2]);
   const Real numer_coeff1 = 2 * (-n[3](0) + n[3](1) + n[5](0) - n[5](1));
   const Real numer_coeff2 =
       n[5](2) * (eigvals[0] - eigvals[1]) + n[3](2) * (-eigvals[0] + eigvals[1]) +
       n[5](1) * (eigvals[0] + eigvals[2]) + (n[3](0) - n[5](0)) * (eigvals[1] + eigvals[2]) -
       n[3](1) * (eigvals[0] + eigvals[2]);
   Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
   const Real denom_const = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (-n[5](0) + n[5](2)) +
                            n[3](0) * (-n[5](1) + n[5](2));
   const Real denom_coeff = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (n[5](0) + n[5](2)) +
                            n[3](0) * (n[5](1) + n[5](2));
   Real denom = denom_const + denom_coeff * sinphi;
   Real residual = -x + numer / denom;
 
   Real deriv_phi;
   Real deriv_coh;
   Real jacobian;
   const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
   unsigned int iter = 0;
   do
   {
     do
     {
       deriv_phi = dphi(intnl_old + dpm[3]);
       deriv_coh = dcohesion(intnl_old + dpm[3]);
       jacobian = -1 +
                  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
                   numer_coeff2 * cosphi * deriv_phi) /
                      denom -
                  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
       x -= residual / jacobian;
       if (iter > _max_iters) // not converging
       {
         nr_converged = false;
         return false;
       }
 
       sinphi = std::sin(phi(intnl_old + x));
       cosphi = std::cos(phi(intnl_old + x));
       coh = cohesion(intnl_old + x);
       cohcos = coh * cosphi;
       numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
       denom = denom_const + denom_coeff * sinphi;
       residual = -x + numer / denom;
       iter++;
     } while (residual * residual > tol);
 
     // now must ensure that yf[3] and yf[5] are both "zero"
     const Real dpm3minusdpm5 =
         (2.0 * (eigvals[0] - eigvals[1]) + x * (n[3](1) - n[3](0) + n[5](1) - n[5](0))) /
         (n[3](0) - n[3](1) + n[5](1) - n[5](0));
     dpm[3] = (x + dpm3minusdpm5) / 2.0;
     dpm[5] = (x - dpm3minusdpm5) / 2.0;
 
     for (unsigned i = 0; i < 3; ++i)
       returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i) - dpm[5] * n[5](i);
     yieldFunctionEigvals(
         returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
   } while (yf[3] * yf[3] > _f_tol * _f_tol && yf[5] * yf[5] > _f_tol * _f_tol);
 
   // so the NR process converged, but we must
   // check Kuhn-Tucker
   nr_converged = true;
 
   if (dpm[3] < -ep_plastic_tolerance || dpm[5] < -ep_plastic_tolerance)
     // Kuhn-Tucker failure.    This is a potential weak-point: if the user has set _f_tol quite
     // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
     // the solution may be wrong and violate Kuhn-Tucker.
     return false;
 
   if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol)
     // Kuhn-Tucker failure
     return false;
 
   // success
   dpm[0] = dpm[1] = dpm[2] = dpm[4] = 0;
   return true;
 }

◆ returnEdge010100()

bool SolidMechanicsPlasticMohrCoulombMulti::returnEdge010100	(	const std::vector< Real > &	eigvals,
		const std::vector< RealVectorValue > &	n,
		std::vector< Real > &	dpm,
		RankTwoTensor &	returned_stress,
		Real	intnl_old,
		Real &	sinphi,
		Real &	cohcos,
		Real	initial_guess,
		Real	mag_E,
		bool &	nr_converged,
		Real	ep_plastic_tolerance,
		std::vector< Real > &	yf
	)		const

private

Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters

eigvals	The three stress eigenvalues, sorted in ascending order
n	The six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]	The six plastic multipliers resulting from the return-map to the edge
returned_stress[out]	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]	The new value of sin(friction angle) after returning.
cohcos[out]	The new value of cohesion*cos(friction angle) after returning.
mag_E	An approximate value for the magnitude of the Young's modulus. This is used to set appropriate tolerances in the Newton-Raphson procedure
initial_guess	A guess for sum(dpm)
nr_converged[out]	Whether the internal Newton-Raphson process converged
ep_plastic_tolerance	The user-set tolerance on the plastic strain.
yf[out]	The yield functions after return

Definition at line 1028 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

 {
   // This returns to the Mohr-Coulomb edge
   // with 010100 being active.  This means that
   // f1=f3=0.  Denoting the returned stress by "a"
   // (in principal stress space), this means that
   // a1=a2 = (2Ccosphi + a0(1-sinphi))/(sinphi+1)
   //
   // Also, a = eigvals - dpm1*n1 - dpm3*n3,
   // where the dpm are the plastic multipliers
   // and the n are the flow directions.
   //
   // Hence we have 5 equations and 5 unknowns (a,
   // dpm1 and dpm3).  Eliminating all unknowns
   // except for x = dpm1+dpm3 gives the residual below
   // (I used mathematica)
 
   Real x = initial_guess;
   sinphi = std::sin(phi(intnl_old + x));
   Real cosphi = std::cos(phi(intnl_old + x));
   Real coh = cohesion(intnl_old + x);
   cohcos = coh * cosphi;
 
   const Real numer_const = -n[1](2) * eigvals[0] - n[3](1) * eigvals[0] + n[3](2) * eigvals[0] -
                            n[1](0) * eigvals[1] + n[1](2) * eigvals[1] + n[3](0) * eigvals[1] -
                            n[3](2) * eigvals[1] + n[1](0) * eigvals[2] - n[3](0) * eigvals[2] +
                            n[3](1) * eigvals[2] - n[1](1) * (-eigvals[0] + eigvals[2]);
   const Real numer_coeff1 = 2 * (n[1](1) - n[1](2) - n[3](1) + n[3](2));
   const Real numer_coeff2 =
       n[3](2) * (-eigvals[0] - eigvals[1]) + n[1](2) * (eigvals[0] + eigvals[1]) +
       n[3](1) * (eigvals[0] + eigvals[2]) + (n[1](0) - n[3](0)) * (eigvals[1] - eigvals[2]) -
       n[1](1) * (eigvals[0] + eigvals[2]);
   Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
   const Real denom_const = -n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (-n[3](0) + n[3](1)) +
                            n[1](1) * (-n[3](2) + n[3](0));
   const Real denom_coeff =
       n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (n[3](0) + n[3](1)) - n[1](1) * (n[3](0) + n[3](2));
   Real denom = denom_const + denom_coeff * sinphi;
   Real residual = -x + numer / denom;
 
   Real deriv_phi;
   Real deriv_coh;
   Real jacobian;
   const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
   unsigned int iter = 0;
   do
   {
     do
     {
       deriv_phi = dphi(intnl_old + dpm[3]);
       deriv_coh = dcohesion(intnl_old + dpm[3]);
       jacobian = -1 +
                  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
                   numer_coeff2 * cosphi * deriv_phi) /
                      denom -
                  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
       x -= residual / jacobian;
       if (iter > _max_iters) // not converging
       {
         nr_converged = false;
         return false;
       }
 
       sinphi = std::sin(phi(intnl_old + x));
       cosphi = std::cos(phi(intnl_old + x));
       coh = cohesion(intnl_old + x);
       cohcos = coh * cosphi;
       numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
       denom = denom_const + denom_coeff * sinphi;
       residual = -x + numer / denom;
       iter++;
     } while (residual * residual > tol);
 
     // now must ensure that yf[1] and yf[3] are both "zero"
     Real dpm1minusdpm3 =
         (2 * (eigvals[1] - eigvals[2]) + x * (n[1](2) - n[1](1) + n[3](2) - n[3](1))) /
         (n[1](1) - n[1](2) + n[3](2) - n[3](1));
     dpm[1] = (x + dpm1minusdpm3) / 2.0;
     dpm[3] = (x - dpm1minusdpm3) / 2.0;
 
     for (unsigned i = 0; i < 3; ++i)
       returned_stress(i, i) = eigvals[i] - dpm[1] * n[1](i) - dpm[3] * n[3](i);
     yieldFunctionEigvals(
         returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
   } while (yf[1] * yf[1] > _f_tol * _f_tol && yf[3] * yf[3] > _f_tol * _f_tol);
 
   // so the NR process converged, but we must
   // check Kuhn-Tucker
   nr_converged = true;
 
   if (dpm[1] < -ep_plastic_tolerance || dpm[3] < -ep_plastic_tolerance)
     // Kuhn-Tucker failure.    This is a potential weak-point: if the user has set _f_tol quite
     // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
     // the solution may be wrong and violate Kuhn-Tucker
     return false;
 
   if (yf[0] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
     // Kuhn-Tucker failure
     return false;
 
   // success
   dpm[0] = dpm[2] = dpm[4] = dpm[5] = 0;
   return true;
 }

◆ returnMap()

bool SolidMechanicsPlasticMohrCoulombMulti::returnMap	(	const RankTwoTensor &	trial_stress,
		Real	intnl_old,
		const RankFourTensor &	E_ijkl,
		Real	ep_plastic_tolerance,
		RankTwoTensor &	returned_stress,
		Real &	returned_intnl,
		std::vector< Real > &	dpm,
		RankTwoTensor &	delta_dp,
		std::vector< Real > &	yf,
		bool &	trial_stress_inadmissible
	)		const

overridevirtual

Performs a custom return-map.

You may choose to over-ride this in your derived SolidMechanicsPlasticXXXX class, and you may implement the return-map algorithm in any way that suits you. Eg, using a Newton-Raphson approach, or a radial-return, etc. This may also be used as a quick way of ascertaining whether (trial_stress, intnl_old) is in fact admissible.

For over-riding this function, please note the following.

(1) Denoting the return value of the function by "successful_return", the only possible output values should be: (A) trial_stress_inadmissible=false, successful_return=true. That is, (trial_stress, intnl_old) is in fact admissible (in the elastic domain). (B) trial_stress_inadmissible=true, successful_return=false. That is (trial_stress, intnl_old) is inadmissible (outside the yield surface), and you didn't return to the yield surface. (C) trial_stress_inadmissible=true, successful_return=true. That is (trial_stress, intnl_old) is inadmissible (outside the yield surface), but you did return to the yield surface. The default implementation only handles case (A) and (B): it does not attempt to do a return-map algorithm.

(2) you must correctly signal "successful_return" using the return value of this function. Don't assume the calling function will do Kuhn-Tucker checking and so forth!

(3) In cases (A) and (B) you needn't set returned_stress, returned_intnl, delta_dp, or dpm. This is for computational efficiency.

(4) In cases (A) and (B), you MUST place the yield function values at (trial_stress, intnl_old) into yf so the calling function can use this information optimally. You will have already calculated these yield function values, which can be quite expensive, and it's not very optimal for the calling function to have to re-calculate them.

(5) In case (C), you need to set: returned_stress (the returned value of stress) returned_intnl (the returned value of the internal variable) delta_dp (the change in plastic strain) dpm (the plastic multipliers needed to bring about the return) yf (yield function values at the returned configuration)

(Note, if you over-ride returnMap, you will probably want to override consistentTangentOpertor too, otherwise it will default to E_ijkl.)

Parameters

	trial_stress	The trial stress
	intnl_old	Value of the internal parameter
	E_ijkl	Elasticity tensor
	ep_plastic_tolerance	Tolerance defined by the user for the plastic strain
[out]	returned_stress	In case (C): lies on the yield surface after returning and produces the correct plastic strain (normality condition). Otherwise: not defined
[out]	returned_intnl	In case (C): the value of the internal parameter after returning. Otherwise: not defined
[out]	dpm	In case (C): the plastic multipliers needed to bring about the return. Otherwise: not defined
[out]	delta_dp	In case (C): The change in plastic strain induced by the return process. Otherwise: not defined
[out]	yf	In case (C): the yield function at (returned_stress, returned_intnl). Otherwise: the yield function at (trial_stress, intnl_old)
[out]	trial_stress_inadmissible	Should be set to false if the trial_stress is admissible, and true if the trial_stress is inadmissible. This can be used by the calling prorgram

Returns: true if a successful return (or a return-map not needed), false if the trial_stress is inadmissible but the return process failed

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 327 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   if (!_use_custom_returnMap)
     return SolidMechanicsPlasticModel::returnMap(trial_stress,
                                                  intnl_old,
                                                  E_ijkl,
                                                  ep_plastic_tolerance,
                                                  returned_stress,
                                                  returned_intnl,
                                                  dpm,
                                                  delta_dp,
                                                  yf,
                                                  trial_stress_inadmissible);
 
   return doReturnMap(trial_stress,
                      intnl_old,
                      E_ijkl,
                      ep_plastic_tolerance,
                      returned_stress,
                      returned_intnl,
                      dpm,
                      delta_dp,
                      yf,
                      trial_stress_inadmissible);
 }

◆ returnPlane()

bool SolidMechanicsPlasticMohrCoulombMulti::returnPlane	(	const std::vector< Real > &	eigvals,
		const std::vector< RealVectorValue > &	n,
		std::vector< Real > &	dpm,
		RankTwoTensor &	returned_stress,
		Real	intnl_old,
		Real &	sinphi,
		Real &	cohcos,
		Real	initial_guess,
		bool &	nr_converged,
		Real	ep_plastic_tolerance,
		std::vector< Real > &	yf
	)		const

private

Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress) will be junk.

Parameters

eigvals	The three stress eigenvalues, sorted in ascending order
n	The six return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are vectors in principal stress space
dpm[out]	The six plastic multipliers resulting from the return-map to the plane
returned_stress[out]	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]	The new value of sin(friction angle) after returning.
cohcos[out]	The new value of cohesion*cos(friction angle) after returning.
initial_guess	A guess for sum(dpm)
nr_converged[out]	Whether the internal Newton-Raphson process converged
ep_plastic_tolerance	The user-set tolerance on the plastic strain.
yf[out]	The yield functions after return

Definition at line 806 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

 {
   // This returns to the Mohr-Coulomb plane using n[3] (ie 000100)
   //
   // The returned point is defined by the f[3]=0 and
   //    a = eigvals - dpm[3]*n[3]
   // where "a" is the returned point and dpm[3] is the plastic multiplier.
   // This equation is a vector equation in principal stress space.
   // (Kuhn-Tucker also demands that dpm[3]>=0, but we leave checking
   // that condition for the end.)
   // Since f[3]=0, we must have
   // a[2]*(1+sinphi) + a[0]*(-1+sinphi) - 2*coh*cosphi = 0
   // which gives dpm[3] as the solution of
   //     alpha*dpm[3] + eigvals[2] - eigvals[0] + beta*sinphi - 2*coh*cosphi = 0
   // with alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0))*sinphi
   //      beta = eigvals[2] + eigvals[0]
 
   mooseAssert(n.size() == 6,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
               "supplied with n of size 6, whereas yours is "
                   << n.size());
   mooseAssert(dpm.size() == 6,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
               "supplied with dpm of size 6, whereas yours is "
                   << dpm.size());
   mooseAssert(yf.size() == 6,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
               "supplied with yf of size 6, whereas yours is "
                   << yf.size());
 
   dpm[3] = initial_guess;
   sinphi = std::sin(phi(intnl_old + dpm[3]));
   Real cosphi = std::cos(phi(intnl_old + dpm[3]));
   Real coh = cohesion(intnl_old + dpm[3]);
   cohcos = coh * cosphi;
 
   Real alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
   Real deriv_phi;
   Real dalpha;
   const Real beta = eigvals[2] + eigvals[0];
   Real deriv_coh;
 
   Real residual =
       alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos; // this is 2*yf[3]
   Real jacobian;
 
   const Real f_tol2 = Utility::pow<2>(_f_tol);
   unsigned int iter = 0;
   do
   {
     deriv_phi = dphi(intnl_old + dpm[3]);
     dalpha = -(n[3](2) + n[3](0)) * cosphi * deriv_phi;
     deriv_coh = dcohesion(intnl_old + dpm[3]);
     jacobian = alpha + dalpha * dpm[3] + beta * cosphi * deriv_phi - 2.0 * deriv_coh * cosphi +
                2.0 * coh * sinphi * deriv_phi;
 
     dpm[3] -= residual / jacobian;
     if (iter > _max_iters) // not converging
     {
       nr_converged = false;
       return false;
     }
 
     sinphi = std::sin(phi(intnl_old + dpm[3]));
     cosphi = std::cos(phi(intnl_old + dpm[3]));
     coh = cohesion(intnl_old + dpm[3]);
     cohcos = coh * cosphi;
     alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
     residual = alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos;
     iter++;
   } while (residual * residual > f_tol2);
 
   // so the NR process converged, but we must
   // check Kuhn-Tucker
   nr_converged = true;
   if (dpm[3] < -ep_plastic_tolerance)
     // Kuhn-Tucker failure
     return false;
 
   for (unsigned i = 0; i < 3; ++i)
     returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i);
   yieldFunctionEigvals(
       returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
 
   // by construction abs(yf[3]) = abs(residual/2) < _f_tol/2
   if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
     // Kuhn-Tucker failure
     return false;
 
   // success!
   dpm[0] = dpm[1] = dpm[2] = dpm[4] = dpm[5] = 0;
   return true;
 }

◆ returnTip()

bool SolidMechanicsPlasticMohrCoulombMulti::returnTip	(	const std::vector< Real > &	eigvals,
		const std::vector< RealVectorValue > &	n,
		std::vector< Real > &	dpm,
		RankTwoTensor &	returned_stress,
		Real	intnl_old,
		Real &	sinphi,
		Real &	cohcos,
		Real	initial_guess,
		bool &	nr_converged,
		Real	ep_plastic_tolerance,
		std::vector< Real > &	yf
	)		const

private

Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are returned.

Note that you must supply THREE suitale n vectors out of the total of SIX flow directions, and then interpret the THREE dpm values appropriately. The return value is true if the internal Newton-Raphson process has converged and Kuhn-Tucker is satisfied, otherwise it is false If the return value is false and/or nr_converged=false then the "out" parameters (sinphi, cohcos, yf, returned_stress, dpm) will be junk.

Parameters

eigvals	The three stress eigenvalues, sorted in ascending order
n	The three return directions, n=E_ijkl*r. Note this algorithm assumes isotropic elasticity, so these are 3 vectors in principal stress space
dpm[out]	The three plastic multipliers resulting from the return-map to the tip.
returned_stress[out]	The returned stress. This will be diagonal, with the return-mapped eigenvalues in the diagonal positions, sorted in ascending order
intnl_old	The internal parameter at stress=eigvals. This algorithm doesn't form the plastic strain, so you will have to use intnl=intnl_old+sum(dpm) if you need the new internal-parameter value at the returned point.
sinphi[out]	The new value of sin(friction angle) after returning.
cohcos[out]	The new value of cohesion*cos(friction angle) after returning.
initial_guess	A guess for sum(dpm)
nr_converged[out]	Whether the internal Newton-Raphson process converged
ep_plastic_tolerance	The user-set tolerance on the plastic strain.
yf[out]	The yield functions after return

Definition at line 638 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap().

 {
   // This returns to the Mohr-Coulomb tip using the THREE directions
   // given in n, and yields the THREE dpm values.  Note that you
   // must supply THREE suitable n vectors out of the total of SIX
   // flow directions, and then interpret the THREE dpm values appropriately.
   //
   // Eg1.  You supply the flow directions n[0], n[1] and n[3] as
   //       the "n" vectors.  This is return-to-the-tip via 110100.
   //       Then the three returned dpm values will be dpm[0], dpm[1] and dpm[3].
 
   // Eg2.  You supply the flow directions n[1], n[3] and n[5] as
   //       the "n" vectors.  This is return-to-the-tip via 010101.
   //       Then the three returned dpm values will be dpm[1], dpm[3] and dpm[5].
 
   // The returned point is defined by the three yield functions (corresonding
   // to the three supplied flow directions) all being zero.
   // that is, returned_stress = diag(cohcot, cohcot, cohcot), where
   // cohcot = cohesion*cosphi/sinphi
   // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
   // The 3 plastic multipliers, dpm, are defiend by the normality condition
   //     eigvals - cohcot = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
   // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
   //  that checking for the end.)
   // (Remember that these "n" vectors and "dpm" values must be interpreted
   //  differently depending on what you pass into this function.)
   // This is a vector equation with solution (A):
   //   dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip;
   //   dpm[1] = triple(eigvals - cohcot, n[2], n[0])/trip;
   //   dpm[2] = triple(eigvals - cohcot, n[0], n[1])/trip;
   // where trip = triple(n[0], n[1], n[2]).
   // By adding the three components of that solution together
   // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
   // and then our Newton-Raphson only involves one variable (x).
   // In the following, i specialise to the isotropic situation.
 
   mooseAssert(n.size() == 3,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
               "supplied with n of size 3, whereas yours is "
                   << n.size());
   mooseAssert(dpm.size() == 3,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
               "supplied with dpm of size 3, whereas yours is "
                   << dpm.size());
   mooseAssert(yf.size() == 6,
               "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
               "supplied with yf of size 6, whereas yours is "
                   << yf.size());
 
   Real x = initial_guess;
   const Real trip = triple_product(n[0], n[1], n[2]);
   sinphi = std::sin(phi(intnl_old + x));
   Real cosphi = std::cos(phi(intnl_old + x));
   Real coh = cohesion(intnl_old + x);
   cohcos = coh * cosphi;
   Real cohcot = cohcos / sinphi;
 
   if (_cohesion.modelName().compare("Constant") != 0 || _phi.modelName().compare("Constant") != 0)
   {
     // Finding x is expensive.  Therefore
     // although x!=0 for Constant Hardening, solution (A)
     // demonstrates that we don't
     // actually need to know x to find the dpm for
     // Constant Hardening.
     //
     // However, for nontrivial Hardening, the following
     // is necessary
     // cohcot_coeff = [1,1,1].(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
     Real cohcot_coeff =
         (n[0](0) * (n[1](1) - n[1](2) - n[2](1)) + (n[1](2) - n[1](1)) * n[2](0) +
          (n[1](0) - n[1](2)) * n[2](1) + n[0](2) * (n[1](0) - n[1](1) - n[2](0) + n[2](1)) +
          n[0](1) * (n[1](2) - n[1](0) + n[2](0) - n[2](2)) +
          (n[0](0) - n[1](0) + n[1](1)) * n[2](2)) /
         trip;
     // eig_term = eigvals.(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
     Real eig_term = eigvals[0] *
                     (-n[0](2) * n[1](1) + n[0](1) * n[1](2) + n[0](2) * n[2](1) -
                      n[1](2) * n[2](1) - n[0](1) * n[2](2) + n[1](1) * n[2](2)) /
                     trip;
     eig_term += eigvals[1] *
                 (n[0](2) * n[1](0) - n[0](0) * n[1](2) - n[0](2) * n[2](0) + n[1](2) * n[2](0) +
                  n[0](0) * n[2](2) - n[1](0) * n[2](2)) /
                 trip;
     eig_term += eigvals[2] *
                 (n[0](0) * n[1](1) - n[1](1) * n[2](0) + n[0](1) * n[2](0) - n[0](1) * n[1](0) -
                  n[0](0) * n[2](1) + n[1](0) * n[2](1)) /
                 trip;
     // and finally, the equation we want to solve is:
     // x - eig_term + cohcot*cohcot_coeff = 0
     // but i divide by cohcot_coeff so the result has the units of
     // stress, so using _f_tol as a convergence check is reasonable
     eig_term /= cohcot_coeff;
     Real residual = x / cohcot_coeff - eig_term + cohcot;
     Real jacobian;
     Real deriv_phi;
     Real deriv_coh;
     unsigned int iter = 0;
     do
     {
       deriv_phi = dphi(intnl_old + x);
       deriv_coh = dcohesion(intnl_old + x);
       jacobian = 1.0 / cohcot_coeff + deriv_coh * cosphi / sinphi -
                  coh * deriv_phi / Utility::pow<2>(sinphi);
       x += -residual / jacobian;
 
       if (iter > _max_iters) // not converging
       {
         nr_converged = false;
         return false;
       }
 
       sinphi = std::sin(phi(intnl_old + x));
       cosphi = std::cos(phi(intnl_old + x));
       coh = cohesion(intnl_old + x);
       cohcos = coh * cosphi;
       cohcot = cohcos / sinphi;
       residual = x / cohcot_coeff - eig_term + cohcot;
       iter++;
     } while (residual * residual > _f_tol * _f_tol / 100);
   }
 
   // so the NR process converged, but we must
   // calculate the individual dpm values and
   // check Kuhn-Tucker
   nr_converged = true;
   if (x < -3 * ep_plastic_tolerance)
     // obviously at least one of the dpm are < -ep_plastic_tolerance.  No point in proceeding.  This
     // is a potential weak-point: if the user has set _f_tol quite large, and ep_plastic_tolerance
     // quite small, the above NR process will quickly converge, but the solution may be wrong and
     // violate Kuhn-Tucker.
     return false;
 
   // The following is the solution (A) written above
   // (dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip, etc)
   // in the isotropic situation
   RealVectorValue v;
   v(0) = eigvals[0] - cohcot;
   v(1) = eigvals[1] - cohcot;
   v(2) = eigvals[2] - cohcot;
   dpm[0] = triple_product(v, n[1], n[2]) / trip;
   dpm[1] = triple_product(v, n[2], n[0]) / trip;
   dpm[2] = triple_product(v, n[0], n[1]) / trip;
 
   if (dpm[0] < -ep_plastic_tolerance || dpm[1] < -ep_plastic_tolerance ||
       dpm[2] < -ep_plastic_tolerance)
     // Kuhn-Tucker failure.  No point in proceeding
     return false;
 
   // Kuhn-Tucker has succeeded: just need returned_stress and yf values
   // I do not use the dpm to calculate returned_stress, because that
   // might add error (and computational effort), simply:
   returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = cohcot;
   // So by construction the yield functions are all zero
   yf[0] = yf[1] = yf[2] = yf[3] = yf[4] = yf[5] = 0;
   return true;
 }

◆ useCustomCTO()

bool SolidMechanicsPlasticModel::useCustomCTO ( ) const

virtualinherited

Returns false. You will want to override this in your derived class if you write a custom consistent tangent operator function.

Reimplemented in SolidMechanicsPlasticTensileMulti, SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticDruckerPragerHyperbolic, and SolidMechanicsPlasticJ2.

Definition at line 213 of file SolidMechanicsPlasticModel.C.

 {
   return false;
 }

◆ useCustomReturnMap()

bool SolidMechanicsPlasticMohrCoulombMulti::useCustomReturnMap ( ) const

overridevirtual

Returns false. You will want to override this in your derived class if you write a custom returnMap function.

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 1164 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   return _use_custom_returnMap;
 }

◆ validParams()

InputParameters SolidMechanicsPlasticMohrCoulombMulti::validParams ( )

static

Definition at line 24 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   InputParameters params = SolidMechanicsPlasticModel::validParams();
   params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
   params.addRequiredParam<UserObjectName>(
       "cohesion", "A SolidMechanicsHardening UserObject that defines hardening of the cohesion");
   params.addRequiredParam<UserObjectName>("friction_angle",
                                           "A SolidMechanicsHardening UserObject "
                                           "that defines hardening of the "
                                           "friction angle (in radians)");
   params.addRequiredParam<UserObjectName>("dilation_angle",
                                           "A SolidMechanicsHardening UserObject "
                                           "that defines hardening of the "
                                           "dilation angle (in radians)");
   params.addParam<unsigned int>("max_iterations",
                                 10,
                                 "Maximum number of Newton-Raphson iterations "
                                 "allowed in the custom return-map algorithm. "
                                 " For highly nonlinear hardening this may "
                                 "need to be higher than 10.");
   params.addParam<Real>("shift",
                         "Yield surface is shifted by this amount to avoid problems with "
                         "defining derivatives when eigenvalues are equal.  If this is "
                         "larger than f_tol, a warning will be issued.  This may be set "
                         "very small when using the custom returnMap.  Default = f_tol.");
   params.addParam<bool>("use_custom_returnMap",
                         true,
                         "Use a custom return-map algorithm for this "
                         "plasticity model, which may speed up "
                         "computations considerably.  Set to true "
                         "only for isotropic elasticity with no "
                         "hardening of the dilation angle.  In this "
                         "case you may set 'shift' very small.");
 
   return params;
 }

◆ yieldFunction()

Real SolidMechanicsPlasticModel::yieldFunction	(	const RankTwoTensor &	stress,
		Real	intnl
	)		const

protectedvirtualinherited

The following functions are what you should override when building single-plasticity models.

The yield function

Parameters

stress	the stress at which to calculate the yield function
intnl	internal parameter

Returns: the yield function

Reimplemented in SolidMechanicsPlasticMeanCapTC, SolidMechanicsPlasticIsotropicSD, SolidMechanicsPlasticDruckerPrager, SolidMechanicsPlasticJ2, SolidMechanicsPlasticOrthotropic, SolidMechanicsPlasticWeakPlaneShear, SolidMechanicsPlasticWeakPlaneTensile, SolidMechanicsPlasticMohrCoulomb, SolidMechanicsPlasticTensile, SolidMechanicsPlasticDruckerPragerHyperbolic, SolidMechanicsPlasticMeanCap, SolidMechanicsPlasticSimpleTester, and SolidMechanicsPlasticWeakPlaneTensileN.

Definition at line 61 of file SolidMechanicsPlasticModel.C.

Referenced by SolidMechanicsPlasticModel::yieldFunctionV().

 {
   return 0.0;
 }

◆ yieldFunctionEigvals()

void SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionEigvals	(	Real	e0,
		Real	e1,
		Real	e2,
		Real	sinphi,
		Real	cohcos,
		std::vector< Real > &	f
	)		const

private

Calculates the yield functions given the eigenvalues of stress.

Parameters

	e0	Smallest eigenvalue
	e1	Middle eigenvalue
	e2	Largest eigenvalue
	sinphi	sin(friction angle)
	cohcos	cohesion*cos(friction angle)
[out]	f	the yield functions

Definition at line 106 of file SolidMechanicsPlasticMohrCoulombMulti.C.

Referenced by doReturnMap(), returnEdge000101(), returnEdge010100(), returnPlane(), and yieldFunctionV().

 {
   // Naively it seems a shame to have 6 yield functions active instead of just
   // 3.  But 3 won't do.  Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
   // Then to return to the yield surface would require 2 positive plastic multipliers
   // and one negative one.  Boo hoo.
 
   f.resize(6);
   f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
   f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
   f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
   f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
   f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
   f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
 }

◆ yieldFunctionV()

void SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionV	(	const RankTwoTensor &	stress,
		Real	intnl,
		std::vector< Real > &	f
	)		const

overridevirtual

Calculates the yield functions.

Note that for single-surface plasticity you don't want to override this - override the private yieldFunction below

Parameters

	stress	the stress at which to calculate the yield function
	intnl	internal parameter
[out]	f	the yield functions

Reimplemented from SolidMechanicsPlasticModel.

Definition at line 89 of file SolidMechanicsPlasticMohrCoulombMulti.C.

 {
   std::vector<Real> eigvals;
   stress.symmetricEigenvalues(eigvals);
   eigvals[0] += _shift;
   eigvals[2] -= _shift;
 
   const Real sinphi = std::sin(phi(intnl));
   const Real cosphi = std::cos(phi(intnl));
   const Real cohcos = cohesion(intnl) * cosphi;
 
   yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
 }

Member Data Documentation

◆ _cohesion

const SolidMechanicsHardeningModel& SolidMechanicsPlasticMohrCoulombMulti::_cohesion

private

Hardening model for cohesion.

Definition at line 95 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by cohesion(), dcohesion(), and returnTip().

◆ _f_tol

const Real SolidMechanicsPlasticModel::_f_tol

inherited

Tolerance on yield function.

Definition at line 170 of file SolidMechanicsPlasticModel.h.

◆ _ic_tol

const Real SolidMechanicsPlasticModel::_ic_tol

inherited

Tolerance on internal constraint.

Definition at line 173 of file SolidMechanicsPlasticModel.h.

◆ _max_iters

const unsigned int SolidMechanicsPlasticMohrCoulombMulti::_max_iters

private

Maximum Newton-Raphison iterations in the custom returnMap algorithm.

Definition at line 104 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by returnEdge000101(), returnEdge010100(), returnPlane(), and returnTip().

◆ _phi

const SolidMechanicsHardeningModel& SolidMechanicsPlasticMohrCoulombMulti::_phi

private

Hardening model for phi.

Definition at line 98 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by dphi(), phi(), and returnTip().

◆ _psi

const SolidMechanicsHardeningModel& SolidMechanicsPlasticMohrCoulombMulti::_psi

private

Hardening model for psi.

Definition at line 101 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by dpsi(), and psi().

◆ _shift

const Real SolidMechanicsPlasticMohrCoulombMulti::_shift

private

yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalues

Definition at line 107 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by df_dsig(), dflowPotential_dintnlV(), doReturnMap(), dyieldFunction_dintnlV(), perturbStress(), SolidMechanicsPlasticMohrCoulombMulti(), and yieldFunctionV().

◆ _use_custom_returnMap

const bool SolidMechanicsPlasticMohrCoulombMulti::_use_custom_returnMap

private

Whether to use the custom return-map algorithm.

Definition at line 110 of file SolidMechanicsPlasticMohrCoulombMulti.h.

Referenced by returnMap(), and useCustomReturnMap().

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

solid_mechanics/include/userobjects/SolidMechanicsPlasticMohrCoulombMulti.h
solid_mechanics/src/userobjects/SolidMechanicsPlasticMohrCoulombMulti.C

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Static Public Attributes

Protected Member Functions

Static Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ return_type

Constructor & Destructor Documentation

◆ SolidMechanicsPlasticMohrCoulombMulti()

Member Function Documentation

◆ activeConstraints()

◆ cohesion()

◆ consistentTangentOperator()

◆ dcohesion()

◆ df_dsig()

◆ dflowPotential_dintnl()

◆ dflowPotential_dintnlV()

◆ dflowPotential_dstress()

◆ dflowPotential_dstressV()

◆ dhardPotential_dintnl()

◆ dhardPotential_dintnlV()

◆ dhardPotential_dstress()

◆ dhardPotential_dstressV()

◆ doReturnMap()

◆ dphi()

◆ dpsi()

◆ dyieldFunction_dintnl()

◆ dyieldFunction_dintnlV()

◆ dyieldFunction_dstress()

◆ dyieldFunction_dstressV()

◆ execute()

◆ finalize()

◆ flowPotential()

◆ flowPotentialV()

◆ hardPotential()

◆ hardPotentialV()

◆ initialize()

◆ KuhnTuckerOK()

◆ KuhnTuckerSingleSurface()

◆ modelName()

◆ numberSurfaces()

◆ perturbStress()

◆ phi()

◆ psi()

◆ returnEdge000101()

◆ returnEdge010100()

◆ returnMap()

◆ returnPlane()

◆ returnTip()

◆ useCustomCTO()

◆ useCustomReturnMap()

◆ validParams()

◆ yieldFunction()

◆ yieldFunctionEigvals()

◆ yieldFunctionV()

Member Data Documentation

◆ _cohesion

◆ _f_tol

◆ _ic_tol

◆ _max_iters

◆ _phi

◆ _psi

◆ _shift

◆ _use_custom_returnMap