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 "TotalLagrangianStressDivergenceBase.h"
11 :
12 : template <class G>
13 : InputParameters
14 5586 : TotalLagrangianStressDivergenceBase<G>::baseParams()
15 : {
16 5586 : InputParameters params = LagrangianStressDivergenceBase::validParams();
17 : // This kernel requires use_displaced_mesh to be off
18 5586 : params.suppressParameter<bool>("use_displaced_mesh");
19 5586 : return params;
20 0 : }
21 :
22 : template <class G>
23 2793 : TotalLagrangianStressDivergenceBase<G>::TotalLagrangianStressDivergenceBase(
24 : const InputParameters & parameters)
25 : : LagrangianStressDivergenceBase(parameters),
26 2793 : _pk1(getMaterialPropertyByName<RankTwoTensor>(_base_name + "pk1_stress")),
27 5586 : _dpk1(getMaterialPropertyByName<RankFourTensor>(_base_name + "pk1_jacobian")),
28 5586 : _dpk1_d_grad_u(getMaterialPropertyByName<RankFourTensor>(_base_name + "dpk1_d_grad_u")),
29 2793 : _dpk1_bypass_fbar(
30 8379 : getMaterialPropertyByName<RankFourTensor>(_base_name + "pk1_jacobian_bypass_fbar"))
31 : {
32 2793 : }
33 :
34 : template <class G>
35 : RankTwoTensor
36 1440805824 : TotalLagrangianStressDivergenceBase<G>::gradTest(unsigned int component)
37 : {
38 : // F-bar doesn't modify the test function
39 1440805824 : return G::gradOp(component, _grad_test[_i][_qp], _test[_i][_qp], _q_point[_qp]);
40 : }
41 :
42 : template <class G>
43 : RankTwoTensor
44 23673216 : TotalLagrangianStressDivergenceBase<G>::gradTrial(unsigned int component)
45 : {
46 : // After the F_ust-wrap architectural change, pk1_jacobian = dPK1/d(F_ust) already
47 : // contains the local F-bar contribution (via the sigma-chain through
48 : // `_d_F_stab_d_F_ust`). The trial gradient is always unstabilized; the non-local F-bar
49 : // Jacobian contribution is added explicitly in `computeQpJacobianDisplacement`.
50 23673216 : return gradTrialUnstabilized(component);
51 : }
52 :
53 : template <class G>
54 : RankTwoTensor
55 56989152 : TotalLagrangianStressDivergenceBase<G>::gradTrialUnstabilized(unsigned int component)
56 : {
57 : // Without F-bar stabilization, simply return the gradient of the trial functions
58 56989152 : return G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
59 : }
60 :
61 : template <class G>
62 : void
63 188471 : TotalLagrangianStressDivergenceBase<G>::precalculateJacobianDisplacement(unsigned int component)
64 : {
65 : // For total Lagrangian, the averaging is taken on the reference frame regardless of geometric
66 : // nonlinearity. The stored quantity is what the strain calc's `_d_F_stab_d_F_avg` contracts
67 : // with: deltaF_avg in `F_bar_mode = total`, deltaf_avg in `F_bar_mode = incremental`.
68 : // total: avg_q(grad_phi[j][q]) -> deltaF_avg per DOF
69 : // incremental: avg_q(grad_phi[j][q] * F_ust_old[q]^{-1})
70 : // -> deltaf_avg per DOF
71 : // (since f_ust = F_ust * F_ust_old^{-1} and F_ust_old is fixed w.r.t. the
72 : // current Newton iterate's disp, so deltaf_ust = deltaF_ust * F_ust_old^{-1}).
73 188471 : const bool incremental = (_F_bar_mode == FBarMode::Incremental);
74 : // F_ust_old[qp]^{-1} depends only on qp, but the average below runs once per trial DOF j.
75 : // Precompute the per-qp inverse once instead of re-inverting the 3x3 for every (qp, j).
76 : std::vector<RankTwoTensor> f_ust_old_inv;
77 188471 : if (incremental)
78 : {
79 60336 : f_ust_old_inv.resize(_qrule->n_points());
80 543024 : for (const auto qp : make_range(_qrule->n_points()))
81 482688 : f_ust_old_inv[qp] = _F_ust_old[qp].inverse();
82 : }
83 1593199 : for (auto j : make_range(_phi.size()))
84 1404728 : _avg_grad_trial[component][j] = StabilizationUtils::elementAverage(
85 12245944 : [this, component, j, incremental, &f_ust_old_inv](unsigned int qp)
86 : {
87 10841216 : const RankTwoTensor g = G::gradOp(component, _grad_phi[j][qp], _phi[j][qp], _q_point[qp]);
88 10841216 : return incremental ? g * f_ust_old_inv[qp] : g;
89 : },
90 : _JxW,
91 : _coord);
92 188471 : }
93 :
94 : template <class G>
95 : void
96 12515274 : TotalLagrangianStressDivergenceBase<G>::precalculateResidual()
97 : {
98 : LagrangianStressDivergenceBase::precalculateResidual();
99 12515274 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
100 : {
101 314880 : computeAverageGradientSpatialTest();
102 314880 : cachePK1ContractionFUst();
103 : }
104 12515274 : }
105 :
106 : template <class G>
107 : void
108 375216 : TotalLagrangianStressDivergenceBase<G>::cachePK1ContractionFUst()
109 : {
110 375216 : const unsigned int n_qp = _qrule->n_points();
111 375216 : _pk1_ddot_F_ust.assign(n_qp, 0.0);
112 3376944 : for (const auto qp : make_range(n_qp))
113 3001728 : _pk1_ddot_F_ust[qp] = _pk1[qp].doubleContraction(_F_ust[qp]);
114 375216 : }
115 :
116 : template <class G>
117 : void
118 478748 : TotalLagrangianStressDivergenceBase<G>::precalculateJacobian()
119 : {
120 478748 : LagrangianStressDivergenceBase::precalculateJacobian();
121 : // Populate the per-(qp, j) local PK1 chain cache for the diagonal column (beta = _alpha).
122 : // Runs regardless of `_stabilize_strain` because the local chain is always needed; the
123 : // helper itself only touches `_avg_grad_trial` when stab is on.
124 478748 : populateLocalPK1Cache(_alpha);
125 478748 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
126 : {
127 20112 : computeAverageGradientSpatialTest();
128 20112 : computeAverageGradientSpatialPhi(_alpha);
129 20112 : computeAvgTestPhiCross(_alpha);
130 : }
131 478748 : }
132 :
133 : template <class G>
134 : void
135 825766 : TotalLagrangianStressDivergenceBase<G>::precalculateOffDiagJacobian(unsigned int jvar)
136 : {
137 825766 : LagrangianStressDivergenceBase::precalculateOffDiagJacobian(jvar);
138 : // Populate the per-(qp, j) local PK1 chain cache for the matched off-diagonal beta.
139 1624108 : for (auto beta : make_range(_ndisp))
140 1620002 : if (jvar == _disp_nums[beta])
141 : {
142 821660 : populateLocalPK1Cache(beta);
143 821660 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
144 : {
145 40224 : computeAverageGradientSpatialTest();
146 40224 : computeAverageGradientSpatialPhi(beta);
147 40224 : computeAvgTestPhiCross(beta);
148 : }
149 : break;
150 : }
151 825766 : }
152 :
153 : template <class G>
154 : void
155 1300408 : TotalLagrangianStressDivergenceBase<G>::populateLocalPK1Cache(unsigned int beta)
156 : {
157 1300408 : const unsigned int n_qp = _qrule->n_points();
158 1300408 : const unsigned int n_phi = _phi.size();
159 1300408 : _dpk1_grad_trial_cache.assign(n_qp, std::vector<RankTwoTensor>(n_phi));
160 2600816 : _grad_trial_cache.assign(n_qp, std::vector<RankTwoTensor>(n_phi));
161 1300408 : if (_stabilize_strain)
162 : {
163 364186 : _delta_PK1_NL_cache.assign(n_qp, std::vector<RankTwoTensor>(n_phi));
164 364186 : _dpk1_total_cache.assign(n_qp, std::vector<RankTwoTensor>(n_phi));
165 : }
166 :
167 : // B-bar volumetric correction (incremental mode only): cache the column-independent
168 : // `pk1 : F_ust` per qp and the per-(qp, j) `dA/dU` term so the Jacobian inner loop reads
169 : // scalars instead of recomputing three double contractions per (_i, _j, component).
170 1300408 : const bool bbar = _stabilize_strain && _F_bar_mode == FBarMode::Incremental;
171 : if (bbar)
172 : {
173 60336 : cachePK1ContractionFUst();
174 120672 : _dA_dU_cache.assign(n_qp, std::vector<Real>(n_phi));
175 : }
176 :
177 11098468 : for (unsigned int qp = 0; qp < n_qp; ++qp)
178 : {
179 9798060 : const RankFourTensor & dpk1 = _dpk1_d_grad_u[qp];
180 9798060 : const RankFourTensor & dFdGU = _d_F_d_grad_u[qp];
181 85826264 : for (unsigned int j = 0; j < n_phi; ++j)
182 : {
183 : // gradTrial == gradTrialUnstabilized in TL (line ~50); single cache feeds both consumers.
184 76028204 : const RankTwoTensor grad_trial = G::gradOp(beta, _grad_phi[j][qp], _phi[j][qp], _q_point[qp]);
185 76028204 : _grad_trial_cache[qp][j] = grad_trial;
186 76028204 : _dpk1_grad_trial_cache[qp][j] = dpk1 * grad_trial;
187 76028204 : if (_stabilize_strain)
188 : {
189 : // Compose the full non-local F-bar PK1 perturbation:
190 : // deltasigma_NL = D_nl[qp] * (_d_F_d_grad_u[qp] * _avg_grad_trial[beta][j])
191 : // deltaPK1_NL = det(F_ust) * deltasigma_NL * F_ust^{-T} (large kinematics)
192 : // = deltasigma_NL (small kinematics)
193 : // D_nl (`_d_nl_fbar`), det(F_ust), and F_ust^{-1} are all per-qp material properties
194 : // from the stress/strain materials -- no kernel-side cache.
195 10656224 : const RankTwoTensor delta_F_avg = dFdGU * _avg_grad_trial[beta][j];
196 10656224 : const RankTwoTensor delta_sigma_nl = (*_d_nl_fbar)[qp] * delta_F_avg;
197 10656224 : if (_large_kinematics)
198 5362112 : _delta_PK1_NL_cache[qp][j] =
199 5362112 : (*_F_ust_det)[qp] * delta_sigma_nl * (*_F_ust_inv)[qp].transpose();
200 : else
201 5294112 : _delta_PK1_NL_cache[qp][j] = delta_sigma_nl;
202 : // Local + non-local PK1 derivative: the single tensor the disp Jacobian contracts against.
203 10656224 : _dpk1_total_cache[qp][j] = _dpk1_grad_trial_cache[qp][j] + _delta_PK1_NL_cache[qp][j];
204 : }
205 : // dA/dU = [ (dPK1_local + dPK1_NL) : F_ust + pk1 : dF_ust ] / 3, all per (qp, j).
206 76028204 : if (bbar)
207 3861504 : _dA_dU_cache[qp][j] = (_dpk1_total_cache[qp][j].doubleContraction(_F_ust[qp]) +
208 3861504 : _pk1[qp].doubleContraction(_grad_trial_cache[qp][j])) /
209 : 3.0;
210 : }
211 : }
212 1300408 : }
213 :
214 : template <class G>
215 : Real
216 81936384 : TotalLagrangianStressDivergenceBase<G>::gradXTestComponent(unsigned int component) const
217 : {
218 : // (grad_x test_i)_component = (grad_X test_i)_j (F_ust^{-1})_{j, component}
219 : // Push-forward via F_ust (= F_actual at alpha=1) matches the PK1 wrap. F_ust^{-1} is a
220 : // material property (computed once per qp by the strain calc), not recomputed here.
221 81936384 : if (!_large_kinematics)
222 0 : return _grad_test[_i][_qp](component);
223 81936384 : const RankTwoTensor & F_inv = (*_F_ust_inv)[_qp];
224 : Real out = 0.0;
225 327745536 : for (unsigned int j = 0; j < 3; ++j)
226 245809152 : out += _grad_test[_i][_qp](j) * F_inv(j, component);
227 : return out;
228 : }
229 :
230 : template <class G>
231 : Real
232 30892032 : TotalLagrangianStressDivergenceBase<G>::gradXPhiComponent(unsigned int component) const
233 : {
234 30892032 : if (!_large_kinematics)
235 0 : return _grad_phi[_j][_qp](component);
236 30892032 : const RankTwoTensor & F_inv = (*_F_ust_inv)[_qp];
237 : Real out = 0.0;
238 123568128 : for (unsigned int j = 0; j < 3; ++j)
239 92676096 : out += _grad_phi[_j][_qp](j) * F_inv(j, component);
240 : return out;
241 : }
242 :
243 : template <class G>
244 : void
245 375216 : TotalLagrangianStressDivergenceBase<G>::computeAverageGradientSpatialTest()
246 : {
247 : // _avg_grad_spatial_test[i] = (1/V_x) int (grad_x test_i)_alpha dV_x
248 : // = (sum_qp (grad_X test_i)_j (F_ust^{-1})_{j, _alpha} * J_ust * w)
249 : // / (sum_qp J_ust * w), w = JxW * coord
250 : // For small kinematics, F_ust = I, J_ust = 1, reducing to OLD's element-averaged grad_test.
251 375216 : _avg_grad_spatial_test.assign(_test.size(), 0.0);
252 : Real V_x = 0.0;
253 3376944 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
254 : {
255 3001728 : const Real J = _large_kinematics ? (*_F_ust_det)[qp] : 1.0;
256 3001728 : const Real w = J * _JxW[qp] * _coord[qp];
257 3001728 : V_x += w;
258 3001728 : const RankTwoTensor F_inv = _large_kinematics ? (*_F_ust_inv)[qp] : RankTwoTensor::Identity();
259 27015552 : for (unsigned int i = 0; i < _test.size(); ++i)
260 : {
261 : Real g_x = 0.0;
262 96055296 : for (unsigned int j = 0; j < 3; ++j)
263 72041472 : g_x += _grad_test[i][qp](j) * F_inv(j, _alpha);
264 24013824 : _avg_grad_spatial_test[i] += g_x * w;
265 : }
266 : }
267 3376944 : for (unsigned int i = 0; i < _test.size(); ++i)
268 3001728 : _avg_grad_spatial_test[i] /= V_x;
269 375216 : }
270 :
271 : template <class G>
272 : void
273 60336 : TotalLagrangianStressDivergenceBase<G>::computeAverageGradientSpatialPhi(unsigned int beta)
274 : {
275 : // Same structure as computeAverageGradientSpatialTest, but for trial functions and component
276 : // beta.
277 60336 : _avg_grad_spatial_phi[beta].assign(_phi.size(), 0.0);
278 : Real V_x = 0.0;
279 543024 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
280 : {
281 482688 : const Real J = _large_kinematics ? (*_F_ust_det)[qp] : 1.0;
282 482688 : const Real w = J * _JxW[qp] * _coord[qp];
283 482688 : V_x += w;
284 482688 : const RankTwoTensor F_inv = _large_kinematics ? (*_F_ust_inv)[qp] : RankTwoTensor::Identity();
285 4344192 : for (unsigned int j = 0; j < _phi.size(); ++j)
286 : {
287 : Real g_x = 0.0;
288 15446016 : for (unsigned int k = 0; k < 3; ++k)
289 11584512 : g_x += _grad_phi[j][qp](k) * F_inv(k, beta);
290 3861504 : _avg_grad_spatial_phi[beta][j] += g_x * w;
291 : }
292 : }
293 543024 : for (unsigned int j = 0; j < _phi.size(); ++j)
294 482688 : _avg_grad_spatial_phi[beta][j] /= V_x;
295 60336 : }
296 :
297 : template <class G>
298 : void
299 60336 : TotalLagrangianStressDivergenceBase<G>::computeAvgTestPhiCross(unsigned int beta)
300 : {
301 : // Populate the cross-products
302 : // _avg_test_phi_cross[b1][b2][i][j] = (1/V_x) int (grad_x test_i)_{b1} * (grad_x phi_j)_{b2}
303 : // dV_x
304 : // for the 4 combinations of (b1, b2) in {_alpha, beta}. The d(avg_grad_spatial_test)/dU
305 : // term in the B-bar Jacobian needs both (alpha, beta) and (beta, alpha) entries; the
306 : // diagonal (alpha == beta) collapses to a single combination.
307 60336 : const unsigned int n_test = _test.size();
308 60336 : const unsigned int n_phi = _phi.size();
309 60336 : std::vector<unsigned int> bs = {_alpha};
310 60336 : if (beta != _alpha)
311 40224 : bs.push_back(beta);
312 : // Zero out and resize the affected slots
313 160896 : for (auto b1 : bs)
314 281568 : for (auto b2 : bs)
315 : {
316 181008 : _avg_test_phi_cross[b1][b2].assign(n_test, std::vector<Real>(n_phi, 0.0));
317 : }
318 : Real V_x = 0.0;
319 : // Cache (grad_x test_i)_{b1} and (grad_x phi_j)_{b2} per qp to avoid recomputing for
320 : // each (i, j) pair.
321 60336 : std::vector<std::vector<Real>> gx_test(bs.size(), std::vector<Real>(n_test, 0.0));
322 60336 : std::vector<std::vector<Real>> gx_phi(bs.size(), std::vector<Real>(n_phi, 0.0));
323 543024 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
324 : {
325 482688 : const Real J = _large_kinematics ? (*_F_ust_det)[qp] : 1.0;
326 482688 : const Real w = J * _JxW[qp] * _coord[qp];
327 482688 : V_x += w;
328 482688 : const RankTwoTensor F_inv = _large_kinematics ? (*_F_ust_inv)[qp] : RankTwoTensor::Identity();
329 1287168 : for (unsigned int bi = 0; bi < bs.size(); ++bi)
330 : {
331 804480 : const unsigned int b = bs[bi];
332 7240320 : for (unsigned int i = 0; i < n_test; ++i)
333 : {
334 : Real g = 0.0;
335 25743360 : for (unsigned int k = 0; k < 3; ++k)
336 19307520 : g += _grad_test[i][qp](k) * F_inv(k, b);
337 6435840 : gx_test[bi][i] = g;
338 : }
339 7240320 : for (unsigned int j = 0; j < n_phi; ++j)
340 : {
341 : Real g = 0.0;
342 25743360 : for (unsigned int k = 0; k < 3; ++k)
343 19307520 : g += _grad_phi[j][qp](k) * F_inv(k, b);
344 6435840 : gx_phi[bi][j] = g;
345 : }
346 : }
347 1287168 : for (unsigned int b1i = 0; b1i < bs.size(); ++b1i)
348 2252544 : for (unsigned int b2i = 0; b2i < bs.size(); ++b2i)
349 : {
350 1448064 : const unsigned int b1 = bs[b1i];
351 1448064 : const unsigned int b2 = bs[b2i];
352 1448064 : auto & dest = _avg_test_phi_cross[b1][b2];
353 13032576 : for (unsigned int i = 0; i < n_test; ++i)
354 104260608 : for (unsigned int j = 0; j < n_phi; ++j)
355 92676096 : dest[i][j] += gx_test[b1i][i] * gx_phi[b2i][j] * w;
356 : }
357 : }
358 160896 : for (auto b1 : bs)
359 281568 : for (auto b2 : bs)
360 : {
361 181008 : auto & dest = _avg_test_phi_cross[b1][b2];
362 1629072 : for (unsigned int i = 0; i < n_test; ++i)
363 13032576 : for (unsigned int j = 0; j < n_phi; ++j)
364 11584512 : dest[i][j] /= V_x;
365 : }
366 60336 : }
367 :
368 : template <class G>
369 : Real
370 784124980 : TotalLagrangianStressDivergenceBase<G>::computeQpResidual()
371 : {
372 784124980 : Real result = gradTest(_alpha).doubleContraction(_pk1[_qp]);
373 :
374 : // OLD-compat B-bar volumetric correction (only in `F_bar_mode = incremental`, the explicit
375 : // OLD-compat mode). This adds
376 : // int (tr sigma / 3) * (avg_grad_x test - grad_x test)_alpha dV_x
377 : // expressed in reference-frame integration as
378 : // sum_qp ((PK1 : F_ust) / 3) * (avg_grad_spatial_test[_i] - grad_x_test_alpha) * JxW * coord
379 : // (the J = det(F_ust) factor is already absorbed into `PK1 : F_ust = J * tr sigma`).
380 : // Without this term, PK1 = det(F_ust) sigma F_ust^{-T} reproduces OLD's `sigma : grad_x test`
381 : // term but not OLD's vol-locking correction; for uniform sigma that's harmless (the integral
382 : // vanishes per element) but for non-uniform sigma -- e.g. mixed-component Dirichlet BCs producing
383 : // non-affine equilibrium with plasticity + a kinematic decomposition that gives spatially varying
384 : // sigma -- it drives the converged displacement away from OLD's. Scoped to incremental mode
385 : // because (a) that's the documented OLD-compat F-bar mode and (b) the existing `total` mode tests
386 : // are calibrated to the F-bar-in-sigma-only formulation (adding the B-bar term there would be
387 : // a separate behavior change).
388 784124980 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
389 : {
390 20152320 : const Real PK1_F_over_3 = _pk1_ddot_F_ust[_qp] / 3.0;
391 20152320 : result += PK1_F_over_3 * (_avg_grad_spatial_test[_i] - gradXTestComponent(_alpha));
392 : }
393 784124980 : return result;
394 : }
395 :
396 : template <class G>
397 : Real
398 598848620 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianDisplacement(unsigned int alpha,
399 : unsigned int beta)
400 : {
401 : // gradTest(alpha) feeds the local Jacobian contraction below.
402 598848620 : const RankTwoTensor grad_test = gradTest(alpha);
403 :
404 : // Local (+ non-local F-bar) Jacobian: J_{alpha beta} = gradTest_alpha : dPK1_total, where
405 : // dPK1_total = dPK1/d(grad u) * grad_phi_beta (local) plus the wrapped non-local F-bar
406 : // perturbation. Both pieces depend only on (qp, j) and are precomputed and pre-summed into
407 : // `_dpk1_total_cache` by `populateLocalPK1Cache(beta)`, so a single contraction covers both.
408 : // Without stabilization there is no non-local piece, so the local cache is contracted directly.
409 1113750104 : Real J = grad_test.doubleContraction(_stabilize_strain ? _dpk1_total_cache[_qp][_j]
410 514901484 : : _dpk1_grad_trial_cache[_qp][_j]);
411 :
412 : // OLD-compat B-bar volumetric correction Jacobian (only in `F_bar_mode = incremental`).
413 : // The residual added the term R_extra = (PK1:F_ust)/3 * (avg_T - T_alpha) summed over qps.
414 : // Linearizing per (i, j, alpha, beta):
415 : // J_extra_qp = (1/3) * d(PK1:F_ust)/dU * (avg_T - T_alpha)
416 : // + (1/3) * (PK1:F_ust) * (d(avg_T)/dU - dT_alpha/dU)
417 : // with derivatives in U (the trial-DOF for component beta at node j):
418 : // d(PK1:F_ust)/dU|qp = (dPK1_local + dPK1_NL) : F_ust + PK1 : dF_ust_local
419 : // dPK1_local = _dpk1_grad_trial_cache[_qp][_j] (cached above)
420 : // dPK1_NL = d_PK1_NL (cached above)
421 : // dF_ust = _grad_trial_cache[_qp][_j] (cached unstabilized trial)
422 : // dT_alpha/dU|qp = -(grad_x test_i)_beta * (grad_x phi_j)_alpha (local cross)
423 : // d(avg_T)/dU = avgCross(alpha, beta) - avgCross(beta, alpha) - avg_T(i) *
424 : // avg_grad_spatial_phi[beta][j]
425 : // avgCross(b1, b2) = _avg_test_phi_cross[b1][b2][i][j] (precomputed)
426 598848620 : if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
427 : {
428 : // A_qp (= pk1:F_ust/3, per qp) and dA_dU (per qp, j) are precomputed in
429 : // `cachePK1ContractionFUst` / `populateLocalPK1Cache` -- no double contractions here.
430 30892032 : const Real A_qp = _pk1_ddot_F_ust[_qp] / 3.0;
431 30892032 : const Real avg_T = _avg_grad_spatial_test[_i];
432 30892032 : const Real T_alpha = gradXTestComponent(alpha);
433 30892032 : const Real B = avg_T - T_alpha;
434 :
435 30892032 : const Real dA_dU = _dA_dU_cache[_qp][_j];
436 :
437 : // dB/dU: -dT_alpha/dU + d(avg_T)/dU
438 : // dT_alpha/dU = -(grad_x test_i)_beta * (grad_x phi_j)_alpha (local cross at _qp)
439 30892032 : const Real T_beta_i = gradXTestComponent(beta);
440 30892032 : const Real phi_alpha_j = gradXPhiComponent(alpha);
441 30892032 : const Real dT_alpha_dU = -T_beta_i * phi_alpha_j;
442 :
443 30892032 : const Real avg_cross_ab = _avg_test_phi_cross[alpha][beta][_i][_j];
444 30892032 : const Real avg_cross_ba = _avg_test_phi_cross[beta][alpha][_i][_j];
445 30892032 : const Real d_avg_T_dU = avg_cross_ab - avg_cross_ba - avg_T * _avg_grad_spatial_phi[beta][_j];
446 :
447 30892032 : const Real dB_dU = d_avg_T_dU - dT_alpha_dU;
448 30892032 : J += dA_dU * B + A_qp * dB_dU;
449 : }
450 :
451 598848620 : return J;
452 : }
453 :
454 : template <class G>
455 : Real
456 843072 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianTemperature(unsigned int cvar)
457 : {
458 : // Multiple eigenstrains may depend on the same coupled var
459 843072 : RankTwoTensor total_deigen;
460 1686144 : for (const auto deigen_darg : _deigenstrain_dargs[cvar])
461 843072 : total_deigen += (*deigen_darg)[_qp];
462 :
463 : // No eigenstrain -> no temperature coupling. Short-circuit before dereferencing the
464 : // d_sigma/d_eigenstrain property (only fetched when eigenstrains are coupled; see
465 : // LagrangianStressDivergenceBase ctor).
466 843072 : if (total_deigen.L2norm() == 0.0)
467 : return 0.0;
468 :
469 : // Direct chain through the constitutive update. The stress material publishes
470 : // d_sigma/d_eigenstrain (= -Jinv * small_jacobian for the objective-rate path,
471 : // with the sign convention that an eigenstrain increase
472 : // reduces the mechanical strain).
473 : // We then wrap to PK1 the same way the residual does:
474 : // dP/dT = det(F) * dsigma/dT * F^{-T} (large kinematics)
475 : // dP/dT = dsigma/dT (small kinematics, P == sigma by convention)
476 843072 : const RankTwoTensor dsigma_dT = (*_dcauchy_stress_d_eigenstrain)[_qp] * total_deigen;
477 843072 : RankTwoTensor dP_dT;
478 843072 : if (_large_kinematics)
479 791808 : dP_dT = _F[_qp].det() * dsigma_dT * _F_inv[_qp].transpose();
480 : else
481 51264 : dP_dT = dsigma_dT;
482 843072 : return dP_dT.doubleContraction(gradTest(_alpha)) * _temperature->phi()[_j][_qp];
483 : }
484 :
485 : template <class G>
486 : Real
487 85888 : TotalLagrangianStressDivergenceBase<G>::computeQpJacobianOutOfPlaneStrain()
488 : {
489 : // d(R_disp_alpha)/d(strain_zz_j) at qp = gradTest_alpha : d(PK1)/d(strain_zz_j).
490 : // strain_zz feeds `_F[(2,2)]` AFTER F-bar runs in `ComputeLagrangianWPSStrain`, so
491 : // strain_zz perturbations bypass F-bar's chain. Use `_dpk1_bypass_fbar` (= the
492 : // pk1_jacobian variant computed with the F-bar `_d_F_stab_d_F_ust` factor REPLACED
493 : // by identity in the sigma chain) for a consistent Jacobian.
494 85888 : return _dpk1_bypass_fbar[_qp].contractionKl(2, 2, gradTest(_alpha)) *
495 85888 : _out_of_plane_strain->phi()[_j][_qp];
496 : }
497 :
498 : template class TotalLagrangianStressDivergenceBase<GradientOperatorCartesian>;
499 : template class TotalLagrangianStressDivergenceBase<GradientOperatorAxisymmetricCylindrical>;
500 : template class TotalLagrangianStressDivergenceBase<GradientOperatorCentrosymmetricSpherical>;
|
