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