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 "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("SolidMechanicsApp", ADComputeIncrementalShellStrain);
26 :
27 : InputParameters
28 644 : ADComputeIncrementalShellStrain::validParams()
29 : {
30 644 : InputParameters params = Material::validParams();
31 644 : params.addClassDescription("Compute a small strain increment for the shell.");
32 1288 : params.addRequiredCoupledVar(
33 : "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
34 1288 : params.addRequiredCoupledVar(
35 : "displacements",
36 : "The displacements appropriate for the simulation geometry and coordinate system");
37 1288 : params.addRequiredCoupledVar(
38 : "thickness",
39 : "Thickness of the shell. Can be supplied as either a number or a variable name.");
40 1288 : params.addRequiredParam<std::string>("through_thickness_order",
41 : "Quadrature order in out of plane direction");
42 1288 : params.addParam<bool>(
43 1288 : "large_strain", false, "Set to true to turn on finite strain calculations.");
44 644 : params.addParam<RealVectorValue>("reference_first_local_direction",
45 644 : RealVectorValue(1, 0, 0),
46 : "Vector that is projected onto an element to define the "
47 : "direction for the first local coordinate direction e1.");
48 644 : return params;
49 0 : }
50 :
51 483 : ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputParameters & parameters)
52 : : Material(parameters),
53 483 : _nrot(coupledComponents("rotations")),
54 483 : _ndisp(coupledComponents("displacements")),
55 483 : _rot_num(_nrot),
56 483 : _disp_num(_ndisp),
57 483 : _thickness(coupledValue("thickness")),
58 966 : _large_strain(getParam<bool>("large_strain")),
59 : _strain_increment(),
60 : _total_strain(),
61 : _total_strain_old(),
62 483 : _nonlinear_sys(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0)),
63 483 : _soln_disp_index(4),
64 483 : _soln_rot_index(4),
65 483 : _soln_vector(20, 1),
66 483 : _strain_vector(5, 1),
67 483 : _nodes(4),
68 483 : _node_normal(declareADProperty<RealVectorValue>("node_normal")),
69 966 : _node_normal_old(getMaterialPropertyOldByName<RealVectorValue>("node_normal")),
70 : _dxyz_dxi(),
71 : _dxyz_deta(),
72 : _dxyz_dzeta(),
73 : _dxyz_dxi_old(),
74 : _dxyz_deta_old(),
75 : _dxyz_dzeta_old(),
76 483 : _v1(4),
77 483 : _v2(4),
78 : _B(),
79 : _B_old(),
80 : _ge(),
81 : _ge_old(),
82 : _J_map(),
83 : _J_map_old(),
84 : _local_transformation_matrix(),
85 : _local_transformation_matrix_old(),
86 : _covariant_transformation_matrix(),
87 : _covariant_transformation_matrix_old(),
88 : _contravariant_transformation_matrix(),
89 : _contravariant_transformation_matrix_old(),
90 : _total_global_strain(),
91 483 : _sol(_nonlinear_sys.currentSolution()),
92 483 : _sol_old(_nonlinear_sys.solutionOld()),
93 1932 : _ref_first_local_dir(getParam<RealVectorValue>("reference_first_local_direction"))
94 :
95 : {
96 : // Checking for consistency between length of the provided displacements and rotations vector
97 483 : if (_ndisp != 3 || _nrot != 2)
98 0 : mooseError(
99 : "ADComputeIncrementalShellStrain: The number of variables supplied in 'displacements' "
100 : "must be 3 and that in 'rotations' must be 2.");
101 :
102 483 : if (_mesh.hasSecondOrderElements())
103 0 : mooseError(
104 : "ADComputeIncrementalShellStrain: Shell element is implemented only for linear elements.");
105 :
106 : // Checking if the size of the first local vector is within an acceptable range
107 483 : if (MooseUtils::absoluteFuzzyEqual(_ref_first_local_dir.norm(), 0.0, 1e-3))
108 0 : mooseError(
109 : "The norm of the defined first local axis is less than the acceptable telerance (1e-3)");
110 :
111 : // fetch coupled variables and gradients (as stateful properties if necessary)
112 1932 : for (unsigned int i = 0; i < _ndisp; ++i)
113 : {
114 2898 : MooseVariable * disp_variable = getVar("displacements", i);
115 1449 : _disp_num[i] = disp_variable->number();
116 :
117 1449 : if (i < _nrot)
118 : {
119 1932 : MooseVariable * rot_variable = getVar("rotations", i);
120 966 : _rot_num[i] = rot_variable->number();
121 : }
122 : }
123 :
124 483 : _t_qrule = std::make_unique<libMesh::QGauss>(
125 1449 : 1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
126 483 : _t_points = _t_qrule->get_points();
127 483 : _strain_increment.resize(_t_points.size());
128 483 : _total_strain.resize(_t_points.size());
129 483 : _total_strain_old.resize(_t_points.size());
130 483 : _B.resize(_t_points.size());
131 483 : _B_old.resize(_t_points.size());
132 483 : _ge.resize(_t_points.size());
133 483 : _ge_old.resize(_t_points.size());
134 483 : _J_map.resize(_t_points.size());
135 483 : _J_map_old.resize(_t_points.size());
136 483 : _dxyz_dxi.resize(_t_points.size());
137 483 : _dxyz_deta.resize(_t_points.size());
138 483 : _dxyz_dzeta.resize(_t_points.size());
139 483 : _dxyz_dxi_old.resize(_t_points.size());
140 483 : _dxyz_deta_old.resize(_t_points.size());
141 483 : _dxyz_dzeta_old.resize(_t_points.size());
142 483 : _local_transformation_matrix.resize(_t_points.size());
143 483 : _local_transformation_matrix_old.resize(_t_points.size());
144 483 : _covariant_transformation_matrix.resize(_t_points.size());
145 483 : _covariant_transformation_matrix_old.resize(_t_points.size());
146 483 : _contravariant_transformation_matrix.resize(_t_points.size());
147 483 : _contravariant_transformation_matrix_old.resize(_t_points.size());
148 483 : _total_global_strain.resize(_t_points.size());
149 :
150 483 : _transformation_matrix = &declareADProperty<RankTwoTensor>("transformation_matrix_element");
151 :
152 1449 : for (unsigned int i = 0; i < _t_points.size(); ++i)
153 : {
154 966 : _strain_increment[i] =
155 966 : &declareADProperty<RankTwoTensor>("strain_increment_t_points_" + std::to_string(i));
156 966 : _total_strain[i] =
157 966 : &declareADProperty<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
158 966 : _total_strain_old[i] =
159 966 : &getMaterialPropertyOldByName<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
160 966 : _B[i] = &declareADProperty<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
161 1932 : _B_old[i] = &getMaterialPropertyOldByName<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
162 966 : _ge[i] = &declareADProperty<RankTwoTensor>("ge_t_points_" + std::to_string(i));
163 1932 : _ge_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>("ge_t_points_" + std::to_string(i));
164 966 : _J_map[i] = &declareADProperty<Real>("J_mapping_t_points_" + std::to_string(i));
165 1932 : _J_map_old[i] = &getMaterialPropertyOldByName<Real>("J_mapping_t_points_" + std::to_string(i));
166 966 : _dxyz_dxi[i] = &declareADProperty<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
167 966 : _dxyz_dxi_old[i] =
168 966 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
169 966 : _dxyz_deta[i] = &declareADProperty<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
170 966 : _dxyz_deta_old[i] =
171 966 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
172 966 : _dxyz_dzeta[i] =
173 966 : &declareADProperty<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
174 966 : _dxyz_dzeta_old[i] =
175 966 : &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
176 : // Create rotation matrix and total strain global for output purposes only
177 966 : _local_transformation_matrix[i] =
178 966 : &declareProperty<RankTwoTensor>("local_transformation_t_points_" + std::to_string(i));
179 966 : _local_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
180 966 : "local_transformation_t_points_" + std::to_string(i));
181 966 : _covariant_transformation_matrix[i] =
182 966 : &declareProperty<RankTwoTensor>("covariant_transformation_t_points_" + std::to_string(i));
183 966 : _covariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
184 966 : "covariant_transformation_t_points_" + std::to_string(i));
185 966 : _contravariant_transformation_matrix[i] = &declareProperty<RankTwoTensor>(
186 966 : "contravariant_transformation_t_points_" + std::to_string(i));
187 966 : _contravariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
188 966 : "contravariant_transformation_t_points_" + std::to_string(i));
189 966 : _total_global_strain[i] =
190 1932 : &declareProperty<RankTwoTensor>("total_global_strain_t_points_" + std::to_string(i));
191 : }
192 :
193 : // used later for computing local coordinate system
194 483 : _x2(1) = 1;
195 483 : _x3(2) = 1;
196 483 : _e1(0) = 1;
197 483 : _e2(1) = 1;
198 483 : _e3(2) = 1;
199 483 : }
200 :
201 : void
202 58498 : ADComputeIncrementalShellStrain::initQpStatefulProperties()
203 : {
204 58498 : unsigned int dim = _current_elem->dim();
205 58498 : if ((dim != 2))
206 0 : mooseError(
207 : "ADComputeIncrementalShellStrain: Shell element is implemented only for 2D elements");
208 58498 : if (_current_elem->n_nodes() != 4)
209 0 : mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four nodes.");
210 58498 : if (_qrule->get_points().size() != 4)
211 2 : mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four "
212 : "quadrature points.");
213 :
214 58496 : computeGMatrix();
215 58496 : computeBMatrix();
216 58496 : }
217 :
218 : void
219 37737 : ADComputeIncrementalShellStrain::computeProperties()
220 : {
221 : // quadrature points in isoparametric space
222 37737 : _2d_points = _qrule->get_points(); // would be in 2D
223 :
224 : // derivatives of shape functions (dphidxi, dphideta and dphidzeta) evaluated at quadrature points
225 : // (in isoparametric space).
226 37737 : unsigned int dim = _current_elem->dim();
227 37737 : FEType fe_type(Utility::string_to_enum<Order>("First"),
228 37737 : Utility::string_to_enum<FEFamily>("LAGRANGE"));
229 37737 : auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
230 37737 : _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
231 37737 : _dphideta_map = fe->get_fe_map().get_dphideta_map();
232 37737 : _phi_map = fe->get_fe_map().get_phi_map();
233 :
234 188685 : for (unsigned int i = 0; i < 4; ++i)
235 150948 : _nodes[i] = _current_elem->node_ptr(i);
236 :
237 188685 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
238 : {
239 452844 : for (unsigned int j = 0; j < _t_points.size(); ++j)
240 : {
241 301896 : (*_ge[j])[i] = (*_ge_old[j])[i];
242 301896 : (*_J_map[j])[i] = (*_J_map_old[j])[i];
243 301896 : (*_dxyz_dxi[j])[i] = (*_dxyz_dxi_old[j])[i];
244 301896 : (*_dxyz_deta[j])[i] = (*_dxyz_deta_old[j])[i];
245 301896 : (*_dxyz_dzeta[j])[i] = (*_dxyz_dzeta_old[j])[i];
246 301896 : (*_B[j])[i] = (*_B_old[j])[i];
247 301896 : (*_local_transformation_matrix[j])[i] = (*_local_transformation_matrix_old[j])[i];
248 301896 : (*_covariant_transformation_matrix[j])[i] = (*_covariant_transformation_matrix_old[j])[i];
249 301896 : (*_contravariant_transformation_matrix[j])[i] =
250 301896 : (*_contravariant_transformation_matrix_old[j])[i];
251 : }
252 : }
253 :
254 37737 : computeSolnVector();
255 :
256 37737 : computeNodeNormal();
257 :
258 188685 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
259 : {
260 452844 : for (unsigned int j = 0; j < _t_points.size(); ++j)
261 : {
262 : // compute strain increment in covariant coordinate system using B and _soln_vector
263 1811376 : for (unsigned int temp1 = 0; temp1 < 5; ++temp1)
264 : {
265 1509480 : _strain_vector(temp1) = 0.0;
266 31699080 : for (unsigned int temp2 = 0; temp2 < 20; ++temp2)
267 60379200 : _strain_vector(temp1) += (*_B[j])[i](temp1, temp2) * _soln_vector(temp2);
268 : }
269 301896 : (*_strain_increment[j])[i](0, 0) = _strain_vector(0);
270 301896 : (*_strain_increment[j])[i](1, 1) = _strain_vector(1);
271 603792 : (*_strain_increment[j])[i](0, 1) = _strain_vector(2);
272 301896 : (*_strain_increment[j])[i](0, 2) = _strain_vector(3);
273 301896 : (*_strain_increment[j])[i](1, 2) = _strain_vector(4);
274 301896 : (*_strain_increment[j])[i](1, 0) = (*_strain_increment[j])[i](0, 1);
275 301896 : (*_strain_increment[j])[i](2, 0) = (*_strain_increment[j])[i](0, 2);
276 301896 : (*_strain_increment[j])[i](2, 1) = (*_strain_increment[j])[i](1, 2);
277 :
278 301896 : (*_total_strain[j])[i] = (*_total_strain_old[j])[i] + (*_strain_increment[j])[i];
279 :
280 1207584 : for (unsigned int ii = 0; ii < 3; ++ii)
281 3622752 : for (unsigned int jj = 0; jj < 3; ++jj)
282 2717064 : _unrotated_total_strain(ii, jj) = MetaPhysicL::raw_value((*_total_strain[j])[i](ii, jj));
283 301896 : (*_total_global_strain[j])[i] = (*_contravariant_transformation_matrix[j])[i] *
284 301896 : _unrotated_total_strain *
285 301896 : (*_contravariant_transformation_matrix[j])[i].transpose();
286 : }
287 : }
288 37737 : }
289 :
290 : void
291 58896 : ADComputeIncrementalShellStrain::computeGMatrix()
292 : {
293 : using std::sqrt;
294 :
295 : // quadrature points in isoparametric space
296 58896 : _2d_points = _qrule->get_points(); // would be in 2D
297 :
298 58896 : unsigned int dim = _current_elem->dim();
299 :
300 : // derivatives of shape functions (dphidxi, dphideta and dphidzeta) evaluated at quadrature points
301 : // (in isoparametric space).
302 58896 : FEType fe_type(Utility::string_to_enum<Order>("First"),
303 58896 : Utility::string_to_enum<FEFamily>("LAGRANGE"));
304 58896 : auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
305 58896 : _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
306 58896 : _dphideta_map = fe->get_fe_map().get_dphideta_map();
307 58896 : _phi_map = fe->get_fe_map().get_phi_map();
308 :
309 294480 : for (unsigned int i = 0; i < _nodes.size(); ++i)
310 235584 : _nodes[i] = _current_elem->node_ptr(i);
311 :
312 58896 : ADRealVectorValue x = (*_nodes[1] - *_nodes[0]);
313 58896 : ADRealVectorValue y = (*_nodes[3] - *_nodes[0]);
314 58896 : ADRealVectorValue normal = x.cross(y);
315 58896 : normal /= normal.norm();
316 :
317 294480 : for (unsigned int k = 0; k < 4; ++k)
318 235584 : _node_normal[k] = normal;
319 :
320 58896 : ADRankTwoTensor a;
321 58896 : ADDenseMatrix b(5, 20);
322 : ADRealVectorValue c;
323 58896 : RankTwoTensor d;
324 : RealVectorValue f;
325 176688 : for (unsigned int t = 0; t < _t_points.size(); ++t)
326 : {
327 117792 : (*_strain_increment[t])[_qp] = a;
328 117792 : (*_total_strain[t])[_qp] = a;
329 117792 : (*_B[t])[_qp] = b;
330 117792 : (*_ge[t])[_qp] = a;
331 117792 : (*_J_map[t])[_qp] = 0;
332 117792 : (*_dxyz_dxi[t])[_qp] = c;
333 117792 : (*_dxyz_deta[t])[_qp] = c;
334 117792 : (*_dxyz_dzeta[t])[_qp] = c;
335 117792 : (*_local_transformation_matrix[t])[_qp] = d;
336 117792 : (*_covariant_transformation_matrix[t])[_qp] = d;
337 117792 : (*_contravariant_transformation_matrix[t])[_qp] = d;
338 : }
339 :
340 : // calculating derivatives of shape function in physical space (dphi/dx, dphi/dy, dphi/dz) at
341 : // quadrature points these are g_{i} in Dvorkin's paper
342 294480 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
343 : {
344 706752 : for (unsigned int j = 0; j < _t_points.size(); ++j)
345 : {
346 471168 : (*_dxyz_dxi[j])[i].zero();
347 1884672 : for (unsigned int component = 0; component < 3; ++component)
348 : {
349 1413504 : (*_dxyz_dxi[j])[i](component) = 0.0;
350 1413504 : (*_dxyz_deta[j])[i](component) = 0.0;
351 1413504 : (*_dxyz_dzeta[j])[i](component) = 0.0;
352 :
353 7067520 : for (unsigned int k = 0; k < _nodes.size(); ++k)
354 : {
355 11308032 : (*_dxyz_dxi[j])[i](component) += _dphidxi_map[k][i] * ((*_nodes[k])(component)) +
356 5654016 : _t_points[j](0) / 2.0 * _thickness[i] *
357 11308032 : _dphidxi_map[k][i] * _node_normal[k](component);
358 11308032 : (*_dxyz_deta[j])[i](component) += _dphideta_map[k][i] * ((*_nodes[k])(component)) +
359 5654016 : _t_points[j](0) / 2.0 * _thickness[i] *
360 11308032 : _dphideta_map[k][i] * _node_normal[k](component);
361 5654016 : (*_dxyz_dzeta[j])[i](component) +=
362 11308032 : _thickness[i] * _phi_map[k][i] * _node_normal[k](component) / 2.0;
363 : }
364 : }
365 : }
366 : }
367 :
368 294480 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
369 : {
370 706752 : for (unsigned int j = 0; j < _t_points.size(); ++j)
371 : {
372 : // calculate gij for elasticity tensor
373 471168 : ADRankTwoTensor gmn;
374 471168 : RankTwoTensor J;
375 1884672 : for (unsigned int component = 0; component < 3; ++component)
376 : {
377 2827008 : gmn(0, 0) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dxi[j])[i](component);
378 2827008 : gmn(1, 1) += (*_dxyz_deta[j])[i](component) * (*_dxyz_deta[j])[i](component);
379 2827008 : gmn(2, 2) += (*_dxyz_dzeta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
380 2827008 : gmn(0, 1) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_deta[j])[i](component);
381 2827008 : gmn(0, 2) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
382 2827008 : gmn(1, 2) += (*_dxyz_deta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
383 :
384 : // why are we dropping derivatives here?!
385 1413504 : J(0, component) = MetaPhysicL::raw_value((*_dxyz_dxi[j])[i](component));
386 1413504 : J(1, component) = MetaPhysicL::raw_value((*_dxyz_deta[j])[i](component));
387 1413504 : J(2, component) = MetaPhysicL::raw_value((*_dxyz_dzeta[j])[i](component));
388 : }
389 471168 : gmn(1, 0) = gmn(0, 1);
390 471168 : gmn(2, 0) = gmn(0, 2);
391 471168 : gmn(2, 1) = gmn(1, 2);
392 :
393 471168 : ADRankTwoTensor gmninv_temp = gmn.inverse();
394 471168 : (*_J_map[j])[i] = sqrt(gmn.det());
395 471168 : (*_covariant_transformation_matrix[j])[i] = J;
396 :
397 471168 : (*_contravariant_transformation_matrix[j])[i] =
398 471168 : (*_covariant_transformation_matrix[j])[i].inverse();
399 :
400 471168 : Real normx = sqrt(J(0, 0) * J(0, 0) + J(0, 1) * J(0, 1) + J(0, 2) * J(0, 2));
401 471168 : Real normy = sqrt(J(1, 0) * J(1, 0) + J(1, 1) * J(1, 1) + J(1, 2) * J(1, 2));
402 471168 : Real normz = sqrt(J(2, 0) * J(2, 0) + J(2, 1) * J(2, 1) + J(2, 2) * J(2, 2));
403 :
404 471168 : J(0, 0) /= normx;
405 471168 : J(0, 1) /= normx;
406 471168 : J(0, 2) /= normx;
407 :
408 471168 : J(1, 0) /= normy;
409 471168 : J(1, 1) /= normy;
410 471168 : J(1, 2) /= normy;
411 :
412 471168 : J(2, 0) /= normz;
413 471168 : J(2, 1) /= normz;
414 471168 : J(2, 2) /= normz;
415 :
416 : // _transformation_matrix is an AD property, but we're not setting the derivatives!
417 471168 : (*_transformation_matrix)[i] = J;
418 :
419 : // compute element's local coordinate
420 471168 : computeLocalCoordinates();
421 :
422 : // calculate the local transformation matrix to be used to map the global stresses to the
423 : // local element coordinate
424 471168 : const auto local_rotation_mat = ADRankTwoTensor::initializeFromRows(_e1, _e2, _e3);
425 :
426 1884672 : for (const auto ii : make_range(Moose::dim))
427 5654016 : for (const auto jj : make_range(Moose::dim))
428 4240512 : (*_local_transformation_matrix[j])[i](ii, jj) =
429 4240512 : MetaPhysicL::raw_value(local_rotation_mat(ii, jj));
430 :
431 : // Calculate the contravariant base vectors g^0, g^1, g^2
432 : // The base vectors for the strain tensor in the convected coordinate
433 : // are g_0, g_1, g_2 (g_i=dx/dr_i)
434 : // The following contravariant base vectors have the property:
435 : // g^i*g_j= 1 if {i=j} otherwise g^i*g_j= 0
436 :
437 471168 : const auto denom = (*_dxyz_dxi[j])[i] * (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]);
438 471168 : const auto gi0 = (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]) / denom;
439 471168 : const auto gi1 = (*_dxyz_dzeta[j])[i].cross((*_dxyz_dxi[j])[i]) / denom;
440 471168 : const auto gi2 = (*_dxyz_dxi[j])[i].cross((*_dxyz_deta[j])[i]) / denom;
441 :
442 : // Calculate the rotation matrix for the elasticity tensor that maps
443 : // the strain tensor (with g_1, g2_, g_3 base vectors) to
444 : // the stress tensor (with base vectors e1, e2, e3)
445 :
446 471168 : (*_ge[j])[i](0, 0) = gi0 * _e1;
447 471168 : (*_ge[j])[i](0, 1) = gi0 * _e2;
448 471168 : (*_ge[j])[i](0, 2) = gi0 * _e3;
449 471168 : (*_ge[j])[i](1, 0) = gi1 * _e1;
450 471168 : (*_ge[j])[i](1, 1) = gi1 * _e2;
451 471168 : (*_ge[j])[i](1, 2) = gi1 * _e3;
452 471168 : (*_ge[j])[i](2, 0) = gi2 * _e1;
453 471168 : (*_ge[j])[i](2, 1) = gi2 * _e2;
454 471168 : (*_ge[j])[i](2, 2) = gi2 * _e3;
455 : }
456 : }
457 58896 : }
458 :
459 : void
460 531797 : ADComputeIncrementalShellStrain::computeLocalCoordinates()
461 : {
462 : using std::abs;
463 :
464 : // default option for the first local vector:the global X-axis is projected on the shell plane
465 : // alternatively the reference first local vector provided by the user can be used to define the
466 : // 1st local axis
467 :
468 : // All nodes in an element have the same normal vector. Therefore, here, we use only the normal
469 : // vecor of the first node as "the element's normal vector"
470 531797 : _e3 = _node_normal[0];
471 :
472 531797 : _e1 = _ref_first_local_dir;
473 :
474 531797 : _e1 /= _e1.norm();
475 :
476 : // The reference first local axis and the normal are considered parallel if the angle between them
477 : // is less than 0.05 degrees
478 531797 : if (MooseUtils::absoluteFuzzyEqual(abs(_e1 * _e3), 1.0, 0.05 * libMesh::pi / 180.0))
479 : {
480 0 : mooseError("The reference first local axis must not be perpendicular to any of the shell "
481 : "elements ");
482 : }
483 :
484 : // we project the reference first local vector on the shell element and calculate the in-plane e1
485 : // and e2 vectors
486 531797 : _e2 = _e3.cross(_e1);
487 531797 : _e2 /= _e2.norm();
488 :
489 531797 : _e1 = _e2.cross(_e3);
490 531797 : _e1 /= _e1.norm();
491 531797 : }
492 :
493 : void
494 37737 : ADComputeIncrementalShellStrain::computeNodeNormal()
495 : {
496 188685 : for (unsigned int k = 0; k < _nodes.size(); ++k)
497 150948 : _node_normal[k] = _node_normal_old[k];
498 37737 : }
499 :
500 : void
501 60629 : ADComputeIncrementalShellStrain::computeBMatrix()
502 : {
503 : // compute nodal local axis
504 60629 : computeLocalCoordinates();
505 303145 : for (unsigned int k = 0; k < _nodes.size(); ++k)
506 : {
507 : _v1[k] = _e1;
508 : _v2[k] = _e2;
509 : }
510 :
511 : // compute B matrix rows correspond to [ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, uz1, uz2, uz3,
512 : // uz4, a1, a2, a3, a4, b1, b2, b3, b4]
513 303145 : for (unsigned int i = 0; i < _2d_points.size(); ++i)
514 : {
515 727548 : for (unsigned int j = 0; j < _t_points.size(); ++j)
516 : {
517 485032 : (*_B[j])[i].resize(5, 20);
518 485032 : (*_B[j])[i].zero();
519 2425160 : for (unsigned int k = 0; k < _nodes.size(); ++k)
520 : {
521 : // corresponding to strain(0,0)
522 3880256 : (*_B[j])[i](0, k) += _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](0);
523 3880256 : (*_B[j])[i](0, 4 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](1);
524 3880256 : (*_B[j])[i](0, 8 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](2);
525 3880256 : (*_B[j])[i](0, 12 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
526 3880256 : (-_v2[k] * (*_dxyz_dxi[j])[i]);
527 3880256 : (*_B[j])[i](0, 16 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
528 3880256 : (_v1[k] * (*_dxyz_dxi[j])[i]);
529 :
530 : // corresponding to strain(1,1)
531 3880256 : (*_B[j])[i](1, k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](0);
532 3880256 : (*_B[j])[i](1, 4 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](1);
533 3880256 : (*_B[j])[i](1, 8 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](2);
534 3880256 : (*_B[j])[i](1, 12 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
535 3880256 : (-_v2[k] * (*_dxyz_deta[j])[i]);
536 3880256 : (*_B[j])[i](1, 16 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
537 3880256 : (_v1[k] * (*_dxyz_deta[j])[i]);
538 :
539 : // corresponding to strain(2,2) = 0
540 :
541 : // corresponding to strain(0,1)
542 5820384 : (*_B[j])[i](2, k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](0) +
543 3880256 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](0));
544 5820384 : (*_B[j])[i](2, 4 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](1) +
545 3880256 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](1));
546 5820384 : (*_B[j])[i](2, 8 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](2) +
547 3880256 : _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](2));
548 1940128 : (*_B[j])[i](2, 12 + k) =
549 1940128 : 0.25 * _t_points[j](0) * _thickness[i] * -_v2[k] *
550 3880256 : (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i] + _dphidxi_map[k][i] * (*_dxyz_deta[j])[i]);
551 1940128 : (*_B[j])[i](2, 16 + k) =
552 1940128 : 0.25 * _t_points[j](0) * _thickness[i] * _v1[k] *
553 3880256 : ((*_dxyz_deta[j])[i] * _dphidxi_map[k][i] + (*_dxyz_dxi[j])[i] * _dphideta_map[k][i]);
554 : }
555 :
556 485032 : _g3_a = _thickness[i] / 4.0 * (_node_normal[2] + _node_normal[3]);
557 485032 : _g3_c = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[1]);
558 485032 : _g3_b = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[3]);
559 485032 : _g3_d = _thickness[i] / 4.0 * (_node_normal[1] + _node_normal[2]);
560 :
561 485032 : _g1_a = 0.5 * ((*_nodes[2]) - (*_nodes[3])) +
562 485032 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[3]);
563 485032 : _g1_c = 0.5 * ((*_nodes[1]) - (*_nodes[0])) +
564 485032 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[1] - _node_normal[0]);
565 485032 : _g2_b = 0.5 * ((*_nodes[3]) - (*_nodes[0])) +
566 485032 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[3] - _node_normal[0]);
567 485032 : _g2_d = 0.5 * ((*_nodes[2]) - (*_nodes[1])) +
568 485032 : _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[1]);
569 :
570 485032 : updateGVectors(); // for large strain problems
571 :
572 : // corresponding to strain(0,2)
573 1940128 : for (unsigned int component = 0; component < 3; component++)
574 : {
575 2910192 : (*_B[j])[i](3, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * _g3_a(component);
576 2910192 : (*_B[j])[i](3, 3 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * -_g3_a(component);
577 2910192 : (*_B[j])[i](3, 1 + component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * _g3_c(component);
578 2910192 : (*_B[j])[i](3, component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * -_g3_c(component);
579 : }
580 485032 : (*_B[j])[i](3, 14) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[2];
581 485032 : (*_B[j])[i](3, 18) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[2];
582 485032 : (*_B[j])[i](3, 15) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[3];
583 485032 : (*_B[j])[i](3, 19) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[3];
584 :
585 485032 : (*_B[j])[i](3, 13) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[1];
586 485032 : (*_B[j])[i](3, 17) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * _v1[1];
587 485032 : (*_B[j])[i](3, 12) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[0];
588 485032 : (*_B[j])[i](3, 16) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * _v1[0];
589 :
590 : // corresponding to strain(1,2)
591 1940128 : for (unsigned int component = 0; component < 3; component++)
592 : {
593 2910192 : (*_B[j])[i](4, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * _g3_d(component);
594 2910192 : (*_B[j])[i](4, 1 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * -_g3_d(component);
595 2910192 : (*_B[j])[i](4, 3 + component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * _g3_b(component);
596 2910192 : (*_B[j])[i](4, component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * -_g3_b(component);
597 : }
598 485032 : (*_B[j])[i](4, 14) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[2];
599 485032 : (*_B[j])[i](4, 18) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[2];
600 485032 : (*_B[j])[i](4, 13) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[1];
601 485032 : (*_B[j])[i](4, 17) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[1];
602 :
603 485032 : (*_B[j])[i](4, 15) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[3];
604 485032 : (*_B[j])[i](4, 19) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[3];
605 485032 : (*_B[j])[i](4, 12) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[0];
606 485032 : (*_B[j])[i](4, 16) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[0];
607 : }
608 : }
609 60629 : }
610 :
611 : void
612 39870 : ADComputeIncrementalShellStrain::computeSolnVector()
613 : {
614 39870 : _soln_vector.zero();
615 :
616 199350 : for (unsigned int j = 0; j < 4; ++j)
617 : {
618 159480 : _soln_disp_index[j].resize(_ndisp);
619 159480 : _soln_rot_index[j].resize(_nrot);
620 :
621 637920 : for (unsigned int i = 0; i < _ndisp; ++i)
622 : {
623 478440 : _soln_disp_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _disp_num[i], 0);
624 478440 : _soln_vector(j + i * _nodes.size()) =
625 478440 : (*_sol)(_soln_disp_index[j][i]) - _sol_old(_soln_disp_index[j][i]);
626 478440 : if (ADReal::do_derivatives)
627 186300 : Moose::derivInsert(
628 : _soln_vector(j + i * _nodes.size()).derivatives(), _soln_disp_index[j][i], 1.);
629 : }
630 :
631 478440 : for (unsigned int i = 0; i < _nrot; ++i)
632 : {
633 318960 : _soln_rot_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _rot_num[i], 0);
634 318960 : _soln_vector(j + 12 + i * _nodes.size()) =
635 318960 : (*_sol)(_soln_rot_index[j][i]) - _sol_old(_soln_rot_index[j][i]);
636 318960 : if (ADReal::do_derivatives)
637 124200 : Moose::derivInsert(
638 : _soln_vector(j + 12 + i * _nodes.size()).derivatives(), _soln_rot_index[j][i], 1.);
639 : }
640 : }
641 39870 : }
|
