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