Line data Source code
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 :
10 : #include "ComputeMultiPlasticityStress.h"
11 : #include "MultiPlasticityDebugger.h"
12 : #include "MatrixTools.h"
13 :
14 : #include "MooseException.h"
15 : #include "RotationMatrix.h" // for rotVecToZ
16 :
17 : #include "libmesh/utility.h"
18 :
19 : registerMooseObject("SolidMechanicsApp", ComputeMultiPlasticityStress);
20 :
21 : InputParameters
22 2655 : ComputeMultiPlasticityStress::validParams()
23 : {
24 2655 : InputParameters params = ComputeStressBase::validParams();
25 2655 : params += MultiPlasticityDebugger::validParams();
26 2655 : params.addClassDescription("Base class for multi-surface finite-strain plasticity");
27 7965 : params.addRangeCheckedParam<unsigned int>("max_NR_iterations",
28 5310 : 20,
29 : "max_NR_iterations>0",
30 : "Maximum number of Newton-Raphson iterations allowed");
31 5310 : 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 7965 : params.addRangeCheckedParam<Real>(
36 : "min_stepsize",
37 5310 : 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 7965 : params.addRangeCheckedParam<Real>("max_stepsize_for_dumb",
43 5310 : 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 5310 : "optimized");
52 5310 : 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 5310 : 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 5310 : params.addParam<bool>("ignore_failures",
71 5310 : 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 5310 : MooseEnum tangent_operator("elastic linear nonlinear", "nonlinear");
78 5310 : 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 5310 : params.addParam<bool>("perform_finite_strain_rotations",
87 5310 : 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 2655 : params.addClassDescription("Material for multi-surface finite-strain plasticity");
94 2655 : return params;
95 2655 : }
96 :
97 1987 : ComputeMultiPlasticityStress::ComputeMultiPlasticityStress(const InputParameters & parameters)
98 : : ComputeStressBase(parameters),
99 : MultiPlasticityDebugger(this),
100 1987 : _max_iter(getParam<unsigned int>("max_NR_iterations")),
101 3974 : _min_stepsize(getParam<Real>("min_stepsize")),
102 3974 : _max_stepsize_for_dumb(getParam<Real>("max_stepsize_for_dumb")),
103 3974 : _ignore_failures(getParam<bool>("ignore_failures")),
104 :
105 3974 : _tangent_operator_type((TangentOperatorEnum)(int)getParam<MooseEnum>("tangent_operator")),
106 :
107 3974 : _epp_tol(getParam<Real>("ep_plastic_tolerance")),
108 :
109 1987 : _dummy_pm(0),
110 :
111 1987 : _cumulative_pm(0),
112 :
113 3974 : _deactivation_scheme((DeactivationSchemeEnum)(int)getParam<MooseEnum>("deactivation_scheme")),
114 :
115 1987 : _n_supplied(parameters.isParamValid("transverse_direction")),
116 4209 : _n_input(_n_supplied ? getParam<RealVectorValue>("transverse_direction") : RealVectorValue()),
117 : _rot(RealTensorValue()),
118 :
119 3974 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
120 :
121 1987 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
122 1987 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
123 1987 : _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
124 3974 : _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
125 1987 : _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
126 3974 : _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
127 1987 : _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
128 1987 : _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
129 : // for visualisation i convert it to Real
130 1987 : _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
131 : // boolean, but for
132 : // visualisation i
133 : // convert it to Real
134 1987 : _ld_encountered(declareProperty<Real>(
135 : "plastic_linear_dependence_encountered")), // this is really a boolean, but for
136 : // visualisation i convert it to Real
137 1987 : _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
138 : // boolean, but for
139 : // visualisation i
140 : // convert it to Real
141 1987 : _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
142 3974 : _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
143 :
144 3974 : _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
145 3974 : _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
146 1987 : _rotation_increment(
147 1987 : getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
148 :
149 3974 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
150 3974 : _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 1987 : _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
155 1999 : hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
156 1999 : _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : nullptr),
157 1987 : _elastic_flexural_rigidity_tensor(
158 1999 : _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
159 : : nullptr),
160 1987 : _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : nullptr),
161 1999 : _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress")
162 : : nullptr),
163 1987 : _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
164 : : nullptr),
165 :
166 1987 : _my_elasticity_tensor(RankFourTensor()),
167 1987 : _my_strain_increment(RankTwoTensor()),
168 1987 : _my_flexural_rigidity_tensor(RankFourTensor()),
169 3974 : _my_curvature(RankTwoTensor())
170 : {
171 1987 : if (_epp_tol <= 0)
172 0 : mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
173 :
174 1987 : if (_n_supplied)
175 : {
176 : // normalise the inputted transverse_direction
177 235 : if (_n_input.norm() == 0)
178 1 : mooseError(
179 : "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
180 : else
181 234 : _n_input /= _n_input.norm();
182 : }
183 :
184 1986 : if (_num_surfaces == 1)
185 1182 : _deactivation_scheme = safe;
186 1986 : }
187 :
188 : void
189 49616 : ComputeMultiPlasticityStress::initQpStatefulProperties()
190 : {
191 49616 : ComputeStressBase::initQpStatefulProperties();
192 :
193 49616 : _plastic_strain[_qp].zero();
194 :
195 49616 : _intnl[_qp].assign(_num_models, 0);
196 :
197 49616 : _yf[_qp].assign(_num_surfaces, 0);
198 :
199 49616 : _dummy_pm.assign(_num_surfaces, 0);
200 :
201 49616 : _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
202 49616 : _linesearch_needed[_qp] = 0;
203 49616 : _ld_encountered[_qp] = 0;
204 49616 : _constraints_added[_qp] = 0;
205 :
206 49616 : _n[_qp] = _n_input;
207 :
208 49616 : if (_cosserat)
209 4800 : (*_couple_stress)[_qp].zero();
210 :
211 49616 : if (_fspb_debug == "jacobian")
212 : {
213 0 : checkDerivatives();
214 0 : mooseError("Finite-differencing completed. Exiting with no error");
215 : }
216 49616 : }
217 :
218 : void
219 2795336 : ComputeMultiPlasticityStress::computeQpStress()
220 : {
221 : // the following "_my" variables can get rotated by preReturnMap and postReturnMap
222 2795336 : _my_elasticity_tensor = _elasticity_tensor[_qp];
223 2795336 : _my_strain_increment = _strain_increment[_qp];
224 2795336 : if (_cosserat)
225 : {
226 19200 : _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
227 19200 : _my_curvature = (*_curvature)[_qp];
228 : }
229 :
230 2795336 : if (_fspb_debug == "jacobian_and_linear_system")
231 : {
232 : // cannot do this at initQpStatefulProperties level since E_ijkl is not defined
233 0 : checkJacobian(_elasticity_tensor[_qp].invSymm(), _intnl_old[_qp]);
234 0 : checkSolution(_elasticity_tensor[_qp].invSymm());
235 0 : mooseError("Finite-differencing completed. Exiting with no error");
236 : }
237 :
238 2795336 : preReturnMap(); // do rotations to new frame if necessary
239 :
240 : unsigned int number_iterations;
241 2795336 : bool linesearch_needed = false;
242 2795336 : bool ld_encountered = false;
243 2795336 : bool constraints_added = false;
244 :
245 2795336 : _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 2795336 : const bool found_solution = quickStep(rot(_stress_old[_qp]),
249 2795336 : _stress[_qp],
250 2795336 : _intnl_old[_qp],
251 2795336 : _intnl[_qp],
252 2795336 : _dummy_pm,
253 : _cumulative_pm,
254 5590672 : rot(_plastic_strain_old[_qp]),
255 2795336 : _plastic_strain[_qp],
256 2795336 : _my_elasticity_tensor,
257 2795336 : _my_strain_increment,
258 2795336 : _yf[_qp],
259 : number_iterations,
260 2795336 : _Jacobian_mult[_qp],
261 : computeQpStress_function,
262 : true);
263 :
264 : // if not purely elastic or the customised stuff failed, do some plastic return
265 2795336 : if (!found_solution)
266 108128 : plasticStep(rot(_stress_old[_qp]),
267 108128 : _stress[_qp],
268 108128 : _intnl_old[_qp],
269 108128 : _intnl[_qp],
270 108128 : rot(_plastic_strain_old[_qp]),
271 108128 : _plastic_strain[_qp],
272 : _my_elasticity_tensor,
273 : _my_strain_increment,
274 108128 : _yf[_qp],
275 : number_iterations,
276 : linesearch_needed,
277 : ld_encountered,
278 : constraints_added,
279 108128 : _Jacobian_mult[_qp]);
280 :
281 2795336 : if (_cosserat)
282 : {
283 19200 : (*_couple_stress)[_qp] = (*_elastic_flexural_rigidity_tensor)[_qp] * _my_curvature;
284 19200 : (*_Jacobian_mult_couple)[_qp] = _my_flexural_rigidity_tensor;
285 : }
286 :
287 2795336 : postReturnMap(); // rotate back from new frame if necessary
288 :
289 2795336 : _iter[_qp] = 1.0 * number_iterations;
290 2795336 : _linesearch_needed[_qp] = linesearch_needed;
291 2795336 : _ld_encountered[_qp] = ld_encountered;
292 2795336 : _constraints_added[_qp] = constraints_added;
293 :
294 : // Update measures of strain
295 2795336 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + _my_strain_increment -
296 2795336 : (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
297 :
298 : // Rotate the tensors to the current configuration
299 2795336 : if (_perform_finite_strain_rotations)
300 : {
301 2795336 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
302 2795336 : _elastic_strain[_qp] =
303 2795336 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
304 2795336 : _plastic_strain[_qp] =
305 2795336 : _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
306 : }
307 2795336 : }
308 :
309 : RankTwoTensor
310 5806928 : ComputeMultiPlasticityStress::rot(const RankTwoTensor & tens)
311 : {
312 5806928 : if (!_n_supplied)
313 5768544 : return tens;
314 38384 : return tens.rotated(_rot);
315 : }
316 :
317 : void
318 2795336 : ComputeMultiPlasticityStress::preReturnMap()
319 : {
320 2795336 : if (_n_supplied)
321 : {
322 : // this is a rotation matrix that will rotate _n to the "z" axis
323 10440 : _rot = RotationMatrix::rotVecToZ(_n[_qp]);
324 :
325 : // rotate the tensors to this frame
326 10440 : _my_elasticity_tensor.rotate(_rot);
327 10440 : _my_strain_increment.rotate(_rot);
328 10440 : if (_cosserat)
329 : {
330 0 : _my_flexural_rigidity_tensor.rotate(_rot);
331 0 : _my_curvature.rotate(_rot);
332 : }
333 : }
334 2795336 : }
335 :
336 : void
337 2795336 : ComputeMultiPlasticityStress::postReturnMap()
338 : {
339 2795336 : if (_n_supplied)
340 : {
341 : // this is a rotation matrix that will rotate "z" axis back to _n
342 10440 : _rot = _rot.transpose();
343 :
344 : // rotate the tensors back to original frame where _n is correctly oriented
345 10440 : _my_elasticity_tensor.rotate(_rot);
346 10440 : _Jacobian_mult[_qp].rotate(_rot);
347 10440 : _my_strain_increment.rotate(_rot);
348 10440 : _stress[_qp].rotate(_rot);
349 10440 : _plastic_strain[_qp].rotate(_rot);
350 10440 : if (_cosserat)
351 : {
352 0 : _my_flexural_rigidity_tensor.rotate(_rot);
353 0 : (*_Jacobian_mult_couple)[_qp].rotate(_rot);
354 0 : _my_curvature.rotate(_rot);
355 0 : (*_couple_stress)[_qp].rotate(_rot);
356 : }
357 :
358 : // Rotate n by _rotation_increment
359 41760 : for (const auto i : make_range(Moose::dim))
360 : {
361 31320 : _n[_qp](i) = 0;
362 125280 : for (const auto j : make_range(Moose::dim))
363 93960 : _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
364 : }
365 : }
366 2795336 : }
367 :
368 : bool
369 2906940 : ComputeMultiPlasticityStress::quickStep(const RankTwoTensor & stress_old,
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 2906940 : iterations = 0;
386 :
387 : unsigned num_plastic_returns;
388 2906940 : RankTwoTensor delta_dp;
389 :
390 : // the following does the customized returnMap algorithm
391 : // for all the plastic models.
392 2906940 : unsigned custom_model = 0;
393 2906940 : 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 2906940 : 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 259254 : plastic_strain = plastic_strain_old;
416 259254 : if (successful_return && final_step)
417 : {
418 39522 : if (called_from == computeQpStress_function)
419 39522 : 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 0 : consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
424 : }
425 259254 : return successful_return;
426 : }
427 2647686 : 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 2647686 : plastic_strain = plastic_strain_old + delta_dp;
433 2647686 : if (final_step)
434 : {
435 2647686 : if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
436 : {
437 2639566 : if (_tangent_operator_type == elastic)
438 0 : consistent_tangent_operator = E_ijkl;
439 : else
440 : {
441 : std::vector<Real> custom_model_pm;
442 5289212 : for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
443 2649646 : custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
444 : consistent_tangent_operator =
445 2639566 : _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 2639566 : }
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 8120 : consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
458 : }
459 2647686 : return true;
460 : }
461 : else // presumably returnMapAll is incorrectly coded!
462 0 : mooseError("ComputeMultiPlasticityStress::quickStep should not get here!");
463 : }
464 :
465 : bool
466 108128 : ComputeMultiPlasticityStress::plasticStep(const RankTwoTensor & stress_old,
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 : {
481 : /**
482 : * the idea in the following is to potentially subdivide the
483 : * strain increment into smaller fractions, of size "step_size".
484 : * First step_size = 1 is tried, and if that fails then 0.5 is
485 : * tried, then 0.25, etc. As soon as a step is successful, its
486 : * results are put into the "good" variables, which are used
487 : * if a subsequent step fails. If >= 2 consecutive steps are
488 : * successful, the step_size is increased by 1.2
489 : *
490 : * The point of all this is that I hope that the
491 : * time-step for the entire mesh need not be cut if there
492 : * are only a few "bad" quadpoints where the return-map
493 : * is difficult
494 : */
495 :
496 : bool return_successful = false;
497 :
498 : // total number of Newton-Raphson iterations used
499 108128 : unsigned int iter = 0;
500 :
501 108128 : 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 108128 : RankTwoTensor stress_good = stress_old;
507 108128 : RankTwoTensor plastic_strain_good = plastic_strain_old;
508 108128 : std::vector<Real> intnl_good(_num_models);
509 268344 : for (unsigned model = 0; model < _num_models; ++model)
510 160216 : intnl_good[model] = intnl_old[model];
511 108128 : 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 108128 : RankTwoTensor this_strain_increment = strain_increment;
516 :
517 108128 : RankTwoTensor dep = step_size * this_strain_increment;
518 :
519 108128 : _cumulative_pm.assign(_num_surfaces, 0);
520 :
521 : unsigned int num_consecutive_successes = 0;
522 219732 : while (time_simulated < 1.0 && step_size >= _min_stepsize)
523 : {
524 111604 : iter = 0;
525 111604 : 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 111604 : step_size <= _max_stepsize_for_dumb,
536 : linesearch_needed,
537 : ld_encountered,
538 : constraints_added,
539 111604 : time_simulated + step_size >= 1,
540 : consistent_tangent_operator,
541 : _cumulative_pm);
542 111604 : iterations += iter;
543 :
544 111604 : if (return_successful)
545 : {
546 110460 : num_consecutive_successes += 1;
547 : time_simulated += step_size;
548 :
549 110460 : 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 2332 : stress_good = stress;
553 2332 : plastic_strain_good = plastic_strain;
554 10796 : for (unsigned model = 0; model < _num_models; ++model)
555 8464 : intnl_good[model] = intnl[model];
556 10796 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
557 8464 : yf_good[surface] = yf[surface];
558 2332 : if (num_consecutive_successes >= 2)
559 1668 : step_size *= 1.2;
560 : }
561 218732 : step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
562 : }
563 : else
564 : {
565 1144 : step_size *= 0.5;
566 : num_consecutive_successes = 0;
567 1144 : stress = stress_good;
568 1144 : plastic_strain = plastic_strain_good;
569 4856 : for (unsigned model = 0; model < _num_models; ++model)
570 3712 : intnl[model] = intnl_good[model];
571 1144 : yf.resize(_num_surfaces); // might have excited with junk
572 4856 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
573 3712 : yf[surface] = yf_good[surface];
574 1144 : dep = step_size * this_strain_increment;
575 : }
576 : }
577 :
578 108128 : if (!return_successful)
579 : {
580 0 : if (_ignore_failures)
581 : {
582 0 : stress = stress_good;
583 0 : plastic_strain = plastic_strain_good;
584 0 : for (unsigned model = 0; model < _num_models; ++model)
585 0 : intnl[model] = intnl_good[model];
586 0 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
587 0 : yf[surface] = yf_good[surface];
588 : }
589 : else
590 : {
591 : Moose::out << "After reducing the stepsize to " << step_size
592 0 : << " with original strain increment with L2norm " << this_strain_increment.L2norm()
593 : << " the returnMap algorithm failed" << std::endl;
594 :
595 0 : _fspb_debug_stress = stress_good + E_ijkl * dep;
596 0 : _fspb_debug_pm.assign(
597 0 : _num_surfaces,
598 : 1); // this is chosen arbitrarily - please change if a more suitable value occurs to you!
599 0 : _fspb_debug_intnl.resize(_num_models);
600 0 : for (unsigned model = 0; model < _num_models; ++model)
601 0 : _fspb_debug_intnl[model] = intnl_good[model];
602 0 : checkDerivatives();
603 0 : checkJacobian(_my_elasticity_tensor.invSymm(), _intnl_old[_qp]);
604 0 : checkSolution(_my_elasticity_tensor.invSymm());
605 0 : mooseError("Exiting\n");
606 : }
607 : }
608 :
609 108128 : return return_successful;
610 108128 : }
611 :
612 : bool
613 111604 : ComputeMultiPlasticityStress::returnMap(const RankTwoTensor & stress_old,
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 111604 : pm.assign(_num_surfaces, 0.0);
637 :
638 111604 : 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,
651 : returnMap_function,
652 : final_step);
653 :
654 111604 : 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 111604 : 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 111604 : RankTwoTensor initial_stress = stress;
689 :
690 111604 : 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 111604 : buildActiveConstraints(f, stress, intnl, E_ijkl, act);
698 :
699 : // Inverse of E_ijkl (assuming symmetric)
700 111604 : 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 111604 : 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
713 111604 : DeactivationSchemeEnum deact_scheme = _deactivation_scheme;
714 :
715 : // For complicated deactivation schemes we have to record the initial active set
716 : std::vector<bool> initial_act;
717 111604 : initial_act.resize(_num_surfaces);
718 111604 : if (_deactivation_scheme == optimized_to_safe ||
719 110420 : _deactivation_scheme == optimized_to_safe_to_dumb ||
720 : _deactivation_scheme == optimized_to_dumb)
721 : {
722 : // if "optimized" fails we can change the deactivation scheme to "safe", etc
723 1184 : deact_scheme = optimized;
724 8096 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
725 6912 : initial_act[surface] = act[surface];
726 : }
727 :
728 111604 : if (_deactivation_scheme == safe_to_dumb)
729 0 : deact_scheme = safe;
730 :
731 : // For "dumb" deactivation, the active set takes all combinations until a solution is found
732 111604 : int dumb_iteration = 0;
733 : std::vector<unsigned int> dumb_order;
734 :
735 111604 : if (_deactivation_scheme == dumb ||
736 111604 : (_deactivation_scheme == optimized_to_safe_to_dumb && can_revert_to_dumb) ||
737 111604 : (_deactivation_scheme == safe_to_dumb && can_revert_to_dumb) ||
738 0 : (_deactivation_scheme == optimized_to_dumb && can_revert_to_dumb))
739 0 : buildDumbOrder(stress, intnl, dumb_order);
740 :
741 111604 : if (_deactivation_scheme == dumb)
742 : {
743 0 : incrementDumb(dumb_iteration, dumb_order, act);
744 0 : 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 111604 : 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 111604 : Real nr_res2 = 0;
757 369732 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
758 258128 : if (act[surface])
759 154080 : 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 532982 : while (nr_res2 > 0.5 && local_iter++ < _max_iter && single_step_success)
771 414452 : 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 118530 : bool nr_good = (nr_res2 <= 0.5 && local_iter <= _max_iter && single_step_success);
787 :
788 118530 : iter += local_iter;
789 :
790 : // 'act' might have changed due to using deact_scheme = optimized, so
791 118530 : actives_tried.insert(activeCombinationNumber(act));
792 :
793 118530 : 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 444 : change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
803 444 : if (!change_scheme && deact_scheme == dumb)
804 0 : increment_dumb = canIncrementDumb(dumb_iteration);
805 :
806 444 : still_finding_solution = (change_scheme || increment_dumb);
807 :
808 444 : if (change_scheme)
809 0 : 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 444 : if (increment_dumb)
819 0 : incrementDumb(dumb_iteration, dumb_order, act);
820 :
821 444 : 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 118086 : kt_good = checkKuhnTucker(f, pm, act);
834 118086 : if (!kt_good)
835 : {
836 5666 : if (deact_scheme != dumb)
837 : {
838 5666 : applyKuhnTucker(f, pm, act);
839 :
840 : // true if we haven't tried this active set before
841 : still_finding_solution =
842 5666 : (actives_tried.find(activeCombinationNumber(act)) == actives_tried.end());
843 5666 : if (!still_finding_solution)
844 : {
845 : // must have tried turning off the constraints already.
846 : // so try changing the scheme
847 16 : if (canChangeScheme(deact_scheme, can_revert_to_dumb))
848 : {
849 : still_finding_solution = true;
850 0 : 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 0 : increment_dumb = canIncrementDumb(dumb_iteration);
865 : still_finding_solution = increment_dumb;
866 0 : if (increment_dumb)
867 0 : incrementDumb(dumb_iteration, dumb_order, act);
868 : }
869 :
870 5666 : 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 118070 : if (nr_good && kt_good)
882 : {
883 : // check admissible
884 : std::vector<Real> all_f;
885 112420 : 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 42428 : admissible = checkAdmissible(stress, intnl, all_f);
890 :
891 42428 : 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 1960 : bool add_constraints = canAddConstraints(act, all_f);
898 1960 : if (add_constraints)
899 : {
900 1960 : constraints_added = true;
901 1960 : std::vector<bool> act_plus(_num_surfaces,
902 1960 : false); // "act" with the positive constraints added in
903 11592 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
904 9632 : if (act[surface] ||
905 5808 : (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol)))
906 : act_plus[surface] = true;
907 1960 : if (actives_tried.find(activeCombinationNumber(act_plus)) == actives_tried.end())
908 : {
909 : // haven't tried this combination of actives yet
910 924 : constraints_added = true;
911 5516 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
912 4592 : 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 1960 : if (!add_constraints)
922 1036 : change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
923 :
924 1960 : if (!add_constraints && !change_scheme && deact_scheme == dumb)
925 0 : increment_dumb = canIncrementDumb(dumb_iteration);
926 :
927 1960 : still_finding_solution = (add_constraints || change_scheme || increment_dumb);
928 :
929 1960 : if (change_scheme)
930 352 : 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 1960 : if (increment_dumb)
940 0 : incrementDumb(dumb_iteration, dumb_order, act);
941 :
942 1960 : 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 112420 : }
950 :
951 117386 : successful_return = (nr_good && admissible && kt_good);
952 : if (successful_return)
953 : break;
954 :
955 : if (still_finding_solution)
956 : {
957 6926 : stress = initial_stress;
958 6926 : delta_dp = RankTwoTensor(); // back to zero change in plastic strain
959 35518 : for (unsigned model = 0; model < _num_models; ++model)
960 28592 : intnl[model] = intnl_old[model]; // back to old internal params
961 6926 : pm.assign(_num_surfaces, 0.0); // back to zero plastic multipliers
962 :
963 6926 : unsigned num_active = numberActive(act);
964 6926 : if (num_active == 0)
965 : {
966 : successful_return = false;
967 : break; // failure
968 : }
969 :
970 6926 : 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 6926 : yieldFunction(stress, intnl, act, f);
975 :
976 6926 : nr_res2 = 0;
977 : unsigned ind = 0;
978 57106 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
979 50180 : if (act[surface])
980 : {
981 22410 : if (f[ind] > _f[modelNumber(surface)]->_f_tol)
982 22266 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
983 22410 : ind++;
984 : }
985 : }
986 : }
987 :
988 : // returned, with either success or failure
989 110920 : if (successful_return)
990 : {
991 110460 : plastic_strain += delta_dp;
992 :
993 364876 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
994 254416 : cumulative_pm[surface] += pm[surface];
995 :
996 110460 : if (final_step)
997 : consistent_tangent_operator =
998 108128 : consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
999 :
1000 110460 : 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 33948 : yieldFunction(stress, intnl, act, f);
1006 : }
1007 : }
1008 :
1009 : return successful_return;
1010 334812 : }
1011 :
1012 : bool
1013 1960 : ComputeMultiPlasticityStress::canAddConstraints(const std::vector<bool> & act,
1014 : const std::vector<Real> & all_f)
1015 : {
1016 4168 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1017 4168 : if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1018 : return true;
1019 : return false;
1020 : }
1021 :
1022 : bool
1023 1496 : ComputeMultiPlasticityStress::canChangeScheme(DeactivationSchemeEnum current_deactivation_scheme,
1024 : bool can_revert_to_dumb)
1025 : {
1026 1496 : if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
1027 : return true;
1028 :
1029 872 : if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
1030 : return true;
1031 :
1032 1144 : if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
1033 : can_revert_to_dumb)
1034 : return true;
1035 :
1036 304 : if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
1037 : can_revert_to_dumb)
1038 : return true;
1039 :
1040 1144 : 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 352 : 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 352 : if (current_deactivation_scheme == optimized &&
1058 352 : (_deactivation_scheme == optimized_to_safe ||
1059 : _deactivation_scheme == optimized_to_safe_to_dumb))
1060 : {
1061 352 : current_deactivation_scheme = safe;
1062 2656 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1063 2304 : act[surface] = initial_act[surface];
1064 : }
1065 0 : else if ((current_deactivation_scheme == safe &&
1066 0 : (_deactivation_scheme == optimized_to_safe_to_dumb ||
1067 : _deactivation_scheme == safe_to_dumb) &&
1068 0 : can_revert_to_dumb) ||
1069 0 : (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_dumb &&
1070 : can_revert_to_dumb))
1071 : {
1072 0 : current_deactivation_scheme = dumb;
1073 0 : dumb_iteration = 0;
1074 0 : buildDumbOrder(initial_stress, intnl_old, dumb_order);
1075 0 : incrementDumb(dumb_iteration, dumb_order, act);
1076 : }
1077 352 : }
1078 :
1079 : bool
1080 414452 : ComputeMultiPlasticityStress::singleStep(Real & nr_res2,
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 414452 : Real nr_res2_before_step = nr_res2;
1098 414452 : RankTwoTensor stress_before_step;
1099 : std::vector<Real> intnl_before_step;
1100 : std::vector<Real> pm_before_step;
1101 414452 : RankTwoTensor delta_dp_before_step;
1102 :
1103 414452 : if (deactivation_scheme == optimized)
1104 : {
1105 : // we potentially use the "before_step" quantities, so record them here
1106 55330 : stress_before_step = stress;
1107 55330 : intnl_before_step.resize(_num_models);
1108 272747 : for (unsigned model = 0; model < _num_models; ++model)
1109 217417 : intnl_before_step[model] = intnl[model];
1110 55330 : pm_before_step.resize(_num_surfaces);
1111 272747 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1112 217417 : pm_before_step[surface] = pm[surface];
1113 55330 : 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 414452 : 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 414452 : 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 837566 : while (constraints_changing)
1138 : {
1139 : // calculate dstress, dpm and dintnl for one full Newton-Raphson step
1140 423138 : nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
1141 :
1142 1130036 : for (unsigned surface = 0; surface < deact_ld.size(); ++surface)
1143 719988 : if (deact_ld[surface])
1144 : {
1145 13090 : ld_encountered = true;
1146 13090 : break;
1147 : }
1148 :
1149 : // perform a line search
1150 : // The line-search will exit with updated values
1151 423138 : 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 423138 : 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 423114 : if (deactivation_scheme == optimized)
1175 : {
1176 321169 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1177 257177 : if (active[surface] && pm[surface] < 0.0)
1178 : constraints_changing = true;
1179 : }
1180 :
1181 63992 : if (constraints_changing)
1182 : {
1183 8686 : stress = stress_before_step;
1184 8686 : delta_dp = delta_dp_before_step;
1185 8686 : nr_res2 = nr_res2_before_step;
1186 48542 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1187 : {
1188 39856 : if (active[surface] && pm[surface] < 0.0)
1189 : {
1190 : // turn off the constraint forever
1191 : active[surface] = false;
1192 9794 : pm_before_step[surface] = 0.0;
1193 9794 : intnl_before_step[modelNumber(surface)] =
1194 9794 : intnl_old[modelNumber(surface)]; // don't want to muck-up hardening!
1195 : }
1196 39856 : intnl[modelNumber(surface)] = intnl_before_step[modelNumber(surface)];
1197 39856 : pm[surface] = pm_before_step[surface];
1198 : }
1199 8686 : 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 423114 : 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 414428 : if (reinstated_actives)
1214 : {
1215 : bool completely_converged = true;
1216 10164 : 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 49734 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1231 43696 : if (active[surface])
1232 27184 : if (f[ind++] > _f[modelNumber(surface)]->_f_tol)
1233 : completely_converged = false;
1234 : }
1235 : else
1236 : completely_converged = false;
1237 :
1238 6038 : if (!completely_converged)
1239 4126 : nr_res2 = residual2(pm, f, epp, ic, active, deact_ld);
1240 : }
1241 :
1242 : return successful_step;
1243 414452 : }
1244 :
1245 : bool
1246 423114 : ComputeMultiPlasticityStress::reinstateLinearDependentConstraints(
1247 : std::vector<bool> & deactivated_due_to_ld)
1248 : {
1249 : bool reinstated_actives = false;
1250 1197439 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1251 774325 : if (deactivated_due_to_ld[surface])
1252 : reinstated_actives = true;
1253 :
1254 423114 : deactivated_due_to_ld.assign(_num_surfaces, false);
1255 423114 : return reinstated_actives;
1256 : }
1257 :
1258 : unsigned
1259 15612 : ComputeMultiPlasticityStress::numberActive(const std::vector<bool> & active)
1260 : {
1261 : unsigned num_active = 0;
1262 105648 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1263 90036 : if (active[surface])
1264 43524 : num_active++;
1265 15612 : return num_active;
1266 : }
1267 :
1268 : bool
1269 42428 : ComputeMultiPlasticityStress::checkAdmissible(const RankTwoTensor & stress,
1270 : const std::vector<Real> & intnl,
1271 : std::vector<Real> & all_f)
1272 : {
1273 : std::vector<bool> act;
1274 42428 : act.assign(_num_surfaces, true);
1275 :
1276 42428 : yieldFunction(stress, intnl, act, all_f);
1277 :
1278 229060 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1279 188592 : if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1280 : return false;
1281 :
1282 : return true;
1283 : }
1284 :
1285 : bool
1286 118086 : ComputeMultiPlasticityStress::checkKuhnTucker(const std::vector<Real> & f,
1287 : const std::vector<Real> & pm,
1288 : const std::vector<bool> & active)
1289 : {
1290 : unsigned ind = 0;
1291 419102 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1292 : {
1293 303414 : if (active[surface])
1294 : {
1295 163690 : if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1296 9938 : if (pm[surface] != 0)
1297 : return false;
1298 : }
1299 139724 : else if (pm[surface] != 0)
1300 0 : mooseError("Crash due to plastic multiplier not being zero. This occurred because of poor "
1301 : "coding!!");
1302 : }
1303 395100 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1304 282680 : if (pm[surface] < 0)
1305 : return false;
1306 :
1307 : return true;
1308 : }
1309 :
1310 : void
1311 5666 : ComputeMultiPlasticityStress::applyKuhnTucker(const std::vector<Real> & f,
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 49014 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1320 : {
1321 43348 : if (active[surface])
1322 : {
1323 23364 : if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1324 8920 : if (pm[surface] != 0)
1325 : {
1326 : turned_off = true;
1327 : active[surface] = false;
1328 : }
1329 : }
1330 19984 : else if (pm[surface] != 0)
1331 0 : 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 5666 : if (!turned_off)
1337 : {
1338 : int surface_to_turn_off = -1;
1339 : Real min_pm = 0;
1340 31424 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1341 28156 : if (pm[surface] < min_pm)
1342 : {
1343 : min_pm = pm[surface];
1344 3268 : surface_to_turn_off = surface;
1345 : }
1346 3268 : if (surface_to_turn_off >= 0)
1347 : active[surface_to_turn_off] = false;
1348 : }
1349 5666 : }
1350 :
1351 : Real
1352 615454 : 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 1693409 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1363 1077955 : if (active[surface])
1364 : {
1365 800982 : if (!deactivated_due_to_ld[surface])
1366 : {
1367 756316 : if (!(pm[surface] == 0 && f[ind] <= 0))
1368 748486 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1369 : }
1370 44666 : else if (deactivated_due_to_ld[surface] && f[ind] > 0)
1371 15628 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1372 800982 : ind++;
1373 : }
1374 :
1375 615454 : nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
1376 :
1377 615454 : std::vector<bool> active_not_deact(_num_surfaces);
1378 1693409 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1379 1399594 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1380 : ind = 0;
1381 1586085 : for (unsigned model = 0; model < _num_models; ++model)
1382 970631 : if (anyActiveSurfaces(model, active_not_deact))
1383 747404 : nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
1384 :
1385 615454 : return nr_res2;
1386 : }
1387 :
1388 : bool
1389 423138 : ComputeMultiPlasticityStress::lineSearch(Real & nr_res2,
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 423138 : 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 423138 : Real f0 = nr_res2; // initial value of residual2
1417 423138 : 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 423138 : ls_pm.resize(pm.size());
1426 :
1427 : // delta_dp during the line-search
1428 423138 : RankTwoTensor ls_delta_dp;
1429 :
1430 : // internal parameter during the line-search
1431 : std::vector<Real> ls_intnl;
1432 423138 : ls_intnl.resize(intnl.size());
1433 :
1434 : // stress during the line-search
1435 423138 : 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 1662906 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1444 1051602 : ls_pm[alpha] = pm[alpha] + dpm[alpha] * lam;
1445 611304 : ls_delta_dp = delta_dp - E_inv * dstress * lam;
1446 1555582 : for (unsigned a = 0; a < intnl.size(); ++a)
1447 944278 : ls_intnl[a] = intnl[a] + dintnl[a] * lam;
1448 611304 : ls_stress = stress + dstress * lam;
1449 :
1450 : // calculate the new active yield functions, epp and active internal constraints
1451 611304 : 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 611304 : nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1455 :
1456 611304 : if (nr_res2 < f0 + 1E-4 * lam * slope)
1457 : break;
1458 188190 : else if (lam < lam_min)
1459 : {
1460 : success = false;
1461 : // restore plastic multipliers, yield functions, etc to original values
1462 120 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1463 96 : ls_pm[alpha] = pm[alpha];
1464 24 : ls_delta_dp = delta_dp;
1465 120 : for (unsigned a = 0; a < intnl.size(); ++a)
1466 96 : ls_intnl[a] = intnl[a];
1467 24 : ls_stress = stress;
1468 24 : calculateConstraints(
1469 : ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1470 24 : nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1471 24 : break;
1472 : }
1473 188166 : else if (lam == 1.0)
1474 : {
1475 : // model as a quadratic
1476 97970 : tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
1477 : }
1478 : else
1479 : {
1480 : // model as a cubic
1481 90196 : Real rhs1 = nr_res2 - f0 - lam * slope;
1482 90196 : Real rhs2 = f2 - f0 - lam2 * slope;
1483 90196 : Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1484 : Real b =
1485 90196 : (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1486 90196 : if (a == 0)
1487 0 : tmp_lam = -slope / (2.0 * b);
1488 : else
1489 : {
1490 90196 : Real disc = Utility::pow<2>(b) - 3 * a * slope;
1491 90196 : if (disc < 0)
1492 0 : tmp_lam = 0.5 * lam;
1493 90196 : else if (b <= 0)
1494 3924 : tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
1495 : else
1496 86272 : tmp_lam = -slope / (b + std::sqrt(disc));
1497 : }
1498 90196 : if (tmp_lam > 0.5 * lam)
1499 3288 : tmp_lam = 0.5 * lam;
1500 : }
1501 : lam2 = lam;
1502 : f2 = nr_res2;
1503 188166 : lam = std::max(tmp_lam, 0.1 * lam);
1504 188166 : }
1505 :
1506 423138 : if (lam < 1.0)
1507 97970 : linesearch_needed = true;
1508 :
1509 : // assign the quantities found in the line-search
1510 : // back to the originals
1511 1197559 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1512 774421 : pm[alpha] = ls_pm[alpha];
1513 423138 : delta_dp = ls_delta_dp;
1514 1090235 : for (unsigned a = 0; a < intnl.size(); ++a)
1515 667097 : intnl[a] = ls_intnl[a];
1516 423138 : stress = ls_stress;
1517 :
1518 423138 : return success;
1519 423138 : }
1520 :
1521 : void
1522 0 : ComputeMultiPlasticityStress::buildDumbOrder(const RankTwoTensor & stress,
1523 : const std::vector<Real> & intnl,
1524 : std::vector<unsigned int> & dumb_order)
1525 : {
1526 0 : if (dumb_order.size() != 0)
1527 0 : return;
1528 :
1529 : std::vector<bool> act;
1530 0 : act.assign(_num_surfaces, true);
1531 :
1532 : std::vector<Real> f;
1533 0 : yieldFunction(stress, intnl, act, f);
1534 : std::vector<RankTwoTensor> df_dstress;
1535 0 : dyieldFunction_dstress(stress, intnl, act, df_dstress);
1536 :
1537 : typedef std::pair<Real, unsigned> pair_for_sorting;
1538 0 : std::vector<pair_for_sorting> dist(_num_surfaces);
1539 0 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1540 : {
1541 0 : dist[surface].first = f[surface] / df_dstress[surface].L2norm();
1542 0 : dist[surface].second = surface;
1543 : }
1544 0 : std::sort(dist.begin(), dist.end()); // sorted in ascending order of f/df_dstress
1545 :
1546 0 : dumb_order.resize(_num_surfaces);
1547 0 : for (unsigned i = 0; i < _num_surfaces; ++i)
1548 0 : dumb_order[i] = dist[_num_surfaces - 1 - i].second;
1549 : // now dumb_order[0] is the surface with the greatest f/df_dstress
1550 0 : }
1551 :
1552 : void
1553 0 : ComputeMultiPlasticityStress::incrementDumb(int & dumb_iteration,
1554 : const std::vector<unsigned int> & dumb_order,
1555 : std::vector<bool> & act)
1556 : {
1557 0 : dumb_iteration += 1;
1558 0 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1559 0 : act[dumb_order[surface]] =
1560 0 : (dumb_iteration &
1561 0 : (1 << surface)); // returns true if the surface_th bit of dumb_iteration == 1
1562 0 : }
1563 :
1564 : bool
1565 0 : ComputeMultiPlasticityStress::canIncrementDumb(int dumb_iteration)
1566 : {
1567 : // (1 << _num_surfaces) = 2^_num_surfaces
1568 0 : return ((dumb_iteration + 1) < (1 << _num_surfaces));
1569 : }
1570 :
1571 : unsigned int
1572 244686 : ComputeMultiPlasticityStress::activeCombinationNumber(const std::vector<bool> & act)
1573 : {
1574 : unsigned num = 0;
1575 914282 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1576 669596 : if (act[surface])
1577 366540 : num += (1 << surface); // (1 << x) = 2^x
1578 :
1579 244686 : return num;
1580 : }
1581 :
1582 : RankFourTensor
1583 116248 : ComputeMultiPlasticityStress::consistentTangentOperator(const RankTwoTensor & stress,
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 :
1590 116248 : if (_tangent_operator_type == elastic)
1591 5024 : 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 111224 : std::vector<bool> act_at_some_step(_num_surfaces);
1597 399336 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1598 288112 : 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 111224 : std::vector<bool> act_vary(_num_surfaces);
1605 399336 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1606 288112 : act_vary[surface] = (pm_this_step[surface] > 0);
1607 :
1608 : std::vector<RankTwoTensor> df_dstress;
1609 111224 : dyieldFunction_dstress(stress, intnl, act_vary, df_dstress);
1610 : std::vector<Real> df_dintnl;
1611 111224 : dyieldFunction_dintnl(stress, intnl, act_vary, df_dintnl);
1612 : std::vector<RankTwoTensor> r;
1613 111224 : flowPotential(stress, intnl, act_vary, r);
1614 : std::vector<RankFourTensor> dr_dstress_at_some_step;
1615 111224 : dflowPotential_dstress(stress, intnl, act_at_some_step, dr_dstress_at_some_step);
1616 : std::vector<RankTwoTensor> dr_dintnl_at_some_step;
1617 111224 : dflowPotential_dintnl(stress, intnl, act_at_some_step, dr_dintnl_at_some_step);
1618 : std::vector<Real> h;
1619 111224 : 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 399336 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1630 288112 : if (act_vary[surface1])
1631 : {
1632 139104 : r_minus_stuff.push_back(r[ind1]);
1633 : ind2 = 0;
1634 594060 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1635 454956 : if (act_at_some_step[surface2])
1636 : {
1637 207056 : if (modelNumber(surface1) == modelNumber(surface2))
1638 : {
1639 : r_minus_stuff.back() -=
1640 174424 : cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
1641 : }
1642 207056 : ind2++;
1643 : }
1644 139104 : ind1++;
1645 : }
1646 :
1647 : unsigned int num_currently_active = 0;
1648 399336 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1649 288112 : if (act_vary[surface])
1650 139104 : 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 111224 : zzz.assign(num_currently_active * num_currently_active, 0.0);
1660 :
1661 : ind1 = 0;
1662 111224 : RankTwoTensor r2;
1663 399336 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1664 288112 : if (act_vary[surface1])
1665 : {
1666 : ind2 = 0;
1667 594060 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1668 454956 : if (act_vary[surface2])
1669 : {
1670 207056 : r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
1671 207056 : zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
1672 207056 : if (modelNumber(surface1) == modelNumber(surface2))
1673 174424 : zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
1674 207056 : ind2++;
1675 : }
1676 139104 : ind1++;
1677 : }
1678 :
1679 111224 : 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 111224 : MatrixTools::inverse(zzz, num_currently_active);
1685 : }
1686 0 : catch (const MooseException & e)
1687 : {
1688 : // in the very rare case of zzz being singular, just return the "elastic" tangent operator
1689 0 : return E_ijkl;
1690 0 : }
1691 : }
1692 :
1693 111224 : RankFourTensor strain_coeff = E_ijkl;
1694 : ind1 = 0;
1695 399336 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1696 288112 : if (act_vary[surface1])
1697 : {
1698 139104 : RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1699 : ind2 = 0;
1700 594060 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1701 454956 : if (act_vary[surface2])
1702 : {
1703 207056 : RankTwoTensor part2 = E_ijkl * df_dstress[ind2];
1704 828224 : for (unsigned i = 0; i < 3; i++)
1705 2484672 : for (unsigned j = 0; j < 3; j++)
1706 7454016 : for (unsigned k = 0; k < 3; k++)
1707 22362048 : for (unsigned l = 0; l < 3; l++)
1708 16771536 : strain_coeff(i, j, k, l) -=
1709 16771536 : part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1710 207056 : ind2++;
1711 : }
1712 139104 : ind1++;
1713 : }
1714 :
1715 111224 : if (_tangent_operator_type == linear)
1716 4512 : return strain_coeff;
1717 :
1718 106712 : RankFourTensor stress_coeff(RankFourTensor::initIdentitySymmetricFour);
1719 :
1720 106712 : RankFourTensor part3;
1721 : ind1 = 0;
1722 379560 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1723 272848 : if (act_at_some_step[surface1])
1724 : {
1725 265600 : part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
1726 132800 : ind1++;
1727 : }
1728 :
1729 106712 : stress_coeff += part3;
1730 :
1731 106712 : 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 379560 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1736 272848 : if (act_vary[surface1])
1737 : {
1738 132800 : RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1739 : ind2 = 0;
1740 563532 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1741 430732 : if (act_vary[surface2])
1742 : {
1743 197168 : RankTwoTensor part2 = part3 * df_dstress[ind2];
1744 788672 : for (unsigned i = 0; i < 3; i++)
1745 2366016 : for (unsigned j = 0; j < 3; j++)
1746 7098048 : for (unsigned k = 0; k < 3; k++)
1747 21294144 : for (unsigned l = 0; l < 3; l++)
1748 15970608 : stress_coeff(i, j, k, l) -=
1749 15970608 : part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1750 197168 : ind2++;
1751 : }
1752 132800 : 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 106712 : RankFourTensor s_inv;
1764 : try
1765 : {
1766 106712 : s_inv = stress_coeff.invSymm();
1767 : }
1768 219 : catch (const MooseException & e)
1769 : {
1770 219 : return strain_coeff; // when stress_coeff is singular (perhaps for incompressible plasticity?)
1771 : // return the "linear" tangent operator
1772 219 : }
1773 :
1774 106493 : return s_inv * strain_coeff;
1775 111224 : }
|
