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