MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
solid_element_3d.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // MAST includes
21 #include "elasticity/solid_element_3d.h"
22 #include "elasticity/stress_output_base.h"
23 #include "base/boundary_condition_base.h"
24 #include "base/system_initialization.h"
25 #include "base/nonlinear_system.h"
26 #include "base/assembly_base.h"
27 #include "base/field_function_base.h"
28 #include "numerics/fem_operator_matrix.h"
29 #include "mesh/geom_elem.h"
30 #include "mesh/fe_base.h"
31 #include "property_cards/element_property_card_base.h"
32 
33 
34 MAST::StructuralElement3D::
35 StructuralElement3D(MAST::SystemInitialization& sys,
36  const MAST::GeomElem& elem,
37  const MAST::ElementPropertyCardBase& p):
38 MAST::StructuralElementBase(sys, elem, p) {
39 
40 }
41 
42 
43 MAST::StructuralElement3D::~StructuralElement3D() {
44 
45 }
46 
47 
48 bool
49 MAST::StructuralElement3D::inertial_residual (bool request_jacobian,
50  RealVectorX& f,
51  RealMatrixX& jac_xddot,
52  RealMatrixX& jac_xdot,
53  RealMatrixX& jac) {
54 
55  std::unique_ptr<MAST::FEBase>
56  fe(_elem.init_fe(true, false));
57 
58  const std::vector<Real>& JxW = fe->get_JxW();
59  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
60  const std::vector<std::vector<Real> >& phi = fe->get_phi();
61 
62  const unsigned int
63  n_phi = (unsigned int)phi.size(),
64  n1 =3,
65  n2 =3*n_phi;
66 
67  RealMatrixX
68  material_mat,
69  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
70  mat2_n2n2 = RealMatrixX::Zero(n2, n2);
71  RealVectorX
72  phi_vec = RealVectorX::Zero(n_phi),
73  vec1_n1 = RealVectorX::Zero(n1),
74  vec2_n2 = RealVectorX::Zero(n2),
75  local_acc = RealVectorX::Zero(n2);
76 
77  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
78  mat_inertia = _property.inertia_matrix(*this);
79 
80  MAST::FEMOperatorMatrix Bmat;
81 
82  local_acc.topRows(n2) = _local_accel.topRows(n2);
83 
84  if (_property.if_diagonal_mass_matrix()) {
85 
86  (*mat_inertia)(xyz[0], _time, material_mat);
87 
88  Real vol = 0.;
89  const unsigned int nshp = fe->n_shape_functions();
90  for (unsigned int i=0; i<JxW.size(); i++)
91  vol += JxW[i];
92  vol /= (1.* nshp);
93  for (unsigned int i_var=0; i_var<3; i_var++)
94  for (unsigned int i=0; i<nshp; i++)
95  jac_xddot(i_var*nshp+i, i_var*nshp+i) =
96  vol*material_mat(i_var, i_var);
97 
98  f.topRows(n2) = jac_xddot.topLeftCorner(n2, n2) * _local_accel.topRows(n2);
99  }
100  else {
101  libMesh::Point p;
102 
103  for (unsigned int qp=0; qp<JxW.size(); qp++) {
104 
105  (*mat_inertia)(xyz[0], _time, material_mat);
106 
107  // now set the shape function values
108  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
109  phi_vec(i_nd) = phi[i_nd][qp];
110 
111  Bmat.reinit(3, phi_vec);
112 
113  Bmat.left_multiply(mat1_n1n2, material_mat);
114 
115  vec1_n1 = mat1_n1n2 * local_acc;
116  Bmat.vector_mult_transpose(vec2_n2, vec1_n1);
117 
118  f.topRows(n2) += JxW[qp] * vec2_n2;
119 
120  if (request_jacobian) {
121 
122  Bmat.right_multiply_transpose(mat2_n2n2,
123  mat1_n1n2);
124  jac_xddot.topLeftCorner(n2, n2) += JxW[qp]*mat2_n2n2;
125  }
126  }
127  }
128 
129  return request_jacobian;
130 }
131 
132 
133 
134 
135 bool
136 MAST::StructuralElement3D::internal_residual(bool request_jacobian,
137  RealVectorX& f,
138  RealMatrixX& jac) {
139 
140  std::unique_ptr<MAST::FEBase>
141  fe(_elem.init_fe(true, false));
142 
143  const std::vector<Real>& JxW = fe->get_JxW();
144  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
145  const unsigned int
146  n_phi = (unsigned int)fe->n_shape_functions(),
147  n1 =6,
148  n2 =3*n_phi,
149  n3 =30,
150  n_nodes =_elem.get_reference_elem().n_nodes();
151 
152  RealMatrixX
153  material_mat,
154  mat_x = RealMatrixX::Zero(6,3),
155  mat_y = RealMatrixX::Zero(6,3),
156  mat_z = RealMatrixX::Zero(6,3),
157  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
158  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
159  mat3_3n2 = RealMatrixX::Zero(3, n2),
160  mat4_33 = RealMatrixX::Zero(3, 3),
161  mat5_n1n3 = RealMatrixX::Zero(n1, n3),
162  mat6_n2n3 = RealMatrixX::Zero(n2, n3),
163  mat7_3n3 = RealMatrixX::Zero(3, n3),
164  Gmat = RealMatrixX::Zero(6, n3),
165  K_alphaalpha = RealMatrixX::Zero(n3, n3),
166  K_ualpha = RealMatrixX::Zero(n2, n3),
167  K_corr = RealMatrixX::Zero(n2, n2);
168  RealVectorX
169  strain = RealVectorX::Zero(6),
170  stress = RealVectorX::Zero(6),
171  vec1_n1 = RealVectorX::Zero(n1),
172  vec2_n2 = RealVectorX::Zero(n2),
173  vec3_3 = RealVectorX::Zero(3),
174  local_disp= RealVectorX::Zero(n2),
175  f_alpha = RealVectorX::Zero(n3),
176  alpha = RealVectorX::Zero(n3);//*_incompatible_sol;
177 
178  // copy the values from the global to the local element
179  local_disp.topRows(n2) = _local_sol.topRows(n2);
180 
181  std::unique_ptr<MAST::FieldFunction<RealMatrixX> > mat_stiff =
182  _property.stiffness_A_matrix(*this);
183 
184  MAST::FEMOperatorMatrix
185  Bmat_lin,
186  Bmat_nl_x,
187  Bmat_nl_y,
188  Bmat_nl_z,
189  Bmat_nl_u,
190  Bmat_nl_v,
191  Bmat_nl_w,
192  Bmat_inc;
193  // six stress components, related to three displacements
194  Bmat_lin.reinit(n1, 3, n_nodes);
195  Bmat_nl_x.reinit(3, 3, n_nodes);
196  Bmat_nl_y.reinit(3, 3, n_nodes);
197  Bmat_nl_z.reinit(3, 3, n_nodes);
198  Bmat_nl_u.reinit(3, 3, n_nodes);
199  Bmat_nl_v.reinit(3, 3, n_nodes);
200  Bmat_nl_w.reinit(3, 3, n_nodes);
201  Bmat_inc.reinit(n1, n3, 1); // six stress-strain components
202 
203  /*
204  // initialize the incompatible mode mapping at element mid-point
205  _init_incompatible_fe_mapping(_elem);
206 
208  // first for loop to evaluate alpha
209  for (unsigned int qp=0; qp<JxW.size(); qp++) {
210 
211  // get the material matrix
212  (*mat_stiff)(xyz[qp], _time, material_mat);
213 
214  this->initialize_green_lagrange_strain_operator(qp,
215  *_fe,
216  local_disp,
217  strain,
218  mat_x, mat_y, mat_z,
219  Bmat_lin,
220  Bmat_nl_x,
221  Bmat_nl_y,
222  Bmat_nl_z,
223  Bmat_nl_u,
224  Bmat_nl_v,
225  Bmat_nl_w);
226  this->initialize_incompatible_strain_operator(qp, *_fe, Bmat_inc, Gmat);
227 
228  // calculate the incompatible mode matrices
229  // incompatible mode diagonal stiffness matrix
230  mat5_n1n3 = material_mat * Gmat;
231  K_alphaalpha += JxW[qp] * ( Gmat.transpose() * mat5_n1n3);
232 
233  // off-diagonal coupling matrix
234  // linear strain term
235  Bmat_lin.right_multiply_transpose(mat6_n2n3, mat5_n1n3);
236  K_ualpha += JxW[qp] * mat6_n2n3;
237 
238  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
239 
240  // nonlinear component
241  // along x
242  mat7_3n3 = mat_x.transpose() * mat5_n1n3;
243  Bmat_nl_x.right_multiply_transpose(mat6_n2n3, mat7_3n3);
244  K_ualpha += JxW[qp] * mat6_n2n3;
245 
246  // along y
247  mat7_3n3 = mat_y.transpose() * mat5_n1n3;
248  Bmat_nl_y.right_multiply_transpose(mat6_n2n3, mat7_3n3);
249  K_ualpha += JxW[qp] * mat6_n2n3;
250 
251  // along z
252  mat7_3n3 = mat_z.transpose() * mat5_n1n3;
253  Bmat_nl_z.right_multiply_transpose(mat6_n2n3, mat7_3n3);
254  K_ualpha += JxW[qp] * mat6_n2n3;
255  }
256  }
257 
258 
259  // incompatible mode corrections
260  K_alphaalpha = K_alphaalpha.inverse();
261  K_corr = K_ualpha * K_alphaalpha * K_ualpha.transpose();
262 
263  if (request_jacobian)
264  jac.topLeftCorner(n2, n2) -= K_corr;
265  */
266 
268  // second for loop to calculate the residual and stiffness contributions
269  for (unsigned int qp=0; qp<JxW.size(); qp++) {
270 
271  // get the material matrix
272  (*mat_stiff)(xyz[qp], _time, material_mat);
273 
274  this->initialize_green_lagrange_strain_operator(qp,
275  *fe,
276  local_disp,
277  strain,
278  mat_x, mat_y, mat_z,
279  Bmat_lin,
280  Bmat_nl_x,
281  Bmat_nl_y,
282  Bmat_nl_z,
283  Bmat_nl_u,
284  Bmat_nl_v,
285  Bmat_nl_w);
286  //this->initialize_incompatible_strain_operator(qp, *_fe, Bmat_inc, Gmat);
287 
288  // calculate the stress
289  stress = material_mat * (strain + Gmat * alpha);
290 
291  // residual from incompatible modes
292  f_alpha += JxW[qp] * Gmat.transpose() * stress;
293 
294  // calculate contribution to the residual
295  // linear strain operator
296  Bmat_lin.vector_mult_transpose(vec2_n2, stress);
297  f.topRows(n2) += JxW[qp] * vec2_n2;
298 
299  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
300 
301  // nonlinear strain operator
302  // x
303  vec3_3 = mat_x.transpose() * stress;
304  Bmat_nl_x.vector_mult_transpose(vec2_n2, vec3_3);
305  f.topRows(n2) += JxW[qp] * vec2_n2;
306 
307  // y
308  vec3_3 = mat_y.transpose() * stress;
309  Bmat_nl_y.vector_mult_transpose(vec2_n2, vec3_3);
310  f.topRows(n2) += JxW[qp] * vec2_n2;
311 
312  // z
313  vec3_3 = mat_z.transpose() * stress;
314  Bmat_nl_z.vector_mult_transpose(vec2_n2, vec3_3);
315  f.topRows(n2) += JxW[qp] * vec2_n2;
316  }
317 
318  if (request_jacobian) {
319 
320  // the strain includes the following expansion
321  // delta_epsilon = B_lin + mat_x B_x + mat_y B_y + mat_z B_z
322  // Hence, the tangent stiffness matrix will include
323  // components from epsilon^T C epsilon
324 
326  // B_lin^T C B_lin
327  Bmat_lin.left_multiply(mat1_n1n2, material_mat);
328  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
329  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
330 
331  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
332 
333  // B_x^T mat_x^T C B_lin
334  mat3_3n2 = mat_x.transpose() * mat1_n1n2;
335  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2);
336  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
337 
338  // B_y^T mat_y^T C B_lin
339  mat3_3n2 = mat_y.transpose() * mat1_n1n2;
340  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2);
341  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
342 
343  // B_z^T mat_z^T C B_lin
344  mat3_3n2 = mat_z.transpose() * mat1_n1n2;
345  Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2);
346  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
347 
349  for (unsigned int i_dim=0; i_dim<3; i_dim++) {
350  switch (i_dim) {
351  case 0:
352  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
353  break;
354 
355  case 1:
356  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
357  break;
358 
359  case 2:
360  Bmat_nl_z.left_multiply(mat1_n1n2, mat_z);
361  break;
362  }
363 
364  // B_lin^T C mat_x_i B_x_i
365  mat1_n1n2 = material_mat * mat1_n1n2;
366  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
367  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
368 
369  // B_x^T mat_x^T C mat_x B_x
370  mat3_3n2 = mat_x.transpose() * mat1_n1n2;
371  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2);
372  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
373 
374  // B_y^T mat_y^T C mat_x B_x
375  mat3_3n2 = mat_y.transpose() * mat1_n1n2;
376  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2);
377  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
378 
379  // B_z^T mat_z^T C mat_x B_x
380  mat3_3n2 = mat_z.transpose() * mat1_n1n2;
381  Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2);
382  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
383 
384  }
385 
386  // use the stress to calculate the final contribution
387  // to the Jacobian stiffness matrix
388  mat4_33(0,0) = stress(0);
389  mat4_33(1,1) = stress(1);
390  mat4_33(2,2) = stress(2);
391  mat4_33(0,1) = mat4_33(1,0) = stress(3);
392  mat4_33(1,2) = mat4_33(2,1) = stress(4);
393  mat4_33(0,2) = mat4_33(2,0) = stress(5);
394 
395  // u-disp
396  Bmat_nl_u.left_multiply(mat3_3n2, mat4_33);
397  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3_3n2);
398  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
399 
400  // v-disp
401  Bmat_nl_v.left_multiply(mat3_3n2, mat4_33);
402  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3_3n2);
403  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
404 
405  // w-disp
406  Bmat_nl_w.left_multiply(mat3_3n2, mat4_33);
407  Bmat_nl_w.right_multiply_transpose(mat2_n2n2, mat3_3n2);
408  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
409  }
410  }
411  }
412 
413  // if jacobian is requested, add a small diagonal value for the
414  // rotational dofs
415  if (request_jacobian)
416  jac.bottomRightCorner(n2, n2) += RealMatrixX::Identity(n2, n2) *
417  1.0e-20 * jac.diagonal().maxCoeff();
418 
419  // correction to the residual from incompatible mode
420  //f.topRows(n2) -= c * (K_alphaalpha * f_alpha);
421 
422  return request_jacobian;
423 }
424 
425 
426 
427 
428 void
429 MAST::StructuralElement3D::
430 update_incompatible_mode_solution(const RealVectorX& dsol) {
431 
432  std::unique_ptr<MAST::FEBase>
433  fe(_elem.init_fe(true, false));
434 
435  const std::vector<Real>& JxW = fe->get_JxW();
436  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
437  const unsigned int
438  n_phi = (unsigned int)fe->n_shape_functions(),
439  n1 =6,
440  n2 =3*n_phi,
441  n3 =30,
442  n_nodes = _elem.get_reference_elem().n_nodes();
443 
444  RealMatrixX
445  material_mat,
446  mat_x = RealMatrixX::Zero(6,3),
447  mat_y = RealMatrixX::Zero(6,3),
448  mat_z = RealMatrixX::Zero(6,3),
449  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
450  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
451  mat3_3n2 = RealMatrixX::Zero(3, n2),
452  mat4_33 = RealMatrixX::Zero(3, 3),
453  mat5_n1n3 = RealMatrixX::Zero(n1, n3),
454  mat6_n2n3 = RealMatrixX::Zero(n2, n3),
455  mat7_3n3 = RealMatrixX::Zero(3, n3),
456  Gmat = RealMatrixX::Zero(6, n3),
457  K_alphaalpha = RealMatrixX::Zero(n3, n3),
458  K_ualpha = RealMatrixX::Zero(n2, n3),
459  K_corr = RealMatrixX::Zero(n2, n2);
460  RealVectorX
461  strain = RealVectorX::Zero(6),
462  stress = RealVectorX::Zero(6),
463  vec1_n1 = RealVectorX::Zero(n1),
464  vec2_n2 = RealVectorX::Zero(n2),
465  vec3_3 = RealVectorX::Zero(3),
466  local_disp= RealVectorX::Zero(n2),
467  f = RealVectorX::Zero(n3),
468  alpha = RealVectorX::Zero(n3);//*_incompatible_sol;
469 
470  // copy the values from the global to the local element
471  local_disp.topRows(n2) = _local_sol.topRows(n2);
472 
473  std::unique_ptr<MAST::FieldFunction<RealMatrixX> > mat_stiff =
474  _property.stiffness_A_matrix(*this);
475 
476  MAST::FEMOperatorMatrix
477  Bmat_lin,
478  Bmat_nl_x,
479  Bmat_nl_y,
480  Bmat_nl_z,
481  Bmat_nl_u,
482  Bmat_nl_v,
483  Bmat_nl_w,
484  Bmat_inc;
485  // six stress components, related to three displacements
486  Bmat_lin.reinit(n1, 3, n_nodes);
487  Bmat_nl_x.reinit(3, 3, n_nodes);
488  Bmat_nl_y.reinit(3, 3, n_nodes);
489  Bmat_nl_z.reinit(3, 3, n_nodes);
490  Bmat_nl_u.reinit(3, 3, n_nodes);
491  Bmat_nl_v.reinit(3, 3, n_nodes);
492  Bmat_nl_w.reinit(3, 3, n_nodes);
493  Bmat_inc.reinit(n1, n3, 1); // six stress-strain components
494 
495 
496  // initialize the incompatible mode mapping at element mid-point
497  _init_incompatible_fe_mapping(_elem.get_reference_elem());
498 
500  // first for loop to evaluate alpha
501  for (unsigned int qp=0; qp<JxW.size(); qp++) {
502 
503  // get the material matrix
504  (*mat_stiff)(xyz[qp], _time, material_mat);
505 
506  this->initialize_green_lagrange_strain_operator(qp,
507  *fe,
508  local_disp,
509  strain,
510  mat_x, mat_y, mat_z,
511  Bmat_lin,
512  Bmat_nl_x,
513  Bmat_nl_y,
514  Bmat_nl_z,
515  Bmat_nl_u,
516  Bmat_nl_v,
517  Bmat_nl_w);
518  //this->initialize_incompatible_strain_operator(qp, *_fe, Bmat_inc, Gmat);
519 
520  // calculate the stress
521  stress = material_mat * (strain + Gmat * alpha);
522 
523  // residual of the incompatible strains
524  f += JxW[qp] * Gmat.transpose() * stress;
525 
526  // calculate the incompatible mode matrices
527  // incompatible mode diagonal stiffness matrix
528  mat5_n1n3 = material_mat * Gmat;
529  K_alphaalpha += JxW[qp] * ( Gmat.transpose() * mat5_n1n3);
530 
531  // off-diagonal coupling matrix
532  // linear strain term
533  Bmat_lin.right_multiply_transpose(mat6_n2n3, mat5_n1n3);
534  K_ualpha += JxW[qp] * mat6_n2n3;
535 
536  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
537 
538  // nonlinear component
539  // along x
540  mat7_3n3 = mat_x.transpose() * mat5_n1n3;
541  Bmat_nl_x.right_multiply_transpose(mat6_n2n3, mat7_3n3);
542  K_ualpha += JxW[qp] * mat6_n2n3;
543 
544  // along y
545  mat7_3n3 = mat_y.transpose() * mat5_n1n3;
546  Bmat_nl_y.right_multiply_transpose(mat6_n2n3, mat7_3n3);
547  K_ualpha += JxW[qp] * mat6_n2n3;
548 
549  // along z
550  mat7_3n3 = mat_z.transpose() * mat5_n1n3;
551  Bmat_nl_z.right_multiply_transpose(mat6_n2n3, mat7_3n3);
552  K_ualpha += JxW[qp] * mat6_n2n3;
553  }
554  }
555 
556 
557  // incompatible mode Jacobian inverse
558  K_alphaalpha = K_alphaalpha.inverse();
559 
560  // update the alpha values
561  alpha += K_alphaalpha * (-f - K_ualpha.transpose() * dsol.topRows(n2));
562 }
563 
564 
565 
566 bool
567 MAST::StructuralElement3D::internal_residual_sensitivity(const MAST::FunctionBase& p,
568  bool request_jacobian,
569  RealVectorX& f,
570  RealMatrixX& jac) {
571 
572  std::unique_ptr<MAST::FEBase>
573  fe(_elem.init_fe(true, false));
574 
575  const std::vector<Real>& JxW = fe->get_JxW();
576  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
577  const unsigned int
578  n_phi = (unsigned int)fe->n_shape_functions(),
579  n1 =6,
580  n2 =3*n_phi,
581  n3 =30,
582  n_nodes =_elem.get_reference_elem().n_nodes();
583 
584  RealMatrixX
585  material_mat,
586  mat_x = RealMatrixX::Zero(6,3),
587  mat_y = RealMatrixX::Zero(6,3),
588  mat_z = RealMatrixX::Zero(6,3),
589  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
590  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
591  mat3_3n2 = RealMatrixX::Zero(3, n2),
592  mat4_33 = RealMatrixX::Zero(3, 3),
593  mat5_n1n3 = RealMatrixX::Zero(n1, n3),
594  mat6_n2n3 = RealMatrixX::Zero(n2, n3),
595  mat7_3n3 = RealMatrixX::Zero(3, n3),
596  Gmat = RealMatrixX::Zero(6, n3),
597  K_alphaalpha = RealMatrixX::Zero(n3, n3),
598  K_ualpha = RealMatrixX::Zero(n2, n3),
599  K_corr = RealMatrixX::Zero(n2, n2);
600  RealVectorX
601  strain = RealVectorX::Zero(6),
602  stress = RealVectorX::Zero(6),
603  vec1_n1 = RealVectorX::Zero(n1),
604  vec2_n2 = RealVectorX::Zero(n2),
605  vec3_3 = RealVectorX::Zero(3),
606  local_disp= RealVectorX::Zero(n2),
607  f_alpha = RealVectorX::Zero(n3),
608  alpha = RealVectorX::Zero(n3);//*_incompatible_sol;
609 
610  // copy the values from the global to the local element
611  local_disp.topRows(n2) = _local_sol.topRows(n2);
612 
613  std::unique_ptr<MAST::FieldFunction<RealMatrixX> > mat_stiff =
614  _property.stiffness_A_matrix(*this);
615 
616  MAST::FEMOperatorMatrix
617  Bmat_lin,
618  Bmat_nl_x,
619  Bmat_nl_y,
620  Bmat_nl_z,
621  Bmat_nl_u,
622  Bmat_nl_v,
623  Bmat_nl_w,
624  Bmat_inc;
625  // six stress components, related to three displacements
626  Bmat_lin.reinit(n1, 3, n_nodes);
627  Bmat_nl_x.reinit(3, 3, n_nodes);
628  Bmat_nl_y.reinit(3, 3, n_nodes);
629  Bmat_nl_z.reinit(3, 3, n_nodes);
630  Bmat_nl_u.reinit(3, 3, n_nodes);
631  Bmat_nl_v.reinit(3, 3, n_nodes);
632  Bmat_nl_w.reinit(3, 3, n_nodes);
633  Bmat_inc.reinit(n1, n3, 1); // six stress-strain components
634 
636  // second for loop to calculate the residual and stiffness contributions
637  for (unsigned int qp=0; qp<JxW.size(); qp++) {
638 
639  // get the material matrix
640  mat_stiff->derivative(p, xyz[qp], _time, material_mat);
641 
642  this->initialize_green_lagrange_strain_operator(qp,
643  *fe,
644  local_disp,
645  strain,
646  mat_x, mat_y, mat_z,
647  Bmat_lin,
648  Bmat_nl_x,
649  Bmat_nl_y,
650  Bmat_nl_z,
651  Bmat_nl_u,
652  Bmat_nl_v,
653  Bmat_nl_w);
654  //this->initialize_incompatible_strain_operator(qp, *_fe, Bmat_inc, Gmat);
655 
656  // calculate the stress
657  stress = material_mat * (strain + Gmat * alpha);
658 
659  // residual from incompatible modes
660  f_alpha += JxW[qp] * Gmat.transpose() * stress;
661 
662  // calculate contribution to the residual
663  // linear strain operator
664  Bmat_lin.vector_mult_transpose(vec2_n2, stress);
665  f.topRows(n2) += JxW[qp] * vec2_n2;
666 
667  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
668 
669  // nonlinear strain operator
670  // x
671  vec3_3 = mat_x.transpose() * stress;
672  Bmat_nl_x.vector_mult_transpose(vec2_n2, vec3_3);
673  f.topRows(n2) += JxW[qp] * vec2_n2;
674 
675  // y
676  vec3_3 = mat_y.transpose() * stress;
677  Bmat_nl_y.vector_mult_transpose(vec2_n2, vec3_3);
678  f.topRows(n2) += JxW[qp] * vec2_n2;
679 
680  // z
681  vec3_3 = mat_z.transpose() * stress;
682  Bmat_nl_z.vector_mult_transpose(vec2_n2, vec3_3);
683  f.topRows(n2) += JxW[qp] * vec2_n2;
684  }
685 
686  if (request_jacobian) {
687 
688  // the strain includes the following expansion
689  // delta_epsilon = B_lin + mat_x B_x + mat_y B_y + mat_z B_z
690  // Hence, the tangent stiffness matrix will include
691  // components from epsilon^T C epsilon
692 
694  // B_lin^T C B_lin
695  Bmat_lin.left_multiply(mat1_n1n2, material_mat);
696  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
697  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
698 
699  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
700 
701  // B_x^T mat_x^T C B_lin
702  mat3_3n2 = mat_x.transpose() * mat1_n1n2;
703  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2);
704  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
705 
706  // B_y^T mat_y^T C B_lin
707  mat3_3n2 = mat_y.transpose() * mat1_n1n2;
708  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2);
709  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
710 
711  // B_z^T mat_z^T C B_lin
712  mat3_3n2 = mat_z.transpose() * mat1_n1n2;
713  Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2);
714  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
715 
717  for (unsigned int i_dim=0; i_dim<3; i_dim++) {
718  switch (i_dim) {
719  case 0:
720  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
721  break;
722 
723  case 1:
724  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
725  break;
726 
727  case 2:
728  Bmat_nl_z.left_multiply(mat1_n1n2, mat_z);
729  break;
730  }
731 
732  // B_lin^T C mat_x_i B_x_i
733  mat1_n1n2 = material_mat * mat1_n1n2;
734  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
735  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
736 
737  // B_x^T mat_x^T C mat_x B_x
738  mat3_3n2 = mat_x.transpose() * mat1_n1n2;
739  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2);
740  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
741 
742  // B_y^T mat_y^T C mat_x B_x
743  mat3_3n2 = mat_y.transpose() * mat1_n1n2;
744  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2);
745  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
746 
747  // B_z^T mat_z^T C mat_x B_x
748  mat3_3n2 = mat_z.transpose() * mat1_n1n2;
749  Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2);
750  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
751 
752  }
753 
754  // use the stress to calculate the final contribution
755  // to the Jacobian stiffness matrix
756  mat4_33(0,0) = stress(0);
757  mat4_33(1,1) = stress(1);
758  mat4_33(2,2) = stress(2);
759  mat4_33(0,1) = mat4_33(1,0) = stress(3);
760  mat4_33(1,2) = mat4_33(2,1) = stress(4);
761  mat4_33(0,2) = mat4_33(2,0) = stress(5);
762 
763  // u-disp
764  Bmat_nl_u.left_multiply(mat3_3n2, mat4_33);
765  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3_3n2);
766  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
767 
768  // v-disp
769  Bmat_nl_v.left_multiply(mat3_3n2, mat4_33);
770  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3_3n2);
771  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
772 
773  // w-disp
774  Bmat_nl_w.left_multiply(mat3_3n2, mat4_33);
775  Bmat_nl_w.right_multiply_transpose(mat2_n2n2, mat3_3n2);
776  jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
777  }
778  }
779  }
780 
781  // if jacobian is requested, add a small diagonal value for the
782  // rotational dofs
783  if (request_jacobian)
784  jac.bottomRightCorner(n2, n2) += RealMatrixX::Identity(n2, n2) *
785  1.0e-20 * jac.diagonal().maxCoeff();
786 
787  return request_jacobian;
788 }
789 
790 
791 
792 
793 bool
794 MAST::StructuralElement3D::prestress_residual (bool request_jacobian,
795  RealVectorX& f,
796  RealMatrixX& jac) {
797 
798  return request_jacobian;
799 }
800 
801 
802 
803 bool
804 MAST::StructuralElement3D::prestress_residual_sensitivity (const MAST::FunctionBase& p,
805  bool request_jacobian,
806  RealVectorX& f,
807  RealMatrixX& jac) {
808 
809  return request_jacobian;
810 }
811 
812 
813 
814 
815 bool
816 MAST::StructuralElement3D::
817 surface_pressure_residual(bool request_jacobian,
818  RealVectorX &f,
819  RealMatrixX &jac,
820  const unsigned int side,
821  MAST::BoundaryConditionBase& bc) {
822 
823  libmesh_assert(!follower_forces); // not implemented yet for follower forces
824 
825  // prepare the side finite element
826  std::unique_ptr<MAST::FEBase>
827  fe(_elem.init_side_fe(side, false, false));
828 
829  const std::vector<Real> &JxW = fe->get_JxW();
830  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
831  const std::vector<std::vector<Real> >& phi = fe->get_phi();
832  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_reference_coordinate();
833  const unsigned int
834  n_phi = (unsigned int)phi.size(),
835  n1 = 3,
836  n2 = 3*n_phi;
837 
838 
839  // get the function from this boundary condition
840  const MAST::FieldFunction<Real>& func =
841  bc.get<MAST::FieldFunction<Real> >("pressure");
842 
843 
844  FEMOperatorMatrix Bmat;
845  Real press;
846 
847  RealVectorX
848  phi_vec = RealVectorX::Zero(n_phi),
849  force = RealVectorX::Zero(n1),
850  local_f = RealVectorX::Zero(n2),
851  vec_n2 = RealVectorX::Zero(n2);
852 
853  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
854 
855  // now set the shape function values
856  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
857  phi_vec(i_nd) = phi[i_nd][qp];
858 
859  Bmat.reinit(n1, phi_vec);
860 
861  // get pressure value
862  func(qpoint[qp], _time, press);
863 
864  // calculate force
865  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
866  force(i_dim) = press * face_normals[qp](i_dim);
867 
868  Bmat.vector_mult_transpose(vec_n2, force);
869 
870  local_f += JxW[qp] * vec_n2;
871  }
872 
873  f.topRows(n2) -= local_f;
874 
875  return (request_jacobian);
876 }
877 
878 
879 
880 
881 
882 bool
883 MAST::StructuralElement3D::
884 surface_pressure_residual_sensitivity(const MAST::FunctionBase& p,
885  bool request_jacobian,
886  RealVectorX &f,
887  RealMatrixX &jac,
888  const unsigned int side,
889  MAST::BoundaryConditionBase& bc) {
890 
891  libmesh_assert(!follower_forces); // not implemented yet for follower forces
892 
893  // prepare the side finite element
894  std::unique_ptr<MAST::FEBase>
895  fe(_elem.init_side_fe(side, false, false));
896 
897  const std::vector<Real> &JxW = fe->get_JxW();
898  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
899  const std::vector<std::vector<Real> >& phi = fe->get_phi();
900  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_reference_coordinate();
901  const unsigned int
902  n_phi = (unsigned int)phi.size(),
903  n1 = 3,
904  n2 = 6*n_phi;
905 
906 
907  // get the function from this boundary condition
908  const MAST::FieldFunction<Real>& func =
909  bc.get<MAST::FieldFunction<Real> >("pressure");
910 
911 
912  FEMOperatorMatrix Bmat;
913  Real press;
914 
915  RealVectorX
916  phi_vec = RealVectorX::Zero(n_phi),
917  force = RealVectorX::Zero(2*n1),
918  local_f = RealVectorX::Zero(n2),
919  vec_n2 = RealVectorX::Zero(n2);
920 
921  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
922 
923  // now set the shape function values
924  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
925  phi_vec(i_nd) = phi[i_nd][qp];
926 
927  Bmat.reinit(2*n1, phi_vec);
928 
929  // get pressure value
930  func.derivative(p, qpoint[qp], _time, press);
931 
932  // calculate force
933  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
934  force(i_dim) = press * face_normals[qp](i_dim);
935 
936  Bmat.vector_mult_transpose(vec_n2, force);
937 
938  local_f += JxW[qp] * vec_n2;
939  }
940 
941  f -= local_f;
942 
943  return (request_jacobian);
944 }
945 
946 
947 
948 
949 
950 bool
951 MAST::StructuralElement3D::thermal_residual(bool request_jacobian,
952  RealVectorX& f,
953  RealMatrixX& jac,
954  MAST::BoundaryConditionBase& bc) {
955 
956  std::unique_ptr<MAST::FEBase>
957  fe(_elem.init_fe(true, false));
958 
959  const std::vector<Real>& JxW = fe->get_JxW();
960  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
961 
962  const unsigned int
963  n_phi = (unsigned int)fe->get_phi().size(),
964  n1 = 6,
965  n2 = 3*n_phi,
966  n_nodes = _elem.get_reference_elem().n_nodes();
967 
968  RealMatrixX
969  material_exp_A_mat,
970  mat_x = RealMatrixX::Zero(6,3),
971  mat_y = RealMatrixX::Zero(6,3),
972  mat_z = RealMatrixX::Zero(6,3),
973  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
974  mat3_3n2 = RealMatrixX::Zero(3,n2),
975  mat4_33 = RealMatrixX::Zero(3,3);
976  RealVectorX
977  vec1_n1 = RealVectorX::Zero(n1),
978  vec2_3 = RealVectorX::Zero(3),
979  vec3_n2 = RealVectorX::Zero(n2),
980  delta_t = RealVectorX::Zero(1),
981  local_disp= RealVectorX::Zero(n2),
982  strain = RealVectorX::Zero(6);
983 
984  // copy the values from the global to the local element
985  local_disp.topRows(n2) = _local_sol.topRows(n2);
986 
987  MAST::FEMOperatorMatrix
988  Bmat_lin,
989  Bmat_nl_x,
990  Bmat_nl_y,
991  Bmat_nl_z,
992  Bmat_nl_u,
993  Bmat_nl_v,
994  Bmat_nl_w;
995  // six stress components, related to three displacements
996  Bmat_lin.reinit(n1, 3, n_nodes);
997  Bmat_nl_x.reinit(3, 3, n_nodes);
998  Bmat_nl_y.reinit(3, 3, n_nodes);
999  Bmat_nl_z.reinit(3, 3, n_nodes);
1000  Bmat_nl_u.reinit(3, 3, n_nodes);
1001  Bmat_nl_v.reinit(3, 3, n_nodes);
1002  Bmat_nl_w.reinit(3, 3, n_nodes);
1003 
1004  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
1005  mat = _property.thermal_expansion_A_matrix(*this);
1006 
1007  const MAST::FieldFunction<Real>
1008  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
1009  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
1010 
1011  Real t, t0;
1012 
1013  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1014 
1015  (*mat) (xyz[qp], _time, material_exp_A_mat);
1016  temp_func (xyz[qp], _time, t);
1017  ref_temp_func(xyz[qp], _time, t0);
1018  delta_t(0) = t-t0;
1019 
1020  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)}
1021 
1022  this->initialize_green_lagrange_strain_operator(qp,
1023  *fe,
1024  local_disp,
1025  strain,
1026  mat_x, mat_y, mat_z,
1027  Bmat_lin,
1028  Bmat_nl_x,
1029  Bmat_nl_y,
1030  Bmat_nl_z,
1031  Bmat_nl_u,
1032  Bmat_nl_v,
1033  Bmat_nl_w);
1034 
1035  // linear strain operotor
1036  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
1037  f.topRows(n2) -= JxW[qp] * vec3_n2;
1038 
1039  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1040 
1041  // nonlinear strain operotor
1042  // x
1043  vec2_3 = mat_x.transpose() * vec1_n1;
1044  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec2_3);
1045  f.topRows(n2) -= JxW[qp] * vec3_n2;
1046 
1047  // y
1048  vec2_3 = mat_y.transpose() * vec1_n1;
1049  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec2_3);
1050  f.topRows(n2) -= JxW[qp] * vec3_n2;
1051 
1052  // z
1053  vec2_3 = mat_z.transpose() * vec1_n1;
1054  Bmat_nl_z.vector_mult_transpose(vec3_n2, vec2_3);
1055  f.topRows(n2) -= JxW[qp] * vec3_n2;
1056 
1057  // Jacobian for the nonlinear case
1058  if (request_jacobian) {
1059 
1060  // use the stress to calculate the final contribution
1061  // to the Jacobian stiffness matrix
1062  mat4_33(0,0) = vec1_n1(0);
1063  mat4_33(1,1) = vec1_n1(1);
1064  mat4_33(2,2) = vec1_n1(2);
1065  mat4_33(0,1) = mat4_33(1,0) = vec1_n1(3);
1066  mat4_33(1,2) = mat4_33(2,1) = vec1_n1(4);
1067  mat4_33(0,2) = mat4_33(2,0) = vec1_n1(5);
1068 
1069  // u-disp
1070  Bmat_nl_u.left_multiply(mat3_3n2, mat4_33);
1071  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3_3n2);
1072  jac.topLeftCorner(n2, n2) -= JxW[qp] * mat2_n2n2;
1073 
1074  // v-disp
1075  Bmat_nl_v.left_multiply(mat3_3n2, mat4_33);
1076  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3_3n2);
1077  jac.topLeftCorner(n2, n2) -= JxW[qp] * mat2_n2n2;
1078 
1079  // w-disp
1080  Bmat_nl_w.left_multiply(mat3_3n2, mat4_33);
1081  Bmat_nl_w.right_multiply_transpose(mat2_n2n2, mat3_3n2);
1082  jac.topLeftCorner(n2, n2) -= JxW[qp] * mat2_n2n2;
1083  }
1084  }
1085  }
1086 
1087  // Jacobian contribution from von Karman strain
1088  return request_jacobian;
1089 }
1090 
1091 
1092 
1093 bool
1094 MAST::StructuralElement3D::thermal_residual_sensitivity(const MAST::FunctionBase& p,
1095  bool request_jacobian,
1096  RealVectorX& f,
1097  RealMatrixX& jac,
1098  MAST::BoundaryConditionBase& bc) {
1099 
1100  // to be implemented
1101  libmesh_error();
1102 
1103  return false;
1104 }
1105 
1106 
1107 
1108 
1109 bool
1110 MAST::StructuralElement3D::
1111 piston_theory_residual(bool request_jacobian,
1112  RealVectorX &f,
1113  RealMatrixX& jac_xdot,
1114  RealMatrixX& jac,
1115  const unsigned int side,
1116  MAST::BoundaryConditionBase& bc) {
1117 
1118 
1119  libmesh_error(); // to be implemented
1120 
1121  return (request_jacobian);
1122 }
1123 
1124 
1125 
1126 bool
1127 MAST::StructuralElement3D::
1128 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
1129  bool request_jacobian,
1130  RealVectorX &f,
1131  RealMatrixX& jac_xdot,
1132  RealMatrixX& jac,
1133  const unsigned int side,
1134  MAST::BoundaryConditionBase& bc) {
1135 
1136 
1137  libmesh_error(); // to be implemented
1138 
1139  return (request_jacobian);
1140 }
1141 
1142 
1143 
1144 
1145 
1146 bool
1147 MAST::StructuralElement3D::calculate_stress(bool request_derivative,
1148  const MAST::FunctionBase* p,
1149  MAST::StressStrainOutputBase& output) {
1150 
1151  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
1152  std::vector<libMesh::Point> qp_loc = fe->get_qpoints();
1153 
1154 
1155  // now that the FE object has been initialized, evaluate the stress values
1156 
1157 
1158  const std::vector<Real> &JxW = fe->get_JxW();
1159  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1160  const unsigned int
1161  n_phi = (unsigned int)fe->n_shape_functions(),
1162  n1 =6,
1163  n2 =3*n_phi,
1164  n3 =30,
1165  n_nodes = _elem.get_reference_elem().n_nodes() ;
1166 
1167  RealMatrixX
1168  material_mat,
1169  mat_x = RealMatrixX::Zero(6,3),
1170  mat_y = RealMatrixX::Zero(6,3),
1171  mat_z = RealMatrixX::Zero(6,3),
1172  Gmat = RealMatrixX::Zero(6, n3);
1173 
1174  RealVectorX
1175  strain = RealVectorX::Zero(6),
1176  stress = RealVectorX::Zero(6),
1177  local_disp= RealVectorX::Zero(n2),
1178  alpha = RealVectorX::Zero(n3);//*_incompatible_sol;
1179 
1180  // copy the values from the global to the local element
1181  local_disp.topRows(n2) = _local_sol.topRows(n2);
1182 
1183  std::unique_ptr<MAST::FieldFunction<RealMatrixX> > mat_stiff =
1184  _property.stiffness_A_matrix(*this);
1185 
1186  MAST::FEMOperatorMatrix
1187  Bmat_lin,
1188  Bmat_nl_x,
1189  Bmat_nl_y,
1190  Bmat_nl_z,
1191  Bmat_nl_u,
1192  Bmat_nl_v,
1193  Bmat_nl_w,
1194  Bmat_inc;
1195  // six stress components, related to three displacements
1196  Bmat_lin.reinit(n1, 3, n_nodes);
1197  Bmat_nl_x.reinit(3, 3, n_nodes);
1198  Bmat_nl_y.reinit(3, 3, n_nodes);
1199  Bmat_nl_z.reinit(3, 3, n_nodes);
1200  Bmat_nl_u.reinit(3, 3, n_nodes);
1201  Bmat_nl_v.reinit(3, 3, n_nodes);
1202  Bmat_nl_w.reinit(3, 3, n_nodes);
1203  Bmat_inc.reinit(n1, n3, 1); // six stress-strain components
1204 
1205  // a reference to the stress output data structure
1206  MAST::StressStrainOutputBase& stress_output =
1207  dynamic_cast<MAST::StressStrainOutputBase&>(output);
1208 
1209  // initialize the incompatible mode mapping at element mid-point
1210  _init_incompatible_fe_mapping(_elem.get_reference_elem());
1211 
1213  // second for loop to calculate the residual and stiffness contributions
1214  for (unsigned int qp=0; qp<qp_loc.size(); qp++) {
1215 
1216  // get the material matrix
1217  (*mat_stiff)(xyz[qp], _time, material_mat);
1218 
1219  this->initialize_green_lagrange_strain_operator(qp,
1220  *fe,
1221  local_disp,
1222  strain,
1223  mat_x, mat_y, mat_z,
1224  Bmat_lin,
1225  Bmat_nl_x,
1226  Bmat_nl_y,
1227  Bmat_nl_z,
1228  Bmat_nl_u,
1229  Bmat_nl_v,
1230  Bmat_nl_w);
1231  //this->initialize_incompatible_strain_operator(qp, *fe, Bmat_inc, Gmat);
1232 
1233  // calculate the stress
1234  strain += Gmat * alpha;
1235  stress = material_mat * strain;
1236 
1237  // set the stress and strain data
1238  MAST::StressStrainOutputBase::Data*
1239  data = nullptr;
1240 
1241  // if neither the derivative nor sensitivity is requested, then
1242  // we assume that a new data entry is to be provided. Otherwise,
1243  // we assume that the stress at this quantity already
1244  // exists, and we only need to append sensitivity/derivative
1245  // data to it
1246  if (!request_derivative && !p)
1247  data = &(stress_output.add_stress_strain_at_qp_location(_elem,
1248  qp,
1249  qp_loc[qp],
1250  xyz[qp],
1251  stress,
1252  strain,
1253  JxW[qp]));
1254  else
1255  data = &(stress_output.get_stress_strain_data_for_elem_at_qp(_elem, qp));
1256 
1257 
1258  if (request_derivative) {
1259  // to be implemented
1260  libmesh_error();
1261  }
1262  }
1263 
1264  return request_derivative;
1265 }
1266 
1267 
1268 
1269 
1270 
1271 
1272 void
1273 MAST::StructuralElement3D::initialize_strain_operator (const unsigned int qp,
1274  const MAST::FEBase& fe,
1275  FEMOperatorMatrix& Bmat) {
1276 
1277  const std::vector<std::vector<libMesh::RealVectorValue> >&
1278  dphi = fe.get_dphi();
1279 
1280  unsigned int n_phi = (unsigned int)dphi.size();
1281  RealVectorX phi = RealVectorX::Zero(n_phi);
1282 
1283  // now set the shape function values
1284  // dN/dx
1285  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1286  phi(i_nd) = dphi[i_nd][qp](0);
1287  Bmat.set_shape_function(0, 0, phi); // epsilon_xx = du/dx
1288  Bmat.set_shape_function(3, 1, phi); // gamma_xy = dv/dx + ...
1289  Bmat.set_shape_function(5, 2, phi); // gamma_zx = dw/dx + ...
1290 
1291  // dN/dy
1292  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1293  phi(i_nd) = dphi[i_nd][qp](1);
1294  Bmat.set_shape_function(1, 1, phi); // epsilon_yy = dv/dy
1295  Bmat.set_shape_function(3, 0, phi); // gamma_xy = du/dy + ...
1296  Bmat.set_shape_function(4, 2, phi); // gamma_yz = dw/dy + ...
1297 
1298  // dN/dz
1299  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1300  phi(i_nd) = dphi[i_nd][qp](2);
1301  Bmat.set_shape_function(2, 2, phi); // epsilon_xx = dw/dz
1302  Bmat.set_shape_function(4, 1, phi); // gamma_xy = dv/dz + ...
1303  Bmat.set_shape_function(5, 0, phi); // gamma_zx = du/dz + ...
1304 }
1305 
1306 
1307 
1308 void
1309 MAST::StructuralElement3D::
1310 initialize_green_lagrange_strain_operator(const unsigned int qp,
1311  const MAST::FEBase& fe,
1312  const RealVectorX& local_disp,
1313  RealVectorX& epsilon,
1314  RealMatrixX& mat_x,
1315  RealMatrixX& mat_y,
1316  RealMatrixX& mat_z,
1317  MAST::FEMOperatorMatrix& Bmat_lin,
1318  MAST::FEMOperatorMatrix& Bmat_nl_x,
1319  MAST::FEMOperatorMatrix& Bmat_nl_y,
1320  MAST::FEMOperatorMatrix& Bmat_nl_z,
1321  MAST::FEMOperatorMatrix& Bmat_nl_u,
1322  MAST::FEMOperatorMatrix& Bmat_nl_v,
1323  MAST::FEMOperatorMatrix& Bmat_nl_w) {
1324 
1325  epsilon.setZero();
1326  mat_x.setZero();
1327  mat_y.setZero();
1328  mat_z.setZero();
1329 
1330  const std::vector<std::vector<libMesh::RealVectorValue> >&
1331  dphi = fe.get_dphi();
1332 
1333  unsigned int n_phi = (unsigned int)dphi.size();
1334  RealVectorX phi = RealVectorX::Zero(n_phi);
1335 
1336  // make sure all matrices are the right size
1337  libmesh_assert_equal_to(epsilon.size(), 6);
1338  libmesh_assert_equal_to(mat_x.rows(), 6);
1339  libmesh_assert_equal_to(mat_x.cols(), 3);
1340  libmesh_assert_equal_to(mat_y.rows(), 6);
1341  libmesh_assert_equal_to(mat_y.cols(), 3);
1342  libmesh_assert_equal_to(mat_z.rows(), 6);
1343  libmesh_assert_equal_to(mat_z.cols(), 3);
1344  libmesh_assert_equal_to(Bmat_lin.m(), 6);
1345  libmesh_assert_equal_to(Bmat_lin.n(), 3*n_phi);
1346  libmesh_assert_equal_to(Bmat_nl_x.m(), 3);
1347  libmesh_assert_equal_to(Bmat_nl_x.n(), 3*n_phi);
1348  libmesh_assert_equal_to(Bmat_nl_y.m(), 3);
1349  libmesh_assert_equal_to(Bmat_nl_y.n(), 3*n_phi);
1350  libmesh_assert_equal_to(Bmat_nl_z.m(), 3);
1351  libmesh_assert_equal_to(Bmat_nl_z.n(), 3*n_phi);
1352 
1353 
1354  // now set the shape function values
1355  // dN/dx
1356  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1357  phi(i_nd) = dphi[i_nd][qp](0);
1358  // linear strain operator
1359  Bmat_lin.set_shape_function(0, 0, phi); // epsilon_xx = du/dx
1360  Bmat_lin.set_shape_function(3, 1, phi); // gamma_xy = dv/dx + ...
1361  Bmat_lin.set_shape_function(5, 2, phi); // gamma_zx = dw/dx + ...
1362 
1363  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1364 
1365  // nonlinear strain operator in x
1366  Bmat_nl_x.set_shape_function(0, 0, phi); // du/dx
1367  Bmat_nl_x.set_shape_function(1, 1, phi); // dv/dx
1368  Bmat_nl_x.set_shape_function(2, 2, phi); // dw/dx
1369 
1370  // nonlinear strain operator in u
1371  Bmat_nl_u.set_shape_function(0, 0, phi); // du/dx
1372  Bmat_nl_v.set_shape_function(0, 1, phi); // dv/dx
1373  Bmat_nl_w.set_shape_function(0, 2, phi); // dw/dx
1374  }
1375 
1376  // dN/dy
1377  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1378  phi(i_nd) = dphi[i_nd][qp](1);
1379  // linear strain operator
1380  Bmat_lin.set_shape_function(1, 1, phi); // epsilon_yy = dv/dy
1381  Bmat_lin.set_shape_function(3, 0, phi); // gamma_xy = du/dy + ...
1382  Bmat_lin.set_shape_function(4, 2, phi); // gamma_yz = dw/dy + ...
1383 
1384  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1385 
1386  // nonlinear strain operator in y
1387  Bmat_nl_y.set_shape_function(0, 0, phi); // du/dy
1388  Bmat_nl_y.set_shape_function(1, 1, phi); // dv/dy
1389  Bmat_nl_y.set_shape_function(2, 2, phi); // dw/dy
1390 
1391  // nonlinear strain operator in v
1392  Bmat_nl_u.set_shape_function(1, 0, phi); // du/dy
1393  Bmat_nl_v.set_shape_function(1, 1, phi); // dv/dy
1394  Bmat_nl_w.set_shape_function(1, 2, phi); // dw/dy
1395  }
1396 
1397  // dN/dz
1398  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1399  phi(i_nd) = dphi[i_nd][qp](2);
1400  Bmat_lin.set_shape_function(2, 2, phi); // epsilon_xx = dw/dz
1401  Bmat_lin.set_shape_function(4, 1, phi); // gamma_xy = dv/dz + ...
1402  Bmat_lin.set_shape_function(5, 0, phi); // gamma_zx = du/dz + ...
1403 
1404  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1405 
1406  // nonlinear strain operator in z
1407  Bmat_nl_z.set_shape_function(0, 0, phi); // du/dz
1408  Bmat_nl_z.set_shape_function(1, 1, phi); // dv/dz
1409  Bmat_nl_z.set_shape_function(2, 2, phi); // dw/dz
1410 
1411  // nonlinear strain operator in w
1412  Bmat_nl_u.set_shape_function(2, 0, phi); // du/dz
1413  Bmat_nl_v.set_shape_function(2, 1, phi); // dv/dz
1414  Bmat_nl_w.set_shape_function(2, 2, phi); // dw/dz
1415 
1416 
1417  // calculate the displacement gradient to create the GL strain
1418  RealVectorX
1419  ddisp_dx = RealVectorX::Zero(3),
1420  ddisp_dy = RealVectorX::Zero(3),
1421  ddisp_dz = RealVectorX::Zero(3);
1422 
1423  Bmat_nl_x.vector_mult(ddisp_dx, local_disp); // {du/dx, dv/dx, dw/dx}
1424  Bmat_nl_y.vector_mult(ddisp_dy, local_disp); // {du/dy, dv/dy, dw/dy}
1425  Bmat_nl_z.vector_mult(ddisp_dz, local_disp); // {du/dz, dv/dz, dw/dz}
1426 
1427  // prepare the displacement gradient matrix: F = grad(u)
1428  RealMatrixX
1429  F = RealMatrixX::Zero(3,3),
1430  E = RealMatrixX::Zero(3,3);
1431  F.col(0) = ddisp_dx;
1432  F.col(1) = ddisp_dy;
1433  F.col(2) = ddisp_dz;
1434 
1435  // this calculates the Green-Lagrange strain in the reference config
1436  E = 0.5*(F + F.transpose() + F.transpose() * F);
1437 
1438  // now, add this to the strain vector
1439  epsilon(0) = E(0,0);
1440  epsilon(1) = E(1,1);
1441  epsilon(2) = E(2,2);
1442  epsilon(3) = E(0,1) + E(1,0);
1443  epsilon(4) = E(1,2) + E(2,1);
1444  epsilon(5) = E(0,2) + E(2,0);
1445 
1446  // now initialize the matrices with strain components
1447  // that multiply the Bmat_nl terms
1448  mat_x.row(0) = ddisp_dx;
1449  mat_x.row(3) = ddisp_dy;
1450  mat_x.row(5) = ddisp_dz;
1451 
1452  mat_y.row(1) = ddisp_dy;
1453  mat_y.row(3) = ddisp_dx;
1454  mat_y.row(4) = ddisp_dz;
1455 
1456 
1457  mat_z.row(2) = ddisp_dz;
1458  mat_z.row(4) = ddisp_dy;
1459  mat_z.row(5) = ddisp_dx;
1460  }
1461  else
1462  Bmat_lin.vector_mult(epsilon, local_disp);
1463 }
1464 
1465 
1466 
1467 
1468 void
1469 MAST::StructuralElement3D::
1470 initialize_incompatible_strain_operator(const unsigned int qp,
1471  const MAST::FEBase& fe,
1472  FEMOperatorMatrix& Bmat,
1473  RealMatrixX& G_mat) {
1474 
1475  RealVectorX phi_vec = RealVectorX::Zero(1);
1476 
1477  // get the location of element coordinates
1478  const std::vector<libMesh::Point>& q_point = fe.get_qrule().get_points();
1479  const Real
1480  xi = q_point[qp](0),
1481  eta = q_point[qp](1),
1482  phi = q_point[qp](2);
1483 
1484  const std::vector<std::vector<Real> >&
1485  dshapedxi = fe.get_dphidxi(),
1486  dshapedeta = fe.get_dphideta(),
1487  dshapedphi = fe.get_dphidzeta();
1488 
1489  const libMesh::Elem& e = _elem.get_reference_elem();
1490 
1491  // calculate the deformed xyz coordinates
1492  const unsigned int n_nodes = e.n_nodes();
1493  RealVectorX
1494  xdef = RealVectorX::Zero(n_nodes),
1495  ydef = RealVectorX::Zero(n_nodes),
1496  zdef = RealVectorX::Zero(n_nodes),
1497  phivec = RealVectorX::Zero(n_nodes);
1498 
1499  // set the current values of nodal coordinates
1500  for (unsigned int i_node=0; i_node<n_nodes; i_node++) {
1501  xdef(i_node) = e.point(i_node)(0);// + _sol(i_node*3+0);
1502  ydef(i_node) = e.point(i_node)(1);// + _sol(i_node*3+1);
1503  zdef(i_node) = e.point(i_node)(2);// + _sol(i_node*3+2);
1504  }
1505 
1506  // calculate dxyz/dxi
1507  // make sure that the number of shape functions is the same as the number
1508  // of nodes. Meaning that this formulation is limited to Lagrange
1509  // elemnts only.
1510  libmesh_assert_equal_to(dshapedxi.size(), n_nodes);
1511 
1512  RealMatrixX
1513  jac = RealMatrixX::Zero(3,3);
1514 
1515  // first derivatives wrt xi
1516  for (unsigned int i_node=0; i_node<n_nodes; i_node++)
1517  phivec(i_node) = dshapedxi[i_node][qp];
1518 
1519  jac(0,0) = phivec.dot(xdef);
1520  jac(0,1) = phivec.dot(ydef);
1521  jac(0,2) = phivec.dot(zdef);
1522 
1523 
1524  // second, derivatives wrt eta
1525  for (unsigned int i_node=0; i_node<n_nodes; i_node++)
1526  phivec(i_node) = dshapedeta[i_node][qp];
1527 
1528  jac(1,0) = phivec.dot(xdef);
1529  jac(1,1) = phivec.dot(ydef);
1530  jac(1,2) = phivec.dot(zdef);
1531 
1532  // lastly, derivatives wrt phi
1533  for (unsigned int i_node=0; i_node<n_nodes; i_node++)
1534  phivec(i_node) = dshapedphi[i_node][qp];
1535 
1536  jac(2,0) = phivec.dot(xdef);
1537  jac(2,1) = phivec.dot(ydef);
1538  jac(2,2) = phivec.dot(zdef);
1539 
1540 
1541  // now set the shape function values
1542  // epsilon_xx
1543  phi_vec(0) = xi; Bmat.set_shape_function(0, 0, phi_vec);
1544  phi_vec(0) = xi*eta; Bmat.set_shape_function(0, 15, phi_vec);
1545  phi_vec(0) = xi*phi; Bmat.set_shape_function(0, 16, phi_vec);
1546  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(0, 24, phi_vec);
1547 
1548 
1549  // epsilon_yy
1550  phi_vec(0) = eta; Bmat.set_shape_function(1, 1, phi_vec);
1551  phi_vec(0) = xi*eta; Bmat.set_shape_function(1, 17, phi_vec);
1552  phi_vec(0) = eta*phi; Bmat.set_shape_function(1, 18, phi_vec);
1553  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(1, 25, phi_vec);
1554 
1555  // epsilon_zz
1556  phi_vec(0) = phi; Bmat.set_shape_function(2, 2, phi_vec);
1557  phi_vec(0) = xi*phi; Bmat.set_shape_function(2, 19, phi_vec);
1558  phi_vec(0) = eta*phi; Bmat.set_shape_function(2, 20, phi_vec);
1559  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(2, 26, phi_vec);
1560 
1561  // epsilon_xy
1562  phi_vec(0) = xi; Bmat.set_shape_function(3, 3, phi_vec);
1563  phi_vec(0) = eta; Bmat.set_shape_function(3, 4, phi_vec);
1564  phi_vec(0) = xi*phi; Bmat.set_shape_function(3, 9, phi_vec);
1565  phi_vec(0) = eta*phi; Bmat.set_shape_function(3, 10, phi_vec);
1566  phi_vec(0) = xi*eta; Bmat.set_shape_function(3, 21, phi_vec);
1567  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(3, 27, phi_vec);
1568 
1569  // epsilon_yz
1570  phi_vec(0) = eta; Bmat.set_shape_function(4, 5, phi_vec);
1571  phi_vec(0) = phi; Bmat.set_shape_function(4, 6, phi_vec);
1572  phi_vec(0) = xi*eta; Bmat.set_shape_function(4, 11, phi_vec);
1573  phi_vec(0) = xi*phi; Bmat.set_shape_function(4, 12, phi_vec);
1574  phi_vec(0) = eta*phi; Bmat.set_shape_function(4, 22, phi_vec);
1575  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(4, 28, phi_vec);
1576 
1577  // epsilon_xz
1578  phi_vec(0) = xi; Bmat.set_shape_function(5, 7, phi_vec);
1579  phi_vec(0) = phi; Bmat.set_shape_function(5, 8, phi_vec);
1580  phi_vec(0) = xi*eta; Bmat.set_shape_function(5, 13, phi_vec);
1581  phi_vec(0) = eta*phi; Bmat.set_shape_function(5, 14, phi_vec);
1582  phi_vec(0) = xi*phi; Bmat.set_shape_function(5, 23, phi_vec);
1583  phi_vec(0) = xi*eta*phi; Bmat.set_shape_function(5, 29, phi_vec);
1584 
1585 
1586  Bmat.left_multiply(G_mat, _T0_inv_tr);
1587  G_mat /= jac.determinant();
1588 }
1589 
1590 
1591 
1592 
1593 void
1594 MAST::StructuralElement3D::_init_incompatible_fe_mapping( const libMesh::Elem& e) {
1595 
1596  libmesh_assert(e.type() == libMesh::HEX8);
1597 
1598  unsigned int nv = _system.system().n_vars();
1599 
1600  libmesh_assert (nv);
1601  libMesh::FEType fe_type = _system.system().variable_type(0); // all variables are assumed to be of same type
1602 
1603 
1604  for (unsigned int i=1; i != nv; ++i)
1605  libmesh_assert(fe_type == _system.system().variable_type(i));
1606 
1607  // Create an adequate quadrature rule
1608  std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(e.dim(), fe_type).release());
1609  const std::vector<libMesh::Point> pts(1); // initializes point to (0,0,0)
1610 
1611  fe->get_dphidxi();
1612  fe->get_dphideta();
1613  fe->get_dphidzeta();
1614 
1615  fe->reinit(&e, &pts); // reinit at (0,0,0)
1616 
1617  _T0_inv_tr = RealMatrixX::Zero(6,6);
1618 
1619  const std::vector<std::vector<Real> >&
1620  dshapedxi = fe->get_dphidxi(),
1621  dshapedeta = fe->get_dphideta(),
1622  dshapedphi = fe->get_dphidzeta();
1623 
1624  // calculate the deformed xyz coordinates
1625  RealVectorX
1626  xdef = RealVectorX::Zero(e.n_nodes()),
1627  ydef = RealVectorX::Zero(e.n_nodes()),
1628  zdef = RealVectorX::Zero(e.n_nodes()),
1629  phi = RealVectorX::Zero(e.n_nodes());
1630 
1631  // set the current values of nodal coordinates
1632  for (unsigned int i_node=0; i_node<e.n_nodes(); i_node++) {
1633  xdef(i_node) = e.point(i_node)(0);// + _local_sol(i_node*3+0);
1634  ydef(i_node) = e.point(i_node)(1);// + _local_sol(i_node*3+1);
1635  zdef(i_node) = e.point(i_node)(2);// + _local_sol(i_node*3+2);
1636  }
1637 
1638  // calculate dxyz/dxi
1639  // make sure that the number of shape functions is the same as the number
1640  // of nodes. Meaning that this formulation is limited to Lagrange
1641  // elemnts only.
1642  libmesh_assert_equal_to(dshapedxi.size(), e.n_nodes());
1643 
1644  RealMatrixX
1645  jac = RealMatrixX::Zero(3,3),
1646  T0 = RealMatrixX::Zero(6,6);
1647 
1648  // first derivatives wrt xi
1649  for (unsigned int i_node=0; i_node<e.n_nodes(); i_node++)
1650  phi(i_node) = dshapedxi[i_node][0];
1651 
1652  jac(0,0) = phi.dot(xdef);
1653  jac(0,1) = phi.dot(ydef);
1654  jac(0,2) = phi.dot(zdef);
1655 
1656 
1657  // second, derivatives wrt eta
1658  for (unsigned int i_node=0; i_node<e.n_nodes(); i_node++)
1659  phi(i_node) = dshapedeta[i_node][0];
1660 
1661  jac(1,0) = phi.dot(xdef);
1662  jac(1,1) = phi.dot(ydef);
1663  jac(1,2) = phi.dot(zdef);
1664 
1665  // lastly, derivatives wrt phi
1666  for (unsigned int i_node=0; i_node<e.n_nodes(); i_node++)
1667  phi(i_node) = dshapedphi[i_node][0];
1668 
1669  jac(2,0) = phi.dot(xdef);
1670  jac(2,1) = phi.dot(ydef);
1671  jac(2,2) = phi.dot(zdef);
1672 
1673 
1674  // we first set the values of the T0 matrix and then get its inverse
1675  T0(0,0) = jac(0,0)*jac(0,0);
1676  T0(0,1) = jac(1,0)*jac(1,0);
1677  T0(0,2) = jac(2,0)*jac(2,0);
1678  T0(0,3) = 2*jac(0,0)*jac(1,0);
1679  T0(0,4) = 2*jac(1,0)*jac(2,0);
1680  T0(0,5) = 2*jac(2,0)*jac(0,0);
1681 
1682  T0(1,0) = jac(0,1)*jac(0,1);
1683  T0(1,1) = jac(1,1)*jac(1,1);
1684  T0(1,2) = jac(2,1)*jac(2,1);
1685  T0(1,3) = 2*jac(0,1)*jac(1,1);
1686  T0(1,4) = 2*jac(1,1)*jac(2,1);
1687  T0(1,5) = 2*jac(2,1)*jac(0,1);
1688 
1689  T0(2,0) = jac(0,2)*jac(0,2);
1690  T0(2,1) = jac(1,2)*jac(1,2);
1691  T0(2,2) = jac(2,2)*jac(2,2);
1692  T0(2,3) = 2*jac(0,2)*jac(1,2);
1693  T0(2,4) = 2*jac(1,2)*jac(2,2);
1694  T0(2,5) = 2*jac(2,2)*jac(0,2);
1695 
1696  T0(3,0) = jac(0,0)*jac(0,1);
1697  T0(3,1) = jac(1,0)*jac(1,1);
1698  T0(3,2) = jac(2,0)*jac(2,1);
1699  T0(3,3) = jac(0,0)*jac(1,1)+jac(1,0)*jac(0,1);
1700  T0(3,4) = jac(1,0)*jac(2,1)+jac(2,0)*jac(1,1);
1701  T0(3,5) = jac(2,1)*jac(0,0)+jac(2,0)*jac(0,1);
1702 
1703  T0(4,0) = jac(0,1)*jac(0,2);
1704  T0(4,1) = jac(1,1)*jac(1,2);
1705  T0(4,2) = jac(2,1)*jac(2,2);
1706  T0(4,3) = jac(0,1)*jac(1,2)+jac(1,1)*jac(0,2);
1707  T0(4,4) = jac(1,1)*jac(2,2)+jac(2,1)*jac(1,2);
1708  T0(4,5) = jac(2,2)*jac(0,1)+jac(2,1)*jac(0,2);
1709 
1710  T0(5,0) = jac(0,0)*jac(0,2);
1711  T0(5,1) = jac(1,0)*jac(1,2);
1712  T0(5,2) = jac(2,0)*jac(2,2);
1713  T0(5,3) = jac(0,0)*jac(1,2)+jac(1,0)*jac(0,2);
1714  T0(5,4) = jac(1,0)*jac(2,2)+jac(2,0)*jac(1,2);
1715  T0(5,5) = jac(2,2)*jac(0,0)+jac(2,0)*jac(0,2);
1716 
1717  _T0_inv_tr = jac.determinant() * T0.inverse().transpose();
1718 }
1719 
1720 
MAST::SystemInitialization::system
MAST::NonlinearSystem & system()
Definition: system_initialization.h:73
MAST::ElementBase::_elem
const MAST::GeomElem & _elem
geometric element for which the computations are performed
Definition: elem_base.h:205
MAST::StressStrainOutputBase
Data structure provides the mechanism to store stress and strain output from a structural analysis...
Definition: stress_output_base.h:61
MAST::StructuralElement3D::_T0_inv_tr
RealMatrixX _T0_inv_tr
Jacobian matrix at element center needed for incompatible modes.
Definition: solid_element_3d.h:418
fe_base.h
MAST::GeomElem::get_reference_elem
virtual const libMesh::Elem & get_reference_elem() const
Definition: geom_elem.cpp:54
MAST::FEMOperatorMatrix::reinit
void reinit(unsigned int n_interpolated_vars, unsigned int n_discrete_vars, unsigned int n_discrete_dofs_per_var)
this initializes the operator for number of rows and variables, assuming that all variables has the s...
Definition: fem_operator_matrix.h:210
MAST::StructuralElement3D::calculate_stress
virtual bool calculate_stress(bool request_derivative, const MAST::FunctionBase *p, MAST::StressStrainOutputBase &output)
Calculates the stress tensor.
Definition: solid_element_3d.cpp:1147
MAST::StructuralElement3D::inertial_residual
virtual bool inertial_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xddot, RealMatrixX &jac_xdot, RealMatrixX &jac)
Calculates the inertial force and the Jacobian matrices.
Definition: solid_element_3d.cpp:49
MAST::StructuralElement3D::internal_residual_sensitivity
virtual bool internal_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the sensitivity of internal residual vector and Jacobian due to strain energy...
Definition: solid_element_3d.cpp:567
MAST::StructuralElement3D::surface_pressure_residual_sensitivity
virtual bool surface_pressure_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to surface pressure.
Definition: solid_element_3d.cpp:884
MAST::FEMOperatorMatrix::set_shape_function
void set_shape_function(unsigned int interpolated_var, unsigned int discrete_var, const RealVectorX &shape_func)
sets the shape function values for the block corresponding to interpolated_var and discrete_var...
Definition: fem_operator_matrix.h:228
MAST::StructuralElement3D::piston_theory_residual_sensitivity
virtual bool piston_theory_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian sensitivity due to piston-theory based surface pressure on t...
Definition: solid_element_3d.h:278
fem_operator_matrix.h
MAST::StructuralElementBase::_local_accel
RealVectorX _local_accel
local acceleration
Definition: structural_element_base.h:800
MAST::FieldFunction::derivative
virtual void derivative(const MAST::FunctionBase &f, ValType &v) const
calculates the value of the function derivative and returns it in v.
Definition: field_function_base.h:78
Real
libMesh::Real Real
Definition: mast_data_types.h:32
MAST::ElementBase::_system
MAST::SystemInitialization & _system
SystemInitialization object associated with this element.
Definition: elem_base.h:200
MAST::FEBase::get_dphi
virtual const std::vector< std::vector< libMesh::RealVectorValue > > & get_dphi() const
Definition: fe_base.cpp:255
nonlinear_system.h
MAST::FEBase::get_dphidxi
virtual const std::vector< std::vector< Real > > & get_dphidxi() const
Definition: fe_base.cpp:367
MAST::StructuralElement3D::StructuralElement3D
StructuralElement3D(MAST::SystemInitialization &sys, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
Definition: solid_element_3d.cpp:35
MAST::StructuralElement3D::initialize_strain_operator
void initialize_strain_operator(const unsigned int qp, const MAST::FEBase &fe, FEMOperatorMatrix &Bmat)
initialize strain operator matrix
Definition: solid_element_3d.cpp:1273
MAST::StructuralElement3D::_init_incompatible_fe_mapping
void _init_incompatible_fe_mapping(const libMesh::Elem &e)
initialize the Jacobian needed for incompatible modes
Definition: solid_element_3d.cpp:1594
system_initialization.h
MAST::StructuralElement3D::initialize_incompatible_strain_operator
void initialize_incompatible_strain_operator(const unsigned int qp, const MAST::FEBase &fe, FEMOperatorMatrix &Bmat, RealMatrixX &G_mat)
initialize incompatible strain operator
Definition: solid_element_3d.cpp:1470
solid_element_3d.h
MAST::FEMOperatorMatrix::n
unsigned int n() const
Definition: fem_operator_matrix.h:52
stress_output_base.h
MAST::ElementPropertyCardBase::inertia_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > inertia_matrix(const MAST::ElementBase &e) const =0
MAST::FEBase::get_dphideta
virtual const std::vector< std::vector< Real > > & get_dphideta() const
Definition: fe_base.cpp:375
MAST::GeomElem::init_fe
virtual std::unique_ptr< MAST::FEBase > init_fe(bool init_grads, bool init_second_order_derivative, int extra_quadrature_order=0) const
initializes the finite element shape function and quadrature object with the order of quadrature rule...
Definition: geom_elem.cpp:148
MAST::StructuralElement3D::piston_theory_residual
virtual bool piston_theory_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to piston-theory based surface pressure on the entire el...
Definition: solid_element_3d.h:261
geom_elem.h
MAST::StructuralElementBase::follower_forces
bool follower_forces
flag for follower forces
Definition: structural_element_base.h:439
MAST::StructuralElement3D::thermal_residual_sensitivity
virtual bool thermal_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the sensitivity of residual vector and Jacobian due to thermal stresses.
Definition: solid_element_3d.cpp:1094
RealMatrixX
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Definition: mast_data_types.h:44
boundary_condition_base.h
MAST::StructuralElement3D::internal_residual
virtual bool internal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal residual vector and Jacobian due to strain energy.
Definition: solid_element_3d.cpp:136
MAST::StructuralElement3D::surface_pressure_residual
bool surface_pressure_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to surface pressure.
Definition: solid_element_3d.cpp:817
MAST::StructuralElement3D::thermal_residual
virtual bool thermal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the residual vector and Jacobian due to thermal stresses.
Definition: solid_element_3d.cpp:951
MAST::ElementPropertyCardBase::stiffness_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > stiffness_A_matrix(const MAST::ElementBase &e) const =0
MAST::FEMOperatorMatrix::vector_mult_transpose
void vector_mult_transpose(T &res, const T &v) const
res = v^T * [this]
Definition: fem_operator_matrix.h:303
RealVectorX
Matrix< Real, Dynamic, 1 > RealVectorX
Definition: mast_data_types.h:35
MAST::FEMOperatorMatrix::right_multiply_transpose
void right_multiply_transpose(T &r, const T &m) const
[R] = [this]^T * [M]
Definition: fem_operator_matrix.h:355
assembly_base.h
MAST::StructuralElementBase::_property
const MAST::ElementPropertyCardBase & _property
element property
Definition: structural_element_base.h:746
MAST::StressStrainOutputBase::Data
This class provides a mechanism to store stress/strain values, their derivatives and sensitivity valu...
Definition: stress_output_base.h:72
MAST::StructuralElement3D::prestress_residual
virtual bool prestress_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the prestress residual vector and Jacobian.
Definition: solid_element_3d.cpp:794
MAST::GeomElem
This class acts as a wrapper around libMesh::Elem for the purpose of providing a uniform interface fo...
Definition: geom_elem.h:59
MAST::ElementPropertyCardBase::strain_type
const MAST::StrainType strain_type() const
returns the type of strain to be used for this element
Definition: element_property_card_base.h:192
MAST::ElementPropertyCardBase::if_diagonal_mass_matrix
bool if_diagonal_mass_matrix() const
returns the type of strain to be used for this element
Definition: element_property_card_base.h:217
MAST::FEMOperatorMatrix::vector_mult
void vector_mult(T &res, const T &v) const
res = [this] * v
Definition: fem_operator_matrix.h:280
MAST::StructuralElement3D::initialize_green_lagrange_strain_operator
void initialize_green_lagrange_strain_operator(const unsigned int qp, const MAST::FEBase &fe, const RealVectorX &local_disp, RealVectorX &epsilon, RealMatrixX &mat_x, RealMatrixX &mat_y, RealMatrixX &mat_z, MAST::FEMOperatorMatrix &Bmat_lin, MAST::FEMOperatorMatrix &Bmat_nl_x, MAST::FEMOperatorMatrix &Bmat_nl_y, MAST::FEMOperatorMatrix &Bmat_nl_z, MAST::FEMOperatorMatrix &Bmat_nl_u, MAST::FEMOperatorMatrix &Bmat_nl_v, MAST::FEMOperatorMatrix &Bmat_nl_w)
initialize the strain operator matrices for the Green-Lagrange strain matrices
Definition: solid_element_3d.cpp:1310
MAST::FunctionSetBase::get
const ValType & get(const std::string &nm) const
returns a constant reference to the specified function
Definition: function_set_base.h:64
MAST::GeomElem::init_side_fe
virtual std::unique_ptr< MAST::FEBase > init_side_fe(unsigned int s, bool init_grads, bool init_second_order_derivative, int extra_quadrature_order=0) const
initializes the finite element shape function and quadrature object for the side with the order of qu...
Definition: geom_elem.cpp:165
MAST::ElementPropertyCardBase::thermal_expansion_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > thermal_expansion_A_matrix(const MAST::ElementBase &e) const =0
field_function_base.h
MAST::StressStrainOutputBase::get_stress_strain_data_for_elem_at_qp
virtual MAST::StressStrainOutputBase::Data & get_stress_strain_data_for_elem_at_qp(const MAST::GeomElem &e, const unsigned int qp)
Definition: stress_output_base.cpp:798
MAST::FEMOperatorMatrix::m
unsigned int m() const
Definition: fem_operator_matrix.h:50
MAST::FEMOperatorMatrix::left_multiply
void left_multiply(T &r, const T &m) const
[R] = [M] * [this]
Definition: fem_operator_matrix.h:415
MAST::StructuralElementBase::_local_sol
RealVectorX _local_sol
local solution
Definition: structural_element_base.h:752
element_property_card_base.h
MAST::StressStrainOutputBase::add_stress_strain_at_qp_location
virtual MAST::StressStrainOutputBase::Data & add_stress_strain_at_qp_location(const MAST::GeomElem &e, const unsigned int qp, const libMesh::Point &quadrature_pt, const libMesh::Point &physical_pt, const RealVectorX &stress, const RealVectorX &strain, Real JxW)
add the stress tensor associated with the qp.
Definition: stress_output_base.cpp:613
MAST::ElementBase::_time
const Real & _time
time for which system is being assembled
Definition: elem_base.h:219
MAST::FEBase::get_dphidzeta
virtual const std::vector< std::vector< Real > > & get_dphidzeta() const
Definition: fe_base.cpp:383
MAST::StructuralElement3D::update_incompatible_mode_solution
virtual void update_incompatible_mode_solution(const RealVectorX &dsol)
updates the incompatible solution for this element.
Definition: solid_element_3d.cpp:430
MAST::FEBase::get_qrule
virtual const libMesh::QBase & get_qrule() const
Definition: fe_base.cpp:418
MAST::StructuralElement3D::prestress_residual_sensitivity
virtual bool prestress_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the sensitivity prestress residual vector and Jacobian.
Definition: solid_element_3d.cpp:804