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