11 #include "DenseMatrix.h"
13 #include "libmesh/quadrature.h"
14 #include "libmesh/utility.h"
15 #include "libmesh/enum_quadrature_type.h"
16 #include "libmesh/fe_type.h"
17 #include "libmesh/string_to_enum.h"
18 #include "libmesh/quadrature_gauss.h"
24 params.addClassDescription(
"Compute a large strain increment for the shell."););
26 template <ComputeStage compute_stage>
28 const InputParameters & parameters)
33 for (
unsigned int i = 0; i <
_t_points.size(); ++i)
34 _B_nl[i] = &declareADProperty<DenseMatrix<Real>>(
"B_nl_t_points_" + std::to_string(i));
37 template <ComputeStage compute_stage>
43 ADDenseMatrix b(5, 20);
44 for (
unsigned int t = 0; t < _t_points.size(); ++t)
48 template <ComputeStage compute_stage>
53 _2d_points = _qrule->get_points();
55 unsigned int dim = _current_elem->dim();
59 FEType fe_type(Utility::string_to_enum<Order>(
"First"),
60 Utility::string_to_enum<FEFamily>(
"LAGRANGE"));
61 auto & fe = _fe_problem.assembly(_tid).getFE(fe_type, dim);
62 _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
63 _dphideta_map = fe->get_fe_map().get_dphideta_map();
64 _phi_map = fe->get_fe_map().get_phi_map();
66 for (
unsigned int i = 0; i < 4; ++i)
67 _nodes[i] = _current_elem->node_ptr(i);
69 for (
unsigned int i = 0; i < _2d_points.size(); ++i)
71 for (
unsigned int j = 0; j < _t_points.size(); ++j)
73 (*_ge[j])[i] = (*_ge_old[j])[i];
74 (*_J_map[j])[i] = (*_J_map_old[j])[i];
81 for (
unsigned int k = 0; k < _nodes.size(); ++k)
82 _node_normal[k] = _node_normal_old[k];
90 for (
unsigned int i = 0; i < _2d_points.size(); ++i)
92 for (
unsigned int j = 0; j < _t_points.size(); ++j)
95 for (
unsigned int temp1 = 0; temp1 < 5; ++temp1)
97 _strain_vector(temp1) = 0.0;
98 for (
unsigned int temp2 = 0; temp2 < 20; ++temp2)
99 _strain_vector(temp1) += (*_B[j])[i](temp1, temp2) * _soln_vector(temp2);
102 (*_strain_increment[j])[i](0, 0) = _strain_vector(0);
103 (*_strain_increment[j])[i](1, 1) = _strain_vector(1);
104 (*_strain_increment[j])[i](0, 1) = _strain_vector(2);
105 (*_strain_increment[j])[i](0, 2) = _strain_vector(3);
106 (*_strain_increment[j])[i](1, 2) = _strain_vector(4);
107 (*_strain_increment[j])[i](1, 0) = (*_strain_increment[j])[i](0, 1);
108 (*_strain_increment[j])[i](2, 0) = (*_strain_increment[j])[i](0, 2);
109 (*_strain_increment[j])[i](2, 1) = (*_strain_increment[j])[i](1, 2);
110 (*_total_strain[j])[i] = (*_total_strain_old[j])[i] + (*_strain_increment[j])[i];
113 copyDualNumbersToValues();
116 template <ComputeStage compute_stage>
121 for (
unsigned int k = 0; k < _nodes.size(); ++k)
124 -_v2[k] * _soln_vector(12 + k) + _v1[k] * _soln_vector(16 + k) + _node_normal_old[k];
125 _node_normal[k] /= _node_normal[k].norm();
129 template <ComputeStage compute_stage>
133 for (
unsigned int i = 0; i < _2d_points.size(); ++i)
135 for (
unsigned int j = 0; j < _t_points.size(); ++j)
142 for (
unsigned int k = 0; k < _nodes.size(); ++k)
147 _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] *
150 _dphideta_map[k][i] *
152 _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] *
155 _thickness[i] * _phi_map[k][i] * _node_normal_old[k](
component) / 2.0;
162 template <ComputeStage compute_stage>
169 0.5 * (_sol_old(_soln_disp_index[2][
component]) - _sol_old(_soln_disp_index[3][
component]));
171 0.5 * (_sol_old(_soln_disp_index[1][
component]) - _sol_old(_soln_disp_index[0][
component]));
173 0.5 * (_sol_old(_soln_disp_index[3][
component]) - _sol_old(_soln_disp_index[0][
component]));
175 0.5 * (_sol_old(_soln_disp_index[2][
component]) - _sol_old(_soln_disp_index[1][
component]));
179 template <ComputeStage compute_stage>
186 for (
unsigned int i = 0; i < _2d_points.size(); ++i)
188 for (
unsigned int j = 0; j < _t_points.size(); ++j)
190 (*_B_nl[j])[i].resize(5, 20);
191 (*_B_nl[j])[i].zero();
192 for (
unsigned int k = 0; k < 4; ++k)
194 for (
unsigned int p = 0; p < 4; ++p)
197 (*_B_nl[j])[i](0, k) += _dphidxi_map[k][i] * _dphidxi_map[p][i] *
198 (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
199 (-_soln_vector(p + 12) * _v2[p](0) +
200 _soln_vector(p + 16) * _v1[p](0)));
201 (*_B_nl[j])[i](0, 4 + k) +=
202 _dphidxi_map[k][i] * _dphidxi_map[p][i] *
203 (_soln_vector(p + 4) +
204 _t_points[j](0) / 2.0 * _thickness[i] *
205 (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
206 (*_B_nl[j])[i](0, 8 + k) +=
207 _dphidxi_map[k][i] * _dphidxi_map[p][i] *
208 (_soln_vector(p + 8) +
209 _t_points[j](0) / 2.0 * _thickness[i] *
210 (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
211 (*_B_nl[j])[i](0, 12 + k) +=
212 _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] * _dphidxi_map[p][i] *
213 (-(_v2[p](0) * _soln_vector(p) + _v2[p](1) * _soln_vector(p + 4) +
214 _v2[p](2) * _soln_vector(p + 8)) +
215 _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
216 (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
217 (*_B_nl[j])[i](0, 16 + k) +=
218 _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] * _dphidxi_map[p][i] *
219 ((_v1[p](0) * _soln_vector(p) + _v1[p](1) * _soln_vector(p + 4) +
220 _v1[p](2) * _soln_vector(p + 8)) +
221 _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
222 (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
225 (*_B_nl[j])[i](1, k) += _dphideta_map[k][i] * _dphideta_map[p][i] *
226 (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
227 (-_soln_vector(p + 12) * _v2[p](0) +
228 _soln_vector(p + 16) * _v1[p](0)));
229 (*_B_nl[j])[i](1, 4 + k) +=
230 _dphideta_map[k][i] * _dphideta_map[p][i] *
231 (_soln_vector(p + 4) +
232 _t_points[j](0) / 2.0 * _thickness[i] *
233 (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
234 (*_B_nl[j])[i](1, 8 + k) +=
235 _dphideta_map[k][i] * _dphideta_map[p][i] *
236 (_soln_vector(p + 8) +
237 _t_points[j](0) / 2.0 * _thickness[i] *
238 (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
239 (*_B_nl[j])[i](1, 12 + k) +=
240 _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] * _dphideta_map[p][i] *
241 (-(_v2[p](0) * _soln_vector(p) + _v2[p](1) * _soln_vector(p + 4) +
242 _v2[p](2) * _soln_vector(p + 8)) +
243 _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
244 (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
245 (*_B_nl[j])[i](1, 16 + k) +=
246 _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] * _dphideta_map[p][i] *
247 ((_v1[p](0) * _soln_vector(p) + _v1[p](1) * _soln_vector(p + 4) +
248 _v1[p](2) * _soln_vector(p + 8)) +
249 _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
250 (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
255 (*_B_nl[j])[i](2, k) += 0.5 *
256 (_dphidxi_map[k][i] * _dphideta_map[p][i] +
257 _dphideta_map[k][i] * _dphidxi_map[p][i]) *
258 (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
259 (-_soln_vector(p + 12) * _v2[p](0) +
260 _soln_vector(p + 16) * _v1[p](0)));
261 (*_B_nl[j])[i](2, 4 + k) +=
263 (_dphidxi_map[k][i] * _dphideta_map[p][i] +
264 _dphideta_map[k][i] * _dphidxi_map[p][i]) *
265 (_soln_vector(p + 4) +
266 _t_points[j](0) / 2.0 * _thickness[i] *
267 (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
268 (*_B_nl[j])[i](2, 8 + k) +=
270 (_dphidxi_map[k][i] * _dphideta_map[p][i] +
271 _dphideta_map[k][i] * _dphidxi_map[p][i]) *
272 (_soln_vector(p + 8) +
273 _t_points[j](0) / 2.0 * _thickness[i] *
274 (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
275 (*_B_nl[j])[i](2, 12 + k) +=
276 _t_points[j](0) * 0.25 *
277 (_dphidxi_map[k][i] * _dphideta_map[p][i] +
278 _dphideta_map[k][i] * _dphidxi_map[p][i]) *
280 (-(_v2[k](0) * _soln_vector(p) + _v2[k](1) * _soln_vector(p + 4) +
281 _v2[k](2) * _soln_vector(p + 8)) +
282 _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
283 (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
284 (*_B_nl[j])[i](2, 16 + k) +=
285 _t_points[j](0) * 0.25 *
286 (_dphidxi_map[k][i] * _dphideta_map[p][i] +
287 _dphideta_map[k][i] * _dphidxi_map[p][i]) *
289 ((_v1[k](0) * _soln_vector(p) + _v1[k](1) * _soln_vector(p + 4) +
290 _v1[k](2) * _soln_vector(p + 8)) +
291 _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
292 (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
300 1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] *
301 (-_soln_vector(12 + 2) * _v2[2](
component) + _soln_vector(16 + 2) * _v1[2](
component) -
306 1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] *
307 (-_soln_vector(12 + 1) * _v2[1](
component) + _soln_vector(16 + 1) * _v1[1](
component) -
312 (*_B_nl[j])[i](3, 12 + 2) +=
313 -1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v2[2](
component) *
315 (*_B_nl[j])[i](3, 16 + 2) +=
316 1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v1[2](
component) *
318 (*_B_nl[j])[i](3, 12 + 3) +=
319 -1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v2[3](
component) *
321 (*_B_nl[j])[i](3, 16 + 3) +=
322 1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v1[3](
component) *
326 (*_B_nl[j])[i](3, 12 + 1) +=
327 -1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v2[1](
component) *
329 (*_B_nl[j])[i](3, 16 + 1) +=
330 1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v1[1](
component) *
332 (*_B_nl[j])[i](3, 12 + 0) +=
333 -1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v2[0](
component) *
335 (*_B_nl[j])[i](3, 16 + 0) +=
336 1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v1[0](
component) *
341 1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] *
342 (-_soln_vector(12 + 2) * _v2[2](
component) + _soln_vector(16 + 2) * _v1[2](
component) -
347 1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] *
348 (-_soln_vector(12 + 3) * _v2[3](
component) + _soln_vector(16 + 3) * _v1[3](
component) -
353 (*_B_nl[j])[i](4, 12 + 2) +=
354 -1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v2[2](
component) *
356 (*_B_nl[j])[i](4, 16 + 2) +=
357 1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v1[2](
component) *
359 (*_B_nl[j])[i](4, 12 + 1) +=
360 -1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v2[1](
component) *
362 (*_B_nl[j])[i](4, 16 + 1) +=
363 1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v1[1](
component) *
367 (*_B_nl[j])[i](4, 12 + 3) +=
368 -1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v2[3](
component) *
370 (*_B_nl[j])[i](4, 16 + 3) +=
371 1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v1[3](
component) *
373 (*_B_nl[j])[i](4, 12 + 0) +=
374 -1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v2[0](
component) *
376 (*_B_nl[j])[i](4, 16 + 0) +=
377 1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v1[0](
component) *