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 "ComputeLagrangianStrain.h"
11 : #include "FactorizedRankTwoTensor.h"
12 : #include "MathUtils.h"
13 : #include "PermutationTensor.h"
14 :
15 : template <class G>
16 : InputParameters
17 5669 : ComputeLagrangianStrainBase<G>::baseParams()
18 : {
19 5669 : InputParameters params = Material::validParams();
20 :
21 11338 : params.addRequiredCoupledVar("displacements", "Displacement variables");
22 11338 : params.addParam<bool>(
23 11338 : "large_kinematics", false, "Use large displacement kinematics in the kernel.");
24 11338 : params.addParam<bool>("stabilize_strain", false, "Average the volumetric strains");
25 11338 : MooseEnum F_bar_mode("total incremental", "total");
26 11338 : params.addParam<MooseEnum>(
27 : "F_bar_mode",
28 : F_bar_mode,
29 : "What deformation gradient F-bar averages over (only used when `stabilize_strain = true`). "
30 : "'total' (default) averages the full F at each qp and rescales each qp's F by "
31 : "cbrt(det(F_avg)/det(F_ust)). 'incremental' averages the incremental F "
32 : "(F_ust * F_ust_old^{-1}) at each qp and rescales by cbrt(det(f_avg)/det(f_ust)); this is "
33 : "bit-for-bit compatible with the OLD `ComputeFiniteStrain` + `volumetric_locking_correction "
34 : "= "
35 : "true` formulation. Set to 'incremental' when cross-checking against the old kernel system.");
36 11338 : params.addParam<bool>(
37 : "publish_rotation_increment",
38 11338 : false,
39 : "If true, publish `rotation_increment = exp(vorticity_increment)` (Rodrigues) for "
40 : "downstream consumers that rotate by it (e.g. `ComputeMultiPlasticityStress` with "
41 : "`perform_finite_strain_rotations = true`). Default false keeps `rotation_increment = I` "
42 : "(the historical behavior -- the Lagrangian objective-rate machinery applies rotation "
43 : "externally). Enable when wrapping plasticity that needs its internal stress state to "
44 : "track the rotated Cauchy stress between steps, in tandem with `rotate_old_stress = true` "
45 : "on the objective rate.");
46 17007 : params.addRangeCheckedParam<Real>(
47 : "alpha",
48 11338 : 1.0,
49 : "alpha >= 0.5 & alpha <= 1.0",
50 : "Generalized midpoint weight for the deformation gradient. 1.0 = backward Euler (default), "
51 : "0.5 = midpoint rule (matches Abaqus/Implicit).");
52 11338 : MooseEnum kinematic_approximation("linear quadratic rashid_approximate rashid_eigen", "linear");
53 11338 : params.addParam<MooseEnum>(
54 : "kinematic_approximation",
55 : kinematic_approximation,
56 : "Approximation to the increment in the spatial velocity gradient: 'linear' (default; "
57 : "dL = I - f^{-1}), 'quadratic' (one more Taylor term), 'rashid_approximate' (Rashid's "
58 : "symmetric+skew formulas), or 'rashid_eigen' (exact log f via polar decomposition + "
59 : "matrix logs). Only affects large_kinematics; small kinematics is always linear.");
60 11338 : params.addParam<std::vector<MaterialPropertyName>>(
61 : "eigenstrain_names", {}, "List of eigenstrains to account for");
62 11338 : params.addParam<std::vector<MaterialPropertyName>>(
63 : "homogenization_gradient_names",
64 : {},
65 : "List of homogenization gradients to add to the displacement gradient");
66 :
67 11338 : params.addParam<std::string>("base_name", "Material property base name");
68 :
69 : // We rely on this *not* having use_displaced mesh on
70 5669 : params.suppressParameter<bool>("use_displaced_mesh");
71 :
72 5669 : return params;
73 5669 : }
74 :
75 : template <class G>
76 4251 : ComputeLagrangianStrainBase<G>::ComputeLagrangianStrainBase(const InputParameters & parameters)
77 : : Material(parameters),
78 : GuaranteeProvider(this),
79 4251 : _ndisp(coupledComponents("displacements")),
80 4251 : _disp(coupledValues("displacements")),
81 8502 : _grad_disp(coupledGradients("displacements")),
82 8550 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
83 8502 : _large_kinematics(getParam<bool>("large_kinematics")),
84 8502 : _stabilize_strain(getParam<bool>("stabilize_strain")),
85 8502 : _F_bar_mode(getParam<MooseEnum>("F_bar_mode").template getEnum<FBarMode>()),
86 8502 : _publish_rotation_increment(getParam<bool>("publish_rotation_increment")),
87 8502 : _alpha(getParam<Real>("alpha")),
88 4251 : _kinematic_approximation(
89 4251 : getParam<MooseEnum>("kinematic_approximation").template getEnum<KinematicApproximation>()),
90 8502 : _eigenstrain_names(getParam<std::vector<MaterialPropertyName>>("eigenstrain_names")),
91 4251 : _eigenstrains(_eigenstrain_names.size()),
92 4251 : _eigenstrains_old(_eigenstrain_names.size()),
93 4251 : _total_strain(declareProperty<RankTwoTensor>(_base_name + "total_strain")),
94 8502 : _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
95 4251 : _mechanical_strain(declareProperty<RankTwoTensor>(_base_name + "mechanical_strain")),
96 8502 : _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
97 4251 : _rotated_mechanical_strain(
98 4251 : declareProperty<RankTwoTensor>(_base_name + "rotated_mechanical_strain")),
99 4251 : _rotated_mechanical_strain_old(
100 4251 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "rotated_mechanical_strain")),
101 4251 : _strain_increment(declareProperty<RankTwoTensor>(_base_name + "strain_increment")),
102 4251 : _deformation_gradient_increment(
103 4251 : declareProperty<RankTwoTensor>(_base_name + "spatial_deformation_gradient_increment")),
104 4251 : _vorticity_increment(declareProperty<RankTwoTensor>(_base_name + "vorticity_increment")),
105 4251 : _F_ust(declareProperty<RankTwoTensor>(_base_name + "unstabilized_deformation_gradient")),
106 4251 : _F_ust_old(
107 4251 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "unstabilized_deformation_gradient")),
108 4251 : _F_actual(declareProperty<RankTwoTensor>(_base_name + "actual_deformation_gradient")),
109 4251 : _F_avg(declareProperty<RankTwoTensor>(_base_name + "average_deformation_gradient")),
110 4251 : _F(declareProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
111 8502 : _F_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
112 4251 : _F_inv(declareProperty<RankTwoTensor>(_base_name + "inverse_deformation_gradient")),
113 4251 : _f_inv(declareProperty<RankTwoTensor>(_base_name + "inverse_incremental_deformation_gradient")),
114 4251 : _F_ust_inv(
115 4251 : declareProperty<RankTwoTensor>(_base_name + "inverse_unstabilized_deformation_gradient")),
116 4251 : _F_ust_det(declareProperty<Real>(_base_name + "det_unstabilized_deformation_gradient")),
117 4251 : _d_deformation_gradient_increment_d_F(declareProperty<RankFourTensor>(
118 : _base_name + "d_spatial_deformation_gradient_increment_d_deformation_gradient")),
119 4251 : _d_vorticity_increment_d_F(declareProperty<RankFourTensor>(
120 : _base_name + "d_vorticity_increment_d_deformation_gradient")),
121 4251 : _d_F_d_grad_u(
122 4251 : declareProperty<RankFourTensor>(_base_name + "d_deformation_gradient_d_grad_displacement")),
123 4251 : _rotation(declareProperty<RankTwoTensor>(_base_name + "rotation")),
124 4251 : _stretch(declareProperty<RankTwoTensor>(_base_name + "stretch")),
125 4251 : _d_rotation_d_F(
126 4251 : declareProperty<RankFourTensor>(_base_name + "d_rotation_d_deformation_gradient")),
127 4251 : _d_F_stab_d_F_ust(declareProperty<RankFourTensor>(_base_name + "d_F_stab_d_F_unstabilized")),
128 4251 : _d_F_stab_d_F_avg(declareProperty<RankFourTensor>(_base_name + "d_F_stab_d_F_average")),
129 8502 : _homogenization_gradient_names(
130 : getParam<std::vector<MaterialPropertyName>>("homogenization_gradient_names")),
131 4251 : _homogenization_contributions(_homogenization_gradient_names.size()),
132 8502 : _rotation_increment(declareProperty<RankTwoTensor>(_base_name + "rotation_increment"))
133 : {
134 : // Couple old displacements only when the simulation is transient. With a Steady executioner
135 : // there is no "previous step", and the generalized midpoint rule treats the old state as
136 : // the undeformed reference (u_n = 0, grad u_n = 0, so F_n = I). The (1 - alpha) contribution
137 : // is then identically zero and we skip it in computeQpUnstabilizedDeformationGradient.
138 4251 : if (_fe_problem.isTransient())
139 : {
140 4155 : _disp_old = coupledValuesOld("displacements");
141 8310 : _grad_disp_old = coupledGradientsOld("displacements");
142 : }
143 :
144 : // Setup eigenstrains
145 4824 : for (auto i : make_range(_eigenstrain_names.size()))
146 : {
147 573 : _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_eigenstrain_names[i]);
148 573 : _eigenstrains_old[i] = &getMaterialPropertyOld<RankTwoTensor>(_eigenstrain_names[i]);
149 : }
150 :
151 : // In the future maybe there is a reason to have more than one, but for now
152 4251 : if (_homogenization_gradient_names.size() > 1)
153 0 : mooseError("ComputeLagrangianStrainBase cannot accommodate more than one "
154 : "homogenization gradient");
155 :
156 : // Setup homogenization contributions
157 4623 : for (unsigned int i = 0; i < _homogenization_gradient_names.size(); i++)
158 372 : _homogenization_contributions[i] =
159 372 : &getMaterialProperty<RankTwoTensor>(_homogenization_gradient_names[i]);
160 :
161 : // The strain calculator is the single source of truth for the kinematics regime. Publish a
162 : // LARGE_KINEMATICS guarantee on the deformation gradient (issued in the constructor so it is in
163 : // place before any consumer's initialSetup); the Lagrangian stress calculators and
164 : // stress-divergence kernels derive their own `large_kinematics` from it. Small kinematics leaves
165 : // the guarantee absent, which the consumers read as `large_kinematics = false`.
166 4251 : if (_large_kinematics)
167 8685 : issueGuarantee(_base_name + "deformation_gradient", Guarantee::LARGE_KINEMATICS);
168 4251 : }
169 :
170 : template <class G>
171 : void
172 676458 : ComputeLagrangianStrainBase<G>::initQpStatefulProperties()
173 : {
174 676458 : _total_strain[_qp].zero();
175 676458 : _mechanical_strain[_qp].zero();
176 676458 : _rotated_mechanical_strain[_qp].zero();
177 676458 : _F[_qp].setToIdentity();
178 676458 : _F_ust[_qp].setToIdentity();
179 676458 : _rotation[_qp].setToIdentity();
180 676458 : }
181 :
182 : template <class G>
183 : void
184 6838728 : ComputeLagrangianStrainBase<G>::computeProperties()
185 : {
186 : // Average the volumetric terms, if required
187 6838728 : computeDeformationGradient();
188 :
189 60194901 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
190 53356175 : computeQpProperties();
191 6838726 : }
192 :
193 : template <class G>
194 : void
195 53356175 : ComputeLagrangianStrainBase<G>::computeQpProperties()
196 : {
197 : // Add in the macroscale gradient contribution to both the stabilized `_F` (used by the
198 : // strain chain via `_f_inv`) AND the unstabilized `_F_ust` (used by the new F_ust-wrap
199 : // architecture in the stress materials). The homogenization gradient is a real
200 : // deformation imposed via the scalar constraint, not a stabilization -- it must appear
201 : // in F_ust too or PK1 = det(F_ust) sigma F_ust^{-T} will be missing its contribution.
202 66611395 : for (auto contribution : _homogenization_contributions)
203 : {
204 13255220 : _F[_qp] += (*contribution)[_qp];
205 13255220 : _F_ust[_qp] += (*contribution)[_qp];
206 : }
207 :
208 : // Publish F_ust^{-1} and det(F_ust) for the large-kinematics consumers' spatial push-forward
209 : // (grad_x = F_ust^{-T} grad_X, J_ust = det F_ust) and the Cauchy -> PK1 wrap. Computed once per
210 : // qp here -- shared by all displacement kernels and the stress calculator via the material
211 : // system -- instead of recomputed per test/trial/qp downstream. Gated on `_large_kinematics`
212 : // (every consumer reads these only on the large-kinematics path) and on `isPropertyActive` (so
213 : // consumers that never request them pay nothing).
214 55483058 : if (_large_kinematics && isPropertyActive(_F_ust_inv.id()))
215 : {
216 33624342 : _F_ust_inv[_qp] = _F_ust[_qp].inverse();
217 33624342 : _F_ust_det[_qp] = _F_ust[_qp].det();
218 : }
219 :
220 : // Skip the Jacobian-only RankFourTensor derivative chain when no downstream consumer
221 : // will read it (i.e. we're in a residual-only sweep). All of `_d_F_d_grad_u`,
222 : // `_d_deformation_gradient_increment_d_F`, `_d_vorticity_increment_d_F`, and the
223 : // `_d_rotation_d_F` slot of the polar decomposition feed only `*_jacobian` material
224 : // properties, which the kernel consumes only during Jacobian or
225 : // residual-and-Jacobian-together assembly.
226 53356175 : const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
227 : _fe_problem.currentlyComputingResidualAndJacobian();
228 :
229 : // dF/d(grad u_{n+1}) = alpha * I^{(4)} for the generalized midpoint rule
230 : // (alpha = 1.0 reduces to backward Euler).
231 : if (need_jacobian)
232 1813103 : _d_F_d_grad_u[_qp] = _alpha * RankFourTensor::IdentityFour();
233 :
234 53356175 : if (_large_kinematics)
235 : {
236 35751225 : _F_inv[_qp] = _F[_qp].inverse();
237 : // For `F_bar_mode = incremental`: `_f_inv` must invert the *incremental* F that was
238 : // F-bar'd (= `gamma_inc * F_ust * F_ust_old^{-1}`, matching OLD `ComputeFiniteStrain`'s
239 : // `_Fhat` after its volumetric-locking correction). Since `_F = gamma_inc * F_ust`,
240 : // `(gamma_inc * F_ust * F_ust_old^{-1})^{-1} = F_ust_old * _F^{-1}` -- i.e., pair the
241 : // unstabilized old F with the (cumulative) stabilized current `_F^{-1}`. Using
242 : // `_F_old` (the *cumulative* F-bar'd previous-step `_F`) here would compound F-bar
243 : // across steps and break OLD-compat. `total` mode keeps the existing form (where `_F`
244 : // and `_F_old` are the cumulative full-F F-bar'd values and the ratio gives the
245 : // F-bar'd incremental F directly).
246 35751225 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
247 917120 : _f_inv[_qp] = _F_ust_old[_qp] * _F_inv[_qp];
248 : else
249 34834105 : _f_inv[_qp] = _F_old[_qp] * _F_inv[_qp];
250 :
251 : // Dispatch to the active kinematic-approximation helper. Each helper returns
252 : // dd (= Deltad), dw (= Deltaw), and -- only when `need_jacobian` -- the
253 : // f^{-1}-derivatives of dL = dd + dw and of dw alone. The (dd, dw) outputs are needed
254 : // every iteration; the RankFour derivatives feed only the Jacobian chain below, so on
255 : // residual-only sweeps the helper skips them entirely.
256 35751225 : RankTwoTensor dd, dw;
257 35751225 : RankFourTensor d_dL_d_f_inv, d_dw_d_f_inv;
258 35751225 : computeQpLargeKinematicIncrement(
259 35751225 : _f_inv[_qp], dd, dw, d_dL_d_f_inv, d_dw_d_f_inv, need_jacobian);
260 :
261 35751223 : if (need_jacobian)
262 : {
263 : // Common chain rule: d(f^{-1})_{pq}/dF_{mn} = -f^{-1}_{pm} * F^{-1}_{nq}.
264 : usingTensorIndices(p_, q_, m_, n_);
265 1302080 : const RankFourTensor d_f_inv_d_F = -_f_inv[_qp].template times<p_, m_, n_, q_>(_F_inv[_qp]);
266 1302080 : _d_deformation_gradient_increment_d_F[_qp] = d_dL_d_f_inv * d_f_inv_d_F;
267 1302080 : _d_vorticity_increment_d_F[_qp] = d_dw_d_f_inv * d_f_inv_d_F;
268 : }
269 :
270 35751223 : setQpIncrementalStrains(dd, dw);
271 : // The polar decomposition (a tensor eigensolve + sqrt per qp, on every residual and Jacobian
272 : // eval) feeds only the Green-Naghdi objective rate's `_rotation` / `_d_rotation_d_F` (and its
273 : // internal `_stretch`). Skip it entirely when no active consumer needs it -- Truesdell uses
274 : // the deformation-gradient increment, Jaumann the vorticity increment, Rashid the Rodrigues
275 : // `_rotation_increment`. The active-property set self-corrects for future consumers with no
276 : // coupling flag.
277 104094321 : if (isPropertyActive(_rotation.id()) || isPropertyActive(_stretch.id()))
278 3159348 : computeQpPolarDecomposition(need_jacobian);
279 : }
280 : // For small deformations we just provide the identity (and always linear)
281 : else
282 : {
283 17604950 : _F_inv[_qp] = RankTwoTensor::Identity();
284 17604950 : _f_inv[_qp] = RankTwoTensor::Identity();
285 17604950 : const RankTwoTensor dL = _F[_qp] - _F_old[_qp];
286 :
287 17604950 : if (need_jacobian)
288 : {
289 : // d(dL)/dF = I^{(4)} when dL = F - F_old. d(dW)/dF = the skew projector.
290 511023 : _d_deformation_gradient_increment_d_F[_qp] = RankFourTensor::IdentityFour();
291 : usingTensorIndices(i_, j_, k_, l_);
292 : const auto I2 = RankTwoTensor::Identity();
293 511023 : _d_vorticity_increment_d_F[_qp] =
294 1022046 : 0.5 * (RankFourTensor::IdentityFour() - I2.template times<j_, k_, i_, l_>(I2));
295 :
296 : // Small kinematics: R = I, U = I, dR/dF = 0. Defensive defaults; GN is not used here.
297 511023 : _d_rotation_d_F[_qp].zero();
298 : }
299 17604950 : _rotation[_qp].setToIdentity();
300 17604950 : _stretch[_qp].setToIdentity();
301 :
302 17604950 : setQpIncrementalStrains(0.5 * (dL + dL.transpose()), 0.5 * (dL - dL.transpose()));
303 : }
304 53356173 : }
305 :
306 : template <class G>
307 : void
308 3159348 : ComputeLagrangianStrainBase<G>::computeQpPolarDecomposition(bool need_jacobian)
309 : {
310 : // Polar decomposition F = R * U of the alpha-weighted, F-bar-stabilized deformation
311 : // gradient at this qp. We decompose _F rather than _F_actual because the rest of the
312 : // kernel/rate chain treats _F as the spatial frame; using _F here matches the pre-3.1
313 : // GN rate exactly (it read _def_grad, which is _F via ComputeLagrangianStressCauchy).
314 3159348 : const RankTwoTensor & F = _F[_qp];
315 3159348 : FactorizedRankTwoTensor C(F.transpose() * F);
316 : // Reuse the sqrt factorization (a tensor eigensolve) for both U and U^{-1}.
317 3159348 : const auto sqrt_C = MathUtils::sqrt(C);
318 3159348 : _stretch[_qp] = sqrt_C.get();
319 3159348 : const RankTwoTensor U_inv = sqrt_C.inverse().get();
320 3159348 : _rotation[_qp] = F * U_inv;
321 :
322 3159348 : if (!need_jacobian)
323 : return;
324 :
325 : // dR/dF closed form. See ComputeLagrangianObjectiveStress.C:221-227.
326 : const auto I = RankTwoTensor::Identity();
327 16396 : const RankTwoTensor Y = _stretch[_qp].trace() * I - _stretch[_qp];
328 16396 : const RankTwoTensor Z = _rotation[_qp] * Y;
329 16396 : const RankTwoTensor O = Z * _rotation[_qp].transpose();
330 : usingTensorIndices(i_, j_, k_, l_);
331 16396 : _d_rotation_d_F[_qp] =
332 32792 : (O.template times<i_, k_, l_, j_>(Y) - Z.template times<i_, l_, k_, j_>(Z)) / Y.det();
333 : }
334 :
335 : template <class G>
336 : void
337 35751225 : ComputeLagrangianStrainBase<G>::computeQpLargeKinematicIncrement(const RankTwoTensor & f_inv,
338 : RankTwoTensor & dd,
339 : RankTwoTensor & dw,
340 : RankFourTensor & d_dL_d_f_inv,
341 : RankFourTensor & d_dw_d_f_inv,
342 : bool need_jacobian)
343 : {
344 35751225 : switch (_kinematic_approximation)
345 : {
346 21187823 : case KinematicApproximation::Linear:
347 21187823 : computeLinearIncrement(f_inv, dd, dw, d_dL_d_f_inv, d_dw_d_f_inv, need_jacobian);
348 21187823 : break;
349 4791488 : case KinematicApproximation::Quadratic:
350 4791488 : computeQuadraticIncrement(f_inv, dd, dw, d_dL_d_f_inv, d_dw_d_f_inv, need_jacobian);
351 4791488 : break;
352 4878784 : case KinematicApproximation::RashidApproximate:
353 4878784 : computeRashidApproximateIncrement(f_inv, dd, dw, d_dL_d_f_inv, d_dw_d_f_inv, need_jacobian);
354 4878784 : break;
355 4893130 : case KinematicApproximation::RashidEigen:
356 4893130 : computeRashidEigenIncrement(f_inv, dd, dw, d_dL_d_f_inv, d_dw_d_f_inv, need_jacobian);
357 4893128 : break;
358 : }
359 35751223 : }
360 :
361 : template <class G>
362 : void
363 0 : ComputeLagrangianStrainBase<G>::computeQpIncrementalStrains(const RankTwoTensor & dL)
364 : {
365 : // Backward-compatible entry point: split into sym/skew and delegate.
366 0 : setQpIncrementalStrains(0.5 * (dL + dL.transpose()), 0.5 * (dL - dL.transpose()));
367 0 : }
368 :
369 : template <class G>
370 : void
371 53356173 : ComputeLagrangianStrainBase<G>::setQpIncrementalStrains(const RankTwoTensor & dd,
372 : const RankTwoTensor & dw)
373 : {
374 53356173 : _strain_increment[_qp] = dd;
375 53356173 : _vorticity_increment[_qp] = dw;
376 : // Full kinematic spatial velocity gradient increment, before any eigenstrain subtraction.
377 : // The objective-rate advection in ComputeLagrangianObjectiveStress consumes this.
378 53356173 : _deformation_gradient_increment[_qp] = dd + dw;
379 :
380 : // Increment the total strain
381 53356173 : _total_strain[_qp] = _total_strain_old[_qp] + _strain_increment[_qp];
382 :
383 : // Get rid of the eigenstrains
384 : // Note we currently do not alter the deformation gradient -- this will be
385 : // needed in the future for a "complete" system
386 53356173 : subtractQpEigenstrainIncrement(_strain_increment[_qp]);
387 :
388 : // Increment the mechanical strain
389 53356173 : _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
390 :
391 : // Additionally maintain the rotated mechanical-strain accumulator,
392 : // eps_n+1 = r_hat (eps_n + Deltad) r_hat^T, where r_hat = exp(Deltaw) via Rodrigues. This matches
393 : // the mechanical_strain output convention of `ComputeFiniteStrain` and is consumed only by aux
394 : // variables -- the stress chain still uses the un-rotated `_mechanical_strain` above. In small
395 : // kinematics dw is small and r_hat ~= I + dw + 1/2 dw^2, which contributes only second-order
396 : // corrections (equivalent to no rotation for small strain).
397 53356173 : RankTwoTensor r_hat;
398 53356173 : if (_large_kinematics)
399 : {
400 35751223 : const Real theta2 = 0.5 * dw.doubleContraction(dw);
401 35751223 : const Real theta = std::sqrt(theta2);
402 : Real f, g;
403 : const Real small_theta = 1.0e-7;
404 35751223 : if (theta < small_theta)
405 : {
406 2853513 : f = 1.0 - theta2 / 6.0;
407 2853513 : g = 0.5 - theta2 / 24.0;
408 : }
409 : else
410 : {
411 32897710 : f = std::sin(theta) / theta;
412 32897710 : g = (1.0 - std::cos(theta)) / theta2;
413 : }
414 35751223 : r_hat = RankTwoTensor::Identity() + f * dw + g * dw * dw;
415 : }
416 : else
417 : {
418 17604950 : r_hat = RankTwoTensor::Identity();
419 : }
420 53356173 : _rotated_mechanical_strain[_qp] =
421 53356173 : r_hat * (_rotated_mechanical_strain_old[_qp] + _strain_increment[_qp]) * r_hat.transpose();
422 :
423 : // Published rotation increment for downstream `ComputeStressBase`-style materials. The
424 : // default (identity) preserves the historical Lagrangian pipeline where the objective
425 : // rate is the sole rotation source. With `publish_rotation_increment = true` we dispatch
426 : // to the `kinematic_approximation`-matched formula so wrapped plasticity with
427 : // `perform_finite_strain_rotations = true` rotates its `_stress` bit-for-bit like OLD's
428 : // `ComputeFiniteStrain` would have done (the rate then runs in `rotate_old_stress`
429 : // passthrough mode so we don't double-rotate). Specifically `rashid_approximate` uses
430 : // OLD's C1/C2/C3 polynomial form (not `exp(dw)`) so the wrapped material's accumulated
431 : // stress storage matches OLD's `_stress` byte-for-byte through return mapping; without
432 : // this, the small (~1e-5) per-step rotation drift between `exp(dw)` and OLD's R_incr
433 : // amplifies through plastic flow into ~1e-3 cumulative stress error.
434 103168 : _rotation_increment[_qp] = (_publish_rotation_increment && _large_kinematics)
435 53459341 : ? computeQpRotationIncrement(_f_inv[_qp], dw)
436 : : RankTwoTensor::Identity();
437 53356173 : }
438 :
439 : template <class G>
440 : void
441 21187823 : ComputeLagrangianStrainBase<G>::computeLinearIncrement(const RankTwoTensor & f_inv,
442 : RankTwoTensor & dd,
443 : RankTwoTensor & dw,
444 : RankFourTensor & d_dL_d_f_inv,
445 : RankFourTensor & d_dw_d_f_inv,
446 : bool need_jacobian) const
447 : {
448 : // Deltal = I - f^{-1}, so Deltad = sym(Deltal), Deltaw = skew(Deltal).
449 21187823 : const RankTwoTensor dL = RankTwoTensor::Identity() - f_inv;
450 21187823 : dd = 0.5 * (dL + dL.transpose());
451 21187823 : dw = 0.5 * (dL - dL.transpose());
452 :
453 21187823 : if (!need_jacobian)
454 : return;
455 : // d(Deltal)/d(f^{-1}) = -I^{(4)}.
456 1153792 : d_dL_d_f_inv = -RankFourTensor::IdentityFour();
457 : // d(Deltaw)/d(f^{-1}) = - skew projector on f^{-1} = -(1/2)(I^{ikjl} - I^{iljk}).
458 : usingTensorIndices(i_, j_, m_, n_);
459 : const auto I2 = RankTwoTensor::Identity();
460 3461376 : d_dw_d_f_inv = -0.5 * (RankFourTensor::IdentityFour() - I2.template times<j_, m_, i_, n_>(I2));
461 : }
462 :
463 : template <class G>
464 : void
465 4791488 : ComputeLagrangianStrainBase<G>::computeQuadraticIncrement(const RankTwoTensor & f_inv,
466 : RankTwoTensor & dd,
467 : RankTwoTensor & dw,
468 : RankFourTensor & d_dL_d_f_inv,
469 : RankFourTensor & d_dw_d_f_inv,
470 : bool need_jacobian) const
471 : {
472 : // Deltal = X + (1/2) X^2 with X = I - f^{-1} (one more Taylor term of -log f^{-1}).
473 4791488 : const RankTwoTensor X = RankTwoTensor::Identity() - f_inv;
474 4791488 : const RankTwoTensor dL = X + 0.5 * X * X;
475 4791488 : dd = 0.5 * (dL + dL.transpose());
476 4791488 : dw = 0.5 * (dL - dL.transpose());
477 :
478 4791488 : if (!need_jacobian)
479 : return;
480 : // dX/d(f^{-1}) = -I^{(4)}, and d(X^2)_{ij}/dX_{mn} = delta_{im} X_{nj} + X_{im} delta_{jn}.
481 : // So d(Deltal)/d(f^{-1}) = -I^{(4)} - (1/2) (delta_{im} X_{nj} + X_{im} delta_{jn}).
482 : usingTensorIndices(i_, j_, m_, n_);
483 : const auto I2 = RankTwoTensor::Identity();
484 : const RankFourTensor dXX_dX =
485 28224 : I2.template times<i_, m_, n_, j_>(X) + X.template times<i_, m_, j_, n_>(I2);
486 28224 : d_dL_d_f_inv = -RankFourTensor::IdentityFour() - 0.5 * dXX_dX;
487 : // dw = (1/2)(dL - dL^T), so d(dw)_{ij}/d... = (1/2)(d(dL)_{ij}/d... - d(dL)_{ji}/d...).
488 56448 : d_dw_d_f_inv = 0.5 * (d_dL_d_f_inv - d_dL_d_f_inv.transposeIj());
489 : }
490 :
491 : template <class G>
492 : void
493 4878784 : ComputeLagrangianStrainBase<G>::computeRashidApproximateIncrement(const RankTwoTensor & f_inv,
494 : RankTwoTensor & dd,
495 : RankTwoTensor & dw,
496 : RankFourTensor & d_dL_d_f_inv,
497 : RankFourTensor & d_dw_d_f_inv,
498 : bool need_jacobian) const
499 : {
500 : // See plan_outline.pdf Sec.2.3 (eq 10-15, with the corrected vorticity).
501 : // X = I - f^{-1}. Symmetric part: A = X X^T - X - X^T, Deltad = -A/2 + A^2/4.
502 : // Skew part (from the rotation tensor): alpha_i = eps_ijk (f^{-1})_jk,
503 : // cos theta = (tr(f^{-1}) - 1)/2, sin theta = sqrt(1 - cos^2theta), Q = (1/4) alpha*alpha,
504 : // Deltaw_ij = -(theta / (2sqrtQ)) eps_ijk alpha_k.
505 : usingTensorIndices(i_, j_, m_, n_);
506 : const auto I2 = RankTwoTensor::Identity();
507 : const auto I4 = RankFourTensor::IdentityFour();
508 :
509 : // ---- symmetric part ----
510 4878784 : const RankTwoTensor X = I2 - f_inv;
511 4878784 : const RankTwoTensor Xt = X.transpose();
512 4878784 : const RankTwoTensor A = X * Xt - X - Xt;
513 4878784 : dd = -0.5 * A + 0.25 * A * A;
514 :
515 : // ---- skew part ----
516 : // alpha_i = eps_ijk (f^{-1})_jk (axial vector of f^{-1}'s skew part, doubled).
517 : RealVectorValue alpha;
518 19515136 : for (unsigned int i = 0; i < 3; ++i)
519 : {
520 : Real ai = 0.0;
521 58545408 : for (unsigned int j = 0; j < 3; ++j)
522 175636224 : for (unsigned int k = 0; k < 3; ++k)
523 131727168 : ai += PermutationTensor::eps(i, j, k) * f_inv(j, k);
524 14636352 : alpha(i) = ai;
525 : }
526 : // Derive theta from the axial vector's magnitude rather than from tr(f^{-1}): the trace
527 : // formula assumes f^{-1} is exactly a rotation, which is only true in the limit. Here
528 : // sin theta = sqrtQ (Q = |skew part|^2/4); a true rotation matches both, and an arbitrary f^{-1}
529 : // is projected onto its nearest "rotation-like" interpretation.
530 4878784 : const Real Q_raw = 0.25 * (alpha * alpha);
531 : // Clamp Q just below 1 so cos theta stays strictly positive (theta -> pi/2 is unphysical for one
532 : // step).
533 4878784 : const Real Q = std::min(Q_raw, 1.0 - 1.0e-12);
534 :
535 : // Small-angle fallback: when sin theta -> 0, theta/(2sqrtQ) -> 1/2 (L'Hopital) and Deltaw ->
536 : // -alpha/2, which matches sym/skew of the linear approximation. d(Deltaw)/d(f^{-1}) is the
537 : // antisymmetrizer (1/2)(delta_im delta_jn - delta_jm delta_in).
538 : const Real small_Q = 1.0e-12;
539 4878784 : if (Q < small_Q)
540 : {
541 : // Deltaw_ij = -(1/2) eps_ijk alpha_k
542 267952 : for (unsigned int i = 0; i < 3; ++i)
543 803856 : for (unsigned int j = 0; j < 3; ++j)
544 : {
545 : Real v = 0.0;
546 2411568 : for (unsigned int k = 0; k < 3; ++k)
547 1808676 : v += PermutationTensor::eps(i, j, k) * alpha(k);
548 602892 : dw(i, j) = -0.5 * v;
549 : }
550 : // d(Deltaw)/d(f^{-1}) = -(1/2)(delta_{im} delta_{jn} - delta_{jm} delta_{in})
551 66988 : if (need_jacobian)
552 12904 : d_dw_d_f_inv = -0.5 * (I4 - I2.template times<j_, m_, i_, n_>(I2));
553 : }
554 : else
555 : {
556 4811796 : const Real sin_theta = std::sqrt(Q);
557 4811796 : const Real cos_theta = std::sqrt(1.0 - Q);
558 4811796 : const Real theta = std::asin(sin_theta);
559 4811796 : const Real coeff = -theta / (2.0 * sin_theta);
560 :
561 19247184 : for (unsigned int i = 0; i < 3; ++i)
562 57741552 : for (unsigned int j = 0; j < 3; ++j)
563 : {
564 : Real v = 0.0;
565 173224656 : for (unsigned int k = 0; k < 3; ++k)
566 129918492 : v += PermutationTensor::eps(i, j, k) * alpha(k);
567 43306164 : dw(i, j) = coeff * v;
568 : }
569 :
570 : // Build d(Deltaw)/d(f^{-1}) analytically.
571 : // Deltaw_ij = c(f^{-1}) * E_ij(f^{-1}), E_ij = eps_ijk alpha_k, c = -theta/(2 sin theta).
572 : //
573 : // dalpha_k/d(f^{-1})_{mn} = eps_{kmn} (Levi-Civita is constant).
574 : // dQ/d(f^{-1})_{mn} = (1/2) alpha_k eps_{kmn}.
575 : // With sin theta = sqrtQ : dsin theta/d(f^{-1})_{mn} = alpha_k eps_{kmn} / (4 sin theta).
576 : // dtheta/d(f^{-1})_{mn} = dsin theta/d(f^{-1})_{mn} / cos theta.
577 : // dc/dtheta = -(sin theta - theta cos theta)/(2 sin^2 theta); dc/d(f^{-1})_{mn} = dc/dtheta *
578 : // dtheta/d(f^{-1})_{mn}
579 : // = (theta cos theta - sin theta) alpha_k eps_{kmn} / (8 sin^3 theta cos theta).
580 : // dE_ij/d(f^{-1})_{mn} = eps_{ijk} eps_{kmn}.
581 : //
582 : // Final:
583 : // d(Deltaw)_{ij}/d(f^{-1})_{mn} = (dc/d(f^{-1})_{mn}) * E_ij + c * dE_ij/d(f^{-1})_{mn}.
584 4811796 : if (need_jacobian)
585 : {
586 53516 : const Real dc_pref =
587 53516 : (theta * cos_theta - sin_theta) / (8.0 * sin_theta * sin_theta * sin_theta * cos_theta);
588 214064 : for (unsigned int i = 0; i < 3; ++i)
589 642192 : for (unsigned int j = 0; j < 3; ++j)
590 : {
591 : Real E_ij = 0.0;
592 1926576 : for (unsigned int k = 0; k < 3; ++k)
593 1444932 : E_ij += PermutationTensor::eps(i, j, k) * alpha(k);
594 1926576 : for (unsigned int m = 0; m < 3; ++m)
595 5779728 : for (unsigned int n = 0; n < 3; ++n)
596 : {
597 : Real eps_alpha = 0.0;
598 17339184 : for (unsigned int k = 0; k < 3; ++k)
599 13004388 : eps_alpha += PermutationTensor::eps(k, m, n) * alpha(k);
600 4334796 : const Real dc_dfinv = dc_pref * eps_alpha;
601 : Real dE_dfinv = 0.0;
602 17339184 : for (unsigned int k = 0; k < 3; ++k)
603 13004388 : dE_dfinv += PermutationTensor::eps(i, j, k) * PermutationTensor::eps(k, m, n);
604 4334796 : d_dw_d_f_inv(i, j, m, n) = dc_dfinv * E_ij + coeff * dE_dfinv;
605 : }
606 : }
607 : }
608 : }
609 :
610 59968 : if (!need_jacobian)
611 : return;
612 :
613 : // ---- symmetric-part derivative + assemble d(Deltal)/d(f^{-1}) (Jacobian only) ----
614 : // dX/d(f^{-1}) = -I^{(4)}.
615 : // d(X X^T)_{ij}/dX_{mn} = delta_{im} X_{jn} + X_{in} delta_{jm} (from (X X^T)_{ij} = X_{ik}
616 : // X_{jk}). d(X^T)_{ij}/dX_{mn} = delta_{jm} delta_{in}.
617 : // -> d(A)/d(f^{-1}) = -[d(X X^T)/dX] + I^{(4)} + (transposed I^{(4)}).
618 : const RankFourTensor d_XXt_dX =
619 59968 : I2.template times<i_, m_, j_, n_>(X) + X.template times<i_, n_, j_, m_>(I2);
620 59968 : const RankFourTensor d_Xt_dX = I2.template times<j_, m_, i_, n_>(I2);
621 59968 : const RankFourTensor dA_dfinv = -d_XXt_dX + I4 + d_Xt_dX;
622 :
623 : // d(A^2)_{ij}/dA_{mn} = delta_{im} A_{nj} + A_{im} delta_{jn}.
624 : const RankFourTensor d_AA_dA =
625 59968 : I2.template times<i_, m_, n_, j_>(A) + A.template times<i_, m_, j_, n_>(I2);
626 59968 : const RankFourTensor d_AA_dfinv = d_AA_dA * dA_dfinv;
627 :
628 59968 : const RankFourTensor d_dd_dfinv = -0.5 * dA_dfinv + 0.25 * d_AA_dfinv;
629 59968 : d_dL_d_f_inv = d_dd_dfinv + d_dw_d_f_inv;
630 : }
631 :
632 : template <class G>
633 : void
634 4893130 : ComputeLagrangianStrainBase<G>::computeRashidEigenIncrement(const RankTwoTensor & f_inv,
635 : RankTwoTensor & dd,
636 : RankTwoTensor & dw,
637 : RankFourTensor & d_dL_d_f_inv,
638 : RankFourTensor & d_dw_d_f_inv,
639 : bool need_jacobian) const
640 : {
641 : // See plan_outline.pdf Sec.2.4. Polar-decompose f^{-1} = r' u', with u' symmetric positive
642 : // definite and r' a proper rotation. The PDF's identity `log f^{-1} = -log d - log w`
643 : // only holds when log u and log r commute (i.e. for non-rotating deformation); in general
644 : // the right stretch of f^{-1} is u' = R * U^{-1} * R^T where R, U are the right polar of f.
645 : // So -log(u') = R * log U * R^T is the *spatial-frame* log strain, not the co-rotated
646 : // log U. We compute it that way first, then rotate by r' = R^T to land in the n-frame
647 : // so that _strain_increment = log U exactly. This is the strain measure the Rashid
648 : // objective stress update consumes (eq. 22, sigma_{n+1} = r_hat (sigma_n + Deltasigma) r_hat^T
649 : // expects Deltasigma in the n-frame).
650 : usingTensorIndices(a_, b_, m_, n_);
651 : const auto I2 = RankTwoTensor::Identity();
652 :
653 : // ---- polar decomposition of f^{-1} ----
654 4893130 : FactorizedRankTwoTensor cprime(f_inv.transpose() * f_inv);
655 : // Reuse the sqrt factorization (a tensor eigensolve) for both u and u^{-1}.
656 4893128 : const auto sqrt_cprime = MathUtils::sqrt(cprime);
657 : const RankTwoTensor u = sqrt_cprime.get();
658 4893128 : const RankTwoTensor u_inv = sqrt_cprime.inverse().get();
659 4893128 : const RankTwoTensor r = f_inv * u_inv;
660 :
661 : // ---- intermediate Deltad_spatial = -log u' = R * log U * R^T (n+1 frame) ----
662 : // u' = sqrt(c'), so log(u') = (1/2) log(c'). Reuse c''s already-computed factorization
663 : // (`sqrt_cprime` shares c''s eigenvectors) instead of re-decomposing u -- one fewer 3x3
664 : // eigensolve per qp, on both residual and Jacobian sweeps.
665 9786256 : const RankTwoTensor dd_spatial = -0.5 * MathUtils::log(cprime).get();
666 :
667 : // ---- Deltaw = -log r via Rodrigues. ----
668 : // log r = phi(theta) (r - r^T) with phi(theta) = theta/(2 sin theta), cos theta = (tr r - 1)/2.
669 4893128 : const Real cos_theta = MathUtils::clamp(0.5 * (r.trace() - 1.0), -1.0, 1.0);
670 4893128 : const Real sin2 = std::max(1.0 - cos_theta * cos_theta, 0.0);
671 4893128 : const Real sin_theta = std::sqrt(sin2);
672 4893128 : const Real theta = std::acos(cos_theta);
673 4893128 : const RankTwoTensor A = r - r.transpose();
674 :
675 : // d(log r)_{ij}/d(r)_{mn} = (dphi/dr_{mn}) A_{ij} + phi (delta_{im} delta_{jn} - delta_{jm}
676 : // delta_{in}).
677 : // dphi/dr_{mn} = (dphi/dtheta)(dtheta/d cos theta)(d cos theta/dr_{mn})
678 : // = psi delta_{mn}, where psi = (theta cos theta - sin theta)/(4 sin^3 theta).
679 : // Small-angle: phi -> 1/2, psi -> -1/12, so d(log r)/dr -> (1/2)(I^(4) - swap_ij).
680 : // `d_logr_d_r` feeds only the Jacobian, so it is built only when `need_jacobian`.
681 : const Real small_sin = 1.0e-7;
682 4893128 : RankFourTensor d_logr_d_r;
683 4893128 : RankTwoTensor log_r;
684 : usingTensorIndices(i_, j_, k_, l_);
685 4893128 : if (std::abs(sin_theta) < small_sin)
686 : {
687 66554 : log_r = 0.5 * A;
688 66554 : if (need_jacobian)
689 : {
690 6438 : const RankFourTensor swap_ij = I2.template times<j_, m_, i_, n_>(I2);
691 12876 : d_logr_d_r = 0.5 * (RankFourTensor::IdentityFour() - swap_ij);
692 : }
693 : }
694 : else
695 : {
696 4826574 : const Real phi = theta / (2.0 * sin_theta);
697 4826574 : log_r = phi * A;
698 4826574 : if (need_jacobian)
699 : {
700 53658 : const Real psi = (theta * cos_theta - sin_theta) / (4.0 * sin_theta * sin2);
701 53658 : const RankFourTensor swap_ij = I2.template times<j_, m_, i_, n_>(I2);
702 : // (dphi/dr)_{mn} = psi delta_{mn} -> outer with A_{ij} gives A.times<i_, j_, m_, n_>(I2) *
703 : // psi.
704 53658 : const RankFourTensor dphi_outer_A = A.template times<i_, j_, m_, n_>(I2);
705 107316 : d_logr_d_r = psi * dphi_outer_A + phi * (RankFourTensor::IdentityFour() - swap_ij);
706 : }
707 : }
708 4893128 : dw = -log_r;
709 :
710 : // ---- Rotate Deltad from spatial back to co-rotated (n) frame: log U = r' * (R log U R^T) *
711 : // r'^T, using r' = R^T from the polar of f^{-1}. ----
712 4893128 : dd = r * dd_spatial * r.transpose();
713 :
714 : // Everything below is the Jacobian chain (RankFour matrix-log derivatives, the
715 : // polar-decomposition dr/d(f^{-1}) closed form, and the three-piece sandwich chain rule).
716 : // Skip it wholesale on residual-only sweeps.
717 4893128 : if (!need_jacobian)
718 : return;
719 :
720 : // d(log u)/d(c') = (1/2) dlog(c'), and d(c')_{ab}/d(f^{-1})_{mn} = delta_{an} f_inv_{mb}
721 : // + delta_{bn} f_inv_{ma}.
722 60096 : const RankFourTensor dlog_cprime = MathUtils::dlog(cprime);
723 : const RankFourTensor d_cprime_d_finv =
724 60096 : I2.template times<a_, n_, m_, b_>(f_inv) + I2.template times<b_, n_, m_, a_>(f_inv);
725 60096 : const RankFourTensor d_dd_spatial_d_finv = -0.5 * (dlog_cprime * d_cprime_d_finv);
726 :
727 : // ---- d(r)/d(f^{-1}) via the polar-decomposition closed form (same trick as
728 : // ComputeLagrangianObjectiveStress::polarDecomposition). ----
729 60096 : const RankTwoTensor Y = u.trace() * I2 - u;
730 60096 : const RankTwoTensor Z = r * Y;
731 60096 : const RankTwoTensor O = Z * r.transpose();
732 : const RankFourTensor d_r_d_finv =
733 60096 : (O.template times<i_, k_, l_, j_>(Y) - Z.template times<i_, l_, k_, j_>(Z)) / Y.det();
734 :
735 60096 : d_dw_d_f_inv = -(d_logr_d_r * d_r_d_finv);
736 :
737 : // The derivative of the sandwich r'*A*r'^T w.r.t. f^{-1} expands to three rank-4 pieces
738 : // (chain rule on r' AND on A). The 6-argument `times<>(RankFourTensor)` overload uses
739 : // x[0..4] (one dummy at index 4), so we reuse the same dummy label `p2_` in each
740 : // sub-contraction.
741 : {
742 : usingTensorIndices(i2_, j2_, m2_, n2_, p2_);
743 60096 : const RankTwoTensor M = dd_spatial * r.transpose(); // dd * r^T (p, j) shape
744 60096 : const RankTwoTensor N = r * dd_spatial; // r * dd (i, q) shape
745 : // T1_{ijmn} = (dr/d_finv)_{ip,mn} * M_{pj}
746 60096 : const RankFourTensor T1 = M.template times<p2_, j2_, i2_, p2_, m2_, n2_>(d_r_d_finv);
747 : // T2_{ijmn} = r_{ip} * (d_dd_spatial_d_finv)_{pq,mn} * r_{jq}: contract twice over the dummy.
748 60096 : const RankFourTensor mid = r.template times<i2_, p2_, p2_, j2_, m2_, n2_>(d_dd_spatial_d_finv);
749 60096 : const RankFourTensor T2 = r.template times<j2_, p2_, i2_, p2_, m2_, n2_>(mid);
750 : // T3_{ijmn} = N_{iq} * (dr/d_finv)_{jq,mn}
751 60096 : const RankFourTensor T3 = N.template times<i2_, p2_, j2_, p2_, m2_, n2_>(d_r_d_finv);
752 60096 : const RankFourTensor d_dd_d_finv = T1 + T2 + T3;
753 60096 : d_dL_d_f_inv = d_dd_d_finv + d_dw_d_f_inv;
754 : }
755 : }
756 :
757 : template <class G>
758 : void
759 53356173 : ComputeLagrangianStrainBase<G>::subtractQpEigenstrainIncrement(RankTwoTensor & strain)
760 : {
761 54017065 : for (auto i : make_range(_eigenstrain_names.size()))
762 660892 : strain -= (*_eigenstrains[i])[_qp] - (*_eigenstrains_old[i])[_qp];
763 53356173 : }
764 :
765 : template <class G>
766 : void
767 53356181 : ComputeLagrangianStrainBase<G>::computeQpUnstabilizedDeformationGradient()
768 : {
769 : // Generalized midpoint: F^alpha_{n+1} = I + alpha * (grad u_{n+1}, u_{n+1})
770 : // + (1 - alpha) * (grad u_n, u_n).
771 : // alpha = 1.0 reduces to backward Euler (no old contribution). With a Steady executioner
772 : // the old displacement is treated as identically zero (F_n = I), so we skip the old call
773 : // entirely - this lets a user run with alpha != 1 in steady mode as well.
774 53356181 : _F_ust[_qp].setToIdentity();
775 53356181 : const bool include_old = _alpha != 1.0 && _fe_problem.isTransient();
776 212132634 : for (auto component : make_range(_ndisp))
777 : {
778 159508841 : G::addGradOp(_F_ust[_qp],
779 : component,
780 158776453 : _alpha * (*_grad_disp[component])[_qp],
781 158776453 : _alpha * (*_disp[component])[_qp],
782 781516 : _q_point[_qp]);
783 158776453 : if (include_old)
784 28453848 : G::addGradOp(_F_ust[_qp],
785 : component,
786 28473816 : (1.0 - _alpha) * (*_grad_disp_old[component])[_qp],
787 28473816 : (1.0 - _alpha) * (*_disp_old[component])[_qp],
788 : _q_point[_qp]);
789 : }
790 53356181 : }
791 :
792 : template <class G>
793 : void
794 53356181 : ComputeLagrangianStrainBase<G>::computeQpActualDeformationGradient()
795 : {
796 : // The literal deformation gradient at n+1 (no alpha weighting, no F-bar). For alpha = 1
797 : // this is identical to _F_ust.
798 53356181 : _F_actual[_qp].setToIdentity();
799 212132634 : for (auto component : make_range(_ndisp))
800 158776453 : G::addGradOp(_F_actual[_qp],
801 : component,
802 : (*_grad_disp[component])[_qp],
803 158776453 : (*_disp[component])[_qp],
804 781516 : _q_point[_qp]);
805 53356181 : }
806 :
807 : template <class G>
808 : void
809 6838728 : ComputeLagrangianStrainBase<G>::computeDeformationGradient()
810 : {
811 : // First calculate the unstabilized deformation gradient at each qp
812 60194909 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
813 : {
814 53356181 : computeQpUnstabilizedDeformationGradient();
815 53356181 : computeQpActualDeformationGradient();
816 53356181 : _F[_qp] = _F_ust[_qp];
817 : }
818 :
819 : usingTensorIndices(i_, j_, k_, l_);
820 : const auto I2 = RankTwoTensor::Identity();
821 :
822 : // The F-bar local/non-local derivative material properties feed only the Jacobian path
823 : // (the stress materials chain them into `_pk1_jacobian`/`_cauchy_jacobian`, and the
824 : // TL kernel's non-local F-bar term contracts `_d_F_stab_d_F_avg`). The residual sweep
825 : // never reads them, so skip the R4 algebra when assembling a residual alone.
826 6838728 : const bool need_jacobian = _fe_problem.currentlyComputingJacobian() ||
827 : _fe_problem.currentlyComputingResidualAndJacobian();
828 :
829 : // If stabilization is on do the volumetric correction
830 6838728 : if (_stabilize_strain)
831 : {
832 : // `_F_avg` consistently stores the element average of `F_ust` regardless of mode;
833 : // it's consumed by the UL kernel for the spatial-frame push-forward, which is a
834 : // purely-geometric quantity. In `F_bar_mode = incremental` the F-bar chain itself
835 : // operates on the *incremental* averaged tensor `f_avg = avg(F_ust * F_ust_old^{-1})`,
836 : // which we compute locally -- the kernel matches by weighting grad_phi by F_ust_old^{-1}
837 : // in its element average. Incremental mode only makes sense for large kinematics
838 : // (OLD's Fhat F-bar lives in the finite-strain code path).
839 1865038 : const bool incremental = (_F_bar_mode == FBarMode::Incremental);
840 1865038 : if (incremental && !_large_kinematics)
841 0 : mooseError("`F_bar_mode = incremental` requires `large_kinematics = true`. The "
842 : "incremental F-bar formulation is the multiplicative correction to the "
843 : "incremental F, which is only defined for large kinematics. Use "
844 : "`F_bar_mode = total` (the default) with small kinematics.");
845 1865038 : const auto F_avg = StabilizationUtils::elementAverage(
846 16482378 : [this](unsigned int qp) { return _F_ust[qp]; }, _JxW, _coord);
847 1750398 : const auto f_avg =
848 : incremental
849 2782158 : ? StabilizationUtils::elementAverage([this](unsigned int qp)
850 917120 : { return _F_ust[qp] * _F_ust_old[qp].inverse(); },
851 : _JxW,
852 : _coord)
853 : : RankTwoTensor();
854 : // Always publish avg(F_ust) -- UL kernel uses this for the push-forward.
855 1865038 : _F_avg.set().setAllValues(F_avg);
856 : // What the F-bar gamma and `_d_F_stab_d_F_avg` are built against (also what the
857 : // kernel-side `_avg_grad_trial` will represent the delta of):
858 1865038 : const auto & avg_for_chain = incremental ? f_avg : F_avg;
859 : // Make the appropriate modification, depending on small or large deformations
860 16482378 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
861 : {
862 14617340 : if (_large_kinematics)
863 : {
864 : // Multiplicative F-bar: F_stab = gamma * F_ust where
865 : // gamma_total = cbrt(det(F_avg) / det(F_ust)) (total mode)
866 : // gamma_inc = cbrt(det(f_avg) / det(f_ust)), f_ust = F_ust*F_ust_old^{-1}
867 : // (incremental mode)
868 : // For incremental: det(f_ust) = det(F_ust)/det(F_ust_old).
869 : // d log det(f_ust)/d F_ust = F_ust^{-T} (same as the total-mode chain)
870 : // so `_d_F_stab_d_F_ust` shares the total-mode shape. The non-local chain
871 : // contracts `_d_F_stab_d_F_avg` with delta(averaged-quantity); kernel computes
872 : // deltaf_avg in incremental mode by weighting grad_phi by F_ust_old^{-1} in its
873 : // element average.
874 : const Real det_ust_local =
875 7829836 : incremental ? _F_ust[_qp].det() / _F_ust_old[_qp].det() : _F[_qp].det();
876 7829836 : const Real gamma = std::pow(avg_for_chain.det() / det_ust_local, 1.0 / 3.0);
877 7829836 : if (need_jacobian)
878 : {
879 122328 : const auto Fust_invT = _F_ust[_qp].inverse().transpose();
880 122328 : const auto avg_invT = avg_for_chain.inverse().transpose();
881 122328 : _d_F_stab_d_F_ust[_qp] =
882 122328 : gamma * RankFourTensor::IdentityFour() -
883 122328 : (gamma / 3.0) * _F_ust[_qp].template times<i_, j_, k_, l_>(Fust_invT);
884 122328 : _d_F_stab_d_F_avg[_qp] =
885 244656 : (gamma / 3.0) * _F_ust[_qp].template times<i_, j_, k_, l_>(avg_invT);
886 : }
887 7829836 : _F[_qp] *= gamma;
888 : }
889 : else
890 : {
891 6787504 : if (need_jacobian)
892 : {
893 : // Additive (trace) F-bar: F_stab = F_ust + (tr(F_avg - F_ust)/3) * I.
894 : // dF_stab/dF_ust = I^(4) - (1/3) * I2 (x) I2 (each diagonal component pulled out).
895 : // dF_stab/dF_avg = (1/3) * I2 (x) I2.
896 168624 : const auto outer = I2.template times<i_, j_, k_, l_>(I2);
897 168624 : _d_F_stab_d_F_ust[_qp] = RankFourTensor::IdentityFour() - (1.0 / 3.0) * outer;
898 168624 : _d_F_stab_d_F_avg[_qp] = (1.0 / 3.0) * outer;
899 : }
900 6787504 : _F[_qp] += (F_avg.trace() - _F[_qp].trace()) * I2 / 3.0;
901 : }
902 : }
903 : }
904 4973690 : else if (need_jacobian)
905 : {
906 : // F-bar off: dF_stab/dF_ust = I^(4), dF_stab/dF_avg = 0.
907 1736981 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
908 : {
909 1522151 : _d_F_stab_d_F_ust[_qp] = RankFourTensor::IdentityFour();
910 1522151 : _d_F_stab_d_F_avg[_qp].zero();
911 : }
912 : }
913 6838728 : }
914 :
915 : template <class G>
916 : RankTwoTensor
917 103168 : ComputeLagrangianStrainBase<G>::computeQpRotationIncrement(const RankTwoTensor & f_inv,
918 : const RankTwoTensor & dw) const
919 : {
920 : // For `rashid_approximate` port OLD `ComputeFiniteStrain`'s C1/C2/C3 polynomial form
921 : // exactly (Rashid 1993). Pairing this rotation with the wrapped material's FSR (via
922 : // `_rotation_increment`) lets the constitutive's `_stress` evolution match OLD's bit-
923 : // for-bit at the same converged displacement state.
924 103168 : if (_kinematic_approximation == KinematicApproximation::RashidApproximate)
925 : {
926 51584 : const Real a[3] = {
927 51584 : f_inv(1, 2) - f_inv(2, 1), f_inv(2, 0) - f_inv(0, 2), f_inv(0, 1) - f_inv(1, 0)};
928 51584 : const Real q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
929 51584 : const Real trFhatinv_1 = f_inv.trace() - 1.0;
930 51584 : const Real p = trFhatinv_1 * trFhatinv_1 / 4.0;
931 51584 : const Real C1_squared = p +
932 51584 : 3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
933 51584 : 2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q);
934 51584 : if (C1_squared <= 0.0)
935 0 : mooseException(
936 : "Cannot take square root of a number less than or equal to zero in the calculation of "
937 : "C1 for the Rashid approximation for the rotation tensor.");
938 51584 : const Real C1 = std::sqrt(C1_squared);
939 : Real C2;
940 51584 : if (q > 0.01)
941 0 : C2 = (1.0 - C1) / (4.0 * q);
942 : else
943 51584 : C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
944 51584 : Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
945 : Utility::pow<3>(p) +
946 51584 : Utility::pow<3>(q) *
947 51584 : (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
948 51584 : 5.0 * Utility::pow<4>(p)) /
949 51584 : (512.0 * Utility::pow<4>(p));
950 51584 : const Real C3_test =
951 51584 : (p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) / Utility::pow<3>(p + q);
952 51584 : if (C3_test <= 0.0)
953 0 : mooseException(
954 : "Cannot take square root of a number less than or equal to zero in the calculation of "
955 : "C3_test for the Rashid approximation for the rotation tensor.");
956 51584 : const Real C3 = 0.5 * std::sqrt(C3_test);
957 51584 : RankTwoTensor R_incr;
958 51584 : R_incr.addIa(C1);
959 206336 : for (unsigned int i = 0; i < 3; ++i)
960 619008 : for (unsigned int j = 0; j < 3; ++j)
961 464256 : R_incr(i, j) += C2 * a[i] * a[j];
962 51584 : R_incr(0, 1) += C3 * a[2];
963 51584 : R_incr(0, 2) -= C3 * a[1];
964 51584 : R_incr(1, 0) -= C3 * a[2];
965 51584 : R_incr(1, 2) += C3 * a[0];
966 51584 : R_incr(2, 0) += C3 * a[1];
967 51584 : R_incr(2, 1) -= C3 * a[0];
968 51584 : return R_incr.transpose();
969 : }
970 :
971 : // For `rashid_eigen`, `linear`, `quadratic`: r_hat = exp(dw) via Rodrigues. For RashidEigen,
972 : // dw is the matrix log of the polar-decomposition R, so exp(dw) recovers that R bit-for-bit
973 : // -- equivalent to OLD `ComputeFiniteStrain`'s EigenSolution rotation.
974 51584 : const Real theta2 = 0.5 * dw.doubleContraction(dw);
975 51584 : const Real theta = std::sqrt(theta2);
976 : Real f, g;
977 : const Real small_theta = 1.0e-7;
978 51584 : if (theta < small_theta)
979 : {
980 3968 : f = 1.0 - theta2 / 6.0;
981 3968 : g = 0.5 - theta2 / 24.0;
982 : }
983 : else
984 : {
985 47616 : f = std::sin(theta) / theta;
986 47616 : g = (1.0 - std::cos(theta)) / theta2;
987 : }
988 51584 : return RankTwoTensor::Identity() + f * dw + g * dw * dw;
989 : }
990 :
991 : template class ComputeLagrangianStrainBase<GradientOperatorCartesian>;
992 : template class ComputeLagrangianStrainBase<GradientOperatorAxisymmetricCylindrical>;
993 : template class ComputeLagrangianStrainBase<GradientOperatorCentrosymmetricSpherical>;
|
