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

MultiPlasticityLinearSystem computes the linear system and handles linear-dependence removal for use in FiniteStrainMultiPlasticity. More...

#include <MultiPlasticityLinearSystem.h>

Inheritance diagram for MultiPlasticityLinearSystem:
[legend]

Public Member Functions

 MultiPlasticityLinearSystem (const MooseObject *moose_object)
 

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

Protected Attributes

Real _svd_tol
 Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent. More...
 
Real _min_f_tol
 Minimum value of the _f_tol parameters for the Yield Function User Objects. More...
 
const InputParameters & _params
 
unsigned int _num_models
 Number of plastic models for this material. More...
 
unsigned int _num_surfaces
 Number of surfaces within the plastic models. More...
 
std::vector< std::vector< unsigned int > > _surfaces_given_model
 _surfaces_given_model[model_number] = vector of surface numbers for this model More...
 
MooseEnum _specialIC
 Allows initial set of active constraints to be chosen optimally. More...
 
std::vector< const TensorMechanicsPlasticModel * > _f
 User objects that define the yield functions, flow potentials, etc. More...
 

Private Member Functions

virtual int singularValuesOfR (const std::vector< RankTwoTensor > &r, std::vector< Real > &s)
 Performs a singular-value decomposition of r and returns the singular values. More...
 
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. More...
 

Detailed Description

MultiPlasticityLinearSystem computes the linear system and handles linear-dependence removal for use in FiniteStrainMultiPlasticity.

Note that if run in debug mode you might have to use the –no-trap-fpe flag because PETSc-LAPACK-BLAS explicitly compute 0/0 and 1/0, and this causes Libmesh to trap the floating-point exceptions

These routines are quite complicated, so here is an extended explanation.

SURFACES AND MODELS

Each plasticity model can have multiple surfaces (eg, Mohr-Coulomb has 6 surfaces), and one internal parameter. This is also described in MultiPlasticityRawComponentAssembler. The

VARIABLE NAMES

_num_surfaces = total number of surfaces _num_models = total number of plasticity models pm = plasticity multiplier. pm.size() = _num_surfaces intnl = internal variable. intnl.size() = _num_models

DEGREES OF FREEDOM

The degrees of freedom are: the 6 components of stress (it is assumed to be symmetric) the plasticity multipliers, pm. the internal parameters, intnl.

Note that in any single Newton-Raphson (NR) iteration, the number of pm and intnl may be different from any other NR iteration. This is because of deactivating surfaces because they are deemed unimportant by the calling program (eg, their yield function is < 0), or because their flow direction is linearly dependent on other surfaces. Therefore:

Hence, more exactly, the degrees of freedom, whose changes will be provided in the NR step are: the 6 components of stress (it is assumed to be symmetric) the plasticity multipliers, pm, belonging to linearly independent surfaces the internal parameters, intnl, belonging to models with at least one linearly-independent surface

THE CONSTRAINTS AND RHS

The variables calculated by calculateConstraints and calculateRHS are: epp = pm*r - E_inv*(trial_stress - stress) = pm*r - delta_dp f = yield function [all the active constraints, including the deactivated_due_to_ld. The latter ones will not be put into the linear system] ic = intnl - intnl_old + pm*h [only for models that contain active surfaces]

Here pm*r = sum_{active_alpha} pm[alpha]*r[alpha]. Note that this contains all the "active" surfaces, even the ones that have been deactivated_due_to_ld. in calculateConstraints, etc, r is a std::vector containing only all the active flow directions (including deactivated_due_to_ld, but not the "not active"). f = all the "active" surfaces, even the ones that have been deactivated_due_to_ld. However, the latter are not put into the RHS pm*h = sum_{active_alpha} pm[alpha]*h[alpha]. Note that this only contains the "active" hardening potentials, even the ones that have been deactivated_due_to_ld. In calculateConstraints, calculateRHS and calculateJacobian, h is a std::vector containing only these "active" ones. Hence, the sum_{active_alpha} contains deactivated_due_to_ld contributions. HOWEVER, if all the surfaces belonging to a model are either "not active" or deactivated_due_to_ld, then this ic is not included in the RHS

