beamel2d.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

#include "beamel2d.h"
#include "global.h"
#include "globmat.h"
#include "genfile.h"
#include "element.h"
#include "node.h"
#include "intpoints.h"
#include <math.h>


beamel2d::beamel2d (void)
{
  long i;
  
  //  number nodes on element
  nne=2;
  //  number of DOFs on element
  ndofe=6;
  //  number of strain/stress components
  tncomp=3;
  //  number of functions approximated
  napfun=3;
  //  order of numerical integration of mass matrix
  intordmm=4;
  //  order of numerical integration of initial stiffness matrix
  intordism=3;
  ssst=plbeam;

  //  number of blocks (parts of geometric matrix)
  nb=1;

  //  number of strain/stress components
  ncomp = new long [nb];
  ncomp[0]=1;
  
  //  cumulative number of components approximated
  cncomp = new long [nb];
  cncomp[0]=0;

  //  number of integration points
  //  order of numerical integration of stiffness matrix
  nip = new long* [nb];
  intordsm = new long* [nb];
  for (i=0;i<nb;i++){
    nip[i] = new long [nb];
    intordsm[i] = new long [nb];
  }
  nip[0][0]=3;
  intordsm[0][0]=3;

  //  total number of integration points
  tnip=3;
}

beamel2d::~beamel2d (void)
{
  long i;
  
  for (i=0;i<nb;i++){
    delete [] nip[i];
    delete [] intordsm[i];
  }
  delete [] nip;
  delete [] intordsm;

  delete [] cncomp;
  delete [] ncomp;
}

void beamel2d::eleminit (long eid)
{
  long ii,jj;
  Mt->elements[eid].nb=nb;
  Mt->elements[eid].intordsm = new long* [nb];
  Mt->elements[eid].nip = new long* [nb];

  for (ii=0;ii<nb;ii++){
    Mt->elements[eid].intordsm[ii] = new long [nb];
    Mt->elements[eid].nip[ii] = new long [nb];
    for (jj=0;jj<nb;jj++){
      Mt->elements[eid].intordsm[ii][jj]=intordsm[ii][jj];
      Mt->elements[eid].nip[ii][jj]=nip[ii][jj];
    }
  }
}

/**
   function assembles %matrix of approximation functions

   ordering of approximation functions
   u - displacement along x
   w - deflection - displacement along z
   phi - rotation around y

   @param n - array containing matrix of approximation functions
   @param xi - natural coordinate from segment <0;1>
   @param l - length of the beam
   @param kappa - 6.0*E*I/k/G/A/l/l
   
   JK, 20.2.2002
*/
void beamel2d::bf_matrix (matrix &n,double xi,double l,double kappa)
{
  vector u(2),w(4),r(4);

  dilat (u.a,xi);
  defl2D_fun (w.a,xi,l,kappa);
  roty_fun (r.a,xi,l,kappa);

  fillm (0.0,n);
  
  n[0][0]=u[0];
  n[0][3]=u[1];
  
  n[1][1]=w[0];
  n[1][2]=w[1];
  n[1][4]=w[2];
  n[1][5]=w[3];
  
  n[2][1]=r[0];
  n[2][2]=r[1];
  n[2][4]=r[2];
  n[2][5]=r[3];
}

/**
   function assembles %matrix of derivatives of approximation function of deflection
   
   @param n - %matrix of derivatives
   @param xi - natural coordinate from segment <0;1>
   @param l - length of the beam
   @param kappa - shear coefficient 6.0*E*I/k/G/A/l/l
   
   JK
*/
void beamel2d::dbf_matrix (matrix &n,double xi,double l,double kappa)
{
  vector dw(4);
  
  der_defl2D_fun (dw.a,xi,l,kappa);
  
  fillm (0.0,n);
  
  n[0][1]=dw[0];
  n[0][2]=dw[1];
  n[0][4]=dw[2];
  n[0][5]=dw[3];
  
}

