MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
dkt_bending_operator.h
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 #ifndef mast_dkt_bending_operator_h
21 #define mast_dkt_bending_operator_h
22 
23 // MAST includes
25 #include "mesh/fe_base.h"
26 #include "mesh/geom_elem.h"
27 
28 // libMesh includes
29 #include "libmesh/elem.h"
30 #include "libmesh/face_tri3.h"
31 #include "libmesh/face_tri6.h"
32 #include "libmesh/node.h"
33 #include "libmesh/fe.h"
34 #include "libmesh/quadrature.h"
35 
36 
37 namespace MAST
38 {
41  public:
43  const std::vector<libMesh::Point>& pts):
44  MAST::BendingOperator2D (elem),
45  _fe (nullptr),
46  _tri6 (nullptr),
47  _nodes (3)
48  {
49  const libMesh::Elem
51 
52  libmesh_assert(e.type() == libMesh::TRI3);
53 
54  _tri6 = new libMesh::Tri6;
55 
56  // next three nodes are created using the corner nodes
57  for (unsigned int i=0; i<3; i++)
58  _nodes[i] = new libMesh::Node;
59 
60  // first edge
61  libMesh::Point p;
62  p = *e.node_ptr(0); p += e.point(1); p*= 0.5;
63  (*_nodes[0]) = p;
64  _nodes[0]->set_id(3);
65 
66  // second edge
67  p = *e.node_ptr(1); p += e.point(2); p*= 0.5;
68  (*_nodes[1]) = p;
69  _nodes[1]->set_id(4);
70 
71  // third edge
72  p = *e.node_ptr(2); p += e.point(0); p*= 0.5;
73  (*_nodes[2]) = p;
74  _nodes[2]->set_id(5);
75 
76  // first three nodes are same as that of the original element
77  for (unsigned int i=0; i<3; i++) {
78  _tri6->set_node( i) = const_cast<libMesh::Node*>(e.node_ptr(i));
79  _tri6->set_node(i+3) = _nodes[i];
80  }
81 
82  // now setup the shape functions
83  _fe = libMesh::FEBase::build(2, libMesh::FEType(libMesh::SECOND,
84  libMesh::LAGRANGE)).release();
85  _fe->get_phi();
86  _fe->get_dphi();
87  _fe->reinit(_tri6, &pts);
88  }
89 
90 
94  virtual ~DKTBendingOperator() {
95  delete _fe;
96  delete _tri6;
97  for (unsigned int i=0; i<3; i++)
98  delete _nodes[i];
99  }
100 
104  virtual bool include_transverse_shear_energy() const {
105  return false;
106  }
107 
108 
115  virtual void
117  const unsigned int qp,
118  FEMOperatorMatrix& Bmat);
119 
124  void
126  const unsigned int qp,
127  const Real z,
128  FEMOperatorMatrix& Bmat_bend);
129 
130  protected:
131 
135  Real _get_edge_length(unsigned int i, unsigned int j);
136 
140  void _get_edge_normal_sine_cosine(unsigned int i,
141  unsigned int j,
142  Real& sine,
143  Real& cosine);
144 
145 
147  RealVectorX& betax,
148  RealVectorX& betay);
149 
150 
154  libMesh::FEBase* _fe;
155 
159  libMesh::Elem* _tri6;
160 
164  std::vector<libMesh::Node*> _nodes;
165  };
166 }
167 
168 
169 
170 inline
171 void
174  const unsigned int qp,
175  FEMOperatorMatrix& Bmat) {
176 
177  this->initialize_bending_strain_operator_for_z(fe, qp, 1., Bmat);
178 }
179 
180 
181 
182 inline
183 void
186  const unsigned int qp,
187  const Real z,
188  FEMOperatorMatrix& Bmat) {
189 
190  // the DKT operator uses its own fe object, and not the one from this argument
191  const std::vector<std::vector<libMesh::RealVectorValue> >& dphi = _fe->get_dphi();
192  const unsigned int n_phi = (unsigned int)dphi.size();
193 
195  phi = RealVectorX::Zero(n_phi),
196  dbetaxdx = RealVectorX::Zero(9),
197  dbetaxdy = RealVectorX::Zero(9),
198  dbetaydx = RealVectorX::Zero(9),
199  dbetaydy = RealVectorX::Zero(9),
200  w = RealVectorX::Zero(3),
201  thetax = RealVectorX::Zero(3),
202  thetay = RealVectorX::Zero(3);
203 
204 
205  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
206  phi(i_nd) = dphi[i_nd][qp](0); // dphi/dx
207  _calculate_dkt_shape_functions(phi, dbetaxdx, dbetaydx);
208 
209  for ( unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
210  phi(i_nd) = dphi[i_nd][qp](1); // dphi/dy
211  _calculate_dkt_shape_functions(phi, dbetaxdy, dbetaydy);
212 
213  // add the shear components together
214  dbetaxdy += dbetaydx;
215 
216 
217  // now set values for individual dofs
218  for (unsigned int i=0; i<3; i++) {
219  w(i) = dbetaxdx( i);
220  thetax(i) = dbetaxdx(3+i);
221  thetay(i) = dbetaxdx(6+i);
222  }
223 
224  w *= z;
225  thetax *= z;
226  thetay *= z;
227 
228  Bmat.set_shape_function(0, 2, w); // epsilon-x: w
229  Bmat.set_shape_function(0, 3, thetax); // epsilon-x: tx
230  Bmat.set_shape_function(0, 4, thetay); // epsilon-x: ty
231 
232 
233  for (unsigned int i=0; i<3; i++) {
234  w(i) = dbetaydy( i);
235  thetax(i) = dbetaydy(3+i);
236  thetay(i) = dbetaydy(6+i);
237  }
238 
239  w *= z;
240  thetax *= z;
241  thetay *= z;
242 
243  Bmat.set_shape_function(1, 2, w); // epsilon-x: w
244  Bmat.set_shape_function(1, 3, thetax); // epsilon-x: tx
245  Bmat.set_shape_function(1, 4, thetay); // epsilon-x: ty
246 
247  for (unsigned int i=0; i<3; i++) {
248  w(i) = dbetaxdy( i);
249  thetax(i) = dbetaxdy(3+i);
250  thetay(i) = dbetaxdy(6+i);
251  }
252 
253  w *= z;
254  thetax *= z;
255  thetay *= z;
256 
257  Bmat.set_shape_function(2, 2, w); // epsilon-x: w
258  Bmat.set_shape_function(2, 3, thetax); // epsilon-x: tx
259  Bmat.set_shape_function(2, 4, thetay); // epsilon-x: ty
260 }
261 
262 
263 
264 inline
265 Real
266 MAST::DKTBendingOperator::_get_edge_length(unsigned int i, unsigned int j)
267 {
268  libMesh::Point l = *_tri6->node_ptr(j);
269  l -= *_tri6->node_ptr(i);
270 
271  return l.norm();
272 }
273 
274 
275 
276 
277 inline
278 void
280  unsigned int j,
281  Real& sine,
282  Real& cosine)
283 {
284  const libMesh::Elem
285  &e = _elem.get_reference_elem();
286 
287  libMesh::Point vec0, vec1, vec2, vec3;
288 
289  // calculate the normal to the element
290  vec0 = *e.node_ptr(0);
291  vec1 = *e.node_ptr(1);
292  vec2 = *e.node_ptr(2);
293 
294  vec1 -= vec0;
295  vec2 -= vec0;
296  vec3 = vec1.cross(vec2);
297  // this is the unit vector normal to the plane of the triangle
298  vec1 = vec3;
299  vec1 /= vec1.norm();
300 
301  // cross product of the length vector with the surface normal will
302  // give the needed vector
303  vec3 = *e.node_ptr(i);
304  vec2 = *e.node_ptr(j);
305 
306  vec2 -= vec3;
307  vec3 = vec2.cross(vec1);
308  vec3 /= vec3.norm(); // this is the unit vector needed
309 
310  // cos of angle between this and the x-axis is simply the
311  // 0th component of this vector
312  cosine = vec3(0);
313  sine = vec3(1);
314 }
315 
316 
317 
318 
319 inline
320 void
322  RealVectorX& betax,
323  RealVectorX& betay)
324 {
325  // -- keep in mind that the index numbers for the elems start at 0.
326  // -- also, the mid side node numbers in the Batoz's paper are different from
327  // the ones used in this library. Hence, use the following association
328  // BATOZ TRI6 node # MAST TRI6 node #
329  // 1 0
330  // 2 1
331  // 3 2
332  // 4 (on edge 2,3) 4 (on edge 1,2)
333  // 5 (on edge 3,1) 5 (on edge 2,0)
334  // 6 (on edge 1,2) 3 (on edge 0,1)
335  // -- all shape functions for this element are in a variable major format. Hence,
336  // they follow Hx for w1, w2, w3, thetax1, thetax2, thetax3, thetay1, thetay2, thetay3.
337  // And then the same this is followed for Hy (from dofs 9-17)
338 
339  // local variables for shape functions
340  Real N1, N2, N3, N4, N5, N6;
341 
342  // local variables for edge lengths and sine/cosines
343  Real l12, l23, l31, cos4, cos5, cos6, sin4, sin5, sin6;
344 
345  N1 = phi(0);
346  N2 = phi(1);
347  N3 = phi(2);
348  N4 = phi(4);
349  N5 = phi(5);
350  N6 = phi(3);
351 
352  _get_edge_normal_sine_cosine(1, 2, sin4, cos4);
353  _get_edge_normal_sine_cosine(2, 0, sin5, cos5);
354  _get_edge_normal_sine_cosine(0, 1, sin6, cos6);
355 
356  l12 = _get_edge_length(0, 1);
357  l23 = _get_edge_length(1, 2);
358  l31 = _get_edge_length(2, 0);
359 
360  betax(0) = 1.5 * (N5 * sin5 / l31 - N6 * sin6 / l12); // Hx, w1
361  betax(1) = 1.5 * (N6 * sin6 / l12 - N4 * sin4 / l23); // Hx, w2
362  betax(2) = 1.5 * (N4 * sin4 / l23 - N5 * sin5 / l31); // Hx, w3
363  betax(3) = (-.75) * (N5 * sin5 * cos5 + N6 * sin6 * cos6); // Hx, thetax1
364  betax(4) = (-.75) * (N4 * sin4 * cos4 + N6 * sin6 * cos6); // Hx, thetax2
365  betax(5) = (-.75) * (N5 * sin5 * cos5 + N4 * sin4 * cos4); // Hx, thetax3
366  betax(6) = (N1 + N5 * (0.5 * cos5 * cos5 - 0.25 * sin5 * sin5) + N6 * (0.5 * cos6 * cos6 - 0.25 * sin6 * sin6)); // Hx, thetay1
367  betax(7) = (N2 + N4 * (0.5 * cos4 * cos4 - 0.25 * sin4 * sin4) + N6 * (0.5 * cos6 * cos6 - 0.25 * sin6 * sin6)); // Hx, thetay2
368  betax(8) = (N3 + N5 * (0.5 * cos5 * cos5 - 0.25 * sin5 * sin5) + N4 * (0.5 * cos4 * cos4 - 0.25 * sin4 * sin4)); // Hx, thetay3
369 
370  betay(0) = 1.5 * (-N5 * cos5 / l31 + N6 * cos6 / l12); // Hy, w1
371  betay(1) = 1.5 * (-N6 * cos6 / l12 + N4 * cos4 / l23); // Hy, w2
372  betay(2) = 1.5 * (-N4 * cos4 / l23 + N5 * cos5 / l31); // Hy, w3
373  betay(3) = (-N1 + N5 * (0.25 * cos5 * cos5 - 0.5 * sin5 * sin5) + N6 * (0.25 * cos6 * cos6 - 0.5 * sin6 * sin6)); // Hy, thetax1
374  betay(4) = (-N2 + N4 * (0.25 * cos4 * cos4 - 0.5 * sin4 * sin4) + N6 * (0.25 * cos6 * cos6 - 0.5 * sin6 * sin6)); // Hy, thetax2
375  betay(5) = (-N3 + N5 * (0.25 * cos5 * cos5 - 0.5 * sin5 * sin5) + N4 * (0.25 * cos4 * cos4 - 0.5 * sin4 * sin4)); // Hy, thetax3
376  betay(6) = (-1.0) * betax(3); // Hy, thetay1
377  betay(7) = (-1.0) * betax(4); // Hy, thetay2
378  betay(8) = (-1.0) * betax(5); // Hy, thetay3
379 }
380 
381 
382 #endif
virtual const libMesh::Elem & get_reference_elem() const
Definition: geom_elem.cpp:54
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...
libMesh::Real Real
Bending strain operator for 1D element.
Real _get_edge_length(unsigned int i, unsigned int j)
returns the length of the side defined by vector from node i to j
Matrix< Real, Dynamic, 1 > RealVectorX
libMesh::Elem * _tri6
6 noded triangle that is used to calculate the shape functions
void initialize_bending_strain_operator_for_z(const MAST::FEBase &fe, const unsigned int qp, const Real z, FEMOperatorMatrix &Bmat_bend)
initializes the bending strain operator for the specified quadrature point and z-location.
void _get_edge_normal_sine_cosine(unsigned int i, unsigned int j, Real &sine, Real &cosine)
returns the cos of normal to the side defined by vector from node i to j
virtual void initialize_bending_strain_operator(const MAST::FEBase &fe, const unsigned int qp, FEMOperatorMatrix &Bmat)
initialze the bending strain operator for DKt element, withouth the z-location.
DKTBendingOperator(MAST::StructuralElementBase &elem, const std::vector< libMesh::Point > &pts)
void _calculate_dkt_shape_functions(const RealVectorX &phi, RealVectorX &betax, RealVectorX &betay)
std::vector< libMesh::Node * > _nodes
vector to store node pointers for the mid-sized nodes
const MAST::GeomElem & _elem
element for which bending operator is created
virtual ~DKTBendingOperator()
destructor
virtual bool include_transverse_shear_energy() const
returns true if this bending operator supports a transverse shear component
libMesh::FEBase * _fe
FE object to get shape functions for this element.