20 #ifndef mast_dkt_bending_operator_h 21 #define mast_dkt_bending_operator_h 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" 43 const std::vector<libMesh::Point>& pts):
52 libmesh_assert(e.type() == libMesh::TRI3);
54 _tri6 =
new libMesh::Tri6;
57 for (
unsigned int i=0; i<3; i++)
58 _nodes[i] =
new libMesh::Node;
62 p = *e.node_ptr(0); p += e.point(1); p*= 0.5;
67 p = *e.node_ptr(1); p += e.point(2); p*= 0.5;
72 p = *e.node_ptr(2); p += e.point(0); p*= 0.5;
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];
83 _fe = libMesh::FEBase::build(2, libMesh::FEType(libMesh::SECOND,
84 libMesh::LAGRANGE)).release();
87 _fe->reinit(_tri6, &pts);
97 for (
unsigned int i=0; i<3; i++)
117 const unsigned int qp,
126 const unsigned int qp,
174 const unsigned int qp,
186 const unsigned int qp,
191 const std::vector<std::vector<libMesh::RealVectorValue> >& dphi =
_fe->get_dphi();
192 const unsigned int n_phi = (
unsigned int)dphi.size();
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);
205 for (
unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
206 phi(i_nd) = dphi[i_nd][qp](0);
209 for (
unsigned int i_nd=0; i_nd<n_phi; i_nd++ )
210 phi(i_nd) = dphi[i_nd][qp](1);
214 dbetaxdy += dbetaydx;
218 for (
unsigned int i=0; i<3; i++) {
220 thetax(i) = dbetaxdx(3+i);
221 thetay(i) = dbetaxdx(6+i);
233 for (
unsigned int i=0; i<3; i++) {
235 thetax(i) = dbetaydy(3+i);
236 thetay(i) = dbetaydy(6+i);
247 for (
unsigned int i=0; i<3; i++) {
249 thetax(i) = dbetaxdy(3+i);
250 thetay(i) = dbetaxdy(6+i);
268 libMesh::Point l = *
_tri6->node_ptr(j);
269 l -= *
_tri6->node_ptr(i);
287 libMesh::Point vec0, vec1, vec2, vec3;
290 vec0 = *e.node_ptr(0);
291 vec1 = *e.node_ptr(1);
292 vec2 = *e.node_ptr(2);
296 vec3 = vec1.cross(vec2);
303 vec3 = *e.node_ptr(i);
304 vec2 = *e.node_ptr(j);
307 vec3 = vec2.cross(vec1);
340 Real N1, N2, N3, N4, N5, N6;
343 Real l12, l23, l31, cos4, cos5, cos6, sin4, sin5, sin6;
360 betax(0) = 1.5 * (N5 * sin5 / l31 - N6 * sin6 / l12);
361 betax(1) = 1.5 * (N6 * sin6 / l12 - N4 * sin4 / l23);
362 betax(2) = 1.5 * (N4 * sin4 / l23 - N5 * sin5 / l31);
363 betax(3) = (-.75) * (N5 * sin5 * cos5 + N6 * sin6 * cos6);
364 betax(4) = (-.75) * (N4 * sin4 * cos4 + N6 * sin6 * cos6);
365 betax(5) = (-.75) * (N5 * sin5 * cos5 + N4 * sin4 * cos4);
366 betax(6) = (N1 + N5 * (0.5 * cos5 * cos5 - 0.25 * sin5 * sin5) + N6 * (0.5 * cos6 * cos6 - 0.25 * sin6 * sin6));
367 betax(7) = (N2 + N4 * (0.5 * cos4 * cos4 - 0.25 * sin4 * sin4) + N6 * (0.5 * cos6 * cos6 - 0.25 * sin6 * sin6));
368 betax(8) = (N3 + N5 * (0.5 * cos5 * cos5 - 0.25 * sin5 * sin5) + N4 * (0.5 * cos4 * cos4 - 0.25 * sin4 * sin4));
370 betay(0) = 1.5 * (-N5 * cos5 / l31 + N6 * cos6 / l12);
371 betay(1) = 1.5 * (-N6 * cos6 / l12 + N4 * cos4 / l23);
372 betay(2) = 1.5 * (-N4 * cos4 / l23 + N5 * cos5 / l31);
373 betay(3) = (-N1 + N5 * (0.25 * cos5 * cos5 - 0.5 * sin5 * sin5) + N6 * (0.25 * cos6 * cos6 - 0.5 * sin6 * sin6));
374 betay(4) = (-N2 + N4 * (0.25 * cos4 * cos4 - 0.5 * sin4 * sin4) + N6 * (0.25 * cos6 * cos6 - 0.5 * sin6 * sin6));
375 betay(5) = (-N3 + N5 * (0.25 * cos5 * cos5 - 0.5 * sin5 * sin5) + N4 * (0.25 * cos4 * cos4 - 0.5 * sin4 * sin4));
376 betay(6) = (-1.0) * betax(3);
377 betay(7) = (-1.0) * betax(4);
378 betay(8) = (-1.0) * betax(5);
virtual bool include_transverse_shear_energy() const
returns true if this bending operator supports a transverse shear component
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...
virtual const libMesh::Elem & get_reference_elem() const
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
libMesh::FEBase * _fe
FE object to get shape functions for this element.