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