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