MAST
mindlin_bending_operator.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 
21 // MAST includes
26 #include "base/nonlinear_system.h"
27 #include "mesh/fe_base.h"
28 #include "mesh/geom_elem.h"
29 #include "base/assembly_base.h"
31 
32 
33 
36 MAST::BendingOperator2D(elem),
37 _shear_quadrature_reduction(2)
38 { }
39 
40 
41 
43 
44 
45 
46 void
49  const unsigned int qp,
50  MAST::FEMOperatorMatrix& Bmat_bend) {
51 
52  this->initialize_bending_strain_operator_for_z(fe, qp, 1., Bmat_bend);
53 }
54 
55 
56 
57 void
60  const unsigned int qp,
61  const Real z,
62  MAST::FEMOperatorMatrix& Bmat_bend) {
63 
64  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe.get_dphi();
65  const std::vector<std::vector<Real> >& phi = fe.get_phi();
66 
67  const unsigned int n_phi = (unsigned int)phi.size();
68 
69  RealVectorX phi_vec = RealVectorX::Zero(n_phi);
70  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
71  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
72 
73  phi_vec *= z;
74  Bmat_bend.set_shape_function(0, 4, phi_vec); // epsilon-x: thetay
75  phi_vec *= -1.0;
76  Bmat_bend.set_shape_function(2, 3, phi_vec); // gamma-xy : thetax
77 
78 
79  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
80  phi_vec(i_nd) = dphi[i_nd][qp](1); // dphi/dy
81 
82  phi_vec *= z;
83  Bmat_bend.set_shape_function(2, 4, phi_vec); // gamma-xy : thetay
84  //Bmat_trans.set_shape_function(1, 2, phi_vec); // gamma-yz : w
85  phi_vec *= -1.0;
86  Bmat_bend.set_shape_function(1, 3, phi_vec); // epsilon-y: thetax
87 }
88 
89 
90 
91 
92 void
95  RealVectorX& local_f,
96  RealMatrixX& local_jac) {
97 
99 
100  // make an fe and quadrature object for the requested order for integrating
101  // transverse shear
102 
103  std::unique_ptr<MAST::FEBase>
104  fe(_elem.init_fe(true, false, -_shear_quadrature_reduction));
105 
106  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe->get_dphi();
107  const std::vector<std::vector<Real> >& phi = fe->get_phi();
108  const std::vector<Real>& JxW = fe->get_JxW();
109  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
110 
111  const unsigned int
112  n_phi = (unsigned int)phi.size(),
113  n2 = 6*n_phi;
114 
116  phi_vec = RealVectorX::Zero(n_phi),
117  vec_n2 = RealVectorX::Zero(n2),
118  vec_2 = RealVectorX::Zero(2);
120  material_trans_shear_mat,
121  mat_n2n2 = RealMatrixX::Zero(n2,n2),
122  mat_2n2 = RealMatrixX::Zero(2,n2);
123 
124 
125  FEMOperatorMatrix Bmat_trans;
126  Bmat_trans.reinit(2, 6, n_phi); // only two shear stresses
127 
128  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
129  mat_stiff = property.transverse_shear_stiffness_matrix(_structural_elem);
130 
131  for (unsigned int qp=0; qp<JxW.size(); qp++) {
132 
133  (*mat_stiff)(xyz[qp],
134  _structural_elem.system().time,
135  material_trans_shear_mat);
136 
138  dphi,
139  JxW,
140  qp,
141  material_trans_shear_mat,
142  Bmat_trans,
143  phi_vec,
144  vec_n2,
145  vec_2,
146  mat_n2n2,
147  mat_2n2,
148  request_jacobian,
149  local_f,
150  local_jac);
151 
152  }
153 }
154 
155 
156 
157 
158 
159 
160 void
163  bool request_jacobian,
164  RealVectorX& local_f,
165  RealMatrixX& local_jac) {
166 
168 
169  // make an fe and quadrature object for the requested order for integrating
170  // transverse shear
171 
172  std::unique_ptr<MAST::FEBase>
173  fe(_elem.init_fe(true, false, -_shear_quadrature_reduction));
174 
175  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe->get_dphi();
176  const std::vector<std::vector<Real> >& phi = fe->get_phi();
177  const std::vector<Real>& JxW = fe->get_JxW();
178  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
179 
180  const unsigned int
181  n_phi = (unsigned int)phi.size(),
182  n2 = 6*n_phi;
183 
185  phi_vec = RealVectorX::Zero(n_phi),
186  vec_n2 = RealVectorX::Zero(n2),
187  vec_2 = RealVectorX::Zero(2);
189  material_trans_shear_mat,
190  mat_n2n2 = RealMatrixX::Zero(n2,n2),
191  mat_2n2 = RealMatrixX::Zero(2,n2);
192 
193 
194  FEMOperatorMatrix Bmat_trans;
195  Bmat_trans.reinit(2, 6, n_phi); // only two shear stresses
196 
197  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
198  mat_stiff = property.transverse_shear_stiffness_matrix(_structural_elem);
199 
200  for (unsigned int qp=0; qp<JxW.size(); qp++) {
201 
202  mat_stiff->derivative(p,
203  xyz[qp],
204  _structural_elem.system().time,
205  material_trans_shear_mat);
206 
208  dphi,
209  JxW,
210  qp,
211  material_trans_shear_mat,
212  Bmat_trans,
213  phi_vec,
214  vec_n2,
215  vec_2,
216  mat_n2n2,
217  mat_2n2,
218  request_jacobian,
219  local_f,
220  local_jac);
221  }
222 }
223 
224 
225 
226 void
230  const unsigned int s,
232  bool request_jacobian,
233  RealVectorX& local_f,
234  RealMatrixX& local_jac) {
235 
237 
238  // make an fe and quadrature object for the requested order for integrating
239  // transverse shear
240 
241  std::unique_ptr<MAST::FEBase>
243 
244  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = fe->get_dphi();
245  const std::vector<std::vector<Real> >& phi = fe->get_phi();
246  const std::vector<libMesh::Point>& xyz = fe->get_xyz();
247  const std::vector<libMesh::Point>& face_normals = fe->get_normals_for_local_coordinate();
248  std::vector<Real> JxW_Vn = fe->get_JxW();
249 
250  const unsigned int
251  n_phi = (unsigned int)phi.size(),
252  n2 = 6*n_phi,
253  dim = 2;
254 
256  phi_vec = RealVectorX::Zero(n_phi),
257  vec_n2 = RealVectorX::Zero(n2),
258  vec_2 = RealVectorX::Zero(2),
259  vel = RealVectorX::Zero(dim);
261  material_trans_shear_mat,
262  dmaterial_trans_shear_mat_dp,
263  mat_n2n2 = RealMatrixX::Zero(n2,n2),
264  mat_2n2 = RealMatrixX::Zero(2,n2);
265 
266 
267  FEMOperatorMatrix Bmat_trans;
268  Bmat_trans.reinit(2, 6, n_phi); // only two shear stresses
269 
270  std::unique_ptr<MAST::FieldFunction<RealMatrixX > >
271  mat_stiff = property.transverse_shear_stiffness_matrix(_structural_elem);
272 
273 
274  Real
275  vn = 0.;
276 
277  // modify the JxW_Vn by multiplying the normal velocity to it
278  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
279 
280  vel_f(xyz[qp], _structural_elem.system().time, vel);
281  vn = 0.;
282  for (unsigned int i=0; i<dim; i++)
283  vn += vel(i)*face_normals[qp](i);
284  JxW_Vn[qp] *= vn;
285  }
286 
287  for (unsigned int qp=0; qp<JxW_Vn.size(); qp++) {
288 
289  (*mat_stiff)(xyz[qp],
290  _structural_elem.system().time,
291  material_trans_shear_mat);
292 
294  dphi,
295  JxW_Vn,
296  qp,
297  material_trans_shear_mat,
298  Bmat_trans,
299  phi_vec,
300  vec_n2,
301  vec_2,
302  mat_n2n2,
303  mat_2n2,
304  request_jacobian,
305  local_f,
306  local_jac);
307  }
308 }
309 
310 
311 
312 void
314 _transverse_shear_operations(const std::vector<std::vector<Real> >& phi,
315  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi,
316  const std::vector<Real>& JxW,
317  const unsigned int qp,
318  const RealMatrixX& material,
319  FEMOperatorMatrix& Bmat,
320  RealVectorX& phi_vec,
321  RealVectorX& vec_n2,
322  RealVectorX& vec_2,
323  RealMatrixX& mat_n2n2,
324  RealMatrixX& mat_2n2,
325  bool request_jacobian,
326  RealVectorX& local_f,
327  RealMatrixX& local_jac) {
328 
329  // initialize the strain operator
330  for ( unsigned int i_nd=0; i_nd<phi.size(); i_nd++ )
331  phi_vec(i_nd) = dphi[i_nd][qp](0); // dphi/dx
332 
333  Bmat.set_shape_function(0, 2, phi_vec); // gamma-xz: w
334 
335  for ( unsigned int i_nd=0; i_nd<phi.size(); i_nd++ )
336  phi_vec(i_nd) = dphi[i_nd][qp](1); // dphi/dy
337 
338  Bmat.set_shape_function(1, 2, phi_vec); // gamma-yz : w
339 
340  for ( unsigned int i_nd=0; i_nd<phi.size(); i_nd++ )
341  phi_vec(i_nd) = phi[i_nd][qp]; // phi
342 
343  Bmat.set_shape_function(0, 4, phi_vec); // gamma-xz: thetay
344  phi_vec *= -1.0;
345  Bmat.set_shape_function(1, 3, phi_vec); // gamma-yz : thetax
346 
347 
348  // now add the transverse shear component
350  vec_2 = material * vec_2;
351  Bmat.vector_mult_transpose(vec_n2, vec_2);
352  local_f += JxW[qp] * vec_n2;
353 
354  if (request_jacobian) {
355 
356  // now add the transverse shear component
357  Bmat.left_multiply(mat_2n2, material);
358  Bmat.right_multiply_transpose(mat_n2n2, mat_2n2);
359  local_jac += JxW[qp] * mat_n2n2;
360  }
361 }
void vector_mult_transpose(T &res, const T &v) const
res = v^T * [this]
void initialize_bending_strain_operator_for_z(const MAST::FEBase &fe, const unsigned int qp, const Real z, MAST::FEMOperatorMatrix &Bmat_bend)
initializes the bending strain operator for the specified quadrature point and z-location.
virtual void initialize_bending_strain_operator(const MAST::FEBase &fe, const unsigned int qp, MAST::FEMOperatorMatrix &Bmat)
initialze the bending strain operator for Mindlin element, withouth the z-location.
void _transverse_shear_operations(const std::vector< std::vector< Real > > &phi, const std::vector< std::vector< libMesh::RealVectorValue > > &dphi, const std::vector< Real > &JxW, const unsigned int qp, const RealMatrixX &material, FEMOperatorMatrix &Bmat_trans, RealVectorX &phi_vec, RealVectorX &vec_n2, RealVectorX &vec_2, RealMatrixX &mat_n2n2, RealMatrixX &mat_2n2, bool request_jacobian, RealVectorX &local_f, RealMatrixX &local_jac)
MAST::NonlinearSystem & system()
Definition: elem_base.cpp:45
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 calculate_transverse_shear_residual(bool request_jacobian, RealVectorX &local_f, RealMatrixX &local_jac)
calculate the transverse shear component for the element
const RealVectorX & local_solution(bool if_sens=false) const
virtual void calculate_transverse_shear_residual_sensitivity(const MAST::FunctionBase &p, bool request_jacobian, RealVectorX &local_f, RealMatrixX &local_jac)
calculate the transverse shear component for the element
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...
MAST::StructuralElementBase & _structural_elem
structural element associated with this
unsigned int _shear_quadrature_reduction
reduction in quadrature for shear energy
libMesh::Real Real
const MAST::ElementPropertyCardBase & elem_property()
returns a constant reference to the finite element object
Matrix< Real, Dynamic, Dynamic > RealMatrixX
void vector_mult(T &res, const T &v) const
res = [this] * v
Bending strain operator for 1D element.
Matrix< Real, Dynamic, 1 > RealVectorX
virtual std::unique_ptr< MAST::FEBase > init_side_fe(unsigned int s, bool init_grads, int extra_quadrature_order=0) const
initializes the finite element shape function and quadrature object for the side with the order of qu...
Definition: geom_elem.cpp:159
virtual const std::vector< std::vector< libMesh::RealVectorValue > > & get_dphi() const
Definition: fe_base.cpp:255
void right_multiply_transpose(T &r, const T &m) const
[R] = [this]^T * [M]
void left_multiply(T &r, const T &m) const
[R] = [M] * [this]
MindlinBendingOperator(MAST::StructuralElementBase &elem)
virtual const std::vector< std::vector< Real > > & get_phi() const
Definition: fe_base.cpp:247
const MAST::GeomElem & _elem
element for which bending operator is created
virtual std::unique_ptr< MAST::FEBase > init_fe(bool init_grads, bool init_second_order_derivative, int extra_quadrature_order=0) const
initializes the finite element shape function and quadrature object with the order of quadrature rule...
Definition: geom_elem.cpp:142
virtual void calculate_transverse_shear_residual_boundary_velocity(const MAST::FunctionBase &p, const unsigned int s, const MAST::FieldFunction< RealVectorX > &vel_f, bool request_jacobian, RealVectorX &local_f, RealMatrixX &local_jac)
calculate the transverse shear component for the element