The RHS 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_active_f], ic[0], ic[1], ..., ic[num_active_ic]) Notice the appearance of only the i>=j "epp" components.

THE JACOBIAN

This is d(-rhs)/d(dof). Remember that the dofs are dependent on what is deactivated_due_to_ld, as specified above. In matrix form, the Jacobian is: ( depp_dstress depp_dpm depp_dintnl ) ( df_dstress 0 df_dintnl ) ( dic_dstress dic_dpm dic_dintnl ) For the "epp" terms, only the i>=j components are kept in the RHS, so only these terms are kept here too

Definition at line 125 of file MultiPlasticityLinearSystem.h.

Constructor & Destructor Documentation

◆ MultiPlasticityLinearSystem()

MultiPlasticityLinearSystem::MultiPlasticityLinearSystem ( const MooseObject *  moose_object)

Definition at line 33 of file MultiPlasticityLinearSystem.C.

35  _svd_tol(_params.get<Real>("linear_dependent")),
36  _min_f_tol(-1.0)
37 {
38  for (unsigned model = 0; model < _num_models; ++model)
39  if (_min_f_tol == -1.0 || _min_f_tol > _f[model]->_f_tol)
40  _min_f_tol = _f[model]->_f_tol;
41 
42  MooseRandom::seed(0);
43 }
Real _svd_tol
Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependen...
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.
MultiPlasticityRawComponentAssembler(const MooseObject *moose_object)
Real _min_f_tol
Minimum value of the _f_tol parameters for the Yield Function User Objects.

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
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfaces_of_modelthe output

Definition at line 810 of file MultiPlasticityRawComponentAssembler.C.

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

