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 800 : HomogenizedTotalLagrangianStressDivergence::validParams()
20 : {
21 800 : InputParameters params = HomogenizationInterface<TotalLagrangianStressDivergence>::validParams();
22 800 : params.addClassDescription("Total Lagrangian stress equilibrium kernel with "
23 : "homogenization constraint Jacobian terms");
24 1600 : params.renameCoupledVar(
25 : "scalar_variable", "macro_var", "Optional scalar field with the macro gradient");
26 :
27 800 : return params;
28 0 : }
29 :
30 400 : HomogenizedTotalLagrangianStressDivergence::HomogenizedTotalLagrangianStressDivergence(
31 400 : const InputParameters & parameters)
32 400 : : HomogenizationInterface<TotalLagrangianStressDivergence>(parameters)
33 : {
34 400 : }
35 :
36 : std::set<std::string>
37 400 : 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 400 : vars.insert(_kappa_var_ptr->name());
42 400 : return vars;
43 : }
44 :
45 : void
46 4577434 : HomogenizedTotalLagrangianStressDivergence::computeScalarResidual()
47 : {
48 4577434 : if (_alpha != 0)
49 3035016 : return;
50 :
51 1542418 : std::vector<Real> scalar_residuals(_k_order);
52 :
53 : // only assemble scalar residual once; i.e. when handling the first displacement component
54 13690802 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
55 : {
56 12148384 : initScalarQpResidual();
57 12148384 : const auto dV = _JxW[_qp] * _coord[_qp];
58 :
59 : // index for residual vector
60 : unsigned int h = 0;
61 :
62 100382016 : for (const auto & [indices, constraint] : cmap())
63 : {
64 88233632 : const auto [i, j] = indices;
65 88233632 : const auto [ctype, ctarget] = constraint;
66 88233632 : const auto cval = ctarget->value(_t, _q_point[_qp]);
67 :
68 : // value to be constrained
69 : Real val;
70 88233632 : if (_large_kinematics)
71 : {
72 47878944 : if (ctype == Homogenization::ConstraintType::Stress)
73 15941168 : val = _pk1[_qp](i, j);
74 31937776 : else if (ctype == Homogenization::ConstraintType::Strain)
75 55830048 : 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 88233632 : scalar_residuals[h++] += (val - cval) * dV;
90 : }
91 : }
92 :
93 1542418 : addResiduals(
94 1542418 : _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
95 1542418 : }
96 :
97 : void
98 192290 : HomogenizedTotalLagrangianStressDivergence::computeScalarJacobian()
99 : {
100 192290 : if (_alpha != 0)
101 : return;
102 :
103 74825 : _local_ke.resize(_k_order, _k_order);
104 :
105 : // only assemble scalar residual once; i.e. when handling the first displacement component
106 549645 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
107 : {
108 474820 : initScalarQpJacobian(_kappa_var);
109 474820 : const auto dV = _JxW[_qp] * _coord[_qp];
110 :
111 : // index for Jacobian row
112 : unsigned int h = 0;
113 :
114 3870344 : for (const auto & [indices1, constraint1] : cmap())
115 : {
116 3395524 : const auto [i, j] = indices1;
117 3395524 : const auto ctype = constraint1.first;
118 :
119 : // index for Jacobian col
120 : unsigned int m = 0;
121 :
122 30684104 : for (const auto & [indices2, constraint2] : cmap())
123 : {
124 27288580 : const auto [k, l] = indices2;
125 27288580 : if (ctype == Homogenization::ConstraintType::Stress)
126 9624080 : _local_ke(h, m++) += dV * (_dpk1[_qp](i, j, k, l));
127 17664500 : else if (ctype == Homogenization::ConstraintType::Strain)
128 : {
129 17664500 : if (_large_kinematics)
130 30781276 : _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 3395524 : h++;
138 : }
139 : }
140 :
141 74825 : addJacobian(_assembly,
142 : _local_ke,
143 74825 : _kappa_var_ptr->dofIndices(),
144 74825 : _kappa_var_ptr->dofIndices(),
145 74825 : _kappa_var_ptr->scalingFactor());
146 : }
147 :
148 : void
149 517060 : 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 517060 : if (jvar_num != _var.number())
155 : return;
156 :
157 : const auto & jvar = getVariable(jvar_num);
158 191890 : const auto jvar_size = jvar.phiSize();
159 191890 : _local_ke.resize(_k_order, jvar_size);
160 :
161 1496730 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
162 : {
163 1304840 : const auto dV = _JxW[_qp] * _coord[_qp];
164 :
165 : // index for Jacobian row
166 : unsigned int h = 0;
167 :
168 11122192 : for (const auto & [indices, constraint] : cmap())
169 : {
170 9817352 : std::tie(_m, _n) = indices;
171 9817352 : _ctype = constraint.first;
172 9817352 : initScalarQpOffDiagJacobian(jvar);
173 85463208 : for (_j = 0; _j < jvar_size; _j++)
174 75645856 : _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num);
175 9817352 : h++;
176 : }
177 : }
178 :
179 191890 : addJacobian(_assembly,
180 : _local_ke,
181 191890 : _kappa_var_ptr->dofIndices(),
182 191890 : jvar.dofIndices(),
183 191890 : _kappa_var_ptr->scalingFactor());
184 : }
185 :
186 : void
187 191890 : HomogenizedTotalLagrangianStressDivergence::computeOffDiagJacobianScalarLocal(
188 : const unsigned int svar_num)
189 : {
190 : // Just in case, skip any other scalar variables
191 191890 : if (svar_num != _kappa_var)
192 : return;
193 :
194 191890 : _local_ke.resize(_test.size(), _k_order);
195 :
196 1496730 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
197 : {
198 : unsigned int l = 0;
199 1304840 : const auto dV = _JxW[_qp] * _coord[_qp];
200 11122192 : for (const auto & [indices, constraint] : cmap())
201 : {
202 : // copy constraint indices to protected variables to pass to Qp routine
203 9817352 : std::tie(_m, _n) = indices;
204 9817352 : _ctype = constraint.first;
205 9817352 : initScalarQpJacobian(svar_num);
206 85463208 : for (_i = 0; _i < _test.size(); _i++)
207 75645856 : _local_ke(_i, l) += dV * computeQpOffDiagJacobianScalar(svar_num);
208 9817352 : l++;
209 : }
210 : }
211 :
212 191890 : addJacobian(
213 191890 : _assembly, _local_ke, _var.dofIndices(), _kappa_var_ptr->dofIndices(), _var.scalingFactor());
214 : }
215 :
216 : Real
217 75645856 : HomogenizedTotalLagrangianStressDivergence::computeQpOffDiagJacobianScalar(
218 : unsigned int /*svar_num*/)
219 : {
220 75645856 : return _dpk1[_qp].contractionKl(_m, _n, gradTest(_alpha));
221 : }
222 :
223 : Real
224 75645856 : HomogenizedTotalLagrangianStressDivergence::computeScalarQpOffDiagJacobian(
225 : unsigned int /*jvar_num*/)
226 : {
227 75645856 : if (_ctype == Homogenization::ConstraintType::Stress)
228 27387712 : return _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha));
229 48258144 : else if (_ctype == Homogenization::ConstraintType::Strain)
230 48258144 : if (_large_kinematics)
231 71797600 : 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 : }
|
