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-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 // MAST includes
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 
36 
37 MAST::StructuralElement2D::
38 StructuralElement2D(MAST::SystemInitialization& sys,
39  MAST::AssemblyBase& assembly,
40  const MAST::GeomElem& elem,
41  const MAST::ElementPropertyCardBase& p):
42 MAST::BendingStructuralElem(sys, assembly, 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 bool
253 MAST::StructuralElement2D::calculate_stress(bool request_derivative,
254  const MAST::FunctionBase* p,
255  MAST::StressStrainOutputBase& output) {
256 
257  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
258 
259  const unsigned int
260  qp_loc_fe_size = (unsigned int)fe->get_qpoints().size(),
261  n_added_qp = 2;
262 
263  std::vector<libMesh::Point>
264  qp_loc_fe = fe->get_qpoints(),
265  qp_loc(qp_loc_fe.size()*n_added_qp);
266  std::vector<Real> JxW = fe->get_JxW();
267  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
268 
269  // we will evaluate the stress at upper and lower layers of element,
270  // so we will add two new points for each qp_loc.
271  // TODO: this will need to be moved to element property section card
272  // to handle composite elements at a later date
273  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
274 
275  qp_loc[i*2] = qp_loc_fe[i];
276  qp_loc[i*2](2) = +1.; // upper skin
277 
278  qp_loc[i*2+1] = qp_loc_fe[i];
279  qp_loc[i*2+1](2) = -1.; // lower skin
280 
282  // scale the JxW of each QP by (1/n_added_qp) since we are
283  // adding extra points and the total volume should add up to be
284  // the same.
286  JxW[i] /= (1.*n_added_qp);
287  }
288 
289 
290 
291  MAST::BendingOperatorType bending_model =
292  _property.bending_model(_elem);
293 
294  std::unique_ptr<MAST::BendingOperator2D>
295  bend(MAST::build_bending_operator_2D(bending_model,
296  *this,
297  qp_loc_fe).release());
298 
299  // now that the FE object has been initialized, evaluate the stress values
300 
301  const unsigned int
302  n_phi = (unsigned int)fe->n_shape_functions(),
303  n1 = this->n_direct_strain_components(),
304  n2 = 6*n_phi,
305  n3 = this->n_von_karman_strain_components();
306 
307  Real
308  z = 0.,
309  z_off = 0.,
310  temp = 0.,
311  ref_t = 0.,
312  alpha = 0.,
313  dtemp = 0.,
314  dref_t= 0.,
315  dalpha= 0.;
316 
317  RealMatrixX
318  material_mat,
319  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
320  dstrain_dX = RealMatrixX::Zero(n1,n2),
321  dstress_dX = RealMatrixX::Zero(n1,n2),
322  mat_n1n2 = RealMatrixX::Zero(n1,n2),
323  eye = RealMatrixX::Identity(n1, n1),
324  mat_x = RealMatrixX::Zero(3,2),
325  mat_y = RealMatrixX::Zero(3,2),
326  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
327  dstress_dX_3D= RealMatrixX::Zero(6,n2);
328 
329  RealVectorX
330  strain = RealVectorX::Zero(n1),
331  stress = RealVectorX::Zero(n1),
332  strain_vk = RealVectorX::Zero(n1),
333  strain_bend = RealVectorX::Zero(n1),
334  strain_3D = RealVectorX::Zero(6),
335  stress_3D = RealVectorX::Zero(6),
336  dstrain_dp = RealVectorX::Zero(n1),
337  dstress_dp = RealVectorX::Zero(n1);
338 
339 
340  FEMOperatorMatrix
341  Bmat_lin,
342  Bmat_nl_x,
343  Bmat_nl_y,
344  Bmat_nl_u,
345  Bmat_nl_v,
346  Bmat_bend,
347  Bmat_vk;
348 
349  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
350  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
351  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
352  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
353  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
354  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
355  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
356 
357  // TODO: remove this const-cast, which may need change in API of
358  // material card
359  const MAST::FieldFunction<RealMatrixX >&
360  mat_stiff =
361  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
362  stiffness_matrix(2);
363 
364  // get the thickness values for the bending strain calculation
365  const MAST::FieldFunction<Real>
366  &h = _property.get<MAST::FieldFunction<Real> >("h"),
367  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
368 
369 
370  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
371  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
372 
373  // a reference to the stress output data structure
374  MAST::StressStrainOutputBase& stress_output =
375  dynamic_cast<MAST::StressStrainOutputBase&>(output);
376 
377  // check to see if the element has any thermal loads specified
378  // The object returns null
379  MAST::BoundaryConditionBase *thermal_load =
380  stress_output.get_thermal_load_for_elem(_elem);
381 
382  const MAST::FieldFunction<Real>
383  *temp_func = nullptr,
384  *ref_temp_func = nullptr,
385  *alpha_func = nullptr;
386 
387  // get pointers to the temperature, if thermal load is specified
388  if (thermal_load) {
389  temp_func =
390  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
391  ref_temp_func =
392  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
393  alpha_func =
394  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
395  }
396 
397 
398 
400  unsigned int
401  qp = 0;
402  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
403  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
404  {
405  qp = qp_loc_index*n_added_qp + section_qp_index;
406 
407  // get the material matrix
408  mat_stiff(xyz[qp_loc_index], _time, material_mat);
409 
410  this->initialize_green_lagrange_strain_operator(qp_loc_index,
411  *fe,
412  _local_sol,
413  strain,
414  mat_x,
415  mat_y,
416  Bmat_lin,
417  Bmat_nl_x,
418  Bmat_nl_y,
419  Bmat_nl_u,
420  Bmat_nl_v);
421 
422  // if thermal load was specified, then set the thermal strain
423  // component of the total strain
424  if (thermal_load) {
425  (*temp_func) (xyz[qp_loc_index], _time, temp);
426  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
427  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
428  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
429  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
430  }
431 
432  if (if_bending) {
433 
434  // von Karman strain
435  if (if_vk) { // get the vonKarman strain operator if needed
436 
437  this->initialize_von_karman_strain_operator(qp_loc_index,
438  *fe,
439  strain_vk,
440  vk_dwdxi_mat,
441  Bmat_vk);
442  strain += strain_vk;
443  }
444 
445  // add to this the bending strain
446  // TODO: add coupling due to h_offset
447  h (xyz[qp_loc_index], _time, z);
448  h_off(xyz[qp_loc_index], _time, z_off);
449  // TODO: this assumes isotropic section. Multilayered sections need
450  // special considerations
451  bend->initialize_bending_strain_operator_for_z(*fe,
452  qp_loc_index,
453  qp_loc[qp](2) * z/2.+z_off,
454  Bmat_bend);
455  Bmat_bend.vector_mult(strain_bend, _local_sol);
456 
457 
458  // add stress due to bending.
459  strain += strain_bend;
460  }
461 
462  // note that this assumes linear material laws
463  stress = material_mat * strain;
464 
465 
466  // now set the data for the 3D stress-strain vector
467  // this is using only the direct strain/stress.
468  // this can be improved by estimating the shear stresses from
469  // torsion and shear flow from bending.
470  stress_3D(0) = stress(0); // sigma-xx
471  stress_3D(1) = stress(1); // sigma-yy
472  stress_3D(3) = stress(2); // tau-xy
473  strain_3D(0) = strain(0); // epsilon-xx
474  strain_3D(1) = strain(1); // epsilon-yy
475  strain_3D(3) = strain(2); // gamma-xy
476 
477  // set the stress and strain data
478  MAST::StressStrainOutputBase::Data*
479  data = nullptr;
480 
481  // if neither the derivative nor sensitivity is requested, then
482  // we assume that a new data entry is to be provided. Otherwise,
483  // we assume that the stress at this quantity already
484  // exists, and we only need to append sensitivity/derivative
485  // data to it
486  if (!request_derivative && !p)
487  data = &(stress_output.add_stress_strain_at_qp_location(_elem,
488  qp,
489  qp_loc[qp],
490  xyz[qp_loc_index],
491  stress_3D,
492  strain_3D,
493  JxW[qp_loc_index]));
494  else
495  data = &(stress_output.get_stress_strain_data_for_elem_at_qp(_elem, qp));
496 
497 
498  // calculate the derivative if requested
499  if (request_derivative || p) {
500 
501  Bmat_lin.left_multiply(dstrain_dX, eye);
502 
503  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
504 
505  Bmat_nl_x.left_multiply(mat_n1n2, mat_x);
506  dstrain_dX += mat_n1n2;
507  Bmat_nl_y.left_multiply(mat_n1n2, mat_y);
508  dstrain_dX += mat_n1n2;
509  }
510 
511  if (if_bending) {
512 
513  // von Karman strain
514  if (if_vk) {
515 
516  Bmat_vk.left_multiply(mat_n1n2, vk_dwdxi_mat);
517  dstrain_dX += mat_n1n2;
518  }
519 
520  // bending strain
521  Bmat_bend.left_multiply(mat_n1n2, eye);
522  dstrain_dX += mat_n1n2;
523  }
524 
525  // note: this assumes linear material laws
526  dstress_dX = material_mat * dstrain_dX;
527 
528  // copy to the 3D structure
529  dstress_dX_3D.row(0) = dstress_dX.row(0); // sigma-xx
530  dstress_dX_3D.row(1) = dstress_dX.row(1); // sigma-yy
531  dstress_dX_3D.row(3) = dstress_dX.row(2); // tau-xy
532  dstrain_dX_3D.row(0) = dstrain_dX.row(0); // epsilon-xx
533  dstrain_dX_3D.row(1) = dstrain_dX.row(1); // epsilon-yy
534  dstrain_dX_3D.row(3) = dstrain_dX.row(2); // gamma-xy
535 
536  if (request_derivative)
537  data->set_derivatives(dstress_dX_3D, dstrain_dX_3D);
538 
539 
540  if (p) {
541  // sensitivity of the response, s, is
542  // ds/dp = partial s/partial p +
543  // partial s/partial X dX/dp
544  // the first part of the sensitivity is obtained from
545  //
546  // the first term includes direct sensitivity of the stress
547  // with respect to the parameter, while holding the solution
548  // constant. This should include influence of shape changes,
549  // if the parameter is shape-dependent.
550  // TODO: include shape sensitivity.
551  // presently, only material parameter is included
552 
553 
554  dstrain_dp = RealVectorX::Zero(n1);
555 
556  // if thermal load was specified, then set the thermal strain
557  // component of the total strain
558  if (thermal_load) {
559  temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp);
560  ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t);
561  alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha);
562  dstrain_dp(0) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t); // epsilon-xx
563  dstrain_dp(1) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t); // epsilon-yy
564  }
565 
566 
567 
568  if (if_bending) {
569 
570  // add to this the bending strain
571  h.derivative (*p,
572  xyz[qp_loc_index], _time, z);
573  h_off.derivative(*p,
574  xyz[qp_loc_index], _time, z_off);
575  // TODO: this assumes isotropic section. Multilayered sections need
576  // special considerations
577  bend->initialize_bending_strain_operator_for_z(*fe,
578  qp_loc_index,
579  qp_loc[qp](2) * z/2.+z_off,
580  Bmat_bend);
581  Bmat_bend.vector_mult(strain_bend, _local_sol);
582 
583 
584  // add stress due to bending.
585  dstrain_dp += strain_bend;
586  }
587 
588 
589  // now use this to calculate the stress sensitivity.
590  dstress_dp = material_mat * dstrain_dp;
591 
592  // get the material matrix sensitivity
593  mat_stiff.derivative(*p,
594  xyz[qp_loc_index],
595  _time,
596  material_mat);
597 
598  // partial sensitivity of strain is zero unless it is a
599  // shape parameter.
600  // TODO: shape sensitivity of strain operator
601 
602  // now use this to calculate the stress sensitivity.
603  dstress_dp += material_mat * strain;
604 
605  //
606  // use the derivative data to evaluate the second term in the
607  // sensitivity
608  //
609  dstress_dp += dstress_dX * _local_sol_sens;
610  dstrain_dp += dstrain_dX * _local_sol_sens;
611 
612  // copy the 3D object
613  stress_3D(0) = dstress_dp(0); // sigma-xx
614  stress_3D(1) = dstress_dp(1); // sigma-yy
615  stress_3D(3) = dstress_dp(2); // tau-xy
616  strain_3D(0) = dstrain_dp(0); // epsilon-xx
617  strain_3D(1) = dstrain_dp(1); // epsilon-yy
618  strain_3D(3) = dstrain_dp(2); // gamma-xy
619 
620  // tell the data object about the sensitivity values
621  data->set_sensitivity(*p,
622  stress_3D,
623  strain_3D);
624  }
625  }
626  }
627 
628  // make sure that the number of data points for this element is
629  // the same as the number of requested points
630  libmesh_assert(qp_loc.size() ==
631  stress_output.n_stress_strain_data_for_elem(_elem));
632 
633  // if either derivative or sensitivity was requested, it was provided
634  // by this routine
635  return request_derivative || p;
636 }
637 
638 
639 
640 void
641 MAST::StructuralElement2D::
642 calculate_stress_boundary_velocity(const MAST::FunctionBase& p,
643  MAST::StressStrainOutputBase& output,
644  const unsigned int s,
645  const MAST::FieldFunction<RealVectorX>& vel_f) {
646 
647  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s, true));
648 
649  const unsigned int
650  qp_loc_fe_size = (unsigned int)fe->get_qpoints().size(),
651  n_added_qp = 2;
652 
653  std::vector<libMesh::Point>
654  qp_loc_fe = fe->get_qpoints(),
655  qp_loc(qp_loc_fe.size()*n_added_qp);
656  std::vector<Real> JxW_Vn = fe->get_JxW();
657  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
658  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
659 
660  // we will evaluate the stress at upper and lower layers of element,
661  // so we will add two new points for each qp_loc.
662  // TODO: this will need to be moved to element property section card
663  // to handle composite elements at a later date
664  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
665 
666  qp_loc[i*2] = qp_loc_fe[i];
667  qp_loc[i*2](2) = +1.; // upper skin
668 
669  qp_loc[i*2+1] = qp_loc_fe[i];
670  qp_loc[i*2+1](2) = -1.; // lower skin
671  }
672 
673  MAST::BendingOperatorType bending_model =
674  _property.bending_model(_elem);
675 
676  std::unique_ptr<MAST::BendingOperator2D>
677  bend(MAST::build_bending_operator_2D(bending_model,
678  *this,
679  qp_loc_fe).release());
680 
681  const unsigned int
682  n_phi = (unsigned int)fe->n_shape_functions(),
683  n1 = this->n_direct_strain_components(),
684  n2 = 6*n_phi,
685  n3 = this->n_von_karman_strain_components(),
686  dim = 2;
687 
688  Real
689  z = 0.,
690  z_off = 0.,
691  temp = 0.,
692  ref_t = 0.,
693  alpha = 0.,
694  vn = 0.;
695 
696  RealMatrixX
697  material_mat,
698  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
699  mat_n1n2 = RealMatrixX::Zero(n1,n2),
700  mat_x = RealMatrixX::Zero(3,2),
701  mat_y = RealMatrixX::Zero(3,2),
702  eye = RealMatrixX::Identity(n1, n1);
703 
704  RealVectorX
705  strain = RealVectorX::Zero(n1),
706  stress = RealVectorX::Zero(n1),
707  strain_vk = RealVectorX::Zero(n1),
708  strain_bend = RealVectorX::Zero(n1),
709  strain_3D = RealVectorX::Zero(6),
710  stress_3D = RealVectorX::Zero(6),
711  vel = RealVectorX::Zero(dim);
712 
713  // modify the JxW_Vn by multiplying the normal velocity to it
714  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
715 
716  vel_f(xyz[qp], _time, vel);
717  vn = 0.;
718  for (unsigned int i=0; i<dim; i++)
719  vn += vel(i)*face_normals[qp](i);
720  JxW_Vn[qp] *= vn/(1.*n_added_qp);
721  }
722 
723  FEMOperatorMatrix
724  Bmat_lin,
725  Bmat_nl_x,
726  Bmat_nl_y,
727  Bmat_nl_u,
728  Bmat_nl_v,
729  Bmat_bend,
730  Bmat_vk;
731 
732  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
733  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
734  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
735  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
736  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
737  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
738  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
739 
740  // TODO: remove this const-cast, which may need change in API of
741  // material card
742  const MAST::FieldFunction<RealMatrixX >&
743  mat_stiff =
744  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
745  stiffness_matrix(2);
746 
747  // get the thickness values for the bending strain calculation
748  const MAST::FieldFunction<Real>
749  &h = _property.get<MAST::FieldFunction<Real> >("h"),
750  &h_off = _property.get<MAST::FieldFunction<Real> >("off");
751 
752 
753  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
754  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
755 
756  // a reference to the stress output data structure
757  MAST::StressStrainOutputBase& stress_output =
758  dynamic_cast<MAST::StressStrainOutputBase&>(output);
759 
760  // check to see if the element has any thermal loads specified
761  // The object returns null
762  MAST::BoundaryConditionBase *thermal_load =
763  stress_output.get_thermal_load_for_elem(_elem);
764 
765  const MAST::FieldFunction<Real>
766  *temp_func = nullptr,
767  *ref_temp_func = nullptr,
768  *alpha_func = nullptr;
769 
770  // get pointers to the temperature, if thermal load is specified
771  if (thermal_load) {
772  temp_func =
773  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
774  ref_temp_func =
775  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
776  alpha_func =
777  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
778  }
779 
781  unsigned int
782  qp = 0;
783  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
784  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
785  {
786  qp = qp_loc_index*n_added_qp + section_qp_index;
787 
788  // get the material matrix
789  mat_stiff(xyz[qp_loc_index], _time, material_mat);
790 
791  this->initialize_green_lagrange_strain_operator(qp_loc_index,
792  *fe,
793  _local_sol,
794  strain,
795  mat_x,
796  mat_y,
797  Bmat_lin,
798  Bmat_nl_x,
799  Bmat_nl_y,
800  Bmat_nl_u,
801  Bmat_nl_v);
802 
803 
804  // if thermal load was specified, then set the thermal strain
805  // component of the total strain
806  if (thermal_load) {
807  (*temp_func) (xyz[qp_loc_index], _time, temp);
808  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
809  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
810  strain(0) -= alpha*(temp-ref_t); // epsilon-xx
811  strain(1) -= alpha*(temp-ref_t); // epsilon-yy
812  }
813 
814  if (if_bending) {
815 
816  // von Karman strain
817  if (if_vk) { // get the vonKarman strain operator if needed
818 
819  this->initialize_von_karman_strain_operator(qp_loc_index,
820  *fe,
821  strain_vk,
822  vk_dwdxi_mat,
823  Bmat_vk);
824  strain += strain_vk;
825  }
826 
827  // add to this the bending strain
828  // TODO: add coupling due to h_offset
829  h (xyz[qp_loc_index], _time, z);
830  h_off(xyz[qp_loc_index], _time, z_off);
831  // TODO: this assumes isotropic section. Multilayered sections need
832  // special considerations
833  bend->initialize_bending_strain_operator_for_z(*fe,
834  qp_loc_index,
835  qp_loc[qp](2) * z/2.+z_off,
836  Bmat_bend);
837  Bmat_bend.vector_mult(strain_bend, _local_sol);
838 
839 
840  // add stress due to bending.
841  strain += strain_bend;
842  }
843 
844  // note that this assumes linear material laws
845  stress = material_mat * strain;
846 
847 
848  // now set the data for the 3D stress-strain vector
849  // this is using only the direct strain/stress.
850  // this can be improved by estimating the shear stresses from
851  // torsion and shear flow from bending.
852  stress_3D(0) = stress(0); // sigma-xx
853  stress_3D(1) = stress(1); // sigma-yy
854  stress_3D(3) = stress(2); // tau-xy
855  strain_3D(0) = strain(0); // epsilon-xx
856  strain_3D(1) = strain(1); // epsilon-yy
857  strain_3D(3) = strain(2); // gamma-xy
858 
859  // if neither the derivative nor sensitivity is requested, then
860  // we assume that a new data entry is to be provided. Otherwise,
861  // we assume that the stress at this quantity already
862  // exists, and we only need to append sensitivity/derivative
863  // data to it
864  stress_output.add_stress_strain_at_boundary_qp_location(_elem,
865  s,
866  qp,
867  qp_loc[qp],
868  xyz[qp_loc_index],
869  stress_3D,
870  strain_3D,
871  JxW_Vn[qp_loc_index]);
872  }
873 
874  // make sure that the number of data points for this element is
875  // the same as the number of requested points
876  libmesh_assert(qp_loc.size() ==
877  stress_output.n_boundary_stress_strain_data_for_elem(_elem));
878 }
879 
880 
881 
882 void
883 MAST::StructuralElement2D::
884 calculate_stress_temperature_derivative(MAST::FEBase& fe_thermal,
885  MAST::StressStrainOutputBase& output) {
886 
887  std::unique_ptr<MAST::FEBase> fe_str(_elem.init_fe(true, false));
888 
889  // make sure that the structural and thermal FE have the same types and
890  // locations
891  libmesh_assert(fe_str->get_fe_type() == fe_thermal.get_fe_type());
892  // this is a weak assertion. We really want that the same qpoints be used,
893  // but we are assuming that if the elem type is the same, and if the
894  // number fo qpoints is the same, then the location of qpoints will also
895  // be the same.
896  libmesh_assert_equal_to(fe_str->get_qpoints().size(), fe_thermal.get_qpoints().size());
897 
898  MAST::StressStrainOutputBase&
899  stress_output = dynamic_cast<MAST::StressStrainOutputBase&>(output);
900 
901  libmesh_assert(stress_output.get_thermal_load_for_elem(_elem));
902 
903  const unsigned int
904  qp_loc_fe_size = (unsigned int)fe_str->get_qpoints().size(),
905  n_added_qp = 2;
906 
907  std::vector<libMesh::Point>
908  qp_loc_fe = fe_str->get_qpoints(),
909  qp_loc(qp_loc_fe.size()*n_added_qp);
910 
911  const std::vector<std::vector<Real>>& phi_temp = fe_thermal.get_phi();
912  const std::vector<libMesh::Point>& xyz = fe_str->get_xyz();
913 
914  // we will evaluate the stress at upper and lower layers of element,
915  // so we will add two new points for each qp_loc.
916  // TODO: this will need to be moved to element property section card
917  // to handle composite elements at a later date
918  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
919 
920  qp_loc[i*2] = qp_loc_fe[i];
921  qp_loc[i*2](2) = +1.; // upper skin
922 
923  qp_loc[i*2+1] = qp_loc_fe[i];
924  qp_loc[i*2+1](2) = -1.; // lower skin
925  }
926 
927 
928 
929  MAST::BendingOperatorType bending_model = _property.bending_model(_elem);
930 
931  std::unique_ptr<MAST::BendingOperator2D>
932  bend(MAST::build_bending_operator_2D(bending_model,
933  *this,
934  qp_loc_fe).release());
935 
936  // now that the FE object has been initialized, evaluate the stress values
937 
938  const unsigned int
939  nt = (unsigned int)fe_thermal.get_phi().size();
940 
941  Real
942  alpha = 0.;
943 
944  RealMatrixX
945  material_mat,
946  dstress_dX,
947  dT_mat = RealMatrixX::Zero(3,1),
948  dstrain_dX_3D = RealMatrixX::Zero(6,nt),
949  dstress_dX_3D = RealMatrixX::Zero(6,nt),
950  Bmat_temp = RealMatrixX::Zero(1,nt);
951 
952  dT_mat(0) = dT_mat(1) = 1.;
953 
954  // TODO: remove this const-cast, which may need change in API of
955  // material card
956  const MAST::FieldFunction<RealMatrixX >&
957  mat_stiff =
958  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).
959  stiffness_matrix(2);
960 
961  const MAST::FieldFunction<Real>
962  &alpha_func =
963  _property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion");
964 
965 
967  unsigned int
968  qp = 0;
969  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
970  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
971  {
972  qp = qp_loc_index*n_added_qp + section_qp_index;
973 
974  // shape function values for the temperature FE
975  for ( unsigned int i_nd=0; i_nd<nt; i_nd++ )
976  Bmat_temp(0, i_nd) = phi_temp[i_nd][qp_loc_index];
977 
978  // get the material matrix
979  mat_stiff(xyz[qp_loc_index], _time, material_mat);
980 
981  // if thermal load was specified, then set the thermal strain
982  // component of the total strain
983  alpha_func(xyz[qp_loc_index], _time, alpha);
984 
985  dstress_dX = alpha * material_mat * dT_mat * Bmat_temp;
986 
987  dstress_dX_3D.row(0) = -dstress_dX.row(0); // sigma-xx
988  dstress_dX_3D.row(1) = -dstress_dX.row(1); // sigma-yy
989  dstress_dX_3D.row(3) = -dstress_dX.row(2); // sigma-xy
990 
991  // set the stress and strain data
992  MAST::StressStrainOutputBase::Data
993  &data = stress_output.get_stress_strain_data_for_elem_at_qp(_elem, qp);
994  data.set_derivatives(dstress_dX_3D, dstrain_dX_3D);
995  }
996 }
997 
998 
999 
1000 
1001 bool
1002 MAST::StructuralElement2D::internal_residual (bool request_jacobian,
1003  RealVectorX& f,
1004  RealMatrixX& jac)
1005 {
1006  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1007  false,
1008  _property.extra_quadrature_order(_elem)));
1009 
1010  const std::vector<Real>& JxW = fe->get_JxW();
1011  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1012 
1013  const unsigned int
1014  n_phi = (unsigned int)fe->get_phi().size(),
1015  n1 = this->n_direct_strain_components(),
1016  n2 =6*n_phi,
1017  n3 = this->n_von_karman_strain_components();
1018 
1019  RealMatrixX
1020  material_A_mat,
1021  material_B_mat,
1022  material_D_mat,
1023  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1024  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1025  mat3,
1026  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1027  mat5_3n2 = RealMatrixX::Zero(3,n2),
1028  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1029  stress = RealMatrixX::Zero(2,2),
1030  mat_x = RealMatrixX::Zero(3,2),
1031  mat_y = RealMatrixX::Zero(3,2),
1032  local_jac = RealMatrixX::Zero(n2,n2);
1033 
1034  RealVectorX
1035  vec1_n1 = RealVectorX::Zero(n1),
1036  vec2_n1 = RealVectorX::Zero(n1),
1037  vec3_n2 = RealVectorX::Zero(n2),
1038  vec4_n3 = RealVectorX::Zero(n3),
1039  vec5_n3 = RealVectorX::Zero(n3),
1040  vec6_n2 = RealVectorX::Zero(n2),
1041  strain = RealVectorX::Zero(3),
1042  local_f = RealVectorX::Zero(n2);
1043 
1044  FEMOperatorMatrix
1045  Bmat_lin,
1046  Bmat_nl_x,
1047  Bmat_nl_y,
1048  Bmat_nl_u,
1049  Bmat_nl_v,
1050  Bmat_bend,
1051  Bmat_vk;
1052 
1053  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1054  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1055  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1056  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1057  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1058  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1059  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1060 
1061  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1062 
1063  MAST::BendingOperatorType
1064  bending_model = _property.bending_model(_elem);
1065 
1066  std::unique_ptr<MAST::BendingOperator2D> bend;
1067 
1068  if (bending_model != MAST::NO_BENDING)
1069  bend.reset(MAST::build_bending_operator_2D(bending_model,
1070  *this,
1071  fe->get_qpoints()).release());
1072 
1073  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1074  mat_stiff_A = _property.stiffness_A_matrix(*this),
1075  mat_stiff_B = _property.stiffness_B_matrix(*this),
1076  mat_stiff_D = _property.stiffness_D_matrix(*this);
1077 
1078  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1079 
1080  // get the material matrix
1081  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1082 
1083  if (bend.get()) {
1084  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1085  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1086  }
1087 
1088  // now calculte the quantity for these matrices
1089  _internal_residual_operation(if_vk,
1090  n2,
1091  qp,
1092  *fe,
1093  JxW,
1094  request_jacobian,
1095  local_f,
1096  local_jac,
1097  _local_sol,
1098  strain,
1099  bend.get(),
1100  Bmat_lin,
1101  Bmat_nl_x,
1102  Bmat_nl_y,
1103  Bmat_nl_u,
1104  Bmat_nl_v,
1105  Bmat_bend,
1106  Bmat_vk,
1107  mat_x,
1108  mat_y,
1109  stress,
1110  vk_dwdxi_mat,
1111  material_A_mat,
1112  material_B_mat,
1113  material_D_mat,
1114  vec1_n1,
1115  vec2_n1,
1116  vec3_n2,
1117  vec4_n3,
1118  vec5_n3,
1119  vec6_n2,
1120  mat1_n1n2,
1121  mat2_n2n2,
1122  mat3,
1123  mat4_n3n2,
1124  mat5_3n2);
1125  }
1126 
1127 
1128  // now calculate the transverse shear contribution if appropriate for the
1129  // element
1130  if (bend.get() &&
1131  bend->include_transverse_shear_energy())
1132  bend->calculate_transverse_shear_residual(request_jacobian,
1133  local_f,
1134  local_jac);
1135 
1136 
1137  // now transform to the global coorodinate system
1138  transform_vector_to_global_system(local_f, vec3_n2);
1139  f += vec3_n2;
1140 
1141  if (request_jacobian) {
1142  // for 2D elements
1143  if (_elem.dim() == 2) {
1144  // add small values to the diagonal of the theta_z dofs
1145  for (unsigned int i=0; i<n_phi; i++)
1146  local_jac(5*n_phi+i, 5*n_phi+i) = 1.0e-8;
1147  }
1148  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1149  jac += mat2_n2n2;
1150  }
1151 
1152  return request_jacobian;
1153 }
1154 
1155 
1156 
1157 
1158 
1159 bool
1160 MAST::StructuralElement2D::internal_residual_sensitivity (const MAST::FunctionBase& p,
1161  bool request_jacobian,
1162  RealVectorX& f,
1163  RealMatrixX& jac)
1164 {
1165  // this should be true if the function is called
1166  libmesh_assert(!p.is_shape_parameter()); // this is not implemented for now
1167 
1168 
1169  // check if the material property or the provided exterior
1170  // values, like temperature, are functions of the sensitivity parameter
1171  bool calculate = false;
1172  calculate = calculate || _property.depends_on(p);
1173 
1174  // nothing to be calculated if the element does not depend on the
1175  // sensitivity parameter.
1176  if (!calculate)
1177  return false;
1178 
1179  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1180  false,
1181  _property.extra_quadrature_order(_elem)));
1182 
1183  const std::vector<Real>& JxW = fe->get_JxW();
1184  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1185 
1186  const unsigned int
1187  n_phi = (unsigned int)fe->get_phi().size(),
1188  n1 = this->n_direct_strain_components(),
1189  n2 =6*n_phi,
1190  n3 = this->n_von_karman_strain_components();
1191 
1192  RealMatrixX
1193  material_A_mat,
1194  material_B_mat,
1195  material_D_mat,
1196  material_trans_shear_mat,
1197  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1198  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1199  mat3,
1200  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1201  mat5_3n2 = RealMatrixX::Zero(3,n2),
1202  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1203  stress = RealMatrixX::Zero(2,2),
1204  mat_x = RealMatrixX::Zero(3,2),
1205  mat_y = RealMatrixX::Zero(3,2),
1206  local_jac = RealMatrixX::Zero(n2,n2);
1207  RealVectorX
1208  vec1_n1 = RealVectorX::Zero(n1),
1209  vec2_n1 = RealVectorX::Zero(n1),
1210  vec3_n2 = RealVectorX::Zero(n2),
1211  vec4_n3 = RealVectorX::Zero(n3),
1212  vec5_n3 = RealVectorX::Zero(n3),
1213  vec6_n2 = RealVectorX::Zero(n2),
1214  strain = RealVectorX::Zero(3),
1215  local_f = RealVectorX::Zero(n2);
1216 
1217  FEMOperatorMatrix
1218  Bmat_lin,
1219  Bmat_nl_x,
1220  Bmat_nl_y,
1221  Bmat_nl_u,
1222  Bmat_nl_v,
1223  Bmat_bend,
1224  Bmat_vk;
1225 
1226  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1227  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1228  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1229  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1230  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1231  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1232  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1233 
1234  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1235 
1236  MAST::BendingOperatorType
1237  bending_model = _property.bending_model(_elem);
1238 
1239  std::unique_ptr<MAST::BendingOperator2D> bend;
1240 
1241  if (bending_model != MAST::NO_BENDING)
1242  bend.reset(MAST::build_bending_operator_2D(bending_model,
1243  *this,
1244  fe->get_qpoints()).release());
1245 
1246  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1247  mat_stiff_A = _property.stiffness_A_matrix(*this),
1248  mat_stiff_B = _property.stiffness_B_matrix(*this),
1249  mat_stiff_D = _property.stiffness_D_matrix(*this);
1250 
1251  // first calculate the sensitivity due to the parameter
1252  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1253 
1254  // get the material matrix
1255  mat_stiff_A->derivative(p, xyz[qp], _time, material_A_mat);
1256 
1257  if (bend.get()) {
1258 
1259  mat_stiff_B->derivative(p, xyz[qp], _time, material_B_mat);
1260  mat_stiff_D->derivative(p, xyz[qp], _time, material_D_mat);
1261  }
1262 
1263  // now calculte the quantity for these matrices
1264  // this accounts for the sensitivity of the material property matrices
1265  _internal_residual_operation(if_vk,
1266  n2,
1267  qp,
1268  *fe,
1269  JxW,
1270  request_jacobian,
1271  local_f,
1272  local_jac,
1273  _local_sol,
1274  strain,
1275  bend.get(),
1276  Bmat_lin,
1277  Bmat_nl_x,
1278  Bmat_nl_y,
1279  Bmat_nl_u,
1280  Bmat_nl_v,
1281  Bmat_bend,
1282  Bmat_vk,
1283  mat_x,
1284  mat_y,
1285  stress,
1286  vk_dwdxi_mat,
1287  material_A_mat,
1288  material_B_mat,
1289  material_D_mat,
1290  vec1_n1,
1291  vec2_n1,
1292  vec3_n2,
1293  vec4_n3,
1294  vec5_n3,
1295  vec6_n2,
1296  mat1_n1n2,
1297  mat2_n2n2,
1298  mat3,
1299  mat4_n3n2,
1300  mat5_3n2);
1301  }
1302 
1303  // now calculate the transverse shear contribution if appropriate for the
1304  // element
1305  if (bend.get() &&
1306  bend->include_transverse_shear_energy())
1307  bend->calculate_transverse_shear_residual_sensitivity(p,
1308  request_jacobian,
1309  local_f,
1310  local_jac);
1311 
1312  // now transform to the global coorodinate system
1313  transform_vector_to_global_system(local_f, vec3_n2);
1314  f += vec3_n2;
1315 
1316  if (request_jacobian) {
1317  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1318  jac += mat2_n2n2;
1319  }
1320 
1321 
1322  return request_jacobian;
1323 }
1324 
1325 
1326 
1327 void
1328 MAST::StructuralElement2D::
1329 internal_residual_boundary_velocity(const MAST::FunctionBase& p,
1330  const unsigned int s,
1331  const MAST::FieldFunction<RealVectorX>& vel_f,
1332  bool request_jacobian,
1333  RealVectorX& f,
1334  RealMatrixX& jac) {
1335 
1336  // prepare the side finite element
1337  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s,
1338  true,
1339  _property.extra_quadrature_order(_elem)));
1340 
1341  std::vector<Real> JxW_Vn = fe->get_JxW();
1342  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1343  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
1344 
1345  const unsigned int
1346  n_phi = (unsigned int)fe->get_phi().size(),
1347  n1 = this->n_direct_strain_components(),
1348  n2 = 6*n_phi,
1349  n3 = this->n_von_karman_strain_components(),
1350  dim = 2;
1351 
1352  RealMatrixX
1353  material_A_mat,
1354  material_B_mat,
1355  material_D_mat,
1356  material_trans_shear_mat,
1357  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1358  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1359  mat3,
1360  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1361  mat5_3n2 = RealMatrixX::Zero(3,n2),
1362  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1363  stress = RealMatrixX::Zero(2,2),
1364  mat_x = RealMatrixX::Zero(3,2),
1365  mat_y = RealMatrixX::Zero(3,2),
1366  local_jac = RealMatrixX::Zero(n2,n2);
1367  RealVectorX
1368  vec1_n1 = RealVectorX::Zero(n1),
1369  vec2_n1 = RealVectorX::Zero(n1),
1370  vec3_n2 = RealVectorX::Zero(n2),
1371  vec4_n3 = RealVectorX::Zero(n3),
1372  vec5_n3 = RealVectorX::Zero(n3),
1373  vec6_n2 = RealVectorX::Zero(n2),
1374  strain = RealVectorX::Zero(3),
1375  local_f = RealVectorX::Zero(n2),
1376  vel = RealVectorX::Zero(dim);
1377 
1378  FEMOperatorMatrix
1379  Bmat_lin,
1380  Bmat_nl_x,
1381  Bmat_nl_y,
1382  Bmat_nl_u,
1383  Bmat_nl_v,
1384  Bmat_bend,
1385  Bmat_vk;
1386 
1387  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1388  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1389  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1390  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1391  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1392  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1393  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1394 
1395  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1396 
1397  MAST::BendingOperatorType
1398  bending_model = _property.bending_model(_elem);
1399 
1400  std::unique_ptr<MAST::BendingOperator2D> bend;
1401 
1402  if (bending_model != MAST::NO_BENDING)
1403  bend.reset(MAST::build_bending_operator_2D(bending_model,
1404  *this,
1405  fe->get_qpoints()).release());
1406 
1407  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1408  mat_stiff_A = _property.stiffness_A_matrix(*this),
1409  mat_stiff_B = _property.stiffness_B_matrix(*this),
1410  mat_stiff_D = _property.stiffness_D_matrix(*this);
1411 
1412  Real
1413  vn = 0.;
1414 
1415  // modify the JxW_Vn by multiplying the normal velocity to it
1416  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
1417 
1418  vel_f(xyz[qp], _time, vel);
1419  vn = 0.;
1420  for (unsigned int i=0; i<dim; i++)
1421  vn += vel(i)*face_normals[qp](i);
1422  JxW_Vn[qp] *= vn;
1423  }
1424 
1425 
1426  // first calculate the sensitivity due to the parameter
1427  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
1428 
1429  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1430 
1431  if (bend.get()) {
1432 
1433  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1434  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1435  }
1436 
1437  // now calculte the quantity for these matrices
1438  // this accounts for the sensitivity of the material property matrices
1439  _internal_residual_operation(if_vk,
1440  n2,
1441  qp,
1442  *fe,
1443  JxW_Vn,
1444  request_jacobian,
1445  local_f,
1446  local_jac,
1447  _local_sol,
1448  strain,
1449  bend.get(),
1450  Bmat_lin,
1451  Bmat_nl_x,
1452  Bmat_nl_y,
1453  Bmat_nl_u,
1454  Bmat_nl_v,
1455  Bmat_bend,
1456  Bmat_vk,
1457  mat_x,
1458  mat_y,
1459  stress,
1460  vk_dwdxi_mat,
1461  material_A_mat,
1462  material_B_mat,
1463  material_D_mat,
1464  vec1_n1,
1465  vec2_n1,
1466  vec3_n2,
1467  vec4_n3,
1468  vec5_n3,
1469  vec6_n2,
1470  mat1_n1n2,
1471  mat2_n2n2,
1472  mat3,
1473  mat4_n3n2,
1474  mat5_3n2);
1475  }
1476 
1477  // now calculate the transverse shear contribution if appropriate for the
1478  // element
1479  if (bend.get() && bend->include_transverse_shear_energy())
1480  bend->calculate_transverse_shear_residual_boundary_velocity
1481  (p, s, vel_f, request_jacobian, local_f, local_jac);
1482 
1483  // now transform to the global coorodinate system
1484  transform_vector_to_global_system(local_f, vec3_n2);
1485  f += vec3_n2;
1486 
1487  if (request_jacobian) {
1488  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1489  jac += mat2_n2n2;
1490  }
1491 }
1492 
1493 
1494 
1495 
1496 bool
1497 MAST::StructuralElement2D::
1498 internal_residual_jac_dot_state_sensitivity (RealMatrixX& jac) {
1499 
1500  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1501  false,
1502  _property.extra_quadrature_order(_elem)));
1503 
1504  const std::vector<Real>& JxW = fe->get_JxW();
1505  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1506  const unsigned int
1507  n_phi = (unsigned int)fe->get_phi().size(),
1508  n1 = this->n_direct_strain_components(),
1509  n2 = 6*n_phi,
1510  n3 = this->n_von_karman_strain_components();
1511 
1512  RealMatrixX
1513  material_A_mat,
1514  material_B_mat,
1515  material_D_mat,
1516  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1517  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1518  mat3,
1519  vk_dwdxi_mat_sens = RealMatrixX::Zero(n1,n3),
1520  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1521  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1522  stress = RealMatrixX::Zero(2,2),
1523  mat_x = RealMatrixX::Zero(3,2),
1524  mat_y = RealMatrixX::Zero(3,2),
1525  local_jac = RealMatrixX::Zero(n2,n2);
1526  RealVectorX
1527  vec1_n1 = RealVectorX::Zero(n1),
1528  vec2_n1 = RealVectorX::Zero(n1),
1529  strain = RealVectorX::Zero(3);
1530 
1531 
1532  FEMOperatorMatrix
1533  Bmat_lin,
1534  Bmat_nl_x,
1535  Bmat_nl_y,
1536  Bmat_nl_u,
1537  Bmat_nl_v,
1538  Bmat_bend,
1539  Bmat_vk;
1540 
1541  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1542  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
1543  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
1544  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
1545  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
1546  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
1547  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1548 
1549  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1550 
1551  MAST::BendingOperatorType
1552  bending_model = _property.bending_model(_elem);
1553 
1554  std::unique_ptr<MAST::BendingOperator2D> bend;
1555 
1556  if (bending_model != MAST::NO_BENDING)
1557  bend.reset(MAST::build_bending_operator_2D(bending_model,
1558  *this,
1559  fe->get_qpoints()).release());
1560 
1561  // without the nonlinear strain, this matrix is zero.
1562  if (!if_vk)
1563  return false;
1564 
1565  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1566  mat_stiff_A = _property.stiffness_A_matrix(*this),
1567  mat_stiff_B = _property.stiffness_B_matrix(*this),
1568  mat_stiff_D = _property.stiffness_D_matrix(*this);
1569 
1570 
1571  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1572 
1573  // get the material matrix
1574  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
1575  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
1576  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
1577 
1578  this->initialize_green_lagrange_strain_operator(qp,
1579  *fe,
1580  _local_sol,
1581  strain,
1582  mat_x,
1583  mat_y,
1584  Bmat_lin,
1585  Bmat_nl_x,
1586  Bmat_nl_y,
1587  Bmat_nl_u,
1588  Bmat_nl_v);
1589 
1590  // first handle constant throught the thickness stresses: membrane and vonKarman
1591  Bmat_lin.vector_mult(vec1_n1, _local_sol_sens);
1592  vec2_n1 = material_A_mat * vec1_n1; // linear direct stress
1593 
1594  // get the bending strain operator
1595  if (bend.get()) {
1596  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
1597 
1598  // evaluate the bending stress and add that to the stress vector
1599  // for evaluation in the nonlinear stress term
1600  Bmat_bend.vector_mult(vec1_n1, _local_sol_sens);
1601  vec2_n1 += material_B_mat * vec1_n1;
1602 
1603  if (if_vk) { // get the vonKarman strain operator if needed
1604 
1605  this->initialize_von_karman_strain_operator(qp,
1606  *fe,
1607  vec1_n1, // epsilon_vk
1608  vk_dwdxi_mat,
1609  Bmat_vk);
1610  this->initialize_von_karman_strain_operator_sensitivity(qp,
1611  *fe,
1612  vk_dwdxi_mat_sens);
1613 
1614  // sensitivity of von Karman strain
1615  vec1_n1(0) = (vk_dwdxi_mat(0,0)*vk_dwdxi_mat_sens(0,0)); // dw/dx dwp/dx
1616  vec1_n1(1) = (vk_dwdxi_mat(1,1)*vk_dwdxi_mat_sens(1,1)); // dw/dy dwp/dy
1617  vec1_n1(2) = (vk_dwdxi_mat(0,0)*vk_dwdxi_mat_sens(1,1) + // dw/dx dwp/dy
1618  vk_dwdxi_mat(1,1)*vk_dwdxi_mat_sens(0,0)); // dw/dy dwp/dx
1619 
1620  vec2_n1 += material_A_mat * vec1_n1;
1621  }
1622  }
1623 
1624  // copy the stress values to a matrix
1625  stress(0,0) = vec2_n1(0); // sigma_xx
1626  stress(1,1) = vec2_n1(1); // sigma_yy
1627  stress(0,1) = vec2_n1(2); // gamma_xy
1628  stress(1,0) = vec2_n1(2); // gamma_xy
1629 
1630 
1631  // copy the stress to use here.
1632  vec2_n1.setZero();
1633 
1634  // now calculate the matrix
1635  // vk - membrane: w-displacement with sens
1636  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1637  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1638  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1639  local_jac += JxW[qp] * mat2_n2n2;
1640 
1641  // vk - bending: w-displacement with stress sens
1642  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
1643  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1644  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1645  local_jac += JxW[qp] * mat2_n2n2;
1646 
1647  // vk - vk: with stress sens and stress
1648  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1649  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1650  mat3 = vk_dwdxi_mat_sens.transpose() * material_A_mat * mat3;
1651  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1652  local_jac += JxW[qp] * mat2_n2n2;
1653 
1654  // vk - vk: with stress and stress sens
1655  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1656  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1657  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1658  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1659  local_jac += JxW[qp] * mat2_n2n2;
1660 
1661  // membrane - vk: w-displacement with sens
1662  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1663  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1664  mat3 = material_A_mat * mat3;
1665  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat3);
1666  local_jac += JxW[qp] * mat2_n2n2;
1667 
1668  // bending - vk: w-displacement with stress sens
1669  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1670  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1671  mat3 = material_B_mat.transpose() * mat3;
1672  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat3);
1673  local_jac += JxW[qp] * mat2_n2n2;
1674 
1675  // vk - vk: w-displacement with stress sens
1676  mat3 = RealMatrixX::Zero(2, n2);
1677  Bmat_vk.left_multiply(mat3, stress);
1678  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1679  local_jac += JxW[qp] * mat2_n2n2;
1680  }
1681 
1682  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1683  jac += mat2_n2n2;
1684 
1685  return true;
1686 }
1687 
1688 
1689 
1690 
1691 void
1692 MAST::StructuralElement2D::_internal_residual_operation
1693 (bool if_vk,
1694  const unsigned int n2,
1695  const unsigned int qp,
1696  const MAST::FEBase& fe,
1697  const std::vector<Real>& JxW,
1698  bool request_jacobian,
1699  RealVectorX& local_f,
1700  RealMatrixX& local_jac,
1701  RealVectorX& local_disp,
1702  RealVectorX& strain_mem,
1703  MAST::BendingOperator2D* bend,
1704  FEMOperatorMatrix& Bmat_lin,
1705  FEMOperatorMatrix& Bmat_nl_x,
1706  FEMOperatorMatrix& Bmat_nl_y,
1707  FEMOperatorMatrix& Bmat_nl_u,
1708  FEMOperatorMatrix& Bmat_nl_v,
1709  MAST::FEMOperatorMatrix& Bmat_bend,
1710  MAST::FEMOperatorMatrix& Bmat_vk,
1711  RealMatrixX& mat_x,
1712  RealMatrixX& mat_y,
1713  RealMatrixX& stress,
1714  RealMatrixX& vk_dwdxi_mat,
1715  RealMatrixX& material_A_mat,
1716  RealMatrixX& material_B_mat,
1717  RealMatrixX& material_D_mat,
1718  RealVectorX& vec1_n1,
1719  RealVectorX& vec2_n1,
1720  RealVectorX& vec3_n2,
1721  RealVectorX& vec4_2,
1722  RealVectorX& vec5_2,
1723  RealVectorX& vec6_n2,
1724  RealMatrixX& mat1_n1n2,
1725  RealMatrixX& mat2_n2n2,
1726  RealMatrixX& mat3,
1727  RealMatrixX& mat4_2n2,
1728  RealMatrixX& mat5_3n2) {
1729 
1730  this->initialize_green_lagrange_strain_operator(qp,
1731  fe,
1732  local_disp,
1733  strain_mem,
1734  mat_x,
1735  mat_y,
1736  Bmat_lin,
1737  Bmat_nl_x,
1738  Bmat_nl_y,
1739  Bmat_nl_u,
1740  Bmat_nl_v);
1741 
1742  vec2_n1 = material_A_mat * strain_mem; // membrane stress
1743 
1744  if (bend) {
1745 
1746  // get the bending strain operator
1747  bend->initialize_bending_strain_operator(fe, qp, Bmat_bend);
1748  Bmat_bend.vector_mult(vec1_n1, local_disp);
1749  vec2_n1 += material_B_mat * vec1_n1;
1750 
1751  if (if_vk) { // get the vonKarman strain operator if needed
1752  this->initialize_von_karman_strain_operator(qp,
1753  fe,
1754  vec1_n1, // epsilon_vk
1755  vk_dwdxi_mat,
1756  Bmat_vk);
1757 
1758  strain_mem += vec1_n1; // epsilon_mem + epsilon_vk
1759  vec2_n1 += material_A_mat * vec1_n1; // stress
1760  }
1761  }
1762 
1763 
1764  // copy the stress values to a matrix
1765  stress(0,0) = vec2_n1(0); // sigma_xx
1766  stress(0,1) = vec2_n1(2); // sigma_xy
1767  stress(1,0) = vec2_n1(2); // sigma_yx
1768  stress(1,1) = vec2_n1(1); // sigma_yy
1769 
1770  // now the internal force vector
1771  // this includes the membrane strain operator with all A and B material operators
1772  Bmat_lin.vector_mult_transpose(vec3_n2, vec2_n1);
1773  local_f.topRows(n2) += JxW[qp] * vec3_n2;
1774 
1775  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1776 
1777  // nonlinear strain operator
1778  // x
1779  vec4_2 = mat_x.transpose() * vec2_n1;
1780  Bmat_nl_x.vector_mult_transpose(vec6_n2, vec4_2);
1781  local_f.topRows(n2) += JxW[qp] * vec6_n2;
1782 
1783  // y
1784  vec4_2 = mat_y.transpose() * vec2_n1;
1785  Bmat_nl_y.vector_mult_transpose(vec6_n2, vec4_2);
1786  local_f.topRows(n2) += JxW[qp] * vec6_n2;
1787  }
1788 
1789 
1790  if (bend) {
1791  if (if_vk) {
1792  // von Karman strain
1793  vec4_2 = vk_dwdxi_mat.transpose() * vec2_n1;
1794  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
1795  local_f += JxW[qp] * vec3_n2;
1796  }
1797 
1798  // now coupling with the bending strain
1799  // B_bend^T [B] B_mem
1800  vec1_n1 = material_B_mat.transpose() * strain_mem;
1801  Bmat_bend.vector_mult_transpose(vec3_n2, vec1_n1);
1802  local_f += JxW[qp] * vec3_n2;
1803 
1804  // now bending stress
1805  Bmat_bend.vector_mult(vec2_n1, local_disp);
1806  vec1_n1 = material_D_mat * vec2_n1;
1807  Bmat_bend.vector_mult_transpose(vec3_n2, vec1_n1);
1808  local_f += JxW[qp] * vec3_n2;
1809  }
1810 
1811  if (request_jacobian) {
1812  // membrane - membrane
1813  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1814  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1815  local_jac += JxW[qp] * mat2_n2n2;
1816 
1817  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
1818 
1819 
1820  // B_lin^T C mat_x B_x
1821  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1822  mat1_n1n2 = material_A_mat * mat1_n1n2;
1823  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1824  local_jac += JxW[qp] * mat2_n2n2;
1825 
1826  // B_lin^T C mat_y B_y
1827  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1828  mat1_n1n2 = material_A_mat * mat1_n1n2;
1829  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1830  local_jac += JxW[qp] * mat2_n2n2;
1831 
1832  // B_x^T mat_x^T C B_lin
1833  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1834  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1835  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1836  local_jac += JxW[qp] * mat2_n2n2;
1837 
1838  // B_x^T mat_x^T C mat_x B_x
1839  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1840  mat1_n1n2 = material_A_mat * mat1_n1n2;
1841  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1842  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1843  local_jac += JxW[qp] * mat2_n2n2;
1844 
1845  // B_x^T mat_x^T C mat_y B_y
1846  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1847  mat1_n1n2 = material_A_mat * mat1_n1n2;
1848  mat5_3n2 = mat_x.transpose() * mat1_n1n2;
1849  Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1850  local_jac += JxW[qp] * mat2_n2n2;
1851 
1852  // B_y^T mat_y^T C B_lin
1853  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1854  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1855  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1856  local_jac += JxW[qp] * mat2_n2n2;
1857 
1858  // B_y^T mat_y^T C mat_x B_x
1859  Bmat_nl_x.left_multiply(mat1_n1n2, mat_x);
1860  mat1_n1n2 = material_A_mat * mat1_n1n2;
1861  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1862  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1863  local_jac += JxW[qp] * mat2_n2n2;
1864 
1865  // B_y^T mat_y^T C mat_y B_y
1866  Bmat_nl_y.left_multiply(mat1_n1n2, mat_y);
1867  mat1_n1n2 = material_A_mat * mat1_n1n2;
1868  mat5_3n2 = mat_y.transpose() * mat1_n1n2;
1869  Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1870  local_jac += JxW[qp] * mat2_n2n2;
1871 
1872  // nonlinear membrane u - nonlinear membrane u
1873  // (geometric / initial stress stiffness)
1874  // u-disp
1875  Bmat_nl_u.left_multiply(mat5_3n2, stress);
1876  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1877  local_jac += JxW[qp] * mat2_n2n2;
1878 
1879  // nonlinear membrane v - nonlinear membrane v
1880  // (geometric / initial stress stiffness)
1881  // v-disp
1882  Bmat_nl_v.left_multiply(mat5_3n2, stress);
1883  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat5_3n2);
1884  local_jac += JxW[qp] * mat2_n2n2;
1885  }
1886 
1887  if (bend) {
1888  if (if_vk) {
1889  // membrane - vk
1890  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1891  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1892  mat3 = material_A_mat * mat3;
1893  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat3);
1894  local_jac += JxW[qp] * mat2_n2n2;
1895 
1896  // vk - membrane
1897  Bmat_lin.left_multiply(mat1_n1n2, material_A_mat);
1898  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1899  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1900  local_jac += JxW[qp] * mat2_n2n2;
1901 
1902  // vk - vk
1903  mat3 = RealMatrixX::Zero(2, n2);
1904  Bmat_vk.left_multiply(mat3, stress);
1905  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1906  local_jac += JxW[qp] * mat2_n2n2;
1907 
1908  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1909  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1910  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1911  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1912  local_jac += JxW[qp] * mat2_n2n2;
1913 
1914  // bending - vk
1915  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1916  Bmat_vk.left_multiply(mat3, vk_dwdxi_mat);
1917  mat3 = material_B_mat.transpose() * mat3;
1918  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat3);
1919  local_jac += JxW[qp] * mat2_n2n2;
1920 
1921  // vk - bending
1922  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
1923  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1924  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
1925  local_jac += JxW[qp] * mat2_n2n2;
1926  }
1927 
1928  // bending - membrane
1929  mat3 = material_B_mat.transpose();
1930  Bmat_lin.left_multiply(mat1_n1n2, mat3);
1931  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1932  local_jac += JxW[qp] * mat2_n2n2;
1933 
1934  // membrane - bending
1935  Bmat_bend.left_multiply(mat1_n1n2, material_B_mat);
1936  Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1937  local_jac += JxW[qp] * mat2_n2n2;
1938 
1939  // bending - bending
1940  Bmat_bend.left_multiply(mat1_n1n2, material_D_mat);
1941  Bmat_bend.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1942  local_jac += JxW[qp] * mat2_n2n2;
1943  }
1944  }
1945 }
1946 
1947 
1948 
1949 
1950 void
1951 MAST::StructuralElement2D::_convert_prestress_A_mat_to_vector(const RealMatrixX& mat,
1952  RealVectorX& vec) const {
1953 
1954  libmesh_assert_equal_to(mat.rows(), 2);
1955  libmesh_assert_equal_to(mat.cols(), 2);
1956  vec = RealVectorX::Zero(3);
1957  vec(0) = mat(0,0); // sigma x
1958  vec(1) = mat(1,1); // sigma y
1959  vec(2) = mat(0,1); // tau xy
1960 }
1961 
1962 
1963 
1964 void
1965 MAST::StructuralElement2D::_convert_prestress_B_mat_to_vector(const RealMatrixX& mat,
1966  RealVectorX& vec) const {
1967 
1968  libmesh_assert_equal_to(mat.rows(), 2);
1969  libmesh_assert_equal_to(mat.cols(), 2);
1970  vec = RealVectorX::Zero(3);
1971  vec(0) = mat(0,0); // sigma x
1972  vec(1) = mat(1,1); // sigma y
1973  vec(2) = mat(0,1); // tau xy
1974 }
1975 
1976 
1977 
1978 
1979 bool
1980 MAST::StructuralElement2D::prestress_residual (bool request_jacobian,
1981  RealVectorX& f,
1982  RealMatrixX& jac)
1983 {
1984  if (!_property.if_prestressed())
1985  return false;
1986 
1987  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1988  false,
1989  _property.extra_quadrature_order(_elem)));
1990 
1991  const std::vector<Real>& JxW = fe->get_JxW();
1992  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1993  const unsigned int
1994  n_phi = (unsigned int)fe->get_phi().size(),
1995  n1 = this->n_direct_strain_components(),
1996  n2 = 6*n_phi,
1997  n3 = this->n_von_karman_strain_components();
1998 
1999  RealMatrixX
2000  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2001  mat3,
2002  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2003  local_jac = RealMatrixX::Zero(n2,n2),
2004  mat_x = RealMatrixX::Zero(3,2),
2005  mat_y = RealMatrixX::Zero(3,2),
2006  prestress_mat_A,
2007  prestress_mat_B;
2008 
2009  RealVectorX
2010  vec2_n1 = RealVectorX::Zero(n1),
2011  vec3_n2 = RealVectorX::Zero(n2),
2012  vec4_n3 = RealVectorX::Zero(n3),
2013  vec5_n3 = RealVectorX::Zero(n3),
2014  local_f = RealVectorX::Zero(n2),
2015  strain = RealVectorX::Zero(3),
2016  prestress_vec_A,
2017  prestress_vec_B;
2018 
2019  FEMOperatorMatrix
2020  Bmat_lin,
2021  Bmat_nl_x,
2022  Bmat_nl_y,
2023  Bmat_nl_u,
2024  Bmat_nl_v,
2025  Bmat_bend,
2026  Bmat_vk;
2027 
2028  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2029  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2030  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2031  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2032  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2033  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2034  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2035 
2036  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2037 
2038  MAST::BendingOperatorType
2039  bending_model = _property.bending_model(_elem);
2040 
2041  std::unique_ptr<MAST::BendingOperator2D> bend;
2042 
2043  if (bending_model != MAST::NO_BENDING)
2044  bend.reset(MAST::build_bending_operator_2D(bending_model,
2045  *this,
2046  fe->get_qpoints()).release());
2047 
2048  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
2049  prestress_A = _property.prestress_A_matrix(*this),
2050  prestress_B = _property.prestress_B_matrix(*this);
2051 
2052  // now calculate the quantity
2053  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2054 
2055  (*prestress_A)(xyz[qp], _time, prestress_mat_A);
2056  (*prestress_B)(xyz[qp], _time, prestress_mat_B);
2057  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
2058  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
2059 
2060  this->initialize_green_lagrange_strain_operator(qp,
2061  *fe,
2062  _local_sol,
2063  strain,
2064  mat_x,
2065  mat_y,
2066  Bmat_lin,
2067  Bmat_nl_x,
2068  Bmat_nl_y,
2069  Bmat_nl_u,
2070  Bmat_nl_v);
2071 
2072  // get the bending strain operator if needed
2073  vec2_n1.setZero(); // used to store vk strain, if applicable
2074  if (bend.get()) {
2075  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2076 
2077  if (if_vk) // get the vonKarman strain operator if needed
2078  this->initialize_von_karman_strain_operator(qp,
2079  *fe,
2080  vec2_n1,
2081  vk_dwdxi_mat,
2082  Bmat_vk);
2083  }
2084 
2085  // first handle constant throught the thickness stresses: membrane and vonKarman
2086  // multiply this with the constant through the thickness strain
2087  // membrane strain
2088  Bmat_lin.vector_mult_transpose(vec3_n2, prestress_vec_A);
2089  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
2090 
2091  if (bend.get()) {
2092  if (if_vk) {
2093  // von Karman strain
2094  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
2095  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_n3);
2096  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
2097  }
2098 
2099  // now coupling with the bending strain
2100  Bmat_bend.vector_mult_transpose(vec3_n2, prestress_vec_B);
2101  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
2102  }
2103 
2104  if (request_jacobian) {
2105  if (bend.get() && if_vk) {
2106  mat3 = RealMatrixX::Zero(2, n2);
2107  Bmat_vk.left_multiply(mat3, prestress_mat_A);
2108  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2109  local_jac += JxW[qp] * mat2_n2n2;
2110  }
2111  }
2112  }
2113 
2114  // now transform to the global coorodinate system
2115  transform_vector_to_global_system(local_f, vec3_n2);
2116  f += vec3_n2;
2117  if (request_jacobian && if_vk) {
2118  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2119  jac += mat2_n2n2;
2120  }
2121 
2122  // only the nonlinear strain returns a Jacobian for prestressing
2123  return (request_jacobian);
2124 }
2125 
2126 
2127 
2128 
2129 
2130 bool
2131 MAST::StructuralElement2D::prestress_residual_sensitivity (const MAST::FunctionBase& p,
2132  bool request_jacobian,
2133  RealVectorX& f,
2134  RealMatrixX& jac)
2135 {
2136  if (!_property.if_prestressed())
2137  return false;
2138 
2139  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2140  false,
2141  _property.extra_quadrature_order(_elem)));
2142 
2143  const std::vector<Real>& JxW = fe->get_JxW();
2144  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2145  const unsigned int
2146  n_phi = (unsigned int)fe->get_phi().size(),
2147  n1 = this->n_direct_strain_components(),
2148  n2 = 6*n_phi,
2149  n3 = this->n_von_karman_strain_components();
2150 
2151  RealMatrixX
2152  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2153  mat3,
2154  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2155  local_jac = RealMatrixX::Zero(n2,n2),
2156  mat_x = RealMatrixX::Zero(3,2),
2157  mat_y = RealMatrixX::Zero(3,2),
2158  prestress_mat_A,
2159  prestress_mat_B;
2160 
2161  RealVectorX
2162  vec2_n1 = RealVectorX::Zero(n1),
2163  vec3_n2 = RealVectorX::Zero(n2),
2164  vec4_n3 = RealVectorX::Zero(n3),
2165  vec5_n3 = RealVectorX::Zero(n3),
2166  local_f = RealVectorX::Zero(n2),
2167  strain = RealVectorX::Zero(3),
2168  prestress_vec_A,
2169  prestress_vec_B;
2170 
2171  FEMOperatorMatrix
2172  Bmat_lin,
2173  Bmat_nl_x,
2174  Bmat_nl_y,
2175  Bmat_nl_u,
2176  Bmat_nl_v,
2177  Bmat_bend,
2178  Bmat_vk;
2179 
2180  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2181  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2182  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2183  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2184  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2185  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2186  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2187 
2188  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2189 
2190  MAST::BendingOperatorType
2191  bending_model = _property.bending_model(_elem);
2192 
2193  std::unique_ptr<MAST::BendingOperator2D> bend;
2194 
2195  if (bending_model != MAST::NO_BENDING)
2196  bend.reset(MAST::build_bending_operator_2D(bending_model,
2197  *this,
2198  fe->get_qpoints()).release());
2199 
2200  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
2201  prestress_A = _property.prestress_A_matrix(*this),
2202  prestress_B = _property.prestress_B_matrix(*this);
2203 
2204  // transform to the local coordinate system
2205  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2206 
2207  prestress_A->derivative(p, xyz[qp], _time, prestress_mat_A);
2208  prestress_B->derivative(p, xyz[qp], _time, prestress_mat_B);
2209  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
2210  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
2211 
2212  this->initialize_green_lagrange_strain_operator(qp,
2213  *fe,
2214  _local_sol,
2215  strain,
2216  mat_x,
2217  mat_y,
2218  Bmat_lin,
2219  Bmat_nl_x,
2220  Bmat_nl_y,
2221  Bmat_nl_u,
2222  Bmat_nl_v);
2223 
2224  // get the bending strain operator if needed
2225  vec2_n1.setZero(); // used to store vk strain, if applicable
2226  if (bend.get()) {
2227  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2228 
2229  if (if_vk) // get the vonKarman strain operator if needed
2230  this->initialize_von_karman_strain_operator(qp,
2231  *fe,
2232  vec2_n1,
2233  vk_dwdxi_mat,
2234  Bmat_vk);
2235  }
2236 
2237  // first handle constant throught the thickness stresses: membrane and vonKarman
2238  // multiply this with the constant through the thickness strain
2239  // membrane strain
2240  Bmat_lin.vector_mult_transpose(vec3_n2, prestress_vec_A);
2241  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
2242 
2243  if (bend.get()) {
2244  if (if_vk) {
2245  // von Karman strain
2246  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
2247  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_n3);
2248  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
2249  }
2250 
2251  // now coupling with the bending strain
2252  Bmat_bend.vector_mult_transpose(vec3_n2, prestress_vec_B);
2253  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
2254  }
2255 
2256  if (request_jacobian) {
2257  if (bend.get() && if_vk) {
2258  mat3 = RealMatrixX::Zero(2, n2);
2259  Bmat_vk.left_multiply(mat3, prestress_mat_A);
2260  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2261  local_jac += JxW[qp] * mat2_n2n2;
2262  }
2263  }
2264  }
2265 
2266  // now transform to the global coorodinate system
2267  transform_vector_to_global_system(local_f, vec3_n2);
2268  f += vec3_n2;
2269  if (request_jacobian && if_vk) {
2270  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2271  jac += mat2_n2n2;
2272  }
2273 
2274  // only the nonlinear strain returns a Jacobian for prestressing
2275  return (request_jacobian);
2276 }
2277 
2278 
2279 
2280 
2281 
2282 bool
2283 MAST::StructuralElement2D::
2284 surface_pressure_residual(bool request_jacobian,
2285  RealVectorX &f,
2286  RealMatrixX &jac,
2287  const unsigned int side,
2288  MAST::BoundaryConditionBase& bc) {
2289 
2290  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2291 
2292  // prepare the side finite element
2293  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false));
2294 
2295  const std::vector<Real> &JxW = fe->get_JxW();
2296  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2297  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2298  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2299  const unsigned int
2300  n_phi = (unsigned int)phi.size(),
2301  n1 = 3,
2302  n2 = 6*n_phi;
2303 
2304 
2305  // get the function from this boundary condition
2306  const MAST::FieldFunction<Real>& p_func =
2307  bc.get<MAST::FieldFunction<Real> >("pressure");
2308 
2309  // get the thickness function to calculate the force
2310  const MAST::FieldFunction<Real>& t_func =
2311  _property.get<MAST::FieldFunction<Real> >("h");
2312 
2313  FEMOperatorMatrix Bmat;
2314  Real
2315  press = 0.,
2316  t_val = 0.;
2317 
2318  RealVectorX
2319  phi_vec = RealVectorX::Zero(n_phi),
2320  force = RealVectorX::Zero(2*n1),
2321  local_f = RealVectorX::Zero(n2),
2322  vec_n2 = RealVectorX::Zero(n2);
2323 
2324  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2325 
2326  // now set the shape function values
2327  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2328  phi_vec(i_nd) = phi[i_nd][qp];
2329 
2330  Bmat.reinit(2*n1, phi_vec);
2331 
2332  // get pressure and thickness values
2333  p_func(qpoint[qp], _time, press);
2334  t_func(qpoint[qp], _time, t_val);
2335 
2336  // calculate force
2337  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2338  force(i_dim) = (press*t_val) * face_normals[qp](i_dim);
2339 
2340  Bmat.vector_mult_transpose(vec_n2, force);
2341 
2342  local_f += JxW[qp] * vec_n2;
2343  }
2344 
2345  // now transform to the global system and add
2346  transform_vector_to_global_system(local_f, vec_n2);
2347  f -= vec_n2;
2348 
2349 
2350  return (request_jacobian);
2351 }
2352 
2353 
2354 
2355 
2356 
2357 bool
2358 MAST::StructuralElement2D::
2359 surface_pressure_residual_sensitivity(const MAST::FunctionBase& p,
2360  bool request_jacobian,
2361  RealVectorX &f,
2362  RealMatrixX &jac,
2363  const unsigned int side,
2364  MAST::BoundaryConditionBase& bc) {
2365 
2366  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2367 
2368  // prepare the side finite element
2369  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false));
2370 
2371  const std::vector<Real> &JxW = fe->get_JxW();
2372  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2373  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2374  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2375  const unsigned int
2376  n_phi = (unsigned int)phi.size(),
2377  n1 = 3,
2378  n2 = 6*n_phi;
2379 
2380 
2381  // get the function from this boundary condition
2382  const MAST::FieldFunction<Real>& p_func =
2383  bc.get<MAST::FieldFunction<Real> >("pressure");
2384 
2385  // get the thickness function to calculate the force
2386  const MAST::FieldFunction<Real>& t_func =
2387  _property.get<MAST::FieldFunction<Real> >("h");
2388 
2389 
2390  FEMOperatorMatrix Bmat;
2391  Real
2392  press = 0.,
2393  dpress = 0.,
2394  t_val = 0.,
2395  dt_val = 0.;
2396 
2397  RealVectorX
2398  phi_vec = RealVectorX::Zero(n_phi),
2399  force = RealVectorX::Zero(2*n1),
2400  local_f = RealVectorX::Zero(n2),
2401  vec_n2 = RealVectorX::Zero(n2);
2402 
2403  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2404 
2405  // now set the shape function values
2406  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2407  phi_vec(i_nd) = phi[i_nd][qp];
2408 
2409  Bmat.reinit(2*n1, phi_vec);
2410 
2411  // get pressure and thickness values and their sensitivities
2412  p_func(qpoint[qp], _time, press);
2413  p_func.derivative(p, qpoint[qp], _time, dpress);
2414  t_func(qpoint[qp], _time, t_val);
2415  t_func.derivative(p, qpoint[qp], _time, dt_val);
2416 
2417  // calculate force
2418  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2419  force(i_dim) = (press*dt_val + dpress*t_val) * face_normals[qp](i_dim);
2420 
2421  Bmat.vector_mult_transpose(vec_n2, force);
2422 
2423  local_f += JxW[qp] * vec_n2;
2424  }
2425 
2426  // now transform to the global system and add
2427  transform_vector_to_global_system(local_f, vec_n2);
2428  f -= vec_n2;
2429 
2430 
2431  return (request_jacobian);
2432 }
2433 
2434 
2435 
2436 
2437 
2438 bool
2439 MAST::StructuralElement2D::thermal_residual (bool request_jacobian,
2440  RealVectorX& f,
2441  RealMatrixX& jac,
2442  MAST::BoundaryConditionBase& bc)
2443 {
2444 
2445  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2446  false,
2447  _property.extra_quadrature_order(_elem)));
2448 
2449  const std::vector<Real>& JxW = fe->get_JxW();
2450  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2451 
2452  const unsigned int
2453  n_phi = (unsigned int)fe->get_phi().size(),
2454  n1 = this->n_direct_strain_components(),
2455  n2 = 6*n_phi,
2456  n3 = this->n_von_karman_strain_components();
2457 
2458  RealMatrixX
2459  material_exp_A_mat,
2460  material_exp_B_mat,
2461  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
2462  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2463  mat3 = RealMatrixX::Zero(2, n2),
2464  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
2465  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2466  stress = RealMatrixX::Zero(2,2),
2467  mat_x = RealMatrixX::Zero(3,2),
2468  mat_y = RealMatrixX::Zero(3,2),
2469  local_jac = RealMatrixX::Zero(n2,n2);
2470 
2471  RealVectorX
2472  vec1_n1 = RealVectorX::Zero(n1),
2473  vec2_n1 = RealVectorX::Zero(n1),
2474  vec3_n2 = RealVectorX::Zero(n2),
2475  vec4_2 = RealVectorX::Zero(2),
2476  vec5_n3 = RealVectorX::Zero(n3),
2477  local_f = RealVectorX::Zero(n2),
2478  strain = RealVectorX::Zero(3),
2479  delta_t = RealVectorX::Zero(1);
2480 
2481  FEMOperatorMatrix
2482  Bmat_lin,
2483  Bmat_nl_x,
2484  Bmat_nl_y,
2485  Bmat_nl_u,
2486  Bmat_nl_v,
2487  Bmat_bend,
2488  Bmat_vk;
2489 
2490  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2491  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2492  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2493  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2494  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2495  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2496  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2497 
2498  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2499 
2500  MAST::BendingOperatorType
2501  bending_model = _property.bending_model(_elem);
2502 
2503  std::unique_ptr<MAST::BendingOperator2D> bend;
2504 
2505  if (bending_model != MAST::NO_BENDING)
2506  bend.reset(MAST::build_bending_operator_2D(bending_model,
2507  *this,
2508  fe->get_qpoints()).release());
2509 
2510  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
2511  expansion_A = _property.thermal_expansion_A_matrix(*this),
2512  expansion_B = _property.thermal_expansion_B_matrix(*this);
2513 
2514  const MAST::FieldFunction<Real>
2515  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
2516  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
2517 
2518  Real
2519  t = 0.,
2520  t0 = 0.,
2521  scaling = 1.;
2522 
2523  if (bc.contains("thermal_jacobian_scaling"))
2524  bc.get<MAST::FieldFunction<Real>>("thermal_jacobian_scaling")(scaling);
2525 
2526  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2527 
2528  // this is moved inside the domain since
2529  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
2530  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
2531 
2532  // get the temperature function
2533  temp_func(xyz[qp], _time, t);
2534  ref_temp_func(xyz[qp], _time, t0);
2535  delta_t(0) = t-t0;
2536 
2537  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)} (with membrane strain)
2538  vec2_n1 = material_exp_B_mat * delta_t; // [C]{alpha (T - T0)} (with bending strain)
2539  stress(0,0) = vec1_n1(0); // sigma_xx
2540  stress(0,1) = vec1_n1(2); // sigma_xy
2541  stress(1,0) = vec1_n1(2); // sigma_yx
2542  stress(1,1) = vec1_n1(1); // sigma_yy
2543 
2544  this->initialize_green_lagrange_strain_operator(qp,
2545  *fe,
2546  _local_sol,
2547  strain,
2548  mat_x,
2549  mat_y,
2550  Bmat_lin,
2551  Bmat_nl_x,
2552  Bmat_nl_y,
2553  Bmat_nl_u,
2554  Bmat_nl_v);
2555 
2556  // membrane strain
2557  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
2558  local_f += JxW[qp] * vec3_n2;
2559 
2560  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
2561 
2562  // nonlinear strain operotor
2563  // x
2564  vec4_2 = mat_x.transpose() * vec1_n1;
2565  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
2566  local_f.topRows(n2) += JxW[qp] * vec3_n2;
2567 
2568  // y
2569  vec4_2 = mat_y.transpose() * vec1_n1;
2570  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
2571  local_f.topRows(n2) += JxW[qp] * vec3_n2;
2572  }
2573 
2574  if (bend.get()) {
2575  // bending strain
2576  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2577  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
2578  local_f += JxW[qp] * vec3_n2;
2579 
2580  // von Karman strain
2581  if (if_vk) {
2582  // get the vonKarman strain operator if needed
2583  this->initialize_von_karman_strain_operator(qp,
2584  *fe,
2585  vec2_n1, // epsilon_vk
2586  vk_dwdxi_mat,
2587  Bmat_vk);
2588  // von Karman strain
2589  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
2590  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
2591  local_f += JxW[qp] * vec3_n2;
2592  }
2593  }
2594 
2595  if (request_jacobian &&
2596  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
2597 
2598  // u-disp
2599  Bmat_nl_u.left_multiply(mat3, stress);
2600  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
2601  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
2602 
2603  // v-disp
2604  Bmat_nl_v.left_multiply(mat3, stress);
2605  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
2606  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
2607  }
2608 
2609  if (request_jacobian && if_vk) {
2610  // vk - vk
2611  mat3 = RealMatrixX::Zero(2, n2);
2612  Bmat_vk.left_multiply(mat3, stress);
2613  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2614  local_jac += JxW[qp] * mat2_n2n2;
2615  }
2616  }
2617 
2618  // now transform to the global coorodinate system
2619  transform_vector_to_global_system(local_f, vec3_n2);
2620  f -= vec3_n2;
2621  if (request_jacobian && if_vk) {
2622  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2623  jac -= scaling * mat2_n2n2;
2624  }
2625 
2626  // Jacobian contribution from von Karman strain
2627  return request_jacobian;
2628 }
2629 
2630 
2631 
2632 
2633 bool
2634 MAST::StructuralElement2D::
2635 thermal_residual_sensitivity (const MAST::FunctionBase& p,
2636  bool request_jacobian,
2637  RealVectorX& f,
2638  RealMatrixX& jac,
2639  MAST::BoundaryConditionBase& bc)
2640 {
2641  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2642  false,
2643  _property.extra_quadrature_order(_elem)));
2644 
2645  const std::vector<Real>& JxW = fe->get_JxW();
2646  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2647 
2648  const unsigned int
2649  n_phi = (unsigned int)fe->get_phi().size(),
2650  n1 = this->n_direct_strain_components(),
2651  n2 = 6*n_phi,
2652  n3 = this->n_von_karman_strain_components();
2653 
2654  RealMatrixX
2655  material_exp_A_mat,
2656  material_exp_B_mat,
2657  material_exp_A_mat_sens,
2658  material_exp_B_mat_sens,
2659  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
2660  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2661  mat3 = RealMatrixX::Zero(2, n2),
2662  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
2663  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2664  stress = RealMatrixX::Zero(2,2),
2665  mat_x = RealMatrixX::Zero(3,2),
2666  mat_y = RealMatrixX::Zero(3,2),
2667  local_jac = RealMatrixX::Zero(n2,n2);
2668 
2669  RealVectorX
2670  vec1_n1 = RealVectorX::Zero(n1),
2671  vec2_n1 = RealVectorX::Zero(n1),
2672  vec3_n2 = RealVectorX::Zero(n2),
2673  vec4_2 = RealVectorX::Zero(2),
2674  vec5_n1 = RealVectorX::Zero(n1),
2675  local_f = RealVectorX::Zero(n2),
2676  strain = RealVectorX::Zero(3),
2677  delta_t = RealVectorX::Zero(1),
2678  delta_t_sens = RealVectorX::Zero(1);
2679 
2680  FEMOperatorMatrix
2681  Bmat_lin,
2682  Bmat_nl_x,
2683  Bmat_nl_y,
2684  Bmat_nl_u,
2685  Bmat_nl_v,
2686  Bmat_bend,
2687  Bmat_vk;
2688 
2689  Bmat_lin.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2690  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2691  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2692  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2693  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2694  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2695  Bmat_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2696 
2697  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2698 
2699  MAST::BendingOperatorType
2700  bending_model = _property.bending_model(_elem);
2701 
2702  std::unique_ptr<MAST::BendingOperator2D> bend;
2703 
2704  if (bending_model != MAST::NO_BENDING)
2705  bend.reset(MAST::build_bending_operator_2D(bending_model,
2706  *this,
2707  fe->get_qpoints()).release());
2708 
2709  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
2710  expansion_A = _property.thermal_expansion_A_matrix(*this),
2711  expansion_B = _property.thermal_expansion_B_matrix(*this);
2712 
2713  // temperature function
2714  const MAST::FieldFunction<Real>
2715  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
2716  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
2717 
2718  Real t, t0, t_sens;
2719 
2720  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2721 
2722  // this is moved inside the domain since
2723  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
2724  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
2725  expansion_A->derivative(p, xyz[qp], _time, material_exp_A_mat_sens);
2726  expansion_B->derivative(p, xyz[qp], _time, material_exp_B_mat_sens);
2727 
2728  // get the temperature function
2729  temp_func(xyz[qp], _time, t);
2730  ref_temp_func(xyz[qp], _time, t0);
2731  delta_t(0) = t-t0;
2732 
2733  // get the temperature function
2734  temp_func(xyz[qp], _time, t);
2735  temp_func.derivative(p, xyz[qp], _time, t_sens);
2736  ref_temp_func(xyz[qp], _time, t0);
2737  delta_t(0) = t-t0;
2738  delta_t_sens(0) = t_sens;
2739 
2740  // now prepare the membrane force sensitivity
2741  vec1_n1 = material_exp_A_mat * delta_t_sens; // [C]{alpha dT/dp} (with membrane strain)
2742  vec2_n1 = material_exp_A_mat_sens * delta_t; // d([C]alpha)/dp (T - T0)} (with membrane strain)
2743  vec1_n1 += vec2_n1;
2744  stress(0,0) = vec1_n1(0); // sigma_xx
2745  stress(0,1) = vec1_n1(2); // sigma_xy
2746  stress(1,0) = vec1_n1(2); // sigma_yx
2747  stress(1,1) = vec1_n1(1); // sigma_yy
2748 
2749  vec2_n1 = material_exp_B_mat * delta_t_sens; // [C]{alpha dT/dp} (with bending strain)
2750  vec5_n1 = material_exp_B_mat_sens * delta_t; // d([C] alpha)/dp (T - T0) (with bending strain)
2751  vec2_n1 += vec5_n1;
2752 
2753 
2754  this->initialize_green_lagrange_strain_operator(qp,
2755  *fe,
2756  _local_sol,
2757  strain,
2758  mat_x,
2759  mat_y,
2760  Bmat_lin,
2761  Bmat_nl_x,
2762  Bmat_nl_y,
2763  Bmat_nl_u,
2764  Bmat_nl_v);
2765 
2766  // membrane strain
2767  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
2768  local_f += JxW[qp] * vec3_n2;
2769 
2770  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
2771 
2772  // nonlinear strain operotor
2773  // x
2774  vec4_2 = mat_x.transpose() * vec1_n1;
2775  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
2776  local_f.topRows(n2) += JxW[qp] * vec3_n2;
2777 
2778  // y
2779  vec4_2 = mat_y.transpose() * vec1_n1;
2780  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
2781  local_f.topRows(n2) += JxW[qp] * vec3_n2;
2782  }
2783 
2784  if (bend.get()) {
2785  // bending strain
2786  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2787  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
2788  local_f += JxW[qp] * vec3_n2;
2789 
2790  // von Karman strain
2791  if (if_vk) {
2792  // get the vonKarman strain operator if needed
2793  this->initialize_von_karman_strain_operator(qp,
2794  *fe,
2795  vec2_n1, // epsilon_vk
2796  vk_dwdxi_mat,
2797  Bmat_vk);
2798  // von Karman strain
2799  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
2800  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
2801  local_f += JxW[qp] * vec3_n2;
2802  }
2803  }
2804 
2805  if (request_jacobian &&
2806  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
2807 
2808  // u-disp
2809  Bmat_nl_u.left_multiply(mat3, stress);
2810  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
2811  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
2812 
2813  // v-disp
2814  Bmat_nl_v.left_multiply(mat3, stress);
2815  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
2816  local_jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2;
2817  }
2818 
2819  if (request_jacobian && if_vk) { // Jacobian only for vk strain
2820 
2821  // vk - vk
2822  mat3 = RealMatrixX::Zero(2, n2);
2823  Bmat_vk.left_multiply(mat3, stress);
2824  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
2825  local_jac += JxW[qp] * mat2_n2n2;
2826  }
2827  }
2828 
2829 
2830  // now transform to the global coorodinate system
2831  transform_vector_to_global_system(local_f, vec3_n2);
2832  f -= vec3_n2;
2833  if (request_jacobian && if_vk) {
2834  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2835  jac -= mat2_n2n2;
2836  }
2837 
2838  // Jacobian contribution from von Karman strain
2839  return request_jacobian;
2840 }
2841 
2842 
2843 
2844 
2845 void
2846 MAST::StructuralElement2D::
2847 thermal_residual_boundary_velocity(const MAST::FunctionBase& p,
2848  const unsigned int s,
2849  const MAST::FieldFunction<RealVectorX>& vel_f,
2850  MAST::BoundaryConditionBase& bc,
2851  bool request_jacobian,
2852  RealVectorX& f,
2853  RealMatrixX& jac) {
2854 
2855  // prepare the side finite element
2856  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s,
2857  true,
2858  _property.extra_quadrature_order(_elem)));
2859 
2860  std::vector<Real> JxW_Vn = fe->get_JxW();
2861  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2862  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2863 
2864  const unsigned int
2865  n_phi = (unsigned int)fe->get_phi().size(),
2866  n1 = this->n_direct_strain_components(), n2=6*n_phi,
2867  n3 = this->n_von_karman_strain_components(),
2868  dim = 2;
2869 
2870  RealMatrixX
2871  material_exp_A_mat,
2872  material_exp_B_mat,
2873  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
2874  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2875  mat3 = RealMatrixX::Zero(2, n2),
2876  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
2877  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
2878  stress = RealMatrixX::Zero(2,2),
2879  mat_x = RealMatrixX::Zero(3,2),
2880  mat_y = RealMatrixX::Zero(3,2),
2881  local_jac = RealMatrixX::Zero(n2,n2);
2882 
2883  RealVectorX
2884  vec1_n1 = RealVectorX::Zero(n1),
2885  vec2_n1 = RealVectorX::Zero(n1),
2886  vec3_n2 = RealVectorX::Zero(n2),
2887  vec4_2 = RealVectorX::Zero(2),
2888  vec5_n3 = RealVectorX::Zero(n3),
2889  local_f = RealVectorX::Zero(n2),
2890  delta_t = RealVectorX::Zero(1),
2891  strain = RealVectorX::Zero(3),
2892  vel = RealVectorX::Zero(dim);
2893 
2894  FEMOperatorMatrix
2895  Bmat_lin,
2896  Bmat_nl_x,
2897  Bmat_nl_y,
2898  Bmat_nl_u,
2899  Bmat_nl_v,
2900  Bmat_bend,
2901  Bmat_vk;
2902 
2903  Bmat_lin.reinit (n1, _system.n_vars(), n_phi); // three stress-strain components
2904  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi);
2905  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi);
2906  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi);
2907  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi);
2908  Bmat_bend.reinit(n1, _system.n_vars(), n_phi);
2909  Bmat_vk.reinit (n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2910 
2911  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2912 
2913  MAST::BendingOperatorType
2914  bending_model = _property.bending_model(_elem);
2915 
2916  std::unique_ptr<MAST::BendingOperator2D> bend;
2917 
2918  if (bending_model != MAST::NO_BENDING)
2919  bend.reset(MAST::build_bending_operator_2D(bending_model,
2920  *this,
2921  fe->get_qpoints()).release());
2922 
2923  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
2924  expansion_A = _property.thermal_expansion_A_matrix(*this),
2925  expansion_B = _property.thermal_expansion_B_matrix(*this);
2926 
2927  const MAST::FieldFunction<Real>
2928  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
2929  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
2930 
2931  Real
2932  vn = 0.,
2933  t = 0.,
2934  t0 = 0.;
2935 
2936  // modify the JxW_Vn by multiplying the normal velocity to it
2937  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
2938 
2939  vel_f(xyz[qp], _time, vel);
2940  vn = 0.;
2941  for (unsigned int i=0; i<dim; i++)
2942  vn += vel(i)*face_normals[qp](i);
2943  JxW_Vn[qp] *= vn;
2944  }
2945 
2946  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
2947 
2948  // this is moved inside the domain since
2949  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
2950  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
2951 
2952  // get the temperature function
2953  temp_func(xyz[qp], _time, t);
2954  ref_temp_func(xyz[qp], _time, t0);
2955  delta_t(0) = t-t0;
2956 
2957  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)} (with membrane strain)
2958  vec2_n1 = material_exp_B_mat * delta_t; // [C]{alpha (T - T0)} (with bending strain)
2959  stress(0,0) = vec1_n1(0); // sigma_xx
2960  stress(0,1) = vec1_n1(2); // sigma_xy
2961  stress(1,0) = vec1_n1(2); // sigma_yx
2962  stress(1,1) = vec1_n1(1); // sigma_yy
2963 
2964  this->initialize_green_lagrange_strain_operator(qp,
2965  *fe,
2966  _local_sol,
2967  strain,
2968  mat_x,
2969  mat_y,
2970  Bmat_lin,
2971  Bmat_nl_x,
2972  Bmat_nl_y,
2973  Bmat_nl_u,
2974  Bmat_nl_v);
2975 
2976  // membrane strain
2977  Bmat_lin.vector_mult_transpose(vec3_n2, vec1_n1);
2978  local_f += JxW_Vn[qp] * vec3_n2;
2979 
2980  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
2981 
2982  // nonlinear strain operotor
2983  // x
2984  vec4_2 = mat_x.transpose() * vec1_n1;
2985  Bmat_nl_x.vector_mult_transpose(vec3_n2, vec4_2);
2986  local_f.topRows(n2) += JxW_Vn[qp] * vec3_n2;
2987 
2988  // y
2989  vec4_2 = mat_y.transpose() * vec1_n1;
2990  Bmat_nl_y.vector_mult_transpose(vec3_n2, vec4_2);
2991  local_f.topRows(n2) += JxW_Vn[qp] * vec3_n2;
2992  }
2993 
2994  if (bend.get()) {
2995  // bending strain
2996  bend->initialize_bending_strain_operator(*fe, qp, Bmat_bend);
2997  Bmat_bend.vector_mult_transpose(vec3_n2, vec2_n1);
2998  local_f += JxW_Vn[qp] * vec3_n2;
2999 
3000  // von Karman strain
3001  if (if_vk) {
3002  // get the vonKarman strain operator if needed
3003  this->initialize_von_karman_strain_operator(qp,
3004  *fe,
3005  vec2_n1, // epsilon_vk
3006  vk_dwdxi_mat,
3007  Bmat_vk);
3008  // von Karman strain
3009  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
3010  Bmat_vk.vector_mult_transpose(vec3_n2, vec4_2);
3011  local_f += JxW_Vn[qp] * vec3_n2;
3012  }
3013  }
3014 
3015  if (request_jacobian &&
3016  _property.strain_type() == MAST::NONLINEAR_STRAIN) {
3017 
3018  // u-disp
3019  Bmat_nl_u.left_multiply(mat3, stress);
3020  Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3);
3021  local_jac.topLeftCorner(n2, n2) += JxW_Vn[qp] * mat2_n2n2;
3022 
3023  // v-disp
3024  Bmat_nl_v.left_multiply(mat3, stress);
3025  Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3);
3026  local_jac.topLeftCorner(n2, n2) += JxW_Vn[qp] * mat2_n2n2;
3027  }
3028 
3029  if (request_jacobian && if_vk) { // Jacobian only for vk strain
3030  // vk - vk
3031  mat3 = RealMatrixX::Zero(2, n2);
3032  Bmat_vk.left_multiply(mat3, stress);
3033  Bmat_vk.right_multiply_transpose(mat2_n2n2, mat3);
3034  local_jac += JxW_Vn[qp] * mat2_n2n2;
3035  }
3036  }
3037 
3038  // now transform to the global coorodinate system
3039  transform_vector_to_global_system(local_f, vec3_n2);
3040  f -= vec3_n2;
3041  if (request_jacobian && if_vk) {
3042  transform_matrix_to_global_system(local_jac, mat2_n2n2);
3043  jac -= mat2_n2n2;
3044  }
3045 }
3046 
3047 
3048 
3049 void
3050 MAST::StructuralElement2D::
3051 thermal_residual_temperature_derivative(const MAST::FEBase& fe_thermal,
3052  RealMatrixX& m) {
3053 
3054 
3055  const std::vector<std::vector<Real>>& phi_temp = fe_thermal.get_phi();
3056  const std::vector<Real>& JxW = fe_thermal.get_JxW();
3057  const std::vector<libMesh::Point>& xyz = fe_thermal.get_xyz();
3058 
3059  std::unique_ptr<MAST::FEBase> fe_str(_elem.init_fe(true, false,
3060  _property.extra_quadrature_order(_elem)));
3061  libmesh_assert(fe_str->get_fe_type() == fe_thermal.get_fe_type());
3062  // this is a weak assertion. We really want that the same qpoints be used,
3063  // but we are assuming that if the elem type is the same, and if the
3064  // number fo qpoints is the same, then the location of qpoints will also
3065  // be the same.
3066  libmesh_assert_equal_to(fe_str->get_qpoints().size(), fe_thermal.get_qpoints().size());
3067 
3068 
3069  const unsigned int
3070  n_phi_str = (unsigned int)fe_str->get_phi().size(),
3071  nt = (unsigned int)fe_thermal.get_phi().size(),
3072  n1 = this->n_direct_strain_components(),
3073  n2 = 6*n_phi_str,
3074  n3 = this->n_von_karman_strain_components();
3075 
3076  RealMatrixX
3077  material_exp_A_mat,
3078  material_exp_B_mat,
3079  mat1_n1nt = RealMatrixX::Zero(n1,nt),
3080  mat2_n1nt = RealMatrixX::Zero(n1,nt),
3081  mat3_n2nt = RealMatrixX::Zero(n2,nt),
3082  mat5,
3083  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
3084  stress = RealMatrixX::Zero(2,2),
3085  mat_x = RealMatrixX::Zero(3,2),
3086  mat_y = RealMatrixX::Zero(3,2),
3087  Bmat_temp = RealMatrixX::Zero(1,nt),
3088  local_m = RealMatrixX::Zero(n2,nt);
3089 
3090  RealVectorX
3091  phi = RealVectorX::Zero(nt),
3092  vec_n1 = RealVectorX::Zero(n1),
3093  strain = RealVectorX::Zero(3);
3094 
3095  FEMOperatorMatrix
3096  Bmat_lin,
3097  Bmat_nl_x,
3098  Bmat_nl_y,
3099  Bmat_nl_u,
3100  Bmat_nl_v,
3101  Bmat_bend,
3102  Bmat_vk;
3103 
3104  Bmat_lin.reinit(n1, _system.n_vars(), n_phi_str); // three stress-strain components
3105  Bmat_nl_x.reinit(2, _system.n_vars(), n_phi_str);
3106  Bmat_nl_y.reinit(2, _system.n_vars(), n_phi_str);
3107  Bmat_nl_u.reinit(2, _system.n_vars(), n_phi_str);
3108  Bmat_nl_v.reinit(2, _system.n_vars(), n_phi_str);
3109  Bmat_bend.reinit(n1, _system.n_vars(), n_phi_str);
3110  Bmat_vk.reinit(n3, _system.n_vars(), n_phi_str); // only dw/dx and dw/dy
3111 
3112  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
3113 
3114  MAST::BendingOperatorType
3115  bending_model = _property.bending_model(_elem);
3116 
3117  std::unique_ptr<MAST::BendingOperator2D> bend;
3118 
3119  if (bending_model != MAST::NO_BENDING)
3120  bend.reset(MAST::build_bending_operator_2D(bending_model,
3121  *this,
3122  fe_str->get_qpoints()).release());
3123 
3124  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
3125  expansion_A = _property.thermal_expansion_A_matrix(*this),
3126  expansion_B = _property.thermal_expansion_B_matrix(*this);
3127 
3128  for (unsigned int qp=0; qp<JxW.size(); qp++) {
3129 
3130  // this is moved inside the domain since
3131  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
3132  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
3133 
3134  // shape function values for the temperature FE
3135  for ( unsigned int i_nd=0; i_nd<nt; i_nd++ )
3136  Bmat_temp(0, i_nd) = phi_temp[i_nd][qp];
3137 
3138  this->initialize_green_lagrange_strain_operator(qp,
3139  *fe_str,
3140  _local_sol,
3141  strain,
3142  mat_x,
3143  mat_y,
3144  Bmat_lin,
3145  Bmat_nl_x,
3146  Bmat_nl_y,
3147  Bmat_nl_u,
3148  Bmat_nl_v);
3149 
3150  mat1_n1nt = material_exp_A_mat * Bmat_temp;
3151  mat2_n1nt = material_exp_B_mat * Bmat_temp;
3152  Bmat_lin.right_multiply_transpose(mat3_n2nt, mat1_n1nt);
3153  local_m += JxW[qp] * mat3_n2nt;
3154 
3155 
3156  if (_property.strain_type() == MAST::NONLINEAR_STRAIN) {
3157 
3158  // nonlinear strain operotor
3159  // x
3160  mat5 = mat_x.transpose() * mat1_n1nt;
3161  Bmat_nl_x.right_multiply_transpose(mat3_n2nt, mat5);
3162  local_m.topRows(n2) += JxW[qp] * mat3_n2nt;
3163 
3164  // y
3165  mat5 = mat_y.transpose() * mat1_n1nt;
3166  Bmat_nl_y.right_multiply_transpose(mat3_n2nt, mat5);
3167  local_m.topRows(n2) += JxW[qp] * mat3_n2nt;
3168  }
3169 
3170  if (bend.get()) {
3171  // bending strain
3172  bend->initialize_bending_strain_operator(*fe_str, qp, Bmat_bend);
3173  Bmat_bend.right_multiply_transpose(mat3_n2nt, mat2_n1nt);
3174  local_m += JxW[qp] * mat3_n2nt;
3175 
3176  // von Karman strain
3177  if (if_vk) {
3178  // get the vonKarman strain operator if needed
3179  this->initialize_von_karman_strain_operator(qp,
3180  *fe_str,
3181  vec_n1, // epsilon_vk
3182  vk_dwdxi_mat,
3183  Bmat_vk);
3184  // von Karman strain
3185  mat5 = vk_dwdxi_mat.transpose() * mat1_n1nt;
3186  Bmat_vk.right_multiply_transpose(mat3_n2nt, mat5);
3187  local_m += JxW[qp] * mat3_n2nt;
3188  }
3189  }
3190  }
3191 
3192  mat3_n2nt.setZero();
3193  phi = RealVectorX::Zero(n2);
3194  vec_n1 = RealVectorX::Zero(n2);
3195  for (unsigned int i=0; i<nt; i++) {
3196  phi = local_m.col(i);
3197  transform_vector_to_global_system(phi, vec_n1);
3198  mat3_n2nt.col(i) = vec_n1;
3199  }
3200  m -= mat3_n2nt;
3201 }
3202 
3203 
3204 
3205 
3206 bool
3207 MAST::StructuralElement2D::
3208 piston_theory_residual(bool request_jacobian,
3209  RealVectorX &f,
3210  RealMatrixX& jac_xdot,
3211  RealMatrixX& jac,
3212  const unsigned int side,
3213  MAST::BoundaryConditionBase& bc) {
3214 
3215  libmesh_assert(false); // to be implemented
3216 
3217  return (request_jacobian);
3218 }
3219 
3220 
3221 
3222 
3223 bool
3224 MAST::StructuralElement2D::
3225 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
3226  bool request_jacobian,
3227  RealVectorX &f,
3228  RealMatrixX& jac_xdot,
3229  RealMatrixX& jac,
3230  const unsigned int side,
3231  MAST::BoundaryConditionBase& bc) {
3232 
3233  libmesh_assert(false); // to be implemented
3234 
3235  return (request_jacobian);
3236 }
3237 
3238 
3239 
3240 
3241 
3242 bool
3243 MAST::StructuralElement2D::
3244 piston_theory_residual(bool request_jacobian,
3245  RealVectorX &f,
3246  RealMatrixX& jac_xdot,
3247  RealMatrixX& jac,
3248  MAST::BoundaryConditionBase& bc) {
3249 
3250  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
3251  libmesh_assert(!follower_forces); // not implemented yet for follower forces
3252 
3253  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
3254 
3255  const std::vector<Real> &JxW = fe->get_JxW();
3256  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
3257  const std::vector<std::vector<Real> >& phi = fe->get_phi();
3258  const unsigned int
3259  n_phi = (unsigned int)phi.size(),
3260  n1 = 2,
3261  n2 = _system.n_vars()*n_phi;
3262 
3263 
3264  // normal for face integration
3265  libMesh::Point normal;
3266  // direction of pressure assumed to be normal (along local z-axis)
3267  // to the element face for 2D and along local y-axis for 1D element.
3268  normal(_elem.dim()) = -1.;
3269 
3270 
3271  // convert to piston theory boundary condition so that the necessary
3272  // flow properties can be obtained
3273  const MAST::PistonTheoryBoundaryCondition& piston_bc =
3274  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
3275 
3276 
3277  // create the constant field functions to pass the dwdx and dwdt values
3278  // to the piston theory pressure functions
3279  MAST::Parameter
3280  dwdx_p ("dwdx", 0.),
3281  dwdt_p ("dwdt", 0.);
3282 
3283  MAST::ConstantFieldFunction
3284  dwdx_f ("dwdx", dwdx_p),
3285  dwdt_f ("dwdx", dwdt_p);
3286 
3287 
3288  std::unique_ptr<MAST::FieldFunction<Real> >
3289  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
3290  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
3291  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
3292 
3293  FEMOperatorMatrix
3294  Bmat_w, // operator matrix for the w-displacement
3295  dBmat; // operator matrix to calculate the derivativ of w wrt x and y
3296 
3297  dBmat.reinit(n1, _system.n_vars(), n_phi);
3298 
3299  RealVectorX
3300  phi_vec = RealVectorX::Zero(n_phi),
3301  force = RealVectorX::Zero(n1),
3302  local_f = RealVectorX::Zero(n2),
3303  vec_n1 = RealVectorX::Zero(n1),
3304  vec_n2 = RealVectorX::Zero(n2),
3305  vel_vec = RealVectorX::Zero(3),
3306  dummy = RealVectorX::Zero(3);
3307 
3308  RealMatrixX
3309  dwdx = RealMatrixX::Zero(3,2),
3310  local_jac_xdot = RealMatrixX::Zero(n2,n2),
3311  local_jac = RealMatrixX::Zero(n2,n2),
3312  mat_n2n2 = RealMatrixX::Zero(n2,n2),
3313  mat_n1n2 = RealMatrixX::Zero(n1,n2),
3314  mat_22 = RealMatrixX::Zero(2,2);
3315 
3316  // we need the velocity vector in the local coordinate system so that
3317  // the appropriate component of the w-derivative can be used
3318  vel_vec = _elem.T_matrix().transpose() * piston_bc.vel_vec();
3319 
3320  Real
3321  dwdt_val = 0.,
3322  dwdx_val = 0.,
3323  p_val = 0.;
3324 
3325 
3326  for (unsigned int qp=0; qp<qpoint.size(); qp++)
3327  {
3328 
3329  // now set the shape function values
3330  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
3331  phi_vec(i_nd) = phi[i_nd][qp];
3332 
3333  // initialize the B matrix for only the w-displacement
3334  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
3335  Bmat_w.set_shape_function(0, 2, phi_vec); // interpolates w-displacement
3336 
3337  // use the Bmat to calculate the velocity vector. Only the w-displacement
3338  // is of interest in the local coordinate, since that is the only
3339  // component normal to the surface.
3340  Bmat_w.right_multiply(vec_n1, _local_vel);
3341  dwdt_val = vec_n1(0);
3342 
3343  // get the operators for dw/dx and dw/dy to calculate the
3344  // normal velocity. We will use the von Karman strain operators
3345  // for this
3346  this->initialize_von_karman_strain_operator(qp,
3347  *fe,
3348  dummy,
3349  dwdx,
3350  dBmat);
3351 
3352  // the diagonal of dwdx matrix stores the
3353  dwdx_val = 0.;
3354  for (unsigned int i=0; i<2; i++)
3355  dwdx_val += dwdx(i,i) * vel_vec(i); // (dw/dx_i)*U_inf . n_i
3356 
3357  // calculate the pressure value
3358  dwdx_p = dwdx_val;
3359  dwdt_p = dwdt_val;
3360  (*pressure)(qpoint[qp], _time, p_val);
3361 
3362  // calculate force
3363  force(0) = p_val * normal(2);
3364 
3365 
3366  Bmat_w.vector_mult_transpose(vec_n2, force);
3367  local_f += JxW[qp] * vec_n2;
3368 
3369 
3370  // calculate the Jacobian if requested
3371  if (request_jacobian) {
3372 
3373  // we need the derivative of cp wrt normal velocity
3374  (*dpressure_dxdot)(qpoint[qp], _time, p_val);
3375 
3376  // calculate the component of Jacobian due to w-velocity
3377  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
3378  local_jac_xdot += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3379 
3380  // now calculate the component of Jacobian
3381  (*dpressure_dx)(qpoint[qp], _time, p_val);
3382 
3383  // derivative wrt x
3384  mat_22.setZero(2,2);
3385  mat_22(0,0) = vel_vec(0);
3386  dBmat.left_multiply(mat_n1n2, mat_22);
3387  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dx
3388  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3389 
3390  // derivative wrt y
3391  mat_22.setZero(2,2);
3392  mat_22(1,1) = vel_vec(1);
3393  dBmat.left_multiply(mat_n1n2, mat_22);
3394  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dy
3395  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3396  }
3397  }
3398 
3399 
3400  // now transform to the global system and add
3401  transform_vector_to_global_system(local_f, vec_n2);
3402  f -= vec_n2;
3403 
3404  // if the Jacobian was requested, then transform it and add to the
3405  // global Jacobian
3406  if (request_jacobian) {
3407  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
3408  jac_xdot -= mat_n2n2;
3409 
3410  transform_matrix_to_global_system(local_jac, mat_n2n2);
3411  jac -= mat_n2n2;
3412  }
3413 
3414  return request_jacobian;
3415 }
3416 
3417 
3418 
3419 
3420 bool
3421 MAST::StructuralElement2D::
3422 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
3423  bool request_jacobian,
3424  RealVectorX &f,
3425  RealMatrixX& jac_xdot,
3426  RealMatrixX& jac,
3427  MAST::BoundaryConditionBase& bc) {
3428 
3429 
3430  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
3431  libmesh_assert(!follower_forces); // not implemented yet for follower forces
3432 
3433  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
3434 
3435  const std::vector<Real> &JxW = fe->get_JxW();
3436  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
3437  const std::vector<std::vector<Real> >& phi = fe->get_phi();
3438  const unsigned int
3439  n_phi = (unsigned int)phi.size(),
3440  n1 = 2,
3441  n2 = _system.n_vars()*n_phi;
3442 
3443 
3444  // normal for face integration
3445  libMesh::Point normal;
3446  // direction of pressure assumed to be normal (along local z-axis)
3447  // to the element face for 2D and along local y-axis for 1D element.
3448  normal(_elem.dim()) = -1.;
3449 
3450 
3451  // convert to piston theory boundary condition so that the necessary
3452  // flow properties can be obtained
3453  const MAST::PistonTheoryBoundaryCondition& piston_bc =
3454  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
3455 
3456 
3457  // create the constant field functions to pass the dwdx and dwdt values
3458  // to the piston theory pressure functions
3459  MAST::Parameter
3460  dwdx_p ("dwdx", 0.),
3461  dwdt_p ("dwdt", 0.);
3462 
3463  MAST::ConstantFieldFunction
3464  dwdx_f ("dwdx", dwdx_p),
3465  dwdt_f ("dwdx", dwdt_p);
3466 
3467 
3468  std::unique_ptr<MAST::FieldFunction<Real> >
3469  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
3470  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
3471  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
3472 
3473  FEMOperatorMatrix
3474  Bmat_w, // operator matrix for the w-displacement
3475  dBmat; // operator matrix to calculate the derivativ of w wrt x and y
3476 
3477  dBmat.reinit(n1, _system.n_vars(), n_phi);
3478 
3479  RealVectorX
3480  phi_vec = RealVectorX::Zero(n_phi),
3481  force = RealVectorX::Zero(n1),
3482  local_f = RealVectorX::Zero(n2),
3483  vec_n1 = RealVectorX::Zero(n1),
3484  vec_n2 = RealVectorX::Zero(n2),
3485  vel_vec = RealVectorX::Zero(3),
3486  dummy = RealVectorX::Zero(3);
3487 
3488  RealMatrixX
3489  dwdx = RealMatrixX::Zero(3,2),
3490  local_jac_xdot = RealMatrixX::Zero(n2,n2),
3491  local_jac = RealMatrixX::Zero(n2,n2),
3492  mat_n2n2 = RealMatrixX::Zero(n2,n2),
3493  mat_n1n2 = RealMatrixX::Zero(n1,n2),
3494  mat_22 = RealMatrixX::Zero(2,2);
3495 
3496  // we need the velocity vector in the local coordinate system so that
3497  // the appropriate component of the w-derivative can be used
3498  vel_vec = _elem.T_matrix().transpose() * piston_bc.vel_vec();
3499 
3500  Real
3501  dwdt_val = 0.,
3502  dwdx_val = 0.,
3503  p_val = 0.;
3504 
3505 
3506  for (unsigned int qp=0; qp<qpoint.size(); qp++)
3507  {
3508 
3509  // now set the shape function values
3510  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
3511  phi_vec(i_nd) = phi[i_nd][qp];
3512 
3513  // initialize the B matrix for only the w-displacement
3514  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
3515  Bmat_w.set_shape_function(0, 2, phi_vec); // interpolates w-displacement
3516 
3517  // use the Bmat to calculate the velocity vector. Only the w-displacement
3518  // is of interest in the local coordinate, since that is the only
3519  // component normal to the surface.
3520  Bmat_w.right_multiply(vec_n1, _local_vel);
3521  dwdt_val = vec_n1(0);
3522 
3523  // get the operators for dw/dx and dw/dy to calculate the
3524  // normal velocity. We will use the von Karman strain operators
3525  // for this
3526  this->initialize_von_karman_strain_operator(qp,
3527  *fe,
3528  dummy,
3529  dwdx,
3530  dBmat);
3531 
3532  // the diagonal of dwdx matrix stores the
3533  dwdx_val = 0.;
3534  for (unsigned int i=0; i<2; i++)
3535  dwdx_val += dwdx(i,i) * vel_vec(i); // (dw/dx_i)*U_inf . n_i
3536 
3537  // calculate the pressure value
3538  dwdx_p = dwdx_val;
3539  dwdt_p = dwdt_val;
3540  pressure->derivative(p, qpoint[qp], _time, p_val);
3541 
3542  // calculate force
3543  force(0) = p_val * normal(2);
3544 
3545 
3546  Bmat_w.vector_mult_transpose(vec_n2, force);
3547  local_f += JxW[qp] * vec_n2;
3548 
3549 
3550  // calculate the Jacobian if requested
3551  if (request_jacobian) {
3552 
3553  // we need the derivative of cp wrt normal velocity
3554  dpressure_dxdot->derivative(p, qpoint[qp], _time, p_val);
3555 
3556  // calculate the component of Jacobian due to w-velocity
3557  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
3558  local_jac_xdot += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3559 
3560  // now calculate the component of Jacobian
3561  dpressure_dx->derivative(p, qpoint[qp], _time, p_val);
3562 
3563  // derivative wrt x
3564  mat_22.setZero(2,2);
3565  mat_22(0,0) = vel_vec(0);
3566  dBmat.left_multiply(mat_n1n2, mat_22);
3567  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dx
3568  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3569 
3570  // derivative wrt y
3571  mat_22.setZero(2,2);
3572  mat_22(1,1) = vel_vec(1);
3573  dBmat.left_multiply(mat_n1n2, mat_22);
3574  Bmat_w.right_multiply_transpose(mat_n2n2, mat_n1n2); // v: B^T dB/dy
3575  local_jac += (JxW[qp] * p_val * normal(2)) * mat_n2n2;
3576  }
3577  }
3578 
3579 
3580  // now transform to the global system and add
3581  transform_vector_to_global_system(local_f, vec_n2);
3582  f -= vec_n2;
3583 
3584  // if the Jacobian was requested, then transform it and add to the
3585  // global Jacobian
3586  if (request_jacobian) {
3587  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
3588  jac_xdot -= mat_n2n2;
3589 
3590  transform_matrix_to_global_system(local_jac, mat_n2n2);
3591  jac -= mat_n2n2;
3592  }
3593 
3594 
3595  return request_jacobian;
3596 }
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:253
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
MAST::ElementPropertyCardBase::thermal_expansion_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > thermal_expansion_B_matrix(const MAST::ElementBase &e) const =0
MAST::ElementBase::_elem
const MAST::GeomElem & _elem
geometric element for which the computations are performed
Definition: elem_base.h:218
MAST::StressStrainOutputBase
Data structure provides the mechanism to store stress and strain output from a structural analysis...
Definition: stress_output_base.h:61
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:57
MAST::ElementPropertyCardBase::prestress_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > prestress_B_matrix(MAST::ElementBase &e) const =0
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:1693
MAST::FEBase::get_JxW
virtual const std::vector< Real > & get_JxW() const
Definition: fe_base.cpp:220
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::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:735
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::StructuralElement2D::StructuralElement2D
StructuralElement2D(MAST::SystemInitialization &sys, MAST::AssemblyBase &assembly, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
Definition: structural_element_2d.cpp:38
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:2131
MAST::FEBase::get_qpoints
virtual const std::vector< libMesh::Point > & get_qpoints() const
Definition: fe_base.cpp:407
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:1498
fem_operator_matrix.h
bending_operator.h
MAST::ElementPropertyCardBase::stiffness_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > stiffness_B_matrix(const MAST::ElementBase &e) const =0
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:719
MAST::ElementBase::_system
MAST::SystemInitialization & _system
SystemInitialization object associated with this element.
Definition: elem_base.h:208
MAST::FEMOperatorMatrix::m
unsigned int m() const
Definition: fem_operator_matrix.h:50
MAST::FEMOperatorMatrix::n
unsigned int n() const
Definition: fem_operator_matrix.h:52
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:3051
stress_output_base.h
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::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::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:1002
MAST::StructuralElement2D::calculate_stress_temperature_derivative
virtual void calculate_stress_temperature_derivative(MAST::FEBase &fe_thermal, MAST::StressStrainOutputBase &output)
Definition: structural_element_2d.cpp:884
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:718
MAST::FunctionBase::is_shape_parameter
virtual bool is_shape_parameter() const
Definition: function_base.h:89
MAST::ElementPropertyCardBase::get_material
virtual const MAST::MaterialPropertyCardBase & get_material() const
return the material property.
Definition: element_property_card_base.h:126
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:3244
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
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:1531
MAST::StructuralElementBase::follower_forces
bool follower_forces
flag for follower forces
Definition: structural_element_base.h:440
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:2847
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::GeomElem::T_matrix
const RealMatrixX & T_matrix() const
O.
Definition: geom_elem.cpp:228
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:149
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:2635
MAST::FEMOperatorMatrix::vector_mult
void vector_mult(T &res, const T &v) const
res = [this] * v
Definition: fem_operator_matrix.h:280
MAST::StressStrainOutputBase::get_thermal_load_for_elem
MAST::BoundaryConditionBase * get_thermal_load_for_elem(const MAST::GeomElem &elem)
Definition: stress_output_base.cpp:787
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:1980
MAST::ElementPropertyCardBase::stiffness_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > stiffness_A_matrix(const MAST::ElementBase &e) const =0
RealVectorX
Matrix< Real, Dynamic, 1 > RealVectorX
Definition: mast_data_types.h:35
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
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:701
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::init_side_fe
virtual std::unique_ptr< MAST::FEBase > init_side_fe(unsigned int s, bool init_grads, int extra_quadrature_order=0) const
initializes the finite element shape function and quadrature object for the side with the order of qu...
Definition: geom_elem.cpp:159
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::FEBase::get_dphi
virtual const std::vector< std::vector< libMesh::RealVectorValue > > & get_dphi() const
Definition: fe_base.cpp:255
MAST::GeomElem::dim
unsigned int dim() const
Definition: geom_elem.cpp:91
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:2284
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:1160
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:1951
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:87
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
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::FEMOperatorMatrix::left_multiply
void left_multiply(T &r, const T &m) const
[R] = [M] * [this]
Definition: fem_operator_matrix.h:415
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:768
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:3422
MAST::FEBase::get_phi
virtual const std::vector< std::vector< Real > > & get_phi() const
Definition: fe_base.cpp:247
MAST::StructuralElementBase::_local_sol
RealVectorX _local_sol
local solution
Definition: structural_element_base.h:707
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:2359
MAST::StructuralElementBase::_local_vel
RealVectorX _local_vel
local velocity
Definition: structural_element_base.h:731
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
piston_theory_boundary_condition.h
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:642
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:583
MAST::BendingOperatorType
BendingOperatorType
Definition: bending_operator.h:48
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:142
MAST::ElementBase::_time
const Real & _time
time for which system is being assembled
Definition: elem_base.h:232
MAST::FEBase::get_fe_type
libMesh::FEType get_fe_type() const
Definition: fe_base.cpp:212
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:1329
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:1607
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:1965
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::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:2439
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::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:64
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::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:122
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:630