814 {
815  active_surfaces_of_model.resize(0);
816  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
817  if (active[_surfaces_given_model[model][model_surface]])
818  active_surfaces_of_model.push_back(model_surface);
819 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ 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
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfacesthe output

Definition at line 799 of file MultiPlasticityRawComponentAssembler.C.

Referenced by calculateConstraints().

802 {
803  active_surfaces.resize(0);
804  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
805  if (active[_surfaces_given_model[model][model_surface]])
806  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
807 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ 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 790 of file MultiPlasticityRawComponentAssembler.C.

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

791 {
792  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
793  if (active[_surfaces_given_model[model][model_surface]])
794  return true;
795  return false;
796 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ 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 TensorMechanicsPlasticModel::activeConstraints to identify the active constraints for each model.

Parameters
fyield functions (should be _num_surfaces of these)
stressstress tensor
intnlinternal parameters
Eijklelasticity tensor (stress = Eijkl*strain)
[out]actthe set of active constraints (will be resized to _num_surfaces)

Definition at line 343 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::returnMap().

348 {
349  mooseAssert(f.size() == _num_surfaces,
350  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
351  << _num_surfaces << " surfaces");
352  mooseAssert(intnl.size() == _num_models,
353  "buildActiveConstraints called with intnl.size = "
354  << intnl.size() << " while there are " << _num_models << " models");
355 
356  if (_specialIC == "rock")
357  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
358  else if (_specialIC == "joint")
359  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
360  else // no specialIC
361  {
362  act.resize(0);
363  unsigned ind = 0;
364  for (unsigned model = 0; model < _num_models; ++model)
365  {
366  std::vector<Real> model_f(0);
367  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
368  model_f.push_back(f[ind++]);
369  std::vector<bool> model_act;
370  RankTwoTensor returned_stress;
371  _f[model]->activeConstraints(
372  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
373  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
374  act.push_back(model_act[model_surface]);
375  }
376  }
377 }
MooseEnum _specialIC
Allows initial set of active constraints to be chosen optimally.
void buildActiveConstraintsRock(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Rock" version Constructs a set of active constraints, given the yield functions, f...
void buildActiveConstraintsJoint(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Joint" version Constructs a set of active constraints, given the yield functions, f.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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 
)
protectedvirtual

The constraints.

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

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]fActive yield function(s)
[out]rActive flow directions
[out]eppPlastic-strain increment constraint
[out]icActive internal-parameter constraint
activeThe active constraints.

Definition at line 227 of file MultiPlasticityLinearSystem.C.

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

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

◆ 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 
)
protectedvirtual

d(rhs)/d(dof)

Definition at line 366 of file MultiPlasticityLinearSystem.C.

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

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

◆ 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 
)
protectedvirtual

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

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

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]rhsthe rhs
activeThe active constraints.
eliminate_ldCheck for linear dependence of constraints and put the results into deactivated_due_to_ld. Usually this should be true, but for certain debug operations it should be false
[out]deactivated_due_to_ldconstraints deactivated due to linear-dependence of flow directions

Definition at line 287 of file MultiPlasticityLinearSystem.C.

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

296 {
297  // see comments at the start of .h file
298 
299  mooseAssert(intnl_old.size() == _num_models,
300  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
301  mooseAssert(intnl.size() == _num_models,
302  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
303  mooseAssert(pm.size() == _num_surfaces,
304  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
305  mooseAssert(active.size() == _num_surfaces,
306  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
307 
308  std::vector<Real> f; // the yield functions
309  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
310  std::vector<Real> ic; // the "internal constraints"
311 
312  std::vector<RankTwoTensor> r;
313  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
314 
315  if (eliminate_ld)
316  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
317  else
318  deactivated_due_to_ld.assign(_num_surfaces, false);
319 
320  std::vector<bool> active_not_deact(_num_surfaces);
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
323 
324  unsigned num_active_f = 0;
325  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
326  if (active_not_deact[surface])
327  num_active_f++;
328 
329  unsigned num_active_ic = 0;
330  for (unsigned model = 0; model < _num_models; ++model)
331  if (anyActiveSurfaces(model, active_not_deact))
332  num_active_ic++;
333 
334  unsigned int dim = 3;
335  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
336  // num_active_f comes from "f",
337  // num_active_f comes from "ic"
338 
339  rhs.resize(system_size);
340 
341  unsigned ind = 0;
342  for (unsigned i = 0; i < dim; ++i)
343  for (unsigned j = 0; j <= i; ++j)
344  rhs[ind++] = -epp(i, j);
345  unsigned active_surface = 0;
346  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
347  if (active[surface])
348  {
349  if (!deactivated_due_to_ld[surface])
350  rhs[ind++] = -f[active_surface];
351  active_surface++;
352  }
353  unsigned active_model = 0;
354  for (unsigned model = 0; model < _num_models; ++model)
355  if (anyActiveSurfaces(model, active))
356  {
357  if (anyActiveSurfaces(model, active_not_deact))
358  rhs[ind++] = -ic[active_model];
359  active_model++;
360  }
361 
362  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
363 }
virtual void eliminateLinearDependence(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs a number of singular-value decompositions to check for linear-dependence of the active direc...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.

◆ 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
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dintnl"
[out]dr_dintnlthe derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 234 of file MultiPlasticityRawComponentAssembler.C.

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

238 {
239  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
240  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
241 
242  dr_dintnl.resize(0);
243  std::vector<unsigned int> active_surfaces_of_model;
244  std::vector<unsigned int>::iterator active_surface;
245  std::vector<RankTwoTensor> model_dr_dintnl;
246  for (unsigned model = 0; model < _num_models; ++model)
247  {
248  activeModelSurfaces(model, active, active_surfaces_of_model);
249  if (active_surfaces_of_model.size() > 0)
250  {
251  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
252  for (active_surface = active_surfaces_of_model.begin();
253  active_surface != active_surfaces_of_model.end();
254  ++active_surface)
255  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
256  }
257  }
258 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dstress"
[out]dr_dstressthe derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 206 of file MultiPlasticityRawComponentAssembler.C.

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

211 {
212  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
213  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
214 
215  dr_dstress.resize(0);
216  std::vector<unsigned int> active_surfaces_of_model;
217  std::vector<unsigned int>::iterator active_surface;
218  std::vector<RankFourTensor> model_dr_dstress;
219  for (unsigned model = 0; model < _num_models; ++model)
220  {
221  activeModelSurfaces(model, active, active_surfaces_of_model);
222  if (active_surfaces_of_model.size() > 0)
223  {
224  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
225  for (active_surface = active_surfaces_of_model.begin();
226  active_surface != active_surfaces_of_model.end();
227  ++active_surface)
228  dr_dstress.push_back(model_dr_dstress[*active_surface]);
229  }
230  }
231 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dintnl"
[out]dh_dintnlthe derivatives. dh_dintnl[a][alpha][b] = dh[a][alpha]/dintnl[b]. Note that the userobjects assume that there is exactly one internal parameter per yield function, so the derivative is only nonzero for a=alpha=b, so that is all we calculate

Definition at line 316 of file MultiPlasticityRawComponentAssembler.C.

Referenced by calculateJacobian().

320 {
321  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
322  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
323 
324  dh_dintnl.resize(0);
325  std::vector<unsigned int> active_surfaces_of_model;
326  std::vector<unsigned int>::iterator active_surface;
327  std::vector<Real> model_dh_dintnl;
328  for (unsigned model = 0; model < _num_models; ++model)
329  {
330  activeModelSurfaces(model, active, active_surfaces_of_model);
331  if (active_surfaces_of_model.size() > 0)
332  {
333  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
334  for (active_surface = active_surfaces_of_model.begin();
335  active_surface != active_surfaces_of_model.end();
336  ++active_surface)
337  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
338  }
339  }
340 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dstress"
[out]dh_dstressthe derivative. dh_dstress[a](i, j) = dh[a]/dstress(k, l)

Definition at line 288 of file MultiPlasticityRawComponentAssembler.C.

Referenced by calculateJacobian().

293 {
294  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
295  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
296 
297  dh_dstress.resize(0);
298  std::vector<unsigned int> active_surfaces_of_model;
299  std::vector<unsigned int>::iterator active_surface;
300  std::vector<RankTwoTensor> model_dh_dstress;
301  for (unsigned model = 0; model < _num_models; ++model)
302  {
303  activeModelSurfaces(model, active, active_surfaces_of_model);
304  if (active_surfaces_of_model.size() > 0)
305  {
306  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
307  for (active_surface = active_surfaces_of_model.begin();
308  active_surface != active_surfaces_of_model.end();
309  ++active_surface)
310  dh_dstress.push_back(model_dh_dstress[*active_surface]);
311  }
312  }
313 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dintnl"
[out]df_dintnlthe derivatives. df_dstress[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 152 of file MultiPlasticityRawComponentAssembler.C.

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

156 {
157  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
158  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
159 
160  df_dintnl.resize(0);
161  std::vector<unsigned int> active_surfaces_of_model;
162  std::vector<unsigned int>::iterator active_surface;
163  std::vector<Real> model_df_dintnl;
164  for (unsigned model = 0; model < _num_models; ++model)
165  {
166  activeModelSurfaces(model, active, active_surfaces_of_model);
167  if (active_surfaces_of_model.size() > 0)
168  {
169  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
170  for (active_surface = active_surfaces_of_model.begin();
171  active_surface != active_surfaces_of_model.end();
172  ++active_surface)
173  df_dintnl.push_back(model_df_dintnl[*active_surface]);
174  }
175  }
176 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dstress"
[out]df_dstressthe derivative (or derivatives in the case of multisurface plasticity). df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Definition at line 124 of file MultiPlasticityRawComponentAssembler.C.

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

129 {
130  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
131  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
132 
133  df_dstress.resize(0);
134  std::vector<unsigned int> active_surfaces_of_model;
135  std::vector<unsigned int>::iterator active_surface;
136  std::vector<RankTwoTensor> model_df_dstress;
137  for (unsigned model = 0; model < _num_models; ++model)
138  {
139  activeModelSurfaces(model, active, active_surfaces_of_model);
140  if (active_surfaces_of_model.size() > 0)
141  {
142  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
143  for (active_surface = active_surfaces_of_model.begin();
144  active_surface != active_surfaces_of_model.end();
145  ++active_surface)
146  df_dstress.push_back(model_df_dstress[*active_surface]);
147  }
148  }
149 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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 
)
privatevirtual

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
stressthe current stress
intnlthe current values of internal parameters
fActive yield function values
rthe flow directions that for those yield functions that are active upon entry to this function
activetrue if active
[out]deactivated_due_to_ldYield functions deactivated due to linearly-dependent flow directions

Definition at line 106 of file MultiPlasticityLinearSystem.C.

Referenced by calculateRHS().

112 {
113  deactivated_due_to_ld.resize(_num_surfaces, false);
114 
115  unsigned num_active = r.size();
116 
117  if (num_active <= 1)
118  return;
119 
120  std::vector<double> s;
121  int info = singularValuesOfR(r, s);
122  if (info != 0)
123  mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
124  "returned with error code ",
125  info);
126 
127  // num_lin_dep are the number of linearly dependent
128  // "r vectors", if num_active <= 6
129  unsigned int num_lin_dep = 0;
130 
131  unsigned i = s.size();
132  while (i-- > 0)
133  if (s[i] < _svd_tol * s[0])
134  num_lin_dep++;
135  else
136  break;
137 
138  if (num_lin_dep == 0 && num_active <= 6)
139  return;
140 
141  // From here on, some flow directions are linearly dependent
142 
143  // Find the signed "distance" of the current (stress, internal) configuration
144  // from the yield surfaces. This distance will not be precise, but
145  // i want to preferentially deactivate yield surfaces that are close
146  // to the current stress point.
147  std::vector<RankTwoTensor> df_dstress;
148  dyieldFunction_dstress(stress, intnl, active, df_dstress);
149 
150  typedef std::pair<Real, unsigned> pair_for_sorting;
151  std::vector<pair_for_sorting> dist(num_active);
152  for (unsigned i = 0; i < num_active; ++i)
153  {
154  dist[i].first = f[i] / df_dstress[i].L2norm();
155  dist[i].second = i;
156  }
157  std::sort(dist.begin(), dist.end()); // sorted in ascending order
158 
159  // There is a potential problem when we have equal f[i], for it can give oscillations
160  bool equals_detected = false;
161  for (unsigned i = 0; i < num_active - 1; ++i)
162  if (std::abs(dist[i].first - dist[i + 1].first) < _min_f_tol)
163  {
164  equals_detected = true;
165  dist[i].first += _min_f_tol * (MooseRandom::rand() - 0.5);
166  }
167  if (equals_detected)
168  std::sort(dist.begin(), dist.end()); // sorted in ascending order
169 
170  std::vector<bool> scheduled_for_deactivation;
171  scheduled_for_deactivation.assign(num_active, false);
172 
173  // In the following loop we go through all the flow directions, from
174  // the one with the largest dist, to the one with the smallest dist,
175  // adding them one-by-one into r_tmp. Upon each addition we check
176  // for linear-dependence. if LD is found, we schedule the most
177  // recently added flow direction for deactivation, and pop it
178  // back off r_tmp
179  unsigned current_yf;
180  current_yf = dist[num_active - 1].second;
181  // the one with largest dist
182  std::vector<RankTwoTensor> r_tmp = {r[current_yf]};
183 
184  unsigned num_kept_active = 1;
185  for (unsigned yf_to_try = 2; yf_to_try <= num_active; ++yf_to_try)
186  {
187  current_yf = dist[num_active - yf_to_try].second;
188  if (num_active == 2) // shortcut to we don't have to singularValuesOfR
189  scheduled_for_deactivation[current_yf] = true;
190  else if (num_kept_active >= 6) // shortcut to we don't have to singularValuesOfR: there can
191  // never be > 6 linearly-independent r vectors
192  scheduled_for_deactivation[current_yf] = true;
193  else
194  {
195  r_tmp.push_back(r[current_yf]);
196  info = singularValuesOfR(r_tmp, s);
197  if (info != 0)
198  mooseError("In finding the SVD in the return-map algorithm, the PETSC LAPACK gesvd routine "
199  "returned with error code ",
200  info);
201  if (s[s.size() - 1] < _svd_tol * s[0])
202  {
203  scheduled_for_deactivation[current_yf] = true;
204  r_tmp.pop_back();
205  num_lin_dep--;
206  }
207  else
208  num_kept_active++;
209  if (num_lin_dep == 0 && num_active <= 6)
210  // have taken out all the vectors that were linearly dependent
211  // so no point continuing
212  break;
213  }
214  }
215 
216  unsigned int old_active_number = 0;
217  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
218  if (active[surface])
219  {
220  if (scheduled_for_deactivation[old_active_number])
221  deactivated_due_to_ld[surface] = true;
222  old_active_number++;
223  }
224 }
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 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.
Real _svd_tol
Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependen...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real _min_f_tol
Minimum value of the _f_tol parameters for the Yield Function User Objects.

◆ 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
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active flow potentials are put into "r"
[out]rthe flow potential (flow potentials in the multi-surface case)

Definition at line 179 of file MultiPlasticityRawComponentAssembler.C.

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

183 {
184  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
185  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
186 
187  r.resize(0);
188  std::vector<unsigned int> active_surfaces_of_model;
189  std::vector<unsigned int>::iterator active_surface;
190  std::vector<RankTwoTensor> model_r;
191  for (unsigned model = 0; model < _num_models; ++model)
192  {
193  activeModelSurfaces(model, active, active_surfaces_of_model);
194  if (active_surfaces_of_model.size() > 0)
195  {
196  _f[model]->flowPotentialV(stress, intnl[model], model_r);
197  for (active_surface = active_surfaces_of_model.begin();
198  active_surface != active_surfaces_of_model.end();
199  ++active_surface)
200  r.push_back(model_r[*active_surface]);
201  }
202  }
203 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ 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
stressthe stress at which to calculate the hardening potential
intnlvector of internal parameters
activeset of active constraints - only the active hardening potentials are put into "h"
[out]hthe hardening potentials. h[alpha] = hardening potential for yield fcn alpha (and, by the above assumption we know which hardening parameter, a, this belongs to)

Definition at line 261 of file MultiPlasticityRawComponentAssembler.C.

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

265 {
266  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
267  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
268 
269  h.resize(0);
270  std::vector<unsigned int> active_surfaces_of_model;
271  std::vector<unsigned int>::iterator active_surface;
272  std::vector<Real> model_h;
273  for (unsigned model = 0; model < _num_models; ++model)
274  {
275  activeModelSurfaces(model, active, active_surfaces_of_model);
276  if (active_surfaces_of_model.size() > 0)
277  {
278  _f[model]->hardPotentialV(stress, intnl[model], model_h);
279  for (active_surface = active_surfaces_of_model.begin();
280  active_surface != active_surfaces_of_model.end();
281  ++active_surface)
282  h.push_back(model_h[*active_surface]);
283  }
284  }
285 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

◆ modelNumber()

unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface)
protectedinherited

◆ 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 
)
protectedvirtual

Performs one Newton-Raphson step.

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

Parameters
stressCurrent value of stress
intnl_oldThe internal variables at the previous "time" step
intnlCurrent value of the internal variables
pmCurrent value of the plasticity multipliers (consistency parameters)
E_invinverse of the elasticity tensor
delta_dpCurrent value of the plastic-strain increment (ie plastic_strain - plastic_strain_old)
[out]dstressThe change in stress for a full Newton step
[out]dpmThe change in all plasticity multipliers for a full Newton step
[out]dintnlThe change in all internal variables for a full Newton step
activeThe active constraints
[out]deactivated_due_to_ldThe constraints deactivated due to linear-dependence of the flow directions

Definition at line 614 of file MultiPlasticityLinearSystem.C.

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

625 {
626  // Calculate RHS and Jacobian
627  std::vector<Real> rhs;
628  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
629 
630  std::vector<std::vector<Real>> jac;
631  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
632 
633  // prepare for LAPACKgesv_ routine provided by PETSc
634  int system_size = rhs.size();
635 
636  std::vector<double> a(system_size * system_size);
637  // Fill in the a "matrix" by going down columns
638  unsigned ind = 0;
639  for (int col = 0; col < system_size; ++col)
640  for (int row = 0; row < system_size; ++row)
641  a[ind++] = jac[row][col];
642 
643  int nrhs = 1;
644  std::vector<int> ipiv(system_size);
645  int info;
646  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
647 
648  if (info != 0)
649  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
650  "routine returned with error code ",
651  info);
652 
653  // Extract the results back to dstress, dpm and dintnl
654  std::vector<bool> active_not_deact(_num_surfaces);
655  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
656  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
657 
658  unsigned int dim = 3;
659  ind = 0;
660 
661  for (unsigned i = 0; i < dim; ++i)
662  for (unsigned j = 0; j <= i; ++j)
663  dstress(i, j) = dstress(j, i) = rhs[ind++];
664  dpm.assign(_num_surfaces, 0);
665  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
666  if (active_not_deact[surface])
667  dpm[surface] = rhs[ind++];
668  dintnl.assign(_num_models, 0);
669  for (unsigned model = 0; model < _num_models; ++model)
670  if (anyActiveSurfaces(model, active_not_deact))
671  dintnl[model] = rhs[ind++];
672 
673  mooseAssert(static_cast<int>(ind) == system_size,
674  "Incorrect extracting of changes from NR solution in nrStep");
675 }
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
unsigned int _num_models
Number of plastic models for this material.

◆ 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_stressthe trial stress
intnl_oldthe old values of the internal parameters
E_ijklthe elasticity tensor
ep_plastic_tolerancethe tolerance on the plastic strain
[out]stressis set to trial_stress in case (A) or (C), and the returned value of stress in case (B).
[out]intnlis set to intnl_old in case (A) or (C), and the returned value of intnl in case (B)
[out]pmZero in case (A) or (C), otherwise the plastic multipliers needed to bring about the returnMap in case (B)
[in/out]cumulative_pm cumulative plastic multipliers, updated in case (B), otherwise left untouched
[out]delta_dpis unchanged in case (A) or (C), and is set to the change in plastic strain in case(B)
[out]yfwill contain the yield function values at (stress, intnl)
[out]num_successful_plastic_returnswill 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 598 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::quickStep().

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

◆ singularValuesOfR()

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

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
rThe flow directions
[out]sThe singular values
Returns
The return value from the PETSc LAPACK gesvd reoutine

Definition at line 46 of file MultiPlasticityLinearSystem.C.

Referenced by eliminateLinearDependence().

48 {
49  int bm = r.size();
50  int bn = 6;
51 
52  s.resize(std::min(bm, bn));
53 
54  // prepare for gesvd or gesdd routine provided by PETSc
55  // Want to find the singular values of matrix
56  // ( r[0](0,0) r[0](0,1) r[0](0,2) r[0](1,1) r[0](1,2) r[0](2,2) )
57  // ( r[1](0,0) r[1](0,1) r[1](0,2) r[1](1,1) r[1](1,2) r[1](2,2) )
58  // 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) )
59  // ( r[3](0,0) r[3](0,1) r[3](0,2) r[3](1,1) r[3](1,2) r[3](2,2) )
60  // ( r[4](0,0) r[4](0,1) r[4](0,2) r[4](1,1) r[4](1,2) r[4](2,2) )
61  // bm = 5
62 
63  std::vector<double> a(bm * 6);
64  // Fill in the a "matrix" by going down columns
65  unsigned ind = 0;
66  for (int col = 0; col < 3; ++col)
67  for (int row = 0; row < bm; ++row)
68  a[ind++] = r[row](0, col);
69  for (int col = 3; col < 5; ++col)
70  for (int row = 0; row < bm; ++row)
71  a[ind++] = r[row](1, col - 2);
72  for (int row = 0; row < bm; ++row)
73  a[ind++] = r[row](2, 2);
74 
75  // u and vt are dummy variables because they won't
76  // get referenced due to the "N" and "N" choices
77  int sizeu = 1;
78  std::vector<double> u(sizeu);
79  int sizevt = 1;
80  std::vector<double> vt(sizevt);
81 
82  int sizework = 16 * (bm + 6); // this is above the lowerbound specified in the LAPACK doco
83  std::vector<double> work(sizework);
84 
85  int info;
86 
87  LAPACKgesvd_("N",
88  "N",
89  &bm,
90  &bn,
91  &a[0],
92  &bm,
93  &s[0],
94  &u[0],
95  &sizeu,
96  &vt[0],
97  &sizevt,
98  &work[0],
99  &sizework,
100  &info);
101 
102  return info;
103 }

◆ 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
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active yield functions are put into "f"
[out]fthe yield function (or functions in the case of multisurface plasticity)

Definition at line 97 of file MultiPlasticityRawComponentAssembler.C.

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

101 {
102  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
103  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
104 
105  f.resize(0);
106  std::vector<unsigned int> active_surfaces_of_model;
107  std::vector<unsigned int>::iterator active_surface;
108  std::vector<Real> model_f;
109  for (unsigned model = 0; model < _num_models; ++model)
110  {
111  activeModelSurfaces(model, active, active_surfaces_of_model);
112  if (active_surfaces_of_model.size() > 0)
113  {
114  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
115  for (active_surface = active_surfaces_of_model.begin();
116  active_surface != active_surfaces_of_model.end();
117  ++active_surface)
118  f.push_back(model_f[*active_surface]);
119  }
120  }
121 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
unsigned int _num_models
Number of plastic models for this material.

Member Data Documentation

◆ _f

std::vector<const TensorMechanicsPlasticModel *> MultiPlasticityRawComponentAssembler::_f
protectedinherited

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

Definition at line 73 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(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), ComputeMultiPlasticityStress::quickStep(), ComputeMultiPlasticityStress::residual2(), ComputeMultiPlasticityStress::returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), ComputeMultiPlasticityStress::singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

◆ _min_f_tol

Real MultiPlasticityLinearSystem::_min_f_tol
protected

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

Definition at line 135 of file MultiPlasticityLinearSystem.h.

Referenced by eliminateLinearDependence(), and MultiPlasticityLinearSystem().

◆ _num_models

unsigned int MultiPlasticityRawComponentAssembler::_num_models
protectedinherited

◆ _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 64 of file MultiPlasticityRawComponentAssembler.h.

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

◆ _params

const InputParameters& MultiPlasticityRawComponentAssembler::_params
protectedinherited

◆ _specialIC

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC
protectedinherited

◆ _surfaces_given_model

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

◆ _svd_tol

Real MultiPlasticityLinearSystem::_svd_tol
protected

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

Definition at line 132 of file MultiPlasticityLinearSystem.h.

Referenced by eliminateLinearDependence().


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