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 "ADComputeIncrementalShellStrain.h"
11 : #include "MooseMesh.h"
12 : #include "Assembly.h"
13 : #include "NonlinearSystem.h"
14 : #include "MooseVariable.h"
15 : #include "ArbitraryQuadrature.h"
16 : #include "DenseMatrix.h"
17 :
18 : #include "libmesh/quadrature.h"
19 : #include "libmesh/utility.h"
20 : #include "libmesh/enum_quadrature_type.h"
21 : #include "libmesh/fe_type.h"
22 : #include "libmesh/string_to_enum.h"
23 : #include "libmesh/quadrature_gauss.h"
24 :
25 : registerMooseObject("TensorMechanicsApp", ADComputeIncrementalShellStrain);
26 :
27 : InputParameters
28 160 : ADComputeIncrementalShellStrain::validParams()
29 : {
30 160 : InputParameters params = Material::validParams();
31 160 : params.addClassDescription("Compute a small strain increment for the shell.");
32 320 : params.addRequiredCoupledVar(
33 : "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
34 320 : params.addRequiredCoupledVar(
35 : "displacements",
36 : "The displacements appropriate for the simulation geometry and coordinate system");
37 320 : params.addRequiredCoupledVar(
38 : "thickness",
39 : "Thickness of the shell. Can be supplied as either a number or a variable name.");
40 320 : params.addRequiredParam<std::string>("through_thickness_order",
41 : "Quadrature order in out of plane direction");
42 320 : params.addParam<bool>(
43 320 : "large_strain", false, "Set to true to turn on finite strain calculations.");
44 160 : return params;
45 0 : }
46 :
47 120 : ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputParameters & parameters)
48 : : Material(parameters),
49 120 : _nrot(coupledComponents("rotations")),
50 120 : _ndisp(coupledComponents("displacements")),
51 120 : _rot_num(_nrot),
52 120 : _disp_num(_ndisp),
53 120 : _thickness(coupledValue("thickness")),
54 240 : _large_strain(getParam<bool>("large_strain")),
55 120 : _strain_increment(),
56 120 : _total_strain(),
57 120 : _total_strain_old(),
58 120 : _nonlinear_sys(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0)),
59 120 : _soln_disp_index(4),
60 120 : _soln_rot_index(4),
61 120 : _soln_vector(20, 1),
62 120 : _strain_vector(5, 1),
63 120 : _nodes(4),
64 120 : _node_normal(declareADProperty<RealVectorValue>("node_normal")),
65 240 : _node_normal_old(getMaterialPropertyOldByName<RealVectorValue>("node_normal")),
66 120 : _dxyz_dxi(),
67 120 : _dxyz_deta(),
68 120 : _dxyz_dzeta(),
69 120 : _dxyz_dxi_old(),
70 120 : _dxyz_deta_old(),
71 120 : _dxyz_dzeta_old(),
72 120 : _v1(4),
73 120 : _v2(4),
74 120 : _B(),
75 120 : _B_old(),
76 120 : _ge(),
77 120 : _ge_old(),
78 120 : _J_map(),
79 120 : _J_map_old(),
80 120 : _covariant_transformation_matrix(),
81 120 : _covariant_transformation_matrix_old(),
82 120 : _contravariant_transformation_matrix(),
83 120 : _contravariant_transformation_matrix_old(),
84 120 : _total_global_strain(),
85 120 : _sol(_nonlinear_sys.currentSolution()),
86 480 : _sol_old(_nonlinear_sys.solutionOld())
87 : {
88 : // Checking for consistency between length of the provided displacements and rotations vector
89 120 : if (_ndisp != 3 || _nrot != 2)
90 0 : mooseError(
91 : "ADComputeIncrementalShellStrain: The number of variables supplied in 'displacements' "
92 : "must be 3 and that in 'rotations' must be 2.");
93 :
94 120 : if (_mesh.hasSecondOrderElements())
95 0 : mooseError(
96 : "ADComputeIncrementalShellStrain: Shell element is implemented only for linear elements.");
97 :
98 : // fetch coupled variables and gradients (as stateful properties if necessary)
99 480 : for (unsigned int i = 0; i < _ndisp; ++i)
100 : {
101 720 : MooseVariable * disp_variable = getVar("displacements", i);
102 360 : _disp_num[i] = disp_variable->number();
103 :
104 360 : if (i < _nrot)
105 : {
106 480 : MooseVariable * rot_variable = getVar("rotations", i);
107 240 : _rot_num[i] = rot_variable->number();
108 : }
109 : }
110 :
111 120 : _t_qrule = std::make_unique<QGauss>(
112 360 : 1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
113 120 : _t_points = _t_qrule->get_points();
114 120 : _strain_increment.resize(_t_points.size());
115 120 : _total_strain.resize(_t_points.size());
116 120 : _total_strain_old.resize(_t_points.size());
117 120 : _B.resize(_t_points.size());
118 120 : _B_old.resize(_t_points.size());
119 120 : _ge.resize(_t_points.size());
120 120 : _ge_old.resize(_t_points.size());
121 120 : _J_map.resize(_t_points.size());
122 120 : _J_map_old.resize(_t_points.size());
123 120 : _dxyz_dxi.resize(_t_points.size());
124 120 : _dxyz_deta.resize(_t_points.size());
125 120 : _dxyz_dzeta.resize(_t_points.size());
126 120 : _dxyz_dxi_old.resize(_t_points.size());
127 120 : _dxyz_deta_old.resize(_t_points.size());
128 120 : _dxyz_dzeta_old.resize(_t_points.size());
129 120 : _covariant_transformation_matrix.resize(_t_points.size());
130 120 : _covariant_transformation_matrix_old.resize(_t_points.size());
131 120 : _contravariant_transformation_matrix.resize(_t_points.size());
132 120 : _contravariant_transformation_matrix_old.resize(_t_points.size());
133 120 : _total_global_strain.resize(_t_points.size());
134 :
135 120 : _transformation_matrix = &declareADProperty<RankTwoTensor>("transformation_matrix_element");
136 :
137 360 : for (unsigned int i = 0; i < _t_points.size(); ++i)
138 : {
139 240 : _strain_increment[i] =
140 240 : &declareADProperty<RankTwoTensor>("strain_increment_t_points_" + std::to_string(i));
141 240 : _total_strain[i] =
142 240 : &declareADProperty<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
143 240 : _total_strain_old[i] =
144 240 : &getMaterialPropertyOldByName<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
145 240 : _B[i] = &declareADProperty<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
146 480 : _B_old[i] = &getMaterialPropertyOldByName<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
147 240 : _ge[i] = &declareADProperty<RankTwoTensor>("ge_t_points_" + std::to_string(i));
148 480 : _ge_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>("ge_t_points_" + std::to_string(i));
149 240 : _J_map[i] = &declareADProperty<Real>("J_mapping_t_points_" + std::to_string(i));
150 480 : _J_map_old[i] = &getMaterialPropertyOldByName<Real>("J_mapping_t_points_" + std::to_string(i));
151 240 : _dxyz_dxi[i] = &declareADProperty<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
152 240 : _dxyz_dxi_old[i] =
153 240 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
154 240 : _dxyz_deta[i] = &declareADProperty<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
155 240 : _dxyz_deta_old[i] =
156 240 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
157 240 : _dxyz_dzeta[i] =
158 240 : &declareADProperty<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
159 240 : _dxyz_dzeta_old[i] =
160 240 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
161 : // Create rotation matrix and total strain global for output purposes only
162 240 : _covariant_transformation_matrix[i] =
163 240 : &declareProperty<RankTwoTensor>("covariant_transformation_t_points_" + std::to_string(i));
164 240 : _covariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
165 240 : "covariant_transformation_t_points_" + std::to_string(i));
166 240 : _contravariant_transformation_matrix[i] = &declareProperty<RankTwoTensor>(
167 240 : "contravariant_transformation_t_points_" + std::to_string(i));
168 240 : _contravariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
169 240 : "contravariant_transformation_t_points_" + std::to_string(i));
170 240 : _total_global_strain[i] =
171 480 : &declareProperty<RankTwoTensor>("total_global_strain_t_points_" + std::to_string(i));
172 : }
173 :
174 : // used later for computing local coordinate system
175 120 : _x2(1) = 1;
176 120 : _x3(2) = 1;
177 120 : }
178 :
179 : void
180 2037 : ADComputeIncrementalShellStrain::initQpStatefulProperties()
181 : {
182 2037 : unsigned int dim = _current_elem->dim();
183 2037 : if ((dim != 2))
184 0 : mooseError(
185 : "ADComputeIncrementalShellStrain: Shell element is implemented only for 2D elements");
186 2037 : if (_current_elem->n_nodes() != 4)
187 0 : mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four nodes.");
188 2037 : if (_qrule->get_points().size() != 4)
189 1 : mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four "
190 : "quadrature points.");
191 2036 : computeGMatrix();
192 2036 : computeBMatrix();
193 2036 : }
194 :
195 : void
196 3118 : ADComputeIncrementalShellStrain::computeProperties()
197 : {
198 : // quadrature points in isoparametric space
199 3118 : _2d_points = _qrule->get_points(); // would be in 2D
200 :
201 : // derivatives of shape functions (dphidxi, dphideta and dphidzeta) evaluated at quadrature points
202 : // (in isoparametric space).
203 3118 : unsigned int dim = _current_elem->dim();
204 3118 : FEType fe_type(Utility::string_to_enum<Order>("First"),
205 3118 : Utility::string_to_enum<FEFamily>("LAGRANGE"));
206 3118 : auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
207 3118 : _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
208 3118 : _dphideta_map = fe->get_fe_map().get_dphideta_map();
209 3118 : _phi_map = fe->get_fe_map().get_phi_map();
210 :
211 15590 : for (unsigned int i = 0; i < 4; ++i)
212 12472 : _nodes[i] = _current_elem->node_ptr(i);
213 :
214 15590 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
215 : {
216 37416 : for (unsigned int j = 0; j < _t_points.size(); ++j)
217 : {
218 24944 : (*_ge[j])[i] = (*_ge_old[j])[i];
219 24944 : (*_J_map[j])[i] = (*_J_map_old[j])[i];
220 24944 : (*_dxyz_dxi[j])[i] = (*_dxyz_dxi_old[j])[i];
221 24944 : (*_dxyz_deta[j])[i] = (*_dxyz_deta_old[j])[i];
222 24944 : (*_dxyz_dzeta[j])[i] = (*_dxyz_dzeta_old[j])[i];
223 24944 : (*_B[j])[i] = (*_B_old[j])[i];
224 24944 : (*_covariant_transformation_matrix[j])[i] = (*_covariant_transformation_matrix_old[j])[i];
225 24944 : (*_contravariant_transformation_matrix[j])[i] =
226 24944 : (*_contravariant_transformation_matrix_old[j])[i];
227 : }
228 : }
229 :
230 3118 : computeSolnVector();
231 :
232 3118 : computeNodeNormal();
233 :
234 15590 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
235 : {
236 37416 : for (unsigned int j = 0; j < _t_points.size(); ++j)
237 : {
238 : // compute strain increment in covariant coordinate system using B and _soln_vector
239 149664 : for (unsigned int temp1 = 0; temp1 < 5; ++temp1)
240 : {
241 124720 : _strain_vector(temp1) = 0.0;
242 2619120 : for (unsigned int temp2 = 0; temp2 < 20; ++temp2)
243 4988800 : _strain_vector(temp1) += (*_B[j])[i](temp1, temp2) * _soln_vector(temp2);
244 : }
245 24944 : (*_strain_increment[j])[i](0, 0) = _strain_vector(0);
246 24944 : (*_strain_increment[j])[i](1, 1) = _strain_vector(1);
247 24944 : (*_strain_increment[j])[i](0, 1) = _strain_vector(2);
248 24944 : (*_strain_increment[j])[i](0, 2) = _strain_vector(3);
249 24944 : (*_strain_increment[j])[i](1, 2) = _strain_vector(4);
250 24944 : (*_strain_increment[j])[i](1, 0) = (*_strain_increment[j])[i](0, 1);
251 24944 : (*_strain_increment[j])[i](2, 0) = (*_strain_increment[j])[i](0, 2);
252 24944 : (*_strain_increment[j])[i](2, 1) = (*_strain_increment[j])[i](1, 2);
253 :
254 24944 : (*_total_strain[j])[i] = (*_total_strain_old[j])[i] + (*_strain_increment[j])[i];
255 :
256 99776 : for (unsigned int ii = 0; ii < 3; ++ii)
257 299328 : for (unsigned int jj = 0; jj < 3; ++jj)
258 224496 : _unrotated_total_strain(ii, jj) = MetaPhysicL::raw_value((*_total_strain[j])[i](ii, jj));
259 24944 : (*_total_global_strain[j])[i] = (*_contravariant_transformation_matrix[j])[i] *
260 49888 : _unrotated_total_strain *
261 24944 : (*_contravariant_transformation_matrix[j])[i].transpose();
262 : }
263 : }
264 3118 : }
265 :
266 : void
267 2212 : ADComputeIncrementalShellStrain::computeGMatrix()
268 : {
269 : // quadrature points in isoparametric space
270 2212 : _2d_points = _qrule->get_points(); // would be in 2D
271 :
272 2212 : unsigned int dim = _current_elem->dim();
273 :
274 : // derivatives of shape functions (dphidxi, dphideta and dphidzeta) evaluated at quadrature points
275 : // (in isoparametric space).
276 2212 : FEType fe_type(Utility::string_to_enum<Order>("First"),
277 2212 : Utility::string_to_enum<FEFamily>("LAGRANGE"));
278 2212 : auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
279 2212 : _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
280 2212 : _dphideta_map = fe->get_fe_map().get_dphideta_map();
281 2212 : _phi_map = fe->get_fe_map().get_phi_map();
282 :
283 11060 : for (unsigned int i = 0; i < 4; ++i)
284 8848 : _nodes[i] = _current_elem->node_ptr(i);
285 :
286 2212 : ADRealVectorValue x = (*_nodes[1] - *_nodes[0]);
287 2212 : ADRealVectorValue y = (*_nodes[3] - *_nodes[0]);
288 2212 : ADRealVectorValue normal = x.cross(y);
289 2212 : normal /= normal.norm();
290 :
291 11060 : for (unsigned int k = 0; k < 4; ++k)
292 8848 : _node_normal[k] = normal;
293 :
294 2212 : ADRankTwoTensor a;
295 2212 : ADDenseMatrix b(5, 20);
296 : ADRealVectorValue c;
297 2212 : RankTwoTensor d;
298 6636 : for (unsigned int t = 0; t < _t_points.size(); ++t)
299 : {
300 4424 : (*_strain_increment[t])[_qp] = a;
301 4424 : (*_total_strain[t])[_qp] = a;
302 4424 : (*_B[t])[_qp] = b;
303 4424 : (*_ge[t])[_qp] = a;
304 4424 : (*_J_map[t])[_qp] = 0;
305 4424 : (*_dxyz_dxi[t])[_qp] = c;
306 4424 : (*_dxyz_deta[t])[_qp] = c;
307 4424 : (*_dxyz_dzeta[t])[_qp] = c;
308 4424 : (*_covariant_transformation_matrix[t])[_qp] = d;
309 4424 : (*_contravariant_transformation_matrix[t])[_qp] = d;
310 : }
311 :
312 : // calculating derivatives of shape function in physical space (dphi/dx, dphi/dy, dphi/dz) at
313 : // quadrature points these are g_{i} in Dvorkin's paper
314 11060 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
315 : {
316 26544 : for (unsigned int j = 0; j < _t_points.size(); ++j)
317 : {
318 17696 : (*_dxyz_dxi[j])[i].zero();
319 70784 : for (unsigned int component = 0; component < 3; ++component)
320 : {
321 53088 : (*_dxyz_dxi[j])[i](component) = 0.0;
322 53088 : (*_dxyz_deta[j])[i](component) = 0.0;
323 53088 : (*_dxyz_dzeta[j])[i](component) = 0.0;
324 :
325 265440 : for (unsigned int k = 0; k < _nodes.size(); ++k)
326 : {
327 424704 : (*_dxyz_dxi[j])[i](component) += _dphidxi_map[k][i] * ((*_nodes[k])(component)) +
328 212352 : _t_points[j](0) / 2.0 * _thickness[i] *
329 424704 : _dphidxi_map[k][i] * _node_normal[k](component);
330 424704 : (*_dxyz_deta[j])[i](component) += _dphideta_map[k][i] * ((*_nodes[k])(component)) +
331 212352 : _t_points[j](0) / 2.0 * _thickness[i] *
332 424704 : _dphideta_map[k][i] * _node_normal[k](component);
333 212352 : (*_dxyz_dzeta[j])[i](component) +=
334 424704 : _thickness[i] * _phi_map[k][i] * _node_normal[k](component) / 2.0;
335 : }
336 : }
337 : }
338 : }
339 :
340 11060 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
341 : {
342 26544 : for (unsigned int j = 0; j < _t_points.size(); ++j)
343 : {
344 : // calculate gij for elasticity tensor
345 17696 : ADRankTwoTensor gmn;
346 17696 : RankTwoTensor J;
347 70784 : for (unsigned int component = 0; component < 3; ++component)
348 : {
349 106176 : gmn(0, 0) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dxi[j])[i](component);
350 106176 : gmn(1, 1) += (*_dxyz_deta[j])[i](component) * (*_dxyz_deta[j])[i](component);
351 106176 : gmn(2, 2) += (*_dxyz_dzeta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
352 106176 : gmn(0, 1) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_deta[j])[i](component);
353 106176 : gmn(0, 2) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
354 106176 : gmn(1, 2) += (*_dxyz_deta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
355 :
356 53088 : J(0, component) = MetaPhysicL::raw_value((*_dxyz_dxi[j])[i](component));
357 53088 : J(1, component) = MetaPhysicL::raw_value((*_dxyz_deta[j])[i](component));
358 53088 : J(2, component) = MetaPhysicL::raw_value((*_dxyz_dzeta[j])[i](component));
359 : }
360 17696 : gmn(1, 0) = gmn(0, 1);
361 17696 : gmn(2, 0) = gmn(0, 2);
362 17696 : gmn(2, 1) = gmn(1, 2);
363 :
364 17696 : ADRankTwoTensor gmninv_temp = gmn.inverse();
365 17696 : (*_J_map[j])[i] = std::sqrt(gmn.det());
366 17696 : (*_covariant_transformation_matrix[j])[i] = J;
367 :
368 17696 : (*_contravariant_transformation_matrix[j])[i] =
369 17696 : (*_covariant_transformation_matrix[j])[i].inverse();
370 :
371 17696 : Real normx = std::sqrt(J(0, 0) * J(0, 0) + J(0, 1) * J(0, 1) + J(0, 2) * J(0, 2));
372 17696 : Real normy = std::sqrt(J(1, 0) * J(1, 0) + J(1, 1) * J(1, 1) + J(1, 2) * J(1, 2));
373 17696 : Real normz = std::sqrt(J(2, 0) * J(2, 0) + J(2, 1) * J(2, 1) + J(2, 2) * J(2, 2));
374 :
375 17696 : J(0, 0) /= normx;
376 17696 : J(0, 1) /= normx;
377 17696 : J(0, 2) /= normx;
378 :
379 17696 : J(1, 0) /= normy;
380 17696 : J(1, 1) /= normy;
381 17696 : J(1, 2) /= normy;
382 :
383 17696 : J(2, 0) /= normz;
384 17696 : J(2, 1) /= normz;
385 17696 : J(2, 2) /= normz;
386 :
387 17696 : (*_transformation_matrix)[i] = J;
388 :
389 : // calculate ge
390 35392 : ADRealVectorValue e3 = (*_dxyz_dzeta[j])[i] / (*_dxyz_dzeta[j])[i].norm();
391 17696 : ADRealVectorValue e1 = (*_dxyz_deta[j])[i].cross(e3);
392 17696 : e1 /= e1.norm();
393 17696 : ADRealVectorValue e2 = e3.cross(e1);
394 17696 : e2 /= e2.norm();
395 :
396 17696 : ADRankTwoTensor local_rotation_mat;
397 17696 : local_rotation_mat(0, 0) = e1(0);
398 17696 : local_rotation_mat(0, 1) = e1(1);
399 17696 : local_rotation_mat(0, 2) = e1(2);
400 17696 : local_rotation_mat(1, 0) = e2(0);
401 17696 : local_rotation_mat(1, 1) = e2(1);
402 17696 : local_rotation_mat(1, 2) = e2(2);
403 17696 : local_rotation_mat(2, 0) = e3(0);
404 17696 : local_rotation_mat(2, 1) = e3(1);
405 17696 : local_rotation_mat(2, 2) = e3(2);
406 :
407 17696 : ADRankTwoTensor gmninv = local_rotation_mat.transpose() * gmninv_temp * local_rotation_mat;
408 :
409 17696 : (*_ge[j])[i](0, 0) = (gmninv * (*_dxyz_dxi[j])[i]) * e1;
410 17696 : (*_ge[j])[i](0, 1) = (gmninv * (*_dxyz_dxi[j])[i]) * e2;
411 17696 : (*_ge[j])[i](0, 2) = (gmninv * (*_dxyz_dxi[j])[i]) * e3;
412 17696 : (*_ge[j])[i](1, 0) = (gmninv * (*_dxyz_deta[j])[i]) * e1;
413 17696 : (*_ge[j])[i](1, 1) = (gmninv * (*_dxyz_deta[j])[i]) * e2;
414 17696 : (*_ge[j])[i](1, 2) = (gmninv * (*_dxyz_deta[j])[i]) * e3;
415 17696 : (*_ge[j])[i](2, 0) = (gmninv * (*_dxyz_dzeta[j])[i]) * e1;
416 17696 : (*_ge[j])[i](2, 1) = (gmninv * (*_dxyz_dzeta[j])[i]) * e2;
417 17696 : (*_ge[j])[i](2, 2) = (gmninv * (*_dxyz_dzeta[j])[i]) * e3;
418 : }
419 : }
420 2212 : }
421 :
422 : void
423 3118 : ADComputeIncrementalShellStrain::computeNodeNormal()
424 : {
425 15590 : for (unsigned int k = 0; k < _nodes.size(); ++k)
426 12472 : _node_normal[k] = _node_normal_old[k];
427 3118 : }
428 :
429 : void
430 3888 : ADComputeIncrementalShellStrain::computeBMatrix()
431 : {
432 : // compute nodal local axis
433 19440 : for (unsigned int k = 0; k < _nodes.size(); ++k)
434 : {
435 15552 : _v1[k] = _x2.cross(_node_normal[k]);
436 31104 : _v1[k] /= _x2.norm() * _node_normal[k].norm();
437 :
438 : // If x2 is parallel to node normal, set V1 to x3
439 15552 : if (MooseUtils::absoluteFuzzyEqual(_v1[k].norm(), 0.0, 1e-6))
440 : _v1[k] = _x3;
441 :
442 31104 : _v2[k] = _node_normal[k].cross(_v1[k]);
443 : }
444 :
445 : // compute B matrix rows correspond to [ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, uz1, uz2, uz3,
446 : // uz4, a1, a2, a3, a4, b1, b2, b3, b4]
447 19440 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
448 : {
449 46656 : for (unsigned int j = 0; j < _t_points.size(); ++j)
450 : {
451 31104 : (*_B[j])[i].resize(5, 20);
452 31104 : (*_B[j])[i].zero();
453 155520 : for (unsigned int k = 0; k < _nodes.size(); ++k)
454 : {
455 : // corresponding to strain(0,0)
456 248832 : (*_B[j])[i](0, k) += _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](0);
457 248832 : (*_B[j])[i](0, 4 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](1);
458 248832 : (*_B[j])[i](0, 8 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](2);
459 248832 : (*_B[j])[i](0, 12 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
460 248832 : (-_v2[k] * (*_dxyz_dxi[j])[i]);
461 248832 : (*_B[j])[i](0, 16 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
462 248832 : (_v1[k] * (*_dxyz_dxi[j])[i]);
463 :
464 : // corresponding to strain(1,1)
465 248832 : (*_B[j])[i](1, k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](0);
466 248832 : (*_B[j])[i](1, 4 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](1);
467 248832 : (*_B[j])[i](1, 8 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](2);
468 248832 : (*_B[j])[i](1, 12 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
469 248832 : (-_v2[k] * (*_dxyz_deta[j])[i]);
470 248832 : (*_B[j])[i](1, 16 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
471 248832 : (_v1[k] * (*_dxyz_deta[j])[i]);
472 :
473 : // corresponding to strain(2,2) = 0
474 :
475 : // corresponding to strain(0,1)
476 373248 : (*_B[j])[i](2, k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](0) +
477 248832 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](0));
478 373248 : (*_B[j])[i](2, 4 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](1) +
479 248832 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](1));
480 373248 : (*_B[j])[i](2, 8 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](2) +
481 248832 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](2));
482 124416 : (*_B[j])[i](2, 12 + k) =
483 124416 : 0.25 * _t_points[j](0) * _thickness[i] * -_v2[k] *
484 248832 : (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i] + _dphidxi_map[k][i] * (*_dxyz_deta[j])[i]);
485 124416 : (*_B[j])[i](2, 16 + k) =
486 124416 : 0.25 * _t_points[j](0) * _thickness[i] * _v1[k] *
487 248832 : ((*_dxyz_deta[j])[i] * _dphidxi_map[k][i] + (*_dxyz_dxi[j])[i] * _dphideta_map[k][i]);
488 : }
489 :
490 31104 : _g3_a = _thickness[i] / 4.0 * (_node_normal[2] + _node_normal[3]);
491 31104 : _g3_c = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[1]);
492 31104 : _g3_b = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[3]);
493 31104 : _g3_d = _thickness[i] / 4.0 * (_node_normal[1] + _node_normal[2]);
494 :
495 31104 : _g1_a = 0.5 * ((*_nodes[2]) - (*_nodes[3])) +
496 31104 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[3]);
497 31104 : _g1_c = 0.5 * ((*_nodes[1]) - (*_nodes[0])) +
498 31104 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[1] - _node_normal[0]);
499 31104 : _g2_b = 0.5 * ((*_nodes[3]) - (*_nodes[0])) +
500 31104 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[3] - _node_normal[0]);
501 31104 : _g2_d = 0.5 * ((*_nodes[2]) - (*_nodes[1])) +
502 31104 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[1]);
503 :
504 31104 : updateGVectors(); // for large strain problems
505 :
506 : // corresponding to strain(0,2)
507 124416 : for (unsigned int component = 0; component < 3; component++)
508 : {
509 186624 : (*_B[j])[i](3, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * _g3_a(component);
510 186624 : (*_B[j])[i](3, 3 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * -_g3_a(component);
511 186624 : (*_B[j])[i](3, 1 + component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * _g3_c(component);
512 186624 : (*_B[j])[i](3, component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * -_g3_c(component);
513 : }
514 31104 : (*_B[j])[i](3, 14) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[2];
515 31104 : (*_B[j])[i](3, 18) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[2];
516 31104 : (*_B[j])[i](3, 15) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[3];
517 31104 : (*_B[j])[i](3, 19) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[3];
518 :
519 31104 : (*_B[j])[i](3, 13) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[1];
520 31104 : (*_B[j])[i](3, 17) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * _v1[1];
521 31104 : (*_B[j])[i](3, 12) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[0];
522 31104 : (*_B[j])[i](3, 16) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * _v1[0];
523 :
524 : // corresponding to strain(1,2)
525 124416 : for (unsigned int component = 0; component < 3; component++)
526 : {
527 186624 : (*_B[j])[i](4, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * _g3_d(component);
528 186624 : (*_B[j])[i](4, 1 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * -_g3_d(component);
529 186624 : (*_B[j])[i](4, 3 + component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * _g3_b(component);
530 186624 : (*_B[j])[i](4, component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * -_g3_b(component);
531 : }
532 31104 : (*_B[j])[i](4, 14) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[2];
533 31104 : (*_B[j])[i](4, 18) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[2];
534 31104 : (*_B[j])[i](4, 13) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[1];
535 31104 : (*_B[j])[i](4, 17) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[1];
536 :
537 31104 : (*_B[j])[i](4, 15) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[3];
538 31104 : (*_B[j])[i](4, 19) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[3];
539 31104 : (*_B[j])[i](4, 12) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[0];
540 31104 : (*_B[j])[i](4, 16) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[0];
541 : }
542 : }
543 3888 : }
544 :
545 : void
546 4970 : ADComputeIncrementalShellStrain::computeSolnVector()
547 : {
548 4970 : _soln_vector.zero();
549 :
550 24850 : for (unsigned int j = 0; j < 4; ++j)
551 : {
552 19880 : _soln_disp_index[j].resize(_ndisp);
553 19880 : _soln_rot_index[j].resize(_nrot);
554 :
555 79520 : for (unsigned int i = 0; i < _ndisp; ++i)
556 : {
557 59640 : _soln_disp_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _disp_num[i], 0);
558 59640 : _soln_vector(j + i * _nodes.size()) =
559 59640 : (*_sol)(_soln_disp_index[j][i]) - _sol_old(_soln_disp_index[j][i]);
560 59640 : if (ADReal::do_derivatives)
561 10872 : Moose::derivInsert(
562 : _soln_vector(j + i * _nodes.size()).derivatives(), _soln_disp_index[j][i], 1.);
563 : }
564 :
565 59640 : for (unsigned int i = 0; i < _nrot; ++i)
566 : {
567 39760 : _soln_rot_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _rot_num[i], 0);
568 39760 : _soln_vector(j + 12 + i * _nodes.size()) =
569 39760 : (*_sol)(_soln_rot_index[j][i]) - _sol_old(_soln_rot_index[j][i]);
570 39760 : if (ADReal::do_derivatives)
571 7248 : Moose::derivInsert(
572 : _soln_vector(j + 12 + i * _nodes.size()).derivatives(), _soln_rot_index[j][i], 1.);
573 : }
574 : }
575 4970 : }
|
