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 "HomogenizedTotalLagrangianStressDivergence.h"
11 :
12 : // MOOSE includes
13 : #include "Function.h"
14 : #include "MooseVariableScalar.h"
15 :
16 : registerMooseObject("SolidMechanicsApp", HomogenizedTotalLagrangianStressDivergence);
17 :
18 : InputParameters
19 624 : HomogenizedTotalLagrangianStressDivergence::validParams()
20 : {
21 624 : InputParameters params = HomogenizationInterface<TotalLagrangianStressDivergence>::validParams();
22 624 : params.addClassDescription("Total Lagrangian stress equilibrium kernel with "
23 : "homogenization constraint Jacobian terms");
24 1248 : params.renameCoupledVar(
25 : "scalar_variable", "macro_var", "Optional scalar field with the macro gradient");
26 :
27 1248 : params.addParam<bool>(
28 1248 : "off_diagonal_jacobian", true, "Include the off diagonal parts of the constraint Jacobian");
29 :
30 624 : return params;
31 0 : }
32 :
33 312 : HomogenizedTotalLagrangianStressDivergence::HomogenizedTotalLagrangianStressDivergence(
34 312 : const InputParameters & parameters)
35 : : HomogenizationInterface<TotalLagrangianStressDivergence>(parameters),
36 624 : _use_off_diagonal(getParam<bool>("off_diagonal_jacobian"))
37 : {
38 312 : }
39 :
40 : std::set<std::string>
41 312 : HomogenizedTotalLagrangianStressDivergence::additionalROVariables()
42 : {
43 : // Add the scalar variable to the list of variables that this kernel contributes to
44 : std::set<std::string> vars = TotalLagrangianStressDivergence::additionalROVariables();
45 312 : vars.insert(_kappa_var_ptr->name());
46 312 : return vars;
47 : }
48 :
49 : void
50 2376640 : HomogenizedTotalLagrangianStressDivergence::computeScalarResidual()
51 : {
52 2376640 : if (_alpha != 0)
53 1574157 : return;
54 :
55 802483 : std::vector<Real> scalar_residuals(_k_order);
56 :
57 : // only assemble scalar residual once; i.e. when handling the first displacement component
58 7104231 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
59 : {
60 6301748 : initScalarQpResidual();
61 6301748 : const auto dV = _JxW[_qp] * _coord[_qp];
62 :
63 : // index for residual vector
64 : unsigned int h = 0;
65 :
66 52122856 : for (const auto & [indices, constraint] : cmap())
67 : {
68 45821108 : const auto [i, j] = indices;
69 45821108 : const auto [ctype, ctarget] = constraint;
70 45821108 : const auto cval = ctarget->value(_t, _q_point[_qp]);
71 :
72 : // value to be constrained
73 : Real val;
74 45821108 : if (_large_kinematics)
75 : {
76 25179316 : if (ctype == Homogenization::ConstraintType::Stress)
77 8266256 : val = _pk1[_qp](i, j);
78 16913060 : else if (ctype == Homogenization::ConstraintType::Strain)
79 29513748 : val = _F[_qp](i, j) - (Real(i == j));
80 : else
81 0 : mooseError("Unknown constraint type in the integral!");
82 : }
83 : else
84 : {
85 20641792 : if (ctype == Homogenization::ConstraintType::Stress)
86 10287360 : val = _pk1[_qp](i, j);
87 10354432 : else if (ctype == Homogenization::ConstraintType::Strain)
88 15506688 : val = 0.5 * (_F[_qp](i, j) + _F[_qp](j, i)) - (Real(i == j));
89 : else
90 0 : mooseError("Unknown constraint type in the integral!");
91 : }
92 :
93 45821108 : scalar_residuals[h++] += (val - cval) * dV;
94 : }
95 : }
96 :
97 802483 : addResiduals(
98 802483 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
99 802483 : }
100 :
101 : void
102 166892 : HomogenizedTotalLagrangianStressDivergence::computeScalarJacobian()
103 : {
104 166892 : if (_alpha != 0)
105 : return;
106 :
107 62270 : _local_ke.resize(_k_order, _k_order);
108 :
109 : // only assemble scalar residual once; i.e. when handling the first displacement component
110 483798 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
111 : {
112 421528 : initScalarQpJacobian(_kappa_var);
113 421528 : const auto dV = _JxW[_qp] * _coord[_qp];
114 :
115 : // index for Jacobian row
116 : unsigned int h = 0;
117 :
118 3532592 : for (const auto & [indices1, constraint1] : cmap())
119 : {
120 3111064 : const auto [i, j] = indices1;
121 3111064 : const auto ctype = constraint1.first;
122 :
123 : // index for Jacobian col
124 : unsigned int m = 0;
125 :
126 28262192 : for (const auto & [indices2, constraint2] : cmap())
127 : {
128 25151128 : const auto [k, l] = indices2;
129 25151128 : if (ctype == Homogenization::ConstraintType::Stress)
130 : // Macro_grad <-> macro_grad: scalar perturbation bypasses F-bar (the
131 : // homogenization material adds to `_F` AFTER F-bar runs), so use the
132 : // bypass variant of pk1_jacobian.
133 8114672 : _local_ke(h, m++) += dV * (_dpk1_bypass_fbar[_qp](i, j, k, l));
134 17036456 : else if (ctype == Homogenization::ConstraintType::Strain)
135 : {
136 17036456 : if (_large_kinematics)
137 28704688 : _local_ke(h, m++) += dV * (Real(i == k && j == l));
138 : else
139 4914168 : _local_ke(h, m++) += dV * (0.5 * Real(i == k && j == l) + 0.5 * Real(i == l && j == k));
140 : }
141 : else
142 0 : mooseError("Unknown constraint type in Jacobian calculator!");
143 : }
144 3111064 : h++;
145 : }
146 : }
147 :
148 62270 : addJacobian(_assembly,
149 : _local_ke,
150 62270 : _kappa_var_ptr->dofIndices(),
151 62270 : _kappa_var_ptr->dofIndices(),
152 62270 : _kappa_var_ptr->scalingFactor());
153 : }
154 :
155 : void
156 463784 : HomogenizedTotalLagrangianStressDivergence::computeScalarOffDiagJacobian(
157 : const unsigned int jvar_num)
158 : {
159 463784 : if (!_use_off_diagonal)
160 : return;
161 :
162 : // ONLY assemble the contribution from _alpha component, which is connected with _var
163 : // The other components are handled by other kernel instances with other _alpha
164 380840 : if (jvar_num != _var.number())
165 : return;
166 :
167 : const auto & jvar = getVariable(jvar_num);
168 138950 : const auto jvar_size = jvar.phiSize();
169 138950 : _local_ke.resize(_k_order, jvar_size);
170 :
171 : // The scalar<->disp Jacobian needs `_avg_grad_trial[_alpha]` populated (for the
172 : // non-local F-bar chain via `_d_F_stab_d_F_avg * deltaF_avg`). The base
173 : // `precalculateOffDiagJacobian` only does this when the off-diag jvar IS a
174 : // displacement; for the scalar-driven path it must be triggered explicitly.
175 138950 : if (_stabilize_strain)
176 : {
177 6378 : _fe_problem.prepareShapes(jvar_num, _tid);
178 6378 : _avg_grad_trial[_alpha].resize(_phi.size());
179 6378 : precalculateJacobianDisplacement(_alpha);
180 : }
181 :
182 1108062 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
183 : {
184 969112 : const auto dV = _JxW[_qp] * _coord[_qp];
185 :
186 : // index for Jacobian row
187 : unsigned int h = 0;
188 :
189 8304944 : for (const auto & [indices, constraint] : cmap())
190 : {
191 7335832 : std::tie(_m, _n) = indices;
192 7335832 : _ctype = constraint.first;
193 7335832 : initScalarQpOffDiagJacobian(jvar);
194 64239096 : for (_j = 0; _j < jvar_size; _j++)
195 56903264 : _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num);
196 7335832 : h++;
197 : }
198 : }
199 :
200 138950 : addJacobian(_assembly,
201 : _local_ke,
202 138950 : _kappa_var_ptr->dofIndices(),
203 138950 : jvar.dofIndices(),
204 138950 : _kappa_var_ptr->scalingFactor());
205 : }
206 :
207 : void
208 166598 : HomogenizedTotalLagrangianStressDivergence::computeOffDiagJacobianScalarLocal(
209 : const unsigned int svar_num)
210 : {
211 166598 : if (!_use_off_diagonal)
212 : return;
213 :
214 : // Just in case, skip any other scalar variables
215 138950 : if (svar_num != _kappa_var)
216 : return;
217 :
218 138950 : _local_ke.resize(_test.size(), _k_order);
219 :
220 1108062 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
221 : {
222 : unsigned int l = 0;
223 969112 : const auto dV = _JxW[_qp] * _coord[_qp];
224 8304944 : for (const auto & [indices, constraint] : cmap())
225 : {
226 : // copy constraint indices to protected variables to pass to Qp routine
227 7335832 : std::tie(_m, _n) = indices;
228 7335832 : _ctype = constraint.first;
229 7335832 : initScalarQpJacobian(svar_num);
230 64239096 : for (_i = 0; _i < _test.size(); _i++)
231 56903264 : _local_ke(_i, l) += dV * computeQpOffDiagJacobianScalar(svar_num);
232 7335832 : l++;
233 : }
234 : }
235 :
236 138950 : addJacobian(
237 138950 : _assembly, _local_ke, _var.dofIndices(), _kappa_var_ptr->dofIndices(), _var.scalingFactor());
238 : }
239 :
240 : Real
241 56903264 : HomogenizedTotalLagrangianStressDivergence::computeQpOffDiagJacobianScalar(
242 : unsigned int /*svar_num*/)
243 : {
244 : // d(disp residual) / d(scalar_{m,n}) = int gradTest_alpha : d(PK1)/d(scalar_{m,n}) dV.
245 : // The macro_gradient adds to `_F` AFTER F-bar runs (in
246 : // `ComputeLagrangianStrainBase::computeQpProperties`), so scalar perturbations
247 : // bypass F-bar's chain -- use `_dpk1_bypass_fbar` (pk1_jacobian with the F-bar
248 : // `_d_F_stab_d_F_ust` factor REPLACED by identity in the sigma chain).
249 56903264 : return _dpk1_bypass_fbar[_qp].contractionKl(_m, _n, gradTest(_alpha));
250 : }
251 :
252 : Real
253 56903264 : HomogenizedTotalLagrangianStressDivergence::computeScalarQpOffDiagJacobian(
254 : unsigned int /*jvar_num*/)
255 : {
256 56903264 : if (_ctype == Homogenization::ConstraintType::Stress)
257 : {
258 : // d(PK1_{m,n})/d(grad u_beta,j) -- local chain via _dpk1 (= dPK1/d(grad u) including
259 : // local F-bar effect via the sigma-chain through `_d_F_stab_d_F_ust`).
260 23587328 : Real J = _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha));
261 :
262 : // Non-local F-bar contribution to PK1 component (m, n) via the shared helper --
263 : // same chain as the regular TL displacement Jacobian but contracted into the single
264 : // (m, n) entry rather than doubled with gradTest. Guarded on `_stabilize_strain`
265 : // because `_avg_grad_trial` is only populated when F-bar is on.
266 23587328 : if (_stabilize_strain)
267 : {
268 331776 : const RankTwoTensor delta_F_avg = _d_F_d_grad_u[_qp] * _avg_grad_trial[_alpha][_j];
269 331776 : J += deltaPK1NonLocalFBar(delta_F_avg)(_m, _n);
270 : }
271 23587328 : return J;
272 : }
273 33315936 : else if (_ctype == Homogenization::ConstraintType::Strain)
274 : {
275 : // d(F_stab_{m,n})/d(disp_alpha_j) -- for F-bar on, the F-bar chain couples F_stab to
276 : // F_ust through both LOCAL (`_d_F_stab_d_F_ust`) and NON-LOCAL
277 : // (`_d_F_stab_d_F_avg * deltaF_avg`) routes. The old `Real(_m == _alpha) *
278 : // gradTrial(_m, _n)` form captured only the F-bar-off case correctly.
279 33315936 : const RankTwoTensor delta_F_ust_local = _d_F_d_grad_u[_qp] * gradTrialUnstabilized(_alpha);
280 33315936 : RankTwoTensor delta_F_stab = _d_F_stab_d_F_ust[_qp] * delta_F_ust_local;
281 33315936 : if (_stabilize_strain)
282 : {
283 516768 : const RankTwoTensor delta_F_avg = _d_F_d_grad_u[_qp] * _avg_grad_trial[_alpha][_j];
284 516768 : delta_F_stab += _d_F_stab_d_F_avg[_qp] * delta_F_avg;
285 : }
286 33315936 : if (_large_kinematics)
287 29755232 : return delta_F_stab(_m, _n);
288 : else
289 3560704 : return 0.5 * (delta_F_stab(_m, _n) + delta_F_stab(_n, _m));
290 : }
291 : else
292 0 : mooseError("Unknown constraint type in kernel calculation!");
293 : }
|
