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 : #pragma once
11 :
12 : #include "StressUpdateBase.h"
13 :
14 : #include <array>
15 :
16 : /**
17 : * MultiParameterPlasticityStressUpdate performs the return-map
18 : * algorithm and associated stress updates for plastic
19 : * models where the yield function and flow directions
20 : * depend on multiple parameters (called "stress_params" in the
21 : * documentation and sp in the code) that
22 : * are themselves functions of stress.
23 : *
24 : * Let the stress_params be S = {S[0], S[1], ... S[N-1]}
25 : * and define _num_sp = N.
26 : *
27 : * For instance, CappedDruckerPrager plasticity has _num_sp = 2
28 : * and defines
29 : * S[0] = p = stress.trace()
30 : * S[1] = q = sqrt(stress.secondInvariant())
31 : *
32 : * The point of the stress_params is to reduce the number of
33 : * degrees of freedom involved in the return-map algorithm
34 : * (for computational efficiency as all the intensive computation
35 : * occurs in "stress_param space") and to allow users to
36 : * write their yield functions, etc, in forms that are
37 : * clear to them to avoid errors.
38 : *
39 : * Derived classes can describe plasticity via multiple yield
40 : * functions. The number of yield function is _num_yf.
41 : * The yield function(s) and flow potential(s)
42 : * must be functions of the stress_params.
43 : *
44 : * Derived classes can use any number of internal parameters.
45 : * The number of internal parameters is _num_intnl. The
46 : * derived classes must define the evolution of these internal
47 : * parameters, which are typically functions of certain
48 : * "distances" in the return-mapping procedure.
49 : *
50 : * For instance, CappedDruckerPrager plasticity has three
51 : * yield functions (a smoothed DruckerPrager cone, a tensile failure
52 : * envelope, and a compressive failure envelope) and two
53 : * internal parameters (plastic shear strain and plastic tensile strain).
54 : * The strength parameters (cohesion, friction angle,
55 : * dilation angle, tensile strength, compressive strength)
56 : * are functions of these internal parameters.
57 : *
58 : * A novel smoothing procedure is used to smooth the derived class's
59 : * yield functions into a single, smooth yield function. The
60 : * smoothing is also performed for the flow potential. All
61 : * return-mapping, etc, processes are performed on this single
62 : * yield function and flow potential.
63 : *
64 : * The return-map algorithm implements the updateState method of
65 : * StressUpdateBase. In particular, the following system of
66 : * equations is solved:
67 : * 0 = _rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j]
68 : * 0 = _rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j]
69 : * ...
70 : * 0 = _rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j]
71 : * 0 = _rhs[N] = f(S, intnl)
72 : * (here there is an implied sum over j)
73 : * for the _num_sp stress_params S, and the single scalar ga
74 : * (actually I solve for gaE = ga * _En, where _En is a normalising
75 : * factor so that gaE is of similar magnitude to the S variables).
76 : * In general E[a, b] = dS[a]/dstress(i, j) E(i, j, k, l) dS[b]/dstress(k, l),
77 : * and in the code below it is assumed that the E[a, b] are independent
78 : * of the stress_params, but adding a dependence would only require
79 : * some Jacobian terms to be modified.
80 : * There are N+1 rhs equations to solve. Here f is
81 : * the smoothed yield function, so the last equation is the admissibility
82 : * condition (the returned stress lies on the yield surface) and g
83 : * is the flow potential so the other conditions implement the normality
84 : * condition. It is up to the derived classes to defined E[i, j]
85 : * so that the rhs really represents the normality condition.
86 : *
87 : * For instance, CappedDruckerPrager has (with Eijkl being the elasticity
88 : * tensor):
89 : * E[0, 0] = _Epp = Eijkl.sum3x3()
90 : * E[0, 1] = 0
91 : * E[1, 0] = 0
92 : * E[1, 1] = Eijkl(0, 1, 0, 1)
93 : */
94 : class MultiParameterPlasticityStressUpdate : public StressUpdateBase
95 : {
96 : public:
97 : static InputParameters validParams();
98 :
99 : MultiParameterPlasticityStressUpdate(const InputParameters & parameters,
100 : unsigned num_sp,
101 : unsigned num_yf,
102 : unsigned num_intnl);
103 :
104 : protected:
105 : virtual void initQpStatefulProperties() override;
106 : using StressUpdateBase::updateState;
107 : virtual void updateState(RankTwoTensor & strain_increment,
108 : RankTwoTensor & inelastic_strain_increment,
109 : const RankTwoTensor & rotation_increment,
110 : RankTwoTensor & stress_new,
111 : const RankTwoTensor & stress_old,
112 : const RankFourTensor & elasticity_tensor,
113 : const RankTwoTensor & elastic_strain_old,
114 : bool compute_full_tangent_operator,
115 : RankFourTensor & tangent_operator) override;
116 :
117 : virtual void propagateQpStatefulProperties() override;
118 :
119 1683 : virtual TangentCalculationMethod getTangentCalculationMethod() override
120 : {
121 1683 : return TangentCalculationMethod::FULL;
122 : }
123 :
124 : /// Internal dimensionality of tensors (currently this is 3 throughout solid mechanics)
125 : constexpr static unsigned _tensor_dimensionality = 3;
126 :
127 : /// Number of stress parameters
128 : const unsigned _num_sp;
129 :
130 : /// An admissible value of stress_params at the initial time
131 : const std::vector<Real> _definitely_ok_sp;
132 :
133 : /// E[i, j] in the system of equations to be solved
134 : std::vector<std::vector<Real>> _Eij;
135 :
136 : /// normalising factor
137 : Real _En;
138 :
139 : /// _Cij[i, j] * _Eij[j, k] = 1 iff j == k
140 : std::vector<std::vector<Real>> _Cij;
141 :
142 : /// Number of yield functions
143 : const unsigned _num_yf;
144 :
145 : /// Number of internal parameters
146 : const unsigned _num_intnl;
147 :
148 : /// Maximum number of Newton-Raphson iterations allowed in the return-map process
149 : const unsigned _max_nr_its;
150 :
151 : /// Whether to perform finite-strain rotations
152 : const bool _perform_finite_strain_rotations;
153 :
154 : /// Smoothing tolerance: edges of the yield surface get smoothed by this amount
155 : const Real _smoothing_tol;
156 :
157 : /// Square of the smoothing tolerance
158 : const Real _smoothing_tol2;
159 :
160 : /// The yield-function tolerance
161 : const Real _f_tol;
162 :
163 : /// Square of the yield-function tolerance
164 : const Real _f_tol2;
165 :
166 : /**
167 : * In order to help the Newton-Raphson procedure, the applied
168 : * strain increment may be applied in sub-increments of size
169 : * greater than this value.
170 : */
171 : const Real _min_step_size;
172 :
173 : /// handles case of initial_stress that is inadmissible being supplied
174 : bool _step_one;
175 :
176 : /// Output a warning message if precision loss is encountered during the return-map process
177 : const bool _warn_about_precision_loss;
178 :
179 : /// plastic strain
180 : MaterialProperty<RankTwoTensor> & _plastic_strain;
181 :
182 : /// Old value of plastic strain
183 : const MaterialProperty<RankTwoTensor> & _plastic_strain_old;
184 :
185 : /// internal parameters
186 : MaterialProperty<std::vector<Real>> & _intnl;
187 :
188 : /// old values of internal parameters
189 : const MaterialProperty<std::vector<Real>> & _intnl_old;
190 :
191 : /// yield functions
192 : MaterialProperty<std::vector<Real>> & _yf;
193 :
194 : /// Number of Newton-Raphson iterations used in the return-map
195 : MaterialProperty<Real> & _iter;
196 :
197 : /// Maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation
198 : MaterialProperty<Real> & _max_iter_used;
199 :
200 : /// Old value of maximum number of Newton-Raphson iterations used in the return-map during the course of the entire simulation
201 : const MaterialProperty<Real> & _max_iter_used_old;
202 :
203 : /// Whether a line-search was needed in the latest Newton-Raphson process (1 if true, 0 otherwise)
204 : MaterialProperty<Real> & _linesearch_needed;
205 :
206 : /**
207 : * Struct designed to hold info about a single yield function
208 : * and its derivatives, as well as the flow directions
209 : */
210 : struct yieldAndFlow
211 : {
212 : Real f; // yield function value
213 : std::vector<Real> df; // df/d(stress_param[i])
214 : std::vector<Real> df_di; // df/d(intnl_variable[a])
215 : std::vector<Real> dg; // d(flow)/d(stress_param[i])
216 : std::vector<std::vector<Real>> d2g; // d^2(flow)/d(sp[i])/d(sp[j])
217 : std::vector<std::vector<Real>> d2g_di; // d^2(flow)/d(sp[i])/d(intnl[a])
218 :
219 : yieldAndFlow() : yieldAndFlow(0, 0) {}
220 :
221 3618 : yieldAndFlow(unsigned num_var, unsigned num_intnl)
222 3618 : : f(0.0),
223 3618 : df(num_var),
224 3618 : df_di(num_intnl),
225 3618 : dg(num_var),
226 3618 : d2g(num_var, std::vector<Real>(num_var, 0.0)),
227 3618 : d2g_di(num_var, std::vector<Real>(num_intnl, 0.0))
228 : {
229 3618 : }
230 :
231 : // Zero everything without resizing/reallocating
232 12560088 : void reset()
233 : {
234 12560088 : f = 0.0;
235 : std::fill(df.begin(), df.end(), 0.0);
236 : std::fill(df_di.begin(), df_di.end(), 0.0);
237 : std::fill(dg.begin(), dg.end(), 0.0);
238 46468368 : for (auto & row : d2g)
239 : std::fill(row.begin(), row.end(), 0.0);
240 46468368 : for (auto & row : d2g_di)
241 : std::fill(row.begin(), row.end(), 0.0);
242 12560088 : }
243 :
244 : // this may be involved in the smoothing of a group of yield functions
245 : bool operator<(const yieldAndFlow & fd) const { return f < fd.f; }
246 : };
247 :
248 : /**
249 : * Computes the smoothed yield function
250 : * @param stress_params The stress parameters (eg stress_params[0] = stress_zz and
251 : * stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
252 : * @param intnl The internal parameters (eg intnl[0] is shear, intnl[1] is tensile)
253 : * @return The smoothed yield function value
254 : */
255 : Real yieldF(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;
256 :
257 : /**
258 : * Computes the smoothed yield function
259 : * @param yfs The values of the individual yield functions
260 : * @return The smoothed yield function value
261 : */
262 : Real yieldF(const std::vector<Real> & yfs) const;
263 :
264 : /**
265 : * Smooths yield functions. The returned value must be zero if abs(f_diff) >= _smoothing_tol
266 : * and otherwise must satisfy, over -_smoothing_tol <= f_diff <= _smoothing_tol:
267 : * (1) C2
268 : * (2) zero at f_diff = +/- _smoothing_tol
269 : * (3) derivative is +/-0.5 at f_diff = +/- _smoothing_tol
270 : * (4) derivative must be in [-0.5, 0.5]
271 : * (5) second derivative is zero at f_diff = +/- _smoothing_tol
272 : * (6) second derivative must be non-negative
273 : * in order to ensure C2 differentiability and convexity of the smoothed yield surface.
274 : */
275 : Real ismoother(Real f_diff) const;
276 :
277 : /**
278 : * Derivative of ismoother
279 : */
280 : Real smoother(Real f_diff) const;
281 :
282 : /**
283 : * Derivative of smoother
284 : */
285 : Real dsmoother(Real f_diff) const;
286 :
287 : /**
288 : * Calculates all yield functions and derivatives, and then performs the smoothing scheme
289 : * @param stress_params[in] The stress parameters (eg stress_params[0] = stress_zz and
290 : * stress_params[1] = sqrt(stress_zx^2 + stress_zy^2))
291 : * @param intnl[in] Internal parameters
292 : * @return The smoothed yield function and derivatives
293 : */
294 : yieldAndFlow smoothAllQuantities(const std::vector<Real> & stress_params,
295 : const std::vector<Real> & intnl) const;
296 :
297 : /**
298 : * Performs a line-search to find stress_params and gaE
299 : * Upon entry:
300 : * - rhs contains the *solution* to the Newton-Raphson (ie nrStep should have been called). If a
301 : * full
302 : * Newton step is used then stress_params[:] += rhs[0:_num_sp-1] and gaE += rhs[_num_sp]
303 : * - res2 contains the residual-squared before applying any of solution
304 : * - stress_params contains the stress_params before applying any of the solution
305 : * - gaE contains gaE before applying any of the solution (that is contained in rhs)
306 : * Upon exit:
307 : * - stress_params will be the stress_params after applying the solution
308 : * - gaE will be the stress_params after applying the solution
309 : * - rhs will contain the updated rhs values (after applying the solution) ready for the next
310 : * Newton-Raphson step,
311 : * - res2 will be the residual-squared after applying the solution
312 : * - intnl will contain the internal variables corresponding to the return from
313 : * trial_stress_params to stress_params
314 : * (and starting from intnl_ok)
315 : * - linesearch_needed will be 1.0 if a linesearch was needed
316 : * - smoothed_q will contain the value of the yield function and its derivatives, etc, at
317 : * (stress_params, intnl)
318 : * @param res2[in,out] the residual-squared, both as an input and output
319 : * @param stress_params[in,out] Upon input the value of the stress_params before the current
320 : * Newton-Raphson process was initiated.
321 : * Upon exit this will hold the values coming from the line search.
322 : * @param trial_stress_params[in] Trial value for the stress_params for this (sub)strain increment
323 : * @param gaE[in,out] Upon input the value of gaE before the current Newton-Raphson iteration was
324 : * initiated. Upon exit this will hold
325 : * the value coming from the line-search
326 : * @param smoothed_q[in,out] Upon input, the value of the smoothed yield function and derivatives
327 : * at the
328 : * prior-to-Newton configuration. Upon exit this is evaluated at the new (stress_params, intnl)
329 : * @param intnl_ok[in] The value of the internal parameters from the start of this (sub)strain
330 : * increment
331 : * @param intnl[in,out] The value of the internal parameters after the line-search has converged
332 : * @param rhs[in,out] Upon entry this contains the solution to the Newton-Raphson. Upon exit this
333 : * contains the updated rhs values
334 : * @return 0 if successful, 1 otherwise
335 : */
336 : int lineSearch(Real & res2,
337 : std::vector<Real> & stress_params,
338 : Real & gaE,
339 : const std::vector<Real> & trial_stress_params,
340 : yieldAndFlow & smoothed_q,
341 : const std::vector<Real> & intnl_ok,
342 : std::vector<Real> & intnl,
343 : std::vector<Real> & rhs,
344 : Real & linesearch_needed) const;
345 :
346 : /**
347 : * Performs a Newton-Raphson step to attempt to zero rhs
348 : * Upon return, rhs will contain the solution.
349 : * @param smoothed_q[in] The value of the smoothed yield function and derivatives prior to this
350 : * Newton-Raphson step
351 : * @param trial_stress_params[in] Trial value for the stress_params for this (sub)strain increment
352 : * @param stress_params[in] The current value of the stress_params
353 : * @param intnl[in] The current value of the internal parameters
354 : * @param gaE[in] The current value of gaE
355 : * @param rhs[in,out] Upon entry, the rhs to zero using Newton-Raphson. Upon exit, the solution
356 : * to the Newton-Raphson problem
357 : * @return 0 if successful, 1 otherwise
358 : */
359 : int nrStep(const yieldAndFlow & smoothed_q,
360 : const std::vector<Real> & trial_stress_params,
361 : const std::vector<Real> & stress_params,
362 : const std::vector<Real> & intnl,
363 : Real gaE,
364 : std::vector<Real> & rhs) const;
365 :
366 : /**
367 : * Calculates the residual-squared for the Newton-Raphson + line-search
368 : * @param rhs[in] The RHS vector
369 : * @return sum_i (rhs[i] * rhs[i])
370 : */
371 : Real calculateRes2(const std::vector<Real> & rhs) const;
372 :
373 : /**
374 : * Calculates the RHS in the following
375 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, j] * dg/dS[j]
376 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, j] * dg/dS[j]
377 : * ...
378 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, j] * dg/dS[j]
379 : * 0 = rhs[N] = f(S, intnl)
380 : * where N = _num_sp
381 : * @param trial_stress_params[in] The trial stress parameters for this (sub)strain increment,
382 : * S[:]^trial
383 : * @param stress_params[in] The current stress parameters, S[:]
384 : * @param gaE[in] ga*_En (the normalisation with _En is so that gaE is of similar magnitude to S)
385 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
386 : * the current stress parameters and the current internal parameters
387 : * @param rhs[out] The result
388 : */
389 : void calculateRHS(const std::vector<Real> & trial_stress_params,
390 : const std::vector<Real> & stress_params,
391 : Real gaE,
392 : const yieldAndFlow & smoothed_q,
393 : std::vector<Real> & rhs) const;
394 :
395 : /**
396 : * Derivative of -RHS with respect to the stress_params and gaE, placed
397 : * into an array ready for solving the linear system using
398 : * LAPACK gsev
399 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
400 : * the
401 : * current values of the stress_params and the internal parameters
402 : * @param dintnl[in] The derivatives of the internal parameters wrt the stress_params
403 : * @param stress_params[in] The current value of the stress_params during the Newton-Raphson
404 : * process
405 : * @param gaE[in] The current value of gaE
406 : * @param jac[out] The outputted derivatives
407 : */
408 : void dnRHSdVar(const yieldAndFlow & smoothed_q,
409 : const std::vector<std::vector<Real>> & dintnl,
410 : const std::vector<Real> & stress_params,
411 : Real gaE,
412 : std::vector<double> & jac) const;
413 :
414 : /**
415 : * Performs any necessary cleaning-up, then throw MooseException(message)
416 : * @param message The message to using in MooseException
417 : */
418 : virtual void errorHandler(const std::string & message) const;
419 :
420 : /**
421 : * Computes the values of the yield functions, given stress_params and intnl parameters.
422 : * Derived classes must override this, to provide the values of the yield functions
423 : * in yf.
424 : * @param stress_params[in] The stress parameters
425 : * @param intnl[in] The internal parameters
426 : * @param[out] yf The yield function values
427 : */
428 : virtual void yieldFunctionValuesV(const std::vector<Real> & stress_params,
429 : const std::vector<Real> & intnl,
430 : std::vector<Real> & yf) const = 0;
431 :
432 : /**
433 : * Completely fills all_q with correct values. These values are:
434 : * (1) the yield function values, yf[i]
435 : * (2) d(yf[i])/d(stress_params[j])
436 : * (3) d(yf[i])/d(intnl[j])
437 : * (4) d(flowPotential[i])/d(stress_params[j])
438 : * (5) d2(flowPotential[i])/d(stress_params[j])/d(stress_params[k])
439 : * (6) d2(flowPotential[i])/d(stress_params[j])/d(intnl[k])
440 : * @param stress_params[in] The stress parameters
441 : * @param intnl[in] The internal parameters
442 : * @param[out] all_q All the desired quantities
443 : */
444 : virtual void computeAllQV(const std::vector<Real> & stress_params,
445 : const std::vector<Real> & intnl,
446 : std::vector<yieldAndFlow> & all_q) const = 0;
447 :
448 : /**
449 : * Derived classes may employ this function to record stuff or do
450 : * other computations prior to the return-mapping algorithm. We
451 : * know that (trial_stress_params, intnl_old) is inadmissible when
452 : * this is called
453 : * @param trial_stress_params[in] The trial values of the stress parameters
454 : * @param stress_trial[in] Trial stress tensor
455 : * @param intnl_old[in] Old value of the internal parameters.
456 : * @param yf[in] The yield functions at (p_trial, q_trial, intnl_old)
457 : * @param Eijkl[in] The elasticity tensor
458 : */
459 : virtual void preReturnMapV(const std::vector<Real> & trial_stress_params,
460 : const RankTwoTensor & stress_trial,
461 : const std::vector<Real> & intnl_old,
462 : const std::vector<Real> & yf,
463 : const RankFourTensor & Eijkl);
464 :
465 : /**
466 : * Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
467 : * The purpose of these "good guesses" is to speed up the Newton-Raphson process by providing
468 : * it with a good initial guess.
469 : * Derived classes may choose to override this if their plastic models are easy
470 : * enough to solve approximately.
471 : * The default values, provided by this class, are simply gaE = 0, stress_params =
472 : * trial_stress_params,
473 : * that is, the "good guess" is just the trial point for this (sub)strain increment.
474 : * @param trial_stress_params[in] The stress_params at the trial point
475 : * @param intnl_old[in] The internal parameters before applying the (sub)strain increment
476 : * @param stress_params[out] The "good guess" value of the stress_params
477 : * @param gaE[out] The "good guess" value of gaE
478 : * @param intnl[out] The "good guess" value of the internal parameters
479 : */
480 : virtual void initializeVarsV(const std::vector<Real> & trial_stress_params,
481 : const std::vector<Real> & intnl_old,
482 : std::vector<Real> & stress_params,
483 : Real & gaE,
484 : std::vector<Real> & intnl) const;
485 :
486 : /**
487 : * Sets the internal parameters based on the trial values of
488 : * stress_params, their current values, and the old values of the
489 : * internal parameters.
490 : * Derived classes must override this.
491 : * @param trial_stress_params[in] The trial stress parameters (eg trial_p and trial_q)
492 : * @param current_stress_params[in] The current stress parameters (eg p and q)
493 : * @param intnl_old[out] Old value of internal parameters
494 : * @param intnl[out] The value of internal parameters to be set
495 : */
496 : virtual void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
497 : const std::vector<Real> & current_stress_params,
498 : const std::vector<Real> & intnl_old,
499 : std::vector<Real> & intnl) const = 0;
500 :
501 : /**
502 : * Sets the derivatives of internal parameters, based on the trial values of
503 : * stress_params, their current values, and the current values of the
504 : * internal parameters.
505 : * Derived classes must override this.
506 : * @param trial_stress_params[in] The trial stress parameters
507 : * @param current_stress_params[in] The current stress parameters
508 : * @param intnl[in] The current value of the internal parameters
509 : * @param dintnl[out] The derivatives dintnl[i][j] = d(intnl[i])/d(stress_param j)
510 : */
511 : virtual void setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
512 : const std::vector<Real> & current_stress_params,
513 : const std::vector<Real> & intnl,
514 : std::vector<std::vector<Real>> & dintnl) const = 0;
515 :
516 : /**
517 : * Computes stress_params, given stress. Derived classes must
518 : * override this
519 : * @param stress[in] Stress tensor
520 : * @param stress_params[out] The compute stress_params
521 : */
522 : virtual void computeStressParams(const RankTwoTensor & stress,
523 : std::vector<Real> & stress_params) const = 0;
524 :
525 : /**
526 : * Derived classes may use this to perform calculations before
527 : * any return-map process is performed, for instance, to initialize
528 : * variables.
529 : * This is called at the very start of updateState, even before
530 : * any checking for admissible stresses, etc, is performed
531 : */
532 : virtual void initializeReturnProcess();
533 :
534 : /**
535 : * Derived classes may use this to perform calculations after the
536 : * return-map process has completed successfully in stress_param space
537 : * but before the returned stress tensor has been calculcated.
538 : * @param rotation_increment[in] The large-strain rotation increment
539 : */
540 : virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment);
541 :
542 : /**
543 : * Sets stress from the admissible parameters.
544 : * This is called after the return-map process has completed
545 : * successfully in stress_param space, just after finalizeReturnProcess
546 : * has been called.
547 : * Derived classes must override this function
548 : * @param stress_trial[in] The trial value of stress
549 : * @param stress_params[in] The value of the stress_params after the return-map process has
550 : * completed successfully
551 : * @param gaE[in] The value of gaE after the return-map process has completed successfully
552 : * @param intnl[in] The value of the internal parameters after the return-map process has
553 : * completed successfully
554 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
555 : * the returned state
556 : * @param Eijkl[in] The elasticity tensor
557 : * @param stress[out] The returned value of the stress tensor
558 : */
559 : virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial,
560 : const std::vector<Real> & stress_params,
561 : Real gaE,
562 : const std::vector<Real> & intnl,
563 : const yieldAndFlow & smoothed_q,
564 : const RankFourTensor & Eijkl,
565 : RankTwoTensor & stress) const = 0;
566 :
567 : /**
568 : * Sets inelastic strain increment from the returned configuration
569 : * This is called after the return-map process has completed
570 : * successfully in stress_param space, just after finalizeReturnProcess
571 : * has been called.
572 : * Derived classes may override this function
573 : * @param stress_trial[in] The trial value of stress
574 : * @param gaE[in] The value of gaE after the return-map process has completed successfully
575 : * @param smoothed_q[in] Holds the current value of yield function and derivatives evaluated at
576 : * the returned configuration
577 : * @param elasticity_tensor[in] The elasticity tensor
578 : * @param returned_stress[in] The stress after the return-map process
579 : * @param inelastic_strain_increment[out] The inelastic strain increment resulting from this
580 : * return-map
581 : */
582 : virtual void setInelasticStrainIncrementAfterReturn(const RankTwoTensor & stress_trial,
583 : Real gaE,
584 : const yieldAndFlow & smoothed_q,
585 : const RankFourTensor & elasticity_tensor,
586 : const RankTwoTensor & returned_stress,
587 : RankTwoTensor & inelastic_strain_increment);
588 :
589 : /**
590 : * Calculates the consistent tangent operator.
591 : * Derived classes may choose to override this for computational
592 : * efficiency. The implementation in this class is quite expensive,
593 : * even though it looks compact and clean, because of all the
594 : * manipulations of RankFourTensors involved.
595 : * @param stress_trial[in] the trial value of the stress tensor for this strain increment
596 : * @param trial_stress_params[in] the trial values of the stress_params for this strain increment
597 : * @param stress[in] the returned value of the stress tensor for this strain increment
598 : * @param stress_params[in] the returned value of the stress_params for this strain increment
599 : * @param gaE[in] the total value of that came from this strain increment
600 : * @param smoothed_q[in] contains the yield function and derivatives evaluated at (p, q)
601 : * @param Eijkl[in] The elasticity tensor
602 : * @param compute_full_tangent_operator[in] true if the full consistent tangent operator is
603 : * needed,
604 : * otherwise false
605 : * @param dvar_dtrial[in] dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])
606 : * for
607 : * this strain increment
608 : * @param[out] cto The consistent tangent operator
609 : */
610 : virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial,
611 : const std::vector<Real> & trial_stress_params,
612 : const RankTwoTensor & stress,
613 : const std::vector<Real> & stress_params,
614 : Real gaE,
615 : const yieldAndFlow & smoothed_q,
616 : const RankFourTensor & Eijkl,
617 : bool compute_full_tangent_operator,
618 : const std::vector<std::vector<Real>> & dvar_dtrial,
619 : RankFourTensor & cto);
620 :
621 : /**
622 : * d(stress_param[i])/d(stress) at given stress. When overriding, ensure to assert that dsp has
623 : * size _num_sp
624 : * @param stress[in] stress tensor
625 : * @param dsp[out] d(stress_param[:])/d(stress)
626 : */
627 : virtual void dstressparam_dstress(const RankTwoTensor & stress,
628 : std::vector<RankTwoTensor> & dsp) const = 0;
629 :
630 : /**
631 : * d2(stress_param[i])/d(stress)/d(stress) at given stress. When overriding, ensure to assert
632 : * that d2sp has size _num_sp
633 : * @param stress[in] stress tensor
634 : * @param d2sp[out] d2(stress_param[:])/d(stress)/d(stress)
635 : */
636 : virtual void d2stressparam_dstress(const RankTwoTensor & stress,
637 : std::vector<RankFourTensor> & d2sp) const = 0;
638 :
639 : /// Sets _Eij and _En and _Cij
640 : virtual void setEffectiveElasticity(const RankFourTensor & Eijkl) = 0;
641 :
642 : /**
643 : * Calculates derivatives of the stress_params and gaE with repect to the
644 : * trial values of the stress_params for the (sub)strain increment.
645 : * After the strain increment has been fully applied, dvar_dtrial will
646 : * contain the result appropriate to the full strain increment. Before
647 : * that time (if applying in sub-strain increments) it will contain the
648 : * result appropriate to the amount of strain increment applied successfully.
649 : * @param elastic_only[in] whether this was an elastic step: if so then the updates to dvar_dtrial
650 : * are fairly trivial
651 : * @param trial_stress_params[in] Trial values of stress_params for this (sub)strain increment
652 : * @param stress_params[in] Returned values of stress_params for this (sub)strain increment
653 : * @param gaE[in] the value of gaE that came from this (sub)strain increment
654 : * @param intnl[in] the value of the internal parameters at the returned position
655 : * @param smoothed_q[in] contains the yield function and derivatives evaluated at (stress_params,
656 : * intnl)
657 : * @param step_size[in] size of this (sub)strain increment
658 : * @param compute_full_tangent_operator[in] true if the full consistent tangent operator is
659 : * needed,
660 : * otherwise false
661 : * @param dvar_dtrial[out] dvar_dtrial[i][j] = d({stress_param[i],gaE})/d(trial_stress_param[j])
662 : */
663 : void dVardTrial(bool elastic_only,
664 : const std::vector<Real> & trial_stress_params,
665 : const std::vector<Real> & stress_params,
666 : Real gaE,
667 : const std::vector<Real> & intnl,
668 : const yieldAndFlow & smoothed_q,
669 : Real step_size,
670 : bool compute_full_tangent_operator,
671 : std::vector<std::vector<Real>> & dvar_dtrial) const;
672 :
673 : /**
674 : * Check whether precision loss has occurred
675 : * @param[in] solution The solution to the Newton-Raphson system
676 : * @param[in] stress_params The currect values of the stress_params for this (sub)strain increment
677 : * @param[in] gaE The currenct value of gaE for this (sub)strain increment
678 : * @return true if precision loss has occurred
679 : */
680 : bool precisionLoss(const std::vector<Real> & solution,
681 : const std::vector<Real> & stress_params,
682 : Real gaE) const;
683 :
684 : private:
685 : /**
686 : * "Trial" value of stress_params that initializes the return-map process
687 : * This is derived from stress = stress_old + Eijkl * strain_increment.
688 : * However, since the return-map process can fail and be restarted by
689 : * applying strain_increment in multiple substeps, _trial_sp can vary
690 : * from substep to substep.
691 : */
692 : std::vector<Real> _trial_sp;
693 :
694 : /**
695 : * "Trial" value of stress that is set at the beginning of the
696 : * return-map process. It is fixed at stress_old + Eijkl * strain_increment
697 : * irrespective of any sub-stepping
698 : */
699 : RankTwoTensor _stress_trial;
700 :
701 : /**
702 : * 0 = rhs[0] = S[0] - S[0]^trial + ga * E[0, i] * dg/dS[i]
703 : * 0 = rhs[1] = S[1] - S[1]^trial + ga * E[1, i] * dg/dS[i]
704 : * ...
705 : * 0 = rhs[N-1] = S[N-1] - S[N-1]^trial + ga * E[N-1, i] * dg/dS[i]
706 : * 0 = rhs[N] = f(S, intnl)
707 : * Here N = num_sp
708 : */
709 : std::vector<Real> _rhs;
710 :
711 : /**
712 : * d({stress_param[i], gaE})/d(trial_stress_param[j])
713 : */
714 : std::vector<std::vector<Real>> _dvar_dtrial;
715 :
716 : /**
717 : * The state (ok_sp, ok_intnl) is known to be admissible, so
718 : * ok_sp are stress_params that are "OK". If the strain_increment
719 : * is applied in substeps then ok_sp is updated after each
720 : * sub strain_increment is applied and the return-map is successful.
721 : * At the end of the entire return-map process _ok_sp will contain
722 : * the stress_params where (_ok_sp, _intnl) is admissible.
723 : */
724 : std::vector<Real> _ok_sp;
725 :
726 : /**
727 : * The state (ok_sp, ok_intnl) is known to be admissible
728 : */
729 : std::vector<Real> _ok_intnl;
730 :
731 : /**
732 : * _del_stress_params = trial_stress_params - ok_sp
733 : * This is fixed at the beginning of the return-map process,
734 : * irrespective of substepping. The return-map problem is:
735 : * apply del_stress_params to stress_prams, and then find
736 : * an admissible (returned) stress_params and gaE
737 : */
738 : std::vector<Real> _del_stress_params;
739 :
740 : /**
741 : * The current values of the stress params during the Newton-Raphson
742 : */
743 : std::vector<Real> _current_sp;
744 :
745 : /**
746 : * The current values of the internal params during the Newton-Raphson
747 : */
748 : std::vector<Real> _current_intnl;
749 :
750 : /**
751 : * Current values of the yield function, derivatives, etc during updateState.
752 : * During Newton-Raphson convergence, this is potentially re-calculated multiple times within the
753 : * loops of updateState and by the lineSearch function that is called within those loops.
754 : */
755 : yieldAndFlow _smoothed_q;
756 :
757 : /**
758 : * d(stress_param[:])/d(stress) used in dstressparam_dstress
759 : * to avoid repeatedly allocating/deallocating this vector
760 : */
761 : mutable std::vector<RankTwoTensor> _dsp_scratch;
762 :
763 : /**
764 : * d(stress_param[:])/d(stress_trial) used in dstressparam_dstress
765 : * to avoid repeatedly allocating/deallocating this vector
766 : */
767 : mutable std::vector<RankTwoTensor> _dsp_trial_scratch;
768 :
769 : /**
770 : * d2(stress_param[:])/d(stress)/d(stress) used in d2stressparam_dstress
771 : * to avoid repeatedly allocating/deallocating this vector
772 : */
773 : mutable std::vector<RankFourTensor> _d2sp_scratch;
774 :
775 : /**
776 : * values of yield functions in yieldF,
777 : * to avoid repeatedly allocating/deallocating this vector
778 : */
779 : mutable std::vector<Real> _yfs_scratch;
780 :
781 : /**
782 : * container to hold old values of old stress_params in lineSearch
783 : * to avoid repeatedly allocating/deallocating this vector
784 : */
785 : mutable std::vector<Real> _sp_params_old_scratch;
786 :
787 : /**
788 : * container to hold initial values of rhs in lineSearch
789 : * to avoid repeatedly allocating/deallocating this vector
790 : */
791 : mutable std::vector<Real> _delta_nr_params_scratch;
792 :
793 : /**
794 : * derivative of internal parameters with respect to stress parameters,
795 : * to avoid repeatedly allocating/deallocating this vector
796 : */
797 : mutable std::vector<std::vector<Real>> _dintnl_scratch;
798 :
799 : /**
800 : * jacobian in the NR process,
801 : * to avoid repeatedly allocating/deallocating this vector
802 : */
803 : mutable std::vector<double> _jac_scratch;
804 :
805 : /**
806 : * container to hold LAPACK ipiv integers,
807 : * to avoid repeatedly allocating/deallocating this vector
808 : */
809 : mutable std::vector<PetscBLASInt> _ipiv_scratch;
810 :
811 : /**
812 : * all the yield function information, for use in smoothAllQuantities,
813 : * to avoid repeatedly allocating/deallocating this vector
814 : */
815 : mutable std::vector<yieldAndFlow> _all_q_scratch;
816 :
817 : private:
818 : /**
819 : * The type of smoother function
820 : */
821 : enum class SmootherFunctionType
822 : {
823 : cos,
824 : poly1,
825 : poly2,
826 : poly3
827 : } _smoother_function_type;
828 : };
|
