MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
structural_element_base.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
28 #include "base/nonlinear_system.h"
29 #include "base/assembly_base.h"
32 #include "mesh/geom_elem.h"
33 #include "mesh/fe_base.h"
35 #include "numerics/utility.h"
36 
37 
39  const MAST::GeomElem& elem,
41 MAST::ElementBase(sys, elem),
42 follower_forces (false),
43 _property (p),
44 _incompatible_sol (nullptr) {
45 
46 }
47 
48 
49 
51 
52 }
53 
54 
55 
56 void
58  bool if_sens) {
59 
60  // convert the vector to the local element coordinate system
61  if (!if_sens) {
62  if (_elem.dim() == 3)
63  _local_sol = vec;
64  else {
65  _local_sol = RealVectorX::Zero(vec.size());
67  }
68  }
69  else {
70 
71  // set the element solution sensitivity.
72  if (_elem.dim() == 3)
73  _local_sol_sens = vec;
74  else {
75  _local_sol_sens = RealVectorX::Zero(vec.size());
77  }
78  }
79 
80  MAST::ElementBase::set_solution(vec, if_sens);
81 }
82 
83 
84 
85 void
87  bool if_sens) {
88 
89  // convert the vector to the local element coordinate system
90  if (!if_sens) {
91  if (_elem.dim() == 3)
92  _local_delta_sol = vec;
93  else {
94  _local_delta_sol = RealVectorX::Zero(vec.size());
96  }
97  }
98  else {
99 
100  // set the element solution sensitivity.
101  if (_elem.dim() == 3)
102  _local_delta_sol_sens = vec;
103  else {
104  _local_delta_sol_sens = RealVectorX::Zero(vec.size());
106  }
107  }
108 
110 }
111 
112 
113 
114 
115 void
117  bool if_sens) {
118 
119  if (!if_sens) {
120  if (_elem.dim() == 3)
121  _local_vel = vec;
122  else {
123  _local_vel = RealVectorX::Zero(vec.size());
125  }
126  }
127  else {
128 
129  if (_elem.dim() == 3)
130  _local_vel_sens = vec;
131  else {
132  _local_vel_sens = RealVectorX::Zero(vec.size());
134  }
135  }
136 
137  MAST::ElementBase::set_velocity(vec, if_sens);
138 }
139 
140 
141 
142 void
144  bool if_sens) {
145 
146  if (!if_sens) {
147  if (_elem.dim() == 3)
148  _local_delta_vel = vec;
149  else {
150  _local_delta_vel = RealVectorX::Zero(vec.size());
152  }
153  }
154  else {
155 
156  if (_elem.dim() == 3)
157  _local_delta_vel_sens = vec;
158  else {
159  _local_delta_vel_sens = RealVectorX::Zero(vec.size());
161  }
162  }
163 
165 }
166 
167 
168 
169 void
171  bool if_sens) {
172 
173  if (!if_sens) {
174  if (_elem.dim() == 3)
175  _local_accel = vec;
176  else {
177  _local_accel = RealVectorX::Zero(vec.size());
179  }
180  }
181  else {
182 
183  if (_elem.dim() == 3)
184  _local_accel_sens = vec;
185  else {
186  _local_accel_sens = RealVectorX::Zero(vec.size());
188  }
189  }
190 
192 }
193 
194 
195 
196 void
198  bool if_sens) {
199 
200  if (!if_sens) {
201  if (_elem.dim() == 3)
202  _local_delta_accel = vec;
203  else {
204  _local_delta_accel = RealVectorX::Zero(vec.size());
206  }
207  }
208  else {
209 
210  if (_elem.dim() == 3)
212  else {
213  _local_delta_accel_sens = RealVectorX::Zero(vec.size());
215  }
216  }
217 
219 }
220 
221 
222 
223 bool
225  RealVectorX& f,
226  RealMatrixX& jac) {
227 
228 
229  // It is assumed that the structural elements implement the Jacobians.
230  // Hence, the residual for the linearized problem will be calculated
231  // using the Jacobian.
232 
233  const unsigned int
234  n_phi = (unsigned int)f.size(),
235  n2 =6*n_phi;
236 
238  v = RealVectorX::Zero(n2);
239 
241  m = RealMatrixX::Zero(n2, n2);
242 
243 
244  this->internal_residual(true, v, m);
245 
246  f += m * _local_delta_sol;
247 
248  if (request_jacobian)
249  jac += m;
250 
251  return request_jacobian;
252 }
253 
254 
255 
256 bool
258  RealVectorX& f,
259  RealMatrixX& jac_xddot,
260  RealMatrixX& jac_xdot,
261  RealMatrixX& jac) {
262 
263  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
264 
265  const std::vector<Real>& JxW = fe->get_JxW();
266  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
267  const std::vector<std::vector<Real> >& phi = fe->get_phi();
268 
269  const unsigned int
270  n_phi = (unsigned int)phi.size(),
271  n_vars = _system.system().n_vars(),
272  n1 =6,
273  n2 =6*n_phi;
274 
276  material_mat,
277  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
278  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
279  local_jac = RealMatrixX::Zero(n2, n2);
281  phi_vec = RealVectorX::Zero(n_phi),
282  vec1_n1 = RealVectorX::Zero(n1),
283  vec2_n2 = RealVectorX::Zero(n2),
284  local_f = RealVectorX::Zero(n2);
285 
286  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
287  mat_inertia = _property.inertia_matrix(*this);
288 
290 
292 
293  (*mat_inertia)(xyz[0], _time, material_mat);
294 
295  Real vol = 0.;
296  const unsigned int nshp = fe->n_shape_functions();
297  for (unsigned int i=0; i<JxW.size(); i++)
298  vol += JxW[i];
299  vol /= (1.* nshp);
300  for (unsigned int i_var=0; i_var<6; i_var++)
301  for (unsigned int i=0; i<nshp; i++)
302  local_jac(i_var*nshp+i, i_var*nshp+i) =
303  vol*material_mat(i_var, i_var);
304 
305  local_f = local_jac * _local_accel;
306  }
307  else {
308 
309  for (unsigned int qp=0; qp<JxW.size(); qp++) {
310 
311  (*mat_inertia)(xyz[qp], _time, material_mat);
312 
313  // now set the shape function values
314  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
315  phi_vec(i_nd) = phi[i_nd][qp];
316 
317  Bmat.reinit(n_vars, phi_vec);
318 
319  Bmat.left_multiply(mat1_n1n2, material_mat);
320 
321  vec1_n1 = mat1_n1n2 * _local_accel;
322  Bmat.vector_mult_transpose(vec2_n2, vec1_n1);
323 
324  local_f += JxW[qp] * vec2_n2;
325 
326  if (request_jacobian) {
327 
328  Bmat.right_multiply_transpose(mat2_n2n2,
329  mat1_n1n2);
330  local_jac += JxW[qp]*mat2_n2n2;
331  }
332  }
333  }
334 
335  // now transform to the global coorodinate system
336  if (_elem.dim() < 3) {
337  transform_vector_to_global_system(local_f, vec2_n2);
338  f += vec2_n2;
339 
340  if (request_jacobian) {
341  transform_matrix_to_global_system(local_jac, mat2_n2n2);
342  jac_xddot += mat2_n2n2;
343  }
344  }
345  else {
346 
347  f += local_f;
348  if (request_jacobian)
349  jac_xddot += local_jac;
350  }
351 
352  return request_jacobian;
353 }
354 
355 
356 
357 
358 bool
360 linearized_inertial_residual (bool request_jacobian,
361  RealVectorX& f,
362  RealMatrixX& jac_xddot,
363  RealMatrixX& jac_xdot,
364  RealMatrixX& jac) {
365 
366  // It is assumed that the structural elements implement the Jacobians.
367  // Hence, the residual for the linearized problem will be calculated
368  // using the Jacobian.
369 
370  const unsigned int
371  n_phi = (unsigned int)f.size(),
372  n2 =6*n_phi;
373 
375  v = RealVectorX::Zero(n2);
376 
378  m_xddot = RealMatrixX::Zero(n2, n2),
379  m_xdot = RealMatrixX::Zero(n2, n2),
380  m_x = RealMatrixX::Zero(n2, n2);
381 
382 
383  this->inertial_residual(true, v, m_xddot, m_xdot, m_x);
384 
385  f += (m_xddot * _local_delta_accel +
386  m_xdot * _local_delta_vel +
387  m_x * _local_delta_sol);
388 
389  if (request_jacobian) {
390 
391  jac_xddot += m_xddot;
392  jac_xdot += m_xdot;
393  jac += m_x;
394  }
395 
396  return request_jacobian;
397 }
398 
399 
400 
401 
402 bool
404  bool request_jacobian,
405  RealVectorX& f,
406  RealMatrixX& jac_xddot,
407  RealMatrixX& jac_xdot,
408  RealMatrixX& jac) {
409 
410  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
411 
412  const std::vector<Real>& JxW = fe->get_JxW();
413  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
414  const std::vector<std::vector<Real> >& phi = fe->get_phi();
415 
416  const unsigned int
417  n_phi = (unsigned int)phi.size(),
418  n_vars = _system.system().n_vars(),
419  n1 =6,
420  n2 =6*n_phi;
421 
423  material_mat,
424  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
425  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
426  local_jac = RealMatrixX::Zero(n2, n2);
428  phi_vec = RealVectorX::Zero(n_phi),
429  vec1_n1 = RealVectorX::Zero(n1),
430  vec2_n2 = RealVectorX::Zero(n2),
431  local_f = RealVectorX::Zero(n2);
432 
433  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
434  mat_inertia = _property.inertia_matrix(*this);
435 
437 
439 
440  // as an approximation, get matrix at the first quadrature point
441  mat_inertia->derivative(p, xyz[0], _time, material_mat);
442 
443  Real vol = 0.;
444  const unsigned int nshp = fe->n_shape_functions();
445  for (unsigned int i=0; i<JxW.size(); i++)
446  vol += JxW[i];
447  vol /= (1.* nshp);
448  for (unsigned int i_var=0; i_var<6; i_var++)
449  for (unsigned int i=0; i<nshp; i++)
450  local_jac(i_var*nshp+i, i_var*nshp+i) =
451  vol*material_mat(i_var, i_var);
452 
453  local_f = local_jac * _local_accel;
454  }
455  else {
456 
457  for (unsigned int qp=0; qp<JxW.size(); qp++) {
458 
459  mat_inertia->derivative(p, xyz[qp], _time, material_mat);
460 
461  // now set the shape function values
462  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
463  phi_vec(i_nd) = phi[i_nd][qp];
464 
465  Bmat.reinit(n_vars, phi_vec);
466 
467  Bmat.left_multiply(mat1_n1n2, material_mat);
468 
469  vec1_n1 = mat1_n1n2 * _local_accel;
470  Bmat.vector_mult_transpose(vec2_n2, vec1_n1);
471 
472  local_f += JxW[qp] * vec2_n2;
473 
474  if (request_jacobian) {
475 
476  Bmat.right_multiply_transpose(mat2_n2n2,
477  mat1_n1n2);
478  local_jac += JxW[qp]*mat2_n2n2;
479  }
480  }
481  }
482 
483  // now transform to the global coorodinate system
484  if (_elem.dim() < 3) {
485  transform_vector_to_global_system(local_f, vec2_n2);
486  f += vec2_n2;
487 
488  if (request_jacobian) {
489  transform_matrix_to_global_system(local_jac, mat2_n2n2);
490  jac_xddot += mat2_n2n2;
491  }
492  }
493  else {
494 
495  f += local_f;
496  if (request_jacobian)
497  jac_xddot += local_jac;
498  }
499 
500  return request_jacobian;
501 }
502 
503 
504 
505 
506 void
509  const unsigned int s,
511  bool request_jacobian,
512  RealVectorX& f,
513  RealMatrixX& jac_xddot,
514  RealMatrixX& jac_xdot,
515  RealMatrixX& jac) {
516 
517  // prepare the side finite element
518  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s, false, false));
519 
520  std::vector<Real> JxW_Vn = fe->get_JxW();
521  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
522  const std::vector<std::vector<Real> >& phi = fe->get_phi();
523  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
524 
525  const unsigned int
526  n_phi = (unsigned int)phi.size(),
527  n_vars = _system.system().n_vars(),
528  n1 =6,
529  n2 =6*n_phi,
530  dim =_elem.dim();
531 
533  material_mat,
534  mat1_n1n2 = RealMatrixX::Zero(n1, n2),
535  mat2_n2n2 = RealMatrixX::Zero(n2, n2),
536  local_jac = RealMatrixX::Zero(n2, n2);
538  phi_vec = RealVectorX::Zero(n_phi),
539  vec1_n1 = RealVectorX::Zero(n1),
540  vec2_n2 = RealVectorX::Zero(n2),
541  local_f = RealVectorX::Zero(n2),
542  vel = RealVectorX::Zero(dim);
543 
544  Real
545  vn = 0.;
546 
547  // modify the JxW_Vn by multiplying the normal velocity to it
548  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
549 
550  vel_f.derivative(p, xyz[qp], _time, vel);
551  vn = 0.;
552  for (unsigned int i=0; i<dim; i++)
553  vn += vel(i)*face_normals[qp](i);
554  JxW_Vn[qp] *= vn;
555  }
556 
557 
558  std::unique_ptr<MAST::FieldFunction<RealMatrixX> >
559  mat_inertia = _property.inertia_matrix(*this);
560 
562 
564 
565  (*mat_inertia)(xyz[0], _time, material_mat);
566 
567  Real vol = 0.;
568  const unsigned int nshp = fe->n_shape_functions();
569  for (unsigned int i=0; i<JxW_Vn.size(); i++)
570  vol += JxW_Vn[i];
571  vol /= (1.* nshp);
572  for (unsigned int i_var=0; i_var<6; i_var++)
573  for (unsigned int i=0; i<nshp; i++)
574  local_jac(i_var*nshp+i, i_var*nshp+i) =
575  vol*material_mat(i_var, i_var);
576 
577  local_f = local_jac * _local_accel;
578  }
579  else {
580 
581  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
582 
583  (*mat_inertia)(xyz[qp], _time, material_mat);
584 
585  // now set the shape function values
586  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
587  phi_vec(i_nd) = phi[i_nd][qp];
588 
589  Bmat.reinit(n_vars, phi_vec);
590 
591  Bmat.left_multiply(mat1_n1n2, material_mat);
592 
593  vec1_n1 = mat1_n1n2 * _local_accel;
594  Bmat.vector_mult_transpose(vec2_n2, vec1_n1);
595 
596  local_f += JxW_Vn[qp] * vec2_n2;
597 
598  if (request_jacobian) {
599 
600  Bmat.right_multiply_transpose(mat2_n2n2,
601  mat1_n1n2);
602  local_jac += JxW_Vn[qp]*mat2_n2n2;
603  }
604  }
605  }
606 
607  // now transform to the global coorodinate system
608  if (_elem.dim() < 3) {
609  transform_vector_to_global_system(local_f, vec2_n2);
610  f += vec2_n2;
611 
612  if (request_jacobian) {
613  transform_matrix_to_global_system(local_jac, mat2_n2n2);
614  jac_xddot += mat2_n2n2;
615  }
616  }
617  else {
618 
619  f += local_f;
620  if (request_jacobian)
621  jac_xddot += local_jac;
622  }
623 }
624 
625 
626 
627 
628 bool
630 side_external_residual(bool request_jacobian,
631  RealVectorX& f,
632  RealMatrixX& jac_xdot,
633  RealMatrixX& jac,
634  std::multimap<libMesh::boundary_id_type, MAST::BoundaryConditionBase*>& bc) {
635 
636  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>> loads;
638 
639  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>>::const_iterator
640  it = loads.begin(),
641  end = loads.end();
642 
643  for ( ; it != end; it++) {
644 
645  std::vector<MAST::BoundaryConditionBase*>::const_iterator
646  bc_it = it->second.begin(),
647  bc_end = it->second.end();
648 
649  for ( ; bc_it != bc_end; bc_it++) {
650 
651  // apply all the types of loading
652  switch ((*bc_it)->type()) {
653 
655  surface_pressure_residual(request_jacobian,
656  f, jac,
657  it->first,
658  **bc_it);
659  break;
660 
661 
663  surface_traction_residual(request_jacobian,
664  f, jac,
665  it->first,
666  **bc_it);
667  break;
668 
669 
672  f, jac,
673  it->first,
674  **bc_it);
675  break;
676 
677 
678  case MAST::PISTON_THEORY:
679  piston_theory_residual(request_jacobian,
680  f,
681  jac_xdot,
682  jac,
683  it->first,
684  **bc_it);
685  break;
686 
687 
689  case MAST::DIRICHLET:
690  // nothing to be done here
691  break;
692 
693  default:
694  // not implemented yet
695  libmesh_error();
696  break;
697  }
698  }
699  }
700  return request_jacobian;
701 }
702 
703 bool
706 (bool request_jacobian,
707  RealVectorX& f,
708  RealMatrixX& jac_xdot,
709  RealMatrixX& jac,
710  std::multimap<libMesh::boundary_id_type, MAST::BoundaryConditionBase*>& bc) {
711 
712  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>> loads;
714 
715  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>>::const_iterator
716  it = loads.begin(),
717  end = loads.end();
718 
719  for ( ; it != end; it++) {
720 
721  std::vector<MAST::BoundaryConditionBase*>::const_iterator
722  bc_it = it->second.begin(),
723  bc_end = it->second.end();
724 
725  for ( ; bc_it != bc_end; bc_it++) {
726 
727  // apply all the types of loading
728  switch ((*bc_it)->type()) {
729 
731  case MAST::DIRICHLET:
732  // nothing to be done here
733  break;
734 
736  case MAST::PISTON_THEORY:
737  default:
738  // not implemented yet
739  libmesh_error();
740  break;
741  }
742  }
743  }
744  return request_jacobian;
745 }
746 
747 
748 
749 bool
752 (bool request_jacobian,
753  ComplexVectorX& f,
754  ComplexMatrixX& jac,
755  std::multimap<libMesh::boundary_id_type, MAST::BoundaryConditionBase*>& bc) {
756 
757  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>> loads;
759 
760  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>>::const_iterator
761  it = loads.begin(),
762  end = loads.end();
763 
764  for ( ; it != end; it++) {
765 
766  std::vector<MAST::BoundaryConditionBase*>::const_iterator
767  bc_it = it->second.begin(),
768  bc_end = it->second.end();
769 
770  for ( ; bc_it != bc_end; bc_it++) {
771 
772  // apply all the types of loading
773  switch ((*bc_it)->type()) {
774 
776 
778  (request_jacobian,
779  f, jac,
780  it->first,
781  **bc_it);
782  break;
783 
784 
786  case MAST::DIRICHLET:
787  // nothing to be done here
788  break;
789 
790  default:
791  // not implemented yet
792  libmesh_error();
793  break;
794  }
795  }
796  }
797  return request_jacobian;
798 }
799 
800 
801 
802 
803 bool
805 volume_external_residual (bool request_jacobian,
806  RealVectorX& f,
807  RealMatrixX& jac_xdot,
808  RealMatrixX& jac,
809  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>& bc) {
810 
811  // iterate over the boundary ids given in the provided force map
812  std::pair<std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator,
813  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator> it;
814 
815  // for each boundary id, check if any of the sides on the element
816  // has the associated boundary
817 
818  libMesh::subdomain_id_type sid = _elem.get_reference_elem().subdomain_id();
819  // find the loads on this boundary and evaluate the f and jac
820  it =bc.equal_range(sid);
821 
822  for ( ; it.first != it.second; it.first++) {
823  // apply all the types of loading
824  switch (it.first->second->type()) {
825 
827  surface_pressure_residual(request_jacobian,
828  f, jac,
829  *it.first->second);
830  break;
831 
832  case MAST::PISTON_THEORY:
833  piston_theory_residual(request_jacobian,
834  f,
835  jac_xdot,
836  jac,
837  *it.first->second);
838  break;
839 
840  case MAST::TEMPERATURE:
841  thermal_residual(request_jacobian,
842  f, jac,
843  *it.first->second);
844  break;
845 
846 
847  default:
848  // not implemented yet
849  libmesh_error();
850  break;
851  }
852  }
853 
854  return request_jacobian;
855 }
856 
857 
858 
859 
860 bool
862 linearized_volume_external_residual (bool request_jacobian,
863  RealVectorX& f,
864  RealMatrixX& jac_xdot,
865  RealMatrixX& jac,
866  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>& bc) {
867 
868  // iterate over the boundary ids given in the provided force map
869  std::pair<std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator,
870  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator> it;
871 
872  // for each boundary id, check if any of the sides on the element
873  // has the associated boundary
874 
875  libMesh::subdomain_id_type sid = _elem.get_reference_elem().subdomain_id();
876  // find the loads on this boundary and evaluate the f and jac
877  it =bc.equal_range(sid);
878 
879  for ( ; it.first != it.second; it.first++) {
880  // apply all the types of loading
881  switch (it.first->second->type()) {
882 
884  linearized_surface_pressure_residual(request_jacobian,
885  f, jac,
886  *it.first->second);
887  break;
888 
889  case MAST::TEMPERATURE: {
890 
891  const unsigned int
892  n = (unsigned int) f.size();
893 
895  local_f = RealVectorX::Zero(n);
897  mat = RealMatrixX::Zero(n, n);
898 
899  // this accounts for the perturbation of displacement.
900  // Perturbation in temperature is not currently handled
901  thermal_residual(true, local_f, mat, *it.first->second);
902  f += mat*_delta_sol;
903  }
904  break;
905 
906  case MAST::PISTON_THEORY:
907  default:
908  // not implemented yet
909  libmesh_error();
910  break;
911  }
912  }
913 
914  return request_jacobian;
915 }
916 
917 
918 
919 
920 bool
923 (bool request_jacobian,
924  ComplexVectorX& f,
925  ComplexMatrixX& jac,
926  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>& bc) {
927 
928  // iterate over the boundary ids given in the provided force map
929  std::pair<std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator,
930  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator> it;
931 
932  // for each boundary id, check if any of the sides on the element
933  // has the associated boundary
934 
935  libMesh::subdomain_id_type sid = _elem.get_reference_elem().subdomain_id();
936  // find the loads on this boundary and evaluate the f and jac
937  it =bc.equal_range(sid);
938 
939  for ( ; it.first != it.second; it.first++) {
940  // apply all the types of loading
941  switch (it.first->second->type()) {
942 
944 
946  (request_jacobian,
947  f, jac,
948  *it.first->second);
949  break;
950 
951  case MAST::TEMPERATURE:
952 
953  case MAST::PISTON_THEORY:
954  default:
955  // not implemented yet
956  libmesh_error();
957  break;
958  }
959  }
960 
961  return request_jacobian;
962 }
963 
964 
965 
966 
967 
968 
969 
970 bool
973  bool request_jacobian,
974  RealVectorX& f,
975  RealMatrixX& jac_xdot,
976  RealMatrixX& jac,
977  std::multimap<libMesh::boundary_id_type, MAST::BoundaryConditionBase*>& bc) {
978 
979  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>> loads;
981 
982  std::map<unsigned int, std::vector<MAST::BoundaryConditionBase*>>::const_iterator
983  it = loads.begin(),
984  end = loads.end();
985 
986  for ( ; it != end; it++) {
987 
988  std::vector<MAST::BoundaryConditionBase*>::const_iterator
989  bc_it = it->second.begin(),
990  bc_end = it->second.end();
991 
992  for ( ; bc_it != bc_end; bc_it++) {
993 
994  // apply all the types of loading
995  switch ((*bc_it)->type()) {
996 
999  request_jacobian,
1000  f, jac,
1001  it->first,
1002  **bc_it);
1003  break;
1004 
1005 
1008  request_jacobian,
1009  f, jac,
1010  it->first,
1011  **bc_it);
1012  break;
1013 
1014 
1017  request_jacobian,
1018  f, jac,
1019  it->first,
1020  **bc_it);
1021  break;
1022 
1023 
1024  case MAST::PISTON_THEORY:
1026  request_jacobian,
1027  f,
1028  jac_xdot,
1029  jac,
1030  it->first,
1031  **bc_it);
1032  break;
1033 
1034 
1036  case MAST::DIRICHLET:
1037  // nothing to be done here
1038  break;
1039 
1040  default:
1041  // not implemented yet
1042  libmesh_error();
1043  break;
1044  }
1045  }
1046  }
1047  return request_jacobian;
1048 }
1049 
1050 
1051 
1052 bool
1055  bool request_jacobian,
1056  RealVectorX& f,
1057  RealMatrixX& jac_xdot,
1058  RealMatrixX& jac,
1059  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>& bc) {
1060 
1061  // iterate over the boundary ids given in the provided force map
1062  std::pair<std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator,
1063  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator> it;
1064 
1065  // for each boundary id, check if any of the sides on the element
1066  // has the associated boundary
1067 
1068  libMesh::subdomain_id_type sid = _elem.get_reference_elem().subdomain_id();
1069  // find the loads on this boundary and evaluate the f and jac
1070  it =bc.equal_range(sid);
1071 
1072  for ( ; it.first != it.second; it.first++) {
1073  // apply all the types of loading
1074  switch (it.first->second->type()) {
1075 
1078  request_jacobian,
1079  f, jac,
1080  *it.first->second);
1081  break;
1082 
1083  case MAST::PISTON_THEORY:
1085  request_jacobian,
1086  f,
1087  jac_xdot,
1088  jac,
1089  *it.first->second);
1090  break;
1091 
1092  case MAST::TEMPERATURE:
1094  request_jacobian,
1095  f, jac,
1096  *it.first->second);
1097  break;
1098 
1099  default:
1100  // not implemented yet
1101  libmesh_error();
1102  break;
1103  }
1104  }
1105 
1106  return request_jacobian;
1107 }
1108 
1109 
1110 
1111 void
1114  const unsigned int s,
1115  const MAST::FieldFunction<RealVectorX>& vel_f,
1116  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>& bc,
1117  bool request_jacobian,
1118  RealVectorX& f,
1119  RealMatrixX& jac) {
1120 
1121  // iterate over the boundary ids given in the provided force map
1122  std::pair<std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator,
1123  std::multimap<libMesh::subdomain_id_type, MAST::BoundaryConditionBase*>::const_iterator> it;
1124 
1125  // for each boundary id, check if any of the sides on the element
1126  // has the associated boundary
1127 
1128  libMesh::subdomain_id_type sid = _elem.get_reference_elem().subdomain_id();
1129  // find the loads on this boundary and evaluate the f and jac
1130  it =bc.equal_range(sid);
1131 
1132  for ( ; it.first != it.second; it.first++) {
1133  // apply all the types of loading
1134  switch (it.first->second->type()) {
1135 
1138  s,
1139  vel_f,
1140  *it.first->second,
1141  request_jacobian,
1142  f,
1143  jac);
1144  break;
1145 
1146  case MAST::TEMPERATURE:
1148  s,
1149  vel_f,
1150  *it.first->second,
1151  request_jacobian,
1152  f,
1153  jac);
1154  break;
1155 
1156  case MAST::PISTON_THEORY:
1157  default:
1158  // not implemented yet
1159  libmesh_error();
1160  break;
1161  }
1162  }
1163 }
1164 
1165 
1166 
1167 bool
1169 surface_pressure_residual(bool request_jacobian,
1170  RealVectorX &f,
1171  RealMatrixX &jac,
1173 
1174  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
1175  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1176 
1177  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
1178 
1179  const std::vector<Real> &JxW = fe->get_JxW();
1180  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1181  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1182  const unsigned int
1183  n_phi = (unsigned int)phi.size(),
1184  n1 =3,
1185  n2 =6*n_phi;
1186 
1187 
1188  // normal for face integration
1189  libMesh::Point normal;
1190  // direction of pressure assumed to be normal (along local z-axis)
1191  // to the element face for 2D and along local y-axis for 1D element.
1192  normal(_elem.dim()) = -1.;
1193 
1194 
1195  // get the function from this boundary condition
1197  bc.get<MAST::FieldFunction<Real> >("pressure");
1198 
1199  Real press;
1200  FEMOperatorMatrix Bmat;
1201 
1202  RealVectorX
1203  phi_vec = RealVectorX::Zero(n_phi),
1204  force = RealVectorX::Zero(2*n1),
1205  local_f = RealVectorX::Zero(n2),
1206  vec_n2 = RealVectorX::Zero(n2);
1207 
1208 
1209  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1210 
1211  // now set the shape function values
1212  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1213  phi_vec(i_nd) = phi[i_nd][qp];
1214 
1215  Bmat.reinit(2*n1, phi_vec);
1216 
1217  // get pressure value
1218  func(qpoint[qp], _time, press);
1219 
1220  // calculate force
1221  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1222  force(i_dim) = press * normal(i_dim);
1223 
1224  Bmat.vector_mult_transpose(vec_n2, force);
1225 
1226  local_f += JxW[qp] * vec_n2;
1227  }
1228 
1229  // now transform to the global system and add
1230  transform_vector_to_global_system(local_f, vec_n2);
1231  f -= vec_n2;
1232 
1233 
1234  return (request_jacobian);
1235 }
1236 
1237 
1238 
1239 
1240 
1241 
1242 bool
1245  bool request_jacobian,
1246  RealVectorX &f,
1247  RealMatrixX &jac,
1249 
1250  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
1251  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1252 
1253  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
1254 
1255  const std::vector<Real> &JxW = fe->get_JxW();
1256  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1257  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1258  const unsigned int
1259  n_phi = (unsigned int)phi.size(),
1260  n1 =3,
1261  n2 =6*n_phi;
1262 
1263 
1264  // normal for face integration
1265  libMesh::Point normal;
1266  // direction of pressure assumed to be normal (along local z-axis)
1267  // to the element face for 2D and along local y-axis for 1D element.
1268  normal(_elem.dim()) = -1.;
1269 
1270 
1271  // get the function from this boundary condition
1273  bc.get<MAST::FieldFunction<Real> >("pressure");
1274 
1275  Real press;
1276  FEMOperatorMatrix Bmat;
1277 
1278  RealVectorX
1279  phi_vec = RealVectorX::Zero(n_phi),
1280  force = RealVectorX::Zero(2*n1),
1281  local_f = RealVectorX::Zero(n2),
1282  vec_n2 = RealVectorX::Zero(n2);
1283 
1284 
1285  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1286 
1287  // now set the shape function values
1288  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1289  phi_vec(i_nd) = phi[i_nd][qp];
1290 
1291  Bmat.reinit(2*n1, phi_vec);
1292 
1293  // get pressure value
1294  func.derivative(p, qpoint[qp], _time, press);
1295 
1296  // calculate force
1297  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1298  force(i_dim) = press * normal(i_dim);
1299 
1300  Bmat.vector_mult_transpose(vec_n2, force);
1301 
1302  local_f += JxW[qp] * vec_n2;
1303  }
1304 
1305  // now transform to the global system and add
1306  transform_vector_to_global_system(local_f, vec_n2);
1307  f -= vec_n2;
1308 
1309 
1310  return (request_jacobian);
1311 }
1312 
1313 
1314 
1315 
1316 void
1319  const unsigned int s,
1320  const MAST::FieldFunction<RealVectorX>& vel_f,
1322  bool request_jacobian,
1323  RealVectorX& f,
1324  RealMatrixX& jac) {
1325 
1326  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements
1327  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1328 
1329  // prepare the side finite element
1330  std::unique_ptr<MAST::FEBase> fe(_elem.init_side_fe(s, false, false));
1331 
1332  std::vector<Real> JxW_Vn = fe->get_JxW();
1333  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
1334  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
1335  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1336 
1337  const unsigned int
1338  n_phi = (unsigned int)phi.size(),
1339  n1 = 3,
1340  n2 = 6*n_phi,
1341  dim = _elem.dim();
1342 
1343 
1344  // normal for face integration
1345  libMesh::Point normal;
1346  // direction of pressure assumed to be normal (along local z-axis)
1347  // to the element face for 2D and along local y-axis for 1D element.
1348  normal(_elem.dim()) = -1.;
1349 
1350 
1351  // get the function from this boundary condition
1353  bc.get<MAST::FieldFunction<Real> >("pressure");
1354 
1355  Real press;
1356  FEMOperatorMatrix Bmat;
1357 
1358  RealVectorX
1359  phi_vec = RealVectorX::Zero(n_phi),
1360  force = RealVectorX::Zero(2*n1),
1361  local_f = RealVectorX::Zero(n2),
1362  vec_n2 = RealVectorX::Zero(n2),
1363  vel = RealVectorX::Zero(dim);
1364 
1365  Real
1366  vn = 0.;
1367 
1368 
1369  // modify the JxW_Vn by multiplying the normal velocity to it
1370  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
1371 
1372  vel_f.derivative(p, xyz[qp], _time, vel);
1373  vn = 0.;
1374  for (unsigned int i=0; i<dim; i++)
1375  vn += vel(i)*face_normals[qp](i);
1376  JxW_Vn[qp] *= vn;
1377  }
1378 
1379 
1380  for (unsigned int qp=0; qp<xyz.size(); qp++) {
1381 
1382  // now set the shape function values
1383  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1384  phi_vec(i_nd) = phi[i_nd][qp];
1385 
1386  Bmat.reinit(2*n1, phi_vec);
1387 
1388  // get pressure value
1389  func(xyz[qp], _time, press);
1390 
1391  // calculate force
1392  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1393  force(i_dim) = press * normal(i_dim);
1394 
1395  Bmat.vector_mult_transpose(vec_n2, force);
1396 
1397  local_f += JxW_Vn[qp] * vec_n2;
1398  }
1399 
1400  // now transform to the global system and add
1401  transform_vector_to_global_system(local_f, vec_n2);
1402  f -= vec_n2;
1403 }
1404 
1405 
1406 
1407 
1408 
1409 bool
1412  RealVectorX &f,
1413  RealMatrixX &jac,
1415 
1416  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements.
1417  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1418  libmesh_assert_equal_to(bc.type(), MAST::SURFACE_PRESSURE);
1419 
1421  press_fn = bc.get<MAST::FieldFunction<Real> >("pressure");
1422 
1423  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
1424 
1425  const std::vector<Real> &JxW = fe->get_JxW();
1426  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1427  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1428  const unsigned int
1429  n_phi = (unsigned int)phi.size(),
1430  n1 = 3,
1431  n2 = 6*n_phi;
1432 
1433  // normal for face integration
1434  libMesh::Point normal;
1435  // direction of pressure assumed to be normal (along local z-axis)
1436  // to the element face for 2D and along local y-axis for 1D element.
1437  normal(_elem.dim()) = -1.;
1438 
1439  RealVectorX
1440  phi_vec = RealVectorX::Zero(n_phi),
1441  w = RealVectorX::Zero(3),
1442  dn_rot = RealVectorX::Zero(3),
1443  force = RealVectorX::Zero(2*n1),
1444  local_f = RealVectorX::Zero(n2),
1445  vec_n2 = RealVectorX::Zero(n2);
1446 
1447  FEMOperatorMatrix Bmat;
1448  Real
1449  press = 0.,
1450  dpress = 0.;
1451 
1452 
1453  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1454 
1455  // now set the shape function values
1456  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1457  phi_vec(i_nd) = phi[i_nd][qp];
1458 
1459  Bmat.reinit(2*n1, phi_vec);
1460 
1461  // get pressure and deformation information
1462  press_fn (qpoint[qp], _time, press);
1463  press_fn.perturbation(qpoint[qp], _time, dpress);
1464 
1465  // calculate force
1466  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1467  force(i_dim) = ( press * dn_rot(i_dim) + // steady pressure
1468  dpress * normal(i_dim)); // unsteady pressure
1469 
1470  Bmat.vector_mult_transpose(vec_n2, force);
1471 
1472  local_f += JxW[qp] * vec_n2;
1473  }
1474 
1475  // now transform to the global system and add
1476  transform_vector_to_global_system(local_f, vec_n2);
1477  f -= vec_n2;
1478 
1479 
1480  return (request_jacobian);
1481 }
1482 
1483 
1484 
1485 
1486 bool
1489 (bool request_jacobian,
1490  ComplexVectorX &f,
1491  ComplexMatrixX &jac,
1493 
1494  libmesh_assert(_elem.dim() < 3); // only applicable for lower dimensional elements.
1495  libmesh_assert(!follower_forces); // not implemented yet for follower forces
1496  libmesh_assert_equal_to(bc.type(), MAST::SURFACE_PRESSURE);
1497 
1499  press_fn = bc.get<MAST::FieldFunction<Real> >("pressure");
1501  dpress_fn = bc.get<MAST::FieldFunction<Complex> >("frequency_domain_pressure");
1502 
1503  std::unique_ptr<MAST::FEBase> fe(_elem.init_fe(false, false));
1504 
1505  const std::vector<Real> &JxW = fe->get_JxW();
1506  const std::vector<libMesh::Point>& qpoint = fe->get_xyz();
1507  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1508  const unsigned int
1509  n_phi = (unsigned int)phi.size(),
1510  n1 = 3,
1511  n2 = 6*n_phi;
1512 
1513  // normal for face integration
1514  libMesh::Point normal;
1515  // direction of pressure assumed to be normal (along local z-axis)
1516  // to the element face for 2D and along local y-axis for 1D element.
1517  normal(_elem.dim()) = -1.;
1518 
1519  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
1521  w = ComplexVectorX::Zero(3),
1522  dn_rot = ComplexVectorX::Zero(3),
1523  force = ComplexVectorX::Zero(2*n1),
1524  local_f = ComplexVectorX::Zero(n2),
1525  vec_n2 = ComplexVectorX::Zero(n2);
1526 
1527  FEMOperatorMatrix Bmat;
1528  Real press = 0.;
1529  Complex dpress = 0.;
1530 
1531 
1532  for (unsigned int qp=0; qp<qpoint.size(); qp++) {
1533 
1534  // now set the shape function values
1535  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
1536  phi_vec(i_nd) = phi[i_nd][qp];
1537 
1538  Bmat.reinit(2*n1, phi_vec);
1539 
1540  // get pressure and deformation information
1541  press_fn (qpoint[qp], _time, press);
1542  dpress_fn(qpoint[qp], _time, dpress);
1543  //dn_rot_fn.freq_domain_motion(pt, normal, w, dn_rot);
1544 
1545  // calculate force
1546  for (unsigned int i_dim=0; i_dim<n1; i_dim++)
1547  force(i_dim) = ( press * dn_rot(i_dim) + // steady pressure
1548  dpress * normal(i_dim)); // unsteady pressure
1549 
1550  Bmat.vector_mult_transpose(vec_n2, force);
1551 
1552  local_f += JxW[qp] * vec_n2;
1553  }
1554 
1555  // now transform to the global system and add
1556  transform_vector_to_global_system(local_f, vec_n2);
1557  f -= vec_n2;
1558 
1559 
1560  return (request_jacobian);
1561 }
1562 
1563 
1564 
1565 template <typename ValType>
1566 void
1568 transform_matrix_to_global_system(const ValType& local_mat,
1569  ValType& global_mat) const {
1570 
1571  // make sure this is only called for 1D and 2D elements
1572  libmesh_assert_less(_elem.dim(), 3);
1573  libmesh_assert_equal_to( local_mat.rows(), local_mat.cols());
1574  libmesh_assert_equal_to(global_mat.rows(), global_mat.cols());
1575  libmesh_assert_equal_to( local_mat.rows(), global_mat.rows());
1576 
1577  if (_elem.use_local_elem()) {
1578 
1579  // assuming same shape function for all 6 dofs.
1580  const unsigned int n_dofs = local_mat.rows()/6;
1581 
1582  ValType mat(6*n_dofs, 6*n_dofs);
1583 
1584  mat.setZero();
1585  global_mat.setZero();
1586 
1587  const RealMatrixX& Tmat = _elem.T_matrix();
1588  // now initialize the global T matrix
1589  for (unsigned int i=0; i<n_dofs; i++)
1590  for (unsigned int j=0; j<3; j++)
1591  for (unsigned int k=0; k<3; k++) {
1592  mat(j*n_dofs+i, k*n_dofs+i) = Tmat(j,k); // for u,v,w
1593  mat((j+3)*n_dofs+i, (k+3)*n_dofs+i) = Tmat(j,k); // for tx,ty,tz
1594  }
1595 
1596  // right multiply with T^tr, and left multiply with T.
1597  global_mat = mat * local_mat * mat.transpose();
1598  }
1599  else
1600  global_mat = local_mat;
1601 }
1602 
1603 
1604 
1605 template <typename ValType>
1606 void
1608 transform_vector_to_local_system(const ValType& global_vec,
1609  ValType& local_vec) const {
1610 
1611  // make sure this is only called for 1D and 2D elements
1612  libmesh_assert_less(_elem.dim(), 3);
1613  libmesh_assert_equal_to( local_vec.size(), global_vec.size());
1614 
1615  if (_elem.use_local_elem()) {
1616 
1617  // assuming same shape function for all 6 dofs.
1618  const unsigned int n_dofs = global_vec.size()/6;
1619  RealMatrixX mat = RealMatrixX::Zero(6*n_dofs, 6*n_dofs);
1620 
1621  local_vec.setZero();
1622 
1623  const RealMatrixX& Tmat = _elem.T_matrix();
1624  // now initialize the global T matrix
1625  for (unsigned int i=0; i<n_dofs; i++)
1626  for (unsigned int j=0; j<3; j++)
1627  for (unsigned int k=0; k<3; k++) {
1628  mat(j*n_dofs+i, k*n_dofs+i) = Tmat(j,k); // for u,v,w
1629  mat((j+3)*n_dofs+i, (k+3)*n_dofs+i) = Tmat(j,k); // for tx,ty,tz
1630  }
1631 
1632  // left multiply with T^tr
1633  local_vec = mat.transpose() * global_vec;
1634  }
1635  else
1636  local_vec = global_vec;
1637 }
1638 
1639 
1640 
1641 template <typename ValType>
1642 void
1644 transform_vector_to_global_system(const ValType& local_vec,
1645  ValType& global_vec) const {
1646 
1647  // make sure this is only called for 1D and 2D elements
1648  libmesh_assert_less(_elem.dim(), 3);
1649  libmesh_assert_equal_to( local_vec.size(), global_vec.size());
1650 
1651  if (_elem.use_local_elem()) {
1652 
1653  // assuming same shape function for all 6 dofs.
1654  const unsigned int n_dofs = local_vec.size()/6;
1655 
1656  RealMatrixX mat = RealMatrixX::Zero(6*n_dofs, 6*n_dofs);
1657 
1658  global_vec.setZero();
1659 
1660  const RealMatrixX& Tmat = _elem.T_matrix();
1661  // now initialize the global T matrix
1662  for (unsigned int i=0; i<n_dofs; i++)
1663  for (unsigned int j=0; j<3; j++)
1664  for (unsigned int k=0; k<3; k++) {
1665  mat(j*n_dofs+i, k*n_dofs+i) = Tmat(j,k); // for u,v,w
1666  mat((j+3)*n_dofs+i, (k+3)*n_dofs+i) = Tmat(j,k); // for tx,ty,tz
1667  }
1668 
1669  // left multiply with T
1670  global_vec = mat* local_vec;
1671  }
1672  else
1673  global_vec = local_vec;
1674 }
1675 
1676 
1677 
1678 
1679 
1680 std::unique_ptr<MAST::StructuralElementBase>
1682  const MAST::GeomElem& elem,
1683  const MAST::ElementPropertyCardBase& p) {
1684 
1685  std::unique_ptr<MAST::StructuralElementBase> e;
1686 
1687  switch (elem.dim()) {
1688  case 1:
1689  e.reset(new MAST::StructuralElement1D(sys, elem, p));
1690  break;
1691 
1692  case 2:
1693  e.reset(new MAST::StructuralElement2D(sys, elem, p));
1694  break;
1695 
1696  case 3:
1697  e.reset(new MAST::StructuralElement3D(sys, elem, p));
1698  break;
1699 
1700  default:
1701  libmesh_error();
1702  break;
1703  }
1704 
1705  return e;
1706 }
1707 
1708 
1709 
1710 Real
1712  const Real vel_U,
1713  const Real gamma,
1714  const Real mach) {
1715 
1716 
1717  Real cp = 0.0;
1718  switch (order)
1719  {
1720  case 3:
1721  cp += (gamma+1.0)/12.0*mach*mach*pow(vel_U,3);
1722  case 2:
1723  cp += (gamma+1.0)/4.0*mach*pow(vel_U,2);
1724  case 1: {
1725  cp += vel_U;
1726  cp *= 2.0/pow(mach*mach-1.,0.5);
1727  }
1728  break;
1729 
1730  default:
1731  libmesh_error_msg("Invalid Piston Theory Order: " << order);
1732  break;
1733  }
1734 
1735  return cp;
1736 }
1737 
1738 
1739 
1740 
1741 Real
1743  const Real vel_U,
1744  const Real gamma,
1745  const Real mach) {
1746 
1747 
1748  Real dcp_dvn = 0.0;
1749  switch (order)
1750  {
1751  case 3:
1752  dcp_dvn += (gamma+1.0)/4.0*mach*mach*pow(vel_U,2);
1753  case 2:
1754  dcp_dvn += (gamma+1.0)/2.0*mach*vel_U;
1755  case 1: {
1756  dcp_dvn += 1.;
1757  dcp_dvn *= 2.0/pow(mach*mach-1.,.5);
1758  }
1759  break;
1760 
1761  default:
1762  libmesh_error_msg("Invalid Piston Theory Order: " << order);
1763  break;
1764  }
1765 
1766  return dcp_dvn;
1767 }
1768 
1769 
1770 
1771 
1772 // template instantiations
1773 template
1774 void
1775 MAST::StructuralElementBase::transform_matrix_to_global_system<RealMatrixX>
1776 (const RealMatrixX& local_mat,
1777  RealMatrixX& global_mat) const;
1778 
1779 
1780 template
1781 void
1782 MAST::StructuralElementBase::transform_vector_to_local_system<RealVectorX>
1783 (const RealVectorX& global_vec,
1784  RealVectorX& local_vec) const;
1785 
1786 
1787 template
1788 void
1789 MAST::StructuralElementBase::transform_vector_to_global_system<RealVectorX>
1790 (const RealVectorX& local_vec,
1791  RealVectorX& global_vec) const;
1792 
1793 
1794 template
1795 void
1796 MAST::StructuralElementBase::transform_matrix_to_global_system<ComplexMatrixX>
1797 (const ComplexMatrixX& local_mat,
1798  ComplexMatrixX& global_mat) const;
1799 
1800 
1801 template
1802 void
1803 MAST::StructuralElementBase::transform_vector_to_local_system<ComplexVectorX>
1804 (const ComplexVectorX& global_vec,
1805  ComplexVectorX& local_vec) const;
1806 
1807 
1808 template
1809 void
1810 MAST::StructuralElementBase::transform_vector_to_global_system<ComplexVectorX>
1811 (const ComplexVectorX& local_vec,
1812  ComplexVectorX& global_vec) const;
1813 
1814 
1815 
RealVectorX _local_accel_sens
local acceleration sensitivity
bool linearized_frequency_domain_volume_external_residual(bool request_jacobian, ComplexVectorX &f, ComplexMatrixX &jac, std::multimap< libMesh::subdomain_id_type, MAST::BoundaryConditionBase *> &bc)
Calculates the frequency domain volume external force contribution to system residual.
MAST::NonlinearSystem & system()
virtual void external_side_loads_for_quadrature_elem(std::multimap< libMesh::boundary_id_type, MAST::BoundaryConditionBase *> &bc, std::map< unsigned int, std::vector< MAST::BoundaryConditionBase *>> &loads) const
From the given list of boundary loads, this identifies the sides of the quadrature element and the lo...
Definition: geom_elem.cpp:183
virtual void surface_pressure_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, MAST::BoundaryConditionBase &bc, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
Calculates the force vector and Jacobian due to surface pressure applied on the entire element domain...
virtual bool thermal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to thermal stresses.
virtual bool linearized_frequency_domain_surface_pressure_residual(bool request_jacobian, ComplexVectorX &f, ComplexMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to small perturbation surface pressure.
virtual bool thermal_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)=0
Calculates the sensitivity of force vector and Jacobian due to thermal stresses.
const MAST::GeomElem & _elem
geometric element for which the computations are performed
Definition: elem_base.h:205
RealVectorX _local_delta_vel_sens
local perturbed velocity sensitivity
virtual bool inertial_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xddot, RealMatrixX &jac_xdot, RealMatrixX &jac)
sensitivity of the inertial force contribution to system residual
virtual const libMesh::Elem & get_reference_elem() const
Definition: geom_elem.cpp:54
virtual bool piston_theory_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to piston-theory based surface pressure on the entire el...
virtual bool linearized_inertial_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xddot, RealMatrixX &jac_xdot, RealMatrixX &jac)
inertial force contribution to system residual of linerized problem
virtual void thermal_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, MAST::BoundaryConditionBase &bc, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)=0
Calculates the sensitivity of force vector and Jacobian due to thermal stresses.
void reinit(unsigned int n_interpolated_vars, unsigned int n_discrete_vars, unsigned int n_discrete_dofs_per_var)
this initializes the operator for number of rows and variables, assuming that all variables has the s...
virtual void set_velocity(const RealVectorX &vec, bool if_sens=false)
stores vec as velocity for element level calculations, or its sensitivity if if_sens is true...
Definition: elem_base.cpp:98
virtual void set_perturbed_acceleration(const RealVectorX &vec, bool if_sens=false)
stores vec as perturbed acceleration for element level calculations, or its sensitivity if if_sens is...
virtual void set_acceleration(const RealVectorX &vec, bool if_sens=false)
stores vec as acceleration for element level calculations, or its sensitivity if if_sens is true...
Matrix< Complex, Dynamic, 1 > ComplexVectorX
RealVectorX _local_delta_sol_sens
local perturbed solution sensitivity
virtual bool inertial_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xddot, RealMatrixX &jac_xdot, RealMatrixX &jac)
inertial force contribution to system residual
RealVectorX _local_delta_accel_sens
local perturbed acceleration sensitivity
RealVectorX _local_accel
local acceleration
bool linearized_volume_external_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::subdomain_id_type, MAST::BoundaryConditionBase *> &bc)
volume external force contribution to system residual.
virtual void derivative(const MAST::FunctionBase &f, ValType &v) const
calculates the value of the function derivative and returns it in v.
libMesh::Real Real
RealVectorX _local_sol_sens
local solution sensitivity
MAST::SystemInitialization & _system
SystemInitialization object associated with this element.
Definition: elem_base.h:200
libMesh::Complex Complex
void transform_vector_to_local_system(const ValType &global_vec, ValType &local_vec) const
Real piston_theory_cp(const unsigned int order, const Real vel_U, const Real gamma, const Real mach)
RealVectorX _local_delta_vel
local perturbed velocity
virtual std::unique_ptr< MAST::FieldFunction< RealMatrixX > > inertia_matrix(const MAST::ElementBase &e) const =0
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
virtual void set_velocity(const RealVectorX &vec, bool if_sens=false)
stores vec as velocity for element level calculations, or its sensitivity if if_sens is true...
RealVectorX _delta_sol
local solution used for linearized analysis
Definition: elem_base.h:249
RealVectorX _local_delta_sol
local perturbed solution
bool follower_forces
flag for follower forces
Matrix< Real, Dynamic, Dynamic > RealMatrixX
Matrix< Complex, Dynamic, Dynamic > ComplexMatrixX
Real piston_theory_dcp_dv(const unsigned int order, const Real vel_U, const Real gamma, const Real mach)
virtual bool surface_traction_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the sensitivity of element vector and matrix quantities for surface traction boundary cond...
bool side_external_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::boundary_id_type, MAST::BoundaryConditionBase *> &bc)
side external force contribution to system residual.
virtual void set_solution(const RealVectorX &vec, bool if_sens=false)
stores vec as solution for element level calculations, or its sensitivity if if_sens is true...
void vector_mult_transpose(T &res, const T &v) const
res = v^T * [this]
Matrix< Real, Dynamic, 1 > RealVectorX
void right_multiply_transpose(T &r, const T &m) const
[R] = [this]^T * [M]
bool volume_external_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::subdomain_id_type, MAST::BoundaryConditionBase *> &bc)
volume external force contribution to system residual.
std::unique_ptr< MAST::StructuralElementBase > build_structural_element(MAST::SystemInitialization &sys, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
builds the structural element for the specified element type
virtual bool piston_theory_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to piston-theory based surface pressure on the entire el...
const MAST::ElementPropertyCardBase & _property
element property
virtual bool surface_traction_residual_shifted_boundary_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the sensitivity of force vector and Jacobian due to surface traction and sensitiity due to...
const RealMatrixX & T_matrix() const
O.
Definition: geom_elem.cpp:236
virtual void set_perturbed_solution(const RealVectorX &vec, bool if_sens=false)
This provides the perturbed solution (or its sensitivity if if_sens is true.) for linearized analysis...
Definition: elem_base.cpp:72
StructuralElementBase(MAST::SystemInitialization &sys, const MAST::GeomElem &elem, const MAST::ElementPropertyCardBase &p)
Constructor.
This class acts as a wrapper around libMesh::Elem for the purpose of providing a uniform interface fo...
Definition: geom_elem.h:59
virtual void set_perturbed_velocity(const RealVectorX &vec, bool if_sens=false)
stores vec as perturbed velocity for element level calculations, or its sensitivity if if_sens is tru...
virtual bool surface_traction_residual_shifted_boundary(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the sensitivity of force vector and Jacobian due to surface traction and sensitiity due to...
unsigned int dim() const
Definition: geom_elem.cpp:91
RealVectorX _local_vel_sens
local velocity sensitivity
virtual bool linearized_surface_pressure_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, MAST::BoundaryConditionBase &bc)
Calculates the force vector and Jacobian due to surface pressure applied on the entire element domain...
bool if_diagonal_mass_matrix() const
returns the type of strain to be used for this element
const ValType & get(const std::string &nm) const
returns a constant reference to the specified function
const MAST::GeomElem & elem() const
Definition: elem_base.h:109
virtual void set_perturbed_solution(const RealVectorX &vec, bool if_sens=false)
stores vec as perturbed solution for element level calculations, or its sensitivity if if_sens is tru...
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::BoundaryConditionType type() const
virtual void set_solution(const RealVectorX &vec, bool if_sens=false)
stores vec as solution for element level calculations, or its sensitivity if if_sens is true...
Definition: elem_base.cpp:60
bool side_external_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::boundary_id_type, MAST::BoundaryConditionBase *> &bc)
sensitivity of the side external force contribution to system residual
bool linearized_side_external_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::boundary_id_type, MAST::BoundaryConditionBase *> &bc)
side external force contribution to system residual.
void transform_matrix_to_global_system(const ValType &local_mat, ValType &global_mat) const
bool linearized_frequency_domain_side_external_residual(bool request_jacobian, ComplexVectorX &f, ComplexMatrixX &jac, std::multimap< libMesh::boundary_id_type, MAST::BoundaryConditionBase *> &bc)
Calculates the external force due to frequency domain side external force contribution to system resi...
void left_multiply(T &r, const T &m) const
[R] = [M] * [this]
virtual bool linearized_internal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
internal force contribution to system residual of the linearized problem
RealVectorX _local_delta_accel
local perturbed acceleration
RealVectorX _local_sol
local solution
virtual bool surface_pressure_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to surface pressure.
bool volume_external_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xdot, RealMatrixX &jac, std::multimap< libMesh::subdomain_id_type, MAST::BoundaryConditionBase *> &bc)
sensitivity of the volume external force contribution to system residual
RealVectorX _local_vel
local velocity
virtual void set_perturbed_velocity(const RealVectorX &vec, bool if_sens=false)
stores vec as perturbed velocity for element level calculations, or its sensitivity if if_sens is tru...
Definition: elem_base.cpp:110
virtual bool surface_traction_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to surface traction.
virtual bool internal_residual(bool request_jacobian, RealVectorX &f, RealMatrixX &jac)=0
internal force contribution to system residual
const Real & _time
time for which system is being assembled
Definition: elem_base.h:219
virtual void set_perturbed_acceleration(const RealVectorX &vec, bool if_sens=false)
stores vec as perturbed acceleration for element level calculations, or its sensitivity if if_sens is...
Definition: elem_base.cpp:134
virtual void inertial_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, bool request_jacobian, RealVectorX &f, RealMatrixX &jac_xddot, RealMatrixX &jac_xdot, RealMatrixX &jac)
sensitivity of the inertial force contribution to system residual
void volume_external_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, std::multimap< libMesh::subdomain_id_type, MAST::BoundaryConditionBase *> &bc, bool request_jacobian, RealVectorX &f, RealMatrixX &jac)
boundary velocity contribution of volume external force.
virtual void set_acceleration(const RealVectorX &vec, bool if_sens=false)
stores vec as acceleration for element level calculations, or its sensitivity if if_sens is true...
Definition: elem_base.cpp:122
void transform_vector_to_global_system(const ValType &local_vec, ValType &global_vec) const
bool use_local_elem() const
Vector and matrix quantities defined on one- and two-dimensional elements that are oriented in two or...
Definition: geom_elem.cpp:322
virtual bool surface_pressure_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &f, RealMatrixX &jac, const unsigned int side, MAST::BoundaryConditionBase &bc)=0
Calculates the force vector and Jacobian due to surface pressure.
This is the base class for elements that implement calculation of finite element quantities over the ...
Definition: elem_base.h:72