MultiPlasticityDebugger Class Reference

MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in their derivatives, etc. More...

#include <MultiPlasticityDebugger.h>

Inheritance diagram for MultiPlasticityDebugger:

Public Member Functions
	MultiPlasticityDebugger (const MooseObject *moose_object)
void	outputAndCheckDebugParameters ()
	Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly.
void	checkDerivatives ()
	Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
void	checkJacobian (const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
	Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.
void	checkSolution (const RankFourTensor &E_inv)
	Checks that Ax does equal b in the NR procedure.
UserObjectName	getUserObjectName (const std::string &param_name) const
const T &	getUserObject (const std::string &param_name, bool is_dependency=true) const
const T &	getUserObjectByName (const UserObjectName &object_name, bool is_dependency=true) const
const UserObjectBase &	getUserObjectBase (const std::string &param_name, bool is_dependency=true) const
const UserObjectBase &	getUserObjectBaseByName (const UserObjectName &object_name, bool is_dependency=true) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObject (const std::string &param_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const
bool	hasUserObjectByName (const UserObjectName &object_name) const

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual void	calculateConstraints (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
	The constraints.
virtual void	calculateRHS (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
	Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void	calculateJacobian (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real > > &jac)
	d(rhs)/d(dof)
virtual void	nrStep (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
	Performs one Newton-Raphson step.
virtual void	yieldFunction (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
	The active yield function(s)
virtual void	dyieldFunction_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
	The derivative of the active yield function(s) with respect to stress.
virtual void	dyieldFunction_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
	The derivative of active yield function(s) with respect to their internal parameters (the user objects assume there is exactly one internal param per yield function)
virtual void	flowPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
	The active flow potential(s) - one for each yield function.
virtual void	dflowPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
	The derivative of the active flow potential(s) with respect to stress.
virtual void	dflowPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
	The derivative of the active flow potentials with respect to the active internal parameters The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha].
virtual void	hardPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
	The active hardening potentials (one for each internal parameter and for each yield function) by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part of model a, so we only calculate those here.
virtual void	dhardPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
	The derivative of the active hardening potentials with respect to stress By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only calculate those here.
virtual void	dhardPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
	The derivative of the active hardening potentials with respect to the active internal parameters.
virtual void	buildActiveConstraints (const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
	Constructs a set of active constraints, given the yield functions, f.
unsigned int	modelNumber (unsigned int surface)
	returns the model number, given the surface number
bool	anyActiveSurfaces (int model, const std::vector< bool > &active)
	returns true if any internal surfaces of the given model are active according to 'active'
void	activeModelSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
	Returns the internal surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.
void	activeSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
	Returns the external surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.
bool	returnMapAll (const RankTwoTensor &trial_stress, const std::vector< Real > &intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &stress, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, RankTwoTensor &delta_dp, std::vector< Real > &yf, unsigned &num_successful_plastic_returns, unsigned &custom_model)
	Performs a returnMap for each plastic model using their inbuilt returnMap functions.
virtual void	addUserObjectDependencyHelper (const UserObjectBase &) const

Protected Attributes
MooseEnum	_fspb_debug
	none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobian_and_linear_system - check entire jacobian and check that Ax=b
RankTwoTensor	_fspb_debug_stress
	Debug the Jacobian entries at this stress.
std::vector< Real >	_fspb_debug_pm
	Debug the Jacobian entires at these plastic multipliers.
std::vector< Real >	_fspb_debug_intnl
	Debug the Jacobian entires at these internal parameters.
Real	_fspb_debug_stress_change
	Debug finite-differencing parameter for the stress.
std::vector< Real >	_fspb_debug_pm_change
	Debug finite-differencing parameters for the plastic multipliers.
std::vector< Real >	_fspb_debug_intnl_change
	Debug finite-differencing parameters for the internal parameters.
Real	_svd_tol
	Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent.
Real	_min_f_tol
	Minimum value of the _f_tol parameters for the Yield Function User Objects.
const InputParameters &	_params
unsigned int	_num_models
	Number of plastic models for this material.
unsigned int	_num_surfaces
	Number of surfaces within the plastic models.
std::vector< std::vector< unsigned int > >	_surfaces_given_model
	_surfaces_given_model[model_number] = vector of surface numbers for this model
MooseEnum	_specialIC
	Allows initial set of active constraints to be chosen optimally.
std::vector< const SolidMechanicsPlasticModel * >	_f
	User objects that define the yield functions, flow potentials, etc.

Private Member Functions
void	fddyieldFunction_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &df_dstress)
	The finite-difference derivative of yield function(s) with respect to stress.
void	fddyieldFunction_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &df_dintnl)
	The finite-difference derivative of yield function(s) with respect to internal parameter(s)
virtual void	fddflowPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankFourTensor > &dr_dstress)
	The finite-difference derivative of the flow potential(s) with respect to stress.
virtual void	fddflowPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &dr_dintnl)
	The finite-difference derivative of the flow potentials with respect to internal parameters.
virtual void	fdJacobian (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real > > &jac)
	The Jacobian calculated using finite differences.
bool	dof_included (unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
virtual int	singularValuesOfR (const std::vector< RankTwoTensor > &r, std::vector< Real > &s)
	Performs a singular-value decomposition of r and returns the singular values.
virtual void	eliminateLinearDependence (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
	Performs a number of singular-value decompositions to check for linear-dependence of the active directions "r" If linear dependence is found, then deactivated_due_to_ld will contain 'true' entries where surfaces need to be deactivated_due_to_ld.
void	buildActiveConstraintsRock (const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
	"Rock" version Constructs a set of active constraints, given the yield functions, f.
void	buildActiveConstraintsJoint (const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
	"Joint" version Constructs a set of active constraints, given the yield functions, f.
const UserObjectBase &	getUserObjectFromFEProblem (const UserObjectName &object_name, const THREAD_ID tid=0) const
const T &	castUserObject (const UserObjectBase &uo_base, const std::string &param_name="") const
void	mooseObjectError (const std::string &param_name, std::stringstream &oss) const
const std::string &	userObjectType (const UserObjectBase &uo) const
const std::string &	userObjectName (const UserObjectBase &uo) const

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

Detailed Description

MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in their derivatives, etc.

Definition at line 19 of file MultiPlasticityDebugger.h.

Constructor & Destructor Documentation

MultiPlasticityDebugger()

MultiPlasticityDebugger::MultiPlasticityDebugger

(

const MooseObject *

moose_object

)

Definition at line 40 of file MultiPlasticityDebugger.C.

  : MultiPlasticityLinearSystem(moose_object),
    _fspb_debug(_params.get<MooseEnum>("debug_fspb")),
    _fspb_debug_stress(_params.get<RealTensorValue>("debug_jac_at_stress")),
    _fspb_debug_pm(_params.get<std::vector<Real>>("debug_jac_at_pm")),
    _fspb_debug_intnl(_params.get<std::vector<Real>>("debug_jac_at_intnl")),
    _fspb_debug_stress_change(_params.get<Real>("debug_stress_change")),
    _fspb_debug_pm_change(_params.get<std::vector<Real>>("debug_pm_change")),
    _fspb_debug_intnl_change(_params.get<std::vector<Real>>("debug_intnl_change"))
{
}

Member Function Documentation

activeModelSurfaces()

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

protectedinherited

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

Parameters

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

Definition at line 809 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), and MultiPlasticityRawComponentAssembler::yieldFunction().

activeSurfaces()

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

protectedinherited

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

Parameters

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

Definition at line 798 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

anyActiveSurfaces()

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

protectedinherited

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

Definition at line 789 of file MultiPlasticityRawComponentAssembler.C.

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

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

buildActiveConstraints()

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

protectedvirtualinherited

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

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

Parameters

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

Definition at line 342 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by ComputeMultiPlasticityStress::returnMap().

buildActiveConstraintsJoint()

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

privateinherited

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

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

Parameters

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

Definition at line 379 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints().

buildActiveConstraintsRock()

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

privateinherited

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

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

Parameters

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

Definition at line 420 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityRawComponentAssembler::buildActiveConstraints().

calculateConstraints()

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

protectedvirtualinherited

The constraints.

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

Parameters

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

Definition at line 227 of file MultiPlasticityLinearSystem.C.

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

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

calculateJacobian()

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

protectedvirtualinherited

d(rhs)/d(dof)

Definition at line 366 of file MultiPlasticityLinearSystem.C.

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

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

calculateRHS()

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

protectedvirtualinherited

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

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

Parameters

	stress	The stress
	intnl_old	old values of the internal parameters
	intnl	internal parameters
	pm	Current value(s) of the plasticity multiplier(s) (consistency parameters)
	delta_dp	Change in plastic strain incurred so far during the return
[out]	rhs	the rhs
	active	The active constraints.
	eliminate_ld	Check for linear dependence of constraints and put the results into deactivated_due_to_ld. Usually this should be true, but for certain debug operations it should be false
[out]	deactivated_due_to_ld	constraints deactivated due to linear-dependence of flow directions

Definition at line 287 of file MultiPlasticityLinearSystem.C.

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

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

checkDerivatives()

void MultiPlasticityDebugger::checkDerivatives

(

)

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

Definition at line 85 of file MultiPlasticityDebugger.C.

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

Referenced by ComputeMultiPlasticityStress::initQpStatefulProperties(), and ComputeMultiPlasticityStress::plasticStep().

checkJacobian()

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

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

Definition at line 163 of file MultiPlasticityDebugger.C.

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

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::plasticStep().

checkSolution()

void MultiPlasticityDebugger::checkSolution

(

const RankFourTensor &

E_inv

)

Checks that Ax does equal b in the NR procedure.

Definition at line 371 of file MultiPlasticityDebugger.C.

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

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::plasticStep().

dflowPotential_dintnl()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 233 of file MultiPlasticityRawComponentAssembler.C.

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

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

dflowPotential_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 205 of file MultiPlasticityRawComponentAssembler.C.

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

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

dhardPotential_dintnl()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the hardening potentials
	intnl	vector of internal parameters
	active	set of active constraints - only the active derivatives are put into "dh_dintnl"
[out]	dh_dintnl	the derivatives. dh_dintnl[a][alpha][b] = dh[a][alpha]/dintnl[b]. Note that the userobjects assume that there is exactly one internal parameter per yield function, so the derivative is only nonzero for a=alpha=b, so that is all we calculate

Definition at line 315 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

dhardPotential_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 287 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

dof_included()

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

private

Definition at line 353 of file MultiPlasticityDebugger.C.

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

Referenced by fdJacobian().

dyieldFunction_dintnl()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 151 of file MultiPlasticityRawComponentAssembler.C.

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

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

dyieldFunction_dstress()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 123 of file MultiPlasticityRawComponentAssembler.C.

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

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

eliminateLinearDependence()

void MultiPlasticityLinearSystem::eliminateLinearDependence ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl, const std::vector< Real > &  f, const std::vector< RankTwoTensor > &  r, const std::vector< bool > &  active, std::vector< bool > &  deactivated_due_to_ld  )

privatevirtualinherited

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

Parameters

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

Definition at line 106 of file MultiPlasticityLinearSystem.C.

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

Referenced by MultiPlasticityLinearSystem::calculateRHS().

fddflowPotential_dintnl()

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

privatevirtual

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

Parameters

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

Definition at line 612 of file MultiPlasticityDebugger.C.

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

Referenced by checkDerivatives().

fddflowPotential_dstress()

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

privatevirtual

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

Parameters

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

Definition at line 583 of file MultiPlasticityDebugger.C.

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

Referenced by checkDerivatives().

fddyieldFunction_dintnl()

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

private

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

Parameters

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

Definition at line 552 of file MultiPlasticityDebugger.C.

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

Referenced by checkDerivatives().

fddyieldFunction_dstress()

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

private

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

Parameters

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

Definition at line 524 of file MultiPlasticityDebugger.C.

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

Referenced by checkDerivatives().

fdJacobian()

void MultiPlasticityDebugger::fdJacobian ( const RankTwoTensor &  stress, const std::vector< Real > &  intnl_old, const std::vector< Real > &  intnl, const std::vector< Real > &  pm, const RankTwoTensor &  delta_dp, const RankFourTensor &  E_inv, bool  eliminate_ld, std::vector< std::vector< Real > > &  jac  )

privatevirtual

The Jacobian calculated using finite differences.

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

Parameters

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

Definition at line 225 of file MultiPlasticityDebugger.C.

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

Referenced by checkJacobian(), and checkSolution().

flowPotential()

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

protectedvirtualinherited

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

Parameters

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

Definition at line 178 of file MultiPlasticityRawComponentAssembler.C.

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

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

hardPotential()

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

protectedvirtualinherited

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

Parameters

	stress	the stress at which to calculate the hardening potential
	intnl	vector of internal parameters
	active	set of active constraints - only the active hardening potentials are put into "h"
[out]	h	the hardening potentials. h[alpha] = hardening potential for yield fcn alpha (and, by the above assumption we know which hardening parameter, a, this belongs to)

Definition at line 260 of file MultiPlasticityRawComponentAssembler.C.

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

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

modelNumber()

unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface )

protectedinherited

returns the model number, given the surface number

Definition at line 783 of file MultiPlasticityRawComponentAssembler.C.

{
  return _model_given_surface[surface];
}

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

nrStep()

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

protectedvirtualinherited

Performs one Newton-Raphson step.

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

Parameters

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

Definition at line 614 of file MultiPlasticityLinearSystem.C.

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

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

outputAndCheckDebugParameters()

void MultiPlasticityDebugger::outputAndCheckDebugParameters

(

)

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

Definition at line 53 of file MultiPlasticityDebugger.C.

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

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

returnMapAll()

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

protectedinherited

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

Performs a returnMap for each plastic model.

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

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

Parameters

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

Returns

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

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

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

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

Definition at line 597 of file MultiPlasticityRawComponentAssembler.C.

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

Referenced by ComputeMultiPlasticityStress::quickStep().

singularValuesOfR()

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

privatevirtualinherited

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

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

Parameters

	r	The flow directions
[out]	s	The singular values

Returns

The return value from the PETSc LAPACK gesvd reoutine

Definition at line 45 of file MultiPlasticityLinearSystem.C.

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

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

validParams()

InputParameters MultiPlasticityDebugger::validParams ( )

static

Definition at line 15 of file MultiPlasticityDebugger.C.

{
  InputParameters params = MultiPlasticityLinearSystem::validParams();
  MooseEnum debug_fspb_type("none crash jacobian jacobian_and_linear_system", "none");
  params.addParam<MooseEnum>("debug_fspb",
                             debug_fspb_type,
                             "Debug types for use by developers when creating new "
                             "plasticity models, not for general use.  2 = debug Jacobian "
                             "entries, 3 = check the entire Jacobian, and check Ax=b");
  params.addParam<RealTensorValue>("debug_jac_at_stress",
                                   RealTensorValue(),
                                   "Debug Jacobian entries at this stress.  For use by developers");
  params.addParam<std::vector<Real>>("debug_jac_at_pm",
                                     "Debug Jacobian entries at these plastic multipliers");
  params.addParam<std::vector<Real>>("debug_jac_at_intnl",
                                     "Debug Jacobian entries at these internal parameters");
  params.addParam<Real>(
      "debug_stress_change", 1.0, "Debug finite differencing parameter for the stress");
  params.addParam<std::vector<Real>>(
      "debug_pm_change", "Debug finite differencing parameters for the plastic multipliers");
  params.addParam<std::vector<Real>>(
      "debug_intnl_change", "Debug finite differencing parameters for the internal parameters");
  return params;
}

Referenced by ComputeMultiPlasticityStress::validParams().

yieldFunction()

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

protectedvirtualinherited

The active yield function(s)

Parameters

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

Definition at line 96 of file MultiPlasticityRawComponentAssembler.C.

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

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

Member Data Documentation

_f

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

protectedinherited

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

Definition at line 71 of file MultiPlasticityRawComponentAssembler.h.

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

_fspb_debug

MooseEnum MultiPlasticityDebugger::_fspb_debug

protected

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

Definition at line 57 of file MultiPlasticityDebugger.h.

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::initQpStatefulProperties().

_fspb_debug_intnl

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl

protected

Debug the Jacobian entires at these internal parameters.

Definition at line 66 of file MultiPlasticityDebugger.h.

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

_fspb_debug_intnl_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change

protected

Debug finite-differencing parameters for the internal parameters.

Definition at line 75 of file MultiPlasticityDebugger.h.

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

_fspb_debug_pm

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm

protected

Debug the Jacobian entires at these plastic multipliers.

Definition at line 63 of file MultiPlasticityDebugger.h.

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

_fspb_debug_pm_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change

protected

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 72 of file MultiPlasticityDebugger.h.

Referenced by fdJacobian(), and outputAndCheckDebugParameters().

_fspb_debug_stress

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress

protected

Debug the Jacobian entries at this stress.

Definition at line 60 of file MultiPlasticityDebugger.h.

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

_fspb_debug_stress_change

Real MultiPlasticityDebugger::_fspb_debug_stress_change

protected

Debug finite-differencing parameter for the stress.

Definition at line 69 of file MultiPlasticityDebugger.h.

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

_min_f_tol

Real MultiPlasticityLinearSystem::_min_f_tol

protectedinherited

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

Definition at line 131 of file MultiPlasticityLinearSystem.h.

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

_model_given_surface

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

privateinherited

given a surface number, this returns the model number

Definition at line 291 of file MultiPlasticityRawComponentAssembler.h.

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

_model_surface_given_surface

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

privateinherited

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

Definition at line 294 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_num_models

unsigned int MultiPlasticityRawComponentAssembler::_num_models

protectedinherited

Number of plastic models for this material.

Definition at line 53 of file MultiPlasticityRawComponentAssembler.h.

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

_num_surfaces

unsigned int MultiPlasticityRawComponentAssembler::_num_surfaces

protectedinherited

Number of surfaces within the plastic models.

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

Definition at line 62 of file MultiPlasticityRawComponentAssembler.h.

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

_params

const InputParameters& MultiPlasticityRawComponentAssembler::_params

protectedinherited

Definition at line 50 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler().

_specialIC

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC

protectedinherited

Allows initial set of active constraints to be chosen optimally.

Definition at line 68 of file MultiPlasticityRawComponentAssembler.h.

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

_surfaces_given_model

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

protectedinherited

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

Definition at line 65 of file MultiPlasticityRawComponentAssembler.h.

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

_svd_tol

Real MultiPlasticityLinearSystem::_svd_tol

protectedinherited

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

Definition at line 128 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().

---

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

• solid_mechanics/include/utils/MultiPlasticityDebugger.h
• solid_mechanics/src/utils/MultiPlasticityDebugger.C