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 576 : HomogenizedTotalLagrangianStressDivergence::validParams()
20 : {
21 576 : InputParameters params = HomogenizationInterface<TotalLagrangianStressDivergence>::validParams();
22 576 : params.addClassDescription("Total Lagrangian stress equilibrium kernel with "
23 : "homogenization constraint Jacobian terms");
24 1152 : params.renameCoupledVar(
25 : "scalar_variable", "macro_var", "Optional scalar field with the macro gradient");
26 :
27 576 : return params;
28 0 : }
29 :
30 288 : HomogenizedTotalLagrangianStressDivergence::HomogenizedTotalLagrangianStressDivergence(
31 288 : const InputParameters & parameters)
32 288 : : HomogenizationInterface<TotalLagrangianStressDivergence>(parameters)
33 : {
34 288 : }
35 :
36 : std::set<std::string>
37 288 : HomogenizedTotalLagrangianStressDivergence::additionalROVariables()
38 : {
39 : // Add the scalar variable to the list of variables that this kernel contributes to
40 : std::set<std::string> vars = TotalLagrangianStressDivergence::additionalROVariables();
41 288 : vars.insert(_kappa_var_ptr->name());
42 288 : return vars;
43 : }
44 :
45 : void
46 2339776 : HomogenizedTotalLagrangianStressDivergence::computeScalarResidual()
47 : {
48 2339776 : if (_alpha != 0)
49 1549581 : return;
50 :
51 790195 : std::vector<Real> scalar_residuals(_k_order);
52 :
53 : // only assemble scalar residual once; i.e. when handling the first displacement component
54 6993639 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
55 : {
56 6203444 : initScalarQpResidual();
57 6203444 : const auto dV = _JxW[_qp] * _coord[_qp];
58 :
59 : // index for residual vector
60 : unsigned int h = 0;
61 :
62 51250408 : for (const auto & [indices, constraint] : cmap())
63 : {
64 45046964 : const auto [i, j] = indices;
65 45046964 : const auto [ctype, ctarget] = constraint;
66 45046964 : const auto cval = ctarget->value(_t, _q_point[_qp]);
67 :
68 : // value to be constrained
69 : Real val;
70 45046964 : if (_large_kinematics)
71 : {
72 24626356 : if (ctype == Homogenization::ConstraintType::Stress)
73 8266256 : val = _pk1[_qp](i, j);
74 16360100 : else if (ctype == Homogenization::ConstraintType::Strain)
75 28592148 : val = _F[_qp](i, j) - (Real(i == j));
76 : else
77 0 : mooseError("Unknown constraint type in the integral!");
78 : }
79 : else
80 : {
81 20420608 : if (ctype == Homogenization::ConstraintType::Stress)
82 10287360 : val = _pk1[_qp](i, j);
83 10133248 : else if (ctype == Homogenization::ConstraintType::Strain)
84 15174912 : val = 0.5 * (_F[_qp](i, j) + _F[_qp](j, i)) - (Real(i == j));
85 : else
86 0 : mooseError("Unknown constraint type in the integral!");
87 : }
88 :
89 45046964 : scalar_residuals[h++] += (val - cval) * dV;
90 : }
91 : }
92 :
93 790195 : addResiduals(
94 790195 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
95 790195 : }
96 :
97 : void
98 139244 : HomogenizedTotalLagrangianStressDivergence::computeScalarJacobian()
99 : {
100 139244 : if (_alpha != 0)
101 : return;
102 :
103 53054 : _local_ke.resize(_k_order, _k_order);
104 :
105 : // only assemble scalar residual once; i.e. when handling the first displacement component
106 400854 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
107 : {
108 347800 : initScalarQpJacobian(_kappa_var);
109 347800 : const auto dV = _JxW[_qp] * _coord[_qp];
110 :
111 : // index for Jacobian row
112 : unsigned int h = 0;
113 :
114 2869040 : for (const auto & [indices1, constraint1] : cmap())
115 : {
116 2521240 : const auto [i, j] = indices1;
117 2521240 : const auto ctype = constraint1.first;
118 :
119 : // index for Jacobian col
120 : unsigned int m = 0;
121 :
122 22806320 : for (const auto & [indices2, constraint2] : cmap())
123 : {
124 20285080 : const auto [k, l] = indices2;
125 20285080 : if (ctype == Homogenization::ConstraintType::Stress)
126 8114672 : _local_ke(h, m++) += dV * (_dpk1[_qp](i, j, k, l));
127 12170408 : else if (ctype == Homogenization::ConstraintType::Strain)
128 : {
129 12170408 : if (_large_kinematics)
130 21184432 : _local_ke(h, m++) += dV * (Real(i == k && j == l));
131 : else
132 2481144 : _local_ke(h, m++) += dV * (0.5 * Real(i == k && j == l) + 0.5 * Real(i == l && j == k));
133 : }
134 : else
135 0 : mooseError("Unknown constraint type in Jacobian calculator!");
136 : }
137 2521240 : h++;
138 : }
139 : }
140 :
141 53054 : addJacobian(_assembly,
142 : _local_ke,
143 53054 : _kappa_var_ptr->dofIndices(),
144 53054 : _kappa_var_ptr->dofIndices(),
145 53054 : _kappa_var_ptr->scalingFactor());
146 : }
147 :
148 : void
149 380840 : HomogenizedTotalLagrangianStressDivergence::computeScalarOffDiagJacobian(
150 : const unsigned int jvar_num)
151 : {
152 : // ONLY assemble the contribution from _alpha component, which is connected with _var
153 : // The other components are handled by other kernel instances with other _alpha
154 380840 : if (jvar_num != _var.number())
155 : return;
156 :
157 : const auto & jvar = getVariable(jvar_num);
158 138950 : const auto jvar_size = jvar.phiSize();
159 138950 : _local_ke.resize(_k_order, jvar_size);
160 :
161 1108062 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
162 : {
163 969112 : const auto dV = _JxW[_qp] * _coord[_qp];
164 :
165 : // index for Jacobian row
166 : unsigned int h = 0;
167 :
168 8304944 : for (const auto & [indices, constraint] : cmap())
169 : {
170 7335832 : std::tie(_m, _n) = indices;
171 7335832 : _ctype = constraint.first;
172 7335832 : initScalarQpOffDiagJacobian(jvar);
173 64239096 : for (_j = 0; _j < jvar_size; _j++)
174 56903264 : _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num);
175 7335832 : h++;
176 : }
177 : }
178 :
179 138950 : addJacobian(_assembly,
180 : _local_ke,
181 138950 : _kappa_var_ptr->dofIndices(),
182 138950 : jvar.dofIndices(),
183 138950 : _kappa_var_ptr->scalingFactor());
184 : }
185 :
186 : void
187 138950 : HomogenizedTotalLagrangianStressDivergence::computeOffDiagJacobianScalarLocal(
188 : const unsigned int svar_num)
189 : {
190 : // Just in case, skip any other scalar variables
191 138950 : if (svar_num != _kappa_var)
192 : return;
193 :
194 138950 : _local_ke.resize(_test.size(), _k_order);
195 :
196 1108062 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
197 : {
198 : unsigned int l = 0;
199 969112 : const auto dV = _JxW[_qp] * _coord[_qp];
200 8304944 : for (const auto & [indices, constraint] : cmap())
201 : {
202 : // copy constraint indices to protected variables to pass to Qp routine
203 7335832 : std::tie(_m, _n) = indices;
204 7335832 : _ctype = constraint.first;
205 7335832 : initScalarQpJacobian(svar_num);
206 64239096 : for (_i = 0; _i < _test.size(); _i++)
207 56903264 : _local_ke(_i, l) += dV * computeQpOffDiagJacobianScalar(svar_num);
208 7335832 : l++;
209 : }
210 : }
211 :
212 138950 : addJacobian(
213 138950 : _assembly, _local_ke, _var.dofIndices(), _kappa_var_ptr->dofIndices(), _var.scalingFactor());
214 : }
215 :
216 : Real
217 56903264 : HomogenizedTotalLagrangianStressDivergence::computeQpOffDiagJacobianScalar(
218 : unsigned int /*svar_num*/)
219 : {
220 56903264 : return _dpk1[_qp].contractionKl(_m, _n, gradTest(_alpha));
221 : }
222 :
223 : Real
224 56903264 : HomogenizedTotalLagrangianStressDivergence::computeScalarQpOffDiagJacobian(
225 : unsigned int /*jvar_num*/)
226 : {
227 56903264 : if (_ctype == Homogenization::ConstraintType::Stress)
228 23587328 : return _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha));
229 33315936 : else if (_ctype == Homogenization::ConstraintType::Strain)
230 33315936 : if (_large_kinematics)
231 49461728 : return Real(_m == _alpha) * gradTrial(_alpha)(_m, _n);
232 : else
233 9460224 : return 0.5 * (Real(_m == _alpha) * gradTrial(_alpha)(_m, _n) +
234 5899520 : Real(_n == _alpha) * gradTrial(_alpha)(_n, _m));
235 : else
236 0 : mooseError("Unknown constraint type in kernel calculation!");
237 : }
|