/**
   function assembles strain-displacement (geometric) %matrix
   
   ordering of strain components
   du/dx - strain e_xx
   phi+dw/dx - strain e_xz
   dphi/dx - curvature
   
   @param gm - geometric %matrix
   @param xi - natural coordinate from segment <0;1>
   @param l - length of the beam
   @param kappa - 6.0*E*I/k/G/A/l/l
   
   JK, 20.2.2002
*/
void beamel2d::geom_matrix (matrix &gm,double xi,double l,double kappa)
{
  vector du(2),dw(4),r(4),dr(4);
  
  der_dilat (du.a,l);
  der_defl2D_fun (dw.a,xi,l,kappa);
  roty_fun (r.a,xi,l,kappa);
  der_roty_fun (dr.a,xi,l,kappa);
  
  fillm (0.0,gm);
  
  gm[0][0] = du[0];
  gm[0][3] = du[1];
  
  gm[1][1] = r[0]+dw[0];
  gm[1][2] = r[1]+dw[1];
  gm[1][4] = r[2]+dw[2];
  gm[1][5] = r[3]+dw[3];
  
  gm[2][1] = dr[0];
  gm[2][2] = dr[1];
  gm[2][4] = dr[2];
  gm[2][5] = dr[3];
}





/**
   function assembles transformation %matrix from local nodal coordinate
   system to the global coordinate system x_g = T x_l
   
   @param nodes - array containing node numbers
   @param tmat - transformation %matrix
   
   JK, 3.2.2002
*/
void beamel2d::transf_matrix (ivector &nodes,matrix &tmat)
{
  long i,n,m;

  fillm (0.0,tmat);

  n=nodes.n;
  m=tmat.m;
  for (i=0;i<m;i++){
    tmat[i][i]=1.0;
  }
  
  for (i=0;i<n;i++){
    if (Mt->nodes[nodes[i]].transf>0){
      tmat[i*3+0][i*3] = Mt->nodes[nodes[i]].e1[0];    tmat[i*3+0][i*3+1] = Mt->nodes[nodes[i]].e2[0];    tmat[i*3+0][i*3+2] = 0.0;
      tmat[i*3+1][i*3] = Mt->nodes[nodes[i]].e1[1];    tmat[i*3+1][i*3+1] = Mt->nodes[nodes[i]].e2[1];    tmat[i*3+1][i*3+2] = 0.0;
      tmat[i*3+2][i*3] = 0.0;                          tmat[i*3+2][i*3+1] = 0.0;                          tmat[i*3+2][i*3+2] = 1.0;
      
    }
  }
}


