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