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 "MultiParameterPlasticityStressUpdate.h"
11 : #include "Conversion.h" // for stringify
12 : #include "MooseEnum.h" // for enum
13 :
14 : // libMesh includes
15 : #include "libmesh/utility.h" // for Utility::pow
16 :
17 : // PETSc includes
18 : #include <petscblaslapack.h> // LAPACKgesv_
19 :
20 : InputParameters
21 3724 : MultiParameterPlasticityStressUpdate::validParams()
22 : {
23 3724 : InputParameters params = StressUpdateBase::validParams();
24 3724 : params.addClassDescription("Return-map and Jacobian algorithms for plastic models where the "
25 : "yield function and flow directions depend on multiple functions of "
26 : "stress");
27 11172 : params.addRangeCheckedParam<unsigned int>(
28 : "max_NR_iterations",
29 7448 : 20,
30 : "max_NR_iterations>0",
31 : "Maximum number of Newton-Raphson iterations allowed during the return-map algorithm");
32 7448 : params.addParam<bool>("perform_finite_strain_rotations",
33 7448 : false,
34 : "Tensors are correctly rotated "
35 : "in finite-strain simulations. "
36 : "For optimal performance you can "
37 : "set this to 'false' if you are "
38 : "only ever using small strains");
39 7448 : params.addRequiredParam<Real>(
40 : "smoothing_tol",
41 : "Intersections of the yield surfaces will be smoothed by this amount (this "
42 : "is measured in units of stress). Often this is related to other physical "
43 : "parameters (eg, 0.1*cohesion) but it is important to set this small enough "
44 : "so that the individual yield surfaces do not mix together in the smoothing "
45 : "process to produce a result where no stress is admissible (for example, "
46 : "mixing together tensile and compressive failure envelopes).");
47 7448 : params.addRequiredParam<Real>("yield_function_tol",
48 : "The return-map process will be deemed to have converged if all "
49 : "yield functions are within yield_function_tol of zero. If this "
50 : "is set very low then precision-loss might be encountered: if the "
51 : "code detects precision loss then it also deems the return-map "
52 : "process has converged.");
53 7448 : MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
54 7448 : params.addParam<Real>("min_step_size",
55 7448 : 1.0,
56 : "In order to help the Newton-Raphson procedure, the applied strain "
57 : "increment may be applied in sub-increments of size greater than this "
58 : "value. Usually it is better for Moose's nonlinear convergence to "
59 : "increase max_NR_iterations rather than decrease this parameter.");
60 7448 : params.addParam<bool>("warn_about_precision_loss",
61 7448 : false,
62 : "Output a message to the console every "
63 : "time precision-loss is encountered "
64 : "during the Newton-Raphson process");
65 7448 : params.addParam<std::vector<Real>>("admissible_stress",
66 : "A single admissible value of the value of the stress "
67 : "parameters for internal parameters = 0. This is used "
68 : "to initialize the return-mapping algorithm during the first "
69 : "nonlinear iteration. If not given then it is assumed that "
70 : "stress parameters = 0 is admissible.");
71 7448 : MooseEnum smoother_fcn_enum("cos poly1 poly2 poly3", "cos");
72 7448 : params.addParam<MooseEnum>("smoother_function_type",
73 : smoother_fcn_enum,
74 : "Type of smoother function to use. 'cos' means (-a/pi)cos(pi x/2/a), "
75 : "'polyN' means a polynomial of degree 2N+2");
76 7448 : params.addParamNamesToGroup("smoother_function_type", "Advanced");
77 3724 : return params;
78 3724 : }
79 :
80 2796 : MultiParameterPlasticityStressUpdate::MultiParameterPlasticityStressUpdate(
81 2796 : const InputParameters & parameters, unsigned num_sp, unsigned num_yf, unsigned num_intnl)
82 : : StressUpdateBase(parameters),
83 2796 : _num_sp(num_sp),
84 2796 : _definitely_ok_sp(isParamValid("admissible_stress")
85 2796 : ? getParam<std::vector<Real>>("admissible_stress")
86 2796 : : std::vector<Real>(_num_sp, 0.0)),
87 2796 : _Eij(num_sp, std::vector<Real>(num_sp)),
88 2796 : _En(1.0),
89 2796 : _Cij(num_sp, std::vector<Real>(num_sp)),
90 2796 : _num_yf(num_yf),
91 2796 : _num_intnl(num_intnl),
92 5592 : _max_nr_its(getParam<unsigned>("max_NR_iterations")),
93 5592 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
94 5592 : _smoothing_tol(getParam<Real>("smoothing_tol")),
95 5592 : _smoothing_tol2(Utility::pow<2>(getParam<Real>("smoothing_tol"))),
96 5592 : _f_tol(getParam<Real>("yield_function_tol")),
97 5592 : _f_tol2(Utility::pow<2>(getParam<Real>("yield_function_tol"))),
98 5592 : _min_step_size(getParam<Real>("min_step_size")),
99 5592 : _step_one(declareRestartableData<bool>("step_one", true)),
100 5592 : _warn_about_precision_loss(getParam<bool>("warn_about_precision_loss")),
101 :
102 2796 : _plastic_strain(declareProperty<RankTwoTensor>(_base_name + "plastic_strain")),
103 5592 : _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "plastic_strain")),
104 2796 : _intnl(declareProperty<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
105 2796 : _intnl_old(
106 2796 : getMaterialPropertyOld<std::vector<Real>>(_base_name + "plastic_internal_parameter")),
107 2796 : _yf(declareProperty<std::vector<Real>>(_base_name + "plastic_yield_function")),
108 2796 : _iter(declareProperty<Real>(_base_name +
109 : "plastic_NR_iterations")), // this is really an unsigned int, but
110 : // for visualisation i convert it to Real
111 2796 : _max_iter_used(declareProperty<Real>(
112 2796 : _base_name + "max_plastic_NR_iterations")), // this is really an unsigned int, but
113 : // for visualisation i convert it to Real
114 5592 : _max_iter_used_old(getMaterialPropertyOld<Real>(_base_name + "max_plastic_NR_iterations")),
115 2796 : _linesearch_needed(
116 2796 : declareProperty<Real>(_base_name + "plastic_linesearch_needed")), // this is really a
117 : // boolean, but for
118 : // visualisation i
119 : // convert it to Real
120 2796 : _trial_sp(num_sp),
121 2796 : _stress_trial(RankTwoTensor()),
122 2796 : _rhs(num_sp + 1),
123 2796 : _dvar_dtrial(num_sp + 1, std::vector<Real>(num_sp, 0.0)),
124 2796 : _ok_sp(num_sp),
125 2796 : _ok_intnl(num_intnl),
126 2796 : _del_stress_params(num_sp),
127 2796 : _current_sp(num_sp),
128 2796 : _current_intnl(num_intnl),
129 2796 : _smoother_function_type(
130 5592 : parameters.get<MooseEnum>("smoother_function_type").getEnum<SmootherFunctionType>())
131 : {
132 2796 : if (_definitely_ok_sp.size() != _num_sp)
133 0 : mooseError("MultiParameterPlasticityStressUpdate: admissible_stress parameter must consist of ",
134 0 : _num_sp,
135 : " numbers");
136 2796 : }
137 :
138 : void
139 204928 : MultiParameterPlasticityStressUpdate::initQpStatefulProperties()
140 : {
141 204928 : _plastic_strain[_qp].zero();
142 204928 : _intnl[_qp].assign(_num_intnl, 0);
143 204928 : _yf[_qp].assign(_num_yf, 0);
144 204928 : _iter[_qp] = 0.0;
145 204928 : _max_iter_used[_qp] = 0.0;
146 204928 : _linesearch_needed[_qp] = 0.0;
147 204928 : }
148 :
149 : void
150 192 : MultiParameterPlasticityStressUpdate::propagateQpStatefulProperties()
151 : {
152 192 : _plastic_strain[_qp] = _plastic_strain_old[_qp];
153 192 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
154 192 : _max_iter_used[_qp] = _max_iter_used_old[_qp];
155 192 : }
156 :
157 : void
158 1428532 : MultiParameterPlasticityStressUpdate::updateState(RankTwoTensor & strain_increment,
159 : RankTwoTensor & inelastic_strain_increment,
160 : const RankTwoTensor & rotation_increment,
161 : RankTwoTensor & stress_new,
162 : const RankTwoTensor & stress_old,
163 : const RankFourTensor & elasticity_tensor,
164 : const RankTwoTensor & /*elastic_strain_old*/,
165 : bool compute_full_tangent_operator,
166 : RankFourTensor & tangent_operator)
167 : {
168 : // Size _yf[_qp] appropriately
169 1428532 : _yf[_qp].assign(_num_yf, 0);
170 : // _plastic_strain and _intnl are usually sized appropriately because they are stateful, but this
171 : // Material may be used from a DiracKernel where stateful materials are not allowed. The best we
172 : // can do is:
173 1428532 : if (_intnl[_qp].size() != _num_intnl)
174 0 : initQpStatefulProperties();
175 :
176 1428532 : initializeReturnProcess();
177 :
178 1428532 : if (_t_step >= 2)
179 1051904 : _step_one = false;
180 :
181 : // initially assume an elastic deformation
182 1428532 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _intnl[_qp].begin());
183 :
184 1428532 : _iter[_qp] = 0.0;
185 1428532 : _max_iter_used[_qp] = std::max(_max_iter_used[_qp], _max_iter_used_old[_qp]);
186 1428532 : _linesearch_needed[_qp] = 0.0;
187 :
188 1428532 : computeStressParams(stress_new, _trial_sp);
189 1428532 : yieldFunctionValuesV(_trial_sp, _intnl[_qp], _yf[_qp]);
190 :
191 1428532 : if (yieldF(_yf[_qp]) <= _f_tol)
192 : {
193 802880 : _plastic_strain[_qp] = _plastic_strain_old[_qp];
194 : inelastic_strain_increment.zero();
195 802880 : if (_fe_problem.currentlyComputingJacobian())
196 198128 : tangent_operator = elasticity_tensor;
197 802880 : return;
198 : }
199 :
200 625652 : _stress_trial = stress_new;
201 : /* The trial stress must be inadmissible
202 : * so we need to return to the yield surface. The following
203 : * equations must be satisfied.
204 : *
205 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
206 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
207 : * ...
208 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
209 : * 0 = rhs[N] = f(S, intnl)
210 : *
211 : * as well as equations defining intnl parameters as functions of
212 : * stress_params, trial_stress_params and intnl_old
213 : *
214 : * The unknowns are S[0], ..., S[N-1], gaE, and the intnl parameters.
215 : * Here gaE = ga * _En (the _En serves to make gaE similar magnitude to S)
216 : * I find it convenient to solve the first N+1 equations for p, q and gaE,
217 : * while substituting the "intnl parameters" equations into these during the solve process
218 : */
219 :
220 2620864 : for (auto & deriv : _dvar_dtrial)
221 1995212 : deriv.assign(_num_sp, 0.0);
222 :
223 625652 : preReturnMapV(_trial_sp, stress_new, _intnl_old[_qp], _yf[_qp], elasticity_tensor);
224 :
225 625652 : setEffectiveElasticity(elasticity_tensor);
226 :
227 625652 : if (_step_one)
228 265828 : std::copy(_definitely_ok_sp.begin(), _definitely_ok_sp.end(), _ok_sp.begin());
229 : else
230 359824 : computeStressParams(stress_old, _ok_sp);
231 625652 : std::copy(_intnl_old[_qp].begin(), _intnl_old[_qp].end(), _ok_intnl.begin());
232 :
233 : // Return-map problem: must apply the following changes in stress_params,
234 : // and find the returned stress_params and gaE
235 1995212 : for (unsigned i = 0; i < _num_sp; ++i)
236 1369560 : _del_stress_params[i] = _trial_sp[i] - _ok_sp[i];
237 :
238 : Real step_taken = 0.0; // amount of del_stress_params that we've applied and the return-map
239 : // problem has succeeded
240 : Real step_size = 1.0; // potentially can apply del_stress_params in substeps
241 : Real gaE_total = 0.0;
242 :
243 : // current values of the yield function, derivatives, etc
244 : yieldAndFlow smoothed_q;
245 :
246 : // In the following sub-stepping procedure it is possible that
247 : // the last step is an elastic step, and therefore smoothed_q won't
248 : // be computed on the last step, so we have to compute it.
249 : bool smoothed_q_calculated = false;
250 :
251 1251304 : while (step_taken < 1.0 && step_size >= _min_step_size)
252 : {
253 625652 : if (1.0 - step_taken < step_size)
254 : // prevent over-shoots of substepping
255 : step_size = 1.0 - step_taken;
256 :
257 : // trial variables in terms of admissible variables
258 1995212 : for (unsigned i = 0; i < _num_sp; ++i)
259 1369560 : _trial_sp[i] = _ok_sp[i] + step_size * _del_stress_params[i];
260 :
261 : // initialize variables that are to be found via Newton-Raphson
262 625652 : _current_sp = _trial_sp;
263 625652 : Real gaE = 0.0;
264 :
265 : // flags indicating failure of Newton-Raphson and line-search
266 : int nr_failure = 0;
267 : int ls_failure = 0;
268 :
269 : // NR iterations taken in this substep
270 : unsigned step_iter = 0;
271 :
272 : // The residual-squared for the line-search
273 625652 : Real res2 = 0.0;
274 :
275 625652 : if (step_size < 1.0 && yieldF(_trial_sp, _ok_intnl) <= _f_tol)
276 : // This is an elastic step
277 : // The "step_size < 1.0" in above condition is for efficiency: we definitely
278 : // know that this is a plastic step if step_size = 1.0
279 : smoothed_q_calculated = false;
280 : else
281 : {
282 : // this is a plastic step
283 :
284 : // initialize current_sp, gaE and current_intnl based on the non-smoothed situation
285 625652 : initializeVarsV(_trial_sp, _ok_intnl, _current_sp, gaE, _current_intnl);
286 : // and find the smoothed yield function, flow potential and derivatives
287 625652 : smoothed_q = smoothAllQuantities(_current_sp, _current_intnl);
288 : smoothed_q_calculated = true;
289 625652 : calculateRHS(_trial_sp, _current_sp, gaE, smoothed_q, _rhs);
290 625652 : res2 = calculateRes2(_rhs);
291 :
292 : // Perform a Newton-Raphson with linesearch to get current_sp, gaE, and also smoothed_q
293 2232986 : while (res2 > _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0)
294 : {
295 : // solve the linear system and store the answer (the "updates") in rhs
296 1607914 : nr_failure = nrStep(smoothed_q, _trial_sp, _current_sp, _current_intnl, gaE, _rhs);
297 1607914 : if (nr_failure != 0)
298 : break;
299 :
300 : // handle precision loss
301 1607914 : if (precisionLoss(_rhs, _current_sp, gaE))
302 : {
303 580 : if (_warn_about_precision_loss)
304 : {
305 : Moose::err << "MultiParameterPlasticityStressUpdate: precision-loss in element "
306 0 : << _current_elem->id() << " quadpoint=" << _qp << ". Stress_params =";
307 0 : for (unsigned i = 0; i < _num_sp; ++i)
308 0 : Moose::err << " " << _current_sp[i];
309 : Moose::err << " gaE = " << gaE << "\n";
310 : }
311 580 : res2 = 0.0;
312 580 : break;
313 : }
314 :
315 : // apply (parts of) the updates, re-calculate smoothed_q, and res2
316 1607334 : ls_failure = lineSearch(res2,
317 : _current_sp,
318 : gaE,
319 : _trial_sp,
320 : smoothed_q,
321 1607334 : _ok_intnl,
322 : _current_intnl,
323 : _rhs,
324 1607334 : _linesearch_needed[_qp]);
325 1607334 : step_iter++;
326 : }
327 : }
328 625652 : if (res2 <= _f_tol2 && step_iter < _max_nr_its && nr_failure == 0 && ls_failure == 0 &&
329 625568 : gaE >= 0.0)
330 : {
331 : // this Newton-Raphson worked fine, or this was an elastic step
332 625568 : std::copy(_current_sp.begin(), _current_sp.end(), _ok_sp.begin());
333 625568 : gaE_total += gaE;
334 625568 : step_taken += step_size;
335 625568 : setIntnlValuesV(_trial_sp, _ok_sp, _ok_intnl, _intnl[_qp]);
336 625568 : std::copy(_intnl[_qp].begin(), _intnl[_qp].end(), _ok_intnl.begin());
337 : // calculate dp/dp_trial, dp/dq_trial, etc, for Jacobian
338 625568 : dVardTrial(!smoothed_q_calculated,
339 : _trial_sp,
340 : _ok_sp,
341 : gaE,
342 : _ok_intnl,
343 : smoothed_q,
344 : step_size,
345 : compute_full_tangent_operator,
346 625568 : _dvar_dtrial);
347 625568 : if (static_cast<Real>(step_iter) > _iter[_qp])
348 436880 : _iter[_qp] = static_cast<Real>(step_iter);
349 625568 : if (static_cast<Real>(step_iter) > _max_iter_used[_qp])
350 22096 : _max_iter_used[_qp] = static_cast<Real>(step_iter);
351 625568 : step_size *= 1.1;
352 : }
353 : else
354 : {
355 : // Newton-Raphson + line-search process failed
356 84 : step_size *= 0.5;
357 : }
358 : }
359 :
360 625652 : if (step_size < _min_step_size)
361 84 : errorHandler("MultiParameterPlasticityStressUpdate: Minimum step-size violated");
362 :
363 : // success!
364 625568 : finalizeReturnProcess(rotation_increment);
365 625568 : yieldFunctionValuesV(_ok_sp, _intnl[_qp], _yf[_qp]);
366 :
367 625568 : if (!smoothed_q_calculated)
368 0 : smoothed_q = smoothAllQuantities(_ok_sp, _intnl[_qp]);
369 :
370 625568 : setStressAfterReturnV(
371 625568 : _stress_trial, _ok_sp, gaE_total, _intnl[_qp], smoothed_q, elasticity_tensor, stress_new);
372 :
373 625568 : setInelasticStrainIncrementAfterReturn(_stress_trial,
374 : gaE_total,
375 : smoothed_q,
376 : elasticity_tensor,
377 : stress_new,
378 : inelastic_strain_increment);
379 :
380 625568 : strain_increment = strain_increment - inelastic_strain_increment;
381 625568 : _plastic_strain[_qp] = _plastic_strain_old[_qp] + inelastic_strain_increment;
382 :
383 625568 : if (_fe_problem.currentlyComputingJacobian())
384 : // for efficiency, do not compute the tangent operator if not currently computing Jacobian
385 97728 : consistentTangentOperatorV(_stress_trial,
386 : _trial_sp,
387 : stress_new,
388 : _ok_sp,
389 : gaE_total,
390 : smoothed_q,
391 : elasticity_tensor,
392 : compute_full_tangent_operator,
393 97728 : _dvar_dtrial,
394 : tangent_operator);
395 625652 : }
396 :
397 : MultiParameterPlasticityStressUpdate::yieldAndFlow
398 2509806 : MultiParameterPlasticityStressUpdate::smoothAllQuantities(const std::vector<Real> & stress_params,
399 : const std::vector<Real> & intnl) const
400 : {
401 2509806 : std::vector<yieldAndFlow> all_q(_num_yf, yieldAndFlow(_num_sp, _num_intnl));
402 2509806 : computeAllQV(stress_params, intnl, all_q);
403 :
404 : /* This routine holds the key to my smoothing strategy. It
405 : * may be proved that this smoothing strategy produces a
406 : * yield surface that is both C2 differentiable and convex,
407 : * assuming the individual yield functions are C2 and
408 : * convex too.
409 : * Of course all the derivatives must also be smoothed.
410 : * Also, I assume that d(flow potential)/dstress gets smoothed
411 : * by the Yield Function (which produces a C2 flow potential).
412 : * See the line identified in the loop below.
413 : * Only time will tell whether this is a good strategy, but it
414 : * works well in all tests so far. Convexity is irrelevant
415 : * for the non-associated case, but at least the return-map
416 : * problem should always have a unique solution.
417 : * For two yield functions+flows, labelled 1 and 2, we
418 : * should have
419 : * d(g1 - g2) . d(f1 - f2) >= 0
420 : * If not then the return-map problem for even the
421 : * multi-surface plasticity with no smoothing won't have a
422 : * unique solution. If the multi-surface plasticity has
423 : * a unique solution then the smoothed version defined
424 : * below will too.
425 : */
426 :
427 : // res_f is the index that contains the smoothed yieldAndFlow
428 : std::size_t res_f = 0;
429 :
430 14917338 : for (std::size_t a = 1; a < all_q.size(); ++a)
431 : {
432 12407532 : if (all_q[res_f].f >= all_q[a].f + _smoothing_tol)
433 : // no smoothing is needed: res_f is already indexes the largest yield function
434 9514184 : continue;
435 2893348 : else if (all_q[a].f >= all_q[res_f].f + _smoothing_tol)
436 : {
437 : // no smoothing is needed, and res_f needs to index to all_q[a]
438 : res_f = a;
439 986340 : continue;
440 : }
441 : else
442 : {
443 : // smoothing is required
444 1907008 : const Real f_diff = all_q[res_f].f - all_q[a].f;
445 1907008 : const Real ism = ismoother(f_diff);
446 1907008 : const Real sm = smoother(f_diff);
447 1907008 : const Real dsm = dsmoother(f_diff);
448 : // we want: all_q[res_f].f = 0.5 * all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism,
449 : // but we have to do the derivatives first
450 6829312 : for (unsigned i = 0; i < _num_sp; ++i)
451 : {
452 18091776 : for (unsigned j = 0; j < _num_sp; ++j)
453 13169472 : all_q[res_f].d2g[i][j] =
454 13169472 : 0.5 * (all_q[res_f].d2g[i][j] + all_q[a].d2g[i][j]) +
455 13169472 : dsm * (all_q[res_f].df[j] - all_q[a].df[j]) * (all_q[res_f].dg[i] - all_q[a].dg[i]) +
456 13169472 : sm * (all_q[res_f].d2g[i][j] - all_q[a].d2g[i][j]);
457 14525760 : for (unsigned j = 0; j < _num_intnl; ++j)
458 9603456 : all_q[res_f].d2g_di[i][j] = 0.5 * (all_q[res_f].d2g_di[i][j] + all_q[a].d2g_di[i][j]) +
459 9603456 : dsm * (all_q[res_f].df_di[j] - all_q[a].df_di[j]) *
460 9603456 : (all_q[res_f].dg[i] - all_q[a].dg[i]) +
461 9603456 : sm * (all_q[res_f].d2g_di[i][j] - all_q[a].d2g_di[i][j]);
462 : }
463 6829312 : for (unsigned i = 0; i < _num_sp; ++i)
464 : {
465 4922304 : all_q[res_f].df[i] = 0.5 * (all_q[res_f].df[i] + all_q[a].df[i]) +
466 4922304 : sm * (all_q[res_f].df[i] - all_q[a].df[i]);
467 : // whether the following (smoothing g with f's smoother) is a good strategy remains to be
468 : // seen...
469 4922304 : all_q[res_f].dg[i] = 0.5 * (all_q[res_f].dg[i] + all_q[a].dg[i]) +
470 4922304 : sm * (all_q[res_f].dg[i] - all_q[a].dg[i]);
471 : }
472 5640640 : for (unsigned i = 0; i < _num_intnl; ++i)
473 3733632 : all_q[res_f].df_di[i] = 0.5 * (all_q[res_f].df_di[i] + all_q[a].df_di[i]) +
474 3733632 : sm * (all_q[res_f].df_di[i] - all_q[a].df_di[i]);
475 1907008 : all_q[res_f].f = 0.5 * (all_q[res_f].f + all_q[a].f + _smoothing_tol) + ism;
476 : }
477 : }
478 5019612 : return all_q[res_f];
479 2509806 : }
480 :
481 : Real
482 2190256 : MultiParameterPlasticityStressUpdate::ismoother(Real f_diff) const
483 : {
484 2190256 : if (std::abs(f_diff) >= _smoothing_tol)
485 : return 0.0;
486 2190256 : switch (_smoother_function_type)
487 : {
488 2190256 : case SmootherFunctionType::cos:
489 2190256 : return -_smoothing_tol / M_PI * std::cos(0.5 * M_PI * f_diff / _smoothing_tol);
490 0 : case SmootherFunctionType::poly1:
491 0 : return 0.75 / _smoothing_tol *
492 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
493 0 : (_smoothing_tol2 / 12.0) * (Utility::pow<4>(f_diff / _smoothing_tol) - 1.0));
494 0 : case SmootherFunctionType::poly2:
495 0 : return 0.625 / _smoothing_tol *
496 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
497 0 : (_smoothing_tol2 / 30.0) * (Utility::pow<6>(f_diff / _smoothing_tol) - 1.0));
498 0 : case SmootherFunctionType::poly3:
499 0 : return (7.0 / 12.0 / _smoothing_tol) *
500 0 : (0.5 * (Utility::pow<2>(f_diff) - _smoothing_tol2) -
501 0 : (_smoothing_tol2 / 56.0) * (Utility::pow<8>(f_diff / _smoothing_tol) - 1.0));
502 : default:
503 : return 0.0;
504 : }
505 : }
506 :
507 : Real
508 1907008 : MultiParameterPlasticityStressUpdate::smoother(Real f_diff) const
509 : {
510 1907008 : if (std::abs(f_diff) >= _smoothing_tol)
511 : return 0.0;
512 1907008 : switch (_smoother_function_type)
513 : {
514 1907008 : case SmootherFunctionType::cos:
515 1907008 : return 0.5 * std::sin(f_diff * M_PI * 0.5 / _smoothing_tol);
516 0 : case SmootherFunctionType::poly1:
517 0 : return 0.75 / _smoothing_tol *
518 0 : (f_diff - (_smoothing_tol / 3.0) * Utility::pow<3>(f_diff / _smoothing_tol));
519 0 : case SmootherFunctionType::poly2:
520 0 : return 0.625 / _smoothing_tol *
521 0 : (f_diff - (_smoothing_tol / 5.0) * Utility::pow<5>(f_diff / _smoothing_tol));
522 0 : case SmootherFunctionType::poly3:
523 0 : return (7.0 / 12.0 / _smoothing_tol) *
524 0 : (f_diff - (_smoothing_tol / 7.0) * Utility::pow<7>(f_diff / _smoothing_tol));
525 : default:
526 : return 0.0;
527 : }
528 : }
529 :
530 : Real
531 1907008 : MultiParameterPlasticityStressUpdate::dsmoother(Real f_diff) const
532 : {
533 1907008 : if (std::abs(f_diff) >= _smoothing_tol)
534 : return 0.0;
535 1907008 : switch (_smoother_function_type)
536 : {
537 1907008 : case SmootherFunctionType::cos:
538 1907008 : return 0.25 * M_PI / _smoothing_tol * std::cos(f_diff * M_PI * 0.5 / _smoothing_tol);
539 0 : case SmootherFunctionType::poly1:
540 0 : return 0.75 / _smoothing_tol * (1.0 - Utility::pow<2>(f_diff / _smoothing_tol));
541 0 : case SmootherFunctionType::poly2:
542 0 : return 0.625 / _smoothing_tol * (1.0 - Utility::pow<4>(f_diff / _smoothing_tol));
543 0 : case SmootherFunctionType::poly3:
544 0 : return (7.0 / 12.0 / _smoothing_tol) * (1.0 - Utility::pow<6>(f_diff / _smoothing_tol));
545 : default:
546 : return 0.0;
547 : }
548 : }
549 :
550 : int
551 1607334 : MultiParameterPlasticityStressUpdate::lineSearch(Real & res2,
552 : std::vector<Real> & stress_params,
553 : Real & gaE,
554 : const std::vector<Real> & trial_stress_params,
555 : yieldAndFlow & smoothed_q,
556 : const std::vector<Real> & intnl_ok,
557 : std::vector<Real> & intnl,
558 : std::vector<Real> & rhs,
559 : Real & linesearch_needed) const
560 : {
561 1607334 : const Real res2_old = res2;
562 1607334 : const std::vector<Real> sp_params_old = stress_params;
563 1607334 : const Real gaE_old = gaE;
564 1607334 : const std::vector<Real> delta_nr_params = rhs;
565 :
566 : Real lam = 1.0; // line-search parameter
567 : const Real lam_min = 1E-10; // minimum value of lam allowed
568 1607334 : const Real slope = -2.0 * res2_old; // "Numerical Recipes" uses -b*A*x, in order to check for
569 : // roundoff, but i hope the nrStep would warn if there were
570 : // problems
571 : Real tmp_lam; // cached value of lam used in quadratic & cubic line search
572 : Real f2 = res2_old; // cached value of f = residual2 used in the cubic in the line search
573 : Real lam2 = lam; // cached value of lam used in the cubic in the line search
574 :
575 : while (true)
576 : {
577 : // update variables using the current line-search parameter
578 6407150 : for (unsigned i = 0; i < _num_sp; ++i)
579 4522996 : stress_params[i] = sp_params_old[i] + lam * delta_nr_params[i];
580 1884154 : gaE = gaE_old + lam * delta_nr_params[_num_sp];
581 :
582 : // and internal parameters
583 1884154 : setIntnlValuesV(trial_stress_params, stress_params, intnl_ok, intnl);
584 :
585 1884154 : smoothed_q = smoothAllQuantities(stress_params, intnl);
586 :
587 : // update rhs for next-time through
588 1884154 : calculateRHS(trial_stress_params, stress_params, gaE, smoothed_q, rhs);
589 1884154 : res2 = calculateRes2(rhs);
590 :
591 : // do the line-search
592 1884154 : if (res2 < res2_old + 1E-4 * lam * slope)
593 : break;
594 276902 : else if (lam < lam_min)
595 : return 1;
596 276820 : else if (lam == 1.0)
597 : {
598 : // model as a quadratic
599 169394 : tmp_lam = -0.5 * slope / (res2 - res2_old - slope);
600 : }
601 : else
602 : {
603 : // model as a cubic
604 107426 : const Real rhs1 = res2 - res2_old - lam * slope;
605 107426 : const Real rhs2 = f2 - res2_old - lam2 * slope;
606 107426 : const Real a = (rhs1 / Utility::pow<2>(lam) - rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
607 : const Real b =
608 107426 : (-lam2 * rhs1 / Utility::pow<2>(lam) + lam * rhs2 / Utility::pow<2>(lam2)) / (lam - lam2);
609 107426 : if (a == 0.0)
610 0 : tmp_lam = -slope / (2.0 * b);
611 : else
612 : {
613 107426 : const Real disc = Utility::pow<2>(b) - 3.0 * a * slope;
614 107426 : if (disc < 0)
615 0 : tmp_lam = 0.5 * lam;
616 107426 : else if (b <= 0)
617 3840 : tmp_lam = (-b + std::sqrt(disc)) / (3.0 * a);
618 : else
619 103586 : tmp_lam = -slope / (b + std::sqrt(disc));
620 : }
621 107426 : if (tmp_lam > 0.5 * lam)
622 2352 : tmp_lam = 0.5 * lam;
623 : }
624 : lam2 = lam;
625 : f2 = res2;
626 276820 : lam = std::max(tmp_lam, 0.1 * lam);
627 276820 : }
628 :
629 1607252 : if (lam < 1.0)
630 169312 : linesearch_needed = 1.0;
631 : return 0;
632 : }
633 :
634 : int
635 1607914 : MultiParameterPlasticityStressUpdate::nrStep(const yieldAndFlow & smoothed_q,
636 : const std::vector<Real> & trial_stress_params,
637 : const std::vector<Real> & stress_params,
638 : const std::vector<Real> & intnl,
639 : Real gaE,
640 : std::vector<Real> & rhs) const
641 : {
642 1607914 : std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
643 1607914 : setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
644 :
645 1607914 : std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
646 1607914 : dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
647 :
648 : // use LAPACK to solve the linear system
649 1607914 : const PetscBLASInt nrhs = 1;
650 1607914 : std::vector<PetscBLASInt> ipiv(_num_sp + 1);
651 : PetscBLASInt info;
652 1607914 : const PetscBLASInt gesv_num_rhs = _num_sp + 1;
653 1607914 : LAPACKgesv_(
654 : &gesv_num_rhs, &nrhs, &jac[0], &gesv_num_rhs, &ipiv[0], &rhs[0], &gesv_num_rhs, &info);
655 3215828 : return info;
656 1607914 : }
657 :
658 : void
659 84 : MultiParameterPlasticityStressUpdate::errorHandler(const std::string & message) const
660 : {
661 168 : throw MooseException(message);
662 : }
663 :
664 : void
665 131712 : MultiParameterPlasticityStressUpdate::initializeReturnProcess()
666 : {
667 131712 : }
668 :
669 : void
670 118256 : MultiParameterPlasticityStressUpdate::finalizeReturnProcess(
671 : const RankTwoTensor & /*rotation_increment*/)
672 : {
673 118256 : }
674 :
675 : void
676 0 : MultiParameterPlasticityStressUpdate::preReturnMapV(
677 : const std::vector<Real> & /*trial_stress_params*/,
678 : const RankTwoTensor & /*stress_trial*/,
679 : const std::vector<Real> & /*intnl_old*/,
680 : const std::vector<Real> & /*yf*/,
681 : const RankFourTensor & /*Eijkl*/)
682 : {
683 0 : }
684 :
685 : Real
686 0 : MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & stress_params,
687 : const std::vector<Real> & intnl) const
688 : {
689 0 : std::vector<Real> yfs(_num_yf);
690 0 : yieldFunctionValuesV(stress_params, intnl, yfs);
691 0 : return yieldF(yfs);
692 : }
693 :
694 : Real
695 1428532 : MultiParameterPlasticityStressUpdate::yieldF(const std::vector<Real> & yfs) const
696 : {
697 1428532 : Real yf = yfs[0];
698 5373948 : for (std::size_t i = 1; i < yfs.size(); ++i)
699 3945416 : if (yf >= yfs[i] + _smoothing_tol)
700 : // no smoothing is needed, and yf is the biggest yield function
701 2946294 : continue;
702 999122 : else if (yfs[i] >= yf + _smoothing_tol)
703 : // no smoothing is needed, and yfs[i] is the biggest yield function
704 : yf = yfs[i];
705 : else
706 283248 : yf = 0.5 * (yf + yfs[i] + _smoothing_tol) + ismoother(yf - yfs[i]);
707 1428532 : return yf;
708 : }
709 :
710 : void
711 0 : MultiParameterPlasticityStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
712 : const std::vector<Real> & intnl_old,
713 : std::vector<Real> & stress_params,
714 : Real & gaE,
715 : std::vector<Real> & intnl) const
716 : {
717 0 : gaE = 0.0;
718 0 : std::copy(trial_stress_params.begin(), trial_stress_params.end(), stress_params.begin());
719 0 : std::copy(intnl_old.begin(), intnl_old.end(), intnl.begin());
720 0 : }
721 :
722 : void
723 0 : MultiParameterPlasticityStressUpdate::consistentTangentOperatorV(
724 : const RankTwoTensor & stress_trial,
725 : const std::vector<Real> & /*trial_stress_params*/,
726 : const RankTwoTensor & stress,
727 : const std::vector<Real> & /*stress_params*/,
728 : Real gaE,
729 : const yieldAndFlow & smoothed_q,
730 : const RankFourTensor & elasticity_tensor,
731 : bool compute_full_tangent_operator,
732 : const std::vector<std::vector<Real>> & dvar_dtrial,
733 : RankFourTensor & cto)
734 : {
735 0 : cto = elasticity_tensor;
736 0 : if (!compute_full_tangent_operator)
737 0 : return;
738 :
739 0 : const Real ga = gaE / _En;
740 :
741 0 : const std::vector<RankTwoTensor> dsp = dstress_param_dstress(stress);
742 0 : const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
743 :
744 0 : for (unsigned a = 0; a < _num_sp; ++a)
745 : {
746 0 : for (unsigned b = 0; b < _num_sp; ++b)
747 : {
748 0 : const RankTwoTensor t = elasticity_tensor * dsp_trial[a];
749 0 : RankTwoTensor s = _Cij[b][a] * dsp[b];
750 0 : for (unsigned c = 0; c < _num_sp; ++c)
751 0 : s -= dsp[b] * _Cij[b][c] * dvar_dtrial[c][a];
752 0 : s = elasticity_tensor * s;
753 0 : cto -= s.outerProduct(t);
754 : }
755 : }
756 :
757 0 : const std::vector<RankFourTensor> d2sp = d2stress_param_dstress(stress);
758 0 : RankFourTensor Tijab;
759 0 : for (unsigned i = 0; i < _num_sp; ++i)
760 0 : Tijab += smoothed_q.dg[i] * d2sp[i];
761 0 : Tijab = ga * elasticity_tensor * Tijab;
762 :
763 0 : RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
764 : try
765 : {
766 0 : inv = inv.transposeMajor().invSymm();
767 : }
768 0 : catch (const MooseException & e)
769 : {
770 : // Cannot form the inverse, so probably at some degenerate place in stress space.
771 : // Just return with the "best estimate" of the cto.
772 0 : mooseWarning("MultiParameterPlasticityStressUpdate: Cannot invert 1+T in consistent tangent "
773 : "operator computation at quadpoint ",
774 0 : _qp,
775 : " of element ",
776 0 : _current_elem->id());
777 : return;
778 0 : }
779 :
780 0 : cto = (cto.transposeMajor() * inv).transposeMajor();
781 : }
782 :
783 : void
784 118256 : MultiParameterPlasticityStressUpdate::setInelasticStrainIncrementAfterReturn(
785 : const RankTwoTensor & /*stress_trial*/,
786 : Real gaE,
787 : const yieldAndFlow & smoothed_q,
788 : const RankFourTensor & /*elasticity_tensor*/,
789 : const RankTwoTensor & returned_stress,
790 : RankTwoTensor & inelastic_strain_increment) const
791 : {
792 118256 : const std::vector<RankTwoTensor> dsp_dstress = dstress_param_dstress(returned_stress);
793 118256 : inelastic_strain_increment = RankTwoTensor();
794 473024 : for (unsigned i = 0; i < _num_sp; ++i)
795 354768 : inelastic_strain_increment += smoothed_q.dg[i] * dsp_dstress[i];
796 118256 : inelastic_strain_increment *= gaE / _En;
797 118256 : }
798 :
799 : Real
800 2509806 : MultiParameterPlasticityStressUpdate::calculateRes2(const std::vector<Real> & rhs) const
801 : {
802 : Real res2 = 0.0;
803 10912168 : for (const auto & r : rhs)
804 8402362 : res2 += r * r;
805 2509806 : return res2;
806 : }
807 :
808 : void
809 2509806 : MultiParameterPlasticityStressUpdate::calculateRHS(const std::vector<Real> & trial_stress_params,
810 : const std::vector<Real> & stress_params,
811 : Real gaE,
812 : const yieldAndFlow & smoothed_q,
813 : std::vector<Real> & rhs) const
814 : {
815 2509806 : const Real ga = gaE / _En;
816 8402362 : for (unsigned i = 0; i < _num_sp; ++i)
817 : {
818 5892556 : rhs[i] = stress_params[i] - trial_stress_params[i];
819 20296500 : for (unsigned j = 0; j < _num_sp; ++j)
820 14403944 : rhs[i] += ga * _Eij[i][j] * smoothed_q.dg[j];
821 : }
822 2509806 : rhs[_num_sp] = smoothed_q.f;
823 2509806 : }
824 :
825 : void
826 1705642 : MultiParameterPlasticityStressUpdate::dnRHSdVar(const yieldAndFlow & smoothed_q,
827 : const std::vector<std::vector<Real>> & dintnl,
828 : const std::vector<Real> & /*stress_params*/,
829 : Real gaE,
830 : std::vector<double> & jac) const
831 : {
832 20563168 : for (auto & jac_entry : jac)
833 18857526 : jac_entry = 0.0;
834 :
835 1705642 : const Real ga = gaE / _En;
836 :
837 : unsigned ind = 0;
838 5617890 : for (unsigned var = 0; var < _num_sp; ++var)
839 : {
840 13239636 : for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
841 : {
842 9327388 : if (var == rhs)
843 3912248 : jac[ind] -= 1.0;
844 32490840 : for (unsigned j = 0; j < _num_sp; ++j)
845 : {
846 23163452 : jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g[j][var];
847 68388756 : for (unsigned k = 0; k < _num_intnl; ++k)
848 45225304 : jac[ind] -= ga * _Eij[rhs][j] * smoothed_q.d2g_di[j][k] * dintnl[k][var];
849 : }
850 9327388 : ind++;
851 : }
852 : // now rhs = _num_sp (that is, the yield function)
853 3912248 : jac[ind] -= smoothed_q.df[var];
854 11614344 : for (unsigned k = 0; k < _num_intnl; ++k)
855 7702096 : jac[ind] -= smoothed_q.df_di[k] * dintnl[k][var];
856 3912248 : ind++;
857 : }
858 :
859 : // now var = _num_sp (that is, gaE)
860 5617890 : for (unsigned rhs = 0; rhs < _num_sp; ++rhs)
861 : {
862 13239636 : for (unsigned j = 0; j < _num_sp; ++j)
863 9327388 : jac[ind] -= (1.0 / _En) * _Eij[rhs][j] * smoothed_q.dg[j];
864 3912248 : ind++;
865 : }
866 : // now rhs = _num_sp (that is, the yield function)
867 1705642 : jac[ind] = 0.0;
868 1705642 : }
869 :
870 : void
871 625568 : MultiParameterPlasticityStressUpdate::dVardTrial(bool elastic_only,
872 : const std::vector<Real> & trial_stress_params,
873 : const std::vector<Real> & stress_params,
874 : Real gaE,
875 : const std::vector<Real> & intnl,
876 : const yieldAndFlow & smoothed_q,
877 : Real step_size,
878 : bool compute_full_tangent_operator,
879 : std::vector<std::vector<Real>> & dvar_dtrial) const
880 : {
881 625568 : if (!_fe_problem.currentlyComputingJacobian())
882 527840 : return;
883 :
884 97728 : if (!compute_full_tangent_operator)
885 : return;
886 :
887 97728 : if (elastic_only)
888 : {
889 : // no change to gaE, and all off-diag stuff remains unchanged from previous step
890 0 : for (unsigned v = 0; v < _num_sp; ++v)
891 0 : dvar_dtrial[v][v] += step_size;
892 : return;
893 : }
894 :
895 97728 : const Real ga = gaE / _En;
896 :
897 97728 : std::vector<std::vector<Real>> dintnl(_num_intnl, std::vector<Real>(_num_sp));
898 97728 : setIntnlDerivativesV(trial_stress_params, stress_params, intnl, dintnl);
899 :
900 : // rhs is described elsewhere, the following are changes in rhs wrt the trial_stress_param
901 : // values
902 : // In the following we use d(intnl)/d(trial variable) = - d(intnl)/d(variable)
903 97728 : std::vector<Real> rhs_cto((_num_sp + 1) * _num_sp);
904 :
905 : unsigned ind = 0;
906 294368 : for (unsigned a = 0; a < _num_sp; ++a)
907 : {
908 : // change in RHS[b] wrt changes in stress_param_trial[a]
909 593472 : for (unsigned b = 0; b < _num_sp; ++b)
910 : {
911 396832 : if (a == b)
912 196640 : rhs_cto[ind] -= 1.0;
913 1201152 : for (unsigned j = 0; j < _num_sp; ++j)
914 2406048 : for (unsigned k = 0; k < _num_intnl; ++k)
915 1601728 : rhs_cto[ind] -= ga * _Eij[b][j] * smoothed_q.d2g_di[j][k] * dintnl[k][a];
916 396832 : ind++;
917 : }
918 : // now b = _num_sp (that is, the yield function)
919 589152 : for (unsigned k = 0; k < _num_intnl; ++k)
920 392512 : rhs_cto[ind] -= smoothed_q.df_di[k] * dintnl[k][a];
921 196640 : ind++;
922 : }
923 :
924 : // jac = d(-rhs)/d(var)
925 97728 : std::vector<double> jac((_num_sp + 1) * (_num_sp + 1));
926 97728 : dnRHSdVar(smoothed_q, dintnl, stress_params, gaE, jac);
927 :
928 97728 : std::vector<PetscBLASInt> ipiv(_num_sp + 1);
929 : PetscBLASInt info;
930 97728 : const PetscBLASInt gesv_num_rhs = _num_sp + 1;
931 97728 : const PetscBLASInt gesv_num_pq = _num_sp;
932 97728 : LAPACKgesv_(&gesv_num_rhs,
933 : &gesv_num_pq,
934 : &jac[0],
935 : &gesv_num_rhs,
936 : &ipiv[0],
937 : &rhs_cto[0],
938 : &gesv_num_rhs,
939 : &info);
940 97728 : if (info != 0)
941 0 : errorHandler("MultiParameterPlasticityStressUpdate: PETSC LAPACK gsev routine returned with "
942 0 : "error code " +
943 0 : Moose::stringify(info));
944 :
945 : ind = 0;
946 97728 : std::vector<std::vector<Real>> dvarn_dtrialn(_num_sp + 1, std::vector<Real>(_num_sp, 0.0));
947 294368 : for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
948 : {
949 593472 : for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
950 : {
951 396832 : dvarn_dtrialn[v][spt] = rhs_cto[ind];
952 396832 : ind++;
953 : }
954 : // the final NR variable is gaE
955 196640 : dvarn_dtrialn[_num_sp][spt] = rhs_cto[ind];
956 196640 : ind++;
957 : }
958 :
959 97728 : const std::vector<std::vector<Real>> dvar_dtrial_old = dvar_dtrial;
960 :
961 294368 : for (unsigned v = 0; v < _num_sp; ++v) // loop over variables in NR procedure
962 : {
963 593472 : for (unsigned spt = 0; spt < _num_sp; ++spt) // loop over trial stress-param variables
964 : {
965 396832 : dvar_dtrial[v][spt] = step_size * dvarn_dtrialn[v][spt];
966 1201152 : for (unsigned a = 0; a < _num_sp; ++a)
967 804320 : dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
968 : }
969 : }
970 : // for gaE the formulae are a little different
971 : const unsigned v = _num_sp;
972 294368 : for (unsigned spt = 0; spt < _num_sp; ++spt)
973 : {
974 196640 : dvar_dtrial[v][spt] += step_size * dvarn_dtrialn[v][spt]; // note +=
975 593472 : for (unsigned a = 0; a < _num_sp; ++a)
976 396832 : dvar_dtrial[v][spt] += dvarn_dtrialn[v][a] * dvar_dtrial_old[a][spt];
977 : }
978 195456 : }
979 :
980 : bool
981 1607914 : MultiParameterPlasticityStressUpdate::precisionLoss(const std::vector<Real> & solution,
982 : const std::vector<Real> & stress_params,
983 : Real gaE) const
984 : {
985 1607914 : if (std::abs(solution[_num_sp]) > 1E-13 * std::abs(gaE))
986 : return false;
987 4432 : for (unsigned i = 0; i < _num_sp; ++i)
988 3852 : if (std::abs(solution[i]) > 1E-13 * std::abs(stress_params[i]))
989 : return false;
990 : return true;
991 : }
|