/**
   function assembles transformation %matrix from local element coordinate
   system to the global problem coordinate system x_g = T x_l
   
   this beam element is formulated in x-z plane, x axis is identical with beam axis
   and it is oriented horizontally, z axis is oriented vertically (downwards)
   
   @param x,z - node coordinates
   @param tmat - transformation %matrix
   
   JK, 3.2.2002
*/
void beamel2d::beam_transf_matrix (long eid,vector &x,vector &z,matrix &tmat)
{
  double c,s,l;
  
  fillm (0.0,tmat);
  
  //  length of the element
  l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
  if (l<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
    fprintf (stderr,"\n in function beamel2d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  //  sine and cosine of element angle
  s=(z[1]-z[0])/l;
  c=(x[1]-x[0])/l;
  
  //  transformation matrix
  tmat[0][0] = c;    tmat[0][1] = -1.0*s;    tmat[0][2] = 0.0;
  tmat[1][0] = s;    tmat[1][1] = c;         tmat[1][2] = 0.0;
  tmat[2][0] = 0.0;  tmat[2][1] = 0.0;       tmat[2][2] = 1.0;

  tmat[3][3] = c;    tmat[3][4] = -1.0*s;    tmat[3][5] = 0.0;
  tmat[4][3] = s;    tmat[4][4] = c;         tmat[4][5] = 0.0;
  tmat[5][3] = 0.0;  tmat[5][4] = 0.0;       tmat[5][5] = 1.0;
}



/**
   function computes stiffness %matrix of two-dimensional beam finite element
   element is based on Mindlin theory (shear stresses are taken into account)
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param sm - stiffness %matrix

   JK, 3.2.2002
*/
void beamel2d::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
{
  long i,ipid,transf;
  double e,g,a,iy,shearcoeff,kappa,l,xi,jac,ww;
  ivector nodes(nne);
  vector x(nne),z(nne),w(intordsm[0][0]),gp(intordsm[0][0]);
  matrix gm(tncomp,ndofe),d(tncomp,tncomp),tmat(ndofe,ndofe);
  
  //  node coordinates in the x-z plane
  Mt->give_node_coord2dxz (x,z,eid);
  //  length of the element
  l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
  if (l<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
    fprintf (stderr,"\n in function beamel2d::stiffness_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  //  coordinates and weights of integration points
  gauss_points (gp.a,w.a,intordsm[0][0]);
  
  
  fillm (0.0,sm);
  ipid=Mt->elements[eid].ipp[ri][ci];
  //  stiffness matrix of the material
  Mm->matstiff (d,ipid);
  //  area of element cross section
  Mc->give_area (eid,a);
  //  moment of inertia of element cross section
  Mc->give_mominer (eid,&iy);
  //  shear coefficient
  Mc->give_shearcoeff (eid,&shearcoeff);
  //  Young's modulus
  e=d[0][0];
  //  shear modulus
  g=d[1][1];
  
  if (shearcoeff<Mp->zero){  kappa=0.0;  shearcoeff=1.0;        }
  else{                      kappa=6.0*e*iy/shearcoeff/g/a/l/l; }
  
  d[0][0]*=a;
  d[1][1]*=a*shearcoeff;
  d[2][2]*=iy;
  
  //  stiffness matrix in local coordinate system
  for (i=0;i<intordsm[0][0];i++){
    xi=(1.0+gp[i])/2.0;  ww=w[i];
    geom_matrix (gm,xi,l,kappa);
    jac=l/2.0*ww;
    bdbj (sm.a,gm.a,d.a,jac,gm.m,gm.n);
  }
  
  //  transformation of stiffness matrix from local element coordinate
  //  system to the global coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  lgmatrixtransf (sm,tmat);

  //  transformation of stiffness matrix from global coordinate system
  //  to the local nodal systems
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
  
}

/**
   function computes stiffness %matrix of two-dimensional beam finite element
   element is based on Mindlin theory (shear stresses are taken into account)
   
   @param eid - number of element
   @param sm - stiffness %matrix

   JK, 3.2.2002
*/
void beamel2d::res_stiffness_matrix (long eid,matrix &sm)
{
  //stiffness_matrix (eid,0,0,sm);
  stiffness_matrix_expl (eid,0,0,sm);
}

/**
   function computes stiffness %matrix of two-dimensional beam finite element
   element is based on Mindlin theory (shear stresses are taken into account)
   stiffness %matrix is formulated explicitly
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param sm - stiffness %matrix

   JK, 3.2.2002
*/
void beamel2d::stiffness_matrix_expl (long eid,long ri,long ci,matrix &sm)
{
  long ipid,transf;
  double e,g,shearcoeff,a,iy,kappa,l;
  ivector nodes(nne);
  vector x(nne),z(nne);
  matrix d(tncomp,tncomp),tmat (ndofe,ndofe);

  
  //  node coordinates in the x-z plane
  Mt->give_node_coord2dxz (x,z,eid);
  
  //  length of the element
  l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
  if (l<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
    fprintf (stderr,"\n in function beamel2d::stiffness_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  
  fillm (0.0,sm);
  ipid=Mt->elements[eid].ipp[ri][ci];
  //  stiffness matrix of the material
  Mm->matstiff (d,ipid);
  //  area of element cross section
  Mc->give_area (eid,a);
  //  moment of inertia of element cross section
  Mc->give_mominer (eid,&iy);
  //  shear coefficient
  Mc->give_shearcoeff (eid,&shearcoeff);
  //  Young's modulus
  e=d[0][0];
  //  shear modulus
  g=d[1][1];
  
  if (shearcoeff<Mp->zero)  kappa=0.0;
  else                      kappa=6.0*e*iy/shearcoeff/g/a/l/l;
  

  //  stiffness matrix in local coordinate system
  sm[0][0]=      e*a/l;
  sm[0][3]= -1.0*e*a/l;
  
  sm[1][1] =  12.0*e*iy/l/l/l;
  sm[1][2] =  -6.0*e*iy/l/l;
  sm[1][4] = -12.0*e*iy/l/l/l;
  sm[1][5] =  -6.0*e*iy/l/l;
  
  sm[2][1] = -6.0*e*iy/l/l;
  sm[2][2] =  4.0*e*iy/l;
  sm[2][4] =  6.0*e*iy/l/l;
  sm[2][5] =  2.0*e*iy/l;
  
  sm[3][0] = -1.0*e*a/l;
  sm[3][3] =      e*a/l;

  sm[4][1] = -12.0*e*iy/l/l/l;
  sm[4][2] =   6.0*e*iy/l/l;
  sm[4][4] =  12.0*e*iy/l/l/l;
  sm[4][5] =   6.0*e*iy/l/l;
  
  sm[5][1] = -6.0*e*iy/l/l;
  sm[5][2] =  2.0*e*iy/l;
  sm[5][4] =  6.0*e*iy/l/l;
  sm[5][5] =  4.0*e*iy/l;
  

  //  transformation of stiffness matrix from local element coordinate
  //  system to the global coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  lgmatrixtransf (sm,tmat);

  //  transformation of stiffness matrix from global coordinate system
  //  to the local nodal systems
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
  
}

/**
   function computes consistent mass %matrix of 2D beam element
   influence of inertial forces from rotations is accounted
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param mm - mass matrix
   
   JK, 3.2.2006
*/
void beamel2d::mass_matrix (long eid,long ri,long ci,matrix &mm)
{
  long i,ipid,transf;
  double xi,ww,jac,e,g,shearcoeff,iy,a,l,kappa,rho;
  ivector nodes(nne);
  vector x(nne),z(nne),w(intordmm),gp(intordmm);
  matrix n(napfun,ndofe),d(tncomp,tncomp),tmat (ndofe,ndofe);


  //  node coordinates in the x-z plane
  Mt->give_node_coord2dxz (x,z,eid);
  gauss_points (gp.a,w.a,intordmm);

  //  length of the element
  l=sqrt ((x[1]-x[0])*(x[1]-x[0])+(z[1]-z[0])*(z[1]-z[0]));
  if (l<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel2d element",eid);
    fprintf (stderr,"\n in function beamel2d::mass_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  
  fillm (0.0,mm);
  ipid=Mt->elements[eid].ipp[ri][ci];

  //  area of element cross section
  Mc->give_area (eid,a);
  //  shear coefficient
  Mc->give_shearcoeff (eid,&shearcoeff);
  //  moment of inertia of element cross section
  Mc->give_mominer (eid,&iy);
  //  density of the material
  Mc->give_densitye (eid,rho);
  //  stiffness matrix of the material
  Mm->matstiff (d,ipid);
  //  Young's modulus
  e=d[0][0];
  //  shear modulus
  g=d[1][1];
  
  if (shearcoeff<Mp->zero)  kappa=0.0;
  else                      kappa=6.0*e*iy/shearcoeff/g/a/l/l;
  
  fillm (0.0,d);
  d[0][0]=1.0;  d[1][1]=1.0;  d[2][2]=iy;
  
  for (i=0;i<intordmm;i++){
    xi=(1.0+gp[i])/2.0;  ww=w[i];
    bf_matrix (n,xi,l,kappa);
    jac=a*l/2.0*ww*rho;
    bdbj (mm.a,n.a,d.a,jac,n.m,n.n);
  }


  //  transformation of mass matrix from local element coordinate
  //  system to the global coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  lgmatrixtransf (mm,tmat);

  //  transformation of mass matrix from global coordinate system
  //  to the local nodal systems
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tmat);
    glmatrixtransf (mm,tmat);
  }
}

/**
   function computes consistent mass %matrix of 2D beam element