https://mooseframework.inl.gov
Loading...
Searching...
No Matches
ADInertialForceShell.C
Go to the documentation of this file.
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 "ADInertialForceShell.h"
11
12#include "SubProblem.h"
13#include "libmesh/utility.h"
14#include "MooseVariable.h"
15#include "Assembly.h"
16#include "NonlinearSystem.h"
17#include "AuxiliarySystem.h"
18#include "MooseMesh.h"
19#include "DenseMatrix.h"
20
21#include "libmesh/utility.h"
22#include "libmesh/quadrature.h"
23#include "libmesh/enum_quadrature_type.h"
24#include "libmesh/fe_type.h"
25#include "libmesh/string_to_enum.h"
26#include "libmesh/quadrature_gauss.h"
27#include "MooseVariable.h"
28
29registerMooseObject("SolidMechanicsApp", ADInertialForceShell);
30
31InputParameters
32ADInertialForceShell::validParams()
33{
34 InputParameters params = ADTimeKernel::validParams();
35
36 params.addClassDescription("Calculates the residual for the inertial force/moment and the "
37 "contribution of mass dependent Rayleigh damping and HHT time "
38 "integration scheme.");
39 params.set<bool>("use_displaced_mesh") = true;
40 params.addRequiredCoupledVar(
41 "rotations",
42 "The rotational variables appropriate for the simulation geometry and coordinate system");
43 params.addRequiredCoupledVar(
44 "displacements",
45 "The displacement variables appropriate for the simulation geometry and coordinate system");
46 params.addCoupledVar("velocities", "Translational velocity variables");
47 params.addCoupledVar("accelerations", "Translational acceleration variables");
48 params.addCoupledVar("rotational_velocities", "Rotational velocity variables");
49 params.addCoupledVar("rotational_accelerations", "Rotational acceleration variables");
50 params.addParam<MaterialPropertyName>(
51 "density",
52 "density",
53 "Name of Material Property or a constant real number defining the density of the beam.");
54 params.addRequiredRangeCheckedParam<unsigned int>(
55 "component",
56 "component<5",
57 "An integer corresponding to the direction "
58 "the variable this kernel acts in. (0 for disp_x, "
59 "1 for disp_y, 2 for disp_z, 3 for alpha, and 4 for beta)");
60 params.addRequiredParam<Real>("thickness", "The kernel's thickness");
61 params.addParam<MaterialPropertyName>("eta",
62 0.0,
63 "Name of material property or a constant real "
64 "number defining the eta parameter for the "
65 "Rayleigh damping.");
66 params.addRangeCheckedParam<Real>("alpha",
67 0.0,
68 "alpha >= -0.3333 & alpha <= 0.0",
69 "Alpha parameter for mass dependent numerical damping induced "
70 "by HHT time integration scheme");
71
72 return params;
73}
74
75ADInertialForceShell::ADInertialForceShell(const InputParameters & parameters)
76 : ADTimeKernel(parameters),
77 _nrot(coupledComponents("rotations")),
78 _ndisp(coupledComponents("displacements")),
79 _rot_num(_nrot),
80 _disp_num(_ndisp),
81 _vel_num(_ndisp),
82 _accel_num(_ndisp),
83 _rot_vel_num(_nrot),
84 _rot_accel_num(_nrot),
85 _component(getParam<unsigned int>("component")),
86 _nodes(4),
87 _v1(4),
88 _v2(4),
89 _node_normal(4),
90 _eta(getMaterialProperty<Real>("eta")),
91 _transformation_matrix(getADMaterialProperty<RankTwoTensor>("transformation_matrix_element")),
92 _J_map(getADMaterialProperty<Real>("J_mapping_t_points_0")),
93 _thickness(getParam<Real>("thickness")),
94 _density(getMaterialProperty<Real>("density")),
95 _alpha(getParam<Real>("alpha"))
96{
97 // Checking for consistency between the length of the provided rotations and displacements vector
98 if (_ndisp != 3 || _nrot != 2)
99 mooseError("InertialForceShell: The number of variables supplied in 'displacements' "
100 "must be 3 and that in 'rotations' must be 2.");
101
102 // fetch coupled displacements and rotations
103 for (unsigned int i = 0; i < _ndisp; ++i)
104 {
105 MooseVariable * disp_variable = getVar("displacements", i);
106 _disp_num[i] = disp_variable->number();
107 }
108 for (unsigned int i = 0; i < _nrot; ++i)
109 {
110 MooseVariable * rot_variable = getVar("rotations", i);
111 _rot_num[i] = rot_variable->number();
112 }
113
114 _x2(1) = 1;
115 _x3(2) = 1;
116}
117
118void
119ADInertialForceShell::computeResidual()
120{
121 prepareVectorTag(_assembly, _var.number());
122
123 precalculateResidual();
124
125 /* ----------------------------------------------------- */
126
127 if (_dt != 0.0)
128 {
129 computeShellInertialForces(_ad_coord, _ad_JxW);
130
131 if (_component < 3)
132 {
133 _global_force[0] = _thickness * _original_local_config.transpose() * _local_force[0];
134 _global_force[1] = _thickness * _original_local_config.transpose() * _local_force[1];
135 _global_force[2] = _thickness * _original_local_config.transpose() * _local_force[2];
136 _global_force[3] = _thickness * _original_local_config.transpose() * _local_force[3];
137
138 _local_re(0) = raw_value(_global_force[0](_component));
139 _local_re(1) = raw_value(_global_force[1](_component));
140 _local_re(2) = raw_value(_global_force[2](_component));
141 _local_re(3) = raw_value(_global_force[3](_component));
142 }
143 else
144 {
145 _global_moment[0] = _original_local_config.transpose() * _local_moment[0];
146 _global_moment[1] = _original_local_config.transpose() * _local_moment[1];
147 _global_moment[2] = _original_local_config.transpose() * _local_moment[2];
148 _global_moment[3] = _original_local_config.transpose() * _local_moment[3];
149
150 _local_re(0) = raw_value(_global_moment[0](_component - 3));
151 _local_re(1) = raw_value(_global_moment[1](_component - 3));
152 _local_re(2) = raw_value(_global_moment[2](_component - 3));
153 _local_re(3) = raw_value(_global_moment[3](_component - 3));
154 }
155 }
156
157 accumulateTaggedLocalResidual();
158
159 if (_has_save_in)
160 {
161 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
162 for (unsigned int i = 0; i < _save_in.size(); i++)
163 _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
164 }
165}
166
167void
168ADInertialForceShell::computeResidualsForJacobian()
169{
170 if (_residuals.size() != _test.size())
171 _residuals.resize(_test.size(), 0);
172 for (auto & r : _residuals)
173 r = 0;
174
175 precalculateResidual();
176
177 mooseAssert(_residuals.size() >= 4,
178 "This is hard coded to index from 0 to 3, so we must have at least four spots in our "
179 "container. I'd prefer to assert that the size == 4, but I don't know what the "
180 "tensor mechanics folks expect.");
181
182 if (_dt != 0.0)
183 {
184 computeShellInertialForces(_ad_coord, _ad_JxW);
185
186 if (_component < 3)
187 {
188 _global_force[0] = _thickness * _original_local_config.transpose() * _local_force[0];
189 _global_force[1] = _thickness * _original_local_config.transpose() * _local_force[1];
190 _global_force[2] = _thickness * _original_local_config.transpose() * _local_force[2];
191 _global_force[3] = _thickness * _original_local_config.transpose() * _local_force[3];
192
193 _residuals[0] = _global_force[0](_component);
194 _residuals[1] = _global_force[1](_component);
195 _residuals[2] = _global_force[2](_component);
196 _residuals[3] = _global_force[3](_component);
197 }
198 else
199 {
200 _global_moment[0] = _original_local_config.transpose() * _local_moment[0];
201 _global_moment[1] = _original_local_config.transpose() * _local_moment[1];
202 _global_moment[2] = _original_local_config.transpose() * _local_moment[2];
203 _global_moment[3] = _original_local_config.transpose() * _local_moment[3];
204
205 _residuals[0] = _global_moment[0](_component - 3);
206 _residuals[1] = _global_moment[1](_component - 3);
207 _residuals[2] = _global_moment[2](_component - 3);
208 _residuals[3] = _global_moment[3](_component - 3);
209 }
210 }
211}
212
213void
214ADInertialForceShell::computeShellInertialForces(const MooseArray<ADReal> & _ad_coord,
215 const MooseArray<ADReal> & _ad_JxW)
216{
217 // Loosely following notation in: "On finite element nonlinear analysis of general shell
218 // structures", PhD thesis by Said Bolourchi (1975).
219
220 _original_local_config = _transformation_matrix[0];
221
222 _nodes[0] = _current_elem->node_ptr(0);
223 _nodes[1] = _current_elem->node_ptr(1);
224 _nodes[2] = _current_elem->node_ptr(2);
225 _nodes[3] = _current_elem->node_ptr(3);
226
227 ADRealVectorValue x = (*_nodes[1] - *_nodes[0]);
228 ADRealVectorValue y = (*_nodes[3] - *_nodes[0]);
229 ADRealVectorValue normal = x.cross(y);
230 normal /= normal.norm();
231
232 _node_normal[0] = normal;
233 _node_normal[1] = normal;
234 _node_normal[2] = normal;
235 _node_normal[3] = normal;
236
237 // compute nodal local axis
238 for (unsigned int k = 0; k < _nodes.size(); ++k)
239 {
240 _v1[k] = _x2.cross(_node_normal[k]);
241
242 // If x2 is parallel to node normal, set V1 to x3
243 if (MooseUtils::absoluteFuzzyEqual(_v1[k].norm(), 0.0, 1e-6))
244 _v1[k] = _x3;
245
246 _v1[k] /= _v1[k].norm();
247
248 _v2[k] = _node_normal[k].cross(_v1[k]);
249 }
250
251 NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(_sys.number());
252
253 if (!nonlinear_sys.solutionUDot())
254 mooseError("InertialForceShell: Time derivative of solution (`u_dot`) is not stored. Please "
255 "set uDotRequested() to true in FEProblemBase before requesting `u_dot`.");
256
257 if (!nonlinear_sys.solutionUDotOld())
258 mooseError("InertialForceShell: Old time derivative of solution (`u_dot_old`) is not "
259 "stored. Please set uDotOldRequested() to true in FEProblemBase before "
260 "requesting `u_dot_old`.");
261
262 if (!nonlinear_sys.solutionUDotDot())
263 mooseError("InertialForceShell: Second time derivative of solution (`u_dotdot`) is not "
264 "stored. Please set uDotDotRequested() to true in FEProblemBase before "
265 "requesting `u_dotdot`.");
266
267 const NumericVector<Number> & vel = *nonlinear_sys.solutionUDot();
268 const NumericVector<Number> & old_vel = *nonlinear_sys.solutionUDotOld();
269 const NumericVector<Number> & accel = *nonlinear_sys.solutionUDotDot();
270
271 for (unsigned int i = 0; i < _ndisp; ++i)
272 {
273 // translational velocities and accelerations
274 unsigned int dof_index_0 = _nodes[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
275 unsigned int dof_index_1 = _nodes[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
276 unsigned int dof_index_2 = _nodes[2]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
277 unsigned int dof_index_3 = _nodes[3]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
278
279 _vel.pos[0](i) = vel(dof_index_0);
280 _vel.pos[1](i) = vel(dof_index_1);
281 _vel.pos[2](i) = vel(dof_index_2);
282 _vel.pos[3](i) = vel(dof_index_3);
283
284 _old_vel.pos[0](i) = old_vel(dof_index_0);
285 _old_vel.pos[1](i) = old_vel(dof_index_1);
286 _old_vel.pos[2](i) = old_vel(dof_index_2);
287 _old_vel.pos[3](i) = old_vel(dof_index_3);
288
289 _accel.pos[0](i) = accel(dof_index_0);
290 _accel.pos[1](i) = accel(dof_index_1);
291 _accel.pos[2](i) = accel(dof_index_2);
292 _accel.pos[3](i) = accel(dof_index_3);
293 }
294
295 for (unsigned int i = 0; i < _nrot; ++i)
296 {
297 // rotational velocities and accelerations
298 unsigned int dof_index_0 = _nodes[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
299 unsigned int dof_index_1 = _nodes[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
300 unsigned int dof_index_2 = _nodes[2]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
301 unsigned int dof_index_3 = _nodes[3]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
302
303 _vel.rot[0](i) = vel(dof_index_0);
304 _vel.rot[1](i) = vel(dof_index_1);
305 _vel.rot[2](i) = vel(dof_index_2);
306 _vel.rot[3](i) = vel(dof_index_3);
307
308 _old_vel.rot[0](i) = old_vel(dof_index_0);
309 _old_vel.rot[1](i) = old_vel(dof_index_1);
310 _old_vel.rot[2](i) = old_vel(dof_index_2);
311 _old_vel.rot[3](i) = old_vel(dof_index_3);
312
313 _accel.rot[0](i) = accel(dof_index_0);
314 _accel.rot[1](i) = accel(dof_index_1);
315 _accel.rot[2](i) = accel(dof_index_2);
316 _accel.rot[3](i) = accel(dof_index_3);
317 }
318 // transform translational and rotational velocities and accelerations to the initial local
319 // configuration of the shell
320 _local_vel.pos[0] = _original_local_config * _vel.pos[0];
321 _local_vel.pos[1] = _original_local_config * _vel.pos[1];
322 _local_vel.pos[2] = _original_local_config * _vel.pos[2];
323 _local_vel.pos[3] = _original_local_config * _vel.pos[3];
324
325 _local_old_vel.pos[0] = _original_local_config * _old_vel.pos[0];
326 _local_old_vel.pos[1] = _original_local_config * _old_vel.pos[1];
327 _local_old_vel.pos[2] = _original_local_config * _old_vel.pos[2];
328 _local_old_vel.pos[3] = _original_local_config * _old_vel.pos[3];
329
330 _local_accel.pos[0] = _original_local_config * _accel.pos[0];
331 _local_accel.pos[1] = _original_local_config * _accel.pos[1];
332 _local_accel.pos[2] = _original_local_config * _accel.pos[2];
333 _local_accel.pos[3] = _original_local_config * _accel.pos[3];
334
335 _local_vel.rot[0] = _original_local_config * _vel.rot[0];
336 _local_vel.rot[1] = _original_local_config * _vel.rot[1];
337 _local_vel.rot[2] = _original_local_config * _vel.rot[2];
338 _local_vel.rot[3] = _original_local_config * _vel.rot[3];
339
340 _local_old_vel.rot[0] = _original_local_config * _old_vel.rot[0];
341 _local_old_vel.rot[1] = _original_local_config * _old_vel.rot[1];
342 _local_old_vel.rot[2] = _original_local_config * _old_vel.rot[2];
343 _local_old_vel.rot[3] = _original_local_config * _old_vel.rot[3];
344
345 _local_accel.rot[0] = _original_local_config * _accel.rot[0];
346 _local_accel.rot[1] = _original_local_config * _accel.rot[1];
347 _local_accel.rot[2] = _original_local_config * _accel.rot[2];
348 _local_accel.rot[3] = _original_local_config * _accel.rot[3];
349
350 // Conversions to ADDenseVector from ADRealVectorValue: Make a method out of this.
351 ADDenseVector local_accel_dv_0(3);
352 ADDenseVector local_accel_dv_1(3);
353 ADDenseVector local_accel_dv_2(3);
354 ADDenseVector local_accel_dv_3(3);
355
356 ADDenseVector local_rot_accel_dv_0(3);
357 ADDenseVector local_rot_accel_dv_1(3);
358 ADDenseVector local_rot_accel_dv_2(3);
359 ADDenseVector local_rot_accel_dv_3(3);
360
361 ADDenseVector local_vel_dv_0(3);
362 ADDenseVector local_vel_dv_1(3);
363 ADDenseVector local_vel_dv_2(3);
364 ADDenseVector local_vel_dv_3(3);
365
366 ADDenseVector local_rot_vel_dv_0(3);
367 ADDenseVector local_rot_vel_dv_1(3);
368 ADDenseVector local_rot_vel_dv_2(3);
369 ADDenseVector local_rot_vel_dv_3(3);
370
371 ADDenseVector local_old_vel_dv_0(3);
372 ADDenseVector local_old_vel_dv_1(3);
373 ADDenseVector local_old_vel_dv_2(3);
374 ADDenseVector local_old_vel_dv_3(3);
375
376 ADDenseVector local_old_rot_vel_dv_0(3);
377 ADDenseVector local_old_rot_vel_dv_1(3);
378 ADDenseVector local_old_rot_vel_dv_2(3);
379 ADDenseVector local_old_rot_vel_dv_3(3);
380
381 for (unsigned int i = 0; i < 3; i++)
382 {
383 local_accel_dv_0(i) = _local_accel.pos[0](i);
384 local_accel_dv_1(i) = _local_accel.pos[1](i);
385 local_accel_dv_2(i) = _local_accel.pos[2](i);
386 local_accel_dv_3(i) = _local_accel.pos[3](i);
387
388 local_rot_accel_dv_0(i) = _local_accel.rot[0](i);
389 local_rot_accel_dv_1(i) = _local_accel.rot[1](i);
390 local_rot_accel_dv_2(i) = _local_accel.rot[2](i);
391 local_rot_accel_dv_3(i) = _local_accel.rot[3](i);
392
393 local_vel_dv_0(i) = _local_vel.pos[0](i);
394 local_vel_dv_1(i) = _local_vel.pos[1](i);
395 local_vel_dv_2(i) = _local_vel.pos[2](i);
396 local_vel_dv_3(i) = _local_vel.pos[3](i);
397
398 local_rot_vel_dv_0(i) = _local_vel.rot[0](i);
399 local_rot_vel_dv_1(i) = _local_vel.rot[1](i);
400 local_rot_vel_dv_2(i) = _local_vel.rot[2](i);
401 local_rot_vel_dv_3(i) = _local_vel.rot[3](i);
402
403 local_old_vel_dv_0(i) = _local_old_vel.pos[0](i);
404 local_old_vel_dv_1(i) = _local_old_vel.pos[1](i);
405 local_old_vel_dv_2(i) = _local_old_vel.pos[2](i);
406 local_old_vel_dv_3(i) = _local_old_vel.pos[3](i);
407
408 local_old_rot_vel_dv_0(i) = _local_old_vel.rot[0](i);
409 local_old_rot_vel_dv_1(i) = _local_old_vel.rot[1](i);
410 local_old_rot_vel_dv_2(i) = _local_old_vel.rot[2](i);
411 local_old_rot_vel_dv_3(i) = _local_old_vel.rot[3](i);
412 }
413 unsigned int dim = _current_elem->dim();
414
415 // Update 0g vectors at plane quadrature points.
416 _0g1_vectors[0] = -0.5 * _thickness * _v1[0];
417 _0g1_vectors[1] = 0.5 * _thickness * _v2[0];
418
419 ADDenseMatrix G1(3, 2);
420 G1(0, 0) = _0g1_vectors[0](0);
421 G1(1, 0) = _0g1_vectors[0](1);
422 G1(2, 0) = _0g1_vectors[0](2);
423 G1(0, 1) = _0g1_vectors[1](0);
424 G1(1, 1) = _0g1_vectors[1](1);
425 G1(2, 1) = _0g1_vectors[1](2);
426 ADDenseMatrix G1T(2, 3);
427 G1.get_transpose(G1T);
428
429 _0g2_vectors[0] = -0.5 * _thickness * _v1[1];
430 _0g2_vectors[1] = 0.5 * _thickness * _v2[1];
431
432 ADDenseMatrix G2(3, 2);
433 G2(0, 0) = _0g2_vectors[0](0);
434 G2(1, 0) = _0g2_vectors[0](1);
435 G2(2, 0) = _0g2_vectors[0](2);
436 G2(0, 1) = _0g2_vectors[1](0);
437 G2(1, 1) = _0g2_vectors[1](1);
438 G2(2, 1) = _0g2_vectors[1](2);
439
440 ADDenseMatrix G2T(2, 3);
441 G2.get_transpose(G2T);
442
443 _0g3_vectors[0] = -0.5 * _thickness * _v1[2];
444 _0g3_vectors[1] = 0.5 * _thickness * _v2[2];
445
446 ADDenseMatrix G3(3, 2);
447 G3(0, 0) = _0g3_vectors[0](0);
448 G3(1, 0) = _0g3_vectors[0](1);
449 G3(2, 0) = _0g3_vectors[0](2);
450 G3(0, 1) = _0g3_vectors[1](0);
451 G3(1, 1) = _0g3_vectors[1](1);
452 G3(2, 1) = _0g3_vectors[1](2);
453
454 ADDenseMatrix G3T(2, 3);
455 G3.get_transpose(G3T);
456
457 _0g4_vectors[0] = -0.5 * _thickness * _v1[3];
458 _0g4_vectors[1] = 0.5 * _thickness * _v2[3];
459
460 ADDenseMatrix G4(3, 2);
461 G4(0, 0) = _0g4_vectors[0](0);
462 G4(1, 0) = _0g4_vectors[0](1);
463 G4(2, 0) = _0g4_vectors[0](2);
464 G4(0, 1) = _0g4_vectors[1](0);
465 G4(1, 1) = _0g4_vectors[1](1);
466 G4(2, 1) = _0g4_vectors[1](2);
467
468 ADDenseMatrix G4T(2, 3);
469 G4.get_transpose(G4T);
470
471 std::vector<ADDenseVector> local_acc;
472 local_acc.resize(4);
473 local_acc[0].resize(3);
474 local_acc[1].resize(3);
475 local_acc[2].resize(3);
476 local_acc[3].resize(3);
477
478 local_acc[0] = local_accel_dv_0;
479 local_acc[1] = local_accel_dv_1;
480 local_acc[2] = local_accel_dv_2;
481 local_acc[3] = local_accel_dv_3;
482
483 std::vector<ADDenseVector> local_rot_acc;
484 local_rot_acc.resize(4);
485 local_rot_acc[0].resize(3);
486 local_rot_acc[1].resize(3);
487 local_rot_acc[2].resize(3);
488 local_rot_acc[3].resize(3);
489 local_rot_acc[0] = local_rot_accel_dv_0;
490 local_rot_acc[1] = local_rot_accel_dv_1;
491 local_rot_acc[2] = local_rot_accel_dv_2;
492 local_rot_acc[3] = local_rot_accel_dv_3;
493
494 // Velocity for Rayleigh damping, including HHT_alpha parameter
495 // {
496 std::vector<ADDenseVector> local_vel;
497 local_vel.resize(4);
498 local_vel[0].resize(3);
499 local_vel[1].resize(3);
500 local_vel[2].resize(3);
501 local_vel[3].resize(3);
502
503 local_vel_dv_0.scale(1 + _alpha);
504 local_old_vel_dv_0.scale(_alpha);
505 local_vel_dv_1.scale(1 + _alpha);
506 local_old_vel_dv_1.scale(_alpha);
507 local_vel_dv_2.scale(1 + _alpha);
508 local_old_vel_dv_2.scale(_alpha);
509 local_vel_dv_3.scale(1 + _alpha);
510 local_old_vel_dv_3.scale(_alpha);
511
512 local_vel_dv_0.add(1.0, local_old_vel_dv_0);
513 local_vel_dv_1.add(1.0, local_old_vel_dv_1);
514 local_vel_dv_2.add(1.0, local_old_vel_dv_2);
515 local_vel_dv_3.add(1.0, local_old_vel_dv_3);
516
517 local_vel[0] = local_vel_dv_0;
518 local_vel[1] = local_vel_dv_1;
519 local_vel[2] = local_vel_dv_2;
520 local_vel[3] = local_vel_dv_3;
521
522 std::vector<ADDenseVector> local_rot_vel;
523 local_rot_vel.resize(4);
524 local_rot_vel[0].resize(3);
525 local_rot_vel[1].resize(3);
526 local_rot_vel[2].resize(3);
527 local_rot_vel[3].resize(3);
528
529 local_rot_vel_dv_0.scale(1 + _alpha);
530 local_old_rot_vel_dv_0.scale(_alpha);
531 local_rot_vel_dv_1.scale(1 + _alpha);
532 local_old_rot_vel_dv_1.scale(_alpha);
533 local_rot_vel_dv_2.scale(1 + _alpha);
534 local_old_rot_vel_dv_2.scale(_alpha);
535 local_rot_vel_dv_3.scale(1 + _alpha);
536 local_old_rot_vel_dv_3.scale(_alpha);
537
538 local_rot_vel_dv_0.add(1.0, local_old_rot_vel_dv_0);
539 local_rot_vel_dv_1.add(1.0, local_old_rot_vel_dv_1);
540 local_rot_vel_dv_2.add(1.0, local_old_rot_vel_dv_2);
541 local_rot_vel_dv_3.add(1.0, local_old_rot_vel_dv_3);
542
543 local_rot_vel[0] = local_rot_vel_dv_0;
544 local_rot_vel[1] = local_rot_vel_dv_1;
545 local_rot_vel[2] = local_rot_vel_dv_2;
546 local_rot_vel[3] = local_rot_vel_dv_3;
547 // }
548
549 FEType fe_type(Utility::string_to_enum<Order>("First"),
550 Utility::string_to_enum<FEFamily>("LAGRANGE"));
551 auto & fe = _fe_problem.assembly(_tid, _sys.number()).getFE(fe_type, dim);
552 _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
553 _dphideta_map = fe->get_fe_map().get_dphideta_map();
554 _phi_map = fe->get_fe_map().get_phi_map();
555
556 // quadrature points in isoparametric space
557 _2d_points = _qrule->get_points(); // would be in 2D
558
559 std::vector<const Node *> nodes;
560 for (unsigned int i = 0; i < 4; ++i)
561 nodes.push_back(_current_elem->node_ptr(i));
562
563 for (unsigned int i = 0; i < _ndisp; i++)
564 for (unsigned int j = 0; j < 4; j++)
565 _local_force[j](i) = 0.0;
566
567 for (unsigned int i = 0; i < _nrot; i++)
568 for (unsigned int j = 0; j < 4; j++)
569 _local_moment[j](i) = 0.0;
570
571 for (unsigned int qp_xy = 0; qp_xy < _2d_points.size(); ++qp_xy)
572 {
573 ADReal factor_qxy = _ad_coord[qp_xy] * _ad_JxW[qp_xy] * _density[qp_xy];
574
575 // Account for inertia on displacement degrees of freedom
576 for (unsigned int dim = 0; dim < 3; dim++)
577 {
578 _local_force[0](dim) +=
579 factor_qxy * (_phi_map[0][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
580 _phi_map[0][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
581 _phi_map[0][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
582 _phi_map[0][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
583
584 if (_eta[0] > TOLERANCE * TOLERANCE)
585 _local_force[0](dim) += factor_qxy * _eta[0] *
586 (_phi_map[0][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
587 _phi_map[0][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
588 _phi_map[0][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
589 _phi_map[0][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
590
591 _local_force[1](dim) +=
592 factor_qxy * (_phi_map[1][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
593 _phi_map[1][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
594 _phi_map[1][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
595 _phi_map[1][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
596
597 if (_eta[0] > TOLERANCE * TOLERANCE)
598 _local_force[1](dim) += factor_qxy * _eta[0] *
599 (_phi_map[1][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
600 _phi_map[1][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
601 _phi_map[1][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
602 _phi_map[1][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
603 _local_force[2](dim) +=
604 factor_qxy * (_phi_map[2][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
605 _phi_map[2][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
606 _phi_map[2][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
607 _phi_map[2][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
608
609 if (_eta[0] > TOLERANCE * TOLERANCE)
610 _local_force[2](dim) += factor_qxy * _eta[0] *
611 (_phi_map[2][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
612 _phi_map[2][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
613 _phi_map[2][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
614 _phi_map[2][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
615
616 _local_force[3](dim) +=
617 factor_qxy * (_phi_map[3][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
618 _phi_map[3][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
619 _phi_map[3][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
620 _phi_map[3][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
621
622 if (_eta[0] > TOLERANCE * TOLERANCE)
623 _local_force[3](dim) += factor_qxy * _eta[0] *
624 (_phi_map[3][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
625 _phi_map[3][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
626 _phi_map[3][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
627 _phi_map[3][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
628 }
629
630 // Account for inertia on rotational degrees of freedom
631 ADReal rot_thickness = _thickness * _thickness * _thickness / 48.0;
632
633 ADDenseVector momentInertia(3);
634 momentInertia(0) = (G1(0, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
635 G1(0, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
636 G2(0, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
637 G2(0, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
638 G3(0, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
639 G3(0, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
640 G4(0, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
641 G4(0, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
642
643 momentInertia(1) = (G1(1, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
644 G1(1, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
645 G2(1, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
646 G2(1, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
647 G3(1, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
648 G3(1, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
649 G4(1, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
650 G4(1, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
651
652 momentInertia(2) = (G1(2, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
653 G1(2, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
654 G2(2, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
655 G2(2, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
656 G3(2, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
657 G3(2, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
658 G4(2, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
659 G4(2, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
660
661 _local_moment[0](0) += factor_qxy * rot_thickness *
662 (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
663 G1T(0, 2) * momentInertia(2));
664
665 _local_moment[0](1) += factor_qxy * rot_thickness *
666 (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
667 G1T(1, 2) * momentInertia(2));
668
669 _local_moment[1](0) += factor_qxy * rot_thickness *
670 (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
671 G2T(0, 2) * momentInertia(2));
672
673 _local_moment[1](1) += factor_qxy * rot_thickness *
674 (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
675 G2T(1, 2) * momentInertia(2));
676
677 _local_moment[2](0) += factor_qxy * rot_thickness *
678 (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
679 G3T(0, 2) * momentInertia(2));
680
681 _local_moment[2](1) += factor_qxy * rot_thickness *
682 (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
683 G3T(1, 2) * momentInertia(2));
684
685 _local_moment[3](0) += factor_qxy * rot_thickness *
686 (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
687 G4T(0, 2) * momentInertia(2));
688
689 _local_moment[3](1) += factor_qxy * rot_thickness *
690 (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
691 G4T(1, 2) * momentInertia(2));
692 }
693}
registerMooseObject
registerMooseObject("SolidMechanicsApp", ADInertialForceShell)
ADReal
DualNumber< Real, DNDerivativeType, true > ADReal
y
const std::vector< double > y
Definition EquilibriumConstantFitTest.C:18
x
const std::vector< double > x
Definition EquilibriumConstantFitTest.C:17
dim
unsigned int dim
int
void ErrorVector unsigned int
ADInertialForceShell::_component
const unsigned int _component
Direction along which residual is calculated.
Definition ADInertialForceShell.h:75
ADInertialForceShell::_dphideta_map
std::vector< std::vector< Real > > _dphideta_map
Derivatives of shape functions w.r.t isoparametric coordinates eta.
Definition ADInertialForceShell.h:111
ADInertialForceShell::_nrot
unsigned int _nrot
Number of coupled rotational variables.
Definition ADInertialForceShell.h:45
ADInertialForceShell::_v1
std::vector< ADRealVectorValue > _v1
First tangential vectors at nodes.
Definition ADInertialForceShell.h:126
ADInertialForceShell::_local_force
std::array< ADRealVectorValue, 4 > _local_force
Forces and moments at the four nodes in the initial local configuration.
Definition ADInertialForceShell.h:99
ADInertialForceShell::computeShellInertialForces
virtual void computeShellInertialForces(const MooseArray< ADReal > &_ad_coord, const MooseArray< ADReal > &_ad_JxW)
Definition ADInertialForceShell.C:214
ADInertialForceShell::_global_moment
std::array< ADRealVectorValue, 4 > _global_moment
Definition ADInertialForceShell.h:105
ADInertialForceShell::_local_accel
PosRotVectors _local_accel
Current shell nodal accelerations in the local frame of reference.
Definition ADInertialForceShell.h:93
ADInertialForceShell::_local_moment
std::array< ADRealVectorValue, 4 > _local_moment
Definition ADInertialForceShell.h:99
ADInertialForceShell::_transformation_matrix
const ADMaterialProperty< RankTwoTensor > & _transformation_matrix
Rotation matrix material property.
Definition ADInertialForceShell.h:156
ADInertialForceShell::computeResidual
virtual void computeResidual() override
Definition ADInertialForceShell.C:119
ADInertialForceShell::_thickness
const ADReal _thickness
Coupled variable for the shell thickness.
Definition ADInertialForceShell.h:162
ADInertialForceShell::computeResidualsForJacobian
virtual void computeResidualsForJacobian() override
Definition ADInertialForceShell.C:168
ADInertialForceShell::_ndisp
unsigned int _ndisp
Number of coupled displacement variables.
Definition ADInertialForceShell.h:48
ADInertialForceShell::_x3
ADRealVectorValue _x3
Helper vector.
Definition ADInertialForceShell.h:135
ADInertialForceShell::_global_force
std::array< ADRealVectorValue, 4 > _global_force
Forces and moments at the four nodes in the global coordinate system.
Definition ADInertialForceShell.h:105
ADInertialForceShell::_density
const MaterialProperty< Real > & _density
Shell material density.
Definition ADInertialForceShell.h:165
ADInertialForceShell::_phi_map
std::vector< std::vector< Real > > _phi_map
Shape function value.
Definition ADInertialForceShell.h:114
ADInertialForceShell::ADInertialForceShell
ADInertialForceShell(const InputParameters &parameters)
Definition ADInertialForceShell.C:75
ADInertialForceShell::_0g2_vectors
std::array< ADRealVectorValue, 2 > _0g2_vectors
Node 2 g vectors.
Definition ADInertialForceShell.h:144
ADInertialForceShell::_2d_points
std::vector< Point > _2d_points
Quadrature points in the in-plane direction in isoparametric coordinate system.
Definition ADInertialForceShell.h:120
ADInertialForceShell::_dphidxi_map
std::vector< std::vector< Real > > _dphidxi_map
Derivatives of shape functions w.r.t isoparametric coordinates xi.
Definition ADInertialForceShell.h:108
ADInertialForceShell::validParams
static InputParameters validParams()
Definition ADInertialForceShell.C:32
ADInertialForceShell::_local_vel
PosRotVectors _local_vel
Current shell nodal velocities in the local frame of reference.
Definition ADInertialForceShell.h:87
ADInertialForceShell::_local_old_vel
PosRotVectors _local_old_vel
Old shell nodal velocities in the local frame of reference.
Definition ADInertialForceShell.h:90
ADInertialForceShell::_node_normal
std::vector< ADRealVectorValue > _node_normal
Normal to the element at the 4 nodes.
Definition ADInertialForceShell.h:138
ADInertialForceShell::_v2
std::vector< ADRealVectorValue > _v2
Second tangential vectors at nodes.
Definition ADInertialForceShell.h:129
ADInertialForceShell::_x2
ADRealVectorValue _x2
Helper vector.
Definition ADInertialForceShell.h:132
ADInertialForceShell::_old_vel
PosRotVectors _old_vel
Old shell nodal velocities in the global frame of reference.
Definition ADInertialForceShell.h:81
ADInertialForceShell::_accel
PosRotVectors _accel
Current shell nodal accelerations in the global frame of reference.
Definition ADInertialForceShell.h:84
ADInertialForceShell::_0g3_vectors
std::array< ADRealVectorValue, 2 > _0g3_vectors
Node 3 g vectors.
Definition ADInertialForceShell.h:147
ADInertialForceShell::_rot_num
std::vector< unsigned int > _rot_num
Variable numbers corresponding to rotational variables.
Definition ADInertialForceShell.h:51
ADInertialForceShell::_original_local_config
ADRankTwoTensor _original_local_config
Rotational transformation from global to initial shell local coordinate system.
Definition ADInertialForceShell.h:72
ADInertialForceShell::_eta
const MaterialProperty< Real > & _eta
Mass proportional Rayleigh damping parameter.
Definition ADInertialForceShell.h:153
ADInertialForceShell::_0g1_vectors
std::array< ADRealVectorValue, 2 > _0g1_vectors
Node 1 g vectors.
Definition ADInertialForceShell.h:141
ADInertialForceShell::_0g4_vectors
std::array< ADRealVectorValue, 2 > _0g4_vectors
Node 4 g vectors.
Definition ADInertialForceShell.h:150
ADInertialForceShell::_vel
PosRotVectors _vel
Current shell nodal velocities in the global frame of reference.
Definition ADInertialForceShell.h:78
ADInertialForceShell::_disp_num
std::vector< unsigned int > _disp_num
Variable numbers corresponding to displacement variables.
Definition ADInertialForceShell.h:54
ADInertialForceShell::_alpha
const Real _alpha
HHT time integration parameter.
Definition ADInertialForceShell.h:168
ADInertialForceShell::_nodes
std::vector< const Node * > _nodes
Vector storing pointers to the nodes of the shell element.
Definition ADInertialForceShell.h:117
ADTimeKernelTempl::_test
const ADTemplateVariableTestValue< T > & _test
ADTimeKernelTempl::_residuals
std::vector< ADReal > _residuals
ADTimeKernelTempl::_var
MooseVariableFE< T > & _var
ADTimeKernelTempl::_ad_JxW
const MooseArray< ADReal > & _ad_JxW
ADTimeKernelTempl::dofIndices
virtual const std::vector< dof_id_type > & dofIndices() const
ADTimeKernelTempl::_ad_coord
const MooseArray< ADReal > & _ad_coord
ADTimeKernelTempl::validParams
static InputParameters validParams()
Assembly::getFE
const FEBase *const & getFE(FEType type, unsigned int dim) const
ADTimeKernelTempl::getVar
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
FEProblemBase::assembly
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override
FEProblemBase::getNonlinearSystemBase
NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num)
InputParameters::addRequiredRangeCheckedParam
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
InputParameters::addRequiredCoupledVar
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
InputParameters::addRequiredParam
void addRequiredParam(const std::string &name, const std::string &doc_string)
InputParameters::addParam
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
InputParameters::addClassDescription
void addClassDescription(const std::string &doc_string)
InputParameters::set
T & set(const std::string &name, bool quiet_mode=false)
InputParameters::addCoupledVar
void addCoupledVar(const std::string &name, const std::string &doc_string)
InputParameters::addRangeCheckedParam
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
ADTimeKernelTempl::_save_in
std::vector< MooseVariableFEBase * > _save_in
ADTimeKernelTempl::_current_elem
const Elem *const & _current_elem
ADTimeKernelTempl::_qrule
const QBase *const & _qrule
ADTimeKernelTempl::mooseError
void mooseError(Args &&... args) const
MooseVariableFE::number
unsigned int number() const
RankTwoTensorTempl::transpose
RankTwoTensorTempl< T > transpose() const
ADTimeKernelTempl::_assembly
Assembly & _assembly
ADTimeKernelTempl::_sys
SystemBase & _sys
ADTimeKernelTempl::precalculateResidual
virtual void precalculateResidual()
ADTimeKernelTempl::_fe_problem
FEProblemBase & _fe_problem
NonlinearSystemBase::solutionUDot
virtual NumericVector< Number > * solutionUDot()
SystemBase::number
unsigned int number() const
NonlinearSystemBase::solutionUDotOld
virtual NumericVector< Number > * solutionUDotOld()
NonlinearSystemBase::solutionUDotDot
virtual NumericVector< Number > * solutionUDotDot()
ADTimeKernelTempl::accumulateTaggedLocalResidual
void accumulateTaggedLocalResidual()
ADTimeKernelTempl::prepareVectorTag
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
ADTimeKernelTempl::_local_re
DenseVector< Number > _local_re
libMesh::DenseMatrix::get_transpose
void get_transpose(DenseMatrix< T > &dest) const
libMesh::DenseVector::resize
void resize(const unsigned int n)
libMesh::DenseVector::add
std::enable_if< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
libMesh::DenseVector::scale
void scale(const T factor)
libMesh::VectorValue::cross
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
PosRotVectors::rot
std::array< ADRealVectorValue, 4 > rot
Definition ADInertialForceShell.h:25
PosRotVectors::pos
std::array< ADRealVectorValue, 4 > pos
Definition ADInertialForceShell.h:24
Generated on Fri Aug 21 2026 13:41:38 for https://mooseframework.inl.gov by doxygen 1.9.8