www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
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:
[legend]

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. More...
 
void checkDerivatives ()
 Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations. More...
 
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. More...
 
void checkSolution (const RankFourTensor &E_inv)
 Checks that Ax does equal b in the NR procedure. More...
 

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

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 More...
 
RankTwoTensor _fspb_debug_stress
 Debug the Jacobian entries at this stress. More...
 
std::vector< Real > _fspb_debug_pm
 Debug the Jacobian entires at these plastic multipliers. More...
 
std::vector< Real > _fspb_debug_intnl
 Debug the Jacobian entires at these internal parameters. More...
 
Real _fspb_debug_stress_change
 Debug finite-differencing parameter for the stress. More...
 
std::vector< Real > _fspb_debug_pm_change
 Debug finite-differencing parameters for the plastic multipliers. More...
 
std::vector< Real > _fspb_debug_intnl_change
 Debug finite-differencing parameters for the internal parameters. More...
 
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

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. More...
 
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) More...
 
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. More...
 
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. More...
 
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. More...
 
bool dof_included (unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
 

Detailed Description

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

Definition at line 22 of file MultiPlasticityDebugger.h.

Constructor & Destructor Documentation

MultiPlasticityDebugger::MultiPlasticityDebugger ( const MooseObject *  moose_object)

Definition at line 37 of file MultiPlasticityDebugger.C.

38  : MultiPlasticityLinearSystem(moose_object),
39  _fspb_debug(_params.get<MooseEnum>("debug_fspb")),
40  _fspb_debug_stress(_params.get<RealTensorValue>("debug_jac_at_stress")),
41  _fspb_debug_pm(_params.get<std::vector<Real>>("debug_jac_at_pm")),
42  _fspb_debug_intnl(_params.get<std::vector<Real>>("debug_jac_at_intnl")),
43  _fspb_debug_stress_change(_params.get<Real>("debug_stress_change")),
44  _fspb_debug_pm_change(_params.get<std::vector<Real>>("debug_pm_change")),
45  _fspb_debug_intnl_change(_params.get<std::vector<Real>>("debug_intnl_change"))
46 {
47 }
MooseEnum _fspb_debug
none - don&#39;t do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
MultiPlasticityLinearSystem(const MooseObject *moose_object)
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

Member Function Documentation

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

813 {
814  active_surfaces_of_model.resize(0);
815  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816  if (active[_surfaces_given_model[model][model_surface]])
817  active_surfaces_of_model.push_back(model_surface);
818 }
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.
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 798 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

801 {
802  active_surfaces.resize(0);
803  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804  if (active[_surfaces_given_model[model][model_surface]])
805  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 }
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.
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.

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

790 {
791  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792  if (active[_surfaces_given_model[model][model_surface]])
793  return true;
794  return false;
795 }
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.
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 339 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::returnMap().

344 {
345  mooseAssert(f.size() == _num_surfaces,
346  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
347  << _num_surfaces
348  << " surfaces");
349  mooseAssert(intnl.size() == _num_models,
350  "buildActiveConstraints called with intnl.size = " << intnl.size()
351  << " while there are "
352  << _num_models
353  << " models");
354 
355  if (_specialIC == "rock")
356  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357  else if (_specialIC == "joint")
358  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359  else // no specialIC
360  {
361  act.resize(0);
362  unsigned ind = 0;
363  for (unsigned model = 0; model < _num_models; ++model)
364  {
365  std::vector<Real> model_f(0);
366  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367  model_f.push_back(f[ind++]);
368  std::vector<bool> model_act;
369  RankTwoTensor returned_stress;
370  _f[model]->activeConstraints(
371  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373  act.push_back(model_act[model_surface]);
374  }
375  }
376 }
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.
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
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 223 of file MultiPlasticityLinearSystem.C.

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

233 {
234  // see comments at the start of .h file
235 
236  mooseAssert(intnl_old.size() == _num_models,
237  "Size of intnl_old is " << intnl_old.size()
238  << " which is incorrect in calculateConstraints");
239  mooseAssert(intnl.size() == _num_models,
240  "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
241  mooseAssert(pm.size() == _num_surfaces,
242  "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
243  mooseAssert(active.size() == _num_surfaces,
244  "Size of active is " << active.size()
245  << " which is incorrect in calculateConstraints");
246 
247  // yield functions
248  yieldFunction(stress, intnl, active, f);
249 
250  // flow directions and "epp"
251  flowPotential(stress, intnl, active, r);
252  epp = RankTwoTensor();
253  unsigned ind = 0;
254  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
255  if (active[surface])
256  epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
257  epp -= delta_dp;
258 
259  // internal constraints
260  std::vector<Real> h;
261  hardPotential(stress, intnl, active, h);
262  ic.resize(0);
263  ind = 0;
264  std::vector<unsigned int> active_surfaces;
265  std::vector<unsigned int>::iterator active_surface;
266  for (unsigned model = 0; model < _num_models; ++model)
267  {
268  activeSurfaces(model, active, active_surfaces);
269  if (active_surfaces.size() > 0)
270  {
271  // some surfaces are active in this model, so must form an internal constraint
272  ic.push_back(intnl[model] - intnl_old[model]);
273  for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
274  ++active_surface)
275  ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
276  // since it was constructed in the same
277  // manner
278  }
279  }
280 }
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=...
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 362 of file MultiPlasticityLinearSystem.C.

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

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

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

292 {
293  // see comments at the start of .h file
294 
295  mooseAssert(intnl_old.size() == _num_models,
296  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
297  mooseAssert(intnl.size() == _num_models,
298  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
299  mooseAssert(pm.size() == _num_surfaces,
300  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
301  mooseAssert(active.size() == _num_surfaces,
302  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
303 
304  std::vector<Real> f; // the yield functions
305  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
306  std::vector<Real> ic; // the "internal constraints"
307 
308  std::vector<RankTwoTensor> r;
309  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
310 
311  if (eliminate_ld)
312  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
313  else
314  deactivated_due_to_ld.assign(_num_surfaces, false);
315 
316  std::vector<bool> active_not_deact(_num_surfaces);
317  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
318  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
319 
320  unsigned num_active_f = 0;
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  if (active_not_deact[surface])
323  num_active_f++;
324 
325  unsigned num_active_ic = 0;
326  for (unsigned model = 0; model < _num_models; ++model)
327  if (anyActiveSurfaces(model, active_not_deact))
328  num_active_ic++;
329 
330  unsigned int dim = 3;
331  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
332  // num_active_f comes from "f",
333  // num_active_f comes from "ic"
334 
335  rhs.resize(system_size);
336 
337  unsigned ind = 0;
338  for (unsigned i = 0; i < dim; ++i)
339  for (unsigned j = 0; j <= i; ++j)
340  rhs[ind++] = -epp(i, j);
341  unsigned active_surface = 0;
342  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
343  if (active[surface])
344  {
345  if (!deactivated_due_to_ld[surface])
346  rhs[ind++] = -f[active_surface];
347  active_surface++;
348  }
349  unsigned active_model = 0;
350  for (unsigned model = 0; model < _num_models; ++model)
351  if (anyActiveSurfaces(model, active))
352  {
353  if (anyActiveSurfaces(model, active_not_deact))
354  rhs[ind++] = -ic[active_model];
355  active_model++;
356  }
357 
358  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
359 }
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.
void MultiPlasticityDebugger::checkDerivatives ( )

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

Definition at line 80 of file MultiPlasticityDebugger.C.

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

81 {
82  Moose::err
83  << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
85 
86  std::vector<bool> act;
87  act.assign(_num_surfaces, true);
88 
89  Moose::err << "\ndyieldFunction_dstress. Relative L2 norms.\n";
90  std::vector<RankTwoTensor> df_dstress;
91  std::vector<RankTwoTensor> fddf_dstress;
94  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
95  {
96  Moose::err << "surface = " << surface << " Relative L2norm = "
97  << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
98  (df_dstress[surface] + fddf_dstress[surface]).L2norm()
99  << "\n";
100  Moose::err << "Coded:\n";
101  df_dstress[surface].print();
102  Moose::err << "Finite difference:\n";
103  fddf_dstress[surface].print();
104  }
105 
106  Moose::err << "\ndyieldFunction_dintnl.\n";
107  std::vector<Real> df_dintnl;
109  Moose::err << "Coded:\n";
110  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
111  Moose::err << df_dintnl[surface] << " ";
112  Moose::err << "\n";
113  std::vector<Real> fddf_dintnl;
115  Moose::err << "Finite difference:\n";
116  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
117  Moose::err << fddf_dintnl[surface] << " ";
118  Moose::err << "\n";
119 
120  Moose::err << "\ndflowPotential_dstress. Relative L2 norms.\n";
121  std::vector<RankFourTensor> dr_dstress;
122  std::vector<RankFourTensor> fddr_dstress;
125  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
126  {
127  Moose::err << "surface = " << surface << " Relative L2norm = "
128  << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
129  (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
130  << "\n";
131  Moose::err << "Coded:\n";
132  dr_dstress[surface].print();
133  Moose::err << "Finite difference:\n";
134  fddr_dstress[surface].print();
135  }
136 
137  Moose::err << "\ndflowPotential_dintnl. Relative L2 norms.\n";
138  std::vector<RankTwoTensor> dr_dintnl;
139  std::vector<RankTwoTensor> fddr_dintnl;
142  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
143  {
144  Moose::err << "surface = " << surface << " Relative L2norm = "
145  << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
146  (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
147  << "\n";
148  Moose::err << "Coded:\n";
149  dr_dintnl[surface].print();
150  Moose::err << "Finite difference:\n";
151  fddr_dintnl[surface].print();
152  }
153 }
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) ...
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.
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 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 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. ...
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
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...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
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.
Real L2norm(const RankTwoTensor &r2tensor)
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...
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 156 of file MultiPlasticityDebugger.C.

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

158 {
159  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
161 
162  std::vector<bool> act;
163  act.assign(_num_surfaces, true);
164  std::vector<bool> deactivated_due_to_ld;
165  deactivated_due_to_ld.assign(_num_surfaces, false);
166 
167  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
168 
169  std::vector<std::vector<Real>> jac;
170  calculateJacobian(_fspb_debug_stress,
173  E_inv,
174  act,
175  deactivated_due_to_ld,
176  jac);
177 
178  std::vector<std::vector<Real>> fdjac;
179  fdJacobian(_fspb_debug_stress,
180  intnl_old,
183  delta_dp,
184  E_inv,
185  false,
186  fdjac);
187 
188  Real L2_numer = 0;
189  Real L2_denom = 0;
190  for (unsigned row = 0; row < jac.size(); ++row)
191  for (unsigned col = 0; col < jac.size(); ++col)
192  {
193  L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
194  L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
195  }
196  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
197 
198  Moose::err << "\nHand-coded Jacobian:\n";
199  for (unsigned row = 0; row < jac.size(); ++row)
200  {
201  for (unsigned col = 0; col < jac.size(); ++col)
202  Moose::err << jac[row][col] << " ";
203  Moose::err << "\n";
204  }
205 
206  Moose::err << "Finite difference Jacobian:\n";
207  for (unsigned row = 0; row < fdjac.size(); ++row)
208  {
209  for (unsigned col = 0; col < fdjac.size(); ++col)
210  Moose::err << fdjac[row][col] << " ";
211  Moose::err << "\n";
212  }
213 }
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.
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)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
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.
void MultiPlasticityDebugger::checkSolution ( const RankFourTensor &  E_inv)

Checks that Ax does equal b in the NR procedure.

Definition at line 362 of file MultiPlasticityDebugger.C.

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

363 {
364  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
365  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
367 
368  std::vector<bool> act;
369  act.assign(_num_surfaces, true);
370  std::vector<bool> deactivated_due_to_ld;
371  deactivated_due_to_ld.assign(_num_surfaces, false);
372 
373  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
374 
375  std::vector<Real> orig_rhs;
376  calculateRHS(_fspb_debug_stress,
380  delta_dp,
381  orig_rhs,
382  act,
383  true,
384  deactivated_due_to_ld);
385 
386  Moose::err << "\nb = ";
387  for (unsigned i = 0; i < orig_rhs.size(); ++i)
388  Moose::err << orig_rhs[i] << " ";
389  Moose::err << "\n\n";
390 
391  std::vector<std::vector<Real>> jac_coded;
392  calculateJacobian(_fspb_debug_stress,
395  E_inv,
396  act,
397  deactivated_due_to_ld,
398  jac_coded);
399 
400  Moose::err
401  << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
402  Moose::err
403  << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
404  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
405  "linear dependence.\n";
406  Moose::err << "Hand-coded Jacobian:\n";
407  for (unsigned row = 0; row < jac_coded.size(); ++row)
408  {
409  for (unsigned col = 0; col < jac_coded.size(); ++col)
410  Moose::err << jac_coded[row][col] << " ";
411  Moose::err << "\n";
412  }
413 
414  deactivated_due_to_ld.assign(_num_surfaces,
415  false); // this potentially gets changed by nrStep, below
416  RankTwoTensor dstress;
417  std::vector<Real> dpm;
418  std::vector<Real> dintnl;
419  nrStep(_fspb_debug_stress,
423  E_inv,
424  delta_dp,
425  dstress,
426  dpm,
427  dintnl,
428  act,
429  deactivated_due_to_ld);
430 
431  std::vector<bool> active_not_deact(_num_surfaces);
432  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
433  active_not_deact[surface] = !deactivated_due_to_ld[surface];
434 
435  std::vector<Real> x;
436  x.assign(orig_rhs.size(), 0);
437  unsigned ind = 0;
438  for (unsigned i = 0; i < 3; ++i)
439  for (unsigned j = 0; j <= i; ++j)
440  x[ind++] = dstress(i, j);
441  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
442  if (active_not_deact[surface])
443  x[ind++] = dpm[surface];
444  for (unsigned model = 0; model < _num_models; ++model)
445  if (anyActiveSurfaces(model, active_not_deact))
446  x[ind++] = dintnl[model];
447 
448  mooseAssert(ind == orig_rhs.size(),
449  "Incorrect extracting of changes from NR solution in the "
450  "finite-difference checking of nrStep");
451 
452  Moose::err << "\nThis yields x =";
453  for (unsigned i = 0; i < orig_rhs.size(); ++i)
454  Moose::err << x[i] << " ";
455  Moose::err << "\n";
456 
457  std::vector<std::vector<Real>> jac_fd;
458  fdJacobian(_fspb_debug_stress,
462  delta_dp,
463  E_inv,
464  true,
465  jac_fd);
466 
467  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
468  Moose::err << "in order to check that the solution is correct\n";
469  Moose::err << "Finite-difference Jacobian:\n";
470  for (unsigned row = 0; row < jac_fd.size(); ++row)
471  {
472  for (unsigned col = 0; col < jac_fd.size(); ++col)
473  Moose::err << jac_fd[row][col] << " ";
474  Moose::err << "\n";
475  }
476 
477  Real L2_numer = 0;
478  Real L2_denom = 0;
479  for (unsigned row = 0; row < jac_coded.size(); ++row)
480  for (unsigned col = 0; col < jac_coded.size(); ++col)
481  {
482  L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
483  L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
484  }
485  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
486  << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
487 
488  std::vector<Real> fd_times_x;
489  fd_times_x.assign(orig_rhs.size(), 0);
490  for (unsigned row = 0; row < orig_rhs.size(); ++row)
491  for (unsigned col = 0; col < orig_rhs.size(); ++col)
492  fd_times_x[row] += jac_fd[row][col] * x[col];
493 
494  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
495  for (unsigned i = 0; i < orig_rhs.size(); ++i)
496  Moose::err << fd_times_x[i] << " ";
497  Moose::err << "\n";
498  Moose::err << "Recall that b = \n";
499  for (unsigned i = 0; i < orig_rhs.size(); ++i)
500  Moose::err << orig_rhs[i] << " ";
501  Moose::err << "\n";
502 
503  L2_numer = 0;
504  L2_denom = 0;
505  for (unsigned i = 0; i < orig_rhs.size(); ++i)
506  {
507  L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
508  L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
509  }
510  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
511 }
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.
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)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
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.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
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.
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 230 of file MultiPlasticityRawComponentAssembler.C.

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

234 {
235  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
236  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
237 
238  dr_dintnl.resize(0);
239  std::vector<unsigned int> active_surfaces_of_model;
240  std::vector<unsigned int>::iterator active_surface;
241  std::vector<RankTwoTensor> model_dr_dintnl;
242  for (unsigned model = 0; model < _num_models; ++model)
243  {
244  activeModelSurfaces(model, active, active_surfaces_of_model);
245  if (active_surfaces_of_model.size() > 0)
246  {
247  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
248  for (active_surface = active_surfaces_of_model.begin();
249  active_surface != active_surfaces_of_model.end();
250  ++active_surface)
251  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
252  }
253  }
254 }
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.
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 202 of file MultiPlasticityRawComponentAssembler.C.

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

207 {
208  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
209  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
210 
211  dr_dstress.resize(0);
212  std::vector<unsigned int> active_surfaces_of_model;
213  std::vector<unsigned int>::iterator active_surface;
214  std::vector<RankFourTensor> model_dr_dstress;
215  for (unsigned model = 0; model < _num_models; ++model)
216  {
217  activeModelSurfaces(model, active, active_surfaces_of_model);
218  if (active_surfaces_of_model.size() > 0)
219  {
220  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
221  for (active_surface = active_surfaces_of_model.begin();
222  active_surface != active_surfaces_of_model.end();
223  ++active_surface)
224  dr_dstress.push_back(model_dr_dstress[*active_surface]);
225  }
226  }
227 }
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.
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 312 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

316 {
317  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
318  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
319 
320  dh_dintnl.resize(0);
321  std::vector<unsigned int> active_surfaces_of_model;
322  std::vector<unsigned int>::iterator active_surface;
323  std::vector<Real> model_dh_dintnl;
324  for (unsigned model = 0; model < _num_models; ++model)
325  {
326  activeModelSurfaces(model, active, active_surfaces_of_model);
327  if (active_surfaces_of_model.size() > 0)
328  {
329  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
330  for (active_surface = active_surfaces_of_model.begin();
331  active_surface != active_surfaces_of_model.end();
332  ++active_surface)
333  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
334  }
335  }
336 }
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.
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 284 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

289 {
290  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
291  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
292 
293  dh_dstress.resize(0);
294  std::vector<unsigned int> active_surfaces_of_model;
295  std::vector<unsigned int>::iterator active_surface;
296  std::vector<RankTwoTensor> model_dh_dstress;
297  for (unsigned model = 0; model < _num_models; ++model)
298  {
299  activeModelSurfaces(model, active, active_surfaces_of_model);
300  if (active_surfaces_of_model.size() > 0)
301  {
302  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
303  for (active_surface = active_surfaces_of_model.begin();
304  active_surface != active_surfaces_of_model.end();
305  ++active_surface)
306  dh_dstress.push_back(model_dh_dstress[*active_surface]);
307  }
308  }
309 }
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.
bool MultiPlasticityDebugger::dof_included ( unsigned int  dof,
const std::vector< bool > &  deactivated_due_to_ld 
)
private

Definition at line 344 of file MultiPlasticityDebugger.C.

Referenced by fdJacobian().

346 {
347  if (dof < unsigned(6))
348  // these are the stress components
349  return true;
350  unsigned eff_dof = dof - 6;
351  if (eff_dof < _num_surfaces)
352  // these are the plastic multipliers, pm
353  return !deactivated_due_to_ld[eff_dof];
354  eff_dof -= _num_surfaces; // now we know the dof is an intnl parameter
355  std::vector<bool> active_surface(_num_surfaces);
356  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
357  active_surface[surface] = !deactivated_due_to_ld[surface];
358  return anyActiveSurfaces(eff_dof, active_surface);
359 }
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; ...
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 148 of file MultiPlasticityRawComponentAssembler.C.

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

152 {
153  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
154  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
155 
156  df_dintnl.resize(0);
157  std::vector<unsigned int> active_surfaces_of_model;
158  std::vector<unsigned int>::iterator active_surface;
159  std::vector<Real> model_df_dintnl;
160  for (unsigned model = 0; model < _num_models; ++model)
161  {
162  activeModelSurfaces(model, active, active_surfaces_of_model);
163  if (active_surfaces_of_model.size() > 0)
164  {
165  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
166  for (active_surface = active_surfaces_of_model.begin();
167  active_surface != active_surfaces_of_model.end();
168  ++active_surface)
169  df_dintnl.push_back(model_df_dintnl[*active_surface]);
170  }
171  }
172 }
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.
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 120 of file MultiPlasticityRawComponentAssembler.C.

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

125 {
126  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
127  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
128 
129  df_dstress.resize(0);
130  std::vector<unsigned int> active_surfaces_of_model;
131  std::vector<unsigned int>::iterator active_surface;
132  std::vector<RankTwoTensor> model_df_dstress;
133  for (unsigned model = 0; model < _num_models; ++model)
134  {
135  activeModelSurfaces(model, active, active_surfaces_of_model);
136  if (active_surfaces_of_model.size() > 0)
137  {
138  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
139  for (active_surface = active_surfaces_of_model.begin();
140  active_surface != active_surfaces_of_model.end();
141  ++active_surface)
142  df_dstress.push_back(model_df_dstress[*active_surface]);
143  }
144  }
145 }
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.
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
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
[out]dr_dintnlthe derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 602 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

605 {
606  dr_dintnl.resize(_num_surfaces);
607 
608  std::vector<bool> act;
609  act.assign(_num_surfaces, true);
610 
611  std::vector<RankTwoTensor> origr;
612  flowPotential(stress, intnl, act, origr);
613 
614  std::vector<Real> intnlep;
615  intnlep.resize(_num_models);
616  for (unsigned model = 0; model < _num_models; ++model)
617  intnlep[model] = intnl[model];
618  Real ep;
619  std::vector<RankTwoTensor> rep;
620  unsigned int model;
621  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
622  {
623  model = modelNumber(surface);
624  ep = _fspb_debug_intnl_change[model];
625  intnlep[model] += ep;
626  flowPotential(stress, intnlep, act, rep);
627  dr_dintnl[surface] = (rep[surface] - origr[surface]) / ep;
628  intnlep[model] -= ep;
629  }
630 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
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.
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
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
[out]dr_dstressthe derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 573 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

576 {
577  dr_dstress.assign(_num_surfaces, RankFourTensor());
578 
579  std::vector<bool> act;
580  act.assign(_num_surfaces, true);
581 
582  Real ep = _fspb_debug_stress_change;
583  RankTwoTensor stressep;
584  std::vector<RankTwoTensor> rep, rep_minus;
585  for (unsigned i = 0; i < 3; ++i)
586  for (unsigned j = 0; j < 3; ++j)
587  {
588  stressep = stress;
589  // do a central difference
590  stressep(i, j) += ep / 2.0;
591  flowPotential(stressep, intnl, act, rep);
592  stressep(i, j) -= ep;
593  flowPotential(stressep, intnl, act, rep_minus);
594  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
595  for (unsigned k = 0; k < 3; ++k)
596  for (unsigned l = 0; l < 3; ++l)
597  dr_dstress[surface](k, l, i, j) = (rep[surface](k, l) - rep_minus[surface](k, l)) / ep;
598  }
599 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
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.
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
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
[out]df_dintnlthe derivative (or derivatives in the case of multisurface plasticity). df_dintnl[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 542 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

545 {
546  df_dintnl.resize(_num_surfaces);
547 
548  std::vector<bool> act;
549  act.assign(_num_surfaces, true);
550 
551  std::vector<Real> origf;
552  yieldFunction(stress, intnl, act, origf);
553 
554  std::vector<Real> intnlep;
555  intnlep.resize(_num_models);
556  for (unsigned model = 0; model < _num_models; ++model)
557  intnlep[model] = intnl[model];
558  Real ep;
559  std::vector<Real> fep;
560  unsigned int model;
561  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
562  {
563  model = modelNumber(surface);
564  ep = _fspb_debug_intnl_change[model];
565  intnlep[model] += ep;
566  yieldFunction(stress, intnlep, act, fep);
567  df_dintnl[surface] = (fep[surface] - origf[surface]) / ep;
568  intnlep[model] -= ep;
569  }
570 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
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.
unsigned int _num_models
Number of plastic models for this material.
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
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
[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 514 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

517 {
518  df_dstress.assign(_num_surfaces, RankTwoTensor());
519 
520  std::vector<bool> act;
521  act.assign(_num_surfaces, true);
522 
523  Real ep = _fspb_debug_stress_change;
524  RankTwoTensor stressep;
525  std::vector<Real> fep, fep_minus;
526  for (unsigned i = 0; i < 3; ++i)
527  for (unsigned j = 0; j < 3; ++j)
528  {
529  stressep = stress;
530  // do a central difference to attempt to capture discontinuities
531  // such as those encountered in tensile and Mohr-Coulomb
532  stressep(i, j) += ep / 2.0;
533  yieldFunction(stressep, intnl, act, fep);
534  stressep(i, j) -= ep;
535  yieldFunction(stressep, intnl, act, fep_minus);
536  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
537  df_dstress[surface](i, j) = (fep[surface] - fep_minus[surface]) / ep;
538  }
539 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
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.
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
stressthe stress at which to calculate the Jacobian
intnl_oldthe old values of internal variables (jacobian is inependent of these, but they are needed to do the finite-differencing cleanly)
intnlthe vector of internal parameters at which to calculate the Jacobian
pmthe plasticity multipliers at which to calculate the Jacobian
delta_dpplastic_strain - plastic_strain_old (Jacobian is independent of this, but it is needed to do the finite-differencing cleanly)
E_invinverse of the elasticity tensor
eliminate_ldonly calculate the Jacobian for the linearly independent constraints
[out]jacthe finite-difference Jacobian

Definition at line 216 of file MultiPlasticityDebugger.C.

Referenced by checkJacobian(), and checkSolution().

224 {
225  std::vector<bool> active;
226  active.assign(_num_surfaces, true);
227 
228  std::vector<bool> deactivated_due_to_ld;
229  std::vector<bool> deactivated_due_to_ld_ep;
230 
231  std::vector<Real> orig_rhs;
232  calculateRHS(stress,
233  intnl_old,
234  intnl,
235  pm,
236  delta_dp,
237  orig_rhs,
238  active,
239  eliminate_ld,
240  deactivated_due_to_ld); // this calculates RHS, and also set deactivated_due_to_ld.
241  // The latter stays fixed for the rest of this routine
242 
243  unsigned int whole_system_size = 6 + _num_surfaces + _num_models;
244  unsigned int system_size =
245  orig_rhs.size(); // will be = whole_system_size if eliminate_ld = false, since all active=true
246  jac.resize(system_size);
247  for (unsigned row = 0; row < system_size; ++row)
248  jac[row].assign(system_size, 0);
249 
250  std::vector<Real> rhs_ep;
251  unsigned col = 0;
252 
253  RankTwoTensor stressep;
254  RankTwoTensor delta_dpep;
255  Real ep = _fspb_debug_stress_change;
256  for (unsigned i = 0; i < 3; ++i)
257  for (unsigned j = 0; j <= i; ++j)
258  {
259  stressep = stress;
260  stressep(i, j) += ep;
261  if (i != j)
262  stressep(j, i) += ep;
263  delta_dpep = delta_dp;
264  for (unsigned k = 0; k < 3; ++k)
265  for (unsigned l = 0; l < 3; ++l)
266  {
267  delta_dpep(k, l) -= E_inv(k, l, i, j) * ep;
268  if (i != j)
269  delta_dpep(k, l) -= E_inv(k, l, j, i) * ep;
270  }
271  active.assign(_num_surfaces, true);
272  calculateRHS(stressep,
273  intnl_old,
274  intnl,
275  pm,
276  delta_dpep,
277  rhs_ep,
278  active,
279  false,
280  deactivated_due_to_ld_ep);
281  unsigned row = 0;
282  for (unsigned dof = 0; dof < whole_system_size; ++dof)
283  if (dof_included(dof, deactivated_due_to_ld))
284  {
285  jac[row][col] =
286  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
287  row++;
288  }
289  col++; // all of the first 6 columns are dof_included since they're stresses
290  }
291 
292  std::vector<Real> pmep;
293  pmep.resize(_num_surfaces);
294  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
295  pmep[surface] = pm[surface];
296  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
297  {
298  if (!dof_included(6 + surface, deactivated_due_to_ld))
299  continue;
300  ep = _fspb_debug_pm_change[surface];
301  pmep[surface] += ep;
302  active.assign(_num_surfaces, true);
303  calculateRHS(
304  stress, intnl_old, intnl, pmep, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
305  unsigned row = 0;
306  for (unsigned dof = 0; dof < whole_system_size; ++dof)
307  if (dof_included(dof, deactivated_due_to_ld))
308  {
309  jac[row][col] =
310  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
311  row++;
312  }
313  pmep[surface] -= ep;
314  col++;
315  }
316 
317  std::vector<Real> intnlep;
318  intnlep.resize(_num_models);
319  for (unsigned model = 0; model < _num_models; ++model)
320  intnlep[model] = intnl[model];
321  for (unsigned model = 0; model < _num_models; ++model)
322  {
323  if (!dof_included(6 + _num_surfaces + model, deactivated_due_to_ld))
324  continue;
325  ep = _fspb_debug_intnl_change[model];
326  intnlep[model] += ep;
327  active.assign(_num_surfaces, true);
328  calculateRHS(
329  stress, intnl_old, intnlep, pm, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
330  unsigned row = 0;
331  for (unsigned dof = 0; dof < whole_system_size; ++dof)
332  if (dof_included(dof, deactivated_due_to_ld))
333  {
334  jac[row][col] =
335  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
336  row++;
337  }
338  intnlep[model] -= ep;
339  col++;
340  }
341 }
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])
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool dof_included(unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
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 175 of file MultiPlasticityRawComponentAssembler.C.

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

179 {
180  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
181  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
182 
183  r.resize(0);
184  std::vector<unsigned int> active_surfaces_of_model;
185  std::vector<unsigned int>::iterator active_surface;
186  std::vector<RankTwoTensor> model_r;
187  for (unsigned model = 0; model < _num_models; ++model)
188  {
189  activeModelSurfaces(model, active, active_surfaces_of_model);
190  if (active_surfaces_of_model.size() > 0)
191  {
192  _f[model]->flowPotentialV(stress, intnl[model], model_r);
193  for (active_surface = active_surfaces_of_model.begin();
194  active_surface != active_surfaces_of_model.end();
195  ++active_surface)
196  r.push_back(model_r[*active_surface]);
197  }
198  }
199 }
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.
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 257 of file MultiPlasticityRawComponentAssembler.C.

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

261 {
262  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
263  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
264 
265  h.resize(0);
266  std::vector<unsigned int> active_surfaces_of_model;
267  std::vector<unsigned int>::iterator active_surface;
268  std::vector<Real> model_h;
269  for (unsigned model = 0; model < _num_models; ++model)
270  {
271  activeModelSurfaces(model, active, active_surfaces_of_model);
272  if (active_surfaces_of_model.size() > 0)
273  {
274  _f[model]->hardPotentialV(stress, intnl[model], model_h);
275  for (active_surface = active_surfaces_of_model.begin();
276  active_surface != active_surfaces_of_model.end();
277  ++active_surface)
278  h.push_back(model_h[*active_surface]);
279  }
280  }
281 }
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.
unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface)
protectedinherited
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
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 610 of file MultiPlasticityLinearSystem.C.

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

621 {
622  // Calculate RHS and Jacobian
623  std::vector<Real> rhs;
624  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
625 
626  std::vector<std::vector<Real>> jac;
627  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
628 
629  // prepare for LAPACKgesv_ routine provided by PETSc
630  int system_size = rhs.size();
631 
632  std::vector<double> a(system_size * system_size);
633  // Fill in the a "matrix" by going down columns
634  unsigned ind = 0;
635  for (int col = 0; col < system_size; ++col)
636  for (int row = 0; row < system_size; ++row)
637  a[ind++] = jac[row][col];
638 
639  int nrhs = 1;
640  std::vector<int> ipiv(system_size);
641  int info;
642  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
643 
644  if (info != 0)
645  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
646  "routine returned with error code ",
647  info);
648 
649  // Extract the results back to dstress, dpm and dintnl
650  std::vector<bool> active_not_deact(_num_surfaces);
651  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
652  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
653 
654  unsigned int dim = 3;
655  ind = 0;
656 
657  for (unsigned i = 0; i < dim; ++i)
658  for (unsigned j = 0; j <= i; ++j)
659  dstress(i, j) = dstress(j, i) = rhs[ind++];
660  dpm.assign(_num_surfaces, 0);
661  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
662  if (active_not_deact[surface])
663  dpm[surface] = rhs[ind++];
664  dintnl.assign(_num_models, 0);
665  for (unsigned model = 0; model < _num_models; ++model)
666  if (anyActiveSurfaces(model, active_not_deact))
667  dintnl[model] = rhs[ind++];
668 
669  mooseAssert(static_cast<int>(ind) == system_size,
670  "Incorrect extracting of changes from NR solution in nrStep");
671 }
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.
void MultiPlasticityDebugger::outputAndCheckDebugParameters ( )

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

Definition at line 50 of file MultiPlasticityDebugger.C.

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

51 {
52  Moose::err << "Debug Parameters are as follows\n";
53  Moose::err << "stress = \n";
54  _fspb_debug_stress.print();
55 
56  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
59  mooseError("The debug parameters have the wrong size\n");
60 
61  Moose::err << "plastic multipliers =\n";
62  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
63  Moose::err << _fspb_debug_pm[surface] << "\n";
64 
65  Moose::err << "internal parameters =\n";
66  for (unsigned model = 0; model < _num_models; ++model)
67  Moose::err << _fspb_debug_intnl[model] << "\n";
68 
69  Moose::err << "finite-differencing parameter for stress-changes:\n"
70  << _fspb_debug_stress_change << "\n";
71  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
72  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
73  Moose::err << _fspb_debug_pm_change[surface] << "\n";
74  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
75  for (unsigned model = 0; model < _num_models; ++model)
76  Moose::err << _fspb_debug_intnl_change[model] << "\n";
77 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
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)

Definition at line 597 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::quickStep().

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

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

97 {
98  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
99  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
100 
101  f.resize(0);
102  std::vector<unsigned int> active_surfaces_of_model;
103  std::vector<unsigned int>::iterator active_surface;
104  std::vector<Real> model_f;
105  for (unsigned model = 0; model < _num_models; ++model)
106  {
107  activeModelSurfaces(model, active, active_surfaces_of_model);
108  if (active_surfaces_of_model.size() > 0)
109  {
110  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
111  for (active_surface = active_surfaces_of_model.begin();
112  active_surface != active_surfaces_of_model.end();
113  ++active_surface)
114  f.push_back(model_f[*active_surface]);
115  }
116  }
117 }
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

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

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

Definition at line 65 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().

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 58 of file MultiPlasticityDebugger.h.

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

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl
protected

Debug the Jacobian entires at these internal parameters.

Definition at line 67 of file MultiPlasticityDebugger.h.

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

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change
protected

Debug finite-differencing parameters for the internal parameters.

Definition at line 76 of file MultiPlasticityDebugger.h.

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

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm
protected

Debug the Jacobian entires at these plastic multipliers.

Definition at line 64 of file MultiPlasticityDebugger.h.

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

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change
protected

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 73 of file MultiPlasticityDebugger.h.

Referenced by fdJacobian(), and outputAndCheckDebugParameters().

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress
protected

Debug the Jacobian entries at this stress.

Definition at line 61 of file MultiPlasticityDebugger.h.

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

Real MultiPlasticityDebugger::_fspb_debug_stress_change
protected

Debug finite-differencing parameter for the stress.

Definition at line 70 of file MultiPlasticityDebugger.h.

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

Real MultiPlasticityLinearSystem::_min_f_tol
protectedinherited

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

Definition at line 132 of file MultiPlasticityLinearSystem.h.

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

unsigned int MultiPlasticityRawComponentAssembler::_num_models
protectedinherited
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 56 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().

const InputParameters& MultiPlasticityRawComponentAssembler::_params
protectedinherited
MooseEnum MultiPlasticityRawComponentAssembler::_specialIC
protectedinherited
std::vector<std::vector<unsigned int> > MultiPlasticityRawComponentAssembler::_surfaces_given_model
protectedinherited
Real MultiPlasticityLinearSystem::_svd_tol
protectedinherited

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

Definition at line 129 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().


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