MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
structural_element_2d.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/structural_element_2d.h"
22 #include "elasticity/piston_theory_boundary_condition.h"
23 #include "elasticity/stress_output_base.h"
24 #include "elasticity/bending_operator.h"
25 #include "property_cards/element_property_card_2D.h"
26 #include "property_cards/material_property_card_base.h"
27 #include "numerics/fem_operator_matrix.h"
28 #include "mesh/fe_base.h"
29 #include "mesh/geom_elem.h"
30 #include "base/system_initialization.h"
31 #include "base/boundary_condition_base.h"
32 #include "base/parameter.h"
33 #include "base/constant_field_function.h"
34 #include "base/assembly_base.h"
35 #include "level_set/level_set_boundary_velocity.h"
36 
37 
38 MAST::StructuralElement2D::
39 StructuralElement2D(MAST::SystemInitialization& sys,
40  const MAST::GeomElem& elem,
41  const MAST::ElementPropertyCardBase& p):
42 MAST::BendingStructuralElem(sys, elem, p) {
43 
44 }
45 
46 
47 
48 MAST::StructuralElement2D::~StructuralElement2D() {
49 
50 }
51 
52 
53 
54 
55 void
56 MAST::StructuralElement2D::
57 initialize_von_karman_strain_operator(const unsigned int qp,
58  const MAST::FEBase& fe,
59  RealVectorX& vk_strain,
60  RealMatrixX& vk_dwdxi_mat,
61  MAST::FEMOperatorMatrix& Bmat_vk) {
62 
63  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
64  const unsigned int n_phi = (unsigned int)dphi.size();
65 
66  libmesh_assert_equal_to(vk_strain.size(), 3);
67  libmesh_assert_equal_to(vk_dwdxi_mat.rows(), 3);
68  libmesh_assert_equal_to(vk_dwdxi_mat.cols(), 2);
69  libmesh_assert_equal_to(Bmat_vk.m(), 2);
70  libmesh_assert_equal_to(Bmat_vk.n(), 6*n_phi);
71  libmesh_assert_less (qp, dphi[0].size());
72 
73  Real dw=0.;
74  vk_strain.setZero();
75  vk_dwdxi_mat.setZero();
76 
77  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
78 
79  dw = 0.;
80  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
81  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
82  dw += phi_vec(i_nd)*_local_sol(2*n_phi+i_nd); // dw/dx
83  }
84  Bmat_vk.set_shape_function(0, 2, phi_vec); // dw/dx
85  vk_dwdxi_mat(0, 0) = dw; // epsilon-xx : dw/dx
86  vk_dwdxi_mat(2, 1) = dw; // gamma-xy : dw/dx
87  vk_strain(0) = 0.5*dw*dw; // 1/2 * (dw/dx)^2
88  vk_strain(2) = dw; // (dw/dx)*(dw/dy) only dw/dx is provided here
89 
90  dw = 0.;
91  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
92  phi_vec(i_nd) = dphi[i_nd][qp](1); // dphi/dy
93  dw += phi_vec(i_nd)*_local_sol(2*n_phi+i_nd); // dw/dy
94  }
95  Bmat_vk.set_shape_function(1, 2, phi_vec); // dw/dy
96  vk_dwdxi_mat(1, 1) = dw; // epsilon-yy : dw/dy
97  vk_dwdxi_mat(2, 0) = dw; // gamma-xy : dw/dy
98  vk_strain(1) = 0.5*dw*dw; // 1/2 * (dw/dy)^2
99  vk_strain(2) *= dw; // (dw/dx)*(dw/dy)
100 }
101 
102 
103 
104 
105 void
106 MAST::StructuralElement2D::
107 initialize_von_karman_strain_operator_sensitivity(const unsigned int qp,
108  const MAST::FEBase& fe,
109  RealMatrixX &vk_dwdxi_mat_sens) {
110 
111  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
112  const unsigned int n_phi = (unsigned int)dphi.size();
113 
114  libmesh_assert_equal_to(vk_dwdxi_mat_sens.rows(), 3);
115  libmesh_assert_equal_to(vk_dwdxi_mat_sens.cols(), 2);
116  libmesh_assert_less (qp, dphi[0].size());
117 
118  Real dw=0.;
119  vk_dwdxi_mat_sens.setZero();
120 
121  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
122 
123  dw = 0.;
124  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
125  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
126  dw += phi_vec(i_nd)*_local_sol_sens(2*n_phi+i_nd); // dw/dx
127  }
128  vk_dwdxi_mat_sens(0, 0) = dw; // epsilon-xx : dw/dx
129  vk_dwdxi_mat_sens(2, 1) = dw; // gamma-xy : dw/dx
130 
131  dw = 0.;
132  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
133  phi_vec(i_nd) = dphi[i_nd][qp](1); // dphi/dy
134  dw += phi_vec(i_nd)*_local_sol_sens(2*n_phi+i_nd); // dw/dy
135  }
136  vk_dwdxi_mat_sens(1, 1) = dw; // epsilon-yy : dw/dy
137  vk_dwdxi_mat_sens(2, 0) = dw; // gamma-xy : dw/dy
138 }
139 
140 
141 void
142 MAST::StructuralElement2D::
143 initialize_green_lagrange_strain_operator(const unsigned int qp,
144  const MAST::FEBase& fe,
145  const RealVectorX& local_disp,
146  RealVectorX& epsilon,
147  RealMatrixX& mat_x,
148  RealMatrixX& mat_y,
149  MAST::FEMOperatorMatrix& Bmat_lin,
150  MAST::FEMOperatorMatrix& Bmat_nl_x,
151  MAST::FEMOperatorMatrix& Bmat_nl_y,
152  MAST::FEMOperatorMatrix& Bmat_nl_u,
153  MAST::FEMOperatorMatrix& Bmat_nl_v) {
154 
155 
156  epsilon.setZero();
157  mat_x.setZero();
158  mat_y.setZero();
159 
160  const std::vector<std::vector<libMesh::RealVectorValue> >&
161  dphi = fe.get_dphi();
162 
163  unsigned int n_phi = (unsigned int)dphi.size();
164  RealVectorX phi = RealVectorX::Zero(n_phi);
165 
166  // make sure all matrices are the right size
167  libmesh_assert_equal_to(epsilon.size(), 3);
168  libmesh_assert_equal_to(mat_x.rows(), 3);
169  libmesh_assert_equal_to(mat_x.cols(), 2);
170  libmesh_assert_equal_to(mat_y.rows(), 3);
171  libmesh_assert_equal_to(mat_y.cols(), 2);
172  libmesh_assert_equal_to(Bmat_lin.m(), 3);
173  libmesh_assert_equal_to(Bmat_lin.n(), 6*n_phi);
174  libmesh_assert_equal_to(Bmat_nl_x.m(), 2);
175  libmesh_assert_equal_to(Bmat_nl_x.n(), 6*n_phi);
176  libmesh_assert_equal_to(Bmat_nl_y.m(), 2);
177  libmesh_assert_equal_to(Bmat_nl_y.n(), 6*n_phi);
178 
179 
180  // now set the shape function values
181  // dN/dx
182  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
183  phi(i_nd) = dphi[i_nd][qp](0);
184  // linear strain operator
185  Bmat_lin.set_shape_function(0, 0, phi); // epsilon_xx = du/dx
186  Bmat_lin.set_shape_function(2, 1, phi); // gamma_xy = dv/dx + ...
187 
188  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
189 
190  // nonlinear strain operator in x
191  Bmat_nl_x.set_shape_function(0, 0, phi); // du/dx
192  Bmat_nl_x.set_shape_function(1, 1, phi); // dv/dx
193 
194  // nonlinear strain operator in u
195  Bmat_nl_u.set_shape_function(0, 0, phi); // du/dx
196  Bmat_nl_v.set_shape_function(0, 1, phi); // dv/dx
197  }
198 
199  // dN/dy
200  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
201  phi(i_nd) = dphi[i_nd][qp](1);
202  // linear strain operator
203  Bmat_lin.set_shape_function(1, 1, phi); // epsilon_yy = dv/dy
204  Bmat_lin.set_shape_function(2, 0, phi); // gamma_xy = du/dy + ...
205 
206  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
207 
208  // nonlinear strain operator in y
209  Bmat_nl_y.set_shape_function(0, 0, phi); // du/dy
210  Bmat_nl_y.set_shape_function(1, 1, phi); // dv/dy
211 
212  // nonlinear strain operator in v
213  Bmat_nl_u.set_shape_function(1, 0, phi); // du/dy
214  Bmat_nl_v.set_shape_function(1, 1, phi); // dv/dy
215 
216  // calculate the displacement gradient to create the
217  RealVectorX
218  ddisp_dx = RealVectorX::Zero(2),
219  ddisp_dy = RealVectorX::Zero(2);
220 
221  Bmat_nl_x.vector_mult(ddisp_dx, local_disp); // {du/dx, dv/dx, dw/dx}
222  Bmat_nl_y.vector_mult(ddisp_dy, local_disp); // {du/dy, dv/dy, dw/dy}
223 
224  // prepare the deformation gradient matrix
225  RealMatrixX
226  F = RealMatrixX::Zero(2,2),
227  E = RealMatrixX::Zero(2,2);
228  F.col(0) = ddisp_dx;
229  F.col(1) = ddisp_dy;
230 
231  // this calculates the Green-Lagrange strain in the reference config
232  E = 0.5*(F + F.transpose() + F.transpose() * F);
233 
234  // now, add this to the strain vector
235  epsilon(0) = E(0,0);
236  epsilon(1) = E(1,1);
237  epsilon(2) = E(0,1) + E(1,0);
238 
239  // now initialize the matrices with strain components
240  // that multiply the Bmat_nl terms
241  mat_x.row(0) = ddisp_dx;
242  mat_x.row(2) = ddisp_dy;
243 
244  mat_y.row(1) = ddisp_dy;
245  mat_y.row(2) = ddisp_dx;
246  }
247  else
248  Bmat_lin.vector_mult(epsilon, local_disp);
249 }
250 
251 
252 
253 void
254 MAST::StructuralElement2D::
255 initialize_strain_operator_gradient(const unsigned int qp,
256  const MAST::FEBase& fe,
257  const RealVectorX& local_disp,
258  RealMatrixX& epsilon_grad,
259  std::vector<MAST::FEMOperatorMatrix>& dBmat_lin) {
260 
261  epsilon_grad.setZero();
262 
263  const std::vector<std::vector<libMesh::RealTensorValue> >&
264  d2phi = fe.get_d2phi();
265 
266  unsigned int n_phi = (unsigned int)d2phi.size();
267  RealVectorX
268  phi = RealVectorX::Zero(n_phi),
269  depsilon = RealVectorX::Zero(3);
270 
271  // make sure all matrices are the right size
272  libmesh_assert_equal_to(epsilon_grad.rows(), 3);
273  libmesh_assert_equal_to(epsilon_grad.cols(), 2);
274  libmesh_assert_equal_to(dBmat_lin.size(), 2);
275 
276 
277  for (unsigned int i=0; i<2; i++) {
278 
279  libmesh_assert_equal_to(dBmat_lin[i].m(), 3);
280  libmesh_assert_equal_to(dBmat_lin[i].n(), 6*n_phi);
281 
282  // now set the shape function values
283  // d2N/dxdxi
284  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
285  phi(i_nd) = d2phi[i_nd][qp](0, i);
286 
287  // linear strain operator
288  dBmat_lin[i].set_shape_function(0, 0, phi); // epsilon_xx = du/dx
289  dBmat_lin[i].set_shape_function(2, 1, phi); // gamma_xy = dv/dx + ...
290 
291  // d2N/dydxi
292  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
293  phi(i_nd) = d2phi[i_nd][qp](1, i);
294 
295  // linear strain operator
296  dBmat_lin[i].set_shape_function(1, 1, phi); // epsilon_yy = dv/dy
297  dBmat_lin[i].set_shape_function(2, 0, phi); // gamma_xy = du/dy + ...
298 
299  dBmat_lin[i].vector_mult(depsilon, local_disp);
300  epsilon_grad.col(i) = depsilon;
301  }
302 }
303 
304 
305 
306 bool
307 MAST::StructuralElement2D::calculate_stress(bool request_derivative,
308  const MAST::FunctionBase* p,
309  MAST::StressStrainOutputBase& output) {
310 
311  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
312 
313  const unsigned int
314  qp_loc_fe_size = (unsigned int)fe->get_qpoints().size(),
315  n_added_qp = 2;
316 
317  std::vector<libMesh::Point>
318  qp_loc_fe = fe->get_qpoints(),
319  qp_loc(qp_loc_fe.size()*n_added_qp);
320  std::vector<Real> JxW = fe->get_JxW();
321  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
322 
323  // we will evaluate the stress at upper and lower layers of element,
324  // so we will add two new points for each qp_loc.
325  // TODO: this will need to be moved to element property section card
326  // to handle composite elements at a later date
327  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
328 
329  qp_loc[i*2] = qp_loc_fe[i];
330  qp_loc[i*2](2) = +1.; // upper skin
331 
332  qp_loc[i*2+1] = qp_loc_fe[i];
333  qp_loc[i*2+1](2) = -1.; // lower skin
334 
336  // scale the JxW of each QP by (1/n_added_qp) since we are
337  // adding extra points and the total volume should add up to be
338  // the same.
340  JxW[i] /= (1.*n_added_qp);
341  }
342 
343 
344 
345  MAST::BendingOperatorType bending_model =
346  _property.bending_model(_elem);
347 
348  std::unique_ptr<MAST::BendingOperator2D>
349  bend(MAST::build_bending_operator_2D(bending_model,
350  *this,
351  qp_loc_fe).release());
352 
353  // now that the FE object has been initialized, evaluate the stress values
354 
355  const unsigned int
356  n_phi = (unsigned int)fe->n_shape_functions(),
357  n1 = this->n_direct_strain_components(),
358  n2 = 6*n_phi,
359  n3 = this->n_von_karman_strain_components();
360 
361  Real
362  z = 0.,
363  z_off = 0.,
364  temp = 0.,
365  ref_t = 0.,
366  alpha = 0.,
367  dtemp = 0.,
368  dref_t= 0.,
369  dalpha= 0.;
370 
371  RealMatrixX
372  material_mat,
373  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
374  dstrain_dX = RealMatrixX::Zero(n1,n2),
375  dstress_dX = RealMatrixX::Zero(n1,n2),
376  mat_n1n2 = RealMatrixX::Zero(n1,n2),
377  eye = RealMatrixX::Identity(n1, n1),
378  mat_x = RealMatrixX::Zero(3,2),
379  mat_y = RealMatrixX::Zero(3,2),
380  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
381  dstress_dX_3D= RealMatrixX::Zero(6,n2);
382 
383  RealVectorX
384  strain = RealVectorX::Zero(n1),
385  stress = RealVectorX::Zero(n1),
386  strain_vk = RealVectorX::Zero(n1),
387  strain_bend = RealVectorX::Zero(n1),
388  strain_3D = RealVectorX::Zero(6),
389  stress_3D = RealVectorX::Zero(6),
390  dstrain_dp = RealVectorX::Zero(n1),
391  dstress_dp = RealVectorX::Zero(n1),
392  vec1 = RealVectorX::Zero(n2),
393  vec2 = RealVectorX::Zero(n2);
394 
395 
396  FEMOperatorMatrix
397  Bmat_lin,
398  Bmat_nl_x,
399  Bmat_nl_y,
400  Bmat_nl_u,
401  Bmat_nl_v,
402  Bmat_bend,
403  Bmat_vk;
404 
405  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
406  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
407  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
408  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
409  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
410  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
411  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
412 
413  // TODO: remove this const-cast, which may need change in API of
414  // material card
415  const MAST::FieldFunction<RealMatrixX >&
416  mat_stiff =
417  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
418  stiffness_matrix(2);
419 
420  // get the thickness values for the bending strain calculation
421  const MAST::FieldFunction<Real>
422  &h = _property.get<MAST::FieldFunction<Real> >("h"),
423  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
424 
425 
426  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
427  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
428 
429  // a reference to the stress output data structure
430  MAST::StressStrainOutputBase& stress_output =
431  dynamic_cast<MAST::StressStrainOutputBase&>(output);
432 
433  // check to see if the element has any thermal loads specified
434  // The object returns null
435  MAST::BoundaryConditionBase *thermal_load =
436  stress_output.get_thermal_load_for_elem(_elem);
437 
438  const MAST::FieldFunction<Real>
439  *temp_func = nullptr,
440  *ref_temp_func = nullptr,
441  *alpha_func = nullptr;
442 
443  // get pointers to the temperature, if thermal load is specified
444  if (thermal_load) {
445  temp_func =
446  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
447  ref_temp_func =
448  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
449  alpha_func =
450  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
451  }
452 
453 
454 
456  unsigned int
457  qp = 0;
458  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
459  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
460  {
461  qp = qp_loc_index*n_added_qp + section_qp_index;
462 
463  // get the material matrix
464  mat_stiff(xyz[qp_loc_index], _time, material_mat);
465 
466  this->initialize_green_lagrange_strain_operator(qp_loc_index,
467  *fe,
468  _local_sol,
469  strain,
470  mat_x,
471  mat_y,
472  Bmat_lin,
473  Bmat_nl_x,
474  Bmat_nl_y,
475  Bmat_nl_u,
476  Bmat_nl_v);
477 
478  // if thermal load was specified, then set the thermal strain
479  // component of the total strain
480  if (thermal_load) {
481  (*temp_func) (xyz[qp_loc_index], _time, temp);
482  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
483  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
484  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
485  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
486  }
487 
488  if (if_bending) {
489 
490  // von Karman strain
491  if (if_vk) { // get the vonKarman strain operator if needed
492 
493  this->initialize_von_karman_strain_operator(qp_loc_index,
494  *fe,
495  strain_vk,
496  vk_dwdxi_mat,
497  Bmat_vk);
498  strain += strain_vk;
499  }
500 
501  // add to this the bending strain
502  // TODO: add coupling due to h_offset
503  h (xyz[qp_loc_index], _time, z);
504  h_off(xyz[qp_loc_index], _time, z_off);
505  // TODO: this assumes isotropic section. Multilayered sections need
506  // special considerations
507  bend->initialize_bending_strain_operator_for_z(*fe,
508  qp_loc_index,
509  qp_loc[qp](2) * z/2.+z_off,
510  Bmat_bend);
511  Bmat_bend.vector_mult(strain_bend, _local_sol);
512 
513 
514  // add stress due to bending.
515  strain += strain_bend;
516  }
517 
518  // note that this assumes linear material laws
519  stress = material_mat * strain;
520 
521 
522  // now set the data for the 3D stress-strain vector
523  // this is using only the direct strain/stress.
524  // this can be improved by estimating the shear stresses from
525  // torsion and shear flow from bending.
526  stress_3D(0) = stress(0); // sigma-xx
527  stress_3D(1) = stress(1); // sigma-yy
528  stress_3D(3) = stress(2); // tau-xy
529  strain_3D(0) = strain(0); // epsilon-xx
530  strain_3D(1) = strain(1); // epsilon-yy
531  strain_3D(3) = strain(2); // gamma-xy
532 
533  // set the stress and strain data
534  MAST::StressStrainOutputBase::Data*
535  data = nullptr;
536 
537  // if neither the derivative nor sensitivity is requested, then
538  // we assume that a new data entry is to be provided. Otherwise,
539  // we assume that the stress at this quantity already
540  // exists, and we only need to append sensitivity/derivative
541  // data to it
542  if (!request_derivative && !p)
543  data = &(stress_output.add_stress_strain_at_qp_location(_elem,
544  qp,
545  qp_loc[qp],
546  xyz[qp_loc_index],
547  stress_3D,
548  strain_3D,
549  JxW[qp_loc_index]));
550  else
551  data = &(stress_output.get_stress_strain_data_for_elem_at_qp(_elem, qp));
552 
553 
554  // calculate the derivative if requested
555  if (request_derivative || p) {
556 
557  Bmat_lin.left_multiply(dstrain_dX, eye);
558 
559  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
560 
561  Bmat_nl_x.left_multiply(mat_n1n2, mat_x);
562  dstrain_dX += mat_n1n2;
563  Bmat_nl_y.left_multiply(mat_n1n2, mat_y);
564  dstrain_dX += mat_n1n2;
565  }
566 
567  if (if_bending) {
568 
569  // von Karman strain
570  if (if_vk) {
571 
572  Bmat_vk.left_multiply(mat_n1n2, vk_dwdxi_mat);
573  dstrain_dX += mat_n1n2;
574  }
575 
576  // bending strain
577  Bmat_bend.left_multiply(mat_n1n2, eye);
578  dstrain_dX += mat_n1n2;
579  }
580 
581  // note: this assumes linear material laws
582  dstress_dX = material_mat * dstrain_dX;
583 
584  // copy to the 3D structure
585  // sigma-xx
586  vec1 = dstress_dX.row(0);
587  this->transform_vector_to_global_system(vec1, vec2);
588  dstress_dX_3D.row(0) = vec2;
589  // sigma-yy
590  vec1 = dstress_dX.row(1);
591  this->transform_vector_to_global_system(vec1, vec2);
592  dstress_dX_3D.row(1) = vec2;
593  // tau-xy
594  vec1 = dstress_dX.row(2);
595  this->transform_vector_to_global_system(vec1, vec2);
596  dstress_dX_3D.row(3) = vec2;
597  // epsilon-xx
598  vec1 = dstrain_dX.row(0);
599  this->transform_vector_to_global_system(vec1, vec2);
600  dstrain_dX_3D.row(0) = vec2;
601  // epsilon-yy
602  vec1 = dstrain_dX.row(1);
603  this->transform_vector_to_global_system(vec1, vec2);
604  dstrain_dX_3D.row(1) = vec2;
605  // gamma-xy
606  vec1 = dstrain_dX.row(2);
607  this->transform_vector_to_global_system(vec1, vec2);
608  dstrain_dX_3D.row(3) = vec2;
609 
610  if (request_derivative)
611  data->set_derivatives(dstress_dX_3D, dstrain_dX_3D);
612 
613 
614  if (p) {
615  // sensitivity of the response, s, is
616  // ds/dp = partial s/partial p +
617  // partial s/partial X dX/dp
618  // the first part of the sensitivity is obtained from
619  //
620  // the first term includes direct sensitivity of the stress
621  // with respect to the parameter, while holding the solution
622  // constant. This should include influence of shape changes,
623  // if the parameter is shape-dependent.
624  // TODO: include shape sensitivity.
625  // presently, only material parameter is included
626 
627 
628  dstrain_dp = RealVectorX::Zero(n1);
629 
630  // if thermal load was specified, then set the thermal strain
631  // component of the total strain
632  if (thermal_load) {
633  temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp);
634  ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t);
635  alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha);
636  dstrain_dp(0) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-xx
637  dstrain_dp(1) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-yy
638  }
639 
640 
641 
642  if (if_bending) {
643 
644  // add to this the bending strain
645  h.derivative (*p,
646  xyz[qp_loc_index], _time, z);
647  h_off.derivative(*p,
648  xyz[qp_loc_index], _time, z_off);
649  // TODO: this assumes isotropic section. Multilayered sections need
650  // special considerations
651  bend->initialize_bending_strain_operator_for_z(*fe,
652  qp_loc_index,
653  qp_loc[qp](2) * z/2.+z_off,
654  Bmat_bend);
655  Bmat_bend.vector_mult(strain_bend, _local_sol);
656 
657 
658  // add stress due to bending.
659  dstrain_dp += strain_bend;
660  }
661 
662 
663  // now use this to calculate the stress sensitivity.
664  dstress_dp = material_mat * dstrain_dp;
665 
666  // get the material matrix sensitivity
667  mat_stiff.derivative(*p,
668  xyz[qp_loc_index],
669  _time,
670  material_mat);
671 
672  // partial sensitivity of strain is zero unless it is a
673  // shape parameter.
674  // TODO: shape sensitivity of strain operator
675 
676  // now use this to calculate the stress sensitivity.
677  dstress_dp += material_mat * strain;
678 
679  //
680  // use the derivative data to evaluate the second term in the
681  // sensitivity
682  //
683  dstress_dp += dstress_dX * _local_sol_sens;
684  dstrain_dp += dstrain_dX * _local_sol_sens;
685 
686  // copy the 3D object
687  stress_3D(0) = dstress_dp(0); // sigma-xx
688  stress_3D(1) = dstress_dp(1); // sigma-yy
689  stress_3D(3) = dstress_dp(2); // tau-xy
690  strain_3D(0) = dstrain_dp(0); // epsilon-xx
691  strain_3D(1) = dstrain_dp(1); // epsilon-yy
692  strain_3D(3) = dstrain_dp(2); // gamma-xy
693 
694  // tell the data object about the sensitivity values
695  data->set_sensitivity(*p,
696  stress_3D,
697  strain_3D);
698  }
699  }
700  }
701 
702  // make sure that the number of data points for this element is
703  // the same as the number of requested points
704  libmesh_assert(qp_loc.size() ==
705  stress_output.n_stress_strain_data_for_elem(_elem));
706 
707  // if either derivative or sensitivity was requested, it was provided
708  // by this routine
709  return request_derivative || p;
710 }
711 
712 
713 
714 void
715 MAST::StructuralElement2D::
716 calculate_stress_boundary_velocity(const MAST::FunctionBase& p,
717  MAST::StressStrainOutputBase& output,
718  const unsigned int s,
719  const MAST::FieldFunction<RealVectorX>& vel_f) {
720 
721  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s, true, false));
722 
723  const unsigned int
724  qp_loc_fe_size = (unsigned int)fe->get_qpoints().size(),
725  n_added_qp = 2;
726 
727  std::vector<libMesh::Point>
728  qp_loc_fe = fe->get_qpoints(),
729  qp_loc(qp_loc_fe.size()*n_added_qp);
730  std::vector<Real> JxW_Vn = fe->get_JxW();
731  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
732  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
733 
734  // we will evaluate the stress at upper and lower layers of element,
735  // so we will add two new points for each qp_loc.
736  // TODO: this will need to be moved to element property section card
737  // to handle composite elements at a later date
738  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
739 
740  qp_loc[i*2] = qp_loc_fe[i];
741  qp_loc[i*2](2) = +1.; // upper skin
742 
743  qp_loc[i*2+1] = qp_loc_fe[i];
744  qp_loc[i*2+1](2) = -1.; // lower skin
745  }
746 
747  MAST::BendingOperatorType bending_model =
748  _property.bending_model(_elem);
749 
750  std::unique_ptr<MAST::BendingOperator2D>
751  bend(MAST::build_bending_operator_2D(bending_model,
752  *this,
753  qp_loc_fe).release());
754 
755  const unsigned int
756  n_phi = (unsigned int)fe->n_shape_functions(),
757  n1 = this->n_direct_strain_components(),
758  n2 = 6*n_phi,
759  n3 = this->n_von_karman_strain_components(),
760  dim = 2;
761 
762  Real
763  z = 0.,
764  z_off = 0.,
765  temp = 0.,
766  ref_t = 0.,
767  alpha = 0.,
768  vn = 0.;
769 
770  RealMatrixX
771  material_mat,
772  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
773  mat_n1n2 = RealMatrixX::Zero(n1,n2),
774  mat_x = RealMatrixX::Zero(3,2),
775  mat_y = RealMatrixX::Zero(3,2),
776  eye = RealMatrixX::Identity(n1, n1);
777 
778  RealVectorX
779  strain = RealVectorX::Zero(n1),
780  stress = RealVectorX::Zero(n1),
781  strain_vk = RealVectorX::Zero(n1),
782  strain_bend = RealVectorX::Zero(n1),
783  strain_3D = RealVectorX::Zero(6),
784  stress_3D = RealVectorX::Zero(6),
785  vel = RealVectorX::Zero(dim);
786 
787  // modify the JxW_Vn by multiplying the normal velocity to it
788  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
789 
790  vel_f.derivative(p, xyz[qp], _time, vel);
791  vn = 0.;
792  for (unsigned int i=0; i<dim; i++)
793  vn += vel(i)*face_normals[qp](i);
794  JxW_Vn[qp] *= vn/(1.*n_added_qp);
795  }
796 
797  FEMOperatorMatrix
798  Bmat_lin,
799  Bmat_nl_x,
800  Bmat_nl_y,
801  Bmat_nl_u,
802  Bmat_nl_v,
803  Bmat_bend,
804  Bmat_vk;
805 
806  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
807  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
808  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
809  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
810  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
811  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
812  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
813 
814  // TODO: remove this const-cast, which may need change in API of
815  // material card
816  const MAST::FieldFunction<RealMatrixX >&
817  mat_stiff =
818  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
819  stiffness_matrix(2);
820 
821  // get the thickness values for the bending strain calculation
822  const MAST::FieldFunction<Real>
823  &h = _property.get<MAST::FieldFunction<Real> >("h"),
824  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
825 
826 
827  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
828  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
829 
830  // a reference to the stress output data structure
831  MAST::StressStrainOutputBase& stress_output =
832  dynamic_cast<MAST::StressStrainOutputBase&>(output);
833 
834  // check to see if the element has any thermal loads specified
835  // The object returns null
836  MAST::BoundaryConditionBase *thermal_load =
837  stress_output.get_thermal_load_for_elem(_elem);
838 
839  const MAST::FieldFunction<Real>
840  *temp_func = nullptr,
841  *ref_temp_func = nullptr,
842  *alpha_func = nullptr;
843 
844  // get pointers to the temperature, if thermal load is specified
845  if (thermal_load) {
846  temp_func =
847  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
848  ref_temp_func =
849  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
850  alpha_func =
851  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
852  }
853 
855  unsigned int
856  qp = 0;
857  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
858  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
859  {
860  qp = qp_loc_index*n_added_qp + section_qp_index;
861 
862  // get the material matrix
863  mat_stiff(xyz[qp_loc_index], _time, material_mat);
864 
865  this->initialize_green_lagrange_strain_operator(qp_loc_index,
866  *fe,
867  _local_sol,
868  strain,
869  mat_x,
870  mat_y,
871  Bmat_lin,
872  Bmat_nl_x,
873  Bmat_nl_y,
874  Bmat_nl_u,
875  Bmat_nl_v);
876 
877 
878  // if thermal load was specified, then set the thermal strain
879  // component of the total strain
880  if (thermal_load) {
881  (*temp_func) (xyz[qp_loc_index], _time, temp);
882  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
883  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
884  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
885  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
886  }
887 
888  if (if_bending) {
889 
890  // von Karman strain
891  if (if_vk) { // get the vonKarman strain operator if needed
892 
893  this->initialize_von_karman_strain_operator(qp_loc_index,
894  *fe,
895  strain_vk,
896  vk_dwdxi_mat,
897  Bmat_vk);
898  strain += strain_vk;
899  }
900 
901  // add to this the bending strain
902  // TODO: add coupling due to h_offset
903  h (xyz[qp_loc_index], _time, z);
904  h_off(xyz[qp_loc_index], _time, z_off);
905  // TODO: this assumes isotropic section. Multilayered sections need
906  // special considerations
907  bend->initialize_bending_strain_operator_for_z(*fe,
908  qp_loc_index,
909  qp_loc[qp](2) * z/2.+z_off,
910  Bmat_bend);
911  Bmat_bend.vector_mult(strain_bend, _local_sol);
912 
913 
914  // add stress due to bending.
915  strain += strain_bend;
916  }
917 
918  // note that this assumes linear material laws
919  stress = material_mat * strain;
920 
921 
922  // now set the data for the 3D stress-strain vector
923  // this is using only the direct strain/stress.
924  // this can be improved by estimating the shear stresses from
925  // torsion and shear flow from bending.
926  stress_3D(0) = stress(0); // sigma-xx
927  stress_3D(1) = stress(1); // sigma-yy
928  stress_3D(3) = stress(2); // tau-xy
929  strain_3D(0) = strain(0); // epsilon-xx
930  strain_3D(1) = strain(1); // epsilon-yy
931  strain_3D(3) = strain(2); // gamma-xy
932 
933  // if neither the derivative nor sensitivity is requested, then
934  // we assume that a new data entry is to be provided. Otherwise,
935  // we assume that the stress at this quantity already
936  // exists, and we only need to append sensitivity/derivative
937  // data to it
938  stress_output.add_stress_strain_at_boundary_qp_location(_elem,
939  s,
940  qp,
941  qp_loc[qp],
942  xyz[qp_loc_index],
943  stress_3D,
944  strain_3D,
945  JxW_Vn[qp_loc_index]);
946  }
947 
948  // make sure that the number of data points for this element is
949  // the same as the number of requested points
950  libmesh_assert(qp_loc.size() ==
951  stress_output.n_boundary_stress_strain_data_for_elem(_elem));
952 }
953 
954 
955 
956 void
957 MAST::StructuralElement2D::
958 calculate_stress_temperature_derivative(MAST::FEBase& fe_thermal,
959  MAST::StressStrainOutputBase& output) {
960 
961  std::unique_ptr<MAST::FEBase> fe_str(_elem.init_fe(true, false));
962 
963  // make sure that the structural and thermal FE have the same types and
964  // locations
965  libmesh_assert(fe_str->get_fe_type() == fe_thermal.get_fe_type());
966  // this is a weak assertion. We really want that the same qpoints be used,
967  // but we are assuming that if the elem type is the same, and if the
968  // number fo qpoints is the same, then the location of qpoints will also
969  // be the same.
970  libmesh_assert_equal_to(fe_str->get_qpoints().size(), fe_thermal.get_qpoints().size());
971 
972  MAST::StressStrainOutputBase&
973  stress_output = dynamic_cast<MAST::StressStrainOutputBase&>(output);
974 
975  libmesh_assert(stress_output.get_thermal_load_for_elem(_elem));
976 
977  const unsigned int
978  qp_loc_fe_size = (unsigned int)fe_str->get_qpoints().size(),
979  n_added_qp = 2;
980 
981  std::vector<libMesh::Point>
982  qp_loc_fe = fe_str->get_qpoints(),
983  qp_loc(qp_loc_fe.size()*n_added_qp);
984 
985  const std::vector<std::vector<Real>>& phi_temp = fe_thermal.get_phi();
986  const std::vector<libMesh::Point>& xyz = fe_str->get_xyz();
987 
988  // we will evaluate the stress at upper and lower layers of element,
989  // so we will add two new points for each qp_loc.
990  // TODO: this will need to be moved to element property section card
991  // to handle composite elements at a later date
992  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
993 
994  qp_loc[i*2] = qp_loc_fe[i];
995  qp_loc[i*2](2) = +1.; // upper skin
996 
997  qp_loc[i*2+1] = qp_loc_fe[i];
998  qp_loc[i*2+1](2) = -1.; // lower skin
999  }
1000 
1001 
1002 
1003  MAST::BendingOperatorType bending_model = _property.bending_model(_elem);
1004 
1005  std::unique_ptr<MAST::BendingOperator2D>
1006  bend(MAST::build_bending_operator_2D(bending_model,
1007  *this,
1008  qp_loc_fe).release());
1009 
1010  // now that the FE object has been initialized, evaluate the stress values
1011 
1012  const unsigned int
1013  nt = (unsigned int)fe_thermal.get_phi().size();
1014 
1015  Real
1016  alpha = 0.;
1017 
1018  RealMatrixX
1019  material_mat,
1020  dstress_dX,
1021  dT_mat = RealMatrixX::Zero(3,1),
1022  dstrain_dX_3D = RealMatrixX::Zero(6,nt),
1023  dstress_dX_3D = RealMatrixX::Zero(6,nt),
1024  Bmat_temp = RealMatrixX::Zero(1,nt);
1025 
1026  dT_mat(0) = dT_mat(1) = 1.;
1027 
1028  // TODO: remove this const-cast, which may need change in API of
1029  // material card
1030  const MAST::FieldFunction<RealMatrixX >&
1031  mat_stiff =
1032  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
1033  stiffness_matrix(2);
1034 
1035  const MAST::FieldFunction<Real>
1036  &alpha_func =
1037  _property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion");
1038 
1039 
1041  unsigned int
1042  qp = 0;
1043  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
1044  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
1045  {
1046  qp = qp_loc_index*n_added_qp + section_qp_index;
1047 
1048  // shape function values for the temperature FE
1049  for ( unsigned int i_nd=0; i_nd<nt; i_nd++ )
1050  Bmat_temp(0, i_nd) = phi_temp[i_nd][qp_loc_index];
1051 
1052  // get the material matrix
1053  mat_stiff(xyz[qp_loc_index], _time, material_mat);
1054 
1055  // if thermal load was specified, then set the thermal strain
1056  // component of the total strain
1057  alpha_func(xyz[qp_loc_index], _time, alpha);
1058 
1059  dstress_dX = alpha * material_mat * dT_mat * Bmat_temp;
1060 
1061  dstress_dX_3D.row(0) = -dstress_dX.row(0); // sigma-xx
1062  dstress_dX_3D.row(1) = -dstress_dX.row(1); // sigma-yy
1063  dstress_dX_3D.row(3) = -dstress_dX.row(2); // sigma-xy
1064 
1065  // set the stress and strain data
1066  MAST::StressStrainOutputBase::Data
1067  &data = stress_output.get_stress_strain_data_for_elem_at_qp(_elem, qp);
1068  data.set_derivatives(dstress_dX_3D, dstrain_dX_3D);
1069  }
1070 }
1071 
1072 
1073 
1074 
1075 bool
1076 MAST::StructuralElement2D::internal_residual (bool request_jacobian,
1077  RealVectorX& f,
1078  RealMatrixX& jac)
1079 {
1080  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1081  false,
1082  _property.extra_quadrature_order(_elem)));
1083 
1084  const std::vector<Real>& JxW = fe->get_JxW();
1085  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1086 
1087  const unsigned int
1088  n_phi = (unsigned int)fe->get_phi().size(),
1089  n1 = this->n_direct_strain_components(),
1090  n2 =6*n_phi,
1091  n3 = this->n_von_karman_strain_components();
1092 
1093  RealMatrixX
1094  material_A_mat,
1095  material_B_mat,
1096  material_D_mat,
1097  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1098  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1099  mat3,
1100  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1101  mat5_3n2 = RealMatrixX::Zero(3,n2),
1102  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1103  stress = RealMatrixX::Zero(2,2),
1104  mat_x = RealMatrixX::Zero(3,2),
1105  mat_y = RealMatrixX::Zero(3,2),
1106  local_jac = RealMatrixX::Zero(n2,n2);
1107 
1108  RealVectorX
1109  vec1_n1 = RealVectorX::Zero(n1),
1110  vec2_n1 = RealVectorX::Zero(n1),
1111  vec3_n2 = RealVectorX::Zero(n2),
1112  vec4_n3 = RealVectorX::Zero(n3),
1113  vec5_n3 = RealVectorX::Zero(n3),
1114  vec6_n2 = RealVectorX::Zero(n2),
1115  strain = RealVectorX::Zero(3),
1116  local_f = RealVectorX::Zero(n2);
1117 
1118  FEMOperatorMatrix
1119  Bmat_lin,
1120  Bmat_nl_x,
1121  Bmat_nl_y,
1122  Bmat_nl_u,
1123  Bmat_nl_v,
1124  Bmat_bend,
1125  Bmat_vk;
1126 
1127  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1128  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1129  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1130  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1131  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1132  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1133  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1134 
1135  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1136 
1137  MAST::BendingOperatorType
1138  bending_model = _property.bending_model(_elem);
1139 
1140  std::unique_ptr<MAST::BendingOperator2D> bend;
1141 
1142  if (bending_model != MAST::NO_BENDING)
1143  bend.reset(MAST::build_bending_operator_2D(bending_model,
1144  *this,
1145  fe->get_qpoints()).release());
1146 
1147  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1148  mat_stiff_A = _property.stiffness_A_matrix(*this),
1149  mat_stiff_B = _property.stiffness_B_matrix(*this),
1150  mat_stiff_D = _property.stiffness_D_matrix(*this);
1151 
1152  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1153 
1154  // get the material matrix
1155  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1156 
1157  if (bend.get()) {
1158  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1159  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1160  }
1161 
1162  // now calculte the quantity for these matrices
1163  _internal_residual_operation(if_vk,
1164  n2,
1165  qp,
1166  *fe,
1167  JxW,
1168  request_jacobian,
1169  local_f,
1170  local_jac,
1171  _local_sol,
1172  strain,
1173  bend.get(),
1174  Bmat_lin,
1175  Bmat_nl_x,
1176  Bmat_nl_y,
1177  Bmat_nl_u,
1178  Bmat_nl_v,
1179  Bmat_bend,
1180  Bmat_vk,
1181  mat_x,
1182  mat_y,
1183  stress,
1184  vk_dwdxi_mat,
1185  material_A_mat,
1186  material_B_mat,
1187  material_D_mat,
1188  vec1_n1,
1189  vec2_n1,
1190  vec3_n2,
1191  vec4_n3,
1192  vec5_n3,
1193  vec6_n2,
1194  mat1_n1n2,
1195  mat2_n2n2,
1196  mat3,
1197  mat4_n3n2,
1198  mat5_3n2);
1199  }
1200 
1201 
1202  // now calculate the transverse shear contribution if appropriate for the
1203  // element
1204  if (bend.get() &&
1205  bend->include_transverse_shear_energy())
1206  bend->calculate_transverse_shear_residual(request_jacobian,
1207  local_f,
1208  local_jac);
1209 
1210 
1211  // now transform to the global coorodinate system
1212  transform_vector_to_global_system(local_f, vec3_n2);
1213  f += vec3_n2;
1214 
1215  if (request_jacobian) {
1216  // for 2D elements
1217  if (_elem.dim() == 2) {
1218  // add small values to the diagonal of the theta_z dofs
1219  for (unsigned int i=0; i<n_phi; i++)
1220  local_jac(5*n_phi+i, 5*n_phi+i) = 1.0e-8;
1221  }
1222  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1223  jac += mat2_n2n2;
1224  }
1225 
1226  return request_jacobian;
1227 }
1228 
1229 
1230 
1231 
1232 
1233 bool
1234 MAST::StructuralElement2D::internal_residual_sensitivity (const MAST::FunctionBase& p,
1235  bool request_jacobian,
1236  RealVectorX& f,
1237  RealMatrixX& jac)
1238 {
1239  // this should be true if the function is called
1240  libmesh_assert(!p.is_shape_parameter()); // this is not implemented for now
1241 
1242 
1243  // check if the material property or the provided exterior
1244  // values, like temperature, are functions of the sensitivity parameter
1245  bool calculate = false;
1246  calculate = calculate || _property.depends_on(p);
1247 
1248  // nothing to be calculated if the element does not depend on the
1249  // sensitivity parameter.
1250  if (!calculate)
1251  return false;
1252 
1253  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1254  false,
1255  _property.extra_quadrature_order(_elem)));
1256 
1257  const std::vector<Real>& JxW = fe->get_JxW();
1258  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1259 
1260  const unsigned int
1261  n_phi = (unsigned int)fe->get_phi().size(),
1262  n1 = this->n_direct_strain_components(),
1263  n2 =6*n_phi,
1264  n3 = this->n_von_karman_strain_components();
1265 
1266  RealMatrixX
1267  material_A_mat,
1268  material_B_mat,
1269  material_D_mat,
1270  material_trans_shear_mat,
1271  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1272  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1273  mat3,
1274  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1275  mat5_3n2 = RealMatrixX::Zero(3,n2),
1276  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1277  stress = RealMatrixX::Zero(2,2),
1278  mat_x = RealMatrixX::Zero(3,2),
1279  mat_y = RealMatrixX::Zero(3,2),
1280  local_jac = RealMatrixX::Zero(n2,n2);
1281  RealVectorX
1282  vec1_n1 = RealVectorX::Zero(n1),
1283  vec2_n1 = RealVectorX::Zero(n1),
1284  vec3_n2 = RealVectorX::Zero(n2),
1285  vec4_n3 = RealVectorX::Zero(n3),
1286  vec5_n3 = RealVectorX::Zero(n3),
1287  vec6_n2 = RealVectorX::Zero(n2),
1288  strain = RealVectorX::Zero(3),
1289  local_f = RealVectorX::Zero(n2);
1290 
1291  FEMOperatorMatrix
1292  Bmat_lin,
1293  Bmat_nl_x,
1294  Bmat_nl_y,
1295  Bmat_nl_u,
1296  Bmat_nl_v,
1297  Bmat_bend,
1298  Bmat_vk;
1299 
1300  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1301  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1302  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1303  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1304  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1305  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1306  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1307 
1308  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1309 
1310  MAST::BendingOperatorType
1311  bending_model = _property.bending_model(_elem);
1312 
1313  std::unique_ptr<MAST::BendingOperator2D> bend;
1314 
1315  if (bending_model != MAST::NO_BENDING)
1316  bend.reset(MAST::build_bending_operator_2D(bending_model,
1317  *this,
1318  fe->get_qpoints()).release());
1319 
1320  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1321  mat_stiff_A = _property.stiffness_A_matrix(*this),
1322  mat_stiff_B = _property.stiffness_B_matrix(*this),
1323  mat_stiff_D = _property.stiffness_D_matrix(*this);
1324 
1325  // first calculate the sensitivity due to the parameter
1326  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1327 
1328  // get the material matrix
1329  mat_stiff_A->derivative(p, xyz[qp], _time, material_A_mat);
1330 
1331  if (bend.get()) {
1332 
1333  mat_stiff_B->derivative(p, xyz[qp], _time, material_B_mat);
1334  mat_stiff_D->derivative(p, xyz[qp], _time, material_D_mat);
1335  }
1336 
1337  // now calculte the quantity for these matrices
1338  // this accounts for the sensitivity of the material property matrices
1339  _internal_residual_operation(if_vk,
1340  n2,
1341  qp,
1342  *fe,
1343  JxW,
1344  request_jacobian,
1345  local_f,
1346  local_jac,
1347  _local_sol,
1348  strain,
1349  bend.get(),
1350  Bmat_lin,
1351  Bmat_nl_x,
1352  Bmat_nl_y,
1353  Bmat_nl_u,
1354  Bmat_nl_v,
1355  Bmat_bend,
1356  Bmat_vk,
1357  mat_x,
1358  mat_y,
1359  stress,
1360  vk_dwdxi_mat,
1361  material_A_mat,
1362  material_B_mat,
1363  material_D_mat,
1364  vec1_n1,
1365  vec2_n1,
1366  vec3_n2,
1367  vec4_n3,
1368  vec5_n3,
1369  vec6_n2,
1370  mat1_n1n2,
1371  mat2_n2n2,
1372  mat3,
1373  mat4_n3n2,
1374  mat5_3n2);
1375  }
1376 
1377  // now calculate the transverse shear contribution if appropriate for the
1378  // element
1379  if (bend.get() &&
1380  bend->include_transverse_shear_energy())
1381  bend->calculate_transverse_shear_residual_sensitivity(p,
1382  request_jacobian,
1383  local_f,
1384  local_jac);
1385 
1386  // now transform to the global coorodinate system
1387  transform_vector_to_global_system(local_f, vec3_n2);
1388  f += vec3_n2;
1389 
1390  if (request_jacobian) {
1391  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1392  jac += mat2_n2n2;
1393  }
1394 
1395 
1396  return request_jacobian;
1397 }
1398 
1399 
1400 
1401 void
1402 MAST::StructuralElement2D::
1403 internal_residual_boundary_velocity(const MAST::FunctionBase& p,
1404  const unsigned int s,
1405  const MAST::FieldFunction<RealVectorX>& vel_f,
1406  bool request_jacobian,
1407  RealVectorX& f,
1408  RealMatrixX& jac) {
1409 
1410  // prepare the side finite element
1411  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s,
1412  true,
1413  _property.extra_quadrature_order(_elem)));
1414 
1415  std::vector<Real> JxW_Vn = fe->get_JxW();
1416  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1417  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
1418 
1419  const unsigned int
1420  n_phi = (unsigned int)fe->get_phi().size(),
1421  n1 = this->n_direct_strain_components(),
1422  n2 = 6*n_phi,
1423  n3 = this->n_von_karman_strain_components(),
1424  dim = 2;
1425 
1426  RealMatrixX
1427  material_A_mat,
1428  material_B_mat,
1429  material_D_mat,
1430  material_trans_shear_mat,
1431  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1432  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1433  mat3,
1434  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1435  mat5_3n2 = RealMatrixX::Zero(3,n2),
1436  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1437  stress = RealMatrixX::Zero(2,2),
1438  mat_x = RealMatrixX::Zero(3,2),
1439  mat_y = RealMatrixX::Zero(3,2),
1440  local_jac = RealMatrixX::Zero(n2,n2);
1441  RealVectorX
1442  vec1_n1 = RealVectorX::Zero(n1),
1443  vec2_n1 = RealVectorX::Zero(n1),
1444  vec3_n2 = RealVectorX::Zero(n2),
1445  vec4_n3 = RealVectorX::Zero(n3),
1446  vec5_n3 = RealVectorX::Zero(n3),
1447  vec6_n2 = RealVectorX::Zero(n2),
1448  strain = RealVectorX::Zero(3),
1449  local_f = RealVectorX::Zero(n2),
1450  vel = RealVectorX::Zero(dim);
1451 
1452  FEMOperatorMatrix
1453  Bmat_lin,
1454  Bmat_nl_x,
1455  Bmat_nl_y,
1456  Bmat_nl_u,
1457  Bmat_nl_v,
1458  Bmat_bend,
1459  Bmat_vk;
1460 
1461  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1462  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1463  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1464  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1465  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1466  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1467  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1468 
1469  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1470 
1471  MAST::BendingOperatorType
1472  bending_model = _property.bending_model(_elem);
1473 
1474  std::unique_ptr<MAST::BendingOperator2D> bend;
1475 
1476  if (bending_model != MAST::NO_BENDING)
1477  bend.reset(MAST::build_bending_operator_2D(bending_model,
1478  *this,
1479  fe->get_qpoints()).release());
1480 
1481  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1482  mat_stiff_A = _property.stiffness_A_matrix(*this),
1483  mat_stiff_B = _property.stiffness_B_matrix(*this),
1484  mat_stiff_D = _property.stiffness_D_matrix(*this);
1485 
1486  Real
1487  vn = 0.;
1488 
1489  // modify the JxW_Vn by multiplying the normal velocity to it
1490  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
1491 
1492  vel_f.derivative(p, xyz[qp], _time, vel);
1493  vn = 0.;
1494  for (unsigned int i=0; i<dim; i++)
1495  vn += vel(i)*face_normals[qp](i);
1496  JxW_Vn[qp] *= vn;
1497  }
1498 
1499 
1500  // first calculate the sensitivity due to the parameter
1501  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
1502 
1503  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1504 
1505  if (bend.get()) {
1506 
1507  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1508  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1509  }
1510 
1511  // now calculte the quantity for these matrices
1512  // this accounts for the sensitivity of the material property matrices
1513  _internal_residual_operation(if_vk,
1514  n2,
1515  qp,
1516  *fe,
1517  JxW_Vn,
1518  request_jacobian,
1519  local_f,
1520  local_jac,
1521  _local_sol,
1522  strain,
1523  bend.get(),
1524  Bmat_lin,
1525  Bmat_nl_x,
1526  Bmat_nl_y,
1527  Bmat_nl_u,
1528  Bmat_nl_v,
1529  Bmat_bend,
1530  Bmat_vk,
1531  mat_x,
1532  mat_y,
1533  stress,
1534  vk_dwdxi_mat,
1535  material_A_mat,
1536  material_B_mat,
1537  material_D_mat,
1538  vec1_n1,
1539  vec2_n1,
1540  vec3_n2,
1541  vec4_n3,
1542  vec5_n3,
1543  vec6_n2,
1544  mat1_n1n2,
1545  mat2_n2n2,
1546  mat3,
1547  mat4_n3n2,
1548  mat5_3n2);
1549  }
1550 
1551  // now calculate the transverse shear contribution if appropriate for the
1552  // element
1553  if (bend.get() && bend->include_transverse_shear_energy())
1554  bend->calculate_transverse_shear_residual_boundary_velocity
1555  (p, s, vel_f, request_jacobian, local_f, local_jac);
1556 
1557  // now transform to the global coorodinate system
1558  transform_vector_to_global_system(local_f, vec3_n2);
1559  f += vec3_n2;
1560 
1561  if (request_jacobian) {
1562  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1563  jac += mat2_n2n2;
1564  }
1565 }
1566 
1567 
1568 
1569 
1570 bool
1571 MAST::StructuralElement2D::
1572 internal_residual_jac_dot_state_sensitivity (RealMatrixX& jac) {
1573 
1574  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1575  false,
1576  _property.extra_quadrature_order(_elem)));
1577 
1578  const std::vector<Real>& JxW = fe->get_JxW();
1579  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1580  const unsigned int
1581  n_phi = (unsigned int)fe->get_phi().size(),
1582  n1 = this->n_direct_strain_components(),
1583  n2 = 6*n_phi,
1584  n3 = this->n_von_karman_strain_components();
1585 
1586  RealMatrixX
1587  material_A_mat,
1588  material_B_mat,
1589  material_D_mat,
1590  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1591  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1592  mat3,
1593  vk_dwdxi_mat_sens = RealMatrixX::Zero(n1,n3),
1594  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1595  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1596  stress = RealMatrixX::Zero(2,2),
1597  mat_x = RealMatrixX::Zero(3,2),
1598  mat_y = RealMatrixX::Zero(3,2),
1599  local_jac = RealMatrixX::Zero(n2,n2);
1600  RealVectorX
1601  vec1_n1 = RealVectorX::Zero(n1),
1602  vec2_n1 = RealVectorX::Zero(n1),
1603  strain = RealVectorX::Zero(3);
1604 
1605 
1606  FEMOperatorMatrix
1607  Bmat_lin,
1608  Bmat_nl_x,
1609  Bmat_nl_y,
1610  Bmat_nl_u,
1611  Bmat_nl_v,
1612  Bmat_bend,
1613  Bmat_vk;
1614 
1615  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1616  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1617  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1618  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1619  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1620  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1621  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1622 
1623  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1624 
1625  MAST::BendingOperatorType
1626  bending_model = _property.bending_model(_elem);
1627 
1628  std::unique_ptr<MAST::BendingOperator2D> bend;
1629 
1630  if (bending_model != MAST::NO_BENDING)
1631  bend.reset(MAST::build_bending_operator_2D(bending_model,
1632  *this,
1633  fe->get_qpoints()).release());
1634 
1635  // without the nonlinear strain, this matrix is zero.
1636  if (!if_vk)
1637  return false;
1638 
1639  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1640  mat_stiff_A = _property.stiffness_A_matrix(*this),
1641  mat_stiff_B = _property.stiffness_B_matrix(*this),
1642  mat_stiff_D = _property.stiffness_D_matrix(*this);
1643 
1644 
1645  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1646 
1647  // get the material matrix
1648  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1649  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1650  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1651 
1652  this->initialize_green_lagrange_strain_operator(qp,
1653  *fe,
1654  _local_sol,
1655  strain,
1656  mat_x,
1657  mat_y,
1658  Bmat_lin,
1659  Bmat_nl_x,
1660  Bmat_nl_y,
1661  Bmat_nl_u,
1662  Bmat_nl_v);
1663 
1664  // first handle constant throught the thickness stresses: membrane and vonKarman
1665  Bmat_lin.vector_mult(vec1_n1, _local_sol_sens);
1666  vec2_n1 = material_A_mat * vec1_n1; // linear direct stress
1667 
1668  // get the bending strain operator
1669  if (bend.get()) {
1670  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
1671 
1672  // evaluate the bending stress and add that to the stress vector
1673  // for evaluation in the nonlinear stress term
1674  Bmat_bend.vector_mult(vec1_n1, _local_sol_sens);
1675  vec2_n1 += material_B_mat * vec1_n1;
1676 
1677  if (if_vk) { // get the vonKarman strain operator if needed
1678 
1679  this->initialize_von_karman_strain_operator(qp,
1680  *fe,
1681  vec1_n1, // epsilon_vk
1682  vk_dwdxi_mat,
1683  Bmat_vk);
1684  this->initialize_von_karman_strain_operator_sensitivity(qp,
1685  *fe,
1686  vk_dwdxi_mat_sens);
1687 
1688  // sensitivity of von Karman strain
1689  vec1_n1(0) = (vk_dwdxi_mat(0,0)*vk_dwdxi_mat_sens(0,0)); // dw/dx dwp/dx
1690  vec1_n1(1) = (vk_dwdxi_mat(1,1)*vk_dwdxi_mat_sens(1,1)); // dw/dy dwp/dy
1691  vec1_n1(2) = (vk_dwdxi_mat(0,0)*vk_dwdxi_mat_sens(1,1) + // dw/dx dwp/dy
1692  vk_dwdxi_mat(1,1)*vk_dwdxi_mat_sens(0,0)); // dw/dy dwp/dx
1693 
1694  vec2_n1 += material_A_mat * vec1_n1;
1695  }
1696  }
1697 
1698  // copy the stress values to a matrix
1699  stress(0,0) = vec2_n1(0); // sigma_xx
1700  stress(1,1) = vec2_n1(1); // sigma_yy
1701  stress(0,1) = vec2_n1(2); // gamma_xy
1702  stress(1,0) = vec2_n1(2); // gamma_xy
1703 
1704 
1705  // copy the stress to use here.
1706  vec2_n1.setZero();
1707 
1708  // now calculate the matrix
1709  // vk - membrane: w-displacement with sens
1710  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1711  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1712  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1713  local_jac += JxW[qp] * mat2_n2n2;
1714 
1715  // vk - bending: w-displacement with stress sens
1716  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
1717  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1718  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1719  local_jac += JxW[qp] * mat2_n2n2;
1720 
1721  // vk - vk: with stress sens and stress
1722  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1723  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1724  mat3 = vk_dwdxi_mat_sens.transpose() * material_A_mat * mat3;
1725  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1726  local_jac += JxW[qp] * mat2_n2n2;
1727 
1728  // vk - vk: with stress and stress sens
1729  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1730  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1731  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1732  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1733  local_jac += JxW[qp] * mat2_n2n2;
1734 
1735  // membrane - vk: w-displacement with sens
1736  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1737  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1738  mat3 = material_A_mat * mat3;
1739  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat3);
1740  local_jac += JxW[qp] * mat2_n2n2;
1741 
1742  // bending - vk: w-displacement with stress sens
1743  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1744  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1745  mat3 = material_B_mat.transpose() * mat3;
1746  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat3);
1747  local_jac += JxW[qp] * mat2_n2n2;
1748 
1749  // vk - vk: w-displacement with stress sens
1750  mat3 = RealMatrixX::Zero(2, n2);
1751  Bmat_vk.left_multiply(mat3, stress);
1752  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1753  local_jac += JxW[qp] * mat2_n2n2;
1754  }
1755 
1756  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1757  jac += mat2_n2n2;
1758 
1759  return true;
1760 }
1761 
1762 
1763 
1764 
1765 void
1766 MAST::StructuralElement2D::_internal_residual_operation
1767 (bool if_vk,
1768  const unsigned int n2,
1769  const unsigned int qp,
1770  const MAST::FEBase& fe,
1771  const std::vector<Real>& JxW,
1772  bool request_jacobian,
1773  RealVectorX& local_f,
1774  RealMatrixX& local_jac,
1775  RealVectorX& local_disp,
1776  RealVectorX& strain_mem,
1777  MAST::BendingOperator2D* bend,
1778  FEMOperatorMatrix& Bmat_lin,
1779  FEMOperatorMatrix& Bmat_nl_x,
1780  FEMOperatorMatrix& Bmat_nl_y,
1781  FEMOperatorMatrix& Bmat_nl_u,
1782  FEMOperatorMatrix& Bmat_nl_v,
1783  MAST::FEMOperatorMatrix& Bmat_bend,
1784  MAST::FEMOperatorMatrix& Bmat_vk,
1785  RealMatrixX& mat_x,
1786  RealMatrixX& mat_y,
1787  RealMatrixX& stress,
1788  RealMatrixX& vk_dwdxi_mat,
1789  RealMatrixX& material_A_mat,
1790  RealMatrixX& material_B_mat,
1791  RealMatrixX& material_D_mat,
1792  RealVectorX& vec1_n1,
1793  RealVectorX& vec2_n1,
1794  RealVectorX& vec3_n2,
1795  RealVectorX& vec4_2,
1796  RealVectorX& vec5_2,
1797  RealVectorX& vec6_n2,
1798  RealMatrixX& mat1_n1n2,
1799  RealMatrixX& mat2_n2n2,
1800  RealMatrixX& mat3,
1801  RealMatrixX& mat4_2n2,
1802  RealMatrixX& mat5_3n2) {
1803 
1804  this->initialize_green_lagrange_strain_operator(qp,
1805  fe,
1806  local_disp,
1807  strain_mem,
1808  mat_x,
1809  mat_y,
1810  Bmat_lin,
1811  Bmat_nl_x,
1812  Bmat_nl_y,
1813  Bmat_nl_u,
1814  Bmat_nl_v);
1815 
1816  vec2_n1 = material_A_mat * strain_mem; // membrane stress
1817 
1818  if (bend) {
1819 
1820  // get the bending strain operator
1821  bend->initialize_bending_strain_operator(fe, qp, Bmat_bend);
1822  Bmat_bend.vector_mult(vec1_n1, local_disp);
1823  vec2_n1 += material_B_mat * vec1_n1;
1824 
1825  if (if_vk) { // get the vonKarman strain operator if needed
1826  this->initialize_von_karman_strain_operator(qp,
1827  fe,
1828  vec1_n1, // epsilon_vk
1829  vk_dwdxi_mat,
1830  Bmat_vk);
1831 
1832  strain_mem += vec1_n1; // epsilon_mem + epsilon_vk
1833  vec2_n1 += material_A_mat * vec1_n1; // stress
1834  }
1835  }
1836 
1837 
1838  // copy the stress values to a matrix
1839  stress(0,0) = vec2_n1(0); // sigma_xx
1840  stress(0,1) = vec2_n1(2); // sigma_xy
1841  stress(1,0) = vec2_n1(2); // sigma_yx
1842  stress(1,1) = vec2_n1(1); // sigma_yy
1843 
1844  // now the internal force vector
1845  // this includes the membrane strain operator with all A and B material operators
1846  Bmat_lin.vector_mult_transpose(vec3_n2, vec2_n1);
1847  local_f.topRows(n2) += JxW[qp] * vec3_n2;
1848 
1849  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1850 
1851  // nonlinear strain operator
1852  // x
1853  vec4_2 = mat_x.transpose() * vec2_n1;
1854  Bmat_nl_x.vector_mult_transpose(vec6_n2, vec4_2);
1855  local_f.topRows(n2) += JxW[qp] * vec6_n2;
1856 
1857  // y
1858  vec4_2 = mat_y.transpose() * vec2_n1;
1859  Bmat_nl_y.vector_mult_transpose(vec6_n2, vec4_2);
1860  local_f.topRows(n2) += JxW[qp] * vec6_n2;
1861  }
1862 
1863 
1864  if (bend) {
1865  if (if_vk) {
1866  // von Karman strain
1867  vec4_2 = vk_dwdxi_mat.transpose() * vec2_n1;
1868  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
1869  local_f += JxW[qp] * vec3_n2;
1870  }
1871 
1872  // now coupling with the bending strain
1873  // B_bend^T [B] B_mem
1874  vec1_n1 = material_B_mat.transpose() * strain_mem;
1875  Bmat_bend.vector_mult_transpose(vec3_n2, vec1_n1);
1876  local_f += JxW[qp] * vec3_n2;
1877 
1878  // now bending stress
1879  Bmat_bend.vector_mult(vec2_n1, local_disp);
1880  vec1_n1 = material_D_mat * vec2_n1;
1881  Bmat_bend.vector_mult_transpose(vec3_n2, vec1_n1);
1882  local_f += JxW[qp] * vec3_n2;
1883  }
1884 
1885  if (request_jacobian) {
1886  // membrane - membrane
1887  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1888  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1889  local_jac += JxW[qp] * mat2_n2n2;
1890 
1891  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1892 
1893 
1894  // B_lin^T C mat_x B_x
1895  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1896  mat1_n1n2 = material_A_mat * mat1_n1n2;
1897  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1898  local_jac += JxW[qp] * mat2_n2n2;
1899 
1900  // B_lin^T C mat_y B_y
1901  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1902  mat1_n1n2 = material_A_mat * mat1_n1n2;
1903  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1904  local_jac += JxW[qp] * mat2_n2n2;
1905 
1906  // B_x^T mat_x^T C B_lin
1907  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1908  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1909  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1910  local_jac += JxW[qp] * mat2_n2n2;
1911 
1912  // B_x^T mat_x^T C mat_x B_x
1913  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1914  mat1_n1n2 = material_A_mat * mat1_n1n2;
1915  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1916  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1917  local_jac += JxW[qp] * mat2_n2n2;
1918 
1919  // B_x^T mat_x^T C mat_y B_y
1920  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1921  mat1_n1n2 = material_A_mat * mat1_n1n2;
1922  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1923  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1924  local_jac += JxW[qp] * mat2_n2n2;
1925 
1926  // B_y^T mat_y^T C B_lin
1927  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1928  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1929  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1930  local_jac += JxW[qp] * mat2_n2n2;
1931 
1932  // B_y^T mat_y^T C mat_x B_x
1933  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1934  mat1_n1n2 = material_A_mat * mat1_n1n2;
1935  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1936  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1937  local_jac += JxW[qp] * mat2_n2n2;
1938 
1939  // B_y^T mat_y^T C mat_y B_y
1940  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1941  mat1_n1n2 = material_A_mat * mat1_n1n2;
1942  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1943  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1944  local_jac += JxW[qp] * mat2_n2n2;
1945 
1946  // nonlinear membrane u - nonlinear membrane u
1947  // (geometric / initial stress stiffness)
1948  // u-disp
1949  Bmat_nl_u.left_multiply(mat5_3n2, stress);
1950  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1951  local_jac += JxW[qp] * mat2_n2n2;
1952 
1953  // nonlinear membrane v - nonlinear membrane v
1954  // (geometric / initial stress stiffness)
1955  // v-disp
1956  Bmat_nl_v.left_multiply(mat5_3n2, stress);
1957  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1958  local_jac += JxW[qp] * mat2_n2n2;
1959  }
1960 
1961  if (bend) {
1962  if (if_vk) {
1963  // membrane - vk
1964  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1965  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1966  mat3 = material_A_mat * mat3;
1967  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat3);
1968  local_jac += JxW[qp] * mat2_n2n2;
1969 
1970  // vk - membrane
1971  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1972  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1973  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1974  local_jac += JxW[qp] * mat2_n2n2;
1975 
1976  // vk - vk
1977  mat3 = RealMatrixX::Zero(2, n2);
1978  Bmat_vk.left_multiply(mat3, stress);
1979  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1980  local_jac += JxW[qp] * mat2_n2n2;
1981 
1982  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1983  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1984  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1985  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1986  local_jac += JxW[qp] * mat2_n2n2;
1987 
1988  // bending - vk
1989  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1990  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1991  mat3 = material_B_mat.transpose() * mat3;
1992  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat3);
1993  local_jac += JxW[qp] * mat2_n2n2;
1994 
1995  // vk - bending
1996  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
1997  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1998  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1999  local_jac += JxW[qp] * mat2_n2n2;
2000  }
2001 
2002  // bending - membrane
2003  mat3 = material_B_mat.transpose();
2004  Bmat_lin.left_multiply(mat1_n1n2, mat3);
2005  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
2006  local_jac += JxW[qp] * mat2_n2n2;
2007 
2008  // membrane - bending
2009  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
2010  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
2011  local_jac += JxW[qp] * mat2_n2n2;
2012 
2013  // bending - bending
2014  Bmat_bend.left_multiply(mat1_n1n2, material_D_mat);
2015  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
2016  local_jac += JxW[qp] * mat2_n2n2;
2017  }
2018  }
2019 }
2020 
2021 
2022 
2023 
2024 void
2025 MAST::StructuralElement2D::_convert_prestress_A_mat_to_vector(const RealMatrixX& mat,
2026  RealVectorX& vec) const {
2027 
2028  libmesh_assert_equal_to(mat.rows(), 2);
2029  libmesh_assert_equal_to(mat.cols(), 2);
2030  vec = RealVectorX::Zero(3);
2031  vec(0) = mat(0,0); // sigma x
2032  vec(1) = mat(1,1); // sigma y
2033  vec(2) = mat(0,1); // tau xy
2034 }
2035 
2036 
2037 
2038 void
2039 MAST::StructuralElement2D::_convert_prestress_B_mat_to_vector(const RealMatrixX& mat,
2040  RealVectorX& vec) const {
2041 
2042  libmesh_assert_equal_to(mat.rows(), 2);
2043  libmesh_assert_equal_to(mat.cols(), 2);
2044  vec = RealVectorX::Zero(3);
2045  vec(0) = mat(0,0); // sigma x
2046  vec(1) = mat(1,1); // sigma y
2047  vec(2) = mat(0,1); // tau xy
2048 }
2049 
2050 
2051 
2052 
2053 bool
2054 MAST::StructuralElement2D::prestress_residual (bool request_jacobian,
2055  RealVectorX& f,
2056  RealMatrixX& jac)
2057 {
2058  if (!_property.if_prestressed())
2059  return false;
2060 
2061  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2062  false,
2063  _property.extra_quadrature_order(_elem)));
2064 
2065  const std::vector<Real>& JxW = fe->get_JxW();
2066  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2067  const unsigned int
2068  n_phi = (unsigned int)fe->get_phi().size(),
2069  n1 = this->n_direct_strain_components(),
2070  n2 = 6*n_phi,
2071  n3 = this->n_von_karman_strain_components();
2072 
2073  RealMatrixX
2074  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2075  mat3,
2076  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2077  local_jac = RealMatrixX::Zero(n2,n2),
2078  mat_x = RealMatrixX::Zero(3,2),
2079  mat_y = RealMatrixX::Zero(3,2),
2080  prestress_mat_A,
2081  prestress_mat_B;
2082 
2083  RealVectorX
2084  vec2_n1 = RealVectorX::Zero(n1),
2085  vec3_n2 = RealVectorX::Zero(n2),
2086  vec4_n3 = RealVectorX::Zero(n3),
2087  vec5_n3 = RealVectorX::Zero(n3),
2088  local_f = RealVectorX::Zero(n2),
2089  strain = RealVectorX::Zero(3),
2090  prestress_vec_A,
2091  prestress_vec_B;
2092 
2093  FEMOperatorMatrix
2094  Bmat_lin,
2095  Bmat_nl_x,
2096  Bmat_nl_y,
2097  Bmat_nl_u,
2098  Bmat_nl_v,
2099  Bmat_bend,
2100  Bmat_vk;
2101 
2102  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2103  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2104  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2105  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2106  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2107  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2108  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2109 
2110  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2111 
2112  MAST::BendingOperatorType
2113  bending_model = _property.bending_model(_elem);
2114 
2115  std::unique_ptr<MAST::BendingOperator2D> bend;
2116 
2117  if (bending_model != MAST::NO_BENDING)
2118  bend.reset(MAST::build_bending_operator_2D(bending_model,
2119  *this,
2120  fe->get_qpoints()).release());
2121 
2122  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
2123  prestress_A = _property.prestress_A_matrix(*this),
2124  prestress_B = _property.prestress_B_matrix(*this);
2125 
2126  // now calculate the quantity
2127  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2128 
2129  (*prestress_A)(xyz[qp], _time, prestress_mat_A);
2130  (*prestress_B)(xyz[qp], _time, prestress_mat_B);
2131  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
2132  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
2133 
2134  this->initialize_green_lagrange_strain_operator(qp,
2135  *fe,
2136  _local_sol,
2137  strain,
2138  mat_x,
2139  mat_y,
2140  Bmat_lin,
2141  Bmat_nl_x,
2142  Bmat_nl_y,
2143  Bmat_nl_u,
2144  Bmat_nl_v);
2145 
2146  // get the bending strain operator if needed
2147  vec2_n1.setZero(); // used to store vk strain, if applicable
2148  if (bend.get()) {
2149  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2150 
2151  if (if_vk) // get the vonKarman strain operator if needed
2152  this->initialize_von_karman_strain_operator(qp,
2153  *fe,
2154  vec2_n1,
2155  vk_dwdxi_mat,
2156  Bmat_vk);
2157  }
2158 
2159  // first handle constant throught the thickness stresses: membrane and vonKarman
2160  // multiply this with the constant through the thickness strain
2161  // membrane strain
2162  Bmat_lin.vector_mult_transpose(vec3_n2, prestress_vec_A);
2163  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
2164 
2165  if (bend.get()) {
2166  if (if_vk) {
2167  // von Karman strain
2168  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
2169  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_n3);
2170  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
2171  }
2172 
2173  // now coupling with the bending strain
2174  Bmat_bend.vector_mult_transpose(vec3_n2, prestress_vec_B);
2175  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
2176  }
2177 
2178  if (request_jacobian) {
2179  if (bend.get() && if_vk) {
2180  mat3 = RealMatrixX::Zero(2, n2);
2181  Bmat_vk.left_multiply(mat3, prestress_mat_A);
2182  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2183  local_jac += JxW[qp] * mat2_n2n2;
2184  }
2185  }
2186  }
2187 
2188  // now transform to the global coorodinate system
2189  transform_vector_to_global_system(local_f, vec3_n2);
2190  f += vec3_n2;
2191  if (request_jacobian && if_vk) {
2192  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2193  jac += mat2_n2n2;
2194  }
2195 
2196  // only the nonlinear strain returns a Jacobian for prestressing
2197  return (request_jacobian);
2198 }
2199 
2200 
2201 
2202 
2203 
2204 bool
2205 MAST::StructuralElement2D::prestress_residual_sensitivity (const MAST::FunctionBase& p,
2206  bool request_jacobian,
2207  RealVectorX& f,
2208  RealMatrixX& jac)
2209 {
2210  if (!_property.if_prestressed())
2211  return false;
2212 
2213  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2214  false,
2215  _property.extra_quadrature_order(_elem)));
2216 
2217  const std::vector<Real>& JxW = fe->get_JxW();
2218  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2219  const unsigned int
2220  n_phi = (unsigned int)fe->get_phi().size(),
2221  n1 = this->n_direct_strain_components(),
2222  n2 = 6*n_phi,
2223  n3 = this->n_von_karman_strain_components();
2224 
2225  RealMatrixX
2226  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2227  mat3,
2228  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2229  local_jac = RealMatrixX::Zero(n2,n2),
2230  mat_x = RealMatrixX::Zero(3,2),
2231  mat_y = RealMatrixX::Zero(3,2),
2232  prestress_mat_A,
2233  prestress_mat_B;
2234 
2235  RealVectorX
2236  vec2_n1 = RealVectorX::Zero(n1),
2237  vec3_n2 = RealVectorX::Zero(n2),
2238  vec4_n3 = RealVectorX::Zero(n3),
2239  vec5_n3 = RealVectorX::Zero(n3),
2240  local_f = RealVectorX::Zero(n2),
2241  strain = RealVectorX::Zero(3),
2242  prestress_vec_A,
2243  prestress_vec_B;
2244 
2245  FEMOperatorMatrix
2246  Bmat_lin,
2247  Bmat_nl_x,
2248  Bmat_nl_y,
2249  Bmat_nl_u,
2250  Bmat_nl_v,
2251  Bmat_bend,
2252  Bmat_vk;
2253 
2254  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2255  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2256  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2257  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2258  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2259  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2260  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2261 
2262  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2263 
2264  MAST::BendingOperatorType
2265  bending_model = _property.bending_model(_elem);
2266 
2267  std::unique_ptr<MAST::BendingOperator2D> bend;
2268 
2269  if (bending_model != MAST::NO_BENDING)
2270  bend.reset(MAST::build_bending_operator_2D(bending_model,
2271  *this,
2272  fe->get_qpoints()).release());
2273 
2274  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
2275  prestress_A = _property.prestress_A_matrix(*this),
2276  prestress_B = _property.prestress_B_matrix(*this);
2277 
2278  // transform to the local coordinate system
2279  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2280 
2281  prestress_A->derivative(p, xyz[qp], _time, prestress_mat_A);
2282  prestress_B->derivative(p, xyz[qp], _time, prestress_mat_B);
2283  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
2284  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
2285 
2286  this->initialize_green_lagrange_strain_operator(qp,
2287  *fe,
2288  _local_sol,
2289  strain,
2290  mat_x,
2291  mat_y,
2292  Bmat_lin,
2293  Bmat_nl_x,
2294  Bmat_nl_y,
2295  Bmat_nl_u,
2296  Bmat_nl_v);
2297 
2298  // get the bending strain operator if needed
2299  vec2_n1.setZero(); // used to store vk strain, if applicable
2300  if (bend.get()) {
2301  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2302 
2303  if (if_vk) // get the vonKarman strain operator if needed
2304  this->initialize_von_karman_strain_operator(qp,
2305  *fe,
2306  vec2_n1,
2307  vk_dwdxi_mat,
2308  Bmat_vk);
2309  }
2310 
2311  // first handle constant throught the thickness stresses: membrane and vonKarman
2312  // multiply this with the constant through the thickness strain
2313  // membrane strain
2314  Bmat_lin.vector_mult_transpose(vec3_n2, prestress_vec_A);
2315  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
2316 
2317  if (bend.get()) {
2318  if (if_vk) {
2319  // von Karman strain
2320  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
2321  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_n3);
2322  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
2323  }
2324 
2325  // now coupling with the bending strain
2326  Bmat_bend.vector_mult_transpose(vec3_n2, prestress_vec_B);
2327  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
2328  }
2329 
2330  if (request_jacobian) {
2331  if (bend.get() && if_vk) {
2332  mat3 = RealMatrixX::Zero(2, n2);
2333  Bmat_vk.left_multiply(mat3, prestress_mat_A);
2334  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2335  local_jac += JxW[qp] * mat2_n2n2;
2336  }
2337  }
2338  }
2339 
2340  // now transform to the global coorodinate system
2341  transform_vector_to_global_system(local_f, vec3_n2);
2342  f += vec3_n2;
2343  if (request_jacobian && if_vk) {
2344  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2345  jac += mat2_n2n2;
2346  }
2347 
2348  // only the nonlinear strain returns a Jacobian for prestressing
2349  return (request_jacobian);
2350 }
2351 
2352 
2353 
2354 
2355 
2356 bool
2357 MAST::StructuralElement2D::
2358 surface_traction_residual(bool request_jacobian,
2359  RealVectorX &f,
2360  RealMatrixX &jac,
2361  const unsigned int side,
2362  MAST::BoundaryConditionBase& bc) {
2363 
2364  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2365 
2366  // prepare the side finite element
2367  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false, false));
2368 
2369  const std::vector<Real> &JxW = fe->get_JxW();
2370  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2371  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2372  const unsigned int
2373  n_phi = (unsigned int)phi.size(),
2374  n1 = 3,
2375  n2 = 6*n_phi;
2376 
2377 
2378  // get the function from this boundary condition
2379  const MAST::FieldFunction<RealVectorX>& trac_func =
2380  bc.get<MAST::FieldFunction<RealVectorX> >("traction");
2381 
2382  // get the thickness function to calculate the force
2383  const MAST::FieldFunction<Real>& t_func =
2384  _property.get<MAST::FieldFunction<Real> >("h");
2385 
2386  FEMOperatorMatrix Bmat;
2387  RealVectorX
2388  trac = RealVectorX::Zero(3);
2389  Real
2390  t_val = 0.;
2391 
2392  RealVectorX
2393  phi_vec = RealVectorX::Zero(n_phi),
2394  force = RealVectorX::Zero(2*n1),
2395  local_f = RealVectorX::Zero(n2),
2396  vec_n2 = RealVectorX::Zero(n2);
2397 
2398  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2399 
2400  // now set the shape function values
2401  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2402  phi_vec(i_nd) = phi[i_nd][qp];
2403 
2404  Bmat.reinit(2*n1, phi_vec);
2405 
2406  // get pressure and thickness values
2407  trac_func(qpoint[qp], _time, trac);
2408  t_func (qpoint[qp], _time, t_val);
2409 
2410  // calculate force
2411  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2412  force(i_dim) = trac(i_dim)*t_val;
2413 
2414  Bmat.vector_mult_transpose(vec_n2, force);
2415 
2416  local_f += JxW[qp] * vec_n2;
2417  }
2418 
2419  // now transform to the global system and add
2420  transform_vector_to_global_system(local_f, vec_n2);
2421  f -= vec_n2;
2422 
2423 
2424  return (request_jacobian);
2425 }
2426 
2427 
2428 
2429 bool
2430 MAST::StructuralElement2D::
2431 surface_traction_residual_sensitivity(const MAST::FunctionBase& p,
2432  bool request_jacobian,
2433  RealVectorX &f,
2434  RealMatrixX &jac,
2435  const unsigned int side,
2436  MAST::BoundaryConditionBase& bc) {
2437 
2438  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2439 
2440  // prepare the side finite element
2441  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false, false));
2442 
2443  const std::vector<Real> &JxW = fe->get_JxW();
2444  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2445  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2446  const unsigned int
2447  n_phi = (unsigned int)phi.size(),
2448  n1 = 3,
2449  n2 = 6*n_phi;
2450 
2451 
2452  // get the function from this boundary condition
2453  const MAST::FieldFunction<RealVectorX>& trac_func =
2454  bc.get<MAST::FieldFunction<RealVectorX> >("traction");
2455 
2456  // get the thickness function to calculate the force
2457  const MAST::FieldFunction<Real>& t_func =
2458  _property.get<MAST::FieldFunction<Real> >("h");
2459 
2460  FEMOperatorMatrix Bmat;
2461  RealVectorX
2462  trac = RealVectorX::Zero(3),
2463  dtrac = RealVectorX::Zero(3);
2464  Real
2465  t_val = 0.,
2466  dt_val = 0.;
2467 
2468  RealVectorX
2469  phi_vec = RealVectorX::Zero(n_phi),
2470  force = RealVectorX::Zero(2*n1),
2471  local_f = RealVectorX::Zero(n2),
2472  vec_n2 = RealVectorX::Zero(n2);
2473 
2474  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2475 
2476  // now set the shape function values
2477  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2478  phi_vec(i_nd) = phi[i_nd][qp];
2479 
2480  Bmat.reinit(2*n1, phi_vec);
2481 
2482  // get pressure and thickness values
2483  trac_func(qpoint[qp], _time, trac);
2484  trac_func.derivative(p, qpoint[qp], _time, dtrac);
2485  t_func (qpoint[qp], _time, t_val);
2486  t_func.derivative(p, qpoint[qp], _time, dt_val);
2487 
2488  // calculate force
2489  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2490  force(i_dim) = (trac(i_dim)*dt_val + dtrac(i_dim)*t_val);
2491 
2492  Bmat.vector_mult_transpose(vec_n2, force);
2493 
2494  local_f += JxW[qp] * vec_n2;
2495  }
2496 
2497  // now transform to the global system and add
2498  transform_vector_to_global_system(local_f, vec_n2);
2499  f -= vec_n2;
2500 
2501 
2502  return (request_jacobian);
2503 }
2504 
2505 
2506 
2507 
2508 bool
2509 MAST::StructuralElement2D::
2510 surface_traction_residual_shifted_boundary(bool request_jacobian,
2511  RealVectorX &f,
2512  RealMatrixX &jac,
2513  const unsigned int side,
2514  MAST::BoundaryConditionBase& bc) {
2515 
2516  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2517 
2518  // prepare the side finite element
2519  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, true, true));
2520 
2521  const std::vector<Real> &JxW = fe->get_JxW();
2522  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2523  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2524  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2525  const unsigned int
2526  n_phi = (unsigned int)phi.size(),
2527  n1 = 3,
2528  n2 = 6*n_phi;
2529 
2530 
2531  // get the function from this boundary condition at the projected point
2532  const MAST::FieldFunction<RealVectorX>& trac_func =
2533  bc.get<MAST::FieldFunction<RealVectorX> >("traction");
2534  const MAST::LevelSetBoundaryVelocity& interface =
2535  bc.get<LevelSetBoundaryVelocity>("phi_vel");
2536 
2537  // get the thickness function to calculate the force
2538  const MAST::FieldFunction<Real>& t_func =
2539  _property.get<MAST::FieldFunction<Real> >("h");
2540 
2541  FEMOperatorMatrix Bmat;
2542 
2543  RealVectorX
2544  trac = RealVectorX::Zero(3),
2545  stress = RealVectorX::Zero(3),
2546  dnormal = RealVectorX::Zero(3),
2547  normal = RealVectorX::Zero(3),
2548  bnd_pt = RealVectorX::Zero(3),
2549  phi_vec = RealVectorX::Zero(n_phi),
2550  force = RealVectorX::Zero(2*n1),
2551  local_f = RealVectorX::Zero(n2),
2552  vec_n2 = RealVectorX::Zero(n2);
2553 
2554  RealMatrixX
2555  normal_mat = RealMatrixX::Zero(2,n1),
2556  stress_grad = RealMatrixX::Zero(n1,2),
2557  mat1_n1n2 = RealMatrixX::Zero(2*n1, n2),
2558  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
2559  local_jac = RealMatrixX::Zero(n2, n2),
2560  dstress_dX = RealMatrixX::Zero(n1,n2),
2561  *dstress_mat= request_jacobian?&dstress_dX:nullptr;
2562 
2563  Real
2564  t_val = 0.,
2565  elem_h = _elem.get_reference_elem().hmax();
2566 
2567  libMesh::Point
2568  pt;
2569 
2570  std::vector<RealMatrixX>
2571  stress_grad_dX(2, RealMatrixX::Zero(n1,n2)),
2572  *stress_grad_mat = request_jacobian?&stress_grad_dX:nullptr;
2573 
2574 
2575  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2576 
2577  // now set the shape function values
2578  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2579  phi_vec(i_nd) = phi[i_nd][qp];
2580 
2581  Bmat.reinit(2*n1, phi_vec);
2582 
2583  // first find the projected point corresponding to this boundary location
2584  interface.search_nearest_interface_point(_elem.get_reference_elem(),
2585  side,
2586  qpoint[qp],
2587  _time,
2588  bnd_pt);
2589  pt(0) = bnd_pt(0); pt(1) = bnd_pt(1); pt(2) = bnd_pt(2);
2590  // now evaluate the traction and normal at this projected point
2591  trac_func (pt, _time, trac);
2592  interface.normal_at_point(pt, _time, normal);
2593  // thickness at the quadrature point
2594  t_func (qpoint[qp], _time, t_val);
2595 
2596  // compute difference in normal vector
2597  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2598  dnormal(i_dim) = face_normals[qp](i_dim) - normal(i_dim);
2599  _surface_normal_voigt_notation(dnormal, normal_mat);
2600 
2601  // compute stress
2602  _compute_stress(*fe, qp, stress, dstress_mat);
2603 
2604  // contribution from dnormal
2605  force.topRows(2) = trac.topRows(2) + normal_mat * stress;
2606 
2607  // contribution from gradient of stress, which uses the
2608  // surface normal at the mapped point
2609  _surface_normal_voigt_notation(normal, normal_mat);
2610  _compute_stress_gradient(*fe, qp, stress_grad, stress_grad_mat);
2611 
2612  for (unsigned int i_dim=0; i_dim<2; i_dim++)
2613  force.topRows(2) -= normal_mat * stress_grad.col(i_dim) * (bnd_pt(i_dim) - qpoint[qp](i_dim));
2614 
2615  Bmat.vector_mult_transpose(vec_n2, force);
2616 
2617  /*libMesh::out
2618  << "id: " << _elem.get_reference_elem().id() << std::endl
2619  << "qp: " << qpoint[qp](0) << " " << qpoint[qp](1) << std::endl
2620  << "x-mapped: " << bnd_pt(0) << " " << bnd_pt(1) << std::endl
2621  << "trac " << trac(0) << " " << trac(1) << std::endl
2622  << "normal: " << normal(0) << " " << normal(1) << std::endl
2623  << "force: " << force.transpose() << std::endl << std::endl;*/
2624 
2625  local_f += t_val * JxW[qp] * vec_n2;
2626 
2627  if (request_jacobian) {
2628 
2629  // contribution from the stress
2630  _surface_normal_voigt_notation(dnormal, normal_mat);
2631  mat1_n1n2.topLeftCorner(2, n2) = normal_mat * dstress_dX;
2632 
2633  // contribution from the stress gradient
2634  _surface_normal_voigt_notation(normal, normal_mat);
2635  for (unsigned int i_dim=0; i_dim<2; i_dim++) {
2636 
2637  mat1_n1n2.topLeftCorner(2, n2) -=
2638  (bnd_pt(i_dim) - qpoint[qp](i_dim)) * normal_mat * stress_grad_dX[i_dim];
2639  }
2640 
2641  Bmat.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
2642  local_jac += t_val * JxW[qp] * mat2_n2n2;
2643  }
2644  }
2645 
2646  // now transform to the global system and add
2647  transform_vector_to_global_system(local_f, vec_n2);
2648  f -= vec_n2;
2649 
2650  if (request_jacobian) {
2651  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2652  jac -= mat2_n2n2;
2653  }
2654 
2655 
2656  return request_jacobian;
2657 }
2658 
2659 
2660 
2661 bool
2662 MAST::StructuralElement2D::
2663 surface_traction_residual_shifted_boundary_sensitivity(const MAST::FunctionBase& p,
2664  bool request_jacobian,
2665  RealVectorX &f,
2666  RealMatrixX &jac,
2667  const unsigned int side,
2668  MAST::BoundaryConditionBase& bc) {
2669 
2670  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2671 
2672  // prepare the side finite element
2673  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, true, true));
2674 
2675  const std::vector<Real> &JxW = fe->get_JxW();
2676  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2677  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2678  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2679  const unsigned int
2680  n_phi = (unsigned int)phi.size(),
2681  n1 = 3,
2682  n2 = 6*n_phi;
2683 
2684 
2685  // get the function from this boundary condition at the projected point
2686  const MAST::FieldFunction<RealVectorX>& trac_func =
2687  bc.get<MAST::FieldFunction<RealVectorX> >("traction");
2688  const MAST::LevelSetBoundaryVelocity& interface =
2689  bc.get<LevelSetBoundaryVelocity>("phi_vel");
2690 
2691  // get the thickness function to calculate the force
2692  const MAST::FieldFunction<Real>& t_func =
2693  _property.get<MAST::FieldFunction<Real> >("h");
2694 
2695  FEMOperatorMatrix Bmat;
2696 
2697  RealVectorX
2698  trac = RealVectorX::Zero(3),
2699  dtrac = RealVectorX::Zero(3),
2700  stress = RealVectorX::Zero(3),
2701  dnormal = RealVectorX::Zero(3),
2702  normal = RealVectorX::Zero(3),
2703  normal_sens = RealVectorX::Zero(3),
2704  bnd_pt = RealVectorX::Zero(3),
2705  bnd_pt_sens = RealVectorX::Zero(3),
2706  phi_vec = RealVectorX::Zero(n_phi),
2707  force = RealVectorX::Zero(2*n1),
2708  local_f = RealVectorX::Zero(n2),
2709  vec_n2 = RealVectorX::Zero(n2);
2710 
2711  RealMatrixX
2712  normal_mat = RealMatrixX::Zero(2,n1),
2713  stress_grad = RealMatrixX::Zero(n1,2),
2714  mat1_n1n2 = RealMatrixX::Zero(2*n1, n2),
2715  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
2716  local_jac = RealMatrixX::Zero(n2, n2),
2717  dstress_dX = RealMatrixX::Zero(n1,n2),
2718  *dstress_mat= request_jacobian?&dstress_dX:nullptr;
2719 
2720  Real
2721  t_val = 0.,
2722  dt_val = 0.,
2723  elem_h = _elem.get_reference_elem().hmax();
2724 
2725  libMesh::Point
2726  pt;
2727 
2728  std::vector<RealMatrixX>
2729  stress_grad_dX(2, RealMatrixX::Zero(n1,n2)),
2730  *stress_grad_mat = request_jacobian?&stress_grad_dX:nullptr;
2731 
2732 
2733  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2734 
2735  // now set the shape function values
2736  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2737  phi_vec(i_nd) = phi[i_nd][qp];
2738 
2739  Bmat.reinit(2*n1, phi_vec);
2740 
2741  // first find the projected point corresponding to this boundary location
2742  interface.search_nearest_interface_point(_elem.get_reference_elem(),
2743  side,
2744  qpoint[qp],
2745  _time,
2746  bnd_pt);
2747  interface.search_nearest_interface_point_derivative(p,
2748  _elem.get_reference_elem(),
2749  side,
2750  qpoint[qp],
2751  _time,
2752  bnd_pt_sens);
2753  pt(0) = bnd_pt(0); pt(1) = bnd_pt(1); pt(2) = bnd_pt(2);
2754 
2755  // now evaluate the traction and normal at this projected point
2756  trac_func(pt, _time, trac);
2757  trac_func.derivative(p, pt, _time, dtrac);
2758  interface.normal_at_point(pt, _time, normal);
2759  interface.normal_derivative_at_point(p, pt, _time, normal_sens);
2760  t_func (qpoint[qp], _time, t_val);
2761  t_func.derivative(p, qpoint[qp], _time, dt_val);
2762 
2763  // compute the contribution of normal and stress-sensitivity
2764  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2765  dnormal(i_dim) = face_normals[qp](i_dim) - normal(i_dim);
2766  _surface_normal_voigt_notation(dnormal, normal_mat);
2767  _compute_stress_sensitivity(p, *fe, qp, stress);
2768  force.topRows(2) = dtrac.topRows(2) + normal_mat * stress;
2769 
2770  // add to this the contribution of the normal sensitivity and stress
2771  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2772  dnormal(i_dim) = normal_sens(i_dim);
2773  _surface_normal_voigt_notation(dnormal, normal_mat);
2774  _compute_stress(*fe, qp, stress, nullptr);
2775  force.topRows(2) -= normal_mat * stress;
2776 
2777  // contribution from sensitivity normal with gradient of stress
2778  _compute_stress_gradient(*fe, qp, stress_grad, stress_grad_mat);
2779  for (unsigned int i_dim=0; i_dim<2; i_dim++)
2780  force.topRows(2) -= normal_mat * stress_grad.col(i_dim) * (bnd_pt(i_dim) - qpoint[qp](i_dim));
2781 
2782  // contribution from sensitivity of mapped point
2783  _surface_normal_voigt_notation(normal, normal_mat);
2784  for (unsigned int i_dim=0; i_dim<2; i_dim++)
2785  force.topRows(2) -= normal_mat * stress_grad.col(i_dim) * bnd_pt_sens(i_dim);
2786 
2787  // contribution from normal with sensitivity of gradient of stress
2788  _compute_stress_gradient_sensitivity(p, *fe, qp, stress_grad);
2789  for (unsigned int i_dim=0; i_dim<2; i_dim++)
2790  force.topRows(2) -= normal_mat * stress_grad.col(i_dim) * (bnd_pt(i_dim) - qpoint[qp](i_dim));
2791 
2792 
2793  Bmat.vector_mult_transpose(vec_n2, force);
2794 
2795  local_f += t_val * JxW[qp] * vec_n2;
2796 
2797 
2798  /*libMesh::out
2799  << "id: " << _elem.get_reference_elem().id() << std::endl
2800  << "qp: " << qpoint[qp](0) << " " << qpoint[qp](1) << std::endl
2801  << "x-mapped: " << bnd_pt(0) << " " << bnd_pt(1) << std::endl
2802  << "x-mapped_sens: " << bnd_pt_sens(0) << " " << bnd_pt_sens(1) << std::endl
2803  << "trac " << trac(0) << " " << trac(1) << std::endl
2804  << "trac_sens: " << dtrac(0) << " " << dtrac(1) << std::endl
2805  << "normal: " << normal(0) << " " << normal(1) << std::endl
2806  << "normal_sens: " << normal_sens(0) << " " << normal_sens(1) << std::endl
2807  << "force: " << force.transpose() << std::endl << std::endl;*/
2808 
2809 
2810  // also include contribution from sensitivity of thickness
2811  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2812  dnormal(i_dim) = face_normals[qp](i_dim) - normal(i_dim);
2813  _surface_normal_voigt_notation(dnormal, normal_mat);
2814 
2815  // compute stress
2816  _compute_stress(*fe, qp, stress, dstress_mat);
2817 
2818  // contribution from dnormal
2819  force.topRows(2) = trac.topRows(2) + normal_mat * stress;
2820 
2821  // contribution from gradient of stress, which uses the
2822  // surface normal at the mapped point
2823  _surface_normal_voigt_notation(normal, normal_mat);
2824  _compute_stress_gradient(*fe, qp, stress_grad, nullptr);
2825 
2826  for (unsigned int i_dim=0; i_dim<2; i_dim++)
2827  force.topRows(2) -= normal_mat * stress_grad.col(i_dim) * (bnd_pt(i_dim) - qpoint[qp](i_dim));
2828 
2829  Bmat.vector_mult_transpose(vec_n2, force);
2830 
2831  local_f += dt_val * JxW[qp] * vec_n2;
2832 
2833 
2834  if (request_jacobian) {
2835 
2836  libmesh_assert(false); // to be implemented
2837  }
2838  }
2839 
2840  // now transform to the global system and add
2841  transform_vector_to_global_system(local_f, vec_n2);
2842  f -= vec_n2;
2843 
2844  if (request_jacobian) {
2845  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2846  jac -= mat2_n2n2;
2847  }
2848 
2849 
2850  return request_jacobian;
2851 }
2852 
2853 
2854 
2855 bool
2856 MAST::StructuralElement2D::
2857 surface_pressure_residual(bool request_jacobian,
2858  RealVectorX &f,
2859  RealMatrixX &jac,
2860  const unsigned int side,
2861  MAST::BoundaryConditionBase& bc) {
2862 
2863  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2864 
2865  // prepare the side finite element
2866  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false, false));
2867 
2868  const std::vector<Real> &JxW = fe->get_JxW();
2869  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2870  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2871  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2872  const unsigned int
2873  n_phi = (unsigned int)phi.size(),
2874  n1 = 3,
2875  n2 = 6*n_phi;
2876 
2877 
2878  // get the function from this boundary condition
2879  const MAST::FieldFunction<Real>& p_func =
2880  bc.get<MAST::FieldFunction<Real> >("pressure");
2881 
2882  // get the thickness function to calculate the force
2883  const MAST::FieldFunction<Real>& t_func =
2884  _property.get<MAST::FieldFunction<Real> >("h");
2885 
2886  FEMOperatorMatrix Bmat;
2887  Real
2888  press = 0.,
2889  t_val = 0.;
2890 
2891  RealVectorX
2892  phi_vec = RealVectorX::Zero(n_phi),
2893  force = RealVectorX::Zero(2*n1),
2894  local_f = RealVectorX::Zero(n2),
2895  vec_n2 = RealVectorX::Zero(n2);
2896 
2897  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2898 
2899  // now set the shape function values
2900  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2901  phi_vec(i_nd) = phi[i_nd][qp];
2902 
2903  Bmat.reinit(2*n1, phi_vec);
2904 
2905  // get pressure and thickness values
2906  p_func(qpoint[qp], _time, press);
2907  t_func(qpoint[qp], _time, t_val);
2908 
2909  // calculate force
2910  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2911  force(i_dim) = (press*t_val) * face_normals[qp](i_dim);
2912 
2913  Bmat.vector_mult_transpose(vec_n2, force);
2914 
2915  local_f += JxW[qp] * vec_n2;
2916  }
2917 
2918  // now transform to the global system and add
2919  transform_vector_to_global_system(local_f, vec_n2);
2920  f -= vec_n2;
2921 
2922 
2923  return (request_jacobian);
2924 }
2925 
2926 
2927 
2928 
2929 
2930 bool
2931 MAST::StructuralElement2D::
2932 surface_pressure_residual_sensitivity(const MAST::FunctionBase& p,
2933  bool request_jacobian,
2934  RealVectorX &f,
2935  RealMatrixX &jac,
2936  const unsigned int side,
2937  MAST::BoundaryConditionBase& bc) {
2938 
2939  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2940 
2941  // prepare the side finite element
2942  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false, false));
2943 
2944  const std::vector<Real> &JxW = fe->get_JxW();
2945  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2946  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2947  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2948  const unsigned int
2949  n_phi = (unsigned int)phi.size(),
2950  n1 = 3,
2951  n2 = 6*n_phi;
2952 
2953 
2954  // get the function from this boundary condition
2955  const MAST::FieldFunction<Real>& p_func =
2956  bc.get<MAST::FieldFunction<Real> >("pressure");
2957 
2958  // get the thickness function to calculate the force
2959  const MAST::FieldFunction<Real>& t_func =
2960  _property.get<MAST::FieldFunction<Real> >("h");
2961 
2962 
2963  FEMOperatorMatrix Bmat;
2964  Real
2965  press = 0.,
2966  dpress = 0.,
2967  t_val = 0.,
2968  dt_val = 0.;
2969 
2970  RealVectorX
2971  phi_vec = RealVectorX::Zero(n_phi),
2972  force = RealVectorX::Zero(2*n1),
2973  local_f = RealVectorX::Zero(n2),
2974  vec_n2 = RealVectorX::Zero(n2);
2975 
2976  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2977 
2978  // now set the shape function values
2979  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2980  phi_vec(i_nd) = phi[i_nd][qp];
2981 
2982  Bmat.reinit(2*n1, phi_vec);
2983 
2984  // get pressure and thickness values and their sensitivities
2985  p_func(qpoint[qp], _time, press);
2986  p_func.derivative(p, qpoint[qp], _time, dpress);
2987  t_func(qpoint[qp], _time, t_val);
2988  t_func.derivative(p, qpoint[qp], _time, dt_val);
2989 
2990  // calculate force
2991  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2992  force(i_dim) = (press*dt_val + dpress*t_val) * face_normals[qp](i_dim);
2993 
2994  Bmat.vector_mult_transpose(vec_n2, force);
2995 
2996  local_f += JxW[qp] * vec_n2;
2997  }
2998 
2999  // now transform to the global system and add
3000  transform_vector_to_global_system(local_f, vec_n2);
3001  f -= vec_n2;
3002 
3003 
3004  return (request_jacobian);
3005 }
3006 
3007 
3008 
3009 
3010 
3011 bool
3012 MAST::StructuralElement2D::thermal_residual (bool request_jacobian,
3013  RealVectorX& f,
3014  RealMatrixX& jac,
3015  MAST::BoundaryConditionBase& bc)
3016 {
3017 
3018  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
3019  false,
3020  _property.extra_quadrature_order(_elem)));
3021 
3022  const std::vector<Real>& JxW = fe->get_JxW();
3023  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
3024 
3025  const unsigned int
3026  n_phi = (unsigned int)fe->get_phi().size(),
3027  n1 = this->n_direct_strain_components(),
3028  n2 = 6*n_phi,
3029  n3 = this->n_von_karman_strain_components();
3030 
3031  RealMatrixX
3032  material_exp_A_mat,
3033  material_exp_B_mat,
3034  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
3035  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
3036  mat3 = RealMatrixX::Zero(2, n2),
3037  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
3038  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
3039  stress = RealMatrixX::Zero(2,2),
3040  mat_x = RealMatrixX::Zero(3,2),
3041  mat_y = RealMatrixX::Zero(3,2),
3042  local_jac = RealMatrixX::Zero(n2,n2);
3043 
3044  RealVectorX
3045  vec1_n1 = RealVectorX::Zero(n1),
3046  vec2_n1 = RealVectorX::Zero(n1),
3047  vec3_n2 = RealVectorX::Zero(n2),
3048  vec4_2 = RealVectorX::Zero(2),
3049  vec5_n3 = RealVectorX::Zero(n3),
3050  local_f = RealVectorX::Zero(n2),
3051  strain = RealVectorX::Zero(3),
3052  delta_t = RealVectorX::Zero(1);
3053 
3054  FEMOperatorMatrix
3055  Bmat_lin,
3056  Bmat_nl_x,
3057  Bmat_nl_y,
3058  Bmat_nl_u,
3059  Bmat_nl_v,
3060  Bmat_bend,
3061  Bmat_vk;
3062 
3063  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
3064  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
3065  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
3066  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
3067  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
3068  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
3069  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
3070 
3071  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
3072 
3073  MAST::BendingOperatorType
3074  bending_model = _property.bending_model(_elem);
3075 
3076  std::unique_ptr<MAST::BendingOperator2D> bend;
3077 
3078  if (bending_model != MAST::NO_BENDING)
3079  bend.reset(MAST::build_bending_operator_2D(bending_model,
3080  *this,
3081  fe->get_qpoints()).release());
3082 
3083  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
3084  expansion_A = _property.thermal_expansion_A_matrix(*this),
3085  expansion_B = _property.thermal_expansion_B_matrix(*this);
3086 
3087  const MAST::FieldFunction<Real>
3088  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
3089  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
3090 
3091  Real
3092  t = 0.,
3093  t0 = 0.,
3094  scaling = 1.;
3095 
3096  if (bc.contains("thermal_jacobian_scaling"))
3097  bc.get<MAST::FieldFunction<Real>>("thermal_jacobian_scaling")(scaling);
3098 
3099  for (unsigned int qp=0; qp<JxW.size(); qp++) {
3100 
3101  // this is moved inside the domain since
3102  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
3103  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
3104 
3105  // get the temperature function
3106  temp_func(xyz[qp], _time, t);
3107  ref_temp_func(xyz[qp], _time, t0);
3108  delta_t(0) = t-t0;
3109 
3110  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)} (with membrane strain)
3111  vec2_n1 = material_exp_B_mat * delta_t; // [C]{alpha (T - T0)} (with bending strain)
3112  stress(0,0) = vec1_n1(0); // sigma_xx
3113  stress(0,1) = vec1_n1(2); // sigma_xy
3114  stress(1,0) = vec1_n1(2); // sigma_yx
3115  stress(1,1) = vec1_n1(1); // sigma_yy
3116 
3117  this->initialize_green_lagrange_strain_operator(qp,
3118  *fe,
3119  _local_sol,
3120  strain,
3121  mat_x,
3122  mat_y,
3123  Bmat_lin,
3124  Bmat_nl_x,
3125  Bmat_nl_y,
3126  Bmat_nl_u,
3127  Bmat_nl_v);
3128 
3129  // membrane strain
3130  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
3131  local_f += JxW[qp] * vec3_n2;
3132 
3133  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
3134 
3135  // nonlinear strain operotor
3136  // x
3137  vec4_2 = mat_x.transpose() * vec1_n1;
3138  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
3139  local_f.topRows(n2) += JxW[qp] * vec3_n2;
3140 
3141  // y
3142  vec4_2 = mat_y.transpose() * vec1_n1;
3143  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
3144  local_f.topRows(n2) += JxW[qp] * vec3_n2;
3145  }
3146 
3147  if (bend.get()) {
3148  // bending strain
3149  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
3150  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
3151  local_f += JxW[qp] * vec3_n2;
3152 
3153  // von Karman strain
3154  if (if_vk) {
3155  // get the vonKarman strain operator if needed
3156  this->initialize_von_karman_strain_operator(qp,
3157  *fe,
3158  vec2_n1, // epsilon_vk
3159  vk_dwdxi_mat,
3160  Bmat_vk);
3161  // von Karman strain
3162  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
3163  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
3164  local_f += JxW[qp] * vec3_n2;
3165  }
3166  }
3167 
3168  if (request_jacobian &&
3169  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
3170 
3171  // u-disp
3172  Bmat_nl_u.left_multiply(mat3, stress);
3173  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
3174  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
3175 
3176  // v-disp
3177  Bmat_nl_v.left_multiply(mat3, stress);
3178  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
3179  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
3180  }
3181 
3182  if (request_jacobian && if_vk) {
3183  // vk - vk
3184  mat3 = RealMatrixX::Zero(2, n2);
3185  Bmat_vk.left_multiply(mat3, stress);
3186  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
3187  local_jac += JxW[qp] * mat2_n2n2;
3188  }
3189  }
3190 
3191  // now transform to the global coorodinate system
3192  transform_vector_to_global_system(local_f, vec3_n2);
3193  f -= vec3_n2;
3194  if (request_jacobian && if_vk) {
3195  transform_matrix_to_global_system(local_jac, mat2_n2n2);
3196  jac -= scaling * mat2_n2n2;
3197  }
3198 
3199  // Jacobian contribution from von Karman strain
3200  return request_jacobian;
3201 }
3202 
3203 
3204 
3205 
3206 bool
3207 MAST::StructuralElement2D::
3208 thermal_residual_sensitivity (const MAST::FunctionBase& p,
3209  bool request_jacobian,
3210  RealVectorX& f,
3211  RealMatrixX& jac,
3212  MAST::BoundaryConditionBase& bc)
3213 {
3214  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
3215  false,
3216  _property.extra_quadrature_order(_elem)));
3217 
3218  const std::vector<Real>& JxW = fe->get_JxW();
3219  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
3220 
3221  const unsigned int
3222  n_phi = (unsigned int)fe->get_phi().size(),
3223  n1 = this->n_direct_strain_components(),
3224  n2 = 6*n_phi,
3225  n3 = this->n_von_karman_strain_components();
3226 
3227  RealMatrixX
3228  material_exp_A_mat,
3229  material_exp_B_mat,
3230  material_exp_A_mat_sens,
3231  material_exp_B_mat_sens,
3232  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
3233  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
3234  mat3 = RealMatrixX::Zero(2, n2),
3235  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
3236  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
3237  stress = RealMatrixX::Zero(2,2),
3238  mat_x = RealMatrixX::Zero(3,2),
3239  mat_y = RealMatrixX::Zero(3,2),
3240  local_jac = RealMatrixX::Zero(n2,n2);
3241 
3242  RealVectorX
3243  vec1_n1 = RealVectorX::Zero(n1),
3244  vec2_n1 = RealVectorX::Zero(n1),
3245  vec3_n2 = RealVectorX::Zero(n2),
3246  vec4_2 = RealVectorX::Zero(2),
3247  vec5_n1 = RealVectorX::Zero(n1),
3248  local_f = RealVectorX::Zero(n2),
3249  strain = RealVectorX::Zero(3),
3250  delta_t = RealVectorX::Zero(1),
3251  delta_t_sens = RealVectorX::Zero(1);
3252 
3253  FEMOperatorMatrix
3254  Bmat_lin,
3255  Bmat_nl_x,
3256  Bmat_nl_y,
3257  Bmat_nl_u,
3258  Bmat_nl_v,
3259  Bmat_bend,
3260  Bmat_vk;
3261 
3262  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
3263  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
3264  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
3265  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
3266  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
3267  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
3268  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
3269 
3270  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
3271 
3272  MAST::BendingOperatorType
3273  bending_model = _property.bending_model(_elem);
3274 
3275  std::unique_ptr<MAST::BendingOperator2D> bend;
3276 
3277  if (bending_model != MAST::NO_BENDING)
3278  bend.reset(MAST::build_bending_operator_2D(bending_model,
3279  *this,
3280  fe->get_qpoints()).release());
3281 
3282  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
3283  expansion_A = _property.thermal_expansion_A_matrix(*this),
3284  expansion_B = _property.thermal_expansion_B_matrix(*this);
3285 
3286  // temperature function
3287  const MAST::FieldFunction<Real>
3288  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
3289  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
3290 
3291  Real t, t0, t_sens, t0_sens;
3292 
3293  for (unsigned int qp=0; qp<JxW.size(); qp++) {
3294 
3295  // this is moved inside the domain since
3296  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
3297  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
3298  expansion_A->derivative(p, xyz[qp], _time, material_exp_A_mat_sens);
3299  expansion_B->derivative(p, xyz[qp], _time, material_exp_B_mat_sens);
3300 
3301  // get the temperature function
3302  temp_func(xyz[qp], _time, t);
3303  ref_temp_func(xyz[qp], _time, t0);
3304  delta_t(0) = t-t0;
3305 
3306  // get the temperature function
3307  temp_func(xyz[qp], _time, t);
3308  temp_func.derivative(p, xyz[qp], _time, t_sens);
3309  ref_temp_func(xyz[qp], _time, t0);
3310  ref_temp_func.derivative(p, xyz[qp], _time, t0_sens);
3311  delta_t(0) = t-t0;
3312  delta_t_sens(0) = t_sens-t0_sens;
3313 
3314  // now prepare the membrane force sensitivity
3315  vec1_n1 = material_exp_A_mat * delta_t_sens; // [C]{alpha dT/dp} (with membrane strain)
3316  vec2_n1 = material_exp_A_mat_sens * delta_t; // d([C]alpha)/dp (T - T0)} (with membrane strain)
3317  vec1_n1 += vec2_n1;
3318  stress(0,0) = vec1_n1(0); // sigma_xx
3319  stress(0,1) = vec1_n1(2); // sigma_xy
3320  stress(1,0) = vec1_n1(2); // sigma_yx
3321  stress(1,1) = vec1_n1(1); // sigma_yy
3322 
3323  vec2_n1 = material_exp_B_mat * delta_t_sens; // [C]{alpha dT/dp} (with bending strain)
3324  vec5_n1 = material_exp_B_mat_sens * delta_t; // d([C] alpha)/dp (T - T0) (with bending strain)
3325  vec2_n1 += vec5_n1;
3326 
3327 
3328  this->initialize_green_lagrange_strain_operator(qp,
3329  *fe,
3330  _local_sol,
3331  strain,
3332  mat_x,
3333  mat_y,
3334  Bmat_lin,
3335  Bmat_nl_x,
3336  Bmat_nl_y,
3337  Bmat_nl_u,
3338  Bmat_nl_v);
3339 
3340  // membrane strain
3341  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
3342  local_f += JxW[qp] * vec3_n2;
3343 
3344  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
3345 
3346  // nonlinear strain operotor
3347  // x
3348  vec4_2 = mat_x.transpose() * vec1_n1;
3349  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
3350  local_f.topRows(n2) += JxW[qp] * vec3_n2;
3351 
3352  // y
3353  vec4_2 = mat_y.transpose() * vec1_n1;
3354  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
3355  local_f.topRows(n2) += JxW[qp] * vec3_n2;
3356  }
3357 
3358  if (bend.get()) {
3359  // bending strain
3360  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
3361  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
3362  local_f += JxW[qp] * vec3_n2;
3363 
3364  // von Karman strain
3365  if (if_vk) {
3366  // get the vonKarman strain operator if needed
3367  this->initialize_von_karman_strain_operator(qp,
3368  *fe,
3369  vec2_n1, // epsilon_vk
3370  vk_dwdxi_mat,
3371  Bmat_vk);
3372  // von Karman strain
3373  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
3374  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
3375  local_f += JxW[qp] * vec3_n2;
3376  }
3377  }
3378 
3379  if (request_jacobian &&
3380  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
3381 
3382  // u-disp
3383  Bmat_nl_u.left_multiply(mat3, stress);
3384  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
3385  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
3386 
3387  // v-disp
3388  Bmat_nl_v.left_multiply(mat3, stress);
3389  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
3390  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
3391  }
3392 
3393  if (request_jacobian && if_vk) { // Jacobian only for vk strain
3394 
3395  // vk - vk
3396  mat3 = RealMatrixX::Zero(2, n2);
3397  Bmat_vk.left_multiply(mat3, stress);
3398  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
3399  local_jac += JxW[qp] * mat2_n2n2;
3400  }
3401  }
3402 
3403 
3404  // now transform to the global coorodinate system
3405  transform_vector_to_global_system(local_f, vec3_n2);
3406  f -= vec3_n2;
3407  if (request_jacobian && if_vk) {
3408  transform_matrix_to_global_system(local_jac, mat2_n2n2);
3409  jac -= mat2_n2n2;
3410  }
3411 
3412  // Jacobian contribution from von Karman strain
3413  return request_jacobian;
3414 }
3415 
3416 
3417 
3418 
3419 void
3420 MAST::StructuralElement2D::
3421 thermal_residual_boundary_velocity(const MAST::FunctionBase& p,
3422  const unsigned int s,
3423  const MAST::FieldFunction<RealVectorX>& vel_f,
3424  MAST::BoundaryConditionBase& bc,
3425  bool request_jacobian,
3426  RealVectorX& f,
3427  RealMatrixX& jac) {
3428 
3429  // prepare the side finite element
3430  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s,
3431  true,
3432  _property.extra_quadrature_order(_elem)));
3433 
3434  std::vector<Real> JxW_Vn = fe->get_JxW();
3435  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
3436  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
3437 
3438  const unsigned int
3439  n_phi = (unsigned int)fe->get_phi().size(),
3440  n1 = this->n_direct_strain_components(), n2=6*n_phi,
3441  n3 = this->n_von_karman_strain_components(),
3442  dim = 2;
3443 
3444  RealMatrixX
3445  material_exp_A_mat,
3446  material_exp_B_mat,
3447  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
3448  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
3449  mat3 = RealMatrixX::Zero(2, n2),
3450  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
3451  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
3452  stress = RealMatrixX::Zero(2,2),
3453  mat_x = RealMatrixX::Zero(3,2),
3454  mat_y = RealMatrixX::Zero(3,2),
3455  local_jac = RealMatrixX::Zero(n2,n2);
3456 
3457  RealVectorX
3458  vec1_n1 = RealVectorX::Zero(n1),
3459  vec2_n1 = RealVectorX::Zero(n1),
3460  vec3_n2 = RealVectorX::Zero(n2),
3461  vec4_2 = RealVectorX::Zero(2),
3462  vec5_n3 = RealVectorX::Zero(n3),
3463  local_f = RealVectorX::Zero(n2),
3464  delta_t = RealVectorX::Zero(1),
3465  strain = RealVectorX::Zero(3),
3466  vel = RealVectorX::Zero(dim);
3467 
3468  FEMOperatorMatrix
3469  Bmat_lin,
3470  Bmat_nl_x,
3471  Bmat_nl_y,
3472  Bmat_nl_u,
3473  Bmat_nl_v,
3474  Bmat_bend,
3475  Bmat_vk;
3476 
3477  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
3478  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
3479  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
3480  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
3481  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
3482  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
3483  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
3484 
3485  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
3486 
3487  MAST::BendingOperatorType
3488  bending_model = _property.bending_model(_elem);
3489 
3490  std::unique_ptr<MAST::BendingOperator2D> bend;
3491 
3492  if (bending_model != MAST::NO_BENDING)
3493  bend.reset(MAST::build_bending_operator_2D(bending_model,
3494  *this,
3495  fe->get_qpoints()).release());
3496 
3497  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
3498  expansion_A = _property.thermal_expansion_A_matrix(*this),
3499  expansion_B = _property.thermal_expansion_B_matrix(*this);
3500 
3501  const MAST::FieldFunction<Real>
3502  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
3503  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
3504 
3505  Real
3506  vn = 0.,
3507  t = 0.,
3508  t0 = 0.;
3509 
3510  // modify the JxW_Vn by multiplying the normal velocity to it
3511  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
3512 
3513  vel_f.derivative(p, xyz[qp], _time, vel);
3514  vn = 0.;
3515  for (unsigned int i=0; i<dim; i++)
3516  vn += vel(i)*face_normals[qp](i);
3517  JxW_Vn[qp] *= vn;
3518  }
3519 
3520  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
3521 
3522  // this is moved inside the domain since
3523  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
3524  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
3525 
3526  // get the temperature function
3527  temp_func(xyz[qp], _time, t);
3528  ref_temp_func(xyz[qp], _time, t0);
3529  delta_t(0) = t-t0;
3530 
3531  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)} (with membrane strain)
3532  vec2_n1 = material_exp_B_mat * delta_t; // [C]{alpha (T - T0)} (with bending strain)
3533  stress(0,0) = vec1_n1(0); // sigma_xx
3534  stress(0,1) = vec1_n1(2); // sigma_xy
3535  stress(1,0) = vec1_n1(2); // sigma_yx
3536  stress(1,1) = vec1_n1(1); // sigma_yy
3537 
3538  this->initialize_green_lagrange_strain_operator(qp,
3539  *fe,
3540  _local_sol,
3541  strain,
3542  mat_x,
3543  mat_y,
3544  Bmat_lin,
3545  Bmat_nl_x,
3546  Bmat_nl_y,
3547  Bmat_nl_u,
3548  Bmat_nl_v);
3549 
3550  // membrane strain
3551  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
3552  local_f += JxW_Vn[qp] * vec3_n2;
3553 
3554  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
3555 
3556  // nonlinear strain operotor
3557  // x
3558  vec4_2 = mat_x.transpose() * vec1_n1;
3559  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
3560  local_f.topRows(n2) += JxW_Vn[qp] * vec3_n2;
3561 
3562  // y
3563  vec4_2 = mat_y.transpose() * vec1_n1;
3564  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
3565  local_f.topRows(n2) += JxW_Vn[qp] * vec3_n2;
3566  }
3567 
3568  if (bend.get()) {
3569  // bending strain
3570  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
3571  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
3572  local_f += JxW_Vn[qp] * vec3_n2;
3573 
3574  // von Karman strain
3575  if (if_vk) {
3576  // get the vonKarman strain operator if needed
3577  this->initialize_von_karman_strain_operator(qp,
3578  *fe,
3579  vec2_n1, // epsilon_vk
3580  vk_dwdxi_mat,
3581  Bmat_vk);
3582  // von Karman strain
3583  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
3584  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
3585  local_f += JxW_Vn[qp] * vec3_n2;
3586  }
3587  }
3588 
3589  if (request_jacobian &&
3590  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
3591 
3592  // u-disp
3593  Bmat_nl_u.left_multiply(mat3, stress);
3594  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
3595  local_jac.topLeftCorner(n2, n2) += JxW_Vn[qp] * mat2_n2n2;
3596 
3597  // v-disp
3598  Bmat_nl_v.left_multiply(mat3, stress);
3599  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
3600  local_jac.topLeftCorner(n2, n2) += JxW_Vn[qp] * mat2_n2n2;
3601  }
3602 
3603  if (request_jacobian && if_vk) { // Jacobian only for vk strain
3604  // vk - vk
3605  mat3 = RealMatrixX::Zero(2, n2);
3606  Bmat_vk.left_multiply(mat3, stress);
3607  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
3608  local_jac += JxW_Vn[qp] * mat2_n2n2;
3609  }
3610  }
3611 
3612  // now transform to the global coorodinate system
3613  transform_vector_to_global_system(local_f, vec3_n2);
3614  f -= vec3_n2;
3615  if (request_jacobian && if_vk) {
3616  transform_matrix_to_global_system(local_jac, mat2_n2n2);
3617  jac -= mat2_n2n2;
3618  }
3619 }
3620 
3621 
3622 
3623 void
3624 MAST::StructuralElement2D::
3625 thermal_residual_temperature_derivative(const MAST::FEBase& fe_thermal,
3626  RealMatrixX& m) {
3627 
3628 
3629  const std::vector<std::vector<Real>>& phi_temp = fe_thermal.get_phi();
3630  const std::vector<Real>& JxW = fe_thermal.get_JxW();
3631  const std::vector<libMesh::Point>& xyz = fe_thermal.get_xyz();
3632 
3633  std::unique_ptr<MAST::FEBase> fe_str(_elem.init_fe(true, false,
3634  _property.extra_quadrature_order(_elem)));
3635  libmesh_assert(fe_str->get_fe_type() == fe_thermal.get_fe_type());
3636  // this is a weak assertion. We really want that the same qpoints be used,
3637  // but we are assuming that if the elem type is the same, and if the
3638  // number fo qpoints is the same, then the location of qpoints will also
3639  // be the same.
3640  libmesh_assert_equal_to(fe_str->get_qpoints().size(), fe_thermal.get_qpoints().size());
3641 
3642 
3643  const unsigned int
3644  n_phi_str = (unsigned int)fe_str->get_phi().size(),
3645  nt = (unsigned int)fe_thermal.get_phi().size(),
3646  n1 = this->n_direct_strain_components(),
3647  n2 = 6*n_phi_str,
3648  n3 = this->n_von_karman_strain_components();
3649 
3650  RealMatrixX
3651  material_exp_A_mat,
3652  material_exp_B_mat,
3653  mat1_n1nt = RealMatrixX::Zero(n1,nt),
3654  mat2_n1nt = RealMatrixX::Zero(n1,nt),
3655  mat3_n2nt = RealMatrixX::Zero(n2,nt),
3656  mat5,
3657  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
3658  stress = RealMatrixX::Zero(2,2),
3659  mat_x = RealMatrixX::Zero(3,2),
3660  mat_y = RealMatrixX::Zero(3,2),
3661  Bmat_temp = RealMatrixX::Zero(1,nt),
3662  local_m = RealMatrixX::Zero(n2,nt);
3663 
3664  RealVectorX
3665  phi = RealVectorX::Zero(nt),
3666  vec_n1 = RealVectorX::Zero(n1),
3667  strain = RealVectorX::Zero(3);
3668 
3669  FEMOperatorMatrix
3670  Bmat_lin,
3671  Bmat_nl_x,
3672  Bmat_nl_y,
3673  Bmat_nl_u,
3674  Bmat_nl_v,
3675  Bmat_bend,
3676  Bmat_vk;
3677 
3678  Bmat_lin.reinit(n1, _system.n_vars(), n_phi_str); // three stress-strain components
3679  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi_str);
3680  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi_str);
3681  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi_str);
3682  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi_str);
3683  Bmat_bend.reinit(n1, _system.n_vars(), n_phi_str);
3684  Bmat_vk.reinit(n3, _system.n_vars(), n_phi_str); // only dw/dx and dw/dy
3685 
3686  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
3687 
3688  MAST::BendingOperatorType
3689  bending_model = _property.bending_model(_elem);
3690 
3691  std::unique_ptr<MAST::BendingOperator2D> bend;
3692 
3693  if (bending_model != MAST::NO_BENDING)
3694  bend.reset(MAST::build_bending_operator_2D(bending_model,
3695  *this,
3696  fe_str->get_qpoints()).release());
3697 
3698  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
3699  expansion_A = _property.thermal_expansion_A_matrix(*this),
3700  expansion_B = _property.thermal_expansion_B_matrix(*this);
3701 
3702  for (unsigned int qp=0; qp<JxW.size(); qp++) {
3703 
3704  // this is moved inside the domain since
3705  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
3706  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
3707 
3708  // shape function values for the temperature FE
3709  for ( unsigned int i_nd=0; i_nd<nt; i_nd++ )
3710  Bmat_temp(0, i_nd) = phi_temp[i_nd][qp];
3711 
3712  this->initialize_green_lagrange_strain_operator(qp,
3713  *fe_str,
3714  _local_sol,
3715  strain,
3716  mat_x,
3717  mat_y,
3718  Bmat_lin,
3719  Bmat_nl_x,
3720  Bmat_nl_y,
3721  Bmat_nl_u,
3722  Bmat_nl_v);
3723 
3724  mat1_n1nt = material_exp_A_mat * Bmat_temp;
3725  mat2_n1nt = material_exp_B_mat * Bmat_temp;
3726  Bmat_lin.right_multiply_transpose(mat3_n2nt, mat1_n1nt);
3727  local_m += JxW[qp] * mat3_n2nt;
3728 
3729 
3730  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
3731 
3732  // nonlinear strain operotor
3733  // x
3734  mat5 = mat_x.transpose() * mat1_n1nt;
3735  Bmat_nl_x.right_multiply_transpose(mat3_n2nt, mat5);
3736  local_m.topRows(n2) += JxW[qp] * mat3_n2nt;
3737 
3738  // y
3739  mat5 = mat_y.transpose() * mat1_n1nt;
3740  Bmat_nl_y.right_multiply_transpose(mat3_n2nt, mat5);
3741  local_m.topRows(n2) += JxW[qp] * mat3_n2nt;
3742  }
3743 
3744  if (bend.get()) {
3745  // bending strain
3746  bend->initialize_bending_strain_operator(*fe_str, qp, Bmat_bend);
3747  Bmat_bend.right_multiply_transpose(mat3_n2nt, mat2_n1nt);
3748  local_m += JxW[qp] * mat3_n2nt;
3749 
3750  // von Karman strain
3751  if (if_vk) {
3752  // get the vonKarman strain operator if needed
3753  this->initialize_von_karman_strain_operator(qp,
3754  *fe_str,
3755  vec_n1, // epsilon_vk
3756  vk_dwdxi_mat,
3757  Bmat_vk);
3758  // von Karman strain
3759  mat5 = vk_dwdxi_mat.transpose() * mat1_n1nt;
3760  Bmat_vk.right_multiply_transpose(mat3_n2nt, mat5);
3761  local_m += JxW[qp] * mat3_n2nt;
3762  }
3763  }
3764  }
3765 
3766  mat3_n2nt.setZero();
3767  phi = RealVectorX::Zero(n2);
3768  vec_n1 = RealVectorX::Zero(n2);
3769  for (unsigned int i=0; i<nt; i++) {
3770  phi = local_m.col(i);
3771  transform_vector_to_global_system(phi, vec_n1);
3772  mat3_n2nt.col(i) = vec_n1;
3773  }
3774  m -= mat3_n2nt;
3775 }
3776 
3777 
3778 
3779 
3780 bool
3781 MAST::StructuralElement2D::
3782 piston_theory_residual(bool request_jacobian,
3783  RealVectorX &f,
3784  RealMatrixX& jac_xdot,
3785  RealMatrixX& jac,
3786  const unsigned int side,
3787  MAST::BoundaryConditionBase& bc) {
3788 
3789  libmesh_assert(false); // to be implemented
3790 
3791  return (request_jacobian);
3792 }
3793 
3794 
3795 
3796 
3797 bool
3798 MAST::StructuralElement2D::
3799 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
3800  bool request_jacobian,
3801  RealVectorX &f,
3802  RealMatrixX& jac_xdot,
3803  RealMatrixX& jac,
3804  const unsigned int side,
3805  MAST::BoundaryConditionBase& bc) {
3806 
3807  libmesh_assert(false); // to be implemented
3808 
3809  return (request_jacobian);
3810 }
3811 
3812 
3813 
3814 
3815 
3816 bool
3817 MAST::StructuralElement2D::
3818 piston_theory_residual(bool request_jacobian,
3819  RealVectorX &f,
3820  RealMatrixX& jac_xdot,
3821  RealMatrixX& jac,
3822  MAST::BoundaryConditionBase& bc) {
3823 
3824  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
3825  libmesh_assert(!follower_forces); // not implemented yet for follower forces
3826 
3827  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
3828 
3829  const std::vector<Real> &JxW = fe->get_JxW();
3830  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
3831  const std::vector<std::vector<Real> >& phi = fe->get_phi();
3832  const unsigned int
3833  n_phi = (unsigned int)phi.size(),
3834  n1 = 2,
3835  n2 = _system.n_vars()*n_phi;
3836 
3837 
3838  // normal for face integration
3839  libMesh::Point normal;
3840  // direction of pressure assumed to be normal (along local z-axis)
3841  // to the element face for 2D and along local y-axis for 1D element.
3842  normal(_elem.dim()) = -1.;
3843 
3844 
3845  // convert to piston theory boundary condition so that the necessary
3846  // flow properties can be obtained
3847  const MAST::PistonTheoryBoundaryCondition& piston_bc =
3848  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
3849 
3850 
3851  // create the constant field functions to pass the dwdx and dwdt values
3852  // to the piston theory pressure functions
3853  MAST::Parameter
3854  dwdx_p ("dwdx", 0.),
3855  dwdt_p ("dwdt", 0.);
3856 
3857  MAST::ConstantFieldFunction
3858  dwdx_f ("dwdx", dwdx_p),
3859  dwdt_f ("dwdx", dwdt_p);
3860 
3861 
3862  std::unique_ptr<MAST::FieldFunction<Real> >
3863  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
3864  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
3865  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
3866 
3867  FEMOperatorMatrix
3868  Bmat_w, // operator matrix for the w-displacement
3869  dBmat; // operator matrix to calculate the derivativ of w wrt x and y
3870 
3871  dBmat.reinit(n1, _system.n_vars(), n_phi);
3872 
3873  RealVectorX
3874  phi_vec = RealVectorX::Zero(n_phi),
3875  force = RealVectorX::Zero(n1),
3876  local_f = RealVectorX::Zero(n2),
3877  vec_n1 = RealVectorX::Zero(n1),
3878  vec_n2 = RealVectorX::Zero(n2),
3879  vel_vec = RealVectorX::Zero(3),
3880  dummy = RealVectorX::Zero(3);
3881 
3882  RealMatrixX
3883  dwdx = RealMatrixX::Zero(3,2),
3884  local_jac_xdot = RealMatrixX::Zero(n2,n2),
3885  local_jac = RealMatrixX::Zero(n2,n2),
3886  mat_n2n2 = RealMatrixX::Zero(n2,n2),
3887  mat_n1n2 = RealMatrixX::Zero(n1,n2),
3888  mat_22 = RealMatrixX::Zero(2,2);
3889 
3890  // we need the velocity vector in the local coordinate system so that
3891  // the appropriate component of the w-derivative can be used
3892  vel_vec = _elem.T_matrix().transpose() * piston_bc.vel_vec();
3893 
3894  Real
3895  dwdt_val = 0.,
3896  dwdx_val = 0.,
3897  p_val = 0.;
3898 
3899 
3900  for (unsigned int qp=0; qp<qpoint.size(); qp++)
3901  {
3902 
3903  // now set the shape function values
3904  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
3905  phi_vec(i_nd) = phi[i_nd][qp];
3906 
3907  // initialize the B matrix for only the w-displacement
3908  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
3909  Bmat_w.set_shape_function(0, 2, phi_vec); // interpolates w-displacement
3910 
3911  // use the Bmat to calculate the velocity vector. Only the w-displacement
3912  // is of interest in the local coordinate, since that is the only
3913  // component normal to the surface.
3914  Bmat_w.right_multiply(vec_n1, _local_vel);
3915  dwdt_val = vec_n1(0);
3916 
3917  // get the operators for dw/dx and dw/dy to calculate the
3918  // normal velocity. We will use the von Karman strain operators
3919  // for this
3920  this->initialize_von_karman_strain_operator(qp,
3921  *fe,
3922  dummy,
3923  dwdx,
3924  dBmat);
3925 
3926  // the diagonal of dwdx matrix stores the
3927  dwdx_val = 0.;
3928  for (unsigned int i=0; i<2; i++)
3929  dwdx_val += dwdx(i,i) * vel_vec(i); // (dw/dx_i)*U_inf . n_i
3930 
3931  // calculate the pressure value
3932  dwdx_p = dwdx_val;
3933  dwdt_p = dwdt_val;
3934  (*pressure)(qpoint[qp], _time, p_val);
3935 
3936  // calculate force
3937  force(0) = p_val * normal(2);
3938 
3939 
3940  Bmat_w.vector_mult_transpose(vec_n2, force);
3941  local_f += JxW[qp] * vec_n2;
3942 
3943 
3944  // calculate the Jacobian if requested
3945  if (request_jacobian) {
3946 
3947  // we need the derivative of cp wrt normal velocity
3948  (*dpressure_dxdot)(qpoint[qp], _time, p_val);
3949 
3950  // calculate the component of Jacobian due to w-velocity
3951  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
3952  local_jac_xdot += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3953 
3954  // now calculate the component of Jacobian
3955  (*dpressure_dx)(qpoint[qp], _time, p_val);
3956 
3957  // derivative wrt x
3958  mat_22.setZero(2,2);
3959  mat_22(0,0) = vel_vec(0);
3960  dBmat.left_multiply(mat_n1n2, mat_22);
3961  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dx
3962  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3963 
3964  // derivative wrt y
3965  mat_22.setZero(2,2);
3966  mat_22(1,1) = vel_vec(1);
3967  dBmat.left_multiply(mat_n1n2, mat_22);
3968  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dy
3969  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3970  }
3971  }
3972 
3973 
3974  // now transform to the global system and add
3975  transform_vector_to_global_system(local_f, vec_n2);
3976  f -= vec_n2;
3977 
3978  // if the Jacobian was requested, then transform it and add to the
3979  // global Jacobian
3980  if (request_jacobian) {
3981  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
3982  jac_xdot -= mat_n2n2;
3983 
3984  transform_matrix_to_global_system(local_jac, mat_n2n2);
3985  jac -= mat_n2n2;
3986  }
3987 
3988  return request_jacobian;
3989 }
3990 
3991 
3992 
3993 
3994 bool
3995 MAST::StructuralElement2D::
3996 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
3997  bool request_jacobian,
3998  RealVectorX &f,
3999  RealMatrixX& jac_xdot,
4000  RealMatrixX& jac,
4001  MAST::BoundaryConditionBase& bc) {
4002 
4003 
4004  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
4005  libmesh_assert(!follower_forces); // not implemented yet for follower forces
4006 
4007  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
4008 
4009  const std::vector<Real> &JxW = fe->get_JxW();
4010  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
4011  const std::vector<std::vector<Real> >& phi = fe->get_phi();
4012  const unsigned int
4013  n_phi = (unsigned int)phi.size(),
4014  n1 = 2,
4015  n2 = _system.n_vars()*n_phi;
4016 
4017 
4018  // normal for face integration
4019  libMesh::Point normal;
4020  // direction of pressure assumed to be normal (along local z-axis)
4021  // to the element face for 2D and along local y-axis for 1D element.
4022  normal(_elem.dim()) = -1.;
4023 
4024 
4025  // convert to piston theory boundary condition so that the necessary
4026  // flow properties can be obtained
4027  const MAST::PistonTheoryBoundaryCondition& piston_bc =
4028  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
4029 
4030 
4031  // create the constant field functions to pass the dwdx and dwdt values
4032  // to the piston theory pressure functions
4033  MAST::Parameter
4034  dwdx_p ("dwdx", 0.),
4035  dwdt_p ("dwdt", 0.);
4036 
4037  MAST::ConstantFieldFunction
4038  dwdx_f ("dwdx", dwdx_p),
4039  dwdt_f ("dwdx", dwdt_p);
4040 
4041 
4042  std::unique_ptr<MAST::FieldFunction<Real> >
4043  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
4044  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
4045  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
4046 
4047  FEMOperatorMatrix
4048  Bmat_w, // operator matrix for the w-displacement
4049  dBmat; // operator matrix to calculate the derivativ of w wrt x and y
4050 
4051  dBmat.reinit(n1, _system.n_vars(), n_phi);
4052 
4053  RealVectorX
4054  phi_vec = RealVectorX::Zero(n_phi),
4055  force = RealVectorX::Zero(n1),
4056  local_f = RealVectorX::Zero(n2),
4057  vec_n1 = RealVectorX::Zero(n1),
4058  vec_n2 = RealVectorX::Zero(n2),
4059  vel_vec = RealVectorX::Zero(3),
4060  dummy = RealVectorX::Zero(3);
4061 
4062  RealMatrixX
4063  dwdx = RealMatrixX::Zero(3,2),
4064  local_jac_xdot = RealMatrixX::Zero(n2,n2),
4065  local_jac = RealMatrixX::Zero(n2,n2),
4066  mat_n2n2 = RealMatrixX::Zero(n2,n2),
4067  mat_n1n2 = RealMatrixX::Zero(n1,n2),
4068  mat_22 = RealMatrixX::Zero(2,2);
4069 
4070  // we need the velocity vector in the local coordinate system so that
4071  // the appropriate component of the w-derivative can be used
4072  vel_vec = _elem.T_matrix().transpose() * piston_bc.vel_vec();
4073 
4074  Real
4075  dwdt_val = 0.,
4076  dwdx_val = 0.,
4077  p_val = 0.;
4078 
4079 
4080  for (unsigned int qp=0; qp<qpoint.size(); qp++)
4081  {
4082 
4083  // now set the shape function values
4084  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
4085  phi_vec(i_nd) = phi[i_nd][qp];
4086 
4087  // initialize the B matrix for only the w-displacement
4088  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
4089  Bmat_w.set_shape_function(0, 2, phi_vec); // interpolates w-displacement
4090 
4091  // use the Bmat to calculate the velocity vector. Only the w-displacement
4092  // is of interest in the local coordinate, since that is the only
4093  // component normal to the surface.
4094  Bmat_w.right_multiply(vec_n1, _local_vel);
4095  dwdt_val = vec_n1(0);
4096 
4097  // get the operators for dw/dx and dw/dy to calculate the
4098  // normal velocity. We will use the von Karman strain operators
4099  // for this
4100  this->initialize_von_karman_strain_operator(qp,
4101  *fe,
4102  dummy,
4103  dwdx,
4104  dBmat);
4105 
4106  // the diagonal of dwdx matrix stores the
4107  dwdx_val = 0.;
4108  for (unsigned int i=0; i<2; i++)
4109  dwdx_val += dwdx(i,i) * vel_vec(i); // (dw/dx_i)*U_inf . n_i
4110 
4111  // calculate the pressure value
4112  dwdx_p = dwdx_val;
4113  dwdt_p = dwdt_val;
4114  pressure->derivative(p, qpoint[qp], _time, p_val);
4115 
4116  // calculate force
4117  force(0) = p_val * normal(2);
4118 
4119 
4120  Bmat_w.vector_mult_transpose(vec_n2, force);
4121  local_f += JxW[qp] * vec_n2;
4122 
4123 
4124  // calculate the Jacobian if requested
4125  if (request_jacobian) {
4126 
4127  // we need the derivative of cp wrt normal velocity
4128  dpressure_dxdot->derivative(p, qpoint[qp], _time, p_val);
4129 
4130  // calculate the component of Jacobian due to w-velocity
4131  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
4132  local_jac_xdot += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
4133 
4134  // now calculate the component of Jacobian
4135  dpressure_dx->derivative(p, qpoint[qp], _time, p_val);
4136 
4137  // derivative wrt x
4138  mat_22.setZero(2,2);
4139  mat_22(0,0) = vel_vec(0);
4140  dBmat.left_multiply(mat_n1n2, mat_22);
4141  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dx
4142  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
4143 
4144  // derivative wrt y
4145  mat_22.setZero(2,2);
4146  mat_22(1,1) = vel_vec(1);
4147  dBmat.left_multiply(mat_n1n2, mat_22);
4148  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dy
4149  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
4150  }
4151  }
4152 
4153 
4154  // now transform to the global system and add
4155  transform_vector_to_global_system(local_f, vec_n2);
4156  f -= vec_n2;
4157 
4158  // if the Jacobian was requested, then transform it and add to the
4159  // global Jacobian
4160  if (request_jacobian) {
4161  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
4162  jac_xdot -= mat_n2n2;
4163 
4164  transform_matrix_to_global_system(local_jac, mat_n2n2);
4165  jac -= mat_n2n2;
4166  }
4167 
4168 
4169  return request_jacobian;
4170 }
4171 
4172 
4173 
4174 
4175 
4176 
4177 void
4178 MAST::StructuralElement2D::_compute_stress(MAST::FEBase& fe,
4179  unsigned int qp,
4180  RealVectorX& stress,
4181  RealMatrixX* dstress_dX) {
4182 
4183  /*
4184  MAST::BendingOperatorType bending_model =
4185  _property.bending_model(_elem);
4186 
4187  std::unique_ptr<MAST::BendingOperator2D>
4188  bend(MAST::build_bending_operator_2D(bending_model,
4189  *this,
4190  qp_loc_fe).release());
4191  */
4192 
4193  std::vector<Real> JxW = fe.get_JxW();
4194  const std::vector<libMesh::Point>& xyz = fe.get_xyz();
4195 
4196  const unsigned int
4197  n_phi = (unsigned int)fe.n_shape_functions(),
4198  n1 = this->n_direct_strain_components(),
4199  n2 = 6*n_phi,
4200  n3 = this->n_von_karman_strain_components();
4201 
4202  /*
4203  Real
4204  z = 0.,
4205  z_off = 0.,
4206  temp = 0.,
4207  ref_t = 0.,
4208  alpha = 0.,
4209  dtemp = 0.,
4210  dref_t= 0.,
4211  dalpha= 0.;
4212  */
4213 
4214  RealMatrixX
4215  material_mat,
4216  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
4217  dstrain_dX = RealMatrixX::Zero(n1,n2),
4218  mat_n1n2 = RealMatrixX::Zero(n1,n2),
4219  eye = RealMatrixX::Identity(n1, n1),
4220  mat_x = RealMatrixX::Zero(3,2),
4221  mat_y = RealMatrixX::Zero(3,2),
4222  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
4223  dstress_dX_3D= RealMatrixX::Zero(6,n2);
4224 
4225  RealVectorX
4226  strain = RealVectorX::Zero(n1),
4227  strain_vk = RealVectorX::Zero(n1),
4228  strain_bend = RealVectorX::Zero(n1),
4229  strain_3D = RealVectorX::Zero(6),
4230  stress_3D = RealVectorX::Zero(6),
4231  dstrain_dp = RealVectorX::Zero(n1),
4232  dstress_dp = RealVectorX::Zero(n1),
4233  vec1 = RealVectorX::Zero(n2),
4234  vec2 = RealVectorX::Zero(n2);
4235 
4236 
4237  FEMOperatorMatrix
4238  Bmat_lin,
4239  Bmat_nl_x,
4240  Bmat_nl_y,
4241  Bmat_nl_u,
4242  Bmat_nl_v,
4243  Bmat_bend,
4244  Bmat_vk;
4245 
4246  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
4247  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
4248  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
4249  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
4250  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
4251  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
4252  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
4253 
4254  // TODO: remove this const-cast, which may need change in API of
4255  // material card
4256  const MAST::FieldFunction<RealMatrixX >&
4257  mat_stiff =
4258  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
4259  stiffness_matrix(2);
4260 
4261  /*
4262  // get the thickness values for the bending strain calculation
4263  const MAST::FieldFunction<Real>
4264  &h = _property.get<MAST::FieldFunction<Real> >("h"),
4265  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
4266 
4267 
4268  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
4269  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
4270 
4271 
4272  // check to see if the element has any thermal loads specified
4273  // The object returns null
4274  MAST::BoundaryConditionBase *thermal_load =
4275  stress_output.get_thermal_load_for_elem(_elem);
4276 
4277  const MAST::FieldFunction<Real>
4278  *temp_func = nullptr,
4279  *ref_temp_func = nullptr,
4280  *alpha_func = nullptr;
4281 
4282  // get pointers to the temperature, if thermal load is specified
4283  if (thermal_load) {
4284  temp_func =
4285  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
4286  ref_temp_func =
4287  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
4288  alpha_func =
4289  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
4290  }
4291  */
4292 
4294  // get the material matrix
4295  mat_stiff(xyz[qp], _time, material_mat);
4296 
4297  this->initialize_green_lagrange_strain_operator(qp,
4298  fe,
4299  _local_sol,
4300  strain,
4301  mat_x,
4302  mat_y,
4303  Bmat_lin,
4304  Bmat_nl_x,
4305  Bmat_nl_y,
4306  Bmat_nl_u,
4307  Bmat_nl_v);
4308 
4309  /*
4310  // if thermal load was specified, then set the thermal strain
4311  // component of the total strain
4312  if (thermal_load) {
4313  (*temp_func) (xyz[qp_loc_index], _time, temp);
4314  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
4315  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
4316  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
4317  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
4318  }
4319 
4320  if (if_bending) {
4321 
4322  // von Karman strain
4323  if (if_vk) { // get the vonKarman strain operator if needed
4324 
4325  this->initialize_von_karman_strain_operator(qp,
4326  fe,
4327  strain_vk,
4328  vk_dwdxi_mat,
4329  Bmat_vk);
4330  strain += strain_vk;
4331  }
4332 
4333  // add to this the bending strain
4334  // TODO: add coupling due to h_offset
4335  h (xyz[qp], _time, z);
4336  h_off(xyz[qp], _time, z_off);
4337  // TODO: this assumes isotropic section. Multilayered sections need
4338  // special considerations
4339  bend->initialize_bending_strain_operator_for_z(fe,
4340  qp,
4341  qp_loc[qp](2) * z/2.+z_off,
4342  Bmat_bend);
4343  Bmat_bend.vector_mult(strain_bend, _local_sol);
4344 
4345 
4346  // add stress due to bending.
4347  strain += strain_bend;
4348  }*/
4349 
4350  // note that this assumes linear material laws
4351  stress = material_mat * strain;
4352 
4353  // calculate the derivative if requested
4354  if (dstress_dX) {
4355 
4356  Bmat_lin.left_multiply(dstrain_dX, eye);
4357 
4358  /*if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
4359 
4360  Bmat_nl_x.left_multiply(mat_n1n2, mat_x);
4361  dstrain_dX += mat_n1n2;
4362  Bmat_nl_y.left_multiply(mat_n1n2, mat_y);
4363  dstrain_dX += mat_n1n2;
4364  }
4365 
4366  if (if_bending) {
4367 
4368  // von Karman strain
4369  if (if_vk) {
4370 
4371  Bmat_vk.left_multiply(mat_n1n2, vk_dwdxi_mat);
4372  dstrain_dX += mat_n1n2;
4373  }
4374 
4375  // bending strain
4376  Bmat_bend.left_multiply(mat_n1n2, eye);
4377  dstrain_dX += mat_n1n2;
4378  }
4379  */
4380 
4381  // note: this assumes linear material laws
4382  *dstress_dX = material_mat * dstrain_dX;
4383 
4384  /*
4385  if (p) {
4386  // sensitivity of the response, s, is
4387  // ds/dp = partial s/partial p +
4388  // partial s/partial X dX/dp
4389  // the first part of the sensitivity is obtained from
4390  //
4391  // the first term includes direct sensitivity of the stress
4392  // with respect to the parameter, while holding the solution
4393  // constant. This should include influence of shape changes,
4394  // if the parameter is shape-dependent.
4395  // TODO: include shape sensitivity.
4396  // presently, only material parameter is included
4397 
4398 
4399  dstrain_dp = RealVectorX::Zero(n1);
4400 
4401  // if thermal load was specified, then set the thermal strain
4402  // component of the total strain
4403  if (thermal_load) {
4404  temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp);
4405  ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t);
4406  alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha);
4407  dstrain_dp(0) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-xx
4408  dstrain_dp(1) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-yy
4409  }
4410 
4411 
4412 
4413  if (if_bending) {
4414 
4415  // add to this the bending strain
4416  h.derivative (*p,
4417  xyz[qp_loc_index], _time, z);
4418  h_off.derivative(*p,
4419  xyz[qp_loc_index], _time, z_off);
4420  // TODO: this assumes isotropic section. Multilayered sections need
4421  // special considerations
4422  bend->initialize_bending_strain_operator_for_z(*fe,
4423  qp_loc_index,
4424  qp_loc[qp](2) * z/2.+z_off,
4425  Bmat_bend);
4426  Bmat_bend.vector_mult(strain_bend, _local_sol);
4427 
4428 
4429  // add stress due to bending.
4430  dstrain_dp += strain_bend;
4431  }
4432 
4433 
4434  // now use this to calculate the stress sensitivity.
4435  dstress_dp = material_mat * dstrain_dp;
4436 
4437  // get the material matrix sensitivity
4438  mat_stiff.derivative(*p,
4439  xyz[qp_loc_index],
4440  _time,
4441  material_mat);
4442 
4443  // partial sensitivity of strain is zero unless it is a
4444  // shape parameter.
4445  // TODO: shape sensitivity of strain operator
4446 
4447  // now use this to calculate the stress sensitivity.
4448  dstress_dp += material_mat * strain;
4449 
4450  //
4451  // use the derivative data to evaluate the second term in the
4452  // sensitivity
4453  //
4454  dstress_dp += dstress_dX * _local_sol_sens;
4455  dstrain_dp += dstrain_dX * _local_sol_sens;
4456 
4457  // copy the 3D object
4458  stress_3D(0) = dstress_dp(0); // sigma-xx
4459  stress_3D(1) = dstress_dp(1); // sigma-yy
4460  stress_3D(3) = dstress_dp(2); // tau-xy
4461  strain_3D(0) = dstrain_dp(0); // epsilon-xx
4462  strain_3D(1) = dstrain_dp(1); // epsilon-yy
4463  strain_3D(3) = dstrain_dp(2); // gamma-xy
4464 
4465  // tell the data object about the sensitivity values
4466  data->set_sensitivity(*p,
4467  stress_3D,
4468  strain_3D);
4469  }
4470  */
4471  }
4472 }
4473 
4474 
4475 
4476 void
4477 MAST::StructuralElement2D::_compute_stress_sensitivity(const MAST::FunctionBase& p,
4478  MAST::FEBase& fe,
4479  unsigned int qp,
4480  RealVectorX& stress) {
4481 
4482  /*
4483  MAST::BendingOperatorType bending_model =
4484  _property.bending_model(_elem);
4485 
4486  std::unique_ptr<MAST::BendingOperator2D>
4487  bend(MAST::build_bending_operator_2D(bending_model,
4488  *this,
4489  qp_loc_fe).release());
4490  */
4491 
4492  std::vector<Real> JxW = fe.get_JxW();
4493  const std::vector<libMesh::Point>& xyz = fe.get_xyz();
4494 
4495  const unsigned int
4496  n_phi = (unsigned int)fe.n_shape_functions(),
4497  n1 = this->n_direct_strain_components(),
4498  n2 = 6*n_phi,
4499  n3 = this->n_von_karman_strain_components();
4500 
4501  /*
4502  Real
4503  z = 0.,
4504  z_off = 0.,
4505  temp = 0.,
4506  ref_t = 0.,
4507  alpha = 0.,
4508  dtemp = 0.,
4509  dref_t= 0.,
4510  dalpha= 0.;
4511  */
4512 
4513  RealMatrixX
4514  material_mat,
4515  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
4516  dstrain_dX = RealMatrixX::Zero(n1,n2),
4517  mat_n1n2 = RealMatrixX::Zero(n1,n2),
4518  eye = RealMatrixX::Identity(n1, n1),
4519  mat_x = RealMatrixX::Zero(3,2),
4520  mat_y = RealMatrixX::Zero(3,2),
4521  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
4522  dstress_dX_3D= RealMatrixX::Zero(6,n2);
4523 
4524  RealVectorX
4525  strain = RealVectorX::Zero(n1),
4526  strain_vk = RealVectorX::Zero(n1),
4527  strain_bend = RealVectorX::Zero(n1),
4528  strain_3D = RealVectorX::Zero(6),
4529  stress_3D = RealVectorX::Zero(6),
4530  dstrain_dp = RealVectorX::Zero(n1),
4531  dstress_dp = RealVectorX::Zero(n1),
4532  vec1 = RealVectorX::Zero(n2),
4533  vec2 = RealVectorX::Zero(n2);
4534 
4535 
4536  FEMOperatorMatrix
4537  Bmat_lin,
4538  Bmat_nl_x,
4539  Bmat_nl_y,
4540  Bmat_nl_u,
4541  Bmat_nl_v,
4542  Bmat_bend,
4543  Bmat_vk;
4544 
4545  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
4546  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
4547  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
4548  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
4549  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
4550  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
4551  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
4552 
4553  // TODO: remove this const-cast, which may need change in API of
4554  // material card
4555  const MAST::FieldFunction<RealMatrixX >&
4556  mat_stiff =
4557  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
4558  stiffness_matrix(2);
4559 
4560  /*
4561  // get the thickness values for the bending strain calculation
4562  const MAST::FieldFunction<Real>
4563  &h = _property.get<MAST::FieldFunction<Real> >("h"),
4564  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
4565 
4566 
4567  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
4568  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
4569 
4570 
4571  // check to see if the element has any thermal loads specified
4572  // The object returns null
4573  MAST::BoundaryConditionBase *thermal_load =
4574  stress_output.get_thermal_load_for_elem(_elem);
4575 
4576  const MAST::FieldFunction<Real>
4577  *temp_func = nullptr,
4578  *ref_temp_func = nullptr,
4579  *alpha_func = nullptr;
4580 
4581  // get pointers to the temperature, if thermal load is specified
4582  if (thermal_load) {
4583  temp_func =
4584  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
4585  ref_temp_func =
4586  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
4587  alpha_func =
4588  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
4589  }
4590  */
4591 
4593  // get the material matrix
4594  mat_stiff.derivative(p, xyz[qp], _time, material_mat);
4595 
4596  this->initialize_green_lagrange_strain_operator(qp,
4597  fe,
4598  _local_sol,
4599  strain,
4600  mat_x,
4601  mat_y,
4602  Bmat_lin,
4603  Bmat_nl_x,
4604  Bmat_nl_y,
4605  Bmat_nl_u,
4606  Bmat_nl_v);
4607 
4608  /*
4609  // if thermal load was specified, then set the thermal strain
4610  // component of the total strain
4611  if (thermal_load) {
4612  (*temp_func) (xyz[qp_loc_index], _time, temp);
4613  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
4614  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
4615  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
4616  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
4617  }
4618 
4619  if (if_bending) {
4620 
4621  // von Karman strain
4622  if (if_vk) { // get the vonKarman strain operator if needed
4623 
4624  this->initialize_von_karman_strain_operator(qp,
4625  fe,
4626  strain_vk,
4627  vk_dwdxi_mat,
4628  Bmat_vk);
4629  strain += strain_vk;
4630  }
4631 
4632  // add to this the bending strain
4633  // TODO: add coupling due to h_offset
4634  h (xyz[qp], _time, z);
4635  h_off(xyz[qp], _time, z_off);
4636  // TODO: this assumes isotropic section. Multilayered sections need
4637  // special considerations
4638  bend->initialize_bending_strain_operator_for_z(fe,
4639  qp,
4640  qp_loc[qp](2) * z/2.+z_off,
4641  Bmat_bend);
4642  Bmat_bend.vector_mult(strain_bend, _local_sol);
4643 
4644 
4645  // add stress due to bending.
4646  strain += strain_bend;
4647  }*/
4648 
4649  // note that this assumes linear material laws
4650  stress = material_mat * strain;
4651 }
4652 
4653 
4654 
4655 void
4656 MAST::StructuralElement2D::
4657 _compute_stress_gradient(MAST::FEBase& fe,
4658  unsigned int qp,
4659  RealMatrixX& stress_grad,
4660  std::vector<RealMatrixX>* dstress_grad_dX) {
4661 
4662  /*
4663  MAST::BendingOperatorType bending_model =
4664  _property.bending_model(_elem);
4665 
4666  std::unique_ptr<MAST::BendingOperator2D>
4667  bend(MAST::build_bending_operator_2D(bending_model,
4668  *this,
4669  qp_loc_fe).release());
4670  */
4671 
4672  std::vector<Real> JxW = fe.get_JxW();
4673  const std::vector<libMesh::Point>& xyz = fe.get_xyz();
4674 
4675  const unsigned int
4676  n_phi = (unsigned int)fe.n_shape_functions(),
4677  n1 = this->n_direct_strain_components(),
4678  n2 = 6*n_phi,
4679  n3 = this->n_von_karman_strain_components();
4680 
4681  /*
4682  Real
4683  z = 0.,
4684  z_off = 0.,
4685  temp = 0.,
4686  ref_t = 0.,
4687  alpha = 0.,
4688  dtemp = 0.,
4689  dref_t= 0.,
4690  dalpha= 0.;
4691  */
4692 
4693  RealMatrixX
4694  material_mat,
4695  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
4696  dstrain = RealMatrixX::Zero(n1,2),
4697  dstrain_dX = RealMatrixX::Zero(n1,n2),
4698  mat_n1n2 = RealMatrixX::Zero(n1,n2),
4699  eye = RealMatrixX::Identity(n1, n1),
4700  mat_x = RealMatrixX::Zero(3,2),
4701  mat_y = RealMatrixX::Zero(3,2),
4702  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
4703  dstress_dX_3D= RealMatrixX::Zero(6,n2);
4704 
4705  RealVectorX
4706  strain_vk = RealVectorX::Zero(n1),
4707  strain_bend = RealVectorX::Zero(n1),
4708  strain_3D = RealVectorX::Zero(6),
4709  stress_3D = RealVectorX::Zero(6),
4710  dstrain_dp = RealVectorX::Zero(n1),
4711  dstress_dp = RealVectorX::Zero(n1),
4712  vec1 = RealVectorX::Zero(n2),
4713  vec2 = RealVectorX::Zero(n2);
4714 
4715 
4716  FEMOperatorMatrix
4717  Bmat_nl_x,
4718  Bmat_nl_y,
4719  Bmat_nl_u,
4720  Bmat_nl_v,
4721  Bmat_bend,
4722  Bmat_vk;
4723 
4724  std::vector<FEMOperatorMatrix>
4725  dBmat_lin(2);
4726 
4727  for (unsigned int i=0; i<2; i++)
4728  dBmat_lin[i].reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
4729 
4730  /*
4731  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
4732  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
4733  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
4734  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
4735  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
4736  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
4737  */
4738 
4739  // TODO: remove this const-cast, which may need change in API of
4740  // material card
4741  const MAST::FieldFunction<RealMatrixX >&
4742  mat_stiff =
4743  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
4744  stiffness_matrix(2);
4745 
4746  /*
4747  // get the thickness values for the bending strain calculation
4748  const MAST::FieldFunction<Real>
4749  &h = _property.get<MAST::FieldFunction<Real> >("h"),
4750  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
4751 
4752 
4753  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
4754  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
4755 
4756 
4757  // check to see if the element has any thermal loads specified
4758  // The object returns null
4759  MAST::BoundaryConditionBase *thermal_load =
4760  stress_output.get_thermal_load_for_elem(_elem);
4761 
4762  const MAST::FieldFunction<Real>
4763  *temp_func = nullptr,
4764  *ref_temp_func = nullptr,
4765  *alpha_func = nullptr;
4766 
4767  // get pointers to the temperature, if thermal load is specified
4768  if (thermal_load) {
4769  temp_func =
4770  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
4771  ref_temp_func =
4772  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
4773  alpha_func =
4774  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
4775  }
4776  */
4777 
4779  // get the material matrix
4780  mat_stiff(xyz[qp], _time, material_mat);
4781 
4782  this->initialize_strain_operator_gradient(qp,
4783  fe,
4784  _local_sol,
4785  dstrain,
4786  dBmat_lin);
4787 
4788  /*
4789  // if thermal load was specified, then set the thermal strain
4790  // component of the total strain
4791  if (thermal_load) {
4792  (*temp_func) (xyz[qp_loc_index], _time, temp);
4793  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
4794  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
4795  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
4796  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
4797  }
4798 
4799  if (if_bending) {
4800 
4801  // von Karman strain
4802  if (if_vk) { // get the vonKarman strain operator if needed
4803 
4804  this->initialize_von_karman_strain_operator(qp,
4805  fe,
4806  strain_vk,
4807  vk_dwdxi_mat,
4808  Bmat_vk);
4809  strain += strain_vk;
4810  }
4811 
4812  // add to this the bending strain
4813  // TODO: add coupling due to h_offset
4814  h (xyz[qp], _time, z);
4815  h_off(xyz[qp], _time, z_off);
4816  // TODO: this assumes isotropic section. Multilayered sections need
4817  // special considerations
4818  bend->initialize_bending_strain_operator_for_z(fe,
4819  qp,
4820  qp_loc[qp](2) * z/2.+z_off,
4821  Bmat_bend);
4822  Bmat_bend.vector_mult(strain_bend, _local_sol);
4823 
4824 
4825  // add stress due to bending.
4826  strain += strain_bend;
4827  }*/
4828 
4829  // note that this assumes linear material laws
4830  stress_grad = material_mat * dstrain;
4831 
4832  // calculate the derivative if requested
4833  if (dstress_grad_dX) {
4834 
4835  for (unsigned int i=0; i<2; i++) {
4836 
4837  dBmat_lin[i].left_multiply(dstrain_dX, eye);
4838 
4839  /*if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
4840 
4841  Bmat_nl_x.left_multiply(mat_n1n2, mat_x);
4842  dstrain_dX += mat_n1n2;
4843  Bmat_nl_y.left_multiply(mat_n1n2, mat_y);
4844  dstrain_dX += mat_n1n2;
4845  }
4846 
4847  if (if_bending) {
4848 
4849  // von Karman strain
4850  if (if_vk) {
4851 
4852  Bmat_vk.left_multiply(mat_n1n2, vk_dwdxi_mat);
4853  dstrain_dX += mat_n1n2;
4854  }
4855 
4856  // bending strain
4857  Bmat_bend.left_multiply(mat_n1n2, eye);
4858  dstrain_dX += mat_n1n2;
4859  }
4860  */
4861 
4862  // note: this assumes linear material laws
4863  (*dstress_grad_dX)[i] = material_mat * dstrain_dX;
4864  }
4865  }
4866 }
4867 
4868 
4869 
4870 
4871 void
4872 MAST::StructuralElement2D::
4873 _compute_stress_gradient_sensitivity(const MAST::FunctionBase& p,
4874  MAST::FEBase& fe,
4875  unsigned int qp,
4876  RealMatrixX& stress_grad) {
4877 
4878  /*
4879  MAST::BendingOperatorType bending_model =
4880  _property.bending_model(_elem);
4881 
4882  std::unique_ptr<MAST::BendingOperator2D>
4883  bend(MAST::build_bending_operator_2D(bending_model,
4884  *this,
4885  qp_loc_fe).release());
4886  */
4887 
4888  stress_grad.setZero();
4889 
4890  std::vector<Real> JxW = fe.get_JxW();
4891  const std::vector<libMesh::Point>& xyz = fe.get_xyz();
4892 
4893  const unsigned int
4894  n_phi = (unsigned int)fe.n_shape_functions(),
4895  n1 = this->n_direct_strain_components(),
4896  n2 = 6*n_phi,
4897  n3 = this->n_von_karman_strain_components();
4898 
4899  /*
4900  Real
4901  z = 0.,
4902  z_off = 0.,
4903  temp = 0.,
4904  ref_t = 0.,
4905  alpha = 0.,
4906  dtemp = 0.,
4907  dref_t= 0.,
4908  dalpha= 0.;
4909  */
4910 
4911  RealMatrixX
4912  material_mat,
4913  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
4914  dstrain = RealMatrixX::Zero(n1,2),
4915  dstrain_dX = RealMatrixX::Zero(n1,n2),
4916  mat_n1n2 = RealMatrixX::Zero(n1,n2),
4917  eye = RealMatrixX::Identity(n1, n1),
4918  mat_x = RealMatrixX::Zero(3,2),
4919  mat_y = RealMatrixX::Zero(3,2),
4920  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
4921  dstress_dX_3D= RealMatrixX::Zero(6,n2);
4922 
4923  RealVectorX
4924  strain_vk = RealVectorX::Zero(n1),
4925  strain_bend = RealVectorX::Zero(n1),
4926  strain_3D = RealVectorX::Zero(6),
4927  stress_3D = RealVectorX::Zero(6),
4928  dstrain_dp = RealVectorX::Zero(n1),
4929  dstress_dp = RealVectorX::Zero(n1),
4930  vec1 = RealVectorX::Zero(n2),
4931  vec2 = RealVectorX::Zero(n2);
4932 
4933 
4934  FEMOperatorMatrix
4935  Bmat_nl_x,
4936  Bmat_nl_y,
4937  Bmat_nl_u,
4938  Bmat_nl_v,
4939  Bmat_bend,
4940  Bmat_vk;
4941 
4942  std::vector<FEMOperatorMatrix>
4943  dBmat_lin(2);
4944 
4945  for (unsigned int i=0; i<2; i++)
4946  dBmat_lin[i].reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
4947 
4948  /*
4949  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
4950  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
4951  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
4952  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
4953  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
4954  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
4955  */
4956 
4957  // TODO: remove this const-cast, which may need change in API of
4958  // material card
4959  const MAST::FieldFunction<RealMatrixX >&
4960  mat_stiff =
4961  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
4962  stiffness_matrix(2);
4963 
4964  /*
4965  // get the thickness values for the bending strain calculation
4966  const MAST::FieldFunction<Real>
4967  &h = _property.get<MAST::FieldFunction<Real> >("h"),
4968  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
4969 
4970 
4971  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
4972  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
4973 
4974 
4975  // check to see if the element has any thermal loads specified
4976  // The object returns null
4977  MAST::BoundaryConditionBase *thermal_load =
4978  stress_output.get_thermal_load_for_elem(_elem);
4979 
4980  const MAST::FieldFunction<Real>
4981  *temp_func = nullptr,
4982  *ref_temp_func = nullptr,
4983  *alpha_func = nullptr;
4984 
4985  // get pointers to the temperature, if thermal load is specified
4986  if (thermal_load) {
4987  temp_func =
4988  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
4989  ref_temp_func =
4990  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
4991  alpha_func =
4992  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
4993  }
4994  */
4995 
4997  // get the material matrix
4998  mat_stiff.derivative(p, xyz[qp], _time, material_mat);
4999 
5000  this->initialize_strain_operator_gradient(qp,
5001  fe,
5002  _local_sol,
5003  dstrain,
5004  dBmat_lin);
5005 
5006  /*
5007  // if thermal load was specified, then set the thermal strain
5008  // component of the total strain
5009  if (thermal_load) {
5010  (*temp_func) (xyz[qp_loc_index], _time, temp);
5011  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
5012  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
5013  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
5014  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
5015  }
5016 
5017  if (if_bending) {
5018 
5019  // von Karman strain
5020  if (if_vk) { // get the vonKarman strain operator if needed
5021 
5022  this->initialize_von_karman_strain_operator(qp,
5023  fe,
5024  strain_vk,
5025  vk_dwdxi_mat,
5026  Bmat_vk);
5027  strain += strain_vk;
5028  }
5029 
5030  // add to this the bending strain
5031  // TODO: add coupling due to h_offset
5032  h (xyz[qp], _time, z);
5033  h_off(xyz[qp], _time, z_off);
5034  // TODO: this assumes isotropic section. Multilayered sections need
5035  // special considerations
5036  bend->initialize_bending_strain_operator_for_z(fe,
5037  qp,
5038  qp_loc[qp](2) * z/2.+z_off,
5039  Bmat_bend);
5040  Bmat_bend.vector_mult(strain_bend, _local_sol);
5041 
5042 
5043  // add stress due to bending.
5044  strain += strain_bend;
5045  }*/
5046 
5047  // note that this assumes linear material laws
5048  stress_grad = material_mat * dstrain;
5049 }
5050 
5051 
5052 void
5053 MAST::StructuralElement2D::_surface_normal_voigt_notation(const RealVectorX& normal,
5054  RealMatrixX& normal_mat) {
5055 
5056  libmesh_assert_equal_to(normal.size(), 3);
5057  libmesh_assert_equal_to(normal_mat.rows(), 2);
5058  libmesh_assert_equal_to(normal_mat.cols(), 3);
5059  normal_mat.setZero();
5060 
5061  /*
5062  // for a 3D stress-tensor the relation will be
5063  normal_mat = RealMatrixX::Zero(3, 6);
5064  // nx
5065  normal_mat(0,0) = normal_mat(1,3) = normal_mat(2,5) = normal(0);
5066  // ny
5067  normal_mat(0,3) = normal_mat(1,1) = normal_mat(2,4) = normal(1);
5068  // nz
5069  normal_mat(0,5) = normal_mat(1,4) = normal_mat(2,2) = normal(2);
5070  */
5071 
5072  // for a 2D stress tensor the relation will be
5073  normal_mat = RealMatrixX::Zero(2, 3);
5074 
5075  // nx
5076  normal_mat(0,0) = normal_mat(1,2) = normal(0);
5077 
5078  // ny
5079  normal_mat(0,2) = normal_mat(1,1) = normal(1);
5080 }
5081 
5082 
MAST::ElementPropertyCardBase::prestress_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > prestress_A_matrix(MAST::ElementBase &e) const =0
MAST::ElementPropertyCardBase::bending_model
virtual MAST::BendingOperatorType bending_model(const MAST::GeomElem &elem) const =0
returns the bending model to be used for the element.
MAST::StructuralElement2D::calculate_stress
virtual bool calculate_stress(bool request_derivative, const MAST::FunctionBase *p, MAST::StressStrainOutputBase &output)
Calculates the stress tensor.
Definition: structural_element_2d.cpp:307
MAST::StructuralElement2D::surface_traction_residual_shifted_boundary_sensitivity
virtual bool surface_traction_residual_shifted_boundary_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the sensitivity of force vector and Jacobian due to surface traction and sensitiity due to...
Definition: structural_element_2d.cpp:2663
MAST::ElementPropertyCardBase::thermal_expansion_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > thermal_expansion_B_matrix(const MAST::ElementBase &e) const =0
MAST::StructuralElement2D::_convert_prestress_B_mat_to_vector
void _convert_prestress_B_mat_to_vector(const RealMatrixX &mat, RealVectorX &vec) const
converts the prestress stress tensor to a vector representation
Definition: structural_element_2d.cpp:2039
MAST::StructuralElement2D::_surface_normal_voigt_notation
void _surface_normal_voigt_notation(const RealVectorX &normal, RealMatrixX &normal_mat)
creates a matrix that can be multiplied with the Voigt notation of stress to compute the surface norm...
Definition: structural_element_2d.cpp:5053
MAST::FunctionBase::is_shape_parameter
virtual bool is_shape_parameter() const
Definition: function_base.h:89
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::StructuralElement2D::surface_traction_residual_sensitivity
virtual bool surface_traction_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the sensitivity of element vector and matrix quantities for surface traction boundary cond...
Definition: structural_element_2d.cpp:2431
fe_base.h
MAST::StructuralElement2D::n_direct_strain_components
virtual unsigned int n_direct_strain_components()
row dimension of the direct strain matrix, also used for the bending operator row dimension ...
Definition: structural_element_2d.h:56
MAST::GeomElem::get_reference_elem
virtual const libMesh::Elem & get_reference_elem() const
Definition: geom_elem.cpp:54
MAST::StructuralElement2D::initialize_strain_operator_gradient
virtual void initialize_strain_operator_gradient(const unsigned int qp, const MAST::FEBase &fe, const RealVectorX &local_disp, RealMatrixX &epsilon_grad, std::vector< MAST::FEMOperatorMatrix > &dBmat_lin)
Definition: structural_element_2d.cpp:255
MAST::ElementPropertyCardBase::prestress_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > prestress_B_matrix(MAST::ElementBase &e) const =0
MAST::FEBase::get_qpoints
virtual const std::vector< libMesh::Point > & get_qpoints() const
Definition: fe_base.cpp:407
MAST::StructuralElement2D::_internal_residual_operation
virtual void _internal_residual_operation(bool if_vk, const unsigned int n2, const unsigned int qp, const MAST::FEBase &fe, const std::vector< Real > &JxW, bool request_jacobian, RealVectorX &local_f, RealMatrixX &local_jac, RealVectorX &local_disp, RealVectorX &strain_mem, MAST::BendingOperator2D *bend, FEMOperatorMatrix &Bmat_lin, FEMOperatorMatrix &Bmat_nl_x, FEMOperatorMatrix &Bmat_nl_y, FEMOperatorMatrix &Bmat_nl_u, FEMOperatorMatrix &Bmat_nl_v, MAST::FEMOperatorMatrix &Bmat_bend, MAST::FEMOperatorMatrix &Bmat_vk, RealMatrixX &mat_x, RealMatrixX &mat_y, RealMatrixX &stress, RealMatrixX &vk_dwdxi_mat, RealMatrixX &material_A_mat, RealMatrixX &material_B_mat, RealMatrixX &material_D_mat, RealVectorX &vec1_n1, RealVectorX &vec2_n1, RealVectorX &vec3_n2, RealVectorX &vec4_2, RealVectorX &vec5_2, RealVectorX &vec6_n2, RealMatrixX &mat1_n1n2, RealMatrixX &mat2_n2n2, RealMatrixX &mat3, RealMatrixX &mat4_2n2, RealMatrixX &mat5_3n2)
performs integration at the quadrature point for the provided matrices.
Definition: structural_element_2d.cpp:1767
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
structural_element_2d.h
MAST::ElementPropertyCardBase::get_material
virtual const MAST::MaterialPropertyCardBase & get_material() const
return the material property.
Definition: element_property_card_base.h:126
MAST::StructuralElement2D::prestress_residual_sensitivity
virtual bool prestress_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal residual vector and Jacobian due to strain energy coming from a prestress...
Definition: structural_element_2d.cpp:2205
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::Parameter
This is a scalar function whose value can be changed and one that can be used as a design variable in...
Definition: parameter.h:35
MAST::StructuralElement2D::internal_residual_jac_dot_state_sensitivity
virtual bool internal_residual_jac_dot_state_sensitivity(RealMatrixX &jac)
calculates d[J]/d{x} .
Definition: structural_element_2d.cpp:1572
fem_operator_matrix.h
MAST::FEBase::get_xyz
virtual const std::vector< libMesh::Point > & get_xyz() const
physical location of the quadrature point in the global coordinate system for the reference element ...
Definition: fe_base.cpp:228
MAST::FEBase::get_d2phi
virtual const std::vector< std::vector< libMesh::RealTensorValue > > & get_d2phi() const
Definition: fe_base.cpp:263
MAST::StructuralElement2D::_compute_stress
void _compute_stress(MAST::FEBase &fe, unsigned int qp, RealVectorX &stress, RealMatrixX *dstress_dX)
Computes the stress at the specified quadrature point of the FE data structure.
Definition: structural_element_2d.cpp:4178
MAST::FEBase::n_shape_functions
virtual unsigned int n_shape_functions() const
Definition: fe_base.cpp:239
bending_operator.h
MAST::ElementPropertyCardBase::stiffness_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > stiffness_B_matrix(const MAST::ElementBase &e) const =0
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::StructuralElementBase::_local_sol_sens
RealVectorX _local_sol_sens
local solution sensitivity
Definition: structural_element_base.h:764
MAST::StructuralElement2D::surface_traction_residual_shifted_boundary
virtual bool surface_traction_residual_shifted_boundary(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the sensitivity of force vector and Jacobian due to surface traction and sensitiity due to...
Definition: structural_element_2d.cpp:2510
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
MAST::StructuralElement2D::_compute_stress_sensitivity
void _compute_stress_sensitivity(const MAST::FunctionBase &f, MAST::FEBase &fe, unsigned int qp, RealVectorX &stress)
Definition: structural_element_2d.cpp:4477
MAST::FEBase::get_phi
virtual const std::vector< std::vector< Real > > & get_phi() const
Definition: fe_base.cpp:247
system_initialization.h
MAST::StructuralElement2D::thermal_residual_temperature_derivative
virtual void thermal_residual_temperature_derivative(const MAST::FEBase &fe_thermal, RealMatrixX &m)
Definition: structural_element_2d.cpp:3625
MAST::FEMOperatorMatrix::n
unsigned int n() const
Definition: fem_operator_matrix.h:52
stress_output_base.h
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
geom_elem.h
MAST::BendingOperator2D::initialize_bending_strain_operator
virtual void initialize_bending_strain_operator(const MAST::FEBase &fe, const unsigned int qp, MAST::FEMOperatorMatrix &Bmat)=0
initialze the bending strain operator for the specified quadrature point
MAST::StructuralElement2D::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: structural_element_2d.cpp:1076
MAST::StructuralElement2D::calculate_stress_temperature_derivative
virtual void calculate_stress_temperature_derivative(MAST::FEBase &fe_thermal, MAST::StressStrainOutputBase &output)
Definition: structural_element_2d.cpp:958
constant_field_function.h
MAST::StructuralElement2D::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: structural_element_2d.cpp:3818
MAST::FunctionSetBase::contains
bool contains(const std::string &nm) const
checks if the card contains the specified property value
Definition: function_set_base.cpp:35
MAST::StructuralElementBase::follower_forces
bool follower_forces
flag for follower forces
Definition: structural_element_base.h:439
MAST::StructuralElement2D::thermal_residual_boundary_velocity
virtual void thermal_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, MAST::BoundaryConditionBase &bc, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
calculates the term on side s: .
Definition: structural_element_2d.cpp:3421
MAST::build_bending_operator_2D
std::unique_ptr< MAST::BendingOperator2D > build_bending_operator_2D(MAST::BendingOperatorType type, MAST::StructuralElementBase &elem, const std::vector< libMesh::Point > &pts)
Definition: bending_operator.cpp:75
RealMatrixX
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Definition: mast_data_types.h:44
MAST::StructuralElement2D::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 the Jacobian due to thermal stresses.
Definition: structural_element_2d.cpp:3208
MAST::StressStrainOutputBase::get_thermal_load_for_elem
MAST::BoundaryConditionBase * get_thermal_load_for_elem(const MAST::GeomElem &elem)
Definition: stress_output_base.cpp:817
MAST::BendingOperator2D
Bending strain operator for 1D element.
Definition: bending_operator.h:157
boundary_condition_base.h
MAST::StructuralElement2D::prestress_residual
virtual bool prestress_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal residual vector and Jacobian due to strain energy coming from a prestress...
Definition: structural_element_2d.cpp:2054
MAST::StressStrainOutputBase::n_stress_strain_data_for_elem
unsigned int n_stress_strain_data_for_elem(const MAST::GeomElem &e) const
Definition: stress_output_base.cpp:748
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::StructuralElement2D::StructuralElement2D
StructuralElement2D(MAST::SystemInitialization &sys, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
Definition: structural_element_2d.cpp:39
MAST::ElementPropertyCardBase::stiffness_D_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > stiffness_D_matrix(const MAST::ElementBase &e) const =0
MAST::StructuralElement2D::initialize_von_karman_strain_operator
virtual void initialize_von_karman_strain_operator(const unsigned int qp, const MAST::FEBase &fe, RealVectorX &vk_strain, RealMatrixX &vk_dwdxi_mat, MAST::FEMOperatorMatrix &Bmat_vk)
initialze the von Karman strain in vK_strain, the operator matrices needed for Jacobian calculation...
Definition: structural_element_2d.cpp:57
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
MAST::StressStrainOutputBase::n_boundary_stress_strain_data_for_elem
unsigned int n_boundary_stress_strain_data_for_elem(const GeomElem &e) const
Definition: stress_output_base.cpp:765
assembly_base.h
MAST::StructuralElement2D::initialize_von_karman_strain_operator_sensitivity
virtual void initialize_von_karman_strain_operator_sensitivity(const unsigned int qp, const MAST::FEBase &fe, RealMatrixX &vk_dwdxi_mat_sens)
initialze the sensitivity of von Karman operator matrices needed for Jacobian calculation.
Definition: structural_element_2d.cpp:107
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::GeomElem::T_matrix
const RealMatrixX & T_matrix() const
O.
Definition: geom_elem.cpp:236
level_set_boundary_velocity.h
MAST::FEBase::get_normals_for_local_coordinate
virtual const std::vector< libMesh::Point > & get_normals_for_local_coordinate() const
normals defined in the coordinate system for the local reference element.
Definition: fe_base.cpp:399
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::GeomElem::dim
unsigned int dim() const
Definition: geom_elem.cpp:91
MAST::FEMOperatorMatrix::vector_mult
void vector_mult(T &res, const T &v) const
res = [this] * v
Definition: fem_operator_matrix.h:280
MAST::StructuralElement2D::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: structural_element_2d.cpp:2857
MAST::StructuralElement2D::internal_residual_sensitivity
virtual bool internal_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the sensitivity internal residual vector and Jacobian due to strain energy.
Definition: structural_element_2d.cpp:1234
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::StressStrainOutputBase::Data::set_derivatives
void set_derivatives(const RealMatrixX &dstress_dX, const RealMatrixX &dstrain_dX)
adds the derivative data
Definition: stress_output_base.cpp:89
material_property_card_base.h
element_property_card_2D.h
MAST::ElementPropertyCardBase::thermal_expansion_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > thermal_expansion_A_matrix(const MAST::ElementBase &e) const =0
MAST::FEBase::get_fe_type
libMesh::FEType get_fe_type() const
Definition: fe_base.cpp:212
MAST::StructuralElement2D::surface_traction_residual
virtual bool surface_traction_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to surface traction.
Definition: structural_element_2d.cpp:2358
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::StructuralElement2D::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 sensitivity of force vector and Jacobian due to piston-theory based surface pressure o...
Definition: structural_element_2d.cpp:3996
MAST::FEBase::get_JxW
virtual const std::vector< Real > & get_JxW() const
Definition: fe_base.cpp:220
MAST::StructuralElementBase::transform_matrix_to_global_system
void transform_matrix_to_global_system(const ValType &local_mat, ValType &global_mat) const
Definition: structural_element_base.cpp:1568
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::StructuralElement2D::_convert_prestress_A_mat_to_vector
void _convert_prestress_A_mat_to_vector(const RealMatrixX &mat, RealVectorX &vec) const
converts the prestress stress tensor to a vector representation
Definition: structural_element_2d.cpp:2025
MAST::StructuralElementBase::_local_sol
RealVectorX _local_sol
local solution
Definition: structural_element_base.h:752
MAST::StructuralElement2D::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: structural_element_2d.cpp:2932
MAST::StructuralElementBase::_local_vel
RealVectorX _local_vel
local velocity
Definition: structural_element_base.h:776
piston_theory_boundary_condition.h
MAST::StructuralElement2D::_compute_stress_gradient
void _compute_stress_gradient(MAST::FEBase &fe, unsigned int qp, RealMatrixX &stress_grad, std::vector< RealMatrixX > *dstress_grad_dX)
computes the gradient of stress in Voigt notation in stress_grad where the three columns are the deri...
Definition: structural_element_2d.cpp:4657
MAST::StructuralElement2D::calculate_stress_boundary_velocity
virtual void calculate_stress_boundary_velocity(const MAST::FunctionBase &p, MAST::StressStrainOutputBase &output, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f)
Calculates the boundary velocity term contributions to the sensitivity of stress at the specified bou...
Definition: structural_element_2d.cpp:716
MAST::StructuralElement2D::initialize_green_lagrange_strain_operator
virtual 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, MAST::FEMOperatorMatrix &Bmat_lin, MAST::FEMOperatorMatrix &Bmat_nl_x, MAST::FEMOperatorMatrix &Bmat_nl_y, MAST::FEMOperatorMatrix &Bmat_nl_u, MAST::FEMOperatorMatrix &Bmat_nl_v)
Definition: structural_element_2d.cpp:143
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::BendingOperatorType
BendingOperatorType
Definition: bending_operator.h:48
MAST::ElementBase::_time
const Real & _time
time for which system is being assembled
Definition: elem_base.h:219
MAST::StructuralElement2D::internal_residual_boundary_velocity
virtual void internal_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
calculates the term on side s: .
Definition: structural_element_2d.cpp:1403
MAST::StructuralElement2D::_compute_stress_gradient_sensitivity
void _compute_stress_gradient_sensitivity(const MAST::FunctionBase &f, MAST::FEBase &fe, unsigned int qp, RealMatrixX &stress_grad)
Definition: structural_element_2d.cpp:4873
MAST::StructuralElement2D::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: structural_element_2d.cpp:3012
MAST::ElementPropertyCardBase::extra_quadrature_order
virtual int extra_quadrature_order(const MAST::GeomElem &elem) const =0
returns the extra quadrature order (on top of the system) that this element should use...
MAST::StructuralElementBase::transform_vector_to_global_system
void transform_vector_to_global_system(const ValType &local_vec, ValType &global_vec) const
Definition: structural_element_base.cpp:1644
MAST::StructuralElement2D::n_von_karman_strain_components
virtual unsigned int n_von_karman_strain_components()
row dimension of the von Karman strain matrix
Definition: structural_element_2d.h:63
MAST::FunctionSetBase::depends_on
virtual bool depends_on(const MAST::FunctionBase &f) const
returns true if the property card depends on the function f
Definition: function_set_base.cpp:61
MAST::StressStrainOutputBase::Data::set_sensitivity
void set_sensitivity(const MAST::FunctionBase &f, const RealVectorX &dstress_df, const RealVectorX &dstrain_df)
sets the sensitivity of the data with respect to a function
Definition: stress_output_base.cpp:124
MAST::StressStrainOutputBase::add_stress_strain_at_boundary_qp_location
virtual MAST::StressStrainOutputBase::Data & add_stress_strain_at_boundary_qp_location(const MAST::GeomElem &e, const unsigned int s, const unsigned int qp, const libMesh::Point &quadrature_pt, const libMesh::Point &physical_pt, const RealVectorX &stress, const RealVectorX &strain, Real JxW_Vn)
add the stress tensor associated with the qp on side s of element e.
Definition: stress_output_base.cpp:660