https://mooseframework.inl.gov
ADComputeIncrementalShellStrain.C
Go to the documentation of this file.
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 
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 
26 
29 {
31  params.addClassDescription("Compute a small strain increment for the shell.");
32  params.addRequiredCoupledVar(
33  "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
34  params.addRequiredCoupledVar(
35  "displacements",
36  "The displacements appropriate for the simulation geometry and coordinate system");
37  params.addRequiredCoupledVar(
38  "thickness",
39  "Thickness of the shell. Can be supplied as either a number or a variable name.");
40  params.addRequiredParam<std::string>("through_thickness_order",
41  "Quadrature order in out of plane direction");
42  params.addParam<bool>(
43  "large_strain", false, "Set to true to turn on finite strain calculations.");
44  params.addParam<RealVectorValue>("reference_first_local_direction",
45  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  return params;
49 }
50 
52  : Material(parameters),
53  _nrot(coupledComponents("rotations")),
54  _ndisp(coupledComponents("displacements")),
55  _rot_num(_nrot),
56  _disp_num(_ndisp),
57  _thickness(coupledValue("thickness")),
58  _large_strain(getParam<bool>("large_strain")),
59  _strain_increment(),
60  _total_strain(),
61  _total_strain_old(),
62  _nonlinear_sys(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0)),
63  _soln_disp_index(4),
64  _soln_rot_index(4),
65  _soln_vector(20, 1),
66  _strain_vector(5, 1),
67  _nodes(4),
68  _node_normal(declareADProperty<RealVectorValue>("node_normal")),
69  _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  _v1(4),
77  _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  _sol(_nonlinear_sys.currentSolution()),
92  _sol_old(_nonlinear_sys.solutionOld()),
93  _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  if (_ndisp != 3 || _nrot != 2)
98  mooseError(
99  "ADComputeIncrementalShellStrain: The number of variables supplied in 'displacements' "
100  "must be 3 and that in 'rotations' must be 2.");
101 
103  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
108  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  for (unsigned int i = 0; i < _ndisp; ++i)
113  {
114  MooseVariable * disp_variable = getVar("displacements", i);
115  _disp_num[i] = disp_variable->number();
116 
117  if (i < _nrot)
118  {
119  MooseVariable * rot_variable = getVar("rotations", i);
120  _rot_num[i] = rot_variable->number();
121  }
122  }
123 
124  _t_qrule = std::make_unique<libMesh::QGauss>(
125  1, Utility::string_to_enum<Order>(getParam<std::string>("through_thickness_order")));
126  _t_points = _t_qrule->get_points();
127  _strain_increment.resize(_t_points.size());
128  _total_strain.resize(_t_points.size());
129  _total_strain_old.resize(_t_points.size());
130  _B.resize(_t_points.size());
131  _B_old.resize(_t_points.size());
132  _ge.resize(_t_points.size());
133  _ge_old.resize(_t_points.size());
134  _J_map.resize(_t_points.size());
135  _J_map_old.resize(_t_points.size());
136  _dxyz_dxi.resize(_t_points.size());
137  _dxyz_deta.resize(_t_points.size());
138  _dxyz_dzeta.resize(_t_points.size());
139  _dxyz_dxi_old.resize(_t_points.size());
140  _dxyz_deta_old.resize(_t_points.size());
141  _dxyz_dzeta_old.resize(_t_points.size());
148  _total_global_strain.resize(_t_points.size());
149 
150  _transformation_matrix = &declareADProperty<RankTwoTensor>("transformation_matrix_element");
151 
152  for (unsigned int i = 0; i < _t_points.size(); ++i)
153  {
154  _strain_increment[i] =
155  &declareADProperty<RankTwoTensor>("strain_increment_t_points_" + std::to_string(i));
156  _total_strain[i] =
157  &declareADProperty<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
158  _total_strain_old[i] =
159  &getMaterialPropertyOldByName<RankTwoTensor>("total_strain_t_points_" + std::to_string(i));
160  _B[i] = &declareADProperty<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
161  _B_old[i] = &getMaterialPropertyOldByName<DenseMatrix<Real>>("B_t_points_" + std::to_string(i));
162  _ge[i] = &declareADProperty<RankTwoTensor>("ge_t_points_" + std::to_string(i));
163  _ge_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>("ge_t_points_" + std::to_string(i));
164  _J_map[i] = &declareADProperty<Real>("J_mapping_t_points_" + std::to_string(i));
165  _J_map_old[i] = &getMaterialPropertyOldByName<Real>("J_mapping_t_points_" + std::to_string(i));
166  _dxyz_dxi[i] = &declareADProperty<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
167  _dxyz_dxi_old[i] =
168  &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dxi_t_points_" + std::to_string(i));
169  _dxyz_deta[i] = &declareADProperty<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
170  _dxyz_deta_old[i] =
171  &getMaterialPropertyOldByName<RealVectorValue>("dxyz_deta_t_points_" + std::to_string(i));
172  _dxyz_dzeta[i] =
173  &declareADProperty<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
174  _dxyz_dzeta_old[i] =
175  &getMaterialPropertyOldByName<RealVectorValue>("dxyz_dzeta_t_points_" + std::to_string(i));
176  // Create rotation matrix and total strain global for output purposes only
178  &declareProperty<RankTwoTensor>("local_transformation_t_points_" + std::to_string(i));
179  _local_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
180  "local_transformation_t_points_" + std::to_string(i));
182  &declareProperty<RankTwoTensor>("covariant_transformation_t_points_" + std::to_string(i));
183  _covariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
184  "covariant_transformation_t_points_" + std::to_string(i));
185  _contravariant_transformation_matrix[i] = &declareProperty<RankTwoTensor>(
186  "contravariant_transformation_t_points_" + std::to_string(i));
187  _contravariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName<RankTwoTensor>(
188  "contravariant_transformation_t_points_" + std::to_string(i));
190  &declareProperty<RankTwoTensor>("total_global_strain_t_points_" + std::to_string(i));
191  }
192 
193  // used later for computing local coordinate system
194  _x2(1) = 1;
195  _x3(2) = 1;
196  _e1(0) = 1;
197  _e2(1) = 1;
198  _e3(2) = 1;
199 }
200 
201 void
203 {
204  unsigned int dim = _current_elem->dim();
205  if ((dim != 2))
206  mooseError(
207  "ADComputeIncrementalShellStrain: Shell element is implemented only for 2D elements");
208  if (_current_elem->n_nodes() != 4)
209  mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four nodes.");
210  if (_qrule->get_points().size() != 4)
211  mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four "
212  "quadrature points.");
213 
214  computeGMatrix();
215  computeBMatrix();
216 }
217 
218 void
220 {
221  // quadrature points in isoparametric space
222  _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  unsigned int dim = _current_elem->dim();
227  FEType fe_type(Utility::string_to_enum<Order>("First"),
228  Utility::string_to_enum<FEFamily>("LAGRANGE"));
229  auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
230  _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
231  _dphideta_map = fe->get_fe_map().get_dphideta_map();
232  _phi_map = fe->get_fe_map().get_phi_map();
233 
234  for (unsigned int i = 0; i < 4; ++i)
235  _nodes[i] = _current_elem->node_ptr(i);
236 
237  for (unsigned int i = 0; i < _2d_points.size(); ++i)
238  {
239  for (unsigned int j = 0; j < _t_points.size(); ++j)
240  {
241  (*_ge[j])[i] = (*_ge_old[j])[i];
242  (*_J_map[j])[i] = (*_J_map_old[j])[i];
243  (*_dxyz_dxi[j])[i] = (*_dxyz_dxi_old[j])[i];
244  (*_dxyz_deta[j])[i] = (*_dxyz_deta_old[j])[i];
245  (*_dxyz_dzeta[j])[i] = (*_dxyz_dzeta_old[j])[i];
246  (*_B[j])[i] = (*_B_old[j])[i];
251  }
252  }
253 
255 
257 
258  for (unsigned int i = 0; i < _2d_points.size(); ++i)
259  {
260  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  for (unsigned int temp1 = 0; temp1 < 5; ++temp1)
264  {
265  _strain_vector(temp1) = 0.0;
266  for (unsigned int temp2 = 0; temp2 < 20; ++temp2)
267  _strain_vector(temp1) += (*_B[j])[i](temp1, temp2) * _soln_vector(temp2);
268  }
269  (*_strain_increment[j])[i](0, 0) = _strain_vector(0);
270  (*_strain_increment[j])[i](1, 1) = _strain_vector(1);
271  (*_strain_increment[j])[i](0, 1) = _strain_vector(2);
272  (*_strain_increment[j])[i](0, 2) = _strain_vector(3);
273  (*_strain_increment[j])[i](1, 2) = _strain_vector(4);
274  (*_strain_increment[j])[i](1, 0) = (*_strain_increment[j])[i](0, 1);
275  (*_strain_increment[j])[i](2, 0) = (*_strain_increment[j])[i](0, 2);
276  (*_strain_increment[j])[i](2, 1) = (*_strain_increment[j])[i](1, 2);
277 
278  (*_total_strain[j])[i] = (*_total_strain_old[j])[i] + (*_strain_increment[j])[i];
279 
280  for (unsigned int ii = 0; ii < 3; ++ii)
281  for (unsigned int jj = 0; jj < 3; ++jj)
285  (*_contravariant_transformation_matrix[j])[i].transpose();
286  }
287  }
288 }
289 
290 void
292 {
293  // quadrature points in isoparametric space
294  _2d_points = _qrule->get_points(); // would be in 2D
295 
296  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  FEType fe_type(Utility::string_to_enum<Order>("First"),
301  Utility::string_to_enum<FEFamily>("LAGRANGE"));
302  auto & fe = _fe_problem.assembly(_tid, _nonlinear_sys.number()).getFE(fe_type, dim);
303  _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
304  _dphideta_map = fe->get_fe_map().get_dphideta_map();
305  _phi_map = fe->get_fe_map().get_phi_map();
306 
307  for (unsigned int i = 0; i < _nodes.size(); ++i)
308  _nodes[i] = _current_elem->node_ptr(i);
309 
310  ADRealVectorValue x = (*_nodes[1] - *_nodes[0]);
311  ADRealVectorValue y = (*_nodes[3] - *_nodes[0]);
312  ADRealVectorValue normal = x.cross(y);
313  normal /= normal.norm();
314 
315  for (unsigned int k = 0; k < 4; ++k)
316  _node_normal[k] = normal;
317 
319  ADDenseMatrix b(5, 20);
323  for (unsigned int t = 0; t < _t_points.size(); ++t)
324  {
325  (*_strain_increment[t])[_qp] = a;
326  (*_total_strain[t])[_qp] = a;
327  (*_B[t])[_qp] = b;
328  (*_ge[t])[_qp] = a;
329  (*_J_map[t])[_qp] = 0;
330  (*_dxyz_dxi[t])[_qp] = c;
331  (*_dxyz_deta[t])[_qp] = c;
332  (*_dxyz_dzeta[t])[_qp] = c;
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  for (unsigned int i = 0; i < _2d_points.size(); ++i)
341  {
342  for (unsigned int j = 0; j < _t_points.size(); ++j)
343  {
344  (*_dxyz_dxi[j])[i].zero();
345  for (unsigned int component = 0; component < 3; ++component)
346  {
347  (*_dxyz_dxi[j])[i](component) = 0.0;
348  (*_dxyz_deta[j])[i](component) = 0.0;
349  (*_dxyz_dzeta[j])[i](component) = 0.0;
350 
351  for (unsigned int k = 0; k < _nodes.size(); ++k)
352  {
353  (*_dxyz_dxi[j])[i](component) += _dphidxi_map[k][i] * ((*_nodes[k])(component)) +
354  _t_points[j](0) / 2.0 * _thickness[i] *
356  (*_dxyz_deta[j])[i](component) += _dphideta_map[k][i] * ((*_nodes[k])(component)) +
357  _t_points[j](0) / 2.0 * _thickness[i] *
359  (*_dxyz_dzeta[j])[i](component) +=
360  _thickness[i] * _phi_map[k][i] * _node_normal[k](component) / 2.0;
361  }
362  }
363  }
364  }
365 
366  for (unsigned int i = 0; i < _2d_points.size(); ++i)
367  {
368  for (unsigned int j = 0; j < _t_points.size(); ++j)
369  {
370  // calculate gij for elasticity tensor
371  ADRankTwoTensor gmn;
372  RankTwoTensor J;
373  for (unsigned int component = 0; component < 3; ++component)
374  {
375  gmn(0, 0) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dxi[j])[i](component);
376  gmn(1, 1) += (*_dxyz_deta[j])[i](component) * (*_dxyz_deta[j])[i](component);
377  gmn(2, 2) += (*_dxyz_dzeta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
378  gmn(0, 1) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_deta[j])[i](component);
379  gmn(0, 2) += (*_dxyz_dxi[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
380  gmn(1, 2) += (*_dxyz_deta[j])[i](component) * (*_dxyz_dzeta[j])[i](component);
381 
382  // why are we dropping derivatives here?!
386  }
387  gmn(1, 0) = gmn(0, 1);
388  gmn(2, 0) = gmn(0, 2);
389  gmn(2, 1) = gmn(1, 2);
390 
391  ADRankTwoTensor gmninv_temp = gmn.inverse();
392  (*_J_map[j])[i] = std::sqrt(gmn.det());
394 
397 
398  Real normx = std::sqrt(J(0, 0) * J(0, 0) + J(0, 1) * J(0, 1) + J(0, 2) * J(0, 2));
399  Real normy = std::sqrt(J(1, 0) * J(1, 0) + J(1, 1) * J(1, 1) + J(1, 2) * J(1, 2));
400  Real normz = std::sqrt(J(2, 0) * J(2, 0) + J(2, 1) * J(2, 1) + J(2, 2) * J(2, 2));
401 
402  J(0, 0) /= normx;
403  J(0, 1) /= normx;
404  J(0, 2) /= normx;
405 
406  J(1, 0) /= normy;
407  J(1, 1) /= normy;
408  J(1, 2) /= normy;
409 
410  J(2, 0) /= normz;
411  J(2, 1) /= normz;
412  J(2, 2) /= normz;
413 
414  // _transformation_matrix is an AD property, but we're not setting the derivatives!
415  (*_transformation_matrix)[i] = J;
416 
417  // compute element's local coordinate
419 
420  // calculate the local transformation matrix to be used to map the global stresses to the
421  // local element coordinate
422  const auto local_rotation_mat = ADRankTwoTensor::initializeFromRows(_e1, _e2, _e3);
423 
424  for (const auto ii : make_range(Moose::dim))
425  for (const auto jj : make_range(Moose::dim))
426  (*_local_transformation_matrix[j])[i](ii, jj) =
427  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  const auto denom = (*_dxyz_dxi[j])[i] * (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]);
436  const auto gi0 = (*_dxyz_deta[j])[i].cross((*_dxyz_dzeta[j])[i]) / denom;
437  const auto gi1 = (*_dxyz_dzeta[j])[i].cross((*_dxyz_dxi[j])[i]) / denom;
438  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  (*_ge[j])[i](0, 0) = gi0 * _e1;
445  (*_ge[j])[i](0, 1) = gi0 * _e2;
446  (*_ge[j])[i](0, 2) = gi0 * _e3;
447  (*_ge[j])[i](1, 0) = gi1 * _e1;
448  (*_ge[j])[i](1, 1) = gi1 * _e2;
449  (*_ge[j])[i](1, 2) = gi1 * _e3;
450  (*_ge[j])[i](2, 0) = gi2 * _e1;
451  (*_ge[j])[i](2, 1) = gi2 * _e2;
452  (*_ge[j])[i](2, 2) = gi2 * _e3;
453  }
454  }
455 }
456 
457 void
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  _e3 = _node_normal[0];
467 
469 
470  _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  if (MooseUtils::absoluteFuzzyEqual(std::abs(_e1 * _e3), 1.0, 0.05 * libMesh::pi / 180.0))
475  {
476  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  _e2 = _e3.cross(_e1);
483  _e2 /= _e2.norm();
484 
485  _e1 = _e2.cross(_e3);
486  _e1 /= _e1.norm();
487 }
488 
489 void
491 {
492  for (unsigned int k = 0; k < _nodes.size(); ++k)
494 }
495 
496 void
498 {
499  // compute nodal local axis
501  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  for (unsigned int i = 0; i < _2d_points.size(); ++i)
510  {
511  for (unsigned int j = 0; j < _t_points.size(); ++j)
512  {
513  (*_B[j])[i].resize(5, 20);
514  (*_B[j])[i].zero();
515  for (unsigned int k = 0; k < _nodes.size(); ++k)
516  {
517  // corresponding to strain(0,0)
518  (*_B[j])[i](0, k) += _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](0);
519  (*_B[j])[i](0, 4 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](1);
520  (*_B[j])[i](0, 8 + k) = _dphidxi_map[k][i] * (*_dxyz_dxi[j])[i](2);
521  (*_B[j])[i](0, 12 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
522  (-_v2[k] * (*_dxyz_dxi[j])[i]);
523  (*_B[j])[i](0, 16 + k) = _dphidxi_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
524  (_v1[k] * (*_dxyz_dxi[j])[i]);
525 
526  // corresponding to strain(1,1)
527  (*_B[j])[i](1, k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](0);
528  (*_B[j])[i](1, 4 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](1);
529  (*_B[j])[i](1, 8 + k) = _dphideta_map[k][i] * (*_dxyz_deta[j])[i](2);
530  (*_B[j])[i](1, 12 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
531  (-_v2[k] * (*_dxyz_deta[j])[i]);
532  (*_B[j])[i](1, 16 + k) = _dphideta_map[k][i] * _t_points[j](0) / 2.0 * _thickness[i] *
533  (_v1[k] * (*_dxyz_deta[j])[i]);
534 
535  // corresponding to strain(2,2) = 0
536 
537  // corresponding to strain(0,1)
538  (*_B[j])[i](2, k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](0) +
539  _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](0));
540  (*_B[j])[i](2, 4 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](1) +
541  _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](1));
542  (*_B[j])[i](2, 8 + k) = 0.5 * (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i](2) +
543  _dphidxi_map[k][i] * (*_dxyz_deta[j])[i](2));
544  (*_B[j])[i](2, 12 + k) =
545  0.25 * _t_points[j](0) * _thickness[i] * -_v2[k] *
546  (_dphideta_map[k][i] * (*_dxyz_dxi[j])[i] + _dphidxi_map[k][i] * (*_dxyz_deta[j])[i]);
547  (*_B[j])[i](2, 16 + k) =
548  0.25 * _t_points[j](0) * _thickness[i] * _v1[k] *
549  ((*_dxyz_deta[j])[i] * _dphidxi_map[k][i] + (*_dxyz_dxi[j])[i] * _dphideta_map[k][i]);
550  }
551 
552  _g3_a = _thickness[i] / 4.0 * (_node_normal[2] + _node_normal[3]);
553  _g3_c = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[1]);
554  _g3_b = _thickness[i] / 4.0 * (_node_normal[0] + _node_normal[3]);
555  _g3_d = _thickness[i] / 4.0 * (_node_normal[1] + _node_normal[2]);
556 
557  _g1_a = 0.5 * ((*_nodes[2]) - (*_nodes[3])) +
558  _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[3]);
559  _g1_c = 0.5 * ((*_nodes[1]) - (*_nodes[0])) +
560  _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[1] - _node_normal[0]);
561  _g2_b = 0.5 * ((*_nodes[3]) - (*_nodes[0])) +
562  _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[3] - _node_normal[0]);
563  _g2_d = 0.5 * ((*_nodes[2]) - (*_nodes[1])) +
564  _t_points[j](0) / 4.0 * _thickness[i] * (_node_normal[2] - _node_normal[1]);
565 
566  updateGVectors(); // for large strain problems
567 
568  // corresponding to strain(0,2)
569  for (unsigned int component = 0; component < 3; component++)
570  {
571  (*_B[j])[i](3, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * _g3_a(component);
572  (*_B[j])[i](3, 3 + component * 4) = 0.125 * (1.0 + _2d_points[i](1)) * -_g3_a(component);
573  (*_B[j])[i](3, 1 + component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * _g3_c(component);
574  (*_B[j])[i](3, component * 4) = 0.125 * (1.0 - _2d_points[i](1)) * -_g3_c(component);
575  }
576  (*_B[j])[i](3, 14) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[2];
577  (*_B[j])[i](3, 18) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[2];
578  (*_B[j])[i](3, 15) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * -_v2[3];
579  (*_B[j])[i](3, 19) = 0.125 * (1.0 + _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_a * _v1[3];
580 
581  (*_B[j])[i](3, 13) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[1];
582  (*_B[j])[i](3, 17) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * _v1[1];
583  (*_B[j])[i](3, 12) = 0.125 * (1.0 - _2d_points[i](1)) * 0.5 * _thickness[i] * _g1_c * -_v2[0];
584  (*_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  for (unsigned int component = 0; component < 3; component++)
588  {
589  (*_B[j])[i](4, 2 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * _g3_d(component);
590  (*_B[j])[i](4, 1 + component * 4) = 0.125 * (1.0 + _2d_points[i](0)) * -_g3_d(component);
591  (*_B[j])[i](4, 3 + component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * _g3_b(component);
592  (*_B[j])[i](4, component * 4) = 0.125 * (1.0 - _2d_points[i](0)) * -_g3_b(component);
593  }
594  (*_B[j])[i](4, 14) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[2];
595  (*_B[j])[i](4, 18) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[2];
596  (*_B[j])[i](4, 13) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * -_v2[1];
597  (*_B[j])[i](4, 17) = 0.125 * (1.0 + _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_d * _v1[1];
598 
599  (*_B[j])[i](4, 15) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[3];
600  (*_B[j])[i](4, 19) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[3];
601  (*_B[j])[i](4, 12) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * -_v2[0];
602  (*_B[j])[i](4, 16) = 0.125 * (1.0 - _2d_points[i](0)) * 0.5 * _thickness[i] * _g2_b * _v1[0];
603  }
604  }
605 }
606 
607 void
609 {
610  _soln_vector.zero();
611 
612  for (unsigned int j = 0; j < 4; ++j)
613  {
614  _soln_disp_index[j].resize(_ndisp);
615  _soln_rot_index[j].resize(_nrot);
616 
617  for (unsigned int i = 0; i < _ndisp; ++i)
618  {
619  _soln_disp_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _disp_num[i], 0);
620  _soln_vector(j + i * _nodes.size()) =
622  if (ADReal::do_derivatives)
624  _soln_vector(j + i * _nodes.size()).derivatives(), _soln_disp_index[j][i], 1.);
625  }
626 
627  for (unsigned int i = 0; i < _nrot; ++i)
628  {
629  _soln_rot_index[j][i] = _nodes[j]->dof_number(_nonlinear_sys.number(), _rot_num[i], 0);
630  _soln_vector(j + 12 + i * _nodes.size()) =
632  if (ADReal::do_derivatives)
634  _soln_vector(j + 12 + i * _nodes.size()).derivatives(), _soln_rot_index[j][i], 1.);
635  }
636  }
637 }
std::vector< ADMaterialProperty< RankTwoTensor > * > _strain_increment
Strain increment in the covariant coordinate system.
NonlinearSystemBase & _nonlinear_sys
Reference to the nonlinear system object.
RankTwoTensorTempl< T > inverse() const
std::vector< const MaterialProperty< DenseMatrix< Real > > * > _B_old
Old B_matrix for small strain.
const MaterialProperty< RealVectorValue > & _node_normal_old
Material property storing the old normal to the element at the 4 nodes.
FEProblemBase & _fe_problem
const QBase *const & _qrule
std::vector< unsigned int > _rot_num
Variable numbers corresponding to the rotational variables.
std::vector< ADRealVectorValue > _v2
First tangential vectors at nodes.
std::vector< std::vector< Real > > _dphidxi_map
Derivatives of shape functions w.r.t isoparametric coordinates xi.
const NumericVector< Number > & _sol_old
ADRealVectorValue _x2
simulation variables
auto norm() const -> decltype(std::norm(Real()))
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
virtual void zero() override final
std::vector< MaterialProperty< RankTwoTensor > * > _contravariant_transformation_matrix
Contravariant base vector matrix material property to transform strain.
std::vector< const MaterialProperty< Real > * > _J_map_old
Old material property containing jacobian of transformation.
std::vector< unsigned int > _disp_num
Variable numbers corresponding to the displacement variables.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_deta_old
Old derivative of global x, y and z w.r.t isoparametric coordinate eta.
unsigned int number() const
std::vector< std::vector< Real > > _dphideta_map
Derivatives of shape functions w.r.t isoparametric coordinates eta.
std::vector< ADMaterialProperty< DenseMatrix< Real > > * > _B
B_matrix for small strain.
unsigned int dim
static const std::string component
Definition: NS.h:153
registerMooseObject("SolidMechanicsApp", ADComputeIncrementalShellStrain)
std::vector< ADMaterialProperty< RankTwoTensor > * > _total_strain
Total strain increment in the covariant coordinate system.
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_dzeta
Derivative of global x, y and z w.r.t isoparametric coordinate zeta.
auto raw_value(const Eigen::Map< T > &in)
static constexpr std::size_t dim
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_dxi_old
Old derivative of global x, y and z w.r.t isoparametric coordinate xi.
MooseVariable * getVar(const std::string &var_name, unsigned int comp)
const std::vector< double > y
const VariableValue & _thickness
Coupled variable for the shell thickness.
const Number zero
unsigned int _nrot
Number of coupled rotational variables.
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual void computeLocalCoordinates()
Computes the element&#39;s local coordinates and store in _e1, _e2 and _e3.
std::vector< const MaterialProperty< RankTwoTensor > * > _covariant_transformation_matrix_old
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override
virtual void computeBMatrix()
Computes the B matrix that connects strains to nodal displacements and rotations. ...
static RankTwoTensorTempl initializeFromRows(const libMesh::TypeVector< ADReal > &row0, const libMesh::TypeVector< ADReal > &row1, const libMesh::TypeVector< ADReal > &row2)
unsigned int _qp
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_deta
Derivative of global x, y and z w.r.t isoparametric coordinate eta.
std::vector< std::vector< unsigned int > > _soln_rot_index
Indices of solution vector corresponding to rotation DOFs in 2 directions at the 4 nodes...
virtual void computeGMatrix()
Computes the transformation matrix from natural coordinates to local cartesian coordinates for elasti...
MooseMesh & _mesh
std::vector< const MaterialProperty< RankTwoTensor > * > _ge_old
Old ge matrix for elasticity tensor conversion.
std::vector< const MaterialProperty< RankTwoTensor > * > _contravariant_transformation_matrix_old
static InputParameters validParams()
const std::vector< double > x
THREAD_ID _tid
Real f(Real x)
Test function for Brents method.
std::vector< const Node * > _nodes
Vector storing pointers to the nodes of the shell element.
std::vector< Point > _t_points
Quadrature points in the out of plane direction in isoparametric coordinate system.
ADDenseVector _strain_vector
Vector that stores the strain in the the 2 axial and 3 shear directions.
unsigned int _ndisp
Number of coupled displacement variables.
std::vector< ADMaterialProperty< RealVectorValue > * > _dxyz_dxi
Derivative of global x, y and z w.r.t isoparametric coordinate xi.
ADDenseVector _soln_vector
Vector that stores the incremental solution at all the 20 DOFs in the 4 noded element.
ADMaterialProperty< RealVectorValue > & _node_normal
Material property storing the normal to the element at the 4 nodes. Stored as a material property for...
unsigned int number() const
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
std::vector< const MaterialProperty< RankTwoTensor > * > _local_transformation_matrix_old
bool hasSecondOrderElements()
std::unique_ptr< libMesh::QGauss > _t_qrule
Quadrature rule in the out of plane direction.
std::vector< ADRealVectorValue > _v1
First tangential vectors at nodes.
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
std::vector< std::vector< unsigned int > > _soln_disp_index
Indices of solution vector corresponding to displacement DOFs in 3 directions at the 4 nodes...
std::vector< MaterialProperty< RankTwoTensor > * > _local_transformation_matrix
Transformation matrix to map stresses from global coordinate to the element&#39;s local coordinate...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void updateGVectors()
Updates the vectors required for shear locking computation for finite rotations.
std::vector< ADMaterialProperty< RankTwoTensor > * > _ge
ge matrix for elasticity tensor conversion
std::vector< std::vector< Real > > _phi_map
Shape function value.
IntRange< T > make_range(T beg, T end)
std::vector< MaterialProperty< RankTwoTensor > * > _total_global_strain
Total strain in global coordinate system.
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
const NumericVector< Number > *const & _sol
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< MaterialProperty< RankTwoTensor > * > _covariant_transformation_matrix
Covariant base vector matrix material property to transform stress.
std::vector< const MaterialProperty< RealVectorValue > * > _dxyz_dzeta_old
Old derivative of global x, y and z w.r.t isoparametric coordinate zeta.
virtual void computeSolnVector()
Computes the 20x1 soln vector and its derivatives for each shell element.
std::vector< Point > _2d_points
Quadrature points in the in-plane direction in isoparametric coordinate system.
std::vector< ADMaterialProperty< Real > * > _J_map
Material property containing jacobian of transformation.
static const std::string k
Definition: NS.h:130
const Elem *const & _current_elem
ADComputeIncrementalShellStrain(const InputParameters &parameters)
ADMaterialProperty< RankTwoTensor > * _transformation_matrix
Rotation matrix material property.
virtual void computeNodeNormal()
Computes the node normal at each node.
const Real pi
std::vector< const MaterialProperty< RankTwoTensor > * > _total_strain_old
Old total strain increment in the covariant coordinate system.