www.mooseframework.org
ADComputeFiniteShellStrain.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "DenseMatrix.h"
12 
13 #include "libmesh/quadrature.h"
14 #include "libmesh/utility.h"
15 #include "libmesh/enum_quadrature_type.h"
16 #include "libmesh/fe_type.h"
17 #include "libmesh/string_to_enum.h"
18 #include "libmesh/quadrature_gauss.h"
19 
21 
24  params.addClassDescription("Compute a large strain increment for the shell."););
25 
26 template <ComputeStage compute_stage>
28  const InputParameters & parameters)
29  : ADComputeIncrementalShellStrain<compute_stage>(parameters), _B_nl()
30 {
31  _B_nl.resize(_t_points.size());
32 
33  for (unsigned int i = 0; i < _t_points.size(); ++i)
34  _B_nl[i] = &declareADProperty<DenseMatrix<Real>>("B_nl_t_points_" + std::to_string(i));
35 }
36 
37 template <ComputeStage compute_stage>
38 void
40 {
41  computeGMatrix();
42 
43  ADDenseMatrix b(5, 20);
44  for (unsigned int t = 0; t < _t_points.size(); ++t)
45  (*_B_nl[t])[_qp] = b;
46 }
47 
48 template <ComputeStage compute_stage>
49 void
51 {
52  // quadrature points in isoparametric space
53  _2d_points = _qrule->get_points(); // would be in 2D
54 
55  unsigned int dim = _current_elem->dim();
56 
57  // derivatives of shape functions (dphidxi, dphideta and dphidzeta) evaluated at quadrature points
58  // (in isoparametric space).
59  FEType fe_type(Utility::string_to_enum<Order>("First"),
60  Utility::string_to_enum<FEFamily>("LAGRANGE"));
61  auto & fe = _fe_problem.assembly(_tid).getFE(fe_type, dim);
62  _dphidxi_map = fe->get_fe_map().get_dphidxi_map();
63  _dphideta_map = fe->get_fe_map().get_dphideta_map();
64  _phi_map = fe->get_fe_map().get_phi_map();
65 
66  for (unsigned int i = 0; i < 4; ++i)
67  _nodes[i] = _current_elem->node_ptr(i);
68 
69  for (unsigned int i = 0; i < _2d_points.size(); ++i)
70  {
71  for (unsigned int j = 0; j < _t_points.size(); ++j)
72  {
73  (*_ge[j])[i] = (*_ge_old[j])[i];
74  (*_J_map[j])[i] = (*_J_map_old[j])[i];
75  }
76  }
77 
78  computeSolnVector();
79  updatedxyz();
80 
81  for (unsigned int k = 0; k < _nodes.size(); ++k)
82  _node_normal[k] = _node_normal_old[k];
83 
84  computeBMatrix();
85 
86  computeNodeNormal();
87 
88  computeBNLMatrix();
89 
90  for (unsigned int i = 0; i < _2d_points.size(); ++i)
91  {
92  for (unsigned int j = 0; j < _t_points.size(); ++j)
93  {
94  // compute strain increment in covariant coordinate system using B and _soln_vector
95  for (unsigned int temp1 = 0; temp1 < 5; ++temp1)
96  {
97  _strain_vector(temp1) = 0.0;
98  for (unsigned int temp2 = 0; temp2 < 20; ++temp2)
99  _strain_vector(temp1) += (*_B[j])[i](temp1, temp2) * _soln_vector(temp2);
100  }
101 
102  (*_strain_increment[j])[i](0, 0) = _strain_vector(0);
103  (*_strain_increment[j])[i](1, 1) = _strain_vector(1);
104  (*_strain_increment[j])[i](0, 1) = _strain_vector(2);
105  (*_strain_increment[j])[i](0, 2) = _strain_vector(3);
106  (*_strain_increment[j])[i](1, 2) = _strain_vector(4);
107  (*_strain_increment[j])[i](1, 0) = (*_strain_increment[j])[i](0, 1);
108  (*_strain_increment[j])[i](2, 0) = (*_strain_increment[j])[i](0, 2);
109  (*_strain_increment[j])[i](2, 1) = (*_strain_increment[j])[i](1, 2);
110  (*_total_strain[j])[i] = (*_total_strain_old[j])[i] + (*_strain_increment[j])[i];
111  }
112  }
113  copyDualNumbersToValues();
114 }
115 
116 template <ComputeStage compute_stage>
117 void
119 {
120  // update _node_normal
121  for (unsigned int k = 0; k < _nodes.size(); ++k)
122  {
123  _node_normal[k] =
124  -_v2[k] * _soln_vector(12 + k) + _v1[k] * _soln_vector(16 + k) + _node_normal_old[k];
125  _node_normal[k] /= _node_normal[k].norm();
126  }
127 }
128 
129 template <ComputeStage compute_stage>
130 void
132 {
133  for (unsigned int i = 0; i < _2d_points.size(); ++i)
134  {
135  for (unsigned int j = 0; j < _t_points.size(); ++j)
136  {
137  for (unsigned int component = 0; component < 3; ++component)
138  {
139  (*_dxyz_dxi[j])[i](component) = 0.0;
140  (*_dxyz_deta[j])[i](component) = 0.0;
141  (*_dxyz_dzeta[j])[i](component) = 0.0;
142  for (unsigned int k = 0; k < _nodes.size(); ++k)
143  {
144  (*_dxyz_dxi[j])[i](component) +=
145  _dphidxi_map[k][i] *
146  ((*_nodes[k])(component) + _sol_old(_soln_disp_index[k][component])) +
147  _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] *
148  _node_normal_old[k](component);
149  (*_dxyz_deta[j])[i](component) +=
150  _dphideta_map[k][i] *
151  ((*_nodes[k])(component) + _sol_old(_soln_disp_index[k][component])) +
152  _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] *
153  _node_normal_old[k](component);
154  (*_dxyz_dzeta[j])[i](component) +=
155  _thickness[i] * _phi_map[k][i] * _node_normal_old[k](component) / 2.0;
156  }
157  }
158  }
159  }
160 }
161 
162 template <ComputeStage compute_stage>
163 void
165 {
166  for (unsigned int component = 0; component < 3; ++component)
167  {
168  _g1_a(component) +=
169  0.5 * (_sol_old(_soln_disp_index[2][component]) - _sol_old(_soln_disp_index[3][component]));
170  _g1_c(component) +=
171  0.5 * (_sol_old(_soln_disp_index[1][component]) - _sol_old(_soln_disp_index[0][component]));
172  _g2_b(component) +=
173  0.5 * (_sol_old(_soln_disp_index[3][component]) - _sol_old(_soln_disp_index[0][component]));
174  _g2_d(component) +=
175  0.5 * (_sol_old(_soln_disp_index[2][component]) - _sol_old(_soln_disp_index[1][component]));
176  }
177 }
178 
179 template <ComputeStage compute_stage>
180 void
182 {
183  // compute BNL matrix - rows correspond to [ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, uz1, uz2, uz3,
184  // uz4, a1, a2, a3, a4, b1, b2, b3, b4]
185 
186  for (unsigned int i = 0; i < _2d_points.size(); ++i)
187  {
188  for (unsigned int j = 0; j < _t_points.size(); ++j)
189  {
190  (*_B_nl[j])[i].resize(5, 20);
191  (*_B_nl[j])[i].zero();
192  for (unsigned int k = 0; k < 4; ++k)
193  {
194  for (unsigned int p = 0; p < 4; ++p) // loop over nodes
195  {
196  // corresponding to strain(0,0)
197  (*_B_nl[j])[i](0, k) += _dphidxi_map[k][i] * _dphidxi_map[p][i] *
198  (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
199  (-_soln_vector(p + 12) * _v2[p](0) +
200  _soln_vector(p + 16) * _v1[p](0)));
201  (*_B_nl[j])[i](0, 4 + k) +=
202  _dphidxi_map[k][i] * _dphidxi_map[p][i] *
203  (_soln_vector(p + 4) +
204  _t_points[j](0) / 2.0 * _thickness[i] *
205  (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
206  (*_B_nl[j])[i](0, 8 + k) +=
207  _dphidxi_map[k][i] * _dphidxi_map[p][i] *
208  (_soln_vector(p + 8) +
209  _t_points[j](0) / 2.0 * _thickness[i] *
210  (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
211  (*_B_nl[j])[i](0, 12 + k) +=
212  _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] * _dphidxi_map[p][i] *
213  (-(_v2[p](0) * _soln_vector(p) + _v2[p](1) * _soln_vector(p + 4) +
214  _v2[p](2) * _soln_vector(p + 8)) +
215  _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
216  (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
217  (*_B_nl[j])[i](0, 16 + k) +=
218  _t_points[j](0) / 2.0 * _thickness[i] * _dphidxi_map[k][i] * _dphidxi_map[p][i] *
219  ((_v1[p](0) * _soln_vector(p) + _v1[p](1) * _soln_vector(p + 4) +
220  _v1[p](2) * _soln_vector(p + 8)) +
221  _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
222  (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
223 
224  // corresponding to strain(1,1)
225  (*_B_nl[j])[i](1, k) += _dphideta_map[k][i] * _dphideta_map[p][i] *
226  (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
227  (-_soln_vector(p + 12) * _v2[p](0) +
228  _soln_vector(p + 16) * _v1[p](0)));
229  (*_B_nl[j])[i](1, 4 + k) +=
230  _dphideta_map[k][i] * _dphideta_map[p][i] *
231  (_soln_vector(p + 4) +
232  _t_points[j](0) / 2.0 * _thickness[i] *
233  (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
234  (*_B_nl[j])[i](1, 8 + k) +=
235  _dphideta_map[k][i] * _dphideta_map[p][i] *
236  (_soln_vector(p + 8) +
237  _t_points[j](0) / 2.0 * _thickness[i] *
238  (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
239  (*_B_nl[j])[i](1, 12 + k) +=
240  _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] * _dphideta_map[p][i] *
241  (-(_v2[p](0) * _soln_vector(p) + _v2[p](1) * _soln_vector(p + 4) +
242  _v2[p](2) * _soln_vector(p + 8)) +
243  _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
244  (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
245  (*_B_nl[j])[i](1, 16 + k) +=
246  _t_points[j](0) / 2.0 * _thickness[i] * _dphideta_map[k][i] * _dphideta_map[p][i] *
247  ((_v1[p](0) * _soln_vector(p) + _v1[p](1) * _soln_vector(p + 4) +
248  _v1[p](2) * _soln_vector(p + 8)) +
249  _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
250  (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
251 
252  // terms corresponding to strain(2,2) are 0.
253 
254  // corresponding to strain(0,1)
255  (*_B_nl[j])[i](2, k) += 0.5 *
256  (_dphidxi_map[k][i] * _dphideta_map[p][i] +
257  _dphideta_map[k][i] * _dphidxi_map[p][i]) *
258  (_soln_vector(p) + _t_points[j](0) / 2.0 * _thickness[i] *
259  (-_soln_vector(p + 12) * _v2[p](0) +
260  _soln_vector(p + 16) * _v1[p](0)));
261  (*_B_nl[j])[i](2, 4 + k) +=
262  0.5 *
263  (_dphidxi_map[k][i] * _dphideta_map[p][i] +
264  _dphideta_map[k][i] * _dphidxi_map[p][i]) *
265  (_soln_vector(p + 4) +
266  _t_points[j](0) / 2.0 * _thickness[i] *
267  (-_soln_vector(p + 12) * _v2[p](1) + _soln_vector(p + 16) * _v1[p](1)));
268  (*_B_nl[j])[i](2, 8 + k) +=
269  0.5 *
270  (_dphidxi_map[k][i] * _dphideta_map[p][i] +
271  _dphideta_map[k][i] * _dphidxi_map[p][i]) *
272  (_soln_vector(p + 8) +
273  _t_points[j](0) / 2.0 * _thickness[i] *
274  (-_soln_vector(p + 12) * _v2[p](2) + _soln_vector(p + 16) * _v1[p](2)));
275  (*_B_nl[j])[i](2, 12 + k) +=
276  _t_points[j](0) * 0.25 *
277  (_dphidxi_map[k][i] * _dphideta_map[p][i] +
278  _dphideta_map[k][i] * _dphidxi_map[p][i]) *
279  _thickness[i] *
280  (-(_v2[k](0) * _soln_vector(p) + _v2[k](1) * _soln_vector(p + 4) +
281  _v2[k](2) * _soln_vector(p + 8)) +
282  _t_points[j](0) / 2.0 * _thickness[i] * _v2[k] *
283  (_v2[p] * _soln_vector(p + 12) - _v1[p] * _soln_vector(p + 16)));
284  (*_B_nl[j])[i](2, 16 + k) +=
285  _t_points[j](0) * 0.25 *
286  (_dphidxi_map[k][i] * _dphideta_map[p][i] +
287  _dphideta_map[k][i] * _dphidxi_map[p][i]) *
288  _thickness[i] *
289  ((_v1[k](0) * _soln_vector(p) + _v1[k](1) * _soln_vector(p + 4) +
290  _v1[k](2) * _soln_vector(p + 8)) +
291  _t_points[j](0) / 2.0 * _thickness[i] * _v1[k] *
292  (-_v2[p] * _soln_vector(p + 12) + _v1[p] * _soln_vector(p + 16)));
293  }
294  }
295 
296  for (unsigned int component = 0; component < 3; ++component)
297  {
298  // corresponding to strain(0,2)
299  (*_B_nl[j])[i](3, 2 + component * 4) +=
300  1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] *
301  (-_soln_vector(12 + 2) * _v2[2](component) + _soln_vector(16 + 2) * _v1[2](component) -
302  _soln_vector(12 + 3) * _v2[3](component) + _soln_vector(16 + 3) * _v1[3](component));
303  (*_B_nl[j])[i](3, 3 + component * 4) += -(*_B_nl[j])[i](3, 2 + component * 4);
304 
305  (*_B_nl[j])[i](3, 1 + component * 4) +=
306  1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] *
307  (-_soln_vector(12 + 1) * _v2[1](component) + _soln_vector(16 + 1) * _v1[1](component) -
308  _soln_vector(12 + 0) * _v2[0](component) + _soln_vector(16 + 0) * _v1[0](component));
309  (*_B_nl[j])[i](3, component * 4) += -(*_B_nl[j])[i](3, 1 + component * 4);
310 
311  // adding contributions corresponding to alpha 2 and 3 and beta 2 and 3
312  (*_B_nl[j])[i](3, 12 + 2) +=
313  -1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v2[2](component) *
314  (_soln_vector(2 + component * 4) - _soln_vector(3 + component * 4));
315  (*_B_nl[j])[i](3, 16 + 2) +=
316  1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v1[2](component) *
317  (_soln_vector(2 + component * 4) - _soln_vector(3 + component * 4));
318  (*_B_nl[j])[i](3, 12 + 3) +=
319  -1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v2[3](component) *
320  (_soln_vector(2 + component * 4) - _soln_vector(3 + component * 4));
321  (*_B_nl[j])[i](3, 16 + 3) +=
322  1.0 / 32.0 * (1.0 + _2d_points[i](1)) * _thickness[i] * _v1[3](component) *
323  (_soln_vector(2 + component * 4) - _soln_vector(3 + component * 4));
324 
325  // adding contributions corresponding to alpha 1 and 0 and beta 1 and 0
326  (*_B_nl[j])[i](3, 12 + 1) +=
327  -1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v2[1](component) *
328  (_soln_vector(1 + component * 4) - _soln_vector(component * 4));
329  (*_B_nl[j])[i](3, 16 + 1) +=
330  1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v1[1](component) *
331  (_soln_vector(1 + component * 4) - _soln_vector(component * 4));
332  (*_B_nl[j])[i](3, 12 + 0) +=
333  -1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v2[0](component) *
334  (_soln_vector(1 + component * 4) - _soln_vector(component * 4));
335  (*_B_nl[j])[i](3, 16 + 0) +=
336  1.0 / 32.0 * (1.0 - _2d_points[i](1)) * _thickness[i] * _v1[0](component) *
337  (_soln_vector(1 + component * 4) - _soln_vector(component * 4));
338 
339  // corresponding to strain(1,2)
340  (*_B_nl[j])[i](4, 2 + component * 4) +=
341  1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] *
342  (-_soln_vector(12 + 2) * _v2[2](component) + _soln_vector(16 + 2) * _v1[2](component) -
343  _soln_vector(12 + 1) * _v2[1](component) + _soln_vector(16 + 1) * _v1[1](component));
344  (*_B_nl[j])[i](4, 1 + component * 4) += -(*_B_nl[j])[i](3, 2 + component * 4);
345 
346  (*_B_nl[j])[i](4, 3 + component * 4) +=
347  1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] *
348  (-_soln_vector(12 + 3) * _v2[3](component) + _soln_vector(16 + 3) * _v1[3](component) -
349  _soln_vector(12 + 0) * _v2[0](component) + _soln_vector(16 + 0) * _v1[0](component));
350  (*_B_nl[j])[i](4, component * 4) += -(*_B_nl[j])[i](3, 3 + component * 4);
351 
352  // adding contributions corresponding to alpha 2, 1 and beta 2 , 1
353  (*_B_nl[j])[i](4, 12 + 2) +=
354  -1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v2[2](component) *
355  (_soln_vector(2 + component * 4) - _soln_vector(1 + component * 4));
356  (*_B_nl[j])[i](4, 16 + 2) +=
357  1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v1[2](component) *
358  (_soln_vector(2 + component * 4) - _soln_vector(1 + component * 4));
359  (*_B_nl[j])[i](4, 12 + 1) +=
360  -1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v2[1](component) *
361  (_soln_vector(2 + component * 4) - _soln_vector(1 + component * 4));
362  (*_B_nl[j])[i](4, 16 + 1) +=
363  1.0 / 32.0 * (1.0 + _2d_points[i](0)) * _thickness[i] * _v1[1](component) *
364  (_soln_vector(2 + component * 4) - _soln_vector(1 + component * 4));
365 
366  // adding contributions corresponding to alpha 3, 0 and beta 3 , 0
367  (*_B_nl[j])[i](4, 12 + 3) +=
368  -1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v2[3](component) *
369  (_soln_vector(3 + component * 4) - _soln_vector(component * 4));
370  (*_B_nl[j])[i](4, 16 + 3) +=
371  1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v1[3](component) *
372  (_soln_vector(3 + component * 4) - _soln_vector(component * 4));
373  (*_B_nl[j])[i](4, 12 + 0) +=
374  -1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v2[0](component) *
375  (_soln_vector(3 + component * 4) - _soln_vector(component * 4));
376  (*_B_nl[j])[i](4, 16 + 0) +=
377  1.0 / 32.0 * (1.0 - _2d_points[i](0)) * _thickness[i] * _v1[0](component) *
378  (_soln_vector(3 + component * 4) - _soln_vector(component * 4));
379  }
380  }
381  }
382 }
ADComputeFiniteShellStrain::updatedxyz
virtual void updatedxyz() override
Updates covariant vectors at each qp for finite rotations.
Definition: ADComputeFiniteShellStrain.C:131
ADComputeFiniteShellStrain::_B_nl
std::vector< ADMaterialProperty(DenseMatrix< Real >) * > _B_nl
Material property to store the B_nl matrix at each quadrature point.
Definition: ADComputeFiniteShellStrain.h:45
defineADValidParams
defineADValidParams(ADComputeFiniteShellStrain, ADComputeIncrementalShellStrain, params.addClassDescription("Compute a large strain increment for the shell.");)
ADComputeFiniteShellStrain::initQpStatefulProperties
virtual void initQpStatefulProperties() override
Definition: ADComputeFiniteShellStrain.C:39
ADComputeFiniteShellStrain::computeProperties
virtual void computeProperties() override
Definition: ADComputeFiniteShellStrain.C:50
MaterialTensorCalculatorTools::component
Real component(const SymmTensor &symm_tensor, unsigned int index)
Definition: MaterialTensorCalculatorTools.C:16
ADComputeIncrementalShellStrain::_t_points
std::vector< Point > _t_points
Quadrature points in the out of plane direction in isoparametric coordinate system.
Definition: ADComputeIncrementalShellStrain.h:153
ADComputeFiniteShellStrain::computeBNLMatrix
virtual void computeBNLMatrix()
Computes the B_nl matrix that connects the nonlinear strains to the nodal displacements and rotations...
Definition: ADComputeFiniteShellStrain.C:181
ADComputeFiniteShellStrain::ADComputeFiniteShellStrain
ADComputeFiniteShellStrain(const InputParameters &parameters)
Definition: ADComputeFiniteShellStrain.C:27
ADComputeFiniteShellStrain.h
ADComputeFiniteShellStrain::computeNodeNormal
virtual void computeNodeNormal() override
Computes the node normal at each node.
Definition: ADComputeFiniteShellStrain.C:118
registerADMooseObject
registerADMooseObject("TensorMechanicsApp", ADComputeFiniteShellStrain)
ADComputeFiniteShellStrain::updateGVectors
virtual void updateGVectors() override
Updates the vectors required for shear locking computation for finite rotations.
Definition: ADComputeFiniteShellStrain.C:164
ADComputeFiniteShellStrain
ADComputeFiniteShellStrain computes the strain increment term for shell elements under finite displac...
Definition: ADComputeFiniteShellStrain.h:19
ADComputeIncrementalShellStrain
Definition: ADComputeIncrementalShellStrain.h:56