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 4174 : ComputeMultiPlasticityStress::validParams()
23 : {
24 4174 : InputParameters params = ComputeStressBase::validParams();
25 4174 : params += MultiPlasticityDebugger::validParams();
26 4174 : params.addClassDescription("Base class for multi-surface finite-strain plasticity");
27 12522 : params.addRangeCheckedParam<unsigned int>("max_NR_iterations",
28 8348 : 20,
29 : "max_NR_iterations>0",
30 : "Maximum number of Newton-Raphson iterations allowed");
31 8348 : 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 12522 : params.addRangeCheckedParam<Real>(
36 : "min_stepsize",
37 8348 : 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 12522 : params.addRangeCheckedParam<Real>("max_stepsize_for_dumb",
43 8348 : 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 8348 : "optimized");
52 8348 : 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 8348 : 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 8348 : params.addParam<bool>("ignore_failures",
71 8348 : 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 8348 : MooseEnum tangent_operator("elastic linear nonlinear", "nonlinear");
78 8348 : 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 8348 : params.addParam<bool>("perform_finite_strain_rotations",
87 8348 : 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 4174 : params.addClassDescription("Material for multi-surface finite-strain plasticity");
94 4174 : return params;
95 4174 : }
96 :
97 3122 : ComputeMultiPlasticityStress::ComputeMultiPlasticityStress(const InputParameters & parameters)
98 : : ComputeStressBase(parameters),
99 : MultiPlasticityDebugger(this),
100 3122 : _max_iter(getParam<unsigned int>("max_NR_iterations")),
101 6244 : _min_stepsize(getParam<Real>("min_stepsize")),
102 6244 : _max_stepsize_for_dumb(getParam<Real>("max_stepsize_for_dumb")),
103 6244 : _ignore_failures(getParam<bool>("ignore_failures")),
104 :
105 6244 : _tangent_operator_type((TangentOperatorEnum)(int)getParam<MooseEnum>("tangent_operator")),
106 :
107 6244 : _epp_tol(getParam<Real>("ep_plastic_tolerance")),
108 :
109 3122 : _dummy_pm(0),
110 :
111 3122 : _cumulative_pm(0),
112 :
113 6244 : _deactivation_scheme((DeactivationSchemeEnum)(int)getParam<MooseEnum>("deactivation_scheme")),
114 :
115 3122 : _n_supplied(parameters.isParamValid("transverse_direction")),
116 6612 : _n_input(_n_supplied ? getParam<RealVectorValue>("transverse_direction") : RealVectorValue()),
117 : _rot(RealTensorValue()),
118 :
119 6244 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
120 :
121 3122 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
122 3122 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
123 3122 : _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
124 6244 : _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
125 3122 : _intnl(declareProperty<std::vector<Real>>("plastic_internal_parameter")),
126 6244 : _intnl_old(getMaterialPropertyOld<std::vector<Real>>("plastic_internal_parameter")),
127 3122 : _yf(declareProperty<std::vector<Real>>("plastic_yield_function")),
128 3122 : _iter(declareProperty<Real>("plastic_NR_iterations")), // this is really an unsigned int, but
129 : // for visualisation i convert it to Real
130 3122 : _linesearch_needed(declareProperty<Real>("plastic_linesearch_needed")), // this is really a
131 : // boolean, but for
132 : // visualisation i
133 : // convert it to Real
134 3122 : _ld_encountered(declareProperty<Real>(
135 : "plastic_linear_dependence_encountered")), // this is really a boolean, but for
136 : // visualisation i convert it to Real
137 3122 : _constraints_added(declareProperty<Real>("plastic_constraints_added")), // this is really a
138 : // boolean, but for
139 : // visualisation i
140 : // convert it to Real
141 3122 : _n(declareProperty<RealVectorValue>("plastic_transverse_direction")),
142 6244 : _n_old(getMaterialPropertyOld<RealVectorValue>("plastic_transverse_direction")),
143 :
144 6244 : _strain_increment(getMaterialPropertyByName<RankTwoTensor>(_base_name + "strain_increment")),
145 6244 : _total_strain_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "total_strain")),
146 3122 : _rotation_increment(
147 3122 : getMaterialPropertyByName<RankTwoTensor>(_base_name + "rotation_increment")),
148 :
149 6244 : _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
150 6244 : _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 3122 : _cosserat(hasMaterialProperty<RankTwoTensor>("curvature") &&
155 3140 : hasMaterialProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
156 3140 : _curvature(_cosserat ? &getMaterialPropertyByName<RankTwoTensor>("curvature") : nullptr),
157 3122 : _elastic_flexural_rigidity_tensor(
158 3140 : _cosserat ? &getMaterialPropertyByName<RankFourTensor>("elastic_flexural_rigidity_tensor")
159 : : nullptr),
160 3122 : _couple_stress(_cosserat ? &declareProperty<RankTwoTensor>("couple_stress") : nullptr),
161 3140 : _couple_stress_old(_cosserat ? &getMaterialPropertyOld<RankTwoTensor>("couple_stress")
162 : : nullptr),
163 3122 : _Jacobian_mult_couple(_cosserat ? &declareProperty<RankFourTensor>("couple_Jacobian_mult")
164 : : nullptr),
165 :
166 3122 : _my_elasticity_tensor(RankFourTensor()),
167 3122 : _my_strain_increment(RankTwoTensor()),
168 3122 : _my_flexural_rigidity_tensor(RankFourTensor()),
169 6244 : _my_curvature(RankTwoTensor())
170 : {
171 3122 : if (_epp_tol <= 0)
172 0 : mooseError("ComputeMultiPlasticityStress: ep_plastic_tolerance must be positive");
173 :
174 3122 : if (_n_supplied)
175 : {
176 : // normalise the inputted transverse_direction
177 368 : if (_n_input.norm() == 0)
178 2 : mooseError(
179 : "ComputeMultiPlasticityStress: transverse_direction vector must not have zero length");
180 : else
181 366 : _n_input /= _n_input.norm();
182 : }
183 :
184 3120 : if (_num_surfaces == 1)
185 1842 : _deactivation_scheme = safe;
186 3120 : }
187 :
188 : void
189 87776 : ComputeMultiPlasticityStress::initQpStatefulProperties()
190 : {
191 87776 : ComputeStressBase::initQpStatefulProperties();
192 :
193 87776 : _plastic_strain[_qp].zero();
194 :
195 87776 : _intnl[_qp].assign(_num_models, 0);
196 :
197 87776 : _yf[_qp].assign(_num_surfaces, 0);
198 :
199 87776 : _dummy_pm.assign(_num_surfaces, 0);
200 :
201 87776 : _iter[_qp] = 0.0; // this is really an unsigned int, but for visualisation i convert it to Real
202 87776 : _linesearch_needed[_qp] = 0;
203 87776 : _ld_encountered[_qp] = 0;
204 87776 : _constraints_added[_qp] = 0;
205 :
206 87776 : _n[_qp] = _n_input;
207 :
208 87776 : if (_cosserat)
209 6400 : (*_couple_stress)[_qp].zero();
210 :
211 87776 : if (_fspb_debug == "jacobian")
212 : {
213 0 : checkDerivatives();
214 0 : mooseError("Finite-differencing completed. Exiting with no error");
215 : }
216 87776 : }
217 :
218 : void
219 6658928 : ComputeMultiPlasticityStress::computeQpStress()
220 : {
221 : // the following "_my" variables can get rotated by preReturnMap and postReturnMap
222 6658928 : _my_elasticity_tensor = _elasticity_tensor[_qp];
223 6658928 : _my_strain_increment = _strain_increment[_qp];
224 6658928 : if (_cosserat)
225 : {
226 19200 : _my_flexural_rigidity_tensor = (*_elastic_flexural_rigidity_tensor)[_qp];
227 19200 : _my_curvature = (*_curvature)[_qp];
228 : }
229 :
230 6658928 : 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 6658928 : preReturnMap(); // do rotations to new frame if necessary
239 :
240 : unsigned int number_iterations;
241 6658928 : bool linesearch_needed = false;
242 6658928 : bool ld_encountered = false;
243 6658928 : bool constraints_added = false;
244 :
245 6658928 : _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 6658928 : const bool found_solution = quickStep(rot(_stress_old[_qp]),
249 6658928 : _stress[_qp],
250 6658928 : _intnl_old[_qp],
251 6658928 : _intnl[_qp],
252 6658928 : _dummy_pm,
253 : _cumulative_pm,
254 13317856 : rot(_plastic_strain_old[_qp]),
255 6658928 : _plastic_strain[_qp],
256 6658928 : _my_elasticity_tensor,
257 6658928 : _my_strain_increment,
258 6658928 : _yf[_qp],
259 : number_iterations,
260 6658928 : _Jacobian_mult[_qp],
261 : computeQpStress_function,
262 : true);
263 :
264 : // if not purely elastic or the customised stuff failed, do some plastic return
265 6658928 : if (!found_solution)
266 138032 : plasticStep(rot(_stress_old[_qp]),
267 138032 : _stress[_qp],
268 138032 : _intnl_old[_qp],
269 138032 : _intnl[_qp],
270 138032 : rot(_plastic_strain_old[_qp]),
271 138032 : _plastic_strain[_qp],
272 : _my_elasticity_tensor,
273 : _my_strain_increment,
274 138032 : _yf[_qp],
275 : number_iterations,
276 : linesearch_needed,
277 : ld_encountered,
278 : constraints_added,
279 138032 : _Jacobian_mult[_qp]);
280 :
281 6658928 : 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 6658928 : postReturnMap(); // rotate back from new frame if necessary
288 :
289 6658928 : _iter[_qp] = 1.0 * number_iterations;
290 6658928 : _linesearch_needed[_qp] = linesearch_needed;
291 6658928 : _ld_encountered[_qp] = ld_encountered;
292 6658928 : _constraints_added[_qp] = constraints_added;
293 :
294 : // Update measures of strain
295 6658928 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + _my_strain_increment -
296 6658928 : (_plastic_strain[_qp] - _plastic_strain_old[_qp]);
297 :
298 : // Rotate the tensors to the current configuration
299 6658928 : if (_perform_finite_strain_rotations)
300 : {
301 6658928 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
302 6658928 : _elastic_strain[_qp] =
303 6658928 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
304 6658928 : _plastic_strain[_qp] =
305 6658928 : _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
306 : }
307 6658928 : }
308 :
309 : RankTwoTensor
310 13593920 : ComputeMultiPlasticityStress::rot(const RankTwoTensor & tens)
311 : {
312 13593920 : if (!_n_supplied)
313 13541504 : return tens;
314 52416 : return tens.rotated(_rot);
315 : }
316 :
317 : void
318 6658928 : ComputeMultiPlasticityStress::preReturnMap()
319 : {
320 6658928 : if (_n_supplied)
321 : {
322 : // this is a rotation matrix that will rotate _n to the "z" axis
323 14304 : _rot = RotationMatrix::rotVecToZ(_n[_qp]);
324 :
325 : // rotate the tensors to this frame
326 14304 : _my_elasticity_tensor.rotate(_rot);
327 14304 : _my_strain_increment.rotate(_rot);
328 14304 : if (_cosserat)
329 : {
330 0 : _my_flexural_rigidity_tensor.rotate(_rot);
331 0 : _my_curvature.rotate(_rot);
332 : }
333 : }
334 6658928 : }
335 :
336 : void
337 6658928 : ComputeMultiPlasticityStress::postReturnMap()
338 : {
339 6658928 : if (_n_supplied)
340 : {
341 : // this is a rotation matrix that will rotate "z" axis back to _n
342 14304 : _rot = _rot.transpose();
343 :
344 : // rotate the tensors back to original frame where _n is correctly oriented
345 14304 : _my_elasticity_tensor.rotate(_rot);
346 14304 : _Jacobian_mult[_qp].rotate(_rot);
347 14304 : _my_strain_increment.rotate(_rot);
348 14304 : _stress[_qp].rotate(_rot);
349 14304 : _plastic_strain[_qp].rotate(_rot);
350 14304 : 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 57216 : for (const auto i : make_range(Moose::dim))
360 : {
361 42912 : _n[_qp](i) = 0;
362 171648 : for (const auto j : make_range(Moose::dim))
363 128736 : _n[_qp](i) += _rotation_increment[_qp](i, j) * _n_old[_qp](j);
364 : }
365 : }
366 6658928 : }
367 :
368 : bool
369 6801584 : 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 6801584 : iterations = 0;
386 :
387 : unsigned num_plastic_returns;
388 6801584 : RankTwoTensor delta_dp;
389 :
390 : // the following does the customized returnMap algorithm
391 : // for all the plastic models.
392 6801584 : unsigned custom_model = 0;
393 6801584 : 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 6801584 : 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 326500 : plastic_strain = plastic_strain_old;
416 326500 : if (successful_return && final_step)
417 : {
418 45812 : if (called_from == computeQpStress_function)
419 45812 : 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 326500 : return successful_return;
426 : }
427 6475084 : 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 6475084 : plastic_strain = plastic_strain_old + delta_dp;
433 6475084 : if (final_step)
434 : {
435 6475084 : if (called_from == computeQpStress_function && _f[custom_model]->useCustomCTO())
436 : {
437 6466028 : if (_tangent_operator_type == elastic)
438 0 : consistent_tangent_operator = E_ijkl;
439 : else
440 : {
441 : std::vector<Real> custom_model_pm;
442 12949432 : for (unsigned surface = 0; surface < _f[custom_model]->numberSurfaces(); ++surface)
443 6483404 : custom_model_pm.push_back(cumulative_pm[_surfaces_given_model[custom_model][surface]]);
444 : consistent_tangent_operator =
445 6466028 : _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 9056 : consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
458 : }
459 6475084 : return true;
460 : }
461 : else // presumably returnMapAll is incorrectly coded!
462 0 : mooseError("ComputeMultiPlasticityStress::quickStep should not get here!");
463 : }
464 :
465 : bool
466 138032 : 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 138032 : unsigned int iter = 0;
500 :
501 138032 : 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 138032 : RankTwoTensor stress_good = stress_old;
507 138032 : RankTwoTensor plastic_strain_good = plastic_strain_old;
508 138032 : std::vector<Real> intnl_good(_num_models);
509 348448 : for (unsigned model = 0; model < _num_models; ++model)
510 210416 : intnl_good[model] = intnl_old[model];
511 138032 : 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 138032 : RankTwoTensor this_strain_increment = strain_increment;
516 :
517 138032 : RankTwoTensor dep = step_size * this_strain_increment;
518 :
519 138032 : _cumulative_pm.assign(_num_surfaces, 0);
520 :
521 : unsigned int num_consecutive_successes = 0;
522 280688 : while (time_simulated < 1.0 && step_size >= _min_stepsize)
523 : {
524 142656 : iter = 0;
525 142656 : 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 142656 : step_size <= _max_stepsize_for_dumb,
536 : linesearch_needed,
537 : ld_encountered,
538 : constraints_added,
539 142656 : time_simulated + step_size >= 1,
540 : consistent_tangent_operator,
541 : _cumulative_pm);
542 142656 : iterations += iter;
543 :
544 142656 : if (return_successful)
545 : {
546 141136 : num_consecutive_successes += 1;
547 : time_simulated += step_size;
548 :
549 141136 : 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 3104 : stress_good = stress;
553 3104 : plastic_strain_good = plastic_strain;
554 14368 : for (unsigned model = 0; model < _num_models; ++model)
555 11264 : intnl_good[model] = intnl[model];
556 14368 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
557 11264 : yf_good[surface] = yf[surface];
558 3104 : if (num_consecutive_successes >= 2)
559 2224 : step_size *= 1.2;
560 : }
561 279360 : step_size = std::min(step_size, 1.0 - time_simulated); // avoid overshoots
562 : }
563 : else
564 : {
565 1520 : step_size *= 0.5;
566 : num_consecutive_successes = 0;
567 1520 : stress = stress_good;
568 1520 : plastic_strain = plastic_strain_good;
569 6448 : for (unsigned model = 0; model < _num_models; ++model)
570 4928 : intnl[model] = intnl_good[model];
571 1520 : yf.resize(_num_surfaces); // might have excited with junk
572 6448 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
573 4928 : yf[surface] = yf_good[surface];
574 1520 : dep = step_size * this_strain_increment;
575 : }
576 : }
577 :
578 138032 : 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 138032 : return return_successful;
610 : }
611 :
612 : bool
613 142656 : 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 142656 : pm.assign(_num_surfaces, 0.0);
637 :
638 142656 : 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 142656 : 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 142656 : 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 142656 : RankTwoTensor initial_stress = stress;
689 :
690 142656 : 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 142656 : buildActiveConstraints(f, stress, intnl, E_ijkl, act);
698 :
699 : // Inverse of E_ijkl (assuming symmetric)
700 142656 : 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 142656 : 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 142656 : 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 142656 : initial_act.resize(_num_surfaces);
718 142656 : if (_deactivation_scheme == optimized_to_safe ||
719 141024 : _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 1632 : deact_scheme = optimized;
724 12128 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
725 10496 : initial_act[surface] = act[surface];
726 : }
727 :
728 142656 : 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 142656 : int dumb_iteration = 0;
733 : std::vector<unsigned int> dumb_order;
734 :
735 142656 : if (_deactivation_scheme == dumb ||
736 142656 : (_deactivation_scheme == optimized_to_safe_to_dumb && can_revert_to_dumb) ||
737 142656 : (_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 142656 : 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 142656 : 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 142656 : Real nr_res2 = 0;
757 487984 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
758 345328 : if (act[surface])
759 203552 : 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 683804 : while (nr_res2 > 0.5 && local_iter++ < _max_iter && single_step_success)
771 531248 : 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 152556 : bool nr_good = (nr_res2 <= 0.5 && local_iter <= _max_iter && single_step_success);
787 :
788 152556 : iter += local_iter;
789 :
790 : // 'act' might have changed due to using deact_scheme = optimized, so
791 152556 : actives_tried.insert(activeCombinationNumber(act));
792 :
793 152556 : 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 592 : change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
803 592 : if (!change_scheme && deact_scheme == dumb)
804 0 : increment_dumb = canIncrementDumb(dumb_iteration);
805 :
806 592 : still_finding_solution = (change_scheme || increment_dumb);
807 :
808 592 : 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 592 : if (increment_dumb)
819 0 : incrementDumb(dumb_iteration, dumb_order, act);
820 :
821 592 : 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 151964 : kt_good = checkKuhnTucker(f, pm, act);
834 151964 : if (!kt_good)
835 : {
836 8204 : if (deact_scheme != dumb)
837 : {
838 8204 : applyKuhnTucker(f, pm, act);
839 :
840 : // true if we haven't tried this active set before
841 : still_finding_solution =
842 8204 : (actives_tried.find(activeCombinationNumber(act)) == actives_tried.end());
843 8204 : 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 8204 : 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 151948 : if (nr_good && kt_good)
882 : {
883 : // check admissible
884 : std::vector<Real> all_f;
885 143760 : 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 56160 : admissible = checkAdmissible(stress, intnl, all_f);
890 :
891 56160 : 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 2624 : bool add_constraints = canAddConstraints(act, all_f);
898 2624 : if (add_constraints)
899 : {
900 2624 : constraints_added = true;
901 2624 : std::vector<bool> act_plus(_num_surfaces,
902 2624 : false); // "act" with the positive constraints added in
903 15936 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
904 13312 : if (act[surface] ||
905 8064 : (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol)))
906 : act_plus[surface] = true;
907 2624 : if (actives_tried.find(activeCombinationNumber(act_plus)) == actives_tried.end())
908 : {
909 : // haven't tried this combination of actives yet
910 1232 : constraints_added = true;
911 7568 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
912 6336 : 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 2624 : if (!add_constraints)
922 1392 : change_scheme = canChangeScheme(deact_scheme, can_revert_to_dumb);
923 :
924 2624 : if (!add_constraints && !change_scheme && deact_scheme == dumb)
925 0 : increment_dumb = canIncrementDumb(dumb_iteration);
926 :
927 2624 : still_finding_solution = (add_constraints || change_scheme || increment_dumb);
928 :
929 2624 : if (change_scheme)
930 480 : 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 2624 : if (increment_dumb)
940 0 : incrementDumb(dumb_iteration, dumb_order, act);
941 :
942 2624 : 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 151036 : successful_return = (nr_good && admissible && kt_good);
952 : if (successful_return)
953 : break;
954 :
955 : if (still_finding_solution)
956 : {
957 9900 : stress = initial_stress;
958 9900 : delta_dp = RankTwoTensor(); // back to zero change in plastic strain
959 53468 : for (unsigned model = 0; model < _num_models; ++model)
960 43568 : intnl[model] = intnl_old[model]; // back to old internal params
961 9900 : pm.assign(_num_surfaces, 0.0); // back to zero plastic multipliers
962 :
963 9900 : unsigned num_active = numberActive(act);
964 9900 : if (num_active == 0)
965 : {
966 : successful_return = false;
967 : break; // failure
968 : }
969 :
970 9900 : 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 9900 : yieldFunction(stress, intnl, act, f);
975 :
976 9900 : nr_res2 = 0;
977 : unsigned ind = 0;
978 82252 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
979 72352 : if (act[surface])
980 : {
981 33772 : if (f[ind] > _f[modelNumber(surface)]->_f_tol)
982 33580 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
983 33772 : ind++;
984 : }
985 : }
986 : }
987 :
988 : // returned, with either success or failure
989 141744 : if (successful_return)
990 : {
991 141136 : plastic_strain += delta_dp;
992 :
993 481536 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
994 340400 : cumulative_pm[surface] += pm[surface];
995 :
996 141136 : if (final_step)
997 : consistent_tangent_operator =
998 138032 : consistentTangentOperator(stress, intnl, E_ijkl, pm, cumulative_pm);
999 :
1000 141136 : 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 45168 : yieldFunction(stress, intnl, act, f);
1006 : }
1007 : }
1008 :
1009 : return successful_return;
1010 : }
1011 :
1012 : bool
1013 2624 : ComputeMultiPlasticityStress::canAddConstraints(const std::vector<bool> & act,
1014 : const std::vector<Real> & all_f)
1015 : {
1016 5552 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1017 5552 : if (!act[surface] && (all_f[surface] > _f[modelNumber(surface)]->_f_tol))
1018 : return true;
1019 : return false;
1020 : }
1021 :
1022 : bool
1023 2000 : ComputeMultiPlasticityStress::canChangeScheme(DeactivationSchemeEnum current_deactivation_scheme,
1024 : bool can_revert_to_dumb)
1025 : {
1026 2000 : if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe)
1027 : return true;
1028 :
1029 1152 : if (current_deactivation_scheme == optimized && _deactivation_scheme == optimized_to_safe_to_dumb)
1030 : return true;
1031 :
1032 1520 : if (current_deactivation_scheme == safe && _deactivation_scheme == safe_to_dumb &&
1033 : can_revert_to_dumb)
1034 : return true;
1035 :
1036 400 : if (current_deactivation_scheme == safe && _deactivation_scheme == optimized_to_safe_to_dumb &&
1037 : can_revert_to_dumb)
1038 : return true;
1039 :
1040 1520 : 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 480 : 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 480 : if (current_deactivation_scheme == optimized &&
1058 480 : (_deactivation_scheme == optimized_to_safe ||
1059 : _deactivation_scheme == optimized_to_safe_to_dumb))
1060 : {
1061 480 : current_deactivation_scheme = safe;
1062 3808 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1063 3328 : 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 480 : }
1078 :
1079 : bool
1080 531248 : 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 531248 : Real nr_res2_before_step = nr_res2;
1098 531248 : RankTwoTensor stress_before_step;
1099 : std::vector<Real> intnl_before_step;
1100 : std::vector<Real> pm_before_step;
1101 531248 : RankTwoTensor delta_dp_before_step;
1102 :
1103 531248 : if (deactivation_scheme == optimized)
1104 : {
1105 : // we potentially use the "before_step" quantities, so record them here
1106 74476 : stress_before_step = stress;
1107 74476 : intnl_before_step.resize(_num_models);
1108 376564 : for (unsigned model = 0; model < _num_models; ++model)
1109 302088 : intnl_before_step[model] = intnl[model];
1110 74476 : pm_before_step.resize(_num_surfaces);
1111 376564 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1112 302088 : pm_before_step[surface] = pm[surface];
1113 74476 : 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 531248 : 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 531248 : 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 1075520 : while (constraints_changing)
1138 : {
1139 : // calculate dstress, dpm and dintnl for one full Newton-Raphson step
1140 544304 : nrStep(stress, intnl_old, intnl, pm, E_inv, delta_dp, dstress, dpm, dintnl, active, deact_ld);
1141 :
1142 1469908 : for (unsigned surface = 0; surface < deact_ld.size(); ++surface)
1143 946596 : if (deact_ld[surface])
1144 : {
1145 20992 : ld_encountered = true;
1146 20992 : break;
1147 : }
1148 :
1149 : // perform a line search
1150 : // The line-search will exit with updated values
1151 544304 : 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 544304 : 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 544272 : if (deactivation_scheme == optimized)
1175 : {
1176 451020 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1177 363520 : if (active[surface] && pm[surface] < 0.0)
1178 : constraints_changing = true;
1179 : }
1180 :
1181 87500 : if (constraints_changing)
1182 : {
1183 13056 : stress = stress_before_step;
1184 13056 : delta_dp = delta_dp_before_step;
1185 13056 : nr_res2 = nr_res2_before_step;
1186 74616 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1187 : {
1188 61560 : if (active[surface] && pm[surface] < 0.0)
1189 : {
1190 : // turn off the constraint forever
1191 : active[surface] = false;
1192 14464 : pm_before_step[surface] = 0.0;
1193 14464 : intnl_before_step[modelNumber(surface)] =
1194 14464 : intnl_old[modelNumber(surface)]; // don't want to muck-up hardening!
1195 : }
1196 61560 : intnl[modelNumber(surface)] = intnl_before_step[modelNumber(surface)];
1197 61560 : pm[surface] = pm_before_step[surface];
1198 : }
1199 13056 : 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 544272 : 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 531216 : if (reinstated_actives)
1214 : {
1215 : bool completely_converged = true;
1216 16112 : 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 77276 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1231 67952 : if (active[surface])
1232 43132 : if (f[ind++] > _f[modelNumber(surface)]->_f_tol)
1233 : completely_converged = false;
1234 : }
1235 : else
1236 : completely_converged = false;
1237 :
1238 9324 : if (!completely_converged)
1239 6788 : nr_res2 = residual2(pm, f, epp, ic, active, deact_ld);
1240 : }
1241 :
1242 : return successful_step;
1243 : }
1244 :
1245 : bool
1246 544272 : ComputeMultiPlasticityStress::reinstateLinearDependentConstraints(
1247 : std::vector<bool> & deactivated_due_to_ld)
1248 : {
1249 : bool reinstated_actives = false;
1250 1580768 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1251 1036496 : if (deactivated_due_to_ld[surface])
1252 : reinstated_actives = true;
1253 :
1254 544272 : deactivated_due_to_ld.assign(_num_surfaces, false);
1255 544272 : return reinstated_actives;
1256 : }
1257 :
1258 : unsigned
1259 22956 : ComputeMultiPlasticityStress::numberActive(const std::vector<bool> & active)
1260 : {
1261 : unsigned num_active = 0;
1262 156868 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1263 133912 : if (active[surface])
1264 67136 : num_active++;
1265 22956 : return num_active;
1266 : }
1267 :
1268 : bool
1269 56160 : ComputeMultiPlasticityStress::checkAdmissible(const RankTwoTensor & stress,
1270 : const std::vector<Real> & intnl,
1271 : std::vector<Real> & all_f)
1272 : {
1273 : std::vector<bool> act;
1274 56160 : act.assign(_num_surfaces, true);
1275 :
1276 56160 : yieldFunction(stress, intnl, act, all_f);
1277 :
1278 311888 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1279 258352 : if (all_f[surface] > _f[modelNumber(surface)]->_f_tol)
1280 : return false;
1281 :
1282 : return true;
1283 : }
1284 :
1285 : bool
1286 151964 : 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 557684 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1292 : {
1293 409604 : if (active[surface])
1294 : {
1295 217928 : if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1296 16624 : if (pm[surface] != 0)
1297 : return false;
1298 : }
1299 191676 : 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 522272 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1304 378512 : if (pm[surface] < 0)
1305 : return false;
1306 :
1307 : return true;
1308 : }
1309 :
1310 : void
1311 8204 : 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 70956 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1320 : {
1321 62752 : if (active[surface])
1322 : {
1323 35368 : if (f[ind++] < -_f[modelNumber(surface)]->_f_tol)
1324 14816 : if (pm[surface] != 0)
1325 : {
1326 : turned_off = true;
1327 : active[surface] = false;
1328 : }
1329 : }
1330 27384 : 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 8204 : if (!turned_off)
1337 : {
1338 : int surface_to_turn_off = -1;
1339 : Real min_pm = 0;
1340 41776 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1341 37456 : if (pm[surface] < min_pm)
1342 : {
1343 : min_pm = pm[surface];
1344 4320 : surface_to_turn_off = surface;
1345 : }
1346 4320 : if (surface_to_turn_off >= 0)
1347 : active[surface_to_turn_off] = false;
1348 : }
1349 8204 : }
1350 :
1351 : Real
1352 771384 : 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 2194704 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1363 1423320 : if (active[surface])
1364 : {
1365 1044944 : if (!deactivated_due_to_ld[surface])
1366 : {
1367 968524 : if (!(pm[surface] == 0 && f[ind] <= 0))
1368 954748 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1369 : }
1370 76420 : else if (deactivated_due_to_ld[surface] && f[ind] > 0)
1371 26308 : nr_res2 += 0.5 * Utility::pow<2>(f[ind] / _f[modelNumber(surface)]->_f_tol);
1372 1044944 : ind++;
1373 : }
1374 :
1375 771384 : nr_res2 += 0.5 * Utility::pow<2>(epp.L2norm() / _epp_tol);
1376 :
1377 771384 : std::vector<bool> active_not_deact(_num_surfaces);
1378 2194704 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1379 1878116 : active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
1380 : ind = 0;
1381 2047200 : for (unsigned model = 0; model < _num_models; ++model)
1382 1275816 : if (anyActiveSurfaces(model, active_not_deact))
1383 956012 : nr_res2 += 0.5 * Utility::pow<2>(ic[ind++] / _f[model]->_ic_tol);
1384 :
1385 771384 : return nr_res2;
1386 : }
1387 :
1388 : bool
1389 544304 : 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 544304 : 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 544304 : Real f0 = nr_res2; // initial value of residual2
1417 544304 : 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 544304 : ls_pm.resize(pm.size());
1426 :
1427 : // delta_dp during the line-search
1428 544304 : RankTwoTensor ls_delta_dp;
1429 :
1430 : // internal parameter during the line-search
1431 : std::vector<Real> ls_intnl;
1432 544304 : ls_intnl.resize(intnl.size());
1433 :
1434 : // stress during the line-search
1435 544304 : 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 2143444 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1444 1378880 : ls_pm[alpha] = pm[alpha] + dpm[alpha] * lam;
1445 764564 : ls_delta_dp = delta_dp - E_inv * dstress * lam;
1446 1995940 : for (unsigned a = 0; a < intnl.size(); ++a)
1447 1231376 : ls_intnl[a] = intnl[a] + dintnl[a] * lam;
1448 764564 : ls_stress = stress + dstress * lam;
1449 :
1450 : // calculate the new active yield functions, epp and active internal constraints
1451 764564 : 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 764564 : nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1455 :
1456 764564 : if (nr_res2 < f0 + 1E-4 * lam * slope)
1457 : break;
1458 220292 : else if (lam < lam_min)
1459 : {
1460 : success = false;
1461 : // restore plastic multipliers, yield functions, etc to original values
1462 160 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1463 128 : ls_pm[alpha] = pm[alpha];
1464 32 : ls_delta_dp = delta_dp;
1465 160 : for (unsigned a = 0; a < intnl.size(); ++a)
1466 128 : ls_intnl[a] = intnl[a];
1467 32 : ls_stress = stress;
1468 32 : calculateConstraints(
1469 : ls_stress, intnl_old, ls_intnl, ls_pm, ls_delta_dp, f, r, epp, ic, active);
1470 32 : nr_res2 = residual2(ls_pm, f, epp, ic, active, deactivated_due_to_ld);
1471 32 : break;
1472 : }
1473 220260 : else if (lam == 1.0)
1474 : {
1475 : // model as a quadratic
1476 123348 : tmp_lam = -slope / 2.0 / (nr_res2 - f0 - slope);
1477 : }
1478 : else
1479 : {
1480 : // model as a cubic
1481 96912 : Real rhs1 = nr_res2 - f0 - lam * slope;
1482 96912 : Real rhs2 = f2 - f0 - lam2 * slope;
1483 96912 : Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1484 : Real b =
1485 96912 : (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
1486 96912 : if (a == 0)
1487 0 : tmp_lam = -slope / (2.0 * b);
1488 : else
1489 : {
1490 96912 : Real disc = Utility::pow<2>(b) - 3 * a * slope;
1491 96912 : if (disc < 0)
1492 0 : tmp_lam = 0.5 * lam;
1493 96912 : else if (b <= 0)
1494 4112 : tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
1495 : else
1496 92800 : tmp_lam = -slope / (b + std::sqrt(disc));
1497 : }
1498 96912 : if (tmp_lam > 0.5 * lam)
1499 3536 : tmp_lam = 0.5 * lam;
1500 : }
1501 : lam2 = lam;
1502 : f2 = nr_res2;
1503 220260 : lam = std::max(tmp_lam, 0.1 * lam);
1504 220260 : }
1505 :
1506 544304 : if (lam < 1.0)
1507 123348 : linesearch_needed = true;
1508 :
1509 : // assign the quantities found in the line-search
1510 : // back to the originals
1511 1580928 : for (unsigned alpha = 0; alpha < pm.size(); ++alpha)
1512 1036624 : pm[alpha] = ls_pm[alpha];
1513 544304 : delta_dp = ls_delta_dp;
1514 1433424 : for (unsigned a = 0; a < intnl.size(); ++a)
1515 889120 : intnl[a] = ls_intnl[a];
1516 544304 : stress = ls_stress;
1517 :
1518 544304 : return success;
1519 : }
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 : }
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 315940 : ComputeMultiPlasticityStress::activeCombinationNumber(const std::vector<bool> & act)
1573 : {
1574 : unsigned num = 0;
1575 1227364 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1576 911424 : if (act[surface])
1577 495028 : num += (1 << surface); // (1 << x) = 2^x
1578 :
1579 315940 : return num;
1580 : }
1581 :
1582 : RankFourTensor
1583 147088 : 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 147088 : if (_tangent_operator_type == elastic)
1591 5568 : 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 141520 : std::vector<bool> act_at_some_step(_num_surfaces);
1597 517888 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1598 376368 : 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 141520 : std::vector<bool> act_vary(_num_surfaces);
1605 517888 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1606 376368 : act_vary[surface] = (pm_this_step[surface] > 0);
1607 :
1608 : std::vector<RankTwoTensor> df_dstress;
1609 141520 : dyieldFunction_dstress(stress, intnl, act_vary, df_dstress);
1610 : std::vector<Real> df_dintnl;
1611 141520 : dyieldFunction_dintnl(stress, intnl, act_vary, df_dintnl);
1612 : std::vector<RankTwoTensor> r;
1613 141520 : flowPotential(stress, intnl, act_vary, r);
1614 : std::vector<RankFourTensor> dr_dstress_at_some_step;
1615 141520 : dflowPotential_dstress(stress, intnl, act_at_some_step, dr_dstress_at_some_step);
1616 : std::vector<RankTwoTensor> dr_dintnl_at_some_step;
1617 141520 : dflowPotential_dintnl(stress, intnl, act_at_some_step, dr_dintnl_at_some_step);
1618 : std::vector<Real> h;
1619 141520 : 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 517888 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1630 376368 : if (act_vary[surface1])
1631 : {
1632 179040 : r_minus_stuff.push_back(r[ind1]);
1633 : ind2 = 0;
1634 781824 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1635 602784 : if (act_at_some_step[surface2])
1636 : {
1637 269408 : if (modelNumber(surface1) == modelNumber(surface2))
1638 : {
1639 : r_minus_stuff.back() -=
1640 224384 : cumulative_pm[surface2] * dr_dintnl_at_some_step[ind2] * h[ind1];
1641 : }
1642 269408 : ind2++;
1643 : }
1644 179040 : ind1++;
1645 : }
1646 :
1647 : unsigned int num_currently_active = 0;
1648 517888 : for (unsigned surface = 0; surface < _num_surfaces; ++surface)
1649 376368 : if (act_vary[surface])
1650 179040 : 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 141520 : zzz.assign(num_currently_active * num_currently_active, 0.0);
1660 :
1661 : ind1 = 0;
1662 141520 : RankTwoTensor r2;
1663 517888 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1664 376368 : if (act_vary[surface1])
1665 : {
1666 : ind2 = 0;
1667 781824 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1668 602784 : if (act_vary[surface2])
1669 : {
1670 269408 : r2 = df_dstress[ind1] * (E_ijkl * r_minus_stuff[ind2]);
1671 269408 : zzz[ind1 * num_currently_active + ind2] += r2(0, 0) + r2(1, 1) + r2(2, 2);
1672 269408 : if (modelNumber(surface1) == modelNumber(surface2))
1673 224384 : zzz[ind1 * num_currently_active + ind2] += df_dintnl[ind1] * h[ind2];
1674 269408 : ind2++;
1675 : }
1676 179040 : ind1++;
1677 : }
1678 :
1679 141520 : 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 141520 : 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 141520 : RankFourTensor strain_coeff = E_ijkl;
1694 : ind1 = 0;
1695 517888 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1696 376368 : if (act_vary[surface1])
1697 : {
1698 179040 : RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1699 : ind2 = 0;
1700 781824 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1701 602784 : if (act_vary[surface2])
1702 : {
1703 269408 : RankTwoTensor part2 = E_ijkl * df_dstress[ind2];
1704 1077632 : for (unsigned i = 0; i < 3; i++)
1705 3232896 : for (unsigned j = 0; j < 3; j++)
1706 9698688 : for (unsigned k = 0; k < 3; k++)
1707 29096064 : for (unsigned l = 0; l < 3; l++)
1708 21822048 : strain_coeff(i, j, k, l) -=
1709 21822048 : part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1710 269408 : ind2++;
1711 : }
1712 179040 : ind1++;
1713 : }
1714 :
1715 141520 : if (_tangent_operator_type == linear)
1716 9024 : return strain_coeff;
1717 :
1718 132496 : RankFourTensor stress_coeff(RankFourTensor::initIdentitySymmetricFour);
1719 :
1720 132496 : RankFourTensor part3;
1721 : ind1 = 0;
1722 478336 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1723 345840 : if (act_at_some_step[surface1])
1724 : {
1725 332864 : part3 += cumulative_pm[surface1] * E_ijkl * dr_dstress_at_some_step[ind1];
1726 166432 : ind1++;
1727 : }
1728 :
1729 132496 : stress_coeff += part3;
1730 :
1731 132496 : 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 478336 : for (unsigned surface1 = 0; surface1 < _num_surfaces; ++surface1)
1736 345840 : if (act_vary[surface1])
1737 : {
1738 166432 : RankTwoTensor part1 = E_ijkl * r_minus_stuff[ind1];
1739 : ind2 = 0;
1740 720768 : for (unsigned surface2 = 0; surface2 < _num_surfaces; ++surface2)
1741 554336 : if (act_vary[surface2])
1742 : {
1743 249632 : RankTwoTensor part2 = part3 * df_dstress[ind2];
1744 998528 : for (unsigned i = 0; i < 3; i++)
1745 2995584 : for (unsigned j = 0; j < 3; j++)
1746 8986752 : for (unsigned k = 0; k < 3; k++)
1747 26960256 : for (unsigned l = 0; l < 3; l++)
1748 20220192 : stress_coeff(i, j, k, l) -=
1749 20220192 : part1(i, j) * part2(k, l) * zzz[ind1 * num_currently_active + ind2];
1750 249632 : ind2++;
1751 : }
1752 166432 : 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 132496 : RankFourTensor s_inv;
1764 : try
1765 : {
1766 132496 : s_inv = stress_coeff.invSymm();
1767 : }
1768 212 : catch (const MooseException & e)
1769 : {
1770 212 : return strain_coeff; // when stress_coeff is singular (perhaps for incompressible plasticity?)
1771 : // return the "linear" tangent operator
1772 212 : }
1773 :
1774 132284 : return s_inv * strain_coeff;
1775 : }
|
