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 828 : HomogenizedTotalLagrangianStressDivergence::validParams()
20 : {
21 828 : InputParameters params = HomogenizationInterface<TotalLagrangianStressDivergence>::validParams();
22 828 : params.addClassDescription("Total Lagrangian stress equilibrium kernel with "
23 : "homogenization constraint Jacobian terms");
24 1656 : params.renameCoupledVar(
25 : "scalar_variable", "macro_var", "Optional scalar field with the macro gradient");
26 :
27 828 : return params;
28 0 : }
29 :
30 414 : HomogenizedTotalLagrangianStressDivergence::HomogenizedTotalLagrangianStressDivergence(
31 414 : const InputParameters & parameters)
32 414 : : HomogenizationInterface<TotalLagrangianStressDivergence>(parameters)
33 : {
34 414 : }
35 :
36 : std::set<std::string>
37 414 : 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 414 : vars.insert(_kappa_var_ptr->name());
42 414 : return vars;
43 : }
44 :
45 : void
46 4578028 : HomogenizedTotalLagrangianStressDivergence::computeScalarResidual()
47 : {
48 4578028 : if (_alpha != 0)
49 3035313 : return;
50 :
51 1542715 : std::vector<Real> scalar_residuals(_k_order);
52 :
53 : // only assemble scalar residual once; i.e. when handling the first displacement component
54 13692287 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
55 : {
56 12149572 : initScalarQpResidual();
57 12149572 : const auto dV = _JxW[_qp] * _coord[_qp];
58 :
59 : // index for residual vector
60 : unsigned int h = 0;
61 :
62 100384392 : for (const auto & [indices, constraint] : cmap())
63 : {
64 88234820 : const auto [i, j] = indices;
65 88234820 : const auto [ctype, ctarget] = constraint;
66 88234820 : const auto cval = ctarget->value(_t, _q_point[_qp]);
67 :
68 : // value to be constrained
69 : Real val;
70 88234820 : if (_large_kinematics)
71 : {
72 47880132 : if (ctype == Homogenization::ConstraintType::Stress)
73 15941168 : val = _pk1[_qp](i, j);
74 31938964 : else if (ctype == Homogenization::ConstraintType::Strain)
75 55831236 : val = _F[_qp](i, j) - (Real(i == j));
76 : else
77 0 : mooseError("Unknown constraint type in the integral!");
78 : }
79 : else
80 : {
81 40354688 : if (ctype == Homogenization::ConstraintType::Stress)
82 20227680 : val = _pk1[_qp](i, j);
83 20127008 : else if (ctype == Homogenization::ConstraintType::Strain)
84 30168672 : 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 88234820 : scalar_residuals[h++] += (val - cval) * dV;
90 : }
91 : }
92 :
93 1542715 : addResiduals(
94 1542715 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
95 1542715 : }
96 :
97 : void
98 192776 : HomogenizedTotalLagrangianStressDivergence::computeScalarJacobian()
99 : {
100 192776 : if (_alpha != 0)
101 : return;
102 :
103 75068 : _local_ke.resize(_k_order, _k_order);
104 :
105 : // only assemble scalar residual once; i.e. when handling the first displacement component
106 550860 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
107 : {
108 475792 : initScalarQpJacobian(_kappa_var);
109 475792 : const auto dV = _JxW[_qp] * _coord[_qp];
110 :
111 : // index for Jacobian row
112 : unsigned int h = 0;
113 :
114 3872288 : for (const auto & [indices1, constraint1] : cmap())
115 : {
116 3396496 : const auto [i, j] = indices1;
117 3396496 : const auto ctype = constraint1.first;
118 :
119 : // index for Jacobian col
120 : unsigned int m = 0;
121 :
122 30686048 : for (const auto & [indices2, constraint2] : cmap())
123 : {
124 27289552 : const auto [k, l] = indices2;
125 27289552 : if (ctype == Homogenization::ConstraintType::Stress)
126 9624080 : _local_ke(h, m++) += dV * (_dpk1[_qp](i, j, k, l));
127 17665472 : else if (ctype == Homogenization::ConstraintType::Strain)
128 : {
129 17665472 : if (_large_kinematics)
130 30782248 : _local_ke(h, m++) += dV * (Real(i == k && j == l));
131 : else
132 3546872 : _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 3396496 : h++;
138 : }
139 : }
140 :
141 75068 : addJacobian(_assembly,
142 : _local_ke,
143 75068 : _kappa_var_ptr->dofIndices(),
144 75068 : _kappa_var_ptr->dofIndices(),
145 75068 : _kappa_var_ptr->scalingFactor());
146 : }
147 :
148 : void
149 518248 : 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 518248 : if (jvar_num != _var.number())
155 : return;
156 :
157 : const auto & jvar = getVariable(jvar_num);
158 192286 : const auto jvar_size = jvar.phiSize();
159 192286 : _local_ke.resize(_k_order, jvar_size);
160 :
161 1498710 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
162 : {
163 1306424 : const auto dV = _JxW[_qp] * _coord[_qp];
164 :
165 : // index for Jacobian row
166 : unsigned int h = 0;
167 :
168 11125360 : for (const auto & [indices, constraint] : cmap())
169 : {
170 9818936 : std::tie(_m, _n) = indices;
171 9818936 : _ctype = constraint.first;
172 9818936 : initScalarQpOffDiagJacobian(jvar);
173 85471128 : for (_j = 0; _j < jvar_size; _j++)
174 75652192 : _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num);
175 9818936 : h++;
176 : }
177 : }
178 :
179 192286 : addJacobian(_assembly,
180 : _local_ke,
181 192286 : _kappa_var_ptr->dofIndices(),
182 192286 : jvar.dofIndices(),
183 192286 : _kappa_var_ptr->scalingFactor());
184 : }
185 :
186 : void
187 192286 : HomogenizedTotalLagrangianStressDivergence::computeOffDiagJacobianScalarLocal(
188 : const unsigned int svar_num)
189 : {
190 : // Just in case, skip any other scalar variables
191 192286 : if (svar_num != _kappa_var)
192 : return;
193 :
194 192286 : _local_ke.resize(_test.size(), _k_order);
195 :
196 1498710 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
197 : {
198 : unsigned int l = 0;
199 1306424 : const auto dV = _JxW[_qp] * _coord[_qp];
200 11125360 : for (const auto & [indices, constraint] : cmap())
201 : {
202 : // copy constraint indices to protected variables to pass to Qp routine
203 9818936 : std::tie(_m, _n) = indices;
204 9818936 : _ctype = constraint.first;
205 9818936 : initScalarQpJacobian(svar_num);
206 85471128 : for (_i = 0; _i < _test.size(); _i++)
207 75652192 : _local_ke(_i, l) += dV * computeQpOffDiagJacobianScalar(svar_num);
208 9818936 : l++;
209 : }
210 : }
211 :
212 192286 : addJacobian(
213 192286 : _assembly, _local_ke, _var.dofIndices(), _kappa_var_ptr->dofIndices(), _var.scalingFactor());
214 : }
215 :
216 : Real
217 75652192 : HomogenizedTotalLagrangianStressDivergence::computeQpOffDiagJacobianScalar(
218 : unsigned int /*svar_num*/)
219 : {
220 75652192 : return _dpk1[_qp].contractionKl(_m, _n, gradTest(_alpha));
221 : }
222 :
223 : Real
224 75652192 : HomogenizedTotalLagrangianStressDivergence::computeScalarQpOffDiagJacobian(
225 : unsigned int /*jvar_num*/)
226 : {
227 75652192 : if (_ctype == Homogenization::ConstraintType::Stress)
228 27387712 : return _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha));
229 48264480 : else if (_ctype == Homogenization::ConstraintType::Strain)
230 48264480 : if (_large_kinematics)
231 71807104 : return Real(_m == _alpha) * gradTrial(_alpha)(_m, _n);
232 : else
233 13419072 : return 0.5 * (Real(_m == _alpha) * gradTrial(_alpha)(_m, _n) +
234 8365600 : Real(_n == _alpha) * gradTrial(_alpha)(_n, _m));
235 : else
236 0 : mooseError("Unknown constraint type in kernel calculation!");
237 : }
|
