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 "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 :
29 : registerMooseObject("SolidMechanicsApp", ADInertialForceShell);
30 :
31 : InputParameters
32 70 : ADInertialForceShell::validParams()
33 : {
34 70 : InputParameters params = ADTimeKernel::validParams();
35 :
36 70 : 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 70 : params.set<bool>("use_displaced_mesh") = true;
40 140 : params.addRequiredCoupledVar(
41 : "rotations",
42 : "The rotational variables appropriate for the simulation geometry and coordinate system");
43 140 : params.addRequiredCoupledVar(
44 : "displacements",
45 : "The displacement variables appropriate for the simulation geometry and coordinate system");
46 140 : params.addCoupledVar("velocities", "Translational velocity variables");
47 140 : params.addCoupledVar("accelerations", "Translational acceleration variables");
48 140 : params.addCoupledVar("rotational_velocities", "Rotational velocity variables");
49 140 : params.addCoupledVar("rotational_accelerations", "Rotational acceleration variables");
50 140 : params.addParam<MaterialPropertyName>(
51 : "density",
52 : "density",
53 : "Name of Material Property or a constant real number defining the density of the beam.");
54 140 : 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 140 : params.addRequiredParam<Real>("thickness", "The kernel's thickness");
61 140 : params.addParam<MaterialPropertyName>("eta",
62 140 : 0.0,
63 : "Name of material property or a constant real "
64 : "number defining the eta parameter for the "
65 : "Rayleigh damping.");
66 210 : params.addRangeCheckedParam<Real>("alpha",
67 140 : 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 70 : return params;
73 0 : }
74 :
75 35 : ADInertialForceShell::ADInertialForceShell(const InputParameters & parameters)
76 : : ADTimeKernel(parameters),
77 35 : _nrot(coupledComponents("rotations")),
78 35 : _ndisp(coupledComponents("displacements")),
79 35 : _rot_num(_nrot),
80 35 : _disp_num(_ndisp),
81 35 : _vel_num(_ndisp),
82 35 : _accel_num(_ndisp),
83 35 : _rot_vel_num(_nrot),
84 35 : _rot_accel_num(_nrot),
85 70 : _component(getParam<unsigned int>("component")),
86 35 : _nodes(4),
87 35 : _v1(4),
88 35 : _v2(4),
89 35 : _node_normal(4),
90 70 : _eta(getMaterialProperty<Real>("eta")),
91 70 : _transformation_matrix(getADMaterialProperty<RankTwoTensor>("transformation_matrix_element")),
92 70 : _J_map(getADMaterialProperty<Real>("J_mapping_t_points_0")),
93 70 : _thickness(getParam<Real>("thickness")),
94 70 : _density(getMaterialProperty<Real>("density")),
95 210 : _alpha(getParam<Real>("alpha"))
96 : {
97 : // Checking for consistency between the length of the provided rotations and displacements vector
98 35 : if (_ndisp != 3 || _nrot != 2)
99 0 : 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 140 : for (unsigned int i = 0; i < _ndisp; ++i)
104 : {
105 210 : MooseVariable * disp_variable = getVar("displacements", i);
106 105 : _disp_num[i] = disp_variable->number();
107 : }
108 105 : for (unsigned int i = 0; i < _nrot; ++i)
109 : {
110 140 : MooseVariable * rot_variable = getVar("rotations", i);
111 70 : _rot_num[i] = rot_variable->number();
112 : }
113 :
114 35 : _x2(1) = 1;
115 35 : _x3(2) = 1;
116 35 : }
117 :
118 : void
119 3640 : ADInertialForceShell::computeResidual()
120 : {
121 3640 : prepareVectorTag(_assembly, _var.number());
122 :
123 3640 : precalculateResidual();
124 :
125 : /* ----------------------------------------------------- */
126 :
127 3640 : if (_dt != 0.0)
128 : {
129 3640 : computeShellInertialForces(_ad_coord, _ad_JxW);
130 :
131 3640 : if (_component < 3)
132 : {
133 2184 : _global_force[0] = _thickness * _original_local_config.transpose() * _local_force[0];
134 4368 : _global_force[1] = _thickness * _original_local_config.transpose() * _local_force[1];
135 4368 : _global_force[2] = _thickness * _original_local_config.transpose() * _local_force[2];
136 4368 : _global_force[3] = _thickness * _original_local_config.transpose() * _local_force[3];
137 :
138 2184 : _local_re(0) = raw_value(_global_force[0](_component));
139 2184 : _local_re(1) = raw_value(_global_force[1](_component));
140 2184 : _local_re(2) = raw_value(_global_force[2](_component));
141 2184 : _local_re(3) = raw_value(_global_force[3](_component));
142 : }
143 : else
144 : {
145 1456 : _global_moment[0] = _original_local_config.transpose() * _local_moment[0];
146 1456 : _global_moment[1] = _original_local_config.transpose() * _local_moment[1];
147 1456 : _global_moment[2] = _original_local_config.transpose() * _local_moment[2];
148 1456 : _global_moment[3] = _original_local_config.transpose() * _local_moment[3];
149 :
150 1456 : _local_re(0) = raw_value(_global_moment[0](_component - 3));
151 1456 : _local_re(1) = raw_value(_global_moment[1](_component - 3));
152 1456 : _local_re(2) = raw_value(_global_moment[2](_component - 3));
153 1456 : _local_re(3) = raw_value(_global_moment[3](_component - 3));
154 : }
155 : }
156 :
157 3640 : accumulateTaggedLocalResidual();
158 :
159 3640 : if (_has_save_in)
160 : {
161 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
162 0 : for (unsigned int i = 0; i < _save_in.size(); i++)
163 0 : _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
164 : }
165 3640 : }
166 :
167 : void
168 520 : ADInertialForceShell::computeResidualsForJacobian()
169 : {
170 520 : if (_residuals.size() != _test.size())
171 35 : _residuals.resize(_test.size(), 0);
172 2600 : for (auto & r : _residuals)
173 2080 : r = 0;
174 :
175 520 : 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 520 : if (_dt != 0.0)
183 : {
184 520 : computeShellInertialForces(_ad_coord, _ad_JxW);
185 :
186 520 : if (_component < 3)
187 : {
188 312 : _global_force[0] = _thickness * _original_local_config.transpose() * _local_force[0];
189 624 : _global_force[1] = _thickness * _original_local_config.transpose() * _local_force[1];
190 624 : _global_force[2] = _thickness * _original_local_config.transpose() * _local_force[2];
191 624 : _global_force[3] = _thickness * _original_local_config.transpose() * _local_force[3];
192 :
193 312 : _residuals[0] = _global_force[0](_component);
194 312 : _residuals[1] = _global_force[1](_component);
195 312 : _residuals[2] = _global_force[2](_component);
196 312 : _residuals[3] = _global_force[3](_component);
197 : }
198 : else
199 : {
200 208 : _global_moment[0] = _original_local_config.transpose() * _local_moment[0];
201 208 : _global_moment[1] = _original_local_config.transpose() * _local_moment[1];
202 208 : _global_moment[2] = _original_local_config.transpose() * _local_moment[2];
203 208 : _global_moment[3] = _original_local_config.transpose() * _local_moment[3];
204 :
205 208 : _residuals[0] = _global_moment[0](_component - 3);
206 208 : _residuals[1] = _global_moment[1](_component - 3);
207 208 : _residuals[2] = _global_moment[2](_component - 3);
208 208 : _residuals[3] = _global_moment[3](_component - 3);
209 : }
210 : }
211 520 : }
212 :
213 : void
214 4160 : ADInertialForceShell::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 4160 : _original_local_config = _transformation_matrix[0];
221 :
222 4160 : _nodes[0] = _current_elem->node_ptr(0);
223 4160 : _nodes[1] = _current_elem->node_ptr(1);
224 4160 : _nodes[2] = _current_elem->node_ptr(2);
225 4160 : _nodes[3] = _current_elem->node_ptr(3);
226 :
227 4160 : ADRealVectorValue x = (*_nodes[1] - *_nodes[0]);
228 4160 : ADRealVectorValue y = (*_nodes[3] - *_nodes[0]);
229 4160 : ADRealVectorValue normal = x.cross(y);
230 8320 : 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 20800 : for (unsigned int k = 0; k < _nodes.size(); ++k)
239 : {
240 33280 : _v1[k] = _x2.cross(_node_normal[k]);
241 :
242 : // If x2 is parallel to node normal, set V1 to x3
243 16640 : if (MooseUtils::absoluteFuzzyEqual(_v1[k].norm(), 0.0, 1e-6))
244 : _v1[k] = _x3;
245 :
246 33280 : _v1[k] /= _v1[k].norm();
247 :
248 33280 : _v2[k] = _node_normal[k].cross(_v1[k]);
249 : }
250 :
251 4160 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(_sys.number());
252 :
253 4160 : if (!nonlinear_sys.solutionUDot())
254 0 : 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 4160 : if (!nonlinear_sys.solutionUDotOld())
258 0 : 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 4160 : if (!nonlinear_sys.solutionUDotDot())
263 0 : 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 4160 : const NumericVector<Number> & vel = *nonlinear_sys.solutionUDot();
268 4160 : const NumericVector<Number> & old_vel = *nonlinear_sys.solutionUDotOld();
269 4160 : const NumericVector<Number> & accel = *nonlinear_sys.solutionUDotDot();
270 :
271 16640 : for (unsigned int i = 0; i < _ndisp; ++i)
272 : {
273 : // translational velocities and accelerations
274 12480 : unsigned int dof_index_0 = _nodes[0]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
275 12480 : unsigned int dof_index_1 = _nodes[1]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
276 12480 : unsigned int dof_index_2 = _nodes[2]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
277 12480 : unsigned int dof_index_3 = _nodes[3]->dof_number(nonlinear_sys.number(), _disp_num[i], 0);
278 :
279 12480 : _vel.pos[0](i) = vel(dof_index_0);
280 12480 : _vel.pos[1](i) = vel(dof_index_1);
281 12480 : _vel.pos[2](i) = vel(dof_index_2);
282 12480 : _vel.pos[3](i) = vel(dof_index_3);
283 :
284 12480 : _old_vel.pos[0](i) = old_vel(dof_index_0);
285 12480 : _old_vel.pos[1](i) = old_vel(dof_index_1);
286 12480 : _old_vel.pos[2](i) = old_vel(dof_index_2);
287 12480 : _old_vel.pos[3](i) = old_vel(dof_index_3);
288 :
289 12480 : _accel.pos[0](i) = accel(dof_index_0);
290 12480 : _accel.pos[1](i) = accel(dof_index_1);
291 12480 : _accel.pos[2](i) = accel(dof_index_2);
292 12480 : _accel.pos[3](i) = accel(dof_index_3);
293 : }
294 :
295 12480 : for (unsigned int i = 0; i < _nrot; ++i)
296 : {
297 : // rotational velocities and accelerations
298 8320 : unsigned int dof_index_0 = _nodes[0]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
299 8320 : unsigned int dof_index_1 = _nodes[1]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
300 8320 : unsigned int dof_index_2 = _nodes[2]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
301 8320 : unsigned int dof_index_3 = _nodes[3]->dof_number(nonlinear_sys.number(), _rot_num[i], 0);
302 :
303 8320 : _vel.rot[0](i) = vel(dof_index_0);
304 8320 : _vel.rot[1](i) = vel(dof_index_1);
305 8320 : _vel.rot[2](i) = vel(dof_index_2);
306 8320 : _vel.rot[3](i) = vel(dof_index_3);
307 :
308 8320 : _old_vel.rot[0](i) = old_vel(dof_index_0);
309 8320 : _old_vel.rot[1](i) = old_vel(dof_index_1);
310 8320 : _old_vel.rot[2](i) = old_vel(dof_index_2);
311 8320 : _old_vel.rot[3](i) = old_vel(dof_index_3);
312 :
313 8320 : _accel.rot[0](i) = accel(dof_index_0);
314 8320 : _accel.rot[1](i) = accel(dof_index_1);
315 8320 : _accel.rot[2](i) = accel(dof_index_2);
316 8320 : _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 4160 : _local_vel.pos[0] = _original_local_config * _vel.pos[0];
321 4160 : _local_vel.pos[1] = _original_local_config * _vel.pos[1];
322 4160 : _local_vel.pos[2] = _original_local_config * _vel.pos[2];
323 4160 : _local_vel.pos[3] = _original_local_config * _vel.pos[3];
324 :
325 4160 : _local_old_vel.pos[0] = _original_local_config * _old_vel.pos[0];
326 4160 : _local_old_vel.pos[1] = _original_local_config * _old_vel.pos[1];
327 4160 : _local_old_vel.pos[2] = _original_local_config * _old_vel.pos[2];
328 4160 : _local_old_vel.pos[3] = _original_local_config * _old_vel.pos[3];
329 :
330 4160 : _local_accel.pos[0] = _original_local_config * _accel.pos[0];
331 4160 : _local_accel.pos[1] = _original_local_config * _accel.pos[1];
332 4160 : _local_accel.pos[2] = _original_local_config * _accel.pos[2];
333 4160 : _local_accel.pos[3] = _original_local_config * _accel.pos[3];
334 :
335 4160 : _local_vel.rot[0] = _original_local_config * _vel.rot[0];
336 4160 : _local_vel.rot[1] = _original_local_config * _vel.rot[1];
337 4160 : _local_vel.rot[2] = _original_local_config * _vel.rot[2];
338 4160 : _local_vel.rot[3] = _original_local_config * _vel.rot[3];
339 :
340 4160 : _local_old_vel.rot[0] = _original_local_config * _old_vel.rot[0];
341 4160 : _local_old_vel.rot[1] = _original_local_config * _old_vel.rot[1];
342 4160 : _local_old_vel.rot[2] = _original_local_config * _old_vel.rot[2];
343 4160 : _local_old_vel.rot[3] = _original_local_config * _old_vel.rot[3];
344 :
345 4160 : _local_accel.rot[0] = _original_local_config * _accel.rot[0];
346 4160 : _local_accel.rot[1] = _original_local_config * _accel.rot[1];
347 4160 : _local_accel.rot[2] = _original_local_config * _accel.rot[2];
348 4160 : _local_accel.rot[3] = _original_local_config * _accel.rot[3];
349 :
350 : // Conversions to ADDenseVector from ADRealVectorValue: Make a method out of this.
351 4160 : ADDenseVector local_accel_dv_0(3);
352 4160 : ADDenseVector local_accel_dv_1(3);
353 4160 : ADDenseVector local_accel_dv_2(3);
354 4160 : ADDenseVector local_accel_dv_3(3);
355 :
356 4160 : ADDenseVector local_rot_accel_dv_0(3);
357 4160 : ADDenseVector local_rot_accel_dv_1(3);
358 4160 : ADDenseVector local_rot_accel_dv_2(3);
359 4160 : ADDenseVector local_rot_accel_dv_3(3);
360 :
361 4160 : ADDenseVector local_vel_dv_0(3);
362 4160 : ADDenseVector local_vel_dv_1(3);
363 4160 : ADDenseVector local_vel_dv_2(3);
364 4160 : ADDenseVector local_vel_dv_3(3);
365 :
366 4160 : ADDenseVector local_rot_vel_dv_0(3);
367 4160 : ADDenseVector local_rot_vel_dv_1(3);
368 4160 : ADDenseVector local_rot_vel_dv_2(3);
369 4160 : ADDenseVector local_rot_vel_dv_3(3);
370 :
371 4160 : ADDenseVector local_old_vel_dv_0(3);
372 4160 : ADDenseVector local_old_vel_dv_1(3);
373 4160 : ADDenseVector local_old_vel_dv_2(3);
374 4160 : ADDenseVector local_old_vel_dv_3(3);
375 :
376 4160 : ADDenseVector local_old_rot_vel_dv_0(3);
377 4160 : ADDenseVector local_old_rot_vel_dv_1(3);
378 4160 : ADDenseVector local_old_rot_vel_dv_2(3);
379 4160 : ADDenseVector local_old_rot_vel_dv_3(3);
380 :
381 16640 : for (unsigned int i = 0; i < 3; i++)
382 : {
383 12480 : local_accel_dv_0(i) = _local_accel.pos[0](i);
384 12480 : local_accel_dv_1(i) = _local_accel.pos[1](i);
385 12480 : local_accel_dv_2(i) = _local_accel.pos[2](i);
386 12480 : local_accel_dv_3(i) = _local_accel.pos[3](i);
387 :
388 12480 : local_rot_accel_dv_0(i) = _local_accel.rot[0](i);
389 12480 : local_rot_accel_dv_1(i) = _local_accel.rot[1](i);
390 12480 : local_rot_accel_dv_2(i) = _local_accel.rot[2](i);
391 12480 : local_rot_accel_dv_3(i) = _local_accel.rot[3](i);
392 :
393 12480 : local_vel_dv_0(i) = _local_vel.pos[0](i);
394 12480 : local_vel_dv_1(i) = _local_vel.pos[1](i);
395 12480 : local_vel_dv_2(i) = _local_vel.pos[2](i);
396 12480 : local_vel_dv_3(i) = _local_vel.pos[3](i);
397 :
398 12480 : local_rot_vel_dv_0(i) = _local_vel.rot[0](i);
399 12480 : local_rot_vel_dv_1(i) = _local_vel.rot[1](i);
400 12480 : local_rot_vel_dv_2(i) = _local_vel.rot[2](i);
401 12480 : local_rot_vel_dv_3(i) = _local_vel.rot[3](i);
402 :
403 12480 : local_old_vel_dv_0(i) = _local_old_vel.pos[0](i);
404 12480 : local_old_vel_dv_1(i) = _local_old_vel.pos[1](i);
405 12480 : local_old_vel_dv_2(i) = _local_old_vel.pos[2](i);
406 12480 : local_old_vel_dv_3(i) = _local_old_vel.pos[3](i);
407 :
408 12480 : local_old_rot_vel_dv_0(i) = _local_old_vel.rot[0](i);
409 12480 : local_old_rot_vel_dv_1(i) = _local_old_vel.rot[1](i);
410 12480 : local_old_rot_vel_dv_2(i) = _local_old_vel.rot[2](i);
411 12480 : local_old_rot_vel_dv_3(i) = _local_old_vel.rot[3](i);
412 : }
413 4160 : unsigned int dim = _current_elem->dim();
414 :
415 : // Update 0g vectors at plane quadrature points.
416 8320 : _0g1_vectors[0] = -0.5 * _thickness * _v1[0];
417 4160 : _0g1_vectors[1] = 0.5 * _thickness * _v2[0];
418 :
419 4160 : ADDenseMatrix G1(3, 2);
420 4160 : G1(0, 0) = _0g1_vectors[0](0);
421 4160 : G1(1, 0) = _0g1_vectors[0](1);
422 4160 : G1(2, 0) = _0g1_vectors[0](2);
423 4160 : G1(0, 1) = _0g1_vectors[1](0);
424 4160 : G1(1, 1) = _0g1_vectors[1](1);
425 4160 : G1(2, 1) = _0g1_vectors[1](2);
426 4160 : ADDenseMatrix G1T(2, 3);
427 4160 : G1.get_transpose(G1T);
428 :
429 8320 : _0g2_vectors[0] = -0.5 * _thickness * _v1[1];
430 4160 : _0g2_vectors[1] = 0.5 * _thickness * _v2[1];
431 :
432 4160 : ADDenseMatrix G2(3, 2);
433 4160 : G2(0, 0) = _0g2_vectors[0](0);
434 4160 : G2(1, 0) = _0g2_vectors[0](1);
435 4160 : G2(2, 0) = _0g2_vectors[0](2);
436 4160 : G2(0, 1) = _0g2_vectors[1](0);
437 4160 : G2(1, 1) = _0g2_vectors[1](1);
438 4160 : G2(2, 1) = _0g2_vectors[1](2);
439 :
440 4160 : ADDenseMatrix G2T(2, 3);
441 4160 : G2.get_transpose(G2T);
442 :
443 8320 : _0g3_vectors[0] = -0.5 * _thickness * _v1[2];
444 4160 : _0g3_vectors[1] = 0.5 * _thickness * _v2[2];
445 :
446 4160 : ADDenseMatrix G3(3, 2);
447 4160 : G3(0, 0) = _0g3_vectors[0](0);
448 4160 : G3(1, 0) = _0g3_vectors[0](1);
449 4160 : G3(2, 0) = _0g3_vectors[0](2);
450 4160 : G3(0, 1) = _0g3_vectors[1](0);
451 4160 : G3(1, 1) = _0g3_vectors[1](1);
452 4160 : G3(2, 1) = _0g3_vectors[1](2);
453 :
454 4160 : ADDenseMatrix G3T(2, 3);
455 4160 : G3.get_transpose(G3T);
456 :
457 8320 : _0g4_vectors[0] = -0.5 * _thickness * _v1[3];
458 4160 : _0g4_vectors[1] = 0.5 * _thickness * _v2[3];
459 :
460 4160 : ADDenseMatrix G4(3, 2);
461 4160 : G4(0, 0) = _0g4_vectors[0](0);
462 4160 : G4(1, 0) = _0g4_vectors[0](1);
463 4160 : G4(2, 0) = _0g4_vectors[0](2);
464 4160 : G4(0, 1) = _0g4_vectors[1](0);
465 4160 : G4(1, 1) = _0g4_vectors[1](1);
466 4160 : G4(2, 1) = _0g4_vectors[1](2);
467 :
468 4160 : ADDenseMatrix G4T(2, 3);
469 4160 : G4.get_transpose(G4T);
470 :
471 : std::vector<ADDenseVector> local_acc;
472 4160 : 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 4160 : 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 4160 : 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 4160 : local_vel_dv_0.scale(1 + _alpha);
504 4160 : local_old_vel_dv_0.scale(_alpha);
505 4160 : local_vel_dv_1.scale(1 + _alpha);
506 4160 : local_old_vel_dv_1.scale(_alpha);
507 4160 : local_vel_dv_2.scale(1 + _alpha);
508 4160 : local_old_vel_dv_2.scale(_alpha);
509 4160 : local_vel_dv_3.scale(1 + _alpha);
510 4160 : local_old_vel_dv_3.scale(_alpha);
511 :
512 4160 : local_vel_dv_0.add(1.0, local_old_vel_dv_0);
513 4160 : local_vel_dv_1.add(1.0, local_old_vel_dv_1);
514 4160 : local_vel_dv_2.add(1.0, local_old_vel_dv_2);
515 4160 : 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 4160 : 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 4160 : local_rot_vel_dv_0.scale(1 + _alpha);
530 4160 : local_old_rot_vel_dv_0.scale(_alpha);
531 4160 : local_rot_vel_dv_1.scale(1 + _alpha);
532 4160 : local_old_rot_vel_dv_1.scale(_alpha);
533 4160 : local_rot_vel_dv_2.scale(1 + _alpha);
534 4160 : local_old_rot_vel_dv_2.scale(_alpha);
535 4160 : local_rot_vel_dv_3.scale(1 + _alpha);
536 4160 : local_old_rot_vel_dv_3.scale(_alpha);
537 :
538 4160 : local_rot_vel_dv_0.add(1.0, local_old_rot_vel_dv_0);
539 4160 : local_rot_vel_dv_1.add(1.0, local_old_rot_vel_dv_1);
540 4160 : local_rot_vel_dv_2.add(1.0, local_old_rot_vel_dv_2);
541 4160 : 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 4160 : FEType fe_type(Utility::string_to_enum<Order>("First"),
550 4160 : Utility::string_to_enum<FEFamily>("LAGRANGE"));
551 4160 : auto & fe = _fe_problem.assembly(_tid, _sys.number()).getFE(fe_type, dim);
552 4160 : _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
553 4160 : _dphideta_map = fe->get_fe_map().get_dphideta_map();
554 4160 : _phi_map = fe->get_fe_map().get_phi_map();
555 :
556 : // quadrature points in isoparametric space
557 4160 : _2d_points = _qrule->get_points(); // would be in 2D
558 :
559 : std::vector<const Node *> nodes;
560 20800 : for (unsigned int i = 0; i < 4; ++i)
561 16640 : nodes.push_back(_current_elem->node_ptr(i));
562 :
563 16640 : for (unsigned int i = 0; i < _ndisp; i++)
564 62400 : for (unsigned int j = 0; j < 4; j++)
565 49920 : _local_force[j](i) = 0.0;
566 :
567 12480 : for (unsigned int i = 0; i < _nrot; i++)
568 41600 : for (unsigned int j = 0; j < 4; j++)
569 33280 : _local_moment[j](i) = 0.0;
570 :
571 20800 : for (unsigned int qp_xy = 0; qp_xy < _2d_points.size(); ++qp_xy)
572 : {
573 16640 : 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 66560 : for (unsigned int dim = 0; dim < 3; dim++)
577 : {
578 : _local_force[0](dim) +=
579 49920 : factor_qxy * (_phi_map[0][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
580 49920 : _phi_map[0][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
581 49920 : _phi_map[0][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
582 99840 : _phi_map[0][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
583 :
584 49920 : if (_eta[0] > TOLERANCE * TOLERANCE)
585 0 : _local_force[0](dim) += factor_qxy * _eta[0] *
586 0 : (_phi_map[0][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
587 0 : _phi_map[0][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
588 0 : _phi_map[0][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
589 0 : _phi_map[0][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
590 :
591 : _local_force[1](dim) +=
592 49920 : factor_qxy * (_phi_map[1][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
593 49920 : _phi_map[1][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
594 49920 : _phi_map[1][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
595 99840 : _phi_map[1][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
596 :
597 49920 : if (_eta[0] > TOLERANCE * TOLERANCE)
598 0 : _local_force[1](dim) += factor_qxy * _eta[0] *
599 0 : (_phi_map[1][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
600 0 : _phi_map[1][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
601 0 : _phi_map[1][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
602 0 : _phi_map[1][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
603 : _local_force[2](dim) +=
604 49920 : factor_qxy * (_phi_map[2][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
605 49920 : _phi_map[2][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
606 49920 : _phi_map[2][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
607 99840 : _phi_map[2][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
608 :
609 49920 : if (_eta[0] > TOLERANCE * TOLERANCE)
610 0 : _local_force[2](dim) += factor_qxy * _eta[0] *
611 0 : (_phi_map[2][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
612 0 : _phi_map[2][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
613 0 : _phi_map[2][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
614 0 : _phi_map[2][qp_xy] * _phi_map[3][qp_xy] * local_vel[3](dim));
615 :
616 : _local_force[3](dim) +=
617 49920 : factor_qxy * (_phi_map[3][qp_xy] * _phi_map[0][qp_xy] * local_acc[0](dim) +
618 49920 : _phi_map[3][qp_xy] * _phi_map[1][qp_xy] * local_acc[1](dim) +
619 49920 : _phi_map[3][qp_xy] * _phi_map[2][qp_xy] * local_acc[2](dim) +
620 99840 : _phi_map[3][qp_xy] * _phi_map[3][qp_xy] * local_acc[3](dim));
621 :
622 49920 : if (_eta[0] > TOLERANCE * TOLERANCE)
623 0 : _local_force[3](dim) += factor_qxy * _eta[0] *
624 0 : (_phi_map[3][qp_xy] * _phi_map[0][qp_xy] * local_vel[0](dim) +
625 0 : _phi_map[3][qp_xy] * _phi_map[1][qp_xy] * local_vel[1](dim) +
626 0 : _phi_map[3][qp_xy] * _phi_map[2][qp_xy] * local_vel[2](dim) +
627 0 : _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 16640 : ADReal rot_thickness = _thickness * _thickness * _thickness / 48.0;
632 :
633 16640 : ADDenseVector momentInertia(3);
634 16640 : momentInertia(0) = (G1(0, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
635 16640 : G1(0, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
636 16640 : G2(0, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
637 16640 : G2(0, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
638 16640 : G3(0, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
639 16640 : G3(0, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
640 16640 : G4(0, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
641 33280 : G4(0, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
642 :
643 16640 : momentInertia(1) = (G1(1, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
644 16640 : G1(1, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
645 16640 : G2(1, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
646 16640 : G2(1, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
647 16640 : G3(1, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
648 16640 : G3(1, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
649 16640 : G4(1, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
650 33280 : G4(1, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
651 :
652 16640 : momentInertia(2) = (G1(2, 0) * (local_rot_acc[0](0) + _eta[0] * local_rot_vel[0](0)) +
653 16640 : G1(2, 1) * (local_rot_acc[0](1) + _eta[0] * local_rot_vel[0](1)) +
654 16640 : G2(2, 0) * (local_rot_acc[1](0) + _eta[0] * local_rot_vel[1](0)) +
655 16640 : G2(2, 1) * (local_rot_acc[1](1) + _eta[0] * local_rot_vel[1](1)) +
656 16640 : G3(2, 0) * (local_rot_acc[2](0) + _eta[0] * local_rot_vel[2](0)) +
657 16640 : G3(2, 1) * (local_rot_acc[2](1) + _eta[0] * local_rot_vel[2](1)) +
658 16640 : G4(2, 0) * (local_rot_acc[3](0) + _eta[0] * local_rot_vel[3](0)) +
659 33280 : G4(2, 1) * (local_rot_acc[3](1) + _eta[0] * local_rot_vel[3](1)));
660 :
661 16640 : _local_moment[0](0) += factor_qxy * rot_thickness *
662 16640 : (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
663 16640 : G1T(0, 2) * momentInertia(2));
664 :
665 16640 : _local_moment[0](1) += factor_qxy * rot_thickness *
666 16640 : (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
667 16640 : G1T(1, 2) * momentInertia(2));
668 :
669 16640 : _local_moment[1](0) += factor_qxy * rot_thickness *
670 16640 : (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
671 16640 : G2T(0, 2) * momentInertia(2));
672 :
673 16640 : _local_moment[1](1) += factor_qxy * rot_thickness *
674 16640 : (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
675 16640 : G2T(1, 2) * momentInertia(2));
676 :
677 16640 : _local_moment[2](0) += factor_qxy * rot_thickness *
678 16640 : (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
679 16640 : G3T(0, 2) * momentInertia(2));
680 :
681 16640 : _local_moment[2](1) += factor_qxy * rot_thickness *
682 16640 : (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
683 16640 : G3T(1, 2) * momentInertia(2));
684 :
685 16640 : _local_moment[3](0) += factor_qxy * rot_thickness *
686 16640 : (G1T(0, 0) * momentInertia(0) + G1T(0, 1) * momentInertia(1) +
687 16640 : G4T(0, 2) * momentInertia(2));
688 :
689 16640 : _local_moment[3](1) += factor_qxy * rot_thickness *
690 16640 : (G1T(1, 0) * momentInertia(0) + G1T(1, 1) * momentInertia(1) +
691 16640 : G4T(1, 2) * momentInertia(2));
692 : }
693 8320 : }
|
