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