www.mooseframework.org
ComputeMultiPlasticityStress.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
9 #include "MatrixTools.h"
10 
11 #include "MooseException.h"
12 #include "RotationMatrix.h" // for rotVecToZ
13 
14 #include "libmesh/utility.h"
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<ComputeStressBase>();
22  params.addClassDescription("Base class for multi-surface finite-strain plasticity");
23  params.addRangeCheckedParam<unsigned int>("max_NR_iterations",
24  20,
25  "max_NR_iterations>0",
26  "Maximum number of Newton-Raphson iterations allowed");
27  params.addRequiredParam<Real>("ep_plastic_tolerance",
28  "The Newton-Raphson process is only deemed "
29  "converged if the plastic strain increment "
30  "constraints have L2 norm less than this.");
31  params.addRangeCheckedParam<Real>(
32  "min_stepsize",
33  0.01,
34  "min_stepsize>0 & min_stepsize<=1",
35  "If ordinary Newton-Raphson + line-search fails, then the applied strain increment is "
36  "subdivided, and the return-map is tried again. This parameter is the minimum fraction of "
37  "applied strain increment that may be applied before the algorithm gives up entirely");
38  params.addRangeCheckedParam<Real>("max_stepsize_for_dumb",
39  0.01,
40  "max_stepsize_for_dumb>0 & max_stepsize_for_dumb<=1",
41  "If your deactivation_scheme is 'something_to_dumb', then "
42  "'dumb' will only be used if the stepsize falls below this "
43  "value. This parameter is useful because the 'dumb' scheme is "
44  "computationally expensive");
45  MooseEnum deactivation_scheme("optimized safe dumb optimized_to_safe safe_to_dumb "
46  "optimized_to_safe_to_dumb optimized_to_dumb",
47  "optimized");
48  params.addParam<MooseEnum>(
49  "deactivation_scheme",
50  deactivation_scheme,
51  "Scheme by which constraints are deactivated. (NOTE: This is irrelevant if there is only "
52  "one yield surface.) safe: return to the yield surface and then deactivate constraints with "
53  "negative plasticity multipliers. optimized: deactivate a constraint as soon as its "
54  "plasticity multiplier becomes negative. dumb: iteratively try all combinations of active "
55  "constraints until the solution is found. You may specify fall-back options. Eg "
56  "optimized_to_safe: first use 'optimized', and if that fails, try the return with 'safe'.");
57  params.addParam<RealVectorValue>(
58  "transverse_direction",
59  "If this parameter is provided, before the return-map algorithm is "
60  "called a rotation is performed so that the 'z' axis in the new "
61  "frame lies along the transverse_direction in the original frame. "
62  "After returning, the inverse rotation is performed. The "
63  "transverse_direction will itself rotate with large strains. This "
64  "is so that transversely-isotropic plasticity models may be easily "
65  "defined in the frame where the isotropy holds in the x-y plane.");
66  params.addParam<bool>("ignore_failures",
67  false,
68  "The return-map algorithm will return with the best admissible "
69  "stresses and internal parameters that it can, even if they don't "
70  "fully correspond to the applied strain increment. To speed "
71  "computations, this flag can be set to true, the max_NR_iterations "
72  "set small, and the min_stepsize large.");
73  MooseEnum tangent_operator("elastic linear nonlinear", "nonlinear");
74  params.addParam<MooseEnum>("tangent_operator",
75  tangent_operator,
76  "Type of tangent operator to return. 'elastic': return the "
77  "elasticity tensor. 'linear': return the consistent tangent operator "
78  "that is correct for plasticity with yield functions linear in "
79  "stress. 'nonlinear': return the full, general consistent tangent "
80  "operator. The calculations assume the hardening potentials are "
81  "independent of stress and hardening parameters.");
82  params.addParam<bool>("perform_finite_strain_rotations",
83  true,
84  "Tensors are correctly rotated in "
85  "finite-strain simulations. For "
86  "optimal performance you can set "
87  "this to 'false' if you are only "
88  "ever using small strains");
89  params.addClassDescription("Material for multi-surface finite-strain plasticity");
90  return params;
91 }
92 
94  : ComputeStressBase(parameters),
96  _max_iter(getParam<unsigned int>("max_NR_iterations")),
97  _min_stepsize(getParam<Real>("min_stepsize")),
98  _max_stepsize_for_dumb(getParam<Real>("max_stepsize_for_dumb")),
99  _ignore_failures(getParam<bool>("ignore_failures")),
100 
101  _tangent_operator_type((TangentOperatorEnum)(int)getParam<MooseEnum>("tangent_operator")),
102 
103  _epp_tol(getParam<Real>("ep_plastic_tolerance")),
104 
105  _dummy_pm(0),
106 
107  _cumulative_pm(0),
108 
109  _deactivation_scheme((DeactivationSchemeEnum)(int)getParam<MooseEnum>("deactivation_scheme")),
110 
111  _n_supplied(parameters.isParamValid("transverse_direction")),
112  _n_input(_n_supplied ? getParam<RealVectorValue>("transverse_direction") : RealVectorValue()),
113  _rot(RealTensorValue()),
114 
115  _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
116 
117  _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
118  _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
119  _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
120  _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
121  _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
122  _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
123  // for visualisation i convert it to Real
124  _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
125  // boolean, but for
126  // visualisation i
127  // convert it to Real
128  _ld_encountered(declareProperty<Real>(
129  "plastic_linear_dependence_encountered")), // this is really a boolean, but for
130  // visualisation i convert it to Real
131  _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
132  // boolean, but for
133  // visualisation i
134  // convert it to Real
135  _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
136  _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
137 
138  _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
139  _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
140  _rotation_increment(
141  getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
142 
143  _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
144  _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
145 
146  // TODO: This design does NOT work. It makes these materials construction order dependent and it
147  // disregards block restrictions.
148  _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
149  hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
150  _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : NULL),
151  _elastic_flexural_rigidity_tensor(
152  _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
153  : NULL),
154  _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : NULL),
155  _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress") : NULL),
156  _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
157  : NULL),
158 
159  _my_elasticity_tensor(RankFourTensor()),
160  _my_strain_increment(RankTwoTensor()),
161  _my_flexural_rigidity_tensor(RankFourTensor()),
162  _my_curvature(RankTwoTensor()),
163  _step_one(
164  declareRestartableData<bool>("step_one", true)) // InitialStress Deprecation: remove this
165 {
166  if (_epp_tol <= 0)
167  mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
168 
169  if (_n_supplied)
170  {
171  // normalise the inputted transverse_direction
172  if (_n_input.norm() == 0)
173  mooseError(
174  "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
175  else
176  _n_input /= _n_input.norm();
177  }
178 
179  if (_num_surfaces == 1)
181 }
182 
183 void
185 {
187 
188  _plastic_strain[_qp].zero();
189 
190  _intnl[_qp].assign(_num_models, 0);
191 
192  _yf[_qp].assign(_num_surfaces, 0);
193 
194  _dummy_pm.assign(_num_surfaces, 0);
195 
196  _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
197  _linesearch_needed[_qp] = 0;
198  _ld_encountered[_qp] = 0;
199  _constraints_added[_qp] = 0;
200 
201  _n[_qp] = _n_input;
202 
203  if (_cosserat)
204  (*_couple_stress)[_qp].zero();
205 
206  if (_fspb_debug == "jacobian")
207  {
209  mooseError("Finite-differencing completed. Exiting with no error");
210  }
211 }
212 
213 void
215 {
216  // the following "_my" variables can get rotated by preReturnMap and postReturnMap
219  if (_cosserat)
220  {
221  _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
222  _my_curvature = (*_curvature)[_qp];
223  }
224 
225  if (_fspb_debug == "jacobian_and_linear_system")
226  {
227  // cannot do this at initQpStatefulProperties level since E_ijkl is not defined
228  checkJacobian(_elasticity_tensor[_qp].invSymm(), _intnl_old[_qp]);
229  checkSolution(_elasticity_tensor[_qp].invSymm());
230  mooseError("Finite-differencing completed. Exiting with no error");
231  }
232 
233  preReturnMap(); // do rotations to new frame if necessary
234 
235  unsigned int number_iterations;
236  bool linesearch_needed = false;
237  bool ld_encountered = false;
238  bool constraints_added = false;
239 
240  _cumulative_pm.assign(_num_surfaces, 0);
241  // try a "quick" return first - this can be purely elastic, or a customised plastic return defined
242  // by a TensorMechanicsPlasticXXXX UserObject
243  const bool found_solution = quickStep(rot(_stress_old[_qp]),
244  _stress[_qp],
245  _intnl_old[_qp],
246  _intnl[_qp],
247  _dummy_pm,
249  rot(_plastic_strain_old[_qp]),
250  _plastic_strain[_qp],
253  _yf[_qp],
254  number_iterations,
255  _Jacobian_mult[_qp],
257  true);
258 
259  // if not purely elastic or the customised stuff failed, do some plastic return
260  if (!found_solution)
262  _stress[_qp],
263  _intnl_old[_qp],
264  _intnl[_qp],
265  rot(_plastic_strain_old[_qp]),
266  _plastic_strain[_qp],
269  _yf[_qp],
270  number_iterations,
271  linesearch_needed,
272  ld_encountered,
273  constraints_added,
274  _Jacobian_mult[_qp]);
275 
276  if (_cosserat)
277  {
278  (*_couple_stress)[_qp] = (*_elastic_flexural_rigidity_tensor)[_qp] * _my_curvature;
279  (*_Jacobian_mult_couple)[_qp] = _my_flexural_rigidity_tensor;
280  }
281 
282  postReturnMap(); // rotate back from new frame if necessary
283 
284  _iter[_qp] = 1.0 * number_iterations;
285  _linesearch_needed[_qp] = linesearch_needed;
286  _ld_encountered[_qp] = ld_encountered;
287  _constraints_added[_qp] = constraints_added;
288 
289  // Update measures of strain
291  (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
292 
293  // Rotate the tensors to the current configuration
295  {
296  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
297  _elastic_strain[_qp] =
298  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
299  _plastic_strain[_qp] =
300  _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
301  }
302 }
303 
304 RankTwoTensor
305 ComputeMultiPlasticityStress::rot(const RankTwoTensor & tens)
306 {
307  if (!_n_supplied)
308  return tens;
309  return tens.rotated(_rot);
310 }
311 
312 void
314 {
315  if (_n_supplied)
316  {
317  // this is a rotation matrix that will rotate _n to the "z" axis
318  _rot = RotationMatrix::rotVecToZ(_n[_qp]);
319 
320  // rotate the tensors to this frame
321  _my_elasticity_tensor.rotate(_rot);
322  _my_strain_increment.rotate(_rot);
323  if (_cosserat)
324  {
326  _my_curvature.rotate(_rot);
327  }
328  }
329 }
330 
331 void
333 {
334  if (_n_supplied)
335  {
336  // this is a rotation matrix that will rotate "z" axis back to _n
337  _rot = _rot.transpose();
338 
339  // rotate the tensors back to original frame where _n is correctly oriented
340  _my_elasticity_tensor.rotate(_rot);
341  _Jacobian_mult[_qp].rotate(_rot);
342  _my_strain_increment.rotate(_rot);
343  _stress[_qp].rotate(_rot);
344  _plastic_strain[_qp].rotate(_rot);
345  if (_cosserat)
346  {
348  (*_Jacobian_mult_couple)[_qp].rotate(_rot);
349  _my_curvature.rotate(_rot);
350  (*_couple_stress)[_qp].rotate(_rot);
351  }
352 
353  // Rotate n by _rotation_increment
354  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
355  {
356  _n[_qp](i) = 0;
357  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
358  _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
359  }
360  }
361 }
362 
363 bool
364 ComputeMultiPlasticityStress::quickStep(const RankTwoTensor & stress_old,
365  RankTwoTensor & stress,
366  const std::vector<Real> & intnl_old,
367  std::vector<Real> & intnl,
368  std::vector<Real> & pm,
369  std::vector<Real> & cumulative_pm,
370  const RankTwoTensor & plastic_strain_old,
371  RankTwoTensor & plastic_strain,
372  const RankFourTensor & E_ijkl,
373  const RankTwoTensor & strain_increment,
374  std::vector<Real> & yf,
375  unsigned int & iterations,
376  RankFourTensor & consistent_tangent_operator,
377  const quickStep_called_from_t called_from,
378  bool final_step)
379 {
380  iterations = 0;
381 
382  unsigned num_plastic_returns;
383  RankTwoTensor delta_dp;
384 
385  // the following does the customized returnMap algorithm
386  // for all the plastic models.
387  unsigned custom_model = 0;
388  bool successful_return = returnMapAll(stress_old + E_ijkl * strain_increment,
389  intnl_old,
390  E_ijkl,
391  _epp_tol,
392  stress,
393  intnl,
394  pm,
395  cumulative_pm,
396  delta_dp,
397  yf,
398  num_plastic_returns,
399  custom_model);
400 
401  // the following updates the plastic_strain, when necessary
402  // and calculates the consistent_tangent_operator, when necessary
403  if (num_plastic_returns == 0)
404  {
405  // if successful_return = true, then a purely elastic step
406  // if successful_return = false, then >=1 plastic model is in
407  // inadmissible zone and failed to return using its customized
408  // returnMap function.
409  // In either case:
410  plastic_strain = plastic_strain_old;
411  if (successful_return && final_step)
412  {
413  if (called_from == computeQpStress_function)
414  consistent_tangent_operator = E_ijkl;
415  else // cannot necessarily use E_ijkl since different plastic models may have been active
416  // during other substeps
417  consistent_tangent_operator =
418  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
419  }
420  return successful_return;
421  }
422  else if (num_plastic_returns == 1 && successful_return)
423  {
424  // one model has successfully completed its custom returnMap algorithm
425  // and the other models have signalled they are elastic at
426  // the trial stress
427  plastic_strain = plastic_strain_old + delta_dp;
428  if (final_step)
429  {
430  if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
431  {
433  consistent_tangent_operator = E_ijkl;
434  else
435  {
436  std::vector<Real> custom_model_pm;
437  for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
438  custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
439  consistent_tangent_operator =
440  _f[custom_model]->consistentTangentOperator(stress_old + E_ijkl * strain_increment,
441  intnl_old[custom_model],
442  stress,
443  intnl[custom_model],
444  E_ijkl,
445  custom_model_pm);
446  }
447  }
448  else // cannot necessarily use the custom consistentTangentOperator since different plastic
449  // models may have been active during other substeps or the custom model says not to use
450  // its custom CTO algorithm
451  consistent_tangent_operator =
452  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
453  }
454  return true;
455  }
456  else // presumably returnMapAll is incorrectly coded!
457  mooseError("ComputeMultiPlasticityStress::quickStep should not get here!");
458 }
459 
460 bool
461 ComputeMultiPlasticityStress::plasticStep(const RankTwoTensor & stress_old,
462  RankTwoTensor & stress,
463  const std::vector<Real> & intnl_old,
464  std::vector<Real> & intnl,
465  const RankTwoTensor & plastic_strain_old,
466  RankTwoTensor & plastic_strain,
467  const RankFourTensor & E_ijkl,
468  const RankTwoTensor & strain_increment,
469  std::vector<Real> & yf,
470  unsigned int & iterations,
471  bool & linesearch_needed,
472  bool & ld_encountered,
473  bool & constraints_added,
474  RankFourTensor & consistent_tangent_operator)
475 {
491  bool return_successful = false;
492 
493  // total number of Newton-Raphson iterations used
494  unsigned int iter = 0;
495 
496  Real step_size = 1.0;
497  Real time_simulated = 0.0;
498 
499  // InitialStress Deprecation: remove the following 2 lines
500  if (_t_step >= 2)
501  _step_one = false;
502 
503  // the "good" variables hold the latest admissible stress
504  // and internal parameters.
505  RankTwoTensor stress_good = stress_old;
506  RankTwoTensor plastic_strain_good = plastic_strain_old;
507  std::vector<Real> intnl_good(_num_models);
508  for (unsigned model = 0; model < _num_models; ++model)
509  intnl_good[model] = intnl_old[model];
510  std::vector<Real> yf_good(_num_surfaces);
511 
512  // Following is necessary because I want strain_increment to be "const"
513  // but I also want to be able to subdivide an initial_stress
514  RankTwoTensor this_strain_increment = strain_increment;
515  // InitialStress Deprecation: remove following 6 lines
516  if (isParamValid("initial_stress") && _step_one)
517  {
518  RankFourTensor E_inv = E_ijkl.invSymm();
519  this_strain_increment += E_inv * stress_old;
520  stress_good = RankTwoTensor();
521  }
522 
523  RankTwoTensor dep = step_size * this_strain_increment;
524 
525  _cumulative_pm.assign(_num_surfaces, 0);
526 
527  unsigned int num_consecutive_successes = 0;
528  while (time_simulated < 1.0 && step_size >= _min_stepsize)
529  {
530  iter = 0;
531  return_successful = returnMap(stress_good,
532  stress,
533  intnl_good,
534  intnl,
535  plastic_strain_good,
536  plastic_strain,
537  E_ijkl,
538  dep,
539  yf,
540  iter,
541  step_size <= _max_stepsize_for_dumb,
542  linesearch_needed,
543  ld_encountered,
544  constraints_added,
545  time_simulated + step_size >= 1,
546  consistent_tangent_operator,
548  iterations += iter;
549 
550  if (return_successful)
551  {
552  num_consecutive_successes += 1;
553  time_simulated += step_size;
554 
555  if (time_simulated < 1.0) // this condition is just for optimization: if time_simulated=1 then
556  // the "good" quantities are no longer needed
557  {
558  stress_good = stress;
559  plastic_strain_good = plastic_strain;
560  for (unsigned model = 0; model < _num_models; ++model)
561  intnl_good[model] = intnl[model];
562  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
563  yf_good[surface] = yf[surface];
564  if (num_consecutive_successes >= 2)
565  step_size *= 1.2;
566  }
567  step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
568  }
569  else
570  {
571  step_size *= 0.5;
572  num_consecutive_successes = 0;
573  stress = stress_good;
574  plastic_strain = plastic_strain_good;
575  for (unsigned model = 0; model < _num_models; ++model)
576  intnl[model] = intnl_good[model];
577  yf.resize(_num_surfaces); // might have excited with junk
578  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
579  yf[surface] = yf_good[surface];
580  dep = step_size * this_strain_increment;
581  }
582  }
583 
584  if (!return_successful)
585  {
586  if (_ignore_failures)
587  {
588  stress = stress_good;
589  plastic_strain = plastic_strain_good;
590  for (unsigned model = 0; model < _num_models; ++model)
591  intnl[model] = intnl_good[model];
592  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
593  yf[surface] = yf_good[surface];
594  }
595  else
596  {
597  Moose::out << "After reducing the stepsize to " << step_size
598  << " with original strain increment with L2norm " << this_strain_increment.L2norm()
599  << " the returnMap algorithm failed\n";
600 
601  _fspb_debug_stress = stress_good + E_ijkl * dep;
602  _fspb_debug_pm.assign(
604  1); // this is chosen arbitrarily - please change if a more suitable value occurs to you!
605  _fspb_debug_intnl.resize(_num_models);
606  for (unsigned model = 0; model < _num_models; ++model)
607  _fspb_debug_intnl[model] = intnl_good[model];
611  mooseError("Exiting\n");
612  }
613  }
614 
615  return return_successful;
616 }
617 
618 bool
619 ComputeMultiPlasticityStress::returnMap(const RankTwoTensor & stress_old,
620  RankTwoTensor & stress,
621  const std::vector<Real> & intnl_old,
622  std::vector<Real> & intnl,
623  const RankTwoTensor & plastic_strain_old,
624  RankTwoTensor & plastic_strain,
625  const RankFourTensor & E_ijkl,
626  const RankTwoTensor & strain_increment,
627  std::vector<Real> & f,
628  unsigned int & iter,
629  bool can_revert_to_dumb,
630  bool & linesearch_needed,
631  bool & ld_encountered,
632  bool & constraints_added,
633  bool final_step,
634  RankFourTensor & consistent_tangent_operator,
635  std::vector<Real> & cumulative_pm)
636 {
637 
638  // The "consistency parameters" (plastic multipliers)
639  // Change in plastic strain in this timestep = pm*flowPotential
640  // Each pm must be non-negative
641  std::vector<Real> pm;
642  pm.assign(_num_surfaces, 0.0);
643 
644  bool successful_return = quickStep(stress_old,
645  stress,
646  intnl_old,
647  intnl,
648  pm,
649  cumulative_pm,
650  plastic_strain_old,
651  plastic_strain,
652  E_ijkl,
653  strain_increment,
654  f,
655  iter,
656  consistent_tangent_operator,
658  final_step);
659 
660  if (successful_return)
661  return successful_return;
662 
663  // Here we know that the trial stress and intnl_old
664  // is inadmissible, and we have to return from those values
665  // value to the yield surface. There are three
666  // types of constraints we have to satisfy, listed
667  // below, and calculated in calculateConstraints(...)
668  // f<=0, epp=0, ic=0 (up to tolerances), and these are
669  // guaranteed to hold if nr_res2<=0.5
670  //
671  // Kuhn-Tucker conditions must also be satisfied
672  // These are:
673  // pm>=0, which may not hold upon exit of the NR loops
674  // due to _deactivation_scheme!=optimized;
675  // pm*f=0 (up to tolerances), which may not hold upon exit
676  // of the NR loops if a constraint got deactivated
677  // due to linear dependence, and then f<0, and its pm>0
678 
679  // Plastic strain constraint, L2 norm must be zero (up to a tolerance)
680  RankTwoTensor epp;
681 
682  // Yield function constraint passed to this function as
683  // std::vector<Real> & f
684  // Each yield function must be <= 0 (up to tolerance)
685  // Note that only the constraints that are active will be
686  // contained in f until the final few lines of returnMap
687 
688  // Internal constraint(s), must be zero (up to a tolerance)
689  // Note that only the constraints that are active will be
690  // contained in ic.
691  std::vector<Real> ic;
692 
693  // Record the stress before Newton-Raphson in case of failure-and-restart
694  RankTwoTensor initial_stress = stress;
695 
696  iter = 0;
697 
698  // Initialize the set of active constraints
699  // At this stage, the active constraints are
700  // those that exceed their _f_tol
701  // active constraints.
702  std::vector<bool> act;
703  buildActiveConstraints(f, stress, intnl, E_ijkl, act);
704 
705  // Inverse of E_ijkl (assuming symmetric)
706  RankFourTensor E_inv = E_ijkl.invSymm();
707 
708  // convenience variable that holds the change in plastic strain incurred during the return
709  // delta_dp = plastic_strain - plastic_strain_old
710  // delta_dp = E^{-1}*(initial_stress - stress), where initial_stress = E*(strain -
711  // plastic_strain_old)
712  RankTwoTensor delta_dp = RankTwoTensor();
713 
714  // whether single step was successful (whether line search was successful, and whether turning off
715  // constraints was successful)
716  bool single_step_success = true;
717 
718  // deactivation scheme
720 
721  // For complicated deactivation schemes we have to record the initial active set
722  std::vector<bool> initial_act;
723  initial_act.resize(_num_surfaces);
727  {
728  // if "optimized" fails we can change the deactivation scheme to "safe", etc
729  deact_scheme = optimized;
730  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
731  initial_act[surface] = act[surface];
732  }
733 
735  deact_scheme = safe;
736 
737  // For "dumb" deactivation, the active set takes all combinations until a solution is found
738  int dumb_iteration = 0;
739  std::vector<unsigned int> dumb_order;
740 
741  if (_deactivation_scheme == dumb ||
742  (_deactivation_scheme == optimized_to_safe_to_dumb && can_revert_to_dumb) ||
743  (_deactivation_scheme == safe_to_dumb && can_revert_to_dumb) ||
744  (_deactivation_scheme == optimized_to_dumb && can_revert_to_dumb))
745  buildDumbOrder(stress, intnl, dumb_order);
746 
747  if (_deactivation_scheme == dumb)
748  {
749  incrementDumb(dumb_iteration, dumb_order, act);
750  yieldFunction(stress, intnl, act, f);
751  }
752 
753  // To avoid any re-trials of "act" combinations that
754  // we've already tried and rejected, i record the
755  // combinations in actives_tried
756  std::set<unsigned int> actives_tried;
757  actives_tried.insert(activeCombinationNumber(act));
758 
759  // The residual-squared that the line-search will reduce
760  // Later it will get contributions from epp and ic, but
761  // at present these are zero
762  Real nr_res2 = 0;
763  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
764  if (act[surface])
765  nr_res2 += 0.5 * Utility::pow<2>(f[surface] / _f[modelNumber(surface)]->_f_tol);
766 
767  successful_return = false;
768 
769  bool still_finding_solution = true;
770  while (still_finding_solution)
771  {
772  single_step_success = true;
773  unsigned int local_iter = 0;
774 
775  // The Newton-Raphson loops
776  while (nr_res2 > 0.5 && local_iter++ < _max_iter && single_step_success)
777  single_step_success = singleStep(nr_res2,
778  stress,
779  intnl_old,
780  intnl,
781  pm,
782  delta_dp,
783  E_inv,
784  f,
785  epp,
786  ic,
787  act,
788  deact_scheme,
789  linesearch_needed,
790  ld_encountered);
791 
792  bool nr_good = (nr_res2 <= 0.5 && local_iter <= _max_iter && single_step_success);
793 
794  iter += local_iter;
795 
796  // 'act' might have changed due to using deact_scheme = optimized, so
797  actives_tried.insert(activeCombinationNumber(act));
798 
799  if (!nr_good)
800  {
801  // failure of NR routine.
802  // We might be able to change the deactivation_scheme and
803  // then re-try the NR starting from the initial values
804  // Or, if deact_scheme == "dumb", we just increarse the
805  // dumb_iteration number and re-try
806  bool change_scheme = false;
807  bool increment_dumb = false;
808  change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
809  if (!change_scheme && deact_scheme == dumb)
810  increment_dumb = canIncrementDumb(dumb_iteration);
811 
812  still_finding_solution = (change_scheme || increment_dumb);
813 
814  if (change_scheme)
815  changeScheme(initial_act,
816  can_revert_to_dumb,
817  initial_stress,
818  intnl_old,
819  deact_scheme,
820  act,
821  dumb_iteration,
822  dumb_order);
823 
824  if (increment_dumb)
825  incrementDumb(dumb_iteration, dumb_order, act);
826 
827  if (!still_finding_solution)
828  {
829  // we cannot change the scheme, or have run out of "dumb" options
830  successful_return = false;
831  break;
832  }
833  }
834 
835  bool kt_good = false;
836  if (nr_good)
837  {
838  // check Kuhn-Tucker
839  kt_good = checkKuhnTucker(f, pm, act);
840  if (!kt_good)
841  {
842  if (deact_scheme != dumb)
843  {
844  applyKuhnTucker(f, pm, act);
845 
846  // true if we haven't tried this active set before
847  still_finding_solution =
848  (actives_tried.find(activeCombinationNumber(act)) == actives_tried.end());
849  if (!still_finding_solution)
850  {
851  // must have tried turning off the constraints already.
852  // so try changing the scheme
853  if (canChangeScheme(deact_scheme, can_revert_to_dumb))
854  {
855  still_finding_solution = true;
856  changeScheme(initial_act,
857  can_revert_to_dumb,
858  initial_stress,
859  intnl_old,
860  deact_scheme,
861  act,
862  dumb_iteration,
863  dumb_order);
864  }
865  }
866  }
867  else
868  {
869  bool increment_dumb = false;
870  increment_dumb = canIncrementDumb(dumb_iteration);
871  still_finding_solution = increment_dumb;
872  if (increment_dumb)
873  incrementDumb(dumb_iteration, dumb_order, act);
874  }
875 
876  if (!still_finding_solution)
877  {
878  // have tried turning off the constraints already,
879  // or have run out of "dumb" options
880  successful_return = false;
881  break;
882  }
883  }
884  }
885 
886  bool admissible = false;
887  if (nr_good && kt_good)
888  {
889  // check admissible
890  std::vector<Real> all_f;
891  if (_num_surfaces == 1)
892  admissible = true; // for a single surface if NR has exited successfully then (stress,
893  // intnl) must be admissible
894  else
895  admissible = checkAdmissible(stress, intnl, all_f);
896 
897  if (!admissible)
898  {
899  // Not admissible.
900  // We can try adding constraints back in
901  // We can try changing the deactivation scheme
902  // Or, if deact_scheme == dumb, just increase dumb_iteration
903  bool add_constraints = canAddConstraints(act, all_f);
904  if (add_constraints)
905  {
906  constraints_added = true;
907  std::vector<bool> act_plus(_num_surfaces,
908  false); // "act" with the positive constraints added in
909  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
910  if (act[surface] ||
911  (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol)))
912  act_plus[surface] = true;
913  if (actives_tried.find(activeCombinationNumber(act_plus)) == actives_tried.end())
914  {
915  // haven't tried this combination of actives yet
916  constraints_added = true;
917  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
918  act[surface] = act_plus[surface];
919  }
920  else
921  add_constraints = false; // haven't managed to add a new combination
922  }
923 
924  bool change_scheme = false;
925  bool increment_dumb = false;
926 
927  if (!add_constraints)
928  change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
929 
930  if (!add_constraints && !change_scheme && deact_scheme == dumb)
931  increment_dumb = canIncrementDumb(dumb_iteration);
932 
933  still_finding_solution = (add_constraints || change_scheme || increment_dumb);
934 
935  if (change_scheme)
936  changeScheme(initial_act,
937  can_revert_to_dumb,
938  initial_stress,
939  intnl_old,
940  deact_scheme,
941  act,
942  dumb_iteration,
943  dumb_order);
944 
945  if (increment_dumb)
946  incrementDumb(dumb_iteration, dumb_order, act);
947 
948  if (!still_finding_solution)
949  {
950  // we cannot change the scheme, or have run out of "dumb" options
951  successful_return = false;
952  break;
953  }
954  }
955  }
956 
957  successful_return = (nr_good && admissible && kt_good);
958  if (successful_return)
959  break;
960 
961  if (still_finding_solution)
962  {
963  stress = initial_stress;
964  delta_dp = RankTwoTensor(); // back to zero change in plastic strain
965  for (unsigned model = 0; model < _num_models; ++model)
966  intnl[model] = intnl_old[model]; // back to old internal params
967  pm.assign(_num_surfaces, 0.0); // back to zero plastic multipliers
968 
969  unsigned num_active = numberActive(act);
970  if (num_active == 0)
971  {
972  successful_return = false;
973  break; // failure
974  }
975 
976  actives_tried.insert(activeCombinationNumber(act));
977 
978  // Since "act" set has changed, either by changing deact_scheme, or by KT failing, so need to
979  // re-calculate nr_res2
980  yieldFunction(stress, intnl, act, f);
981 
982  nr_res2 = 0;
983  unsigned ind = 0;
984  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
985  if (act[surface])
986  {
987  if (f[ind] > _f[modelNumber(surface)]->_f_tol)
988  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
989  ind++;
990  }
991  }
992  }
993 
994  // returned, with either success or failure
995  if (successful_return)
996  {
997  plastic_strain += delta_dp;
998 
999  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1000  cumulative_pm[surface] += pm[surface];
1001 
1002  if (final_step)
1003  consistent_tangent_operator =
1004  consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
1005 
1006  if (f.size() != _num_surfaces)
1007  {
1008  // at this stage f.size() = num_active, but we need to return with all the yield functions
1009  // evaluated, so:
1010  act.assign(_num_surfaces, true);
1011  yieldFunction(stress, intnl, act, f);
1012  }
1013  }
1014 
1015  return successful_return;
1016 }
1017 
1018 bool
1020  const std::vector<Real> & all_f)
1021 {
1022  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1023  if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1024  return true;
1025  return false;
1026 }
1027 
1028 bool
1030  bool can_revert_to_dumb)
1031 {
1032  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
1033  return true;
1034 
1035  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
1036  return true;
1037 
1038  if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
1039  can_revert_to_dumb)
1040  return true;
1041 
1042  if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
1043  can_revert_to_dumb)
1044  return true;
1045 
1046  if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1047  can_revert_to_dumb)
1048  return true;
1049 
1050  return false;
1051 }
1052 
1053 void
1054 ComputeMultiPlasticityStress::changeScheme(const std::vector<bool> & initial_act,
1055  bool can_revert_to_dumb,
1056  const RankTwoTensor & initial_stress,
1057  const std::vector<Real> & intnl_old,
1058  DeactivationSchemeEnum & current_deactivation_scheme,
1059  std::vector<bool> & act,
1060  int & dumb_iteration,
1061  std::vector<unsigned int> & dumb_order)
1062 {
1063  if (current_deactivation_scheme == optimized &&
1066  {
1067  current_deactivation_scheme = safe;
1068  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1069  act[surface] = initial_act[surface];
1070  }
1071  else if ((current_deactivation_scheme == safe &&
1074  can_revert_to_dumb) ||
1075  (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1076  can_revert_to_dumb))
1077  {
1078  current_deactivation_scheme = dumb;
1079  dumb_iteration = 0;
1080  buildDumbOrder(initial_stress, intnl_old, dumb_order);
1081  incrementDumb(dumb_iteration, dumb_order, act);
1082  }
1083 }
1084 
1085 bool
1087  RankTwoTensor & stress,
1088  const std::vector<Real> & intnl_old,
1089  std::vector<Real> & intnl,
1090  std::vector<Real> & pm,
1091  RankTwoTensor & delta_dp,
1092  const RankFourTensor & E_inv,
1093  std::vector<Real> & f,
1094  RankTwoTensor & epp,
1095  std::vector<Real> & ic,
1096  std::vector<bool> & active,
1097  DeactivationSchemeEnum deactivation_scheme,
1098  bool & linesearch_needed,
1099  bool & ld_encountered)
1100 {
1101  bool successful_step; // return value
1102 
1103  Real nr_res2_before_step = nr_res2;
1104  RankTwoTensor stress_before_step;
1105  std::vector<Real> intnl_before_step;
1106  std::vector<Real> pm_before_step;
1107  RankTwoTensor delta_dp_before_step;
1108 
1109  if (deactivation_scheme == optimized)
1110  {
1111  // we potentially use the "before_step" quantities, so record them here
1112  stress_before_step = stress;
1113  intnl_before_step.resize(_num_models);
1114  for (unsigned model = 0; model < _num_models; ++model)
1115  intnl_before_step[model] = intnl[model];
1116  pm_before_step.resize(_num_surfaces);
1117  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1118  pm_before_step[surface] = pm[surface];
1119  delta_dp_before_step = delta_dp;
1120  }
1121 
1122  // During the Newton-Raphson procedure, we'll be
1123  // changing the following parameters in order to
1124  // (attempt to) satisfy the constraints.
1125  RankTwoTensor dstress; // change in stress
1126  std::vector<Real> dpm; // change in plasticity multipliers ("consistency parameters"). For ALL
1127  // contraints (active and deactive)
1128  std::vector<Real>
1129  dintnl; // change in internal parameters. For ALL internal params (active and deactive)
1130 
1131  // The constraints that have been deactivated for this NR step
1132  // due to the flow directions being linearly dependent
1133  std::vector<bool> deact_ld;
1134  deact_ld.assign(_num_surfaces, false);
1135 
1136  /* After NR and linesearch, if _deactivation_scheme == "optimized", the
1137  * active plasticity multipliers are checked for non-negativity. If some
1138  * are negative then they are deactivated forever, and the NR step is
1139  * re-done starting from the _before_step quantities recorded above
1140  */
1141  bool constraints_changing = true;
1142  bool reinstated_actives;
1143  while (constraints_changing)
1144  {
1145  // calculate dstress, dpm and dintnl for one full Newton-Raphson step
1146  nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
1147 
1148  for (unsigned surface = 0; surface < deact_ld.size(); ++surface)
1149  if (deact_ld[surface])
1150  {
1151  ld_encountered = true;
1152  break;
1153  }
1154 
1155  // perform a line search
1156  // The line-search will exit with updated values
1157  successful_step = lineSearch(nr_res2,
1158  stress,
1159  intnl_old,
1160  intnl,
1161  pm,
1162  E_inv,
1163  delta_dp,
1164  dstress,
1165  dpm,
1166  dintnl,
1167  f,
1168  epp,
1169  ic,
1170  active,
1171  deact_ld,
1172  linesearch_needed);
1173 
1174  if (!successful_step)
1175  // completely bomb out
1176  return successful_step;
1177 
1178  // See if any active constraints need to be removed, and the step re-done
1179  constraints_changing = false;
1180  if (deactivation_scheme == optimized)
1181  {
1182  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1183  if (active[surface] && pm[surface] < 0.0)
1184  constraints_changing = true;
1185  }
1186 
1187  if (constraints_changing)
1188  {
1189  stress = stress_before_step;
1190  delta_dp = delta_dp_before_step;
1191  nr_res2 = nr_res2_before_step;
1192  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1193  {
1194  if (active[surface] && pm[surface] < 0.0)
1195  {
1196  // turn off the constraint forever
1197  active[surface] = false;
1198  pm_before_step[surface] = 0.0;
1199  intnl_before_step[modelNumber(surface)] =
1200  intnl_old[modelNumber(surface)]; // don't want to muck-up hardening!
1201  }
1202  intnl[modelNumber(surface)] = intnl_before_step[modelNumber(surface)];
1203  pm[surface] = pm_before_step[surface];
1204  }
1205  if (numberActive(active) == 0)
1206  {
1207  // completely bomb out
1208  successful_step = false;
1209  return successful_step;
1210  }
1211  }
1212 
1213  // reinstate any active values that have been turned off due to linear-dependence
1214  reinstated_actives = reinstateLinearDependentConstraints(deact_ld);
1215  } // ends "constraints_changing" loop
1216 
1217  // if active constraints were reinstated then nr_res2 needs to be re-calculated so it is correct
1218  // upson returning
1219  if (reinstated_actives)
1220  {
1221  bool completely_converged = true;
1222  if (successful_step && nr_res2 < 0.5)
1223  {
1224  // Here we have converged to the correct solution if
1225  // all the yield functions are < 0. Excellent!
1226  //
1227  // This is quite tricky - perhaps i can refactor to make it more obvious.
1228  //
1229  // Because actives are now reinstated, the residual2
1230  // calculation below will give nr_res2 > 0.5, because it won't
1231  // realise that we only had to set the active-but-not-deactivated f = 0,
1232  // and not the entire active set. If we pass that nr_res2 back from
1233  // this function then the calling function will not realise we've converged!
1234  // Therefore, check for this case
1235  unsigned ind = 0;
1236  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1237  if (active[surface])
1238  if (f[ind++] > _f[modelNumber(surface)]->_f_tol)
1239  completely_converged = false;
1240  }
1241  else
1242  completely_converged = false;
1243 
1244  if (!completely_converged)
1245  nr_res2 = residual2(pm, f, epp, ic, active, deact_ld);
1246  }
1247 
1248  return successful_step;
1249 }
1250 
1251 bool
1253  std::vector<bool> & deactivated_due_to_ld)
1254 {
1255  bool reinstated_actives = false;
1256  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1257  if (deactivated_due_to_ld[surface])
1258  reinstated_actives = true;
1259 
1260  deactivated_due_to_ld.assign(_num_surfaces, false);
1261  return reinstated_actives;
1262 }
1263 
1264 unsigned
1265 ComputeMultiPlasticityStress::numberActive(const std::vector<bool> & active)
1266 {
1267  unsigned num_active = 0;
1268  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1269  if (active[surface])
1270  num_active++;
1271  return num_active;
1272 }
1273 
1274 bool
1276  const std::vector<Real> & intnl,
1277  std::vector<Real> & all_f)
1278 {
1279  std::vector<bool> act;
1280  act.assign(_num_surfaces, true);
1281 
1282  yieldFunction(stress, intnl, act, all_f);
1283 
1284  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1285  if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1286  return false;
1287 
1288  return true;
1289 }
1290 
1291 bool
1293  const std::vector<Real> & pm,
1294  const std::vector<bool> & active)
1295 {
1296  unsigned ind = 0;
1297  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1298  {
1299  if (active[surface])
1300  {
1301  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1302  if (pm[surface] != 0)
1303  return false;
1304  }
1305  else if (pm[surface] != 0)
1306  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1307  "coding!!");
1308  }
1309  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1310  if (pm[surface] < 0)
1311  return false;
1312 
1313  return true;
1314 }
1315 
1316 void
1318  const std::vector<Real> & pm,
1319  std::vector<bool> & active)
1320 {
1321  bool turned_off = false;
1322  unsigned ind = 0;
1323 
1324  // turn off all active surfaces that have f<0 and pm!=0
1325  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1326  {
1327  if (active[surface])
1328  {
1329  if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1330  if (pm[surface] != 0)
1331  {
1332  turned_off = true;
1333  active[surface] = false;
1334  }
1335  }
1336  else if (pm[surface] != 0)
1337  mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1338  "coding!!");
1339  }
1340 
1341  // if didn't turn off anything yet, turn off surface with minimum pm
1342  if (!turned_off)
1343  {
1344  int surface_to_turn_off = -1;
1345  Real min_pm = 0;
1346  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1347  if (pm[surface] < min_pm)
1348  {
1349  min_pm = pm[surface];
1350  surface_to_turn_off = surface;
1351  }
1352  if (surface_to_turn_off >= 0)
1353  active[surface_to_turn_off] = false;
1354  }
1355 }
1356 
1357 Real
1358 ComputeMultiPlasticityStress::residual2(const std::vector<Real> & pm,
1359  const std::vector<Real> & f,
1360  const RankTwoTensor & epp,
1361  const std::vector<Real> & ic,
1362  const std::vector<bool> & active,
1363  const std::vector<bool> & deactivated_due_to_ld)
1364 {
1365  Real nr_res2 = 0;
1366  unsigned ind = 0;
1367 
1368  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1369  if (active[surface])
1370  {
1371  if (!deactivated_due_to_ld[surface])
1372  {
1373  if (!(pm[surface] == 0 && f[ind] <= 0))
1374  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1375  }
1376  else if (deactivated_due_to_ld[surface] && f[ind] > 0)
1377  nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1378  ind++;
1379  }
1380 
1381  nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
1382 
1383  std::vector<bool> active_not_deact(_num_surfaces);
1384  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1385  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1386  ind = 0;
1387  for (unsigned model = 0; model < _num_models; ++model)
1388  if (anyActiveSurfaces(model, active_not_deact))
1389  nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
1390 
1391  return nr_res2;
1392 }
1393 
1394 bool
1396  RankTwoTensor & stress,
1397  const std::vector<Real> & intnl_old,
1398  std::vector<Real> & intnl,
1399  std::vector<Real> & pm,
1400  const RankFourTensor & E_inv,
1401  RankTwoTensor & delta_dp,
1402  const RankTwoTensor & dstress,
1403  const std::vector<Real> & dpm,
1404  const std::vector<Real> & dintnl,
1405  std::vector<Real> & f,
1406  RankTwoTensor & epp,
1407  std::vector<Real> & ic,
1408  const std::vector<bool> & active,
1409  const std::vector<bool> & deactivated_due_to_ld,
1410  bool & linesearch_needed)
1411 {
1412  // Line search algorithm straight out of "Numerical Recipes"
1413 
1414  bool success =
1415  true; // return value: will be false if linesearch couldn't reduce the residual-squared
1416 
1417  // Aim is to decrease residual2
1418 
1419  Real lam = 1.0; // the line-search parameter: 1.0 is a full Newton step
1420  Real lam_min =
1421  1E-10; // minimum value of lam allowed - perhaps this should be dynamically calculated?
1422  Real f0 = nr_res2; // initial value of residual2
1423  Real slope = -2 * nr_res2; // "Numerical Recipes" uses -b*A*x, in order to check for roundoff, but
1424  // i hope the nrStep would warn if there were problems.
1425  Real tmp_lam; // cached value of lam used in quadratic & cubic line search
1426  Real f2 = nr_res2; // cached value of f = residual2 used in the cubic in the line search
1427  Real lam2 = lam; // cached value of lam used in the cubic in the line search
1428 
1429  // pm during the line-search
1430  std::vector<Real> ls_pm;
1431  ls_pm.resize(pm.size());
1432 
1433  // delta_dp during the line-search
1434  RankTwoTensor ls_delta_dp;
1435 
1436  // internal parameter during the line-search
1437  std::vector<Real> ls_intnl;
1438  ls_intnl.resize(intnl.size());
1439 
1440  // stress during the line-search
1441  RankTwoTensor ls_stress;
1442 
1443  // flow directions (not used in line search, but calculateConstraints returns this parameter)
1444  std::vector<RankTwoTensor> r;
1445 
1446  while (true)
1447  {
1448  // update the variables using this line-search parameter
1449  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1450  ls_pm[alpha] = pm[alpha] + dpm[alpha] * lam;
1451  ls_delta_dp = delta_dp - E_inv * dstress * lam;
1452  for (unsigned a = 0; a < intnl.size(); ++a)
1453  ls_intnl[a] = intnl[a] + dintnl[a] * lam;
1454  ls_stress = stress + dstress * lam;
1455 
1456  // calculate the new active yield functions, epp and active internal constraints
1457  calculateConstraints(ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1458 
1459  // calculate the new residual-squared
1460  nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1461 
1462  if (nr_res2 < f0 + 1E-4 * lam * slope)
1463  break;
1464  else if (lam < lam_min)
1465  {
1466  success = false;
1467  // restore plastic multipliers, yield functions, etc to original values
1468  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1469  ls_pm[alpha] = pm[alpha];
1470  ls_delta_dp = delta_dp;
1471  for (unsigned a = 0; a < intnl.size(); ++a)
1472  ls_intnl[a] = intnl[a];
1473  ls_stress = stress;
1475  ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1476  nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1477  break;
1478  }
1479  else if (lam == 1.0)
1480  {
1481  // model as a quadratic
1482  tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
1483  }
1484  else
1485  {
1486  // model as a cubic
1487  Real rhs1 = nr_res2 - f0 - lam * slope;
1488  Real rhs2 = f2 - f0 - lam2 * slope;
1489  Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1490  Real b =
1491  (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1492  if (a == 0)
1493  tmp_lam = -slope / (2.0 * b);
1494  else
1495  {
1496  Real disc = Utility::pow<2>(b) - 3 * a * slope;
1497  if (disc < 0)
1498  tmp_lam = 0.5 * lam;
1499  else if (b <= 0)
1500  tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
1501  else
1502  tmp_lam = -slope / (b + std::sqrt(disc));
1503  }
1504  if (tmp_lam > 0.5 * lam)
1505  tmp_lam = 0.5 * lam;
1506  }
1507  lam2 = lam;
1508  f2 = nr_res2;
1509  lam = std::max(tmp_lam, 0.1 * lam);
1510  }
1511 
1512  if (lam < 1.0)
1513  linesearch_needed = true;
1514 
1515  // assign the quantities found in the line-search
1516  // back to the originals
1517  for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1518  pm[alpha] = ls_pm[alpha];
1519  delta_dp = ls_delta_dp;
1520  for (unsigned a = 0; a < intnl.size(); ++a)
1521  intnl[a] = ls_intnl[a];
1522  stress = ls_stress;
1523 
1524  return success;
1525 }
1526 
1527 void
1529  const std::vector<Real> & intnl,
1530  std::vector<unsigned int> & dumb_order)
1531 {
1532  if (dumb_order.size() != 0)
1533  return;
1534 
1535  std::vector<bool> act;
1536  act.assign(_num_surfaces, true);
1537 
1538  std::vector<Real> f;
1539  yieldFunction(stress, intnl, act, f);
1540  std::vector<RankTwoTensor> df_dstress;
1541  dyieldFunction_dstress(stress, intnl, act, df_dstress);
1542 
1543  typedef std::pair<Real, unsigned> pair_for_sorting;
1544  std::vector<pair_for_sorting> dist(_num_surfaces);
1545  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1546  {
1547  dist[surface].first = f[surface] / df_dstress[surface].L2norm();
1548  dist[surface].second = surface;
1549  }
1550  std::sort(dist.begin(), dist.end()); // sorted in ascending order of f/df_dstress
1551 
1552  dumb_order.resize(_num_surfaces);
1553  for (unsigned i = 0; i < _num_surfaces; ++i)
1554  dumb_order[i] = dist[_num_surfaces - 1 - i].second;
1555  // now dumb_order[0] is the surface with the greatest f/df_dstress
1556 }
1557 
1558 void
1560  const std::vector<unsigned int> & dumb_order,
1561  std::vector<bool> & act)
1562 {
1563  dumb_iteration += 1;
1564  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1565  act[dumb_order[surface]] =
1566  (dumb_iteration &
1567  (1 << surface)); // returns true if the surface_th bit of dumb_iteration == 1
1568 }
1569 
1570 bool
1572 {
1573  // (1 << _num_surfaces) = 2^_num_surfaces
1574  return ((dumb_iteration + 1) < (1 << _num_surfaces));
1575 }
1576 
1577 unsigned int
1579 {
1580  unsigned num = 0;
1581  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1582  if (act[surface])
1583  num += (1 << surface); // (1 << x) = 2^x
1584 
1585  return num;
1586 }
1587 
1588 RankFourTensor
1590  const std::vector<Real> & intnl,
1591  const RankFourTensor & E_ijkl,
1592  const std::vector<Real> & pm_this_step,
1593  const std::vector<Real> & cumulative_pm)
1594 {
1595 
1597  return E_ijkl;
1598 
1599  // Typically act_at_some_step = act, but it is possible
1600  // that when subdividing a strain increment, a surface
1601  // is only active for one sub-step
1602  std::vector<bool> act_at_some_step(_num_surfaces);
1603  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1604  act_at_some_step[surface] = (cumulative_pm[surface] > 0);
1605 
1606  // "act" might contain surfaces that are linearly dependent
1607  // with others. Only the plastic multipliers that are > 0
1608  // for this strain increment need to be varied to find
1609  // the consistent tangent operator
1610  std::vector<bool> act_vary(_num_surfaces);
1611  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1612  act_vary[surface] = (pm_this_step[surface] > 0);
1613 
1614  std::vector<RankTwoTensor> df_dstress;
1615  dyieldFunction_dstress(stress, intnl, act_vary, df_dstress);
1616  std::vector<Real> df_dintnl;
1617  dyieldFunction_dintnl(stress, intnl, act_vary, df_dintnl);
1618  std::vector<RankTwoTensor> r;
1619  flowPotential(stress, intnl, act_vary, r);
1620  std::vector<RankFourTensor> dr_dstress_at_some_step;
1621  dflowPotential_dstress(stress, intnl, act_at_some_step, dr_dstress_at_some_step);
1622  std::vector<RankTwoTensor> dr_dintnl_at_some_step;
1623  dflowPotential_dintnl(stress, intnl, act_at_some_step, dr_dintnl_at_some_step);
1624  std::vector<Real> h;
1625  hardPotential(stress, intnl, act_vary, h);
1626 
1627  unsigned ind1;
1628  unsigned ind2;
1629 
1630  // r_minus_stuff[alpha] = r[alpha] -
1631  // pm_cumulatve[gamma]*dr[gamma]_dintnl[a]_at_some_step*h[a][alpha], with alpha only being in
1632  // act_vary, but gamma being act_at_some_step
1633  std::vector<RankTwoTensor> r_minus_stuff;
1634  ind1 = 0;
1635  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1636  if (act_vary[surface1])
1637  {
1638  r_minus_stuff.push_back(r[ind1]);
1639  ind2 = 0;
1640  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1641  if (act_at_some_step[surface2])
1642  {
1643  if (modelNumber(surface1) == modelNumber(surface2))
1644  {
1645  r_minus_stuff.back() -=
1646  cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
1647  }
1648  ind2++;
1649  }
1650  ind1++;
1651  }
1652 
1653  unsigned int num_currently_active = 0;
1654  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1655  if (act_vary[surface])
1656  num_currently_active += 1;
1657 
1658  // zzz is a matrix in the form that can be easily
1659  // inverted by MatrixTools::inverse
1660  // Eg for num_currently_active = 3
1661  // (zzz[0] zzz[1] zzz[2])
1662  // (zzz[3] zzz[4] zzz[5])
1663  // (zzz[6] zzz[7] zzz[8])
1664  std::vector<PetscScalar> zzz;
1665  zzz.assign(num_currently_active * num_currently_active, 0.0);
1666 
1667  ind1 = 0;
1668  RankTwoTensor r2;
1669  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1670  if (act_vary[surface1])
1671  {
1672  ind2 = 0;
1673  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1674  if (act_vary[surface2])
1675  {
1676  r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
1677  zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
1678  if (modelNumber(surface1) == modelNumber(surface2))
1679  zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
1680  ind2++;
1681  }
1682  ind1++;
1683  }
1684 
1685  if (num_currently_active > 0)
1686  {
1687  // invert zzz, in place. if num_currently_active = 0 then zzz is not needed.
1688  try
1689  {
1690  MatrixTools::inverse(zzz, num_currently_active);
1691  }
1692  catch (const MooseException & e)
1693  {
1694  // in the very rare case of zzz being singular, just return the "elastic" tangent operator
1695  return E_ijkl;
1696  }
1697  }
1698 
1699  RankFourTensor strain_coeff = E_ijkl;
1700  ind1 = 0;
1701  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1702  if (act_vary[surface1])
1703  {
1704  RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1705  ind2 = 0;
1706  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1707  if (act_vary[surface2])
1708  {
1709  RankTwoTensor part2 = E_ijkl * df_dstress[ind2];
1710  for (unsigned i = 0; i < 3; i++)
1711  for (unsigned j = 0; j < 3; j++)
1712  for (unsigned k = 0; k < 3; k++)
1713  for (unsigned l = 0; l < 3; l++)
1714  strain_coeff(i, j, k, l) -=
1715  part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1716  ind2++;
1717  }
1718  ind1++;
1719  }
1720 
1722  return strain_coeff;
1723 
1724  RankFourTensor stress_coeff(RankFourTensor::initIdentitySymmetricFour);
1725 
1726  RankFourTensor part3;
1727  ind1 = 0;
1728  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1729  if (act_at_some_step[surface1])
1730  {
1731  part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
1732  ind1++;
1733  }
1734 
1735  stress_coeff += part3;
1736 
1737  part3 = part3.transposeMajor(); // this is because below i want df_dstress[ind2]*part3, and this
1738  // equals (part3.transposeMajor())*df_dstress[ind2]
1739 
1740  ind1 = 0;
1741  for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1742  if (act_vary[surface1])
1743  {
1744  RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1745  ind2 = 0;
1746  for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1747  if (act_vary[surface2])
1748  {
1749  RankTwoTensor part2 = part3 * df_dstress[ind2];
1750  for (unsigned i = 0; i < 3; i++)
1751  for (unsigned j = 0; j < 3; j++)
1752  for (unsigned k = 0; k < 3; k++)
1753  for (unsigned l = 0; l < 3; l++)
1754  stress_coeff(i, j, k, l) -=
1755  part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1756  ind2++;
1757  }
1758  ind1++;
1759  }
1760 
1761  // need to find the inverse of stress_coeff, but remember
1762  // stress_coeff does not have the symmetries commonly found
1763  // in tensor mechanics:
1764  // stress_coeff(i, j, k, l) = stress_coeff(j, i, k, l) = stress_coeff(i, j, l, k) !=
1765  // stress_coeff(k, l, i, j)
1766  // (note the final not-equals). We want s_inv, such that
1767  // s_inv(i, j, m, n)*stress_coeff(m, n, k, l) = (de_ik*de_jl + de_il*de_jk)/2
1768  // where de_ij = 1 if i=j, and 0 otherwise.
1769  RankFourTensor s_inv;
1770  try
1771  {
1772  s_inv = stress_coeff.invSymm();
1773  }
1774  catch (const MooseException & e)
1775  {
1776  return strain_coeff; // when stress_coeff is singular (perhaps for incompressible plasticity?)
1777  // return the "linear" tangent operator
1778  }
1779 
1780  return s_inv * strain_coeff;
1781 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
bool _perform_finite_strain_rotations
whether to perform the rotations necessary in finite-strain simulations
virtual bool reinstateLinearDependentConstraints(std::vector< bool > &deactivated_due_to_ld)
makes all deactivated_due_to_ld false, and if >0 of them were initially true, returns true ...
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
ComputeStressBase is the base class for stress tensors.
bool canAddConstraints(const std::vector< bool > &act, const std::vector< Real > &all_f)
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
quickStep_called_from_t
The functions from which quickStep can be called.
MooseEnum _fspb_debug
none - don&#39;t do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
virtual bool singleStep(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, RankTwoTensor &delta_dp, const RankFourTensor &E_inv, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, std::vector< bool > &active, DeactivationSchemeEnum deactivation_scheme, bool &linesearch_needed, bool &ld_encountered)
Performs a single Newton-Raphson + linesearch step Constraints are deactivated and the step is re-don...
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.
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment (coming from ComputeIncrementalSmallStrain, for example)
const MaterialProperty< RealVectorValue > & _n_old
old value of transverse direction
MaterialProperty< Real > & _constraints_added
Whether constraints were added in during the latest Newton-Raphson process (1 if true, 0 otherwise)
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
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.
bool canChangeScheme(DeactivationSchemeEnum current_deactivation_scheme, bool can_revert_to_dumb)
RankFourTensor _my_flexural_rigidity_tensor
Flexual rigidity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankTwoTensor > & _strain_increment
strain increment (coming from ComputeIncrementalSmallStrain, for example)
TangentOperatorEnum
The type of tangent operator to return. tangent operator = d(stress_rate)/d(strain_rate).
const MaterialProperty< RankTwoTensor > & _stress_old
Old value of stress.
unsigned int _max_iter
Maximum number of Newton-Raphson iterations allowed.
RankTwoTensor _my_strain_increment
Strain increment that can be rotated by this class, and split into multiple increments (ie...
unsigned int activeCombinationNumber(const std::vector< bool > &act)
virtual void applyKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
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)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
bool _ignore_failures
Even if the returnMap fails, return the best values found for stress and internal parameters...
RankFourTensor _my_elasticity_tensor
Elasticity tensor that can be rotated by this class (ie, its not const)
MaterialProperty< Real > & _iter
Number of Newton-Raphson iterations used in the return-map.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Old value of elastic strain.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
InputParameters validParams< ComputeMultiPlasticityStress >()
bool _n_supplied
User supplied the transverse direction vector.
MaterialProperty< Real > & _ld_encountered
Whether linear-dependence was encountered in the latest Newton-Raphson process (1 if true...
std::vector< const TensorMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.
void changeScheme(const std::vector< bool > &initial_act, bool can_revert_to_dumb, const RankTwoTensor &initial_stress, const std::vector< Real > &intnl_old, DeactivationSchemeEnum &current_deactivation_scheme, std::vector< bool > &act, int &dumb_iteration, std::vector< unsigned int > &dumb_order)
virtual bool checkKuhnTucker(const std::vector< Real > &f, const std::vector< Real > &pm, const std::vector< bool > &active)
Checks Kuhn-Tucker conditions, and alters "active" if appropriate.
virtual bool lineSearch(Real &nr_res2, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, const RankFourTensor &E_inv, RankTwoTensor &delta_dp, const RankTwoTensor &dstress, const std::vector< Real > &dpm, const std::vector< Real > &dintnl, std::vector< Real > &f, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, bool &linesearch_needed)
Performs a line search.
RankTwoTensor rot(const RankTwoTensor &tens)
MaterialProperty< std::vector< Real > > & _yf
yield functions
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...
void checkSolution(const RankFourTensor &E_inv)
Checks that Ax does equal b in the NR procedure.
void checkDerivatives()
Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.
MaterialProperty< RealVectorValue > & _n
current value of transverse direction
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; ...
std::vector< Real > _dummy_pm
dummy "consistency parameters" (plastic multipliers) used in quickStep when called from computeQpStre...
bool & _step_one
True if this is the first timestep (timestep < 2).
virtual void incrementDumb(int &dumb_iteration, const std::vector< unsigned int > &dumb_order, std::vector< bool > &act)
Increments "dumb_iteration" by 1, and sets "act" appropriately (act[alpha] = true iff alpha_th bit of...
MaterialProperty< Real > & _linesearch_needed
Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise) ...
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.
RealVectorValue _n_input
the supplied transverse direction vector
Real _max_stepsize_for_dumb
"dumb" deactivation will only be used if the stepsize falls below this quantity
const MaterialProperty< std::vector< Real > > & _intnl_old
old values of internal parameters
enum ComputeMultiPlasticityStress::TangentOperatorEnum _tangent_operator_type
virtual bool checkAdmissible(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &all_f)
Checks whether the yield functions are in the admissible region.
virtual bool returnMap(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &f, unsigned int &iter, bool can_revert_to_dumb, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, bool final_step, RankFourTensor &consistent_tangent_operator, std::vector< Real > &cumulative_pm)
Implements the return map.
virtual Real residual2(const std::vector< Real > &pm, const std::vector< Real > &f, const RankTwoTensor &epp, const std::vector< Real > &ic, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld)
The residual-squared.
InputParameters validParams< ComputeStressBase >()
ComputeMultiPlasticityStress(const InputParameters &parameters)
MaterialProperty< RankTwoTensor > & _plastic_strain
plastic strain
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.
Real _epp_tol
Tolerance on the plastic strain increment ("direction") constraint.
virtual void initQpStatefulProperties() override
RankFourTensor consistentTangentOperator(const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &pm_this_step, const std::vector< Real > &cumulative_pm)
Computes the consistent tangent operator (another name for the jacobian = d(stress_rate)/d(strain_rat...
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Real _min_stepsize
Minimum fraction of applied strain that may be applied during adaptive stepsizing.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in t...
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...
RankTwoTensor _my_curvature
Curvature that can be rotated by this class, and split into multiple increments (ie, its not const)
std::vector< Real > _cumulative_pm
the sum of the plastic multipliers over all the sub-steps.
MaterialProperty< RankTwoTensor > & _elastic_strain
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Old value of plastic strain.
virtual unsigned int numberActive(const std::vector< bool > &active)
counts the number of active constraints
RealTensorValue _rot
rotation matrix that takes _n to (0, 0, 1)
virtual bool quickStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, RankFourTensor &consistent_tangent_operator, const quickStep_called_from_t called_from, bool final_step)
Attempts to find an admissible (stress, intnl) by using the customized return-map algorithms defined ...
bool _cosserat
whether Cosserat mechanics should be used
void buildDumbOrder(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< unsigned int > &dumb_order)
Builds the order which "dumb" activation will take.
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...
enum ComputeMultiPlasticityStress::DeactivationSchemeEnum _deactivation_scheme
virtual bool plasticStep(const RankTwoTensor &stress_old, RankTwoTensor &stress, const std::vector< Real > &intnl_old, std::vector< Real > &intnl, const RankTwoTensor &plastic_strain_old, RankTwoTensor &plastic_strain, const RankFourTensor &E_ijkl, const RankTwoTensor &strain_increment, std::vector< Real > &yf, unsigned int &iterations, bool &linesearch_needed, bool &ld_encountered, bool &constraints_added, RankFourTensor &consistent_tangent_operator)
performs a plastic step
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.
virtual void 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 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.
InputParameters validParams< MultiPlasticityDebugger >()