MAST
structural_element_1d.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_1d.h"
22 #include "elasticity/piston_theory_boundary_condition.h"
23 #include "elasticity/stress_output_base.h"
24 #include "elasticity/bending_operator.h"
25 #include "numerics/fem_operator_matrix.h"
26 #include "property_cards/element_property_card_1D.h"
27 #include "property_cards/material_property_card_base.h"
28 #include "base/system_initialization.h"
29 #include "base/boundary_condition_base.h"
30 #include "base/parameter.h"
31 #include "base/constant_field_function.h"
32 #include "base/assembly_base.h"
33 #include "mesh/fe_base.h"
34 #include "mesh/geom_elem.h"
35 
36 
37 MAST::StructuralElement1D::StructuralElement1D(MAST::SystemInitialization& sys,
38  MAST::AssemblyBase& assembly,
39  const MAST::GeomElem& elem,
40  const MAST::ElementPropertyCardBase& p):
41 MAST::BendingStructuralElem(sys, assembly, elem, p) {
42 
43 }
44 
45 
46 
47 MAST::StructuralElement1D::~StructuralElement1D() {
48 
49 }
50 
51 
52 
53 
54 void
55 MAST::StructuralElement1D::
56 initialize_direct_strain_operator(const unsigned int qp,
57  const MAST::FEBase& fe,
58  MAST::FEMOperatorMatrix& Bmat) {
59 
60  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
61 
62  unsigned int n_phi = (unsigned int)dphi.size();
63  RealVectorX phi = RealVectorX::Zero(n_phi);
64 
65  libmesh_assert_equal_to(Bmat.m(), 2);
66  libmesh_assert_equal_to(Bmat.n(), 6*n_phi);
67  libmesh_assert_less (qp, dphi[0].size());
68 
69  // now set the shape function values
70  // dN/dx
71  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
72  phi(i_nd) = dphi[i_nd][qp](0);
73  Bmat.set_shape_function(0, 0, phi); // epsilon_xx = du/dx
74  Bmat.set_shape_function(1, 3, phi); // torsion operator = dtheta_x/dx
75 }
76 
77 
78 
79 void
80 MAST::StructuralElement1D::
81 initialize_von_karman_strain_operator(const unsigned int qp,
82  const MAST::FEBase& fe,
83  RealVectorX& vk_strain,
84  RealMatrixX& vk_dvdxi_mat,
85  RealMatrixX& vk_dwdxi_mat,
86  MAST::FEMOperatorMatrix& Bmat_v_vk,
87  MAST::FEMOperatorMatrix& Bmat_w_vk) {
88 
89  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
90  const unsigned int n_phi = (unsigned int)dphi.size();
91 
92  libmesh_assert_equal_to(vk_strain.size(), 2);
93  libmesh_assert_equal_to(vk_dvdxi_mat.rows(), 2);
94  libmesh_assert_equal_to(vk_dvdxi_mat.cols(), 2);
95  libmesh_assert_equal_to(Bmat_v_vk.m(), 2);
96  libmesh_assert_equal_to(Bmat_v_vk.n(), 6*n_phi);
97  libmesh_assert_equal_to(vk_dwdxi_mat.rows(), 2);
98  libmesh_assert_equal_to(vk_dwdxi_mat.cols(), 2);
99  libmesh_assert_equal_to(Bmat_w_vk.m(), 2);
100  libmesh_assert_equal_to(Bmat_w_vk.n(), 6*n_phi);
101  libmesh_assert_less (qp, dphi[0].size());
102 
103  Real dv=0., dw=0.;
104  vk_strain.setZero();
105  vk_dvdxi_mat.setZero();
106  vk_dwdxi_mat.setZero();
107 
108  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
109 
110  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
111  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
112  dv += phi_vec(i_nd)*_local_sol(n_phi+i_nd); // dv/dx
113  dw += phi_vec(i_nd)*_local_sol(2*n_phi+i_nd); // dw/dx
114  }
115 
116  Bmat_v_vk.set_shape_function(0, 1, phi_vec); // dv/dx
117  Bmat_w_vk.set_shape_function(0, 2, phi_vec); // dw/dx
118  vk_dvdxi_mat(0, 0) = dv; // epsilon-xx : dv/dx
119  vk_dwdxi_mat(0, 0) = dw; // epsilon-xx : dw/dx
120  vk_strain(0) = 0.5*(dv*dv+dw*dw); // 1/2 * [(dv/dx)^2 + (dw/dx)^2]
121 }
122 
123 
124 
125 
126 void
127 MAST::StructuralElement1D::
128 initialize_von_karman_strain_operator_sensitivity(const unsigned int qp,
129  const MAST::FEBase& fe,
130  RealMatrixX& vk_dvdxi_mat_sens,
131  RealMatrixX& vk_dwdxi_mat_sens) {
132 
133  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
134  const unsigned int n_phi = (unsigned int)dphi.size();
135 
136  libmesh_assert_equal_to(vk_dvdxi_mat_sens.rows(), 2);
137  libmesh_assert_equal_to(vk_dvdxi_mat_sens.cols(), 2);
138  libmesh_assert_equal_to(vk_dwdxi_mat_sens.rows(), 2);
139  libmesh_assert_equal_to(vk_dwdxi_mat_sens.cols(), 2);
140  libmesh_assert_less (qp, dphi[0].size());
141 
142  Real dv=0., dw=0.;
143  vk_dvdxi_mat_sens.setZero();
144  vk_dwdxi_mat_sens.setZero();
145 
146  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
147 
148  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ ) {
149  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
150  dv += phi_vec(i_nd)*_local_sol_sens(n_phi+i_nd); // dv/dx
151  dw += phi_vec(i_nd)*_local_sol_sens(2*n_phi+i_nd); // dw/dx
152  }
153 
154  vk_dvdxi_mat_sens(0, 0) = dv; // epsilon-xx : dv/dx
155  vk_dwdxi_mat_sens(0, 0) = dw; // epsilon-xx : dw/dx
156 }
157 
158 
159 
160 
161 bool
162 MAST::StructuralElement1D::calculate_stress(bool request_derivative,
163  const MAST::FunctionBase* p,
164  MAST::StressStrainOutputBase& output) {
165 
166  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
167 
168  const unsigned int
169  qp_loc_fe_size = (unsigned int)fe->get_qpoints().size(),
170  n_added_qp = 4;
171 
172  std::vector<libMesh::Point>
173  qp_loc_fe = fe->get_qpoints(),
174  qp_loc(qp_loc_fe_size*n_added_qp);
175 
176 
177  // we will evaluate the stress at upper and lower layers of element,
178  // so we will add two new points for each qp_loc
179  // TODO: move this to element section property class for composite materials
180  for (unsigned int i=0; i<qp_loc_fe.size(); i++) {
181 
182  qp_loc[i*4] = qp_loc_fe[i];
183  qp_loc[i*4](1) = +1.;
184  qp_loc[i*4](2) = +1.; // upper right
185 
186  qp_loc[i*4+1] = qp_loc_fe[i];
187  qp_loc[i*4+1](1) = -1.;
188  qp_loc[i*4+1](2) = +1.; // upper left
189 
190  qp_loc[i*4+2] = qp_loc_fe[i];
191  qp_loc[i*4+2](1) = +1.;
192  qp_loc[i*4+2](2) = -1.; // lower right
193 
194  qp_loc[i*4+3] = qp_loc_fe[i];
195  qp_loc[i*4+3](1) = -1.;
196  qp_loc[i*4+3](2) = -1.; // lower left
197  }
198 
199 
200  MAST::BendingOperatorType bending_model = _property.bending_model(_elem);
201 
202  std::unique_ptr<MAST::BendingOperator1D>
203  bend(MAST::build_bending_operator_1D(bending_model,
204  *this,
205  qp_loc_fe));
206 
207 
208  // now that the FE object has been initialized, evaluate the stress values
209 
210 
211  const std::vector<Real> &JxW = fe->get_JxW();
212  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
213  const unsigned int
214  n_phi = (unsigned int)fe->n_shape_functions(),
215  n1 = this->n_direct_strain_components(),
216  n2 = 6*n_phi,
217  n3 = this->n_von_karman_strain_components();
218 
219  Real
220  y = 0.,
221  z = 0.,
222  y_off = 0.,
223  z_off = 0.,
224  temp = 0.,
225  ref_t = 0.,
226  alpha = 0.,
227  dtemp = 0.,
228  dref_t= 0.,
229  dalpha= 0.;
230 
231  RealMatrixX
232  material_mat,
233  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
234  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
235  dstrain_dX = RealMatrixX::Zero(n1,n2),
236  dstress_dX = RealMatrixX::Zero(n1,n2),
237  mat_n1n2 = RealMatrixX::Zero(n1,n2),
238  eye = RealMatrixX::Identity(n1, n1),
239  dstrain_dX_3D= RealMatrixX::Zero(6,n2),
240  dstress_dX_3D= RealMatrixX::Zero(6,n2);
241 
242  RealVectorX
243  strain = RealVectorX::Zero(n1),
244  stress = RealVectorX::Zero(n1),
245  strain_bend = RealVectorX::Zero(n1),
246  strain_vk = RealVectorX::Zero(n3),
247  strain_3D = RealVectorX::Zero(6),
248  stress_3D = RealVectorX::Zero(6),
249  dstrain_dp = RealVectorX::Zero(n1),
250  dstress_dp = RealVectorX::Zero(n1);
251 
252  MAST::FEMOperatorMatrix
253  Bmat_mem,
254  Bmat_bend_v,
255  Bmat_bend_w,
256  Bmat_v_vk,
257  Bmat_w_vk;
258 
259  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
260  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
261  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
262  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
263  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
264 
265  // TODO: remove this const-cast, which may need change in API of
266  // material card
267  const MAST::FieldFunction<RealMatrixX >&
268  mat_stiff =
269  const_cast<MAST::MaterialPropertyCardBase&>(_property.get_material()).stiffness_matrix(1);
270 
271  // get the thickness values for the bending strain calculation
272  const MAST::FieldFunction<Real>
273  &hy = _property.get<MAST::FieldFunction<Real> >("hy"),
274  &hz = _property.get<MAST::FieldFunction<Real> >("hz"),
275  &hy_off = _property.get<MAST::FieldFunction<Real> >("hy_off"),
276  &hz_off = _property.get<MAST::FieldFunction<Real> >("hz_off");
277 
278 
279  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN),
280  if_bending = (_property.bending_model(_elem) != MAST::NO_BENDING);
281 
282  // a reference to the stress output data structure
283  MAST::StressStrainOutputBase& stress_output =
284  dynamic_cast<MAST::StressStrainOutputBase&>(output);
285 
286  // check to see if the element has any thermal loads specified
287  // The object returns null
288  MAST::BoundaryConditionBase *thermal_load =
289  stress_output.get_thermal_load_for_elem(_elem);
290 
291  const MAST::FieldFunction<Real>
292  *temp_func = nullptr,
293  *ref_temp_func = nullptr,
294  *alpha_func = nullptr;
295 
296  // get pointers to the temperature, if thermal load is specified
297  if (thermal_load) {
298  temp_func =
299  &(thermal_load->get<MAST::FieldFunction<Real> >("temperature"));
300  ref_temp_func =
301  &(thermal_load->get<MAST::FieldFunction<Real> >("ref_temperature"));
302  alpha_func =
303  &(_property.get_material().get<MAST::FieldFunction<Real> >("alpha_expansion"));
304  }
305 
306  // TODO: improve the stress calculation by including shear stress due
307  // to torsion and shear flow due to beam bending.
308 
310  // second for loop to calculate the residual and stiffness contributions
312  unsigned int
313  qp = 0;
314  for (unsigned int qp_loc_index=0; qp_loc_index<qp_loc_fe_size; qp_loc_index++)
315  for (unsigned int section_qp_index=0; section_qp_index<n_added_qp; section_qp_index++)
316  {
317  qp = qp_loc_index*n_added_qp + section_qp_index;
318 
319  // get the material matrix
320  mat_stiff(xyz[qp_loc_index], _time, material_mat);
321 
322  this->initialize_direct_strain_operator(qp_loc_index, *fe, Bmat_mem);
323 
324  // first handle constant throught the thickness stresses: membrane and vonKarman
325  Bmat_mem.vector_mult(strain, _local_sol);
326 
327  // if thermal load was specified, then set the thermal strain
328  // component of the total strain
329  if (thermal_load) {
330  (*temp_func) (xyz[qp_loc_index], _time, temp);
331  (*ref_temp_func)(xyz[qp_loc_index], _time, ref_t);
332  (*alpha_func) (xyz[qp_loc_index], _time, alpha);
333  strain(0) -= alpha*(temp-ref_t);
334  }
335 
336 
337  if (if_bending) {
338  // von Karman strain
339  if (if_vk) { // get the vonKarman strain operator if needed
340 
341  this->initialize_von_karman_strain_operator(qp_loc_index,
342  *fe,
343  strain_vk,
344  vk_dvdxi_mat,
345  vk_dwdxi_mat,
346  Bmat_v_vk,
347  Bmat_w_vk);
348  strain += strain_vk;
349  }
350 
351  // add to this the bending strain
352  hy (xyz[qp_loc_index], _time, y);
353  hz (xyz[qp_loc_index], _time, z);
354  hy_off(xyz[qp_loc_index], _time, y_off);
355  hz_off(xyz[qp_loc_index], _time, z_off);
356 
357  // TODO: this assumes isotropic section. Multilayered sections need
358  // special considerations
359  // This operator depends on the y and z thickness values. Sensitivity
360  // analysis should include the sensitivity of this operator on
361  // these thickness values
362  bend->initialize_bending_strain_operator_for_yz(*fe,
363  qp_loc_index,
364  qp_loc[qp](1) * y/2.+y_off,
365  qp_loc[qp](2) * z/2.+z_off,
366  Bmat_bend_v,
367  Bmat_bend_w);
368  Bmat_bend_v.vector_mult(strain_bend, _local_sol);
369  strain += strain_bend;
370 
371  Bmat_bend_w.vector_mult(strain_bend, _local_sol);
372  strain += strain_bend;
373  }
374 
375 
376  // note that this assumes linear material laws
377  stress = material_mat * strain;
378 
379  // now set the data for the 3D stress-strain vector
380  // this is using only the direct strain/stress.
381  // this can be improved by estimating the shear stresses from
382  // torsion and shear flow from bending.
383  strain_3D(0) = strain(0);
384  stress_3D(0) = stress(0);
385 
386  // set the stress and strain data
387  MAST::StressStrainOutputBase::Data*
388  data = nullptr;
389 
390  // if neither the derivative nor sensitivity is requested, then
391  // we assume that a new data entry is to be provided. Otherwise,
392  // we assume that the stress at this quantity already
393  // exists, and we only need to append sensitivity/derivative
394  // data to it
395  if (!request_derivative && !p)
396  data = &(stress_output.add_stress_strain_at_qp_location(_elem,
397  qp,
398  qp_loc[qp],
399  xyz[qp_loc_index],
400  stress_3D,
401  strain_3D,
402  JxW[qp_loc_index]));
403  else
404  data = &(stress_output.get_stress_strain_data_for_elem_at_qp(_elem,
405  qp));
406 
407  // calculate the derivative if requested
408  if (request_derivative || p) {
409 
410  Bmat_mem.left_multiply(dstrain_dX, eye); // membrane strain is linear
411 
412  if (if_bending) {
413 
414  // von Karman strain
415  if (if_vk) {
416 
417  Bmat_v_vk.left_multiply(mat_n1n2, vk_dvdxi_mat);
418  dstrain_dX += mat_n1n2;
419 
420  Bmat_w_vk.left_multiply(mat_n1n2, vk_dwdxi_mat);
421  dstrain_dX += mat_n1n2;
422  }
423 
424  // bending strain
425  Bmat_bend_v.left_multiply(mat_n1n2, eye);
426  dstrain_dX += mat_n1n2;
427 
428  Bmat_bend_w.left_multiply(mat_n1n2, eye);
429  dstrain_dX += mat_n1n2;
430  }
431 
432  // note: this assumes linear material laws
433  dstress_dX = material_mat * dstrain_dX;
434 
435  // copy to the 3D structure
436  dstress_dX_3D.row(0) = dstress_dX.row(0);
437  dstrain_dX_3D.row(0) = dstrain_dX.row(0);
438 
439  if (request_derivative)
440  data->set_derivatives(dstress_dX_3D, dstrain_dX_3D);
441 
442 
443  if (p) {
444  // sensitivity of the response, s, is
445  // ds/dp = partial s/partial p +
446  // partial s/partial X dX/dp
447  // the first part of the sensitivity is obtained from
448  //
449  // the first term includes direct sensitivity of the stress
450  // with respect to the parameter, while holding the solution
451  // constant. This should include influence of shape changes,
452  // if the parameter is shape-dependent.
453  // TODO: include shape sensitivity.
454  // presently, only material parameter is included
455 
456  dstrain_dp = RealVectorX::Zero(n1);
457 
458  // if thermal load was specified, then set the thermal strain
459  // component of the total strain
460  if (thermal_load) {
461  temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp);
462  ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t);
463  alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha);
464  dstrain_dp(0) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t);
465  }
466 
467  // include the dependence of strain on the thickness
468  if (if_bending) {
469 
470  // add to this the bending strain
471  hy.derivative(*p, xyz[qp_loc_index], _time, y);
472  hz.derivative(*p, xyz[qp_loc_index], _time, z);
473  hy_off.derivative(*p, xyz[qp_loc_index], _time, y_off);
474  hz_off.derivative(*p, xyz[qp_loc_index], _time, z_off);
475 
476  bend->initialize_bending_strain_operator_for_yz(*fe,
477  qp_loc_index,
478  qp_loc[qp](1) * y/2.+y_off,
479  qp_loc[qp](2) * z/2.+z_off,
480  Bmat_bend_v,
481  Bmat_bend_w);
482  Bmat_bend_v.vector_mult(strain_bend, _local_sol);
483  dstrain_dp += strain_bend;
484 
485  Bmat_bend_w.vector_mult(strain_bend, _local_sol);
486  dstrain_dp += strain_bend;
487  }
488 
489 
490  // now use this to calculate the stress sensitivity.
491  dstress_dp = material_mat * dstrain_dp;
492 
493  // get the material matrix sensitivity
494  mat_stiff.derivative(*p, xyz[qp_loc_index], _time, material_mat);
495 
496  // partial sensitivity of strain is zero unless it is a
497  // shape parameter.
498  // TODO: shape sensitivity of strain operator
499 
500  // now use this to calculate the stress sensitivity.
501  dstress_dp += material_mat * strain;
502 
503 
504  //
505  // use the derivative data to evaluate the second term in the
506  // sensitivity
507  //
508 
509  dstress_dp += dstress_dX * _local_sol_sens;
510  dstrain_dp += dstrain_dX * _local_sol_sens;
511 
512  // copy the 3D object
513  stress_3D(0) = dstress_dp(0);
514  strain_3D(0) = dstrain_dp(0);
515 
516  // tell the data object about the sensitivity values
517  data->set_sensitivity(*p,
518  stress_3D,
519  strain_3D);
520  }
521  }
522  }
523 
524  // make sure that the number of data points for this element is
525  // the same as the number of requested points
526  libmesh_assert(qp_loc.size() ==
527  stress_output.n_stress_strain_data_for_elem(_elem));
528 
529  // if either derivative or sensitivity was requested, it was provided
530  // by this routine
531  return request_derivative || p;
532 }
533 
534 
535 
536 
537 
538 bool
539 MAST::StructuralElement1D::internal_residual (bool request_jacobian,
540  RealVectorX& f,
541  RealMatrixX& jac)
542 {
543  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
544  false,
545  _property.extra_quadrature_order(_elem)));
546 
547  const std::vector<Real>& JxW = fe->get_JxW();
548  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
549  const unsigned int
550  n_phi = (unsigned int)fe->get_phi().size(),
551  n1 = this->n_direct_strain_components(),
552  n2 = 6*n_phi,
553  n3 = this->n_von_karman_strain_components();
554 
555  RealMatrixX
556  material_A_mat,
557  material_B_mat,
558  material_D_mat,
559  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
560  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
561  mat3,
562  mat4_n3n2 = RealMatrixX::Zero(n3,2),
563  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
564  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
565  stress = RealMatrixX::Zero(2,2),
566  stress_l = RealMatrixX::Zero(2,2),
567  local_jac = RealMatrixX::Zero(n2,n2);
568 
569  RealVectorX
570  vec1_n1 = RealVectorX::Zero(n1),
571  vec2_n1 = RealVectorX::Zero(n1),
572  vec3_n2 = RealVectorX::Zero(n2),
573  vec4_n3 = RealVectorX::Zero(n3),
574  vec5_n3 = RealVectorX::Zero(n3),
575  local_f = RealVectorX::Zero(n2);
576 
577  local_f.setZero();
578  local_jac.setZero();
579 
580  MAST::FEMOperatorMatrix
581  Bmat_mem,
582  Bmat_bend_v,
583  Bmat_bend_w,
584  Bmat_v_vk,
585  Bmat_w_vk;
586 
587  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
588  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
589  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
590  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
591  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
592 
593  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
594 
595  MAST::BendingOperatorType
596  bending_model = _property.bending_model(_elem);
597 
598  std::unique_ptr<MAST::BendingOperator1D> bend;
599 
600  if (bending_model != MAST::NO_BENDING)
601  bend.reset(MAST::build_bending_operator_1D(bending_model,
602  *this,
603  fe->get_qpoints()).release());
604 
605 
606  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
607  mat_stiff_A = _property.stiffness_A_matrix(*this),
608  mat_stiff_B = _property.stiffness_B_matrix(*this),
609  mat_stiff_D = _property.stiffness_D_matrix(*this);
610 
611  for (unsigned int qp=0; qp<JxW.size(); qp++) {
612 
613  // get the material matrix
614  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
615 
616  if (bend.get()) {
617  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
618  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
619  }
620 
621  // now calculte the quantity for these matrices
622  _internal_residual_operation(if_vk, n2, qp, *fe, JxW,
623  request_jacobian,
624  local_f, local_jac,
625  bend.get(),
626  Bmat_mem, Bmat_bend_v, Bmat_bend_w,
627  Bmat_v_vk, Bmat_w_vk,
628  stress, stress_l, vk_dvdxi_mat, vk_dwdxi_mat,
629  material_A_mat,
630  material_B_mat, material_D_mat, vec1_n1,
631  vec2_n1, vec3_n2, vec4_n3,
632  vec5_n3, mat1_n1n2, mat2_n2n2,
633  mat3, mat4_n3n2);
634 
635  }
636 
637 
638  // now calculate the transverse shear contribution if appropriate for the
639  // element
640  if (bend.get() && bend->include_transverse_shear_energy())
641  bend->calculate_transverse_shear_residual(request_jacobian,
642  local_f,
643  local_jac);
644 
645 
646  // now transform to the global coorodinate system
647  transform_vector_to_global_system(local_f, vec3_n2);
648  f += vec3_n2;
649 
650  if (request_jacobian) {
651  transform_matrix_to_global_system(local_jac, mat2_n2n2);
652  jac += mat2_n2n2;
653  }
654 
655  return request_jacobian;
656 }
657 
658 
659 
660 
661 
662 bool
663 MAST::StructuralElement1D::internal_residual_sensitivity (const MAST::FunctionBase& p,
664  bool request_jacobian,
665  RealVectorX& f,
666  RealMatrixX& jac)
667 {
668  // this should be true if the function is called
669  libmesh_assert(!p.is_shape_parameter()); // this is not implemented for now
670 
671 
672  // check if the material property or the provided exterior
673  // values, like temperature, are functions of the sensitivity parameter
674  bool calculate = false;
675  calculate = calculate || _property.depends_on(p);
676 
677  // nothing to be calculated if the element does not depend on the
678  // sensitivity parameter.
679  if (!calculate)
680  return false;
681 
682  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
683  false,
684  _property.extra_quadrature_order(_elem)));
685 
686  const std::vector<Real>& JxW = fe->get_JxW();
687  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
688  const unsigned int
689  n_phi = (unsigned int)fe->get_phi().size(),
690  n1 = this->n_direct_strain_components(),
691  n2 = 6*n_phi,
692  n3 = this->n_von_karman_strain_components();
693 
694  RealMatrixX
695  material_A_mat,
696  material_B_mat,
697  material_D_mat,
698  material_trans_shear_mat,
699  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
700  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
701  mat3,
702  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
703  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
704  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
705  stress = RealMatrixX::Zero(2,2),
706  stress_l = RealMatrixX::Zero(2,2),
707  local_jac = RealMatrixX::Zero(n2,n2);
708  RealVectorX
709  vec1_n1 = RealVectorX::Zero(n1),
710  vec2_n1 = RealVectorX::Zero(n1),
711  vec3_n2 = RealVectorX::Zero(n2),
712  vec4_n3 = RealVectorX::Zero(n3),
713  vec5_n3 = RealVectorX::Zero(n3),
714  local_f = RealVectorX::Zero(n2);
715 
716  local_f.setZero();
717  local_jac.setZero();
718 
719  MAST::FEMOperatorMatrix
720  Bmat_mem,
721  Bmat_bend_v,
722  Bmat_bend_w,
723  Bmat_v_vk,
724  Bmat_w_vk;
725 
726  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
727  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
728  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
729  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
730  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
731 
732  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
733 
734  MAST::BendingOperatorType
735  bending_model = _property.bending_model(_elem);
736 
737  std::unique_ptr<MAST::BendingOperator1D> bend;
738 
739  if (bending_model != MAST::NO_BENDING)
740  bend.reset(MAST::build_bending_operator_1D(bending_model,
741  *this,
742  fe->get_qpoints()).release());
743 
744  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
745  mat_stiff_A = _property.stiffness_A_matrix(*this),
746  mat_stiff_B = _property.stiffness_B_matrix(*this),
747  mat_stiff_D = _property.stiffness_D_matrix(*this);
748 
749  // first calculate the sensitivity due to the parameter
750  for (unsigned int qp=0; qp<JxW.size(); qp++) {
751 
752  // get the material matrix
753  mat_stiff_A->derivative(p, xyz[qp], _time, material_A_mat);
754 
755  if (bend.get()) {
756  mat_stiff_B->derivative(p, xyz[qp], _time, material_B_mat);
757  mat_stiff_D->derivative(p, xyz[qp], _time, material_D_mat);
758  }
759 
760  // now calculte the quantity for these matrices
761  _internal_residual_operation(if_vk, n2, qp, *fe, JxW,
762  request_jacobian,
763  local_f, local_jac,
764  bend.get(),
765  Bmat_mem, Bmat_bend_v, Bmat_bend_w,
766  Bmat_v_vk, Bmat_w_vk,
767  stress, stress_l, vk_dvdxi_mat, vk_dwdxi_mat,
768  material_A_mat,
769  material_B_mat, material_D_mat, vec1_n1,
770  vec2_n1, vec3_n2, vec4_n3,
771  vec5_n3, mat1_n1n2, mat2_n2n2,
772  mat3, mat4_n3n2);
773  }
774 
775  // now calculate the transverse shear contribution if appropriate for the
776  // element
777  if (bend.get() && bend->include_transverse_shear_energy())
778  bend->calculate_transverse_shear_residual_sensitivity(p,
779  request_jacobian,
780  local_f,
781  local_jac);
782 
783  // now transform to the global coorodinate system
784  transform_vector_to_global_system(local_f, vec3_n2);
785  f += vec3_n2;
786 
787  if (request_jacobian) {
788  transform_matrix_to_global_system(local_jac, mat2_n2n2);
789  jac += mat2_n2n2;
790  }
791 
792  return request_jacobian;
793 }
794 
795 
796 
797 
798 bool
799 MAST::StructuralElement1D::
800 internal_residual_jac_dot_state_sensitivity (RealMatrixX& jac) {
801 
802  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
803  false,
804  _property.extra_quadrature_order(_elem)));
805 
806  const std::vector<Real>& JxW = fe->get_JxW();
807  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
808  const unsigned int
809  n_phi = (unsigned int)fe->get_phi().size(),
810  n1 = this->n_direct_strain_components(),
811  n2 = 6*n_phi,
812  n3 = this->n_von_karman_strain_components();
813 
814  RealMatrixX
815  material_A_mat,
816  material_B_mat,
817  material_D_mat,
818  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
819  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
820  mat3,
821  vk_dvdxi_mat_sens = RealMatrixX::Zero(n1,n3),
822  vk_dwdxi_mat_sens = RealMatrixX::Zero(n1,n3),
823  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
824  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
825  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
826  stress = RealMatrixX::Zero(2,2),
827  local_jac = RealMatrixX::Zero(n2,n2);
828  RealVectorX
829  vec1_n1 = RealVectorX::Zero(n1),
830  vec2_n1 = RealVectorX::Zero(n1),
831  vec3_n2 = RealVectorX::Zero(n2),
832  vec4_n3 = RealVectorX::Zero(n3),
833  vec5_n3 = RealVectorX::Zero(n3);
834 
835  local_jac.setZero();
836 
837 
838  MAST::FEMOperatorMatrix
839  Bmat_mem,
840  Bmat_bend_v,
841  Bmat_bend_w,
842  Bmat_v_vk,
843  Bmat_w_vk;
844 
845  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
846  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
847  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
848  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
849  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
850 
851  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
852 
853  MAST::BendingOperatorType
854  bending_model = _property.bending_model(_elem);
855 
856  std::unique_ptr<MAST::BendingOperator1D> bend;
857 
858  if (bending_model != MAST::NO_BENDING)
859  bend.reset(MAST::build_bending_operator_1D(bending_model,
860  *this,
861  fe->get_qpoints()).release());
862 
863  // without the nonlinear strain, this matrix is zero.
864  if (!if_vk)
865  return false;
866 
867  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
868  mat_stiff_A = _property.stiffness_A_matrix(*this),
869  mat_stiff_B = _property.stiffness_B_matrix(*this),
870  mat_stiff_D = _property.stiffness_D_matrix(*this);
871 
872 
873  for (unsigned int qp=0; qp<JxW.size(); qp++) {
874 
875  // get the material matrix
876  (*mat_stiff_A)(xyz[qp], _time, material_A_mat);
877 
878  (*mat_stiff_B)(xyz[qp], _time, material_B_mat);
879  (*mat_stiff_D)(xyz[qp], _time, material_D_mat);
880 
881  // now calculte the quantity for these matrices
882  this->initialize_direct_strain_operator(qp, *fe, Bmat_mem);
883 
884  // first handle constant throught the thickness stresses: membrane and vonKarman
885  Bmat_mem.vector_mult(vec1_n1, _local_sol_sens);
886  vec2_n1 = material_A_mat * vec1_n1; // linear direct stress
887 
888  // copy the stress values to a matrix
889  stress(0,0) = vec2_n1(0); // sigma_xx
890 
891  // get the bending strain operator
892  vec2_n1.setZero(); // used to store vk strain, if applicable
893  if (bend.get()) {
894  bend->initialize_bending_strain_operator(*fe, qp,
895  Bmat_bend_v,
896  Bmat_bend_w);
897 
898  // evaluate the bending stress and add that to the stress vector
899  // for evaluation in the nonlinear stress term
900  Bmat_bend_v.vector_mult(vec2_n1, _local_sol_sens);
901  vec1_n1 = material_B_mat * vec2_n1;
902  stress(0,0) += vec1_n1(0);
903 
904  Bmat_bend_w.vector_mult(vec2_n1, _local_sol_sens);
905  vec1_n1 = material_B_mat * vec2_n1;
906  stress(0,0) += vec1_n1(0);
907 
908 
909  if (if_vk) { // get the vonKarman strain operator if needed
910 
911  this->initialize_von_karman_strain_operator(qp,
912  *fe,
913  vec2_n1,
914  vk_dvdxi_mat,
915  vk_dwdxi_mat,
916  Bmat_v_vk,
917  Bmat_w_vk);
918  this->initialize_von_karman_strain_operator_sensitivity(qp,
919  *fe,
920  vk_dvdxi_mat_sens,
921  vk_dwdxi_mat_sens);
922  // sensitivity of von Karman strain
923  vec2_n1.setZero();
924  vec2_n1(0) = (vk_dvdxi_mat(0,0)*vk_dvdxi_mat_sens(0,0) +
925  vk_dwdxi_mat(0,0)*vk_dwdxi_mat_sens(0,0));
926  vec1_n1 = material_A_mat * vec2_n1;
927  stress(0,0) += vec1_n1(0);
928  }
929  }
930 
931  // copy the stress to use here.
932  vec1_n1.setZero();
933 
934  // now calculate the matrix
935  // membrane - vk: v-displacement
936  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
937  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat_sens);
938  mat3 = material_A_mat * mat3;
939  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat3);
940  local_jac += JxW[qp] * mat2_n2n2;
941 
942  // membrane - vk: w-displacement
943  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
944  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
945  mat3 = material_A_mat * mat3;
946  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat3);
947  local_jac += JxW[qp] * mat2_n2n2;
948 
949  // vk - membrane: v-displacement
950  Bmat_mem.left_multiply(mat1_n1n2, material_A_mat);
951  mat3 = vk_dvdxi_mat_sens.transpose() * mat1_n1n2;
952  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
953  local_jac += JxW[qp] * mat2_n2n2;
954 
955  // vk - membrane: w-displacement
956  Bmat_mem.left_multiply(mat1_n1n2, material_A_mat);
957  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
958  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
959  local_jac += JxW[qp] * mat2_n2n2;
960 
961  // vk - vk: v-displacement
962  mat3 = RealMatrixX::Zero(2, n2);
963  Bmat_v_vk.left_multiply(mat3, stress);
964  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
965  local_jac += JxW[qp] * mat2_n2n2;
966 
967  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
968  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
969  mat3 = vk_dvdxi_mat_sens.transpose() * material_A_mat * mat3;
970  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
971  local_jac += JxW[qp] * mat2_n2n2;
972 
973  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
974  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat_sens);
975  mat3 = vk_dvdxi_mat.transpose() * material_A_mat * mat3;
976  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
977  local_jac += JxW[qp] * mat2_n2n2;
978 
979  // vk - vk: w-displacement
980  mat3 = RealMatrixX::Zero(2, n2);
981  Bmat_w_vk.left_multiply(mat3, stress);
982  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
983  local_jac += JxW[qp] * mat2_n2n2;
984 
985  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
986  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
987  mat3 = vk_dwdxi_mat_sens.transpose() * material_A_mat * mat3;
988  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
989  local_jac += JxW[qp] * mat2_n2n2;
990 
991  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
992  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
993  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
994  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
995  local_jac += JxW[qp] * mat2_n2n2;
996 
997  // coupling of v, w-displacements
998  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
999  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1000  mat3 = vk_dvdxi_mat.transpose() * material_A_mat * mat3;
1001  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1002  local_jac += JxW[qp] * mat2_n2n2;
1003 
1004  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1005  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
1006  mat3 = vk_dvdxi_mat_sens.transpose() * material_A_mat * mat3;
1007  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1008  local_jac += JxW[qp] * mat2_n2n2;
1009 
1010  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1011  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat_sens);
1012  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1013  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1014  local_jac += JxW[qp] * mat2_n2n2;
1015 
1016  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1017  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
1018  mat3 = vk_dwdxi_mat_sens.transpose() * material_A_mat * mat3;
1019  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1020  local_jac += JxW[qp] * mat2_n2n2;
1021 
1022  // bending - vk: v-displacement
1023  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1024  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat_sens);
1025  mat3 = material_B_mat.transpose() * mat3;
1026  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat3);
1027  local_jac += JxW[qp] * mat2_n2n2;
1028  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat3);
1029  local_jac += JxW[qp] * mat2_n2n2;
1030 
1031  // bending - vk: w-displacement
1032  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1033  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat_sens);
1034  mat3 = material_B_mat.transpose() * mat3;
1035  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat3);
1036  local_jac += JxW[qp] * mat2_n2n2;
1037  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat3);
1038  local_jac += JxW[qp] * mat2_n2n2;
1039 
1040  // vk - bending: v-displacement
1041  Bmat_bend_v.left_multiply(mat1_n1n2, material_B_mat);
1042  mat3 = vk_dvdxi_mat_sens.transpose() * mat1_n1n2;
1043  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1044  local_jac += JxW[qp] * mat2_n2n2;
1045 
1046  Bmat_bend_v.left_multiply(mat1_n1n2, material_B_mat);
1047  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1048  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1049  local_jac += JxW[qp] * mat2_n2n2;
1050 
1051  // vk - bending: w-displacement
1052  Bmat_bend_w.left_multiply(mat1_n1n2, material_B_mat);
1053  mat3 = vk_dvdxi_mat_sens.transpose() * mat1_n1n2;
1054  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1055  local_jac += JxW[qp] * mat2_n2n2;
1056 
1057  Bmat_bend_w.left_multiply(mat1_n1n2, material_B_mat);
1058  mat3 = vk_dwdxi_mat_sens.transpose() * mat1_n1n2;
1059  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1060  local_jac += JxW[qp] * mat2_n2n2;
1061  }
1062 
1063  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1064  jac += mat2_n2n2;
1065 
1066  return true;
1067 }
1068 
1069 
1070 
1071 
1072 void
1073 MAST::StructuralElement1D::
1074 _internal_residual_operation(bool if_vk,
1075  const unsigned int n2,
1076  const unsigned int qp,
1077  const MAST::FEBase& fe,
1078  const std::vector<Real>& JxW,
1079  bool request_jacobian,
1080  RealVectorX& local_f,
1081  RealMatrixX& local_jac,
1082  MAST::BendingOperator1D* bend,
1083  MAST::FEMOperatorMatrix& Bmat_mem,
1084  MAST::FEMOperatorMatrix& Bmat_bend_v,
1085  MAST::FEMOperatorMatrix& Bmat_bend_w,
1086  MAST::FEMOperatorMatrix& Bmat_v_vk,
1087  MAST::FEMOperatorMatrix& Bmat_w_vk,
1088  RealMatrixX& stress,
1089  RealMatrixX& stress_l,
1090  RealMatrixX& vk_dvdxi_mat,
1091  RealMatrixX& vk_dwdxi_mat,
1092  RealMatrixX& material_A_mat,
1093  RealMatrixX& material_B_mat,
1094  RealMatrixX& material_D_mat,
1095  RealVectorX& vec1_n1,
1096  RealVectorX& vec2_n1,
1097  RealVectorX& vec3_n2,
1098  RealVectorX& vec4_2,
1099  RealVectorX& vec5_2,
1100  RealMatrixX& mat1_n1n2,
1101  RealMatrixX& mat2_n2n2,
1102  RealMatrixX& mat3,
1103  RealMatrixX& mat4_2n2)
1104 {
1105  this->initialize_direct_strain_operator(qp, fe, Bmat_mem);
1106 
1107  // first handle constant throught the thickness stresses: membrane and vonKarman
1108  Bmat_mem.vector_mult(vec1_n1, _local_sol);
1109  vec2_n1 = material_A_mat * vec1_n1; // linear direct stress
1110 
1111  // copy the stress values to a matrix
1112  stress_l(0,0) = vec2_n1(0); // sigma_xx
1113  stress(0,0) = vec2_n1(0);
1114  stress(1,1) = vec1_n1(0); // temporary storage of the axial strain
1115  stress(0,1) = vec2_n1(1); // temporary storage of the torsion force
1116 
1117  // get the bending strain operator
1118  vec2_n1.setZero(); // used to store vk strain, if applicable
1119  if (bend) {
1120  bend->initialize_bending_strain_operator(fe, qp,
1121  Bmat_bend_v,
1122  Bmat_bend_w);
1123 
1124  // evaluate the bending stress and add that to the stress vector
1125  // for evaluation in the nonlinear stress term
1126  Bmat_bend_v.vector_mult(vec2_n1, _local_sol);
1127  vec1_n1 = material_B_mat * vec2_n1;
1128  stress_l(0,0) += vec1_n1(0);
1129  stress(0,0) += vec1_n1(0);
1130 
1131  Bmat_bend_w.vector_mult(vec2_n1, _local_sol);
1132  vec1_n1 = material_B_mat * vec2_n1;
1133  stress_l(0,0) += vec1_n1(0);
1134  stress(0,0) += vec1_n1(0);
1135 
1136  if (if_vk) { // get the vonKarman strain operator if needed
1137 
1138  this->initialize_von_karman_strain_operator(qp,
1139  fe,
1140  vec2_n1, // epsilon_vk
1141  vk_dvdxi_mat,
1142  vk_dwdxi_mat,
1143  Bmat_v_vk,
1144  Bmat_w_vk);
1145  vec1_n1 = material_A_mat * vec2_n1;
1146  stress(0,0) += vec1_n1(0); // total strain that multiplies with the membrane strain
1147  stress(1,1) += vec2_n1(0); // add the two strains to get the direct strain
1148  }
1149  }
1150 
1151  // copy the stress to use here.
1152  vec1_n1.setZero();
1153  vec1_n1(0) = stress(0,0);
1154  vec1_n1(1) = stress(0,1); // use the torsion strain from the temporary location
1155  stress(0, 1) = 0.; // zero out the temporary value storing the torsion strain
1156 
1157  // now the internal force vector
1158  // this includes the membrane strain operator with all A and B material operators
1159  Bmat_mem.vector_mult_transpose(vec3_n2, vec1_n1);
1160  local_f += JxW[qp] * vec3_n2;
1161 
1162  if (bend) {
1163  if (if_vk) {
1164  // von Karman strain: direct stress
1165  vec4_2 = vk_dvdxi_mat.transpose() * vec1_n1;
1166  Bmat_v_vk.vector_mult_transpose(vec3_n2, vec4_2);
1167  local_f += JxW[qp] * vec3_n2;
1168 
1169  // von Karman strain: direct stress
1170  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
1171  Bmat_w_vk.vector_mult_transpose(vec3_n2, vec4_2);
1172  local_f += JxW[qp] * vec3_n2;
1173  }
1174 
1175  // use the direct strain from the temprary storage
1176  vec2_n1(0) = stress(1,1);
1177  stress(1,1) = 0.;
1178  // now coupling with the bending strain
1179  // B_bend^T [B] B_mem
1180  vec1_n1 = material_B_mat.transpose() * vec2_n1;
1181  Bmat_bend_v.vector_mult_transpose(vec3_n2, vec1_n1);
1182  local_f += JxW[qp] * vec3_n2;
1183  Bmat_bend_w.vector_mult_transpose(vec3_n2, vec1_n1);
1184  local_f += JxW[qp] * vec3_n2;
1185 
1186  // now bending stress
1187  Bmat_bend_v.vector_mult(vec2_n1, _local_sol);
1188  vec1_n1 = material_D_mat * vec2_n1;
1189  Bmat_bend_v.vector_mult_transpose(vec3_n2, vec1_n1);
1190  local_f += JxW[qp] * vec3_n2;
1191 
1192  Bmat_bend_w.vector_mult(vec2_n1, _local_sol);
1193  vec1_n1 = material_D_mat * vec2_n1;
1194  Bmat_bend_w.vector_mult_transpose(vec3_n2, vec1_n1);
1195  local_f += JxW[qp] * vec3_n2;
1196  }
1197 
1198  if (request_jacobian) {
1199  // membrane - membrane
1200  Bmat_mem.left_multiply(mat1_n1n2, material_A_mat);
1201  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1202  local_jac += JxW[qp] * mat2_n2n2;
1203 
1204  if (bend) {
1205  if (if_vk) {
1206  // membrane - vk: v-displacement
1207  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1208  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
1209  mat3 = material_A_mat * mat3;
1210  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat3);
1211  local_jac += JxW[qp] * mat2_n2n2;
1212 
1213  // membrane - vk: w-displacement
1214  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1215  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
1216  mat3 = material_A_mat * mat3;
1217  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat3);
1218  local_jac += JxW[qp] * mat2_n2n2;
1219 
1220  // vk - membrane: v-displacement
1221  Bmat_mem.left_multiply(mat1_n1n2, material_A_mat);
1222  mat3 = vk_dvdxi_mat.transpose() * mat1_n1n2;
1223  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1224  local_jac += JxW[qp] * mat2_n2n2;
1225 
1226  // vk - membrane: w-displacement
1227  Bmat_mem.left_multiply(mat1_n1n2, material_A_mat);
1228  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1229  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1230  local_jac += JxW[qp] * mat2_n2n2;
1231 
1232  // if only the first order term of the Jacobian is needed, for
1233  // example for linearized buckling analysis, then the linear
1234  // stress combined with the variation of the von Karman strain
1235  // is included. Otherwise, all terms are included
1236  /*if (if_ignore_ho_jac) {
1237  // vk - vk: v-displacement: first order term
1238  mat3 = RealMatrixX::Zero(2, n2);
1239  Bmat_v_vk.left_multiply(mat3, stress_l);
1240  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1241  local_jac += JxW[qp] * mat2_n2n2;
1242 
1243  // vk - vk: v-displacement: first order term
1244  mat3 = RealMatrixX::Zero(2, n2);
1245  Bmat_w_vk.left_multiply(mat3, stress_l);
1246  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1247  local_jac += JxW[qp] * mat2_n2n2;
1248  }
1249  else*/ {
1250  // vk - vk: v-displacement
1251  mat3 = RealMatrixX::Zero(2, n2);
1252  Bmat_v_vk.left_multiply(mat3, stress);
1253  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1254  local_jac += JxW[qp] * mat2_n2n2;
1255 
1256  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1257  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
1258  mat3 = vk_dvdxi_mat.transpose() * material_A_mat * mat3;
1259  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1260  local_jac += JxW[qp] * mat2_n2n2;
1261 
1262  // vk - vk: w-displacement
1263  mat3 = RealMatrixX::Zero(2, n2);
1264  Bmat_w_vk.left_multiply(mat3, stress);
1265  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1266  local_jac += JxW[qp] * mat2_n2n2;
1267 
1268  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1269  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
1270  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1271  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1272  local_jac += JxW[qp] * mat2_n2n2;
1273 
1274  // coupling of v, w-displacements
1275  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1276  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
1277  mat3 = vk_dvdxi_mat.transpose() * material_A_mat * mat3;
1278  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1279  local_jac += JxW[qp] * mat2_n2n2;
1280 
1281  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1282  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
1283  mat3 = vk_dwdxi_mat.transpose() * material_A_mat * mat3;
1284  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1285  local_jac += JxW[qp] * mat2_n2n2;
1286 
1287  }
1288 
1289  // bending - vk: v-displacement
1290  mat3 = RealMatrixX::Zero(vk_dvdxi_mat.rows(), n2);
1291  Bmat_v_vk.left_multiply(mat3, vk_dvdxi_mat);
1292  mat3 = material_B_mat.transpose() * mat3;
1293  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat3);
1294  local_jac += JxW[qp] * mat2_n2n2;
1295  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat3);
1296  local_jac += JxW[qp] * mat2_n2n2;
1297 
1298  // bending - vk: w-displacement
1299  mat3 = RealMatrixX::Zero(vk_dwdxi_mat.rows(), n2);
1300  Bmat_w_vk.left_multiply(mat3, vk_dwdxi_mat);
1301  mat3 = material_B_mat.transpose() * mat3;
1302  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat3);
1303  local_jac += JxW[qp] * mat2_n2n2;
1304  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat3);
1305  local_jac += JxW[qp] * mat2_n2n2;
1306 
1307  // vk - bending: v-displacement
1308  Bmat_bend_v.left_multiply(mat1_n1n2, material_B_mat);
1309  mat3 = vk_dvdxi_mat.transpose() * mat1_n1n2;
1310  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1311  local_jac += JxW[qp] * mat2_n2n2;
1312 
1313  Bmat_bend_v.left_multiply(mat1_n1n2, material_B_mat);
1314  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1315  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1316  local_jac += JxW[qp] * mat2_n2n2;
1317 
1318  // vk - bending: w-displacement
1319  Bmat_bend_w.left_multiply(mat1_n1n2, material_B_mat);
1320  mat3 = vk_dvdxi_mat.transpose() * mat1_n1n2;
1321  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1322  local_jac += JxW[qp] * mat2_n2n2;
1323 
1324  Bmat_bend_w.left_multiply(mat1_n1n2, material_B_mat);
1325  mat3 = vk_dwdxi_mat.transpose() * mat1_n1n2;
1326  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1327  local_jac += JxW[qp] * mat2_n2n2;
1328  }
1329 
1330  // bending - membrane
1331  mat3 = material_B_mat.transpose();
1332  Bmat_mem.left_multiply(mat1_n1n2, mat3);
1333  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1334  local_jac += JxW[qp] * mat2_n2n2;
1335  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1336  local_jac += JxW[qp] * mat2_n2n2;
1337 
1338  // membrane - bending
1339  Bmat_bend_v.left_multiply(mat1_n1n2, material_B_mat);
1340  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1341  local_jac += JxW[qp] * mat2_n2n2;
1342 
1343  Bmat_bend_w.left_multiply(mat1_n1n2, material_B_mat);
1344  Bmat_mem.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1345  local_jac += JxW[qp] * mat2_n2n2;
1346 
1347  // bending - bending
1348  Bmat_bend_v.left_multiply(mat1_n1n2, material_D_mat);
1349  Bmat_bend_v.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1350  local_jac += JxW[qp] * mat2_n2n2;
1351 
1352  Bmat_bend_w.left_multiply(mat1_n1n2, material_D_mat);
1353  Bmat_bend_w.right_multiply_transpose(mat2_n2n2, mat1_n1n2);
1354  local_jac += JxW[qp] * mat2_n2n2;
1355  }
1356  }
1357 }
1358 
1359 
1360 
1361 void
1362 MAST::StructuralElement1D::_convert_prestress_A_mat_to_vector(const RealMatrixX& mat,
1363  RealVectorX& vec) const {
1364 
1365  libmesh_assert_equal_to(mat.rows(), 2);
1366  libmesh_assert_equal_to(mat.cols(), 2);
1367  vec = RealVectorX::Zero(2);
1368  vec(0) = mat(0,0);
1369 }
1370 
1371 
1372 void
1373 MAST::StructuralElement1D::_convert_prestress_B_mat_to_vector(const RealMatrixX& mat,
1374  RealVectorX& vec) const {
1375 
1376  libmesh_assert_equal_to(mat.rows(), 2);
1377  libmesh_assert_equal_to(mat.cols(), 2);
1378  vec = RealVectorX::Zero(2);
1379  vec(0) = mat(0,0);
1380  vec(1) = mat(0,1);
1381 }
1382 
1383 
1384 
1385 bool
1386 MAST::StructuralElement1D::prestress_residual (bool request_jacobian,
1387  RealVectorX& f,
1388  RealMatrixX& jac)
1389 {
1390  if (!_property.if_prestressed())
1391  return false;
1392 
1393  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1394  false,
1395  _property.extra_quadrature_order(_elem)));
1396 
1397  const std::vector<Real>& JxW = fe->get_JxW();
1398  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1399  const unsigned int
1400  n_phi = (unsigned int)fe->get_phi().size(),
1401  n1 = this->n_direct_strain_components(),
1402  n2 =6*n_phi,
1403  n3 = this->n_von_karman_strain_components();
1404 
1405  RealMatrixX
1406  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
1407  mat3,
1408  vk_dvdxi_mat = RealMatrixX::Zero(n1, n3),
1409  vk_dwdxi_mat = RealMatrixX::Zero(n1, n3),
1410  local_jac = RealMatrixX::Zero(n2, n2),
1411  prestress_mat_A,
1412  prestress_mat_B;
1413  RealVectorX
1414  vec2_n1 = RealVectorX::Zero(n1),
1415  vec3_n2 = RealVectorX::Zero(n2),
1416  vec4_n3 = RealVectorX::Zero(n3),
1417  vec5_n3 = RealVectorX::Zero(n3),
1418  local_f = RealVectorX::Zero(n2),
1419  prestress_vec_A,
1420  prestress_vec_B;
1421 
1422  local_f.setZero();
1423  local_jac.setZero();
1424 
1425  MAST::FEMOperatorMatrix
1426  Bmat_mem,
1427  Bmat_bend_v,
1428  Bmat_bend_w,
1429  Bmat_v_vk,
1430  Bmat_w_vk;
1431 
1432  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1433  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
1434  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
1435  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
1436  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1437 
1438  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1439 
1440  MAST::BendingOperatorType
1441  bending_model = _property.bending_model(_elem);
1442 
1443  std::unique_ptr<MAST::BendingOperator1D> bend;
1444 
1445  if (bending_model != MAST::NO_BENDING)
1446  bend.reset(MAST::build_bending_operator_1D(bending_model,
1447  *this,
1448  fe->get_qpoints()).release());
1449 
1450  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
1451  prestress_A = _property.prestress_A_matrix(*this),
1452  prestress_B = _property.prestress_B_matrix(*this);
1453 
1454  // now calculate the quantity
1455  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1456 
1457  (*prestress_A)(xyz[qp], _time, prestress_mat_A);
1458  (*prestress_B)(xyz[qp], _time, prestress_mat_B);
1459  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
1460  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
1461 
1462  this->initialize_direct_strain_operator(qp, *fe, Bmat_mem);
1463 
1464  // get the bending strain operator if needed
1465  vec2_n1.setZero(); // used to store vk strain, if applicable
1466  if (bend.get()) {
1467  bend->initialize_bending_strain_operator(*fe, qp,
1468  Bmat_bend_v,
1469  Bmat_bend_w);
1470 
1471  if (if_vk) // get the vonKarman strain operator if needed
1472  this->initialize_von_karman_strain_operator(qp,
1473  *fe,
1474  vec2_n1,
1475  vk_dvdxi_mat,
1476  vk_dwdxi_mat,
1477  Bmat_v_vk,
1478  Bmat_w_vk);
1479  }
1480 
1481  // first handle constant throught the thickness stresses: membrane and vonKarman
1482  // multiply this with the constant through the thickness strain
1483  // membrane strain
1484  Bmat_mem.vector_mult_transpose(vec3_n2, prestress_vec_A);
1485  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
1486 
1487  if (bend.get()) {
1488  if (if_vk) {
1489  // von Karman strain: v-displacement
1490  vec4_n3 = vk_dvdxi_mat.transpose() * prestress_vec_A;
1491  Bmat_v_vk.vector_mult_transpose(vec3_n2, vec4_n3);
1492  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
1493 
1494  // von Karman strain: w-displacement
1495  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
1496  Bmat_w_vk.vector_mult_transpose(vec3_n2, vec4_n3);
1497  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
1498  }
1499 
1500  // now coupling with the bending strain
1501  Bmat_bend_v.vector_mult_transpose(vec3_n2, prestress_vec_B);
1502  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
1503  Bmat_bend_w.vector_mult_transpose(vec3_n2, prestress_vec_B);
1504  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
1505  }
1506 
1507  if (request_jacobian) {
1508  if (bend.get() && if_vk) {
1509  // v-displacement
1510  mat3 = RealMatrixX::Zero(2, n2);
1511  Bmat_v_vk.left_multiply(mat3, prestress_mat_A);
1512  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1513  local_jac += JxW[qp] * mat2_n2n2;
1514 
1515  // w-displacement
1516  mat3 = RealMatrixX::Zero(2, n2);
1517  Bmat_w_vk.left_multiply(mat3, prestress_mat_A);
1518  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1519  local_jac += JxW[qp] * mat2_n2n2;
1520  }
1521  }
1522  }
1523 
1524  // now transform to the global coorodinate system
1525  transform_vector_to_global_system(local_f, vec3_n2);
1526  f += vec3_n2;
1527  if (request_jacobian && if_vk) {
1528  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1529  jac += mat2_n2n2;
1530  }
1531 
1532  // only the nonlinear strain returns a Jacobian for prestressing
1533  return (request_jacobian);
1534 }
1535 
1536 
1537 
1538 
1539 
1540 bool
1541 MAST::StructuralElement1D::prestress_residual_sensitivity (const MAST::FunctionBase& p,
1542  bool request_jacobian,
1543  RealVectorX& f,
1544  RealMatrixX& jac)
1545 {
1546  if (!_property.if_prestressed())
1547  return false;
1548 
1549  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1550  false,
1551  _property.extra_quadrature_order(_elem)));
1552 
1553  const std::vector<Real>& JxW = fe->get_JxW();
1554  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1555  const unsigned int
1556  n_phi = (unsigned int)fe->get_phi().size(),
1557  n1 = this->n_direct_strain_components(),
1558  n2 = 6*n_phi,
1559  n3 = this->n_von_karman_strain_components();
1560 
1561  RealMatrixX
1562  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1563  mat3,
1564  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1565  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
1566  local_jac = RealMatrixX::Zero(n2,n2),
1567  prestress_mat_A,
1568  prestress_mat_B;
1569  RealVectorX
1570  vec2_n1 = RealVectorX::Zero(n1),
1571  vec3_n2 = RealVectorX::Zero(n2),
1572  vec4_n3 = RealVectorX::Zero(n3),
1573  vec5_n3 = RealVectorX::Zero(n3),
1574  local_f = RealVectorX::Zero(n2),
1575  prestress_vec_A,
1576  prestress_vec_B;
1577 
1578  local_f.setZero();
1579  local_jac.setZero();
1580 
1581  MAST::FEMOperatorMatrix
1582  Bmat_mem,
1583  Bmat_bend_v,
1584  Bmat_bend_w,
1585  Bmat_v_vk,
1586  Bmat_w_vk;
1587 
1588  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1589  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
1590  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
1591  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
1592  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1593 
1594  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1595 
1596  MAST::BendingOperatorType
1597  bending_model = _property.bending_model(_elem);
1598 
1599  std::unique_ptr<MAST::BendingOperator1D> bend;
1600 
1601  if (bending_model != MAST::NO_BENDING)
1602  bend.reset(MAST::build_bending_operator_1D(bending_model,
1603  *this,
1604  fe->get_qpoints()).release());
1605 
1606  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
1607  prestress_A = _property.prestress_A_matrix(*this),
1608  prestress_B = _property.prestress_B_matrix(*this);
1609 
1610  // transform to the local coordinate system
1611  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1612 
1613  prestress_A->derivative(p, xyz[qp], _time, prestress_mat_A);
1614  prestress_B->derivative(p, xyz[qp], _time, prestress_mat_B);
1615  _convert_prestress_A_mat_to_vector(prestress_mat_A, prestress_vec_A);
1616  _convert_prestress_B_mat_to_vector(prestress_mat_B, prestress_vec_B);
1617 
1618  this->initialize_direct_strain_operator(qp, *fe, Bmat_mem);
1619 
1620  // get the bending strain operator if needed
1621  vec2_n1.setZero(); // used to store vk strain, if applicable
1622  if (bend.get()) {
1623  bend->initialize_bending_strain_operator(*fe, qp,
1624  Bmat_bend_v,
1625  Bmat_bend_w);
1626 
1627  if (if_vk) // get the vonKarman strain operator if needed
1628  this->initialize_von_karman_strain_operator(qp,
1629  *fe,
1630  vec2_n1,
1631  vk_dvdxi_mat,
1632  vk_dwdxi_mat,
1633  Bmat_v_vk,
1634  Bmat_w_vk);
1635  }
1636 
1637  // first handle constant throught the thickness stresses: membrane and vonKarman
1638  // multiply this with the constant through the thickness strain
1639  // membrane strain
1640  Bmat_mem.vector_mult_transpose(vec3_n2, prestress_vec_A);
1641  local_f += JxW[qp] * vec3_n2; // epsilon_mem * sigma_0
1642 
1643  if (bend.get()) {
1644  if (if_vk) {
1645  // von Karman strain: v-displacement
1646  vec4_n3 = vk_dvdxi_mat.transpose() * prestress_vec_A;
1647  Bmat_v_vk.vector_mult_transpose(vec3_n2, vec4_n3);
1648  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
1649 
1650  // von Karman strain: w-displacement
1651  vec4_n3 = vk_dwdxi_mat.transpose() * prestress_vec_A;
1652  Bmat_w_vk.vector_mult_transpose(vec3_n2, vec4_n3);
1653  local_f += JxW[qp] * vec3_n2; // epsilon_vk * sigma_0
1654  }
1655 
1656  // now coupling with the bending strain
1657  Bmat_bend_v.vector_mult_transpose(vec3_n2, prestress_vec_B);
1658  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
1659  Bmat_bend_w.vector_mult_transpose(vec3_n2, prestress_vec_B);
1660  local_f += JxW[qp] * vec3_n2; // epsilon_bend * sigma_0
1661  }
1662 
1663  if (request_jacobian) {
1664  if (bend.get() && if_vk) {
1665  // v-displacement
1666  mat3 = RealMatrixX::Zero(2, n2);
1667  Bmat_v_vk.left_multiply(mat3, prestress_mat_A);
1668  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1669  local_jac += JxW[qp] * mat2_n2n2;
1670 
1671  // w-displacement
1672  mat3 = RealMatrixX::Zero(2, n2);
1673  Bmat_w_vk.left_multiply(mat3, prestress_mat_A);
1674  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
1675  local_jac += JxW[qp] * mat2_n2n2;
1676  }
1677  }
1678  }
1679 
1680  // now transform to the global coorodinate system
1681  transform_vector_to_global_system(local_f, vec3_n2);
1682  f += vec3_n2;
1683  if (request_jacobian && if_vk) {
1684  transform_matrix_to_global_system(local_jac, mat2_n2n2);
1685  jac += mat2_n2n2;
1686  }
1687 
1688  // only the nonlinear strain returns a Jacobian for prestressing
1689  return (request_jacobian);
1690 }
1691 
1692 
1693 
1694 
1695 bool
1696 MAST::StructuralElement1D::
1697 surface_pressure_residual(bool request_jacobian,
1698  RealVectorX &f,
1699  RealMatrixX &jac,
1700  const unsigned int side,
1701  MAST::BoundaryConditionBase& bc) {
1702 
1703  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1704 
1705  // prepare the side finite element
1706  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false));
1707 
1708  const std::vector<Real> &JxW = fe->get_JxW();
1709  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1710  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1711  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
1712  const unsigned int
1713  n_phi = (unsigned int)phi.size(),
1714  n1 = 3,
1715  n2 = 6*n_phi;
1716 
1717 
1718  // get the function from this boundary condition
1719  const MAST::FieldFunction<Real>& p_func =
1720  bc.get<MAST::FieldFunction<Real> >("pressure");
1721 
1722  const MAST::FieldFunction<Real>& A_func =
1723  dynamic_cast<const MAST::ElementPropertyCard1D&>(_property).A();
1724 
1725  FEMOperatorMatrix Bmat;
1726  Real
1727  press = 0.,
1728  A_val = 0.;
1729 
1730  RealVectorX
1731  phi_vec = RealVectorX::Zero(n_phi),
1732  force = RealVectorX::Zero(2*n1),
1733  local_f = RealVectorX::Zero(n2),
1734  vec_n2 = RealVectorX::Zero(n2);
1735 
1736  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1737 
1738  // now set the shape function values
1739  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1740  phi_vec(i_nd) = phi[i_nd][qp];
1741 
1742  Bmat.reinit(2*n1, phi_vec);
1743 
1744  // get pressure cross-sectional area value
1745  p_func(qpoint[qp], _time, press);
1746  A_func(qpoint[qp], _time, A_val);
1747 
1748  // calculate force
1749  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1750  force(i_dim) = (press*A_val) * face_normals[qp](i_dim);
1751 
1752  Bmat.vector_mult_transpose(vec_n2, force);
1753 
1754  local_f += JxW[qp] * vec_n2;
1755  }
1756 
1757  // now transform to the global system and add
1758  if (_elem.dim() < 3) {
1759  transform_vector_to_global_system(local_f, vec_n2);
1760  f -= vec_n2;
1761  }
1762  else
1763  f -= local_f;
1764 
1765 
1766  return (request_jacobian);
1767 }
1768 
1769 
1770 
1771 
1772 
1773 bool
1774 MAST::StructuralElement1D::
1775 surface_pressure_residual_sensitivity(const MAST::FunctionBase& p,
1776  bool request_jacobian,
1777  RealVectorX &f,
1778  RealMatrixX &jac,
1779  const unsigned int side,
1780  MAST::BoundaryConditionBase& bc) {
1781 
1782  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1783 
1784  // prepare the side finite element
1785  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(side, false));
1786 
1787  const std::vector<Real> &JxW = fe->get_JxW();
1788  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1789  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1790  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
1791  const unsigned int
1792  n_phi = (unsigned int)phi.size(),
1793  n1 = 3,
1794  n2 = 6*n_phi;
1795 
1796 
1797  // get the function from this boundary condition
1798  const MAST::FieldFunction<Real>& p_func =
1799  bc.get<MAST::FieldFunction<Real> >("pressure");
1800  const MAST::FieldFunction<Real>& A_func =
1801  dynamic_cast<const MAST::ElementPropertyCard1D&>(_property).A();
1802 
1803 
1804  FEMOperatorMatrix Bmat;
1805  Real
1806  press = 0.,
1807  dpress = 0.,
1808  A_val = 0.,
1809  dA_val = 0.;
1810 
1811  RealVectorX
1812  phi_vec = RealVectorX::Zero(n_phi),
1813  force = RealVectorX::Zero(2*n1),
1814  local_f = RealVectorX::Zero(n2),
1815  vec_n2 = RealVectorX::Zero(n2);
1816 
1817  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1818 
1819  // now set the shape function values
1820  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1821  phi_vec(i_nd) = phi[i_nd][qp];
1822 
1823  Bmat.reinit(2*n1, phi_vec);
1824 
1825  // get pressure and area values and their sensitivities
1826  p_func(qpoint[qp], _time, press);
1827  p_func.derivative(p, qpoint[qp], _time, dpress);
1828  A_func(qpoint[qp], _time, A_val);
1829  A_func.derivative(p, qpoint[qp], _time, dA_val);
1830 
1831  // calculate force
1832  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1833  force(i_dim) = (press*dA_val + dpress*A_val) * face_normals[qp](i_dim);
1834 
1835  Bmat.vector_mult_transpose(vec_n2, force);
1836 
1837  local_f += JxW[qp] * vec_n2;
1838  }
1839 
1840  // now transform to the global system and add
1841  if (_elem.dim() < 3) {
1842  transform_vector_to_global_system(local_f, vec_n2);
1843  f -= vec_n2;
1844  }
1845  else
1846  f -= local_f;
1847 
1848 
1849  return (request_jacobian);
1850 }
1851 
1852 
1853 
1854 
1855 
1856 
1857 bool
1858 MAST::StructuralElement1D::thermal_residual (bool request_jacobian,
1859  RealVectorX& f,
1860  RealMatrixX& jac,
1861  MAST::BoundaryConditionBase& bc)
1862 {
1863  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
1864  false,
1865  _property.extra_quadrature_order(_elem)));
1866 
1867  const std::vector<Real>& JxW = fe->get_JxW();
1868  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1869  const unsigned int
1870  n_phi = (unsigned int)fe->get_phi().size(),
1871  n1 = this->n_direct_strain_components(),
1872  n2 = 6*n_phi,
1873  n3 = this->n_von_karman_strain_components();
1874 
1875  RealMatrixX
1876  material_exp_A_mat,
1877  material_exp_B_mat,
1878  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
1879  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
1880  mat3,
1881  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
1882  vk_dvdxi_mat = RealMatrixX::Zero(n1,n3),
1883  vk_dwdxi_mat = RealMatrixX::Zero(n1,n3),
1884  stress = RealMatrixX::Zero(2,2),
1885  local_jac = RealMatrixX::Zero(n2, n2);
1886  RealVectorX
1887  vec1_n1 = RealVectorX::Zero(n1),
1888  vec2_n1 = RealVectorX::Zero(n1),
1889  vec3_n2 = RealVectorX::Zero(n2),
1890  vec4_2 = RealVectorX::Zero(2),
1891  vec5_n3 = RealVectorX::Zero(n3),
1892  local_f = RealVectorX::Zero(n2),
1893  delta_t = RealVectorX::Zero(1);
1894 
1895  local_f.setZero();
1896  local_jac.setZero();
1897 
1898  MAST::FEMOperatorMatrix
1899  Bmat_mem,
1900  Bmat_bend_v,
1901  Bmat_bend_w,
1902  Bmat_v_vk,
1903  Bmat_w_vk;
1904 
1905  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
1906  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
1907  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
1908  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
1909  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
1910 
1911  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
1912 
1913  MAST::BendingOperatorType
1914  bending_model = _property.bending_model(_elem);
1915 
1916  std::unique_ptr<MAST::BendingOperator1D> bend;
1917 
1918  if (bending_model != MAST::NO_BENDING)
1919  bend.reset(MAST::build_bending_operator_1D(bending_model,
1920  *this,
1921  fe->get_qpoints()).release());
1922 
1923  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
1924  expansion_A = _property.thermal_expansion_A_matrix(*this),
1925  expansion_B = _property.thermal_expansion_B_matrix(*this);
1926 
1927  // temperature function
1928  const MAST::FieldFunction<Real>
1929  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
1930  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
1931 
1932  Real
1933  t = 0.,
1934  t0 = 0.,
1935  scaling = 1.;
1936 
1937  if (bc.contains("thermal_jacobian_scaling"))
1938  bc.get<MAST::FieldFunction<Real>>("thermal_jacobian_scaling")(scaling);
1939 
1940  for (unsigned int qp=0; qp<JxW.size(); qp++) {
1941 
1942  // get the material property
1943  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
1944  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
1945 
1946  // get the temperature function
1947  temp_func (xyz[qp], _time, t);
1948  ref_temp_func(xyz[qp], _time, t0);
1949  delta_t(0) = t-t0;
1950 
1951  vec1_n1 = material_exp_A_mat * delta_t; // [C]{alpha (T - T0)} (with membrane strain)
1952  stress(0,0) = vec1_n1(0); // sigma_xx
1953  vec2_n1 = material_exp_B_mat * delta_t; // [C]{alpha (T - T0)} (with bending strain)
1954 
1955  this->initialize_direct_strain_operator(qp, *fe, Bmat_mem);
1956 
1957  // membrane strain
1958  Bmat_mem.vector_mult_transpose(vec3_n2, vec1_n1);
1959  local_f += JxW[qp] * vec3_n2;
1960 
1961  if (bend.get()) {
1962  // bending strain
1963  bend->initialize_bending_strain_operator(*fe, qp,
1964  Bmat_bend_v,
1965  Bmat_bend_w);
1966  Bmat_bend_v.vector_mult_transpose(vec3_n2, vec2_n1);
1967  local_f += JxW[qp] * vec3_n2;
1968  Bmat_bend_w.vector_mult_transpose(vec3_n2, vec2_n1);
1969  local_f += JxW[qp] * vec3_n2;
1970 
1971  // von Karman strain
1972  if (if_vk) {
1973  // get the vonKarman strain operator if needed
1974  this->initialize_von_karman_strain_operator(qp,
1975  *fe,
1976  vec2_n1, // epsilon_vk
1977  vk_dvdxi_mat,
1978  vk_dwdxi_mat,
1979  Bmat_v_vk,
1980  Bmat_w_vk);
1981  // von Karman strain: v-displacement
1982  vec4_2 = vk_dvdxi_mat.transpose() * vec1_n1;
1983  Bmat_v_vk.vector_mult_transpose(vec3_n2, vec4_2);
1984  local_f += JxW[qp] * vec3_n2;
1985 
1986  // von Karman strain: w-displacement
1987  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
1988  Bmat_w_vk.vector_mult_transpose(vec3_n2, vec4_2);
1989  local_f += JxW[qp] * vec3_n2;
1990  }
1991 
1992  if (request_jacobian && if_vk) { // Jacobian only for vk strain
1993 
1994  // vk - vk: v-displacement
1995  mat3 = RealMatrixX::Zero(2, n2);
1996  Bmat_v_vk.left_multiply(mat3, stress);
1997  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
1998  local_jac += JxW[qp] * mat2_n2n2;
1999 
2000  // vk - vk: w-displacement
2001  mat3 = RealMatrixX::Zero(2, n2);
2002  Bmat_w_vk.left_multiply(mat3, stress);
2003  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
2004  local_jac += JxW[qp] * mat2_n2n2;
2005  }
2006  }
2007  }
2008 
2009 
2010  // now transform to the global coorodinate system
2011  transform_vector_to_global_system(local_f, vec3_n2);
2012  f -= vec3_n2;
2013  if (request_jacobian && if_vk) {
2014  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2015  jac -= scaling * mat2_n2n2;
2016  }
2017 
2018  // Jacobian contribution from von Karman strain
2019  return request_jacobian;
2020 }
2021 
2022 
2023 
2024 
2025 bool
2026 MAST::StructuralElement1D::thermal_residual_sensitivity (const MAST::FunctionBase& p,
2027  bool request_jacobian,
2028  RealVectorX& f,
2029  RealMatrixX& jac,
2030  MAST::BoundaryConditionBase& bc)
2031 {
2032  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true,
2033  false,
2034  _property.extra_quadrature_order(_elem)));
2035 
2036  const std::vector<Real>& JxW = fe->get_JxW();
2037  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
2038  const unsigned int
2039  n_phi = (unsigned int)fe->get_phi().size(),
2040  n1 = this->n_direct_strain_components(),
2041  n2 = 6*n_phi,
2042  n3 = this->n_von_karman_strain_components();
2043 
2044  RealMatrixX
2045  material_exp_A_mat,
2046  material_exp_B_mat,
2047  material_exp_A_mat_sens,
2048  material_exp_B_mat_sens,
2049  mat1_n1n2 = RealMatrixX::Zero(n1,n2),
2050  mat2_n2n2 = RealMatrixX::Zero(n2,n2),
2051  mat3,
2052  mat4_n3n2 = RealMatrixX::Zero(n3,n2),
2053  vk_dvdxi_mat = RealMatrixX::Zero(2,2),
2054  vk_dwdxi_mat = RealMatrixX::Zero(2,2),
2055  stress = RealMatrixX::Zero(2,2),
2056  local_jac = RealMatrixX::Zero(n2,n2);
2057  RealVectorX
2058  vec1_n1 = RealVectorX::Zero(n1),
2059  vec2_n1 = RealVectorX::Zero(n1),
2060  vec3_n2 = RealVectorX::Zero(n2),
2061  vec4_2 = RealVectorX::Zero(2),
2062  vec5_n1 = RealVectorX::Zero(n1),
2063  local_f = RealVectorX::Zero(n2),
2064  delta_t = RealVectorX::Zero(1),
2065  delta_t_sens = RealVectorX::Zero(1);
2066 
2067  MAST::FEMOperatorMatrix
2068  Bmat_mem,
2069  Bmat_bend_v,
2070  Bmat_bend_w,
2071  Bmat_v_vk,
2072  Bmat_w_vk;
2073 
2074  Bmat_mem.reinit(n1, _system.n_vars(), n_phi); // three stress-strain components
2075  Bmat_bend_v.reinit(n1, _system.n_vars(), n_phi);
2076  Bmat_bend_w.reinit(n1, _system.n_vars(), n_phi);
2077  Bmat_v_vk.reinit(n3, _system.n_vars(), n_phi); // only dv/dx and dv/dy
2078  Bmat_w_vk.reinit(n3, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2079 
2080  bool if_vk = (_property.strain_type() == MAST::NONLINEAR_STRAIN);
2081 
2082  MAST::BendingOperatorType
2083  bending_model = _property.bending_model(_elem);
2084 
2085  std::unique_ptr<MAST::BendingOperator1D> bend;
2086 
2087  if (bending_model != MAST::NO_BENDING)
2088  bend.reset(MAST::build_bending_operator_1D(bending_model,
2089  *this,
2090  fe->get_qpoints()).release());
2091 
2092  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
2093  expansion_A = _property.thermal_expansion_A_matrix(*this),
2094  expansion_B = _property.thermal_expansion_B_matrix(*this);
2095 
2096  // temperature function
2097  const MAST::FieldFunction<Real>
2098  &temp_func = bc.get<MAST::FieldFunction<Real> >("temperature"),
2099  &ref_temp_func = bc.get<MAST::FieldFunction<Real> >("ref_temperature");
2100 
2101  Real t, t0, t_sens;
2102 
2103  for (unsigned int qp=0; qp<JxW.size(); qp++) {
2104 
2105  // get the material property
2106  (*expansion_A)(xyz[qp], _time, material_exp_A_mat);
2107  (*expansion_B)(xyz[qp], _time, material_exp_B_mat);
2108  expansion_A->derivative(p, xyz[qp], _time, material_exp_A_mat_sens);
2109  expansion_B->derivative(p, xyz[qp], _time, material_exp_B_mat_sens);
2110 
2111  // get the temperature function
2112  temp_func(xyz[qp], _time, t);
2113  temp_func.derivative(p, xyz[qp], _time, t_sens);
2114  ref_temp_func(xyz[qp], _time, t0);
2115  delta_t(0) = t-t0;
2116  delta_t_sens(0) = t_sens;
2117 
2118  // now prepare the membrane force sensitivity
2119  vec1_n1 = material_exp_A_mat * delta_t_sens; // [C]{alpha (dT/dp)} (with membrane strain)
2120  vec2_n1 = material_exp_A_mat_sens * delta_t; // d([C].{alpha})/dp (T - T0)} (with membrane
2121  vec1_n1 += vec2_n1; // sensitivity of the thermal membrane force
2122  stress(0,0) = vec1_n1(0); // sigma_xx
2123 
2124  // now prepare the membrane-bending coupling force sensitivity
2125  vec2_n1 = material_exp_B_mat * delta_t_sens; // [C]{alpha dT/dp} (with bending strain)
2126  vec5_n1 = material_exp_B_mat_sens * delta_t; // d([C].{alpha})/dp (T - T0)} (with bending strain)
2127  vec2_n1 += vec5_n1; // sensitivity of the thermal membrane force
2128 
2129 
2130  this->initialize_direct_strain_operator(qp, *fe, Bmat_mem);
2131 
2132  // membrane strain
2133  Bmat_mem.vector_mult_transpose(vec3_n2, vec1_n1);
2134  local_f += JxW[qp] * vec3_n2;
2135 
2136  if (bend.get()) {
2137  // bending strain
2138  bend->initialize_bending_strain_operator(*fe, qp,
2139  Bmat_bend_v,
2140  Bmat_bend_w);
2141  Bmat_bend_v.vector_mult_transpose(vec3_n2, vec2_n1);
2142  local_f += JxW[qp] * vec3_n2;
2143  Bmat_bend_w.vector_mult_transpose(vec3_n2, vec2_n1);
2144  local_f += JxW[qp] * vec3_n2;
2145 
2146  // von Karman strain
2147  if (if_vk) {
2148  // get the vonKarman strain operator if needed
2149  this->initialize_von_karman_strain_operator(qp,
2150  *fe,
2151  vec2_n1, // epsilon_vk
2152  vk_dvdxi_mat,
2153  vk_dwdxi_mat,
2154  Bmat_v_vk,
2155  Bmat_w_vk);
2156  // von Karman strain: v-displacement
2157  vec4_2 = vk_dvdxi_mat.transpose() * vec1_n1;
2158  Bmat_v_vk.vector_mult_transpose(vec3_n2, vec4_2);
2159  local_f += JxW[qp] * vec3_n2;
2160 
2161  // von Karman strain: w-displacement
2162  vec4_2 = vk_dwdxi_mat.transpose() * vec1_n1;
2163  Bmat_w_vk.vector_mult_transpose(vec3_n2, vec4_2);
2164  local_f += JxW[qp] * vec3_n2;
2165  }
2166 
2167  if (request_jacobian && if_vk) { // Jacobian only for vk strain
2168  // vk - vk: v-displacement
2169  mat3 = RealMatrixX::Zero(2, n2);
2170  Bmat_v_vk.left_multiply(mat3, stress);
2171  Bmat_v_vk.right_multiply_transpose(mat2_n2n2, mat3);
2172  local_jac += JxW[qp] * mat2_n2n2;
2173 
2174  // vk - vk: w-displacement
2175  mat3 = RealMatrixX::Zero(2, n2);
2176  Bmat_w_vk.left_multiply(mat3, stress);
2177  Bmat_w_vk.right_multiply_transpose(mat2_n2n2, mat3);
2178  local_jac += JxW[qp] * mat2_n2n2;
2179  }
2180  }
2181  }
2182 
2183 
2184  // now transform to the global coorodinate system
2185  transform_vector_to_global_system(local_f, vec3_n2);
2186  f -= vec3_n2;
2187  if (request_jacobian && if_vk) {
2188  transform_matrix_to_global_system(local_jac, mat2_n2n2);
2189  jac -= mat2_n2n2;
2190  }
2191 
2192  // Jacobian contribution from von Karman strain
2193  return request_jacobian;
2194 }
2195 
2196 
2197 
2198 
2199 
2200 
2201 
2202 bool
2203 MAST::StructuralElement1D::
2204 piston_theory_residual(bool request_jacobian,
2205  RealVectorX &f,
2206  RealMatrixX& jac_xdot,
2207  RealMatrixX& jac,
2208  MAST::BoundaryConditionBase& bc) {
2209 
2210  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2211 
2212  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
2213 
2214  const std::vector<Real> &JxW = fe->get_JxW();
2215  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2216  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2217  const unsigned int
2218  n_phi = (unsigned int)phi.size(),
2219  n1 = this->n_direct_strain_components(),
2220  n2 = 6*n_phi;
2221 
2222 
2223 
2224  // convert to piston theory boundary condition so that the necessary
2225  // flow properties can be obtained
2226  const MAST::PistonTheoryBoundaryCondition& piston_bc =
2227  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
2228 
2229  // create the constant field functions to pass the dwdx and dwdt values
2230  // to the piston theory pressure functions
2231  MAST::Parameter
2232  dwdx_p ("dwdx", 0.),
2233  dwdt_p ("dwdt", 0.);
2234 
2235  MAST::ConstantFieldFunction
2236  dwdx_f ("dwdx", dwdx_p),
2237  dwdt_f ("dwdx", dwdt_p);
2238 
2239  std::unique_ptr<MAST::FieldFunction<Real> >
2240  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
2241  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
2242  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
2243 
2244  FEMOperatorMatrix
2245  Bmat_v, // operator matrix for the v-displacement
2246  Bmat_w, // operator matrix for the w-displacement
2247  Bmat_dvdx, // operator matrix to calculate the derivativ of v wrt x
2248  Bmat_dwdx; // operator matrix to calculate the derivativ of w wrt x
2249 
2250  Bmat_dvdx.reinit(2, _system.n_vars(), n_phi); // only dv/dx and dv/dy
2251  Bmat_dwdx.reinit(2, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2252 
2253  RealVectorX
2254  phi_vec = RealVectorX::Zero(n_phi),
2255  force = RealVectorX::Zero(n1),
2256  local_f = RealVectorX::Zero(n2),
2257  vec_n1 = RealVectorX::Zero(n1),
2258  vec_n2 = RealVectorX::Zero(n2),
2259  vel = RealVectorX::Zero(3),
2260  p_val = RealVectorX::Zero(3),
2261  normal = RealVectorX::Zero(3),
2262  dummy = RealVectorX::Zero(2);
2263 
2264 
2265  // direction of pressure assumed to be normal (along local z-axis)
2266  // to the element face for 2D and along local y-axis for 1D element.
2267  normal(_elem.dim()) = -1.;
2268 
2269 
2270  RealMatrixX
2271  dvdx = RealMatrixX::Zero(2,2),
2272  dwdx = RealMatrixX::Zero(2,2),
2273  local_jac = RealMatrixX::Zero(n2,n2),
2274  local_jac_xdot = RealMatrixX::Zero(n2,n2),
2275  mat_n2n2 = RealMatrixX::Zero(n2,n2);
2276 
2277 
2278  for (unsigned int qp=0; qp<qpoint.size(); qp++)
2279  {
2280 
2281  // now set the shape function values
2282  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2283  phi_vec(i_nd) = phi[i_nd][qp];
2284 
2285  Bmat_v.reinit(n1, _system.n_vars(), n_phi);
2286  Bmat_v.set_shape_function(0, 1, phi_vec);
2287  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
2288  Bmat_w.set_shape_function(1, 2, phi_vec);
2289 
2290  // use the Bmat to calculate the velocity vector
2291  Bmat_v.right_multiply(vec_n1, _local_vel);
2292  vel(1) = vec_n1(0); // v_dot
2293  Bmat_w.right_multiply(vec_n1, _local_vel);
2294  vel(2) = vec_n1(1); // w_dot
2295 
2296  // get the operators for dv/dx and dw/dx. We will use the
2297  // von Karman strain operators for this
2298  initialize_von_karman_strain_operator(qp,
2299  *fe,
2300  dummy,
2301  dvdx,
2302  dwdx,
2303  Bmat_dvdx,
2304  Bmat_dwdx);
2305 
2306 
2307  // now use this information to evaluate the normal cp due to
2308  // both the v and w motions
2309  dwdx_p = dvdx(0,0);
2310  dwdt_p = vel(1);
2311  (*pressure)(qpoint[qp], _time, p_val(1));
2312 
2313 
2314  dwdx_p = dwdx(0,0);
2315  dwdt_p = vel(2);
2316  (*pressure)(qpoint[qp], _time, p_val(2));
2317 
2318 
2319  // calculate force and discrete vector for v-disp
2320  force.setZero();
2321  force(0) = p_val(1) * normal(1);
2322  Bmat_v.vector_mult_transpose(vec_n2, force);
2323  local_f += JxW[qp] * vec_n2;
2324 
2325  // now the w-disp
2326  force.setZero();
2327  force(1) = p_val(2) * normal(2);
2328  Bmat_v.vector_mult_transpose(vec_n2, force);
2329  local_f += JxW[qp] * vec_n2;
2330 
2331 
2332  // calculate the Jacobian if requested
2333  if (request_jacobian) {
2334 
2335  // we need the derivative of cp wrt normal velocity
2336  dwdx_p = dvdx(0,0);
2337  dwdt_p = vel(1);
2338  (*dpressure_dxdot)(qpoint[qp], _time, p_val(1));
2339 
2340  // calculate the component of Jacobian due to v-velocity
2341  Bmat_v.right_multiply_transpose(mat_n2n2, Bmat_v);
2342  local_jac_xdot += JxW[qp] * p_val(1) * normal(1) * mat_n2n2;
2343 
2344  // now wrt v
2345  (*dpressure_dx)(qpoint[qp], _time, p_val(1));
2346  // now use calculate the component of Jacobian due to deformation
2347  Bmat_v.right_multiply_transpose(mat_n2n2, Bmat_dvdx); // v: B^T dB/dx
2348  local_jac += (JxW[qp] * p_val(1) * normal(1)) * mat_n2n2;
2349 
2350 
2351  dwdx_p = dwdx(0,0);
2352  dwdt_p = vel(2);
2353  (*dpressure_dxdot)(qpoint[qp], _time, p_val(2));
2354 
2355  // calculate the component of Jacobian due to w-velocity
2356  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
2357  local_jac_xdot += JxW[qp] * p_val(2) * normal(2) * mat_n2n2;
2358 
2359  // now wrt w
2360  (*dpressure_dx)(qpoint[qp], _time, p_val(2));
2361  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_dwdx); // w: B^T dB/dx
2362  local_jac += (JxW[qp] * p_val(2) * normal(2)) * mat_n2n2;
2363  }
2364  }
2365 
2366  // now transform to the global system and add
2367  transform_vector_to_global_system(local_f, vec_n2);
2368  f -= vec_n2;
2369 
2370  // if the Jacobian was requested, then transform it and add to the
2371  // global Jacobian
2372  if (request_jacobian) {
2373  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
2374  jac_xdot -= mat_n2n2;
2375 
2376  transform_matrix_to_global_system(local_jac, mat_n2n2);
2377  jac -= mat_n2n2;
2378  }
2379 
2380  return (request_jacobian);
2381 }
2382 
2383 
2384 
2385 
2386 
2387 
2388 
2389 
2390 
2391 
2392 bool
2393 MAST::StructuralElement1D::
2394 piston_theory_residual_sensitivity(const MAST::FunctionBase& p,
2395  bool request_jacobian,
2396  RealVectorX &f,
2397  RealMatrixX& jac_xdot,
2398  RealMatrixX& jac,
2399  MAST::BoundaryConditionBase& bc) {
2400 
2401 
2402  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2403 
2404  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(true, false));
2405 
2406  const std::vector<Real> &JxW = fe->get_JxW();
2407  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2408  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2409  const unsigned int
2410  n_phi = (unsigned int)phi.size(),
2411  n1 = this->n_direct_strain_components(),
2412  n2 = 6*n_phi;
2413 
2414 
2415 
2416  // convert to piston theory boundary condition so that the necessary
2417  // flow properties can be obtained
2418  const MAST::PistonTheoryBoundaryCondition& piston_bc =
2419  dynamic_cast<MAST::PistonTheoryBoundaryCondition&>(bc);
2420 
2421  // create the constant field functions to pass the dwdx and dwdt values
2422  // to the piston theory pressure functions
2423  MAST::Parameter
2424  dwdx_p ("dwdx", 0.),
2425  dwdt_p ("dwdt", 0.);
2426 
2427  MAST::ConstantFieldFunction
2428  dwdx_f ("dwdx", dwdx_p),
2429  dwdt_f ("dwdx", dwdt_p);
2430 
2431  std::unique_ptr<MAST::FieldFunction<Real> >
2432  pressure (piston_bc.get_pressure_function(dwdx_f, dwdt_f).release()),
2433  dpressure_dx (piston_bc.get_dpdx_function (dwdx_f, dwdt_f).release()),
2434  dpressure_dxdot (piston_bc.get_dpdxdot_function (dwdx_f, dwdt_f).release());
2435 
2436  FEMOperatorMatrix
2437  Bmat_v, // operator matrix for the v-displacement
2438  Bmat_w, // operator matrix for the w-displacement
2439  Bmat_dvdx, // operator matrix to calculate the derivativ of v wrt x
2440  Bmat_dwdx; // operator matrix to calculate the derivativ of w wrt x
2441 
2442  Bmat_dvdx.reinit(2, _system.n_vars(), n_phi); // only dv/dx and dv/dy
2443  Bmat_dwdx.reinit(2, _system.n_vars(), n_phi); // only dw/dx and dw/dy
2444 
2445  RealVectorX
2446  phi_vec = RealVectorX::Zero(n_phi),
2447  force = RealVectorX::Zero(n1),
2448  local_f = RealVectorX::Zero(n2),
2449  vec_n1 = RealVectorX::Zero(n1),
2450  vec_n2 = RealVectorX::Zero(n2),
2451  vel = RealVectorX::Zero(3),
2452  p_val = RealVectorX::Zero(3),
2453  normal = RealVectorX::Zero(3),
2454  dummy = RealVectorX::Zero(2);
2455 
2456 
2457  // direction of pressure assumed to be normal (along local z-axis)
2458  // to the element face for 2D and along local y-axis for 1D element.
2459  normal(_elem.dim()) = -1.;
2460 
2461 
2462  RealMatrixX
2463  dvdx = RealMatrixX::Zero(2,2),
2464  dwdx = RealMatrixX::Zero(2,2),
2465  local_jac = RealMatrixX::Zero(n2,n2),
2466  local_jac_xdot = RealMatrixX::Zero(n2,n2),
2467  mat_n2n2 = RealMatrixX::Zero(n2,n2);
2468 
2469 
2470  for (unsigned int qp=0; qp<qpoint.size(); qp++)
2471  {
2472 
2473  // now set the shape function values
2474  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2475  phi_vec(i_nd) = phi[i_nd][qp];
2476 
2477  Bmat_v.reinit(n1, _system.n_vars(), n_phi);
2478  Bmat_v.set_shape_function(0, 1, phi_vec);
2479  Bmat_w.reinit(n1, _system.n_vars(), n_phi);
2480  Bmat_w.set_shape_function(1, 2, phi_vec);
2481 
2482  // use the Bmat to calculate the velocity vector
2483  Bmat_v.right_multiply(vec_n1, _local_vel);
2484  vel(1) = vec_n1(0); // v_dot
2485  Bmat_w.right_multiply(vec_n1, _local_vel);
2486  vel(2) = vec_n1(1); // w_dot
2487 
2488  // get the operators for dv/dx and dw/dx. We will use the
2489  // von Karman strain operators for this
2490  initialize_von_karman_strain_operator(qp,
2491  *fe,
2492  dummy,
2493  dvdx,
2494  dwdx,
2495  Bmat_dvdx,
2496  Bmat_dwdx);
2497 
2498 
2499  // now use this information to evaluate the normal cp due to
2500  // both the v and w motions
2501  dwdx_p = dvdx(0,0);
2502  dwdt_p = vel(1);
2503  pressure->derivative(p, qpoint[qp], _time, p_val(1));
2504 
2505 
2506  dwdx_p = dwdx(0,0);
2507  dwdt_p = vel(2);
2508  pressure->derivative(p, qpoint[qp], _time, p_val(2));
2509 
2510 
2511  // calculate force and discrete vector for v-disp
2512  force.setZero();
2513  force(0) = p_val(1) * normal(1);
2514  Bmat_v.vector_mult_transpose(vec_n2, force);
2515  local_f += JxW[qp] * vec_n2;
2516 
2517  // now the w-disp
2518  force.setZero();
2519  force(1) = p_val(2) * normal(2);
2520  Bmat_v.vector_mult_transpose(vec_n2, force);
2521  local_f += JxW[qp] * vec_n2;
2522 
2523 
2524  // calculate the Jacobian if requested
2525  if (request_jacobian) {
2526 
2527  // we need the derivative of cp wrt normal velocity
2528  dwdx_p = dvdx(0,0);
2529  dwdt_p = vel(1);
2530  dpressure_dxdot->derivative(p, qpoint[qp], _time, p_val(1));
2531 
2532  // calculate the component of Jacobian due to v-velocity
2533  Bmat_v.right_multiply_transpose(mat_n2n2, Bmat_v);
2534  local_jac_xdot += JxW[qp] * p_val(1) * normal(1) * mat_n2n2;
2535 
2536  // now wrt v
2537  dpressure_dx->derivative(p, qpoint[qp], _time, p_val(1));
2538 
2539  // now use calculate the component of Jacobian due to deformation
2540  Bmat_v.right_multiply_transpose(mat_n2n2, Bmat_dvdx); // v: B^T dB/dx
2541  local_jac += (JxW[qp] * p_val(1) * normal(1)) * mat_n2n2;
2542 
2543 
2544  dwdx_p = dwdx(0,0);
2545  dwdt_p = vel(2);
2546  dpressure_dxdot->derivative(p, qpoint[qp], _time, p_val(2));
2547 
2548  // calculate the component of Jacobian due to w-velocity
2549  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_w);
2550  local_jac_xdot += JxW[qp] * p_val(2) * normal(2) * mat_n2n2;
2551 
2552  // now wrt w
2553  dpressure_dx->derivative(p, qpoint[qp], _time, p_val(2));
2554  Bmat_w.right_multiply_transpose(mat_n2n2, Bmat_dwdx); // w: B^T dB/dx
2555  local_jac += (JxW[qp] * p_val(2) * normal(2)) * mat_n2n2;
2556  }
2557  }
2558 
2559  // now transform to the global system and add
2560  transform_vector_to_global_system(local_f, vec_n2);
2561  f -= vec_n2;
2562 
2563  // if the Jacobian was requested, then transform it and add to the
2564  // global Jacobian
2565  if (request_jacobian) {
2566  transform_matrix_to_global_system(local_jac_xdot, mat_n2n2);
2567  jac_xdot -= mat_n2n2;
2568 
2569  transform_matrix_to_global_system(local_jac, mat_n2n2);
2570  jac -= mat_n2n2;
2571  }
2572 
2573 
2574  // no parametric sensitivity is calculated for piston theory at this point.
2575  return (request_jacobian);
2576 }
2577 
2578 
2579 
2580 
2581 
2582 /*bool
2583 MAST::StructuralElement1D::
2584 linearized_frequency_domain_surface_pressure_residual
2585  (bool request_jacobian,
2586  ComplexVectorX &f,
2587  ComplexMatrixX &jac,
2588  const unsigned int side,
2589  MAST::BoundaryConditionBase& bc) {
2590 
2591  libmesh_assert(!follower_forces); // not implemented yet for follower forces
2592  libmesh_assert_equal_to(bc.type(), MAST::SMALL_DISTURBANCE_MOTION);
2593 
2594 
2595  MAST::FieldFunction<Real>&
2596  press_fn = bc.get<MAST::FieldFunction<Real> > ("pressure");
2597  MAST::FieldFunction<Complex>&
2598  dpress_fn = bc.get<MAST::FieldFunction<Complex> > ("dpressure");
2599  MAST::FieldFunction<ComplexVectorX>&
2600  dn_rot_fn = bc.get<MAST::FieldFunction<ComplexVectorX> >("dnormal");
2601 
2602 
2603  libMesh::FEBase *fe_ptr = nullptr;
2604  libMesh::QBase *qrule_ptr = nullptr;
2605  _get_side_fe_and_qrule(get_elem_for_quadrature(),
2606  side,
2607  &fe_ptr,
2608  &qrule_ptr,
2609  false);
2610  std::unique_ptr<libMesh::FEBase> fe(fe_ptr);
2611  std::unique_ptr<libMesh::QBase> qrule(qrule_ptr);
2612 
2613 
2614  // Physical location of the quadrature points
2615  const std::vector<Real> &JxW = fe->get_JxW();
2616  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
2617  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2618  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
2619 
2620  const unsigned int
2621  n_phi = (unsigned int)phi.size(),
2622  n1 = 3,
2623  n2 = 6*n_phi;
2624 
2625  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
2626  ComplexVectorX
2627  dn_rot = ComplexVectorX::Zero(3),
2628  force = ComplexVectorX::Zero(2*n1),
2629  local_f = ComplexVectorX::Zero(n2),
2630  vec_n2 = ComplexVectorX::Zero(n2);
2631 
2632 
2633  FEMOperatorMatrix Bmat;
2634  Real press;
2635  Complex dpress;
2636 
2637  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
2638 
2639  // now set the shape function values
2640  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
2641  phi_vec(i_nd) = phi[i_nd][qp];
2642 
2643  Bmat.reinit(2*n1, phi_vec);
2644 
2645  // get pressure and deformation information
2646  press_fn (qpoint[qp], _time, press);
2647  dpress_fn(qpoint[qp], _time, dpress);
2648  dn_rot_fn(qpoint[qp], _time, dn_rot);
2649 
2650  // press = 0.;
2651  // dpress = Complex(2./4.*std::real(dn_rot(0)), 2./4./.1*std::imag(utrans(1)));
2652  // libMesh::out << q_point[qp](0)
2653  // << " " << std::real(utrans(1))
2654  // << " " << std::imag(utrans(1))
2655  // << " " << std::real(dn_rot(0))
2656  // << " " << std::imag(dn_rot(0))
2657  // << " " << std::real(press)
2658  // << " " << std::imag(press)
2659  // << " " << std::real(dpress)
2660  // << " " << std::imag(dpress) << std::endl;
2661 
2662  // calculate force
2663  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
2664  force(i_dim) = ( press * dn_rot(i_dim) + // steady pressure
2665  dpress * face_normals[qp](i_dim)); // unsteady pressure
2666 
2667 
2668  Bmat.vector_mult_transpose(vec_n2, force);
2669 
2670  local_f -= JxW[qp] * vec_n2;
2671  }
2672 
2673  // now transform to the global system and add
2674  transform_vector_to_global_system(local_f, vec_n2);
2675  f += vec_n2;
2676 
2677 
2678  return (request_jacobian);
2679 }
2680 */
2681 
2682 
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.
structural_element_1d.h
MAST::StructuralElement1D::initialize_direct_strain_operator
virtual void initialize_direct_strain_operator(const unsigned int qp, const MAST::FEBase &fe, MAST::FEMOperatorMatrix &Bmat)
initialize membrane strain operator matrix
Definition: structural_element_1d.cpp:56
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::build_bending_operator_1D
std::unique_ptr< MAST::BendingOperator1D > build_bending_operator_1D(MAST::BendingOperatorType type, MAST::StructuralElementBase &elem, const std::vector< libMesh::Point > &pts)
builds a bending operator and returns it in a smart-pointer
Definition: bending_operator.cpp:46
MAST::StructuralElement1D::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_1d.h:54
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::StructuralElement1D::_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_1d.cpp:1373
element_property_card_1D.h
MAST::ElementPropertyCardBase::prestress_B_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > prestress_B_matrix(MAST::ElementBase &e) const =0
MAST::FEMOperatorMatrix::reinit
void reinit(unsigned int n_interpolated_vars, unsigned int n_discrete_vars, unsigned int n_discrete_dofs_per_var)
this initializes the operator for number of rows and variables, assuming that all variables has the s...
Definition: fem_operator_matrix.h:210
MAST::StructuralElement1D::prestress_residual
virtual bool prestress_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal force vector and Jacobian due to strain energy coming from a prestress...
Definition: structural_element_1d.cpp:1386
MAST::StructuralElement1D::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_1d.cpp:2394
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::StructuralElement1D::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_1d.cpp:2204
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::StructuralElement1D::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 force vector and the Jacobian due to thermal stresses.
Definition: structural_element_1d.cpp:2026
fem_operator_matrix.h
MAST::StructuralElement1D::_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, MAST::BendingOperator1D *bend_op, MAST::FEMOperatorMatrix &Bmat_mem, MAST::FEMOperatorMatrix &Bmat_bend_v, MAST::FEMOperatorMatrix &Bmat_bend_w, MAST::FEMOperatorMatrix &Bmat_v_vk, MAST::FEMOperatorMatrix &Bmat_w_vk, RealMatrixX &stress, RealMatrixX &stress_l, RealMatrixX &vk_dvdxi_mat, 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, RealMatrixX &mat1_n1n2, RealMatrixX &mat2_n2n2, RealMatrixX &mat3, RealMatrixX &mat4_2n2)
performs integration at the quadrature point for the provided matrices.
Definition: structural_element_1d.cpp:1074
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
stress_output_base.h
geom_elem.h
MAST::FunctionSetBase::depends_on
virtual bool depends_on(const MAST::FunctionBase &f) const
returns true if the property card depends on the function f
Definition: function_set_base.cpp:61
MAST::StressStrainOutputBase::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::StructuralElement1D::initialize_von_karman_strain_operator_sensitivity
void initialize_von_karman_strain_operator_sensitivity(const unsigned int qp, const MAST::FEBase &fe, RealMatrixX &vk_dvdxi_mat_sens, RealMatrixX &vk_dwdxi_mat_sens)
initialze the sensitivity of von Karman operator matrices needed for Jacobian calculation.
Definition: structural_element_1d.cpp:128
MAST::StructuralElement1D::internal_residual_jac_dot_state_sensitivity
virtual bool internal_residual_jac_dot_state_sensitivity(RealMatrixX &jac)
calculates d[J]/d{x} .
Definition: structural_element_1d.cpp:800
MAST::FieldFunction::derivative
virtual void derivative(const MAST::FunctionBase &f, ValType &v) const
calculates the value of the function derivative and returns it in v.
Definition: field_function_base.h:78
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
RealMatrixX
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Definition: mast_data_types.h:44
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::FEMOperatorMatrix::vector_mult
void vector_mult(T &res, const T &v) const
res = [this] * v
Definition: fem_operator_matrix.h:280
MAST::PistonTheoryBoundaryCondition::get_dpdx_function
std::unique_ptr< MAST::FieldFunction< Real > > get_dpdx_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
Definition: piston_theory_boundary_condition.cpp:558
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::BendingOperator1D
Bending strain operator for 1D element.
Definition: bending_operator.h:123
boundary_condition_base.h
MAST::StructuralElement1D::prestress_residual_sensitivity
virtual bool prestress_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal force vector and Jacobian due to strain energy coming from a prestress...
Definition: structural_element_1d.cpp:1541
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
assembly_base.h
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::PistonTheoryBoundaryCondition::get_dpdxdot_function
std::unique_ptr< MAST::FieldFunction< Real > > get_dpdxdot_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
Definition: piston_theory_boundary_condition.cpp:579
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::StructuralElement1D::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_1d.cpp:1775
MAST::FEBase::get_dphi
virtual const std::vector< std::vector< libMesh::RealVectorValue > > & get_dphi() const
Definition: fe_base.cpp:255
MAST::PistonTheoryBoundaryCondition::get_pressure_function
std::unique_ptr< MAST::FieldFunction< Real > > get_pressure_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
Definition: piston_theory_boundary_condition.cpp:537
MAST::BendingOperator1D::initialize_bending_strain_operator
virtual void initialize_bending_strain_operator(const MAST::FEBase &fe, const unsigned int qp, MAST::FEMOperatorMatrix &Bmat_v, MAST::FEMOperatorMatrix &Bmat_w)=0
initialze the bending strain operator for the specified quadrature point
MAST::GeomElem::dim
unsigned int dim() const
Definition: geom_elem.cpp:91
MAST::StructuralElement1D::_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_1d.cpp:1362
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
MAST::ElementPropertyCardBase::thermal_expansion_A_matrix
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > thermal_expansion_A_matrix(const MAST::ElementBase &e) const =0
MAST::StructuralElement1D::n_von_karman_strain_components
virtual unsigned int n_von_karman_strain_components()
row dimension of the von Karman strain matrix
Definition: structural_element_1d.h:61
MAST::StructuralElement1D::thermal_residual
virtual bool thermal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to thermal stresses.
Definition: structural_element_1d.cpp:1858
MAST::StructuralElement1D::calculate_stress
virtual bool calculate_stress(bool request_derivative, const MAST::FunctionBase *p, MAST::StressStrainOutputBase &output)
Calculates the stress tensor.
Definition: structural_element_1d.cpp:162
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::StructuralElementBase::_local_sol
RealVectorX _local_sol
local solution
Definition: structural_element_base.h:707
MAST::StructuralElementBase::_local_vel
RealVectorX _local_vel
local velocity
Definition: structural_element_base.h:731
MAST::StructuralElement1D::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_dvdxi_mat, RealMatrixX &vk_dwdxi_mat, MAST::FEMOperatorMatrix &Bmat_v_vk, MAST::FEMOperatorMatrix &Bmat_w_vk)
initialze the von Karman strain in vK_strain, the operator matrices needed for Jacobian calculation...
Definition: structural_element_1d.cpp:81
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::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::StructuralElement1D::StructuralElement1D
StructuralElement1D(MAST::SystemInitialization &sys, MAST::AssemblyBase &assembly, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
Definition: structural_element_1d.cpp:37
MAST::ElementBase::_time
const Real & _time
time for which system is being assembled
Definition: elem_base.h:232
MAST::StructuralElement1D::internal_residual_sensitivity
virtual bool internal_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the sensitivity internal force vector and Jacobian due to strain energy.
Definition: structural_element_1d.cpp:663
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::StructuralElement1D::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_1d.cpp:1697
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::StructuralElement1D::internal_residual
virtual bool internal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the internal force vector and Jacobian due to strain energy.
Definition: structural_element_1d.cpp:539
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