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