beamel3d.cpp

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

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


beamel3d::beamel3d (void)
{
  long i,j;
  
  //  number of nodes on element
  nne=2;
  //  number of DOFs on element
  ndofe=12;
  //  number of strain/stress components
  //  in this case, displacements/rotations/forces/moments are used
  tncomp=6;
  //  number of functions approximated
  napfun=6;
  //  order of numerical integration of mass matrix
  intordmm=4;

  intordism=2;
  ssst=spacebeam;

  //  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];
  }
  
  //  there are two integration points
  //  each integration point stores nodal forces and moments at one node
  nip[0][0]=2;
  intordsm[0][0]=1;
  
  tnip=0;
  for (i=0;i<nb;i++){
    for (j=0;j<nb;j++){
      tnip+=nip[i][j];
    }
  }
}

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

/**
   function initiates data on element
   
   JK
*/
void beamel3d::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 geometric %matrix
   (function assembles %matrix of derivatives of approximation functions)
   
   
   ordering of approximation functions

   du/dx - strain e_xx
   phi+dv/dx - strain e_yz
   phi+dw/dx - strain e_xz
   dphix/dx - curvature
   dphiy/dx - curvature
   dphiz/dx - curvature

   @param n - array containing geometric %matrix
   @param s - natural coordinate from segment <0;1>
   @param dl - length of the beam
   @param gy - 6.0*E*I_y/k/G/A/l/l
   @param gz - 6.0*E*I_z/k/G/A/l/l
   
   PF, 20.12.2002, revised by JK 2.9.2006
*/
void beamel3d::geom_matrix (matrix &n,double s,double dl,double gy,double gz)
{
  double aj1,ll;

  fillm (0.0,n);
  ll=dl*dl;
  
  //    Vnitrni sily   N, Vy, Vz, Mx, My, Mz
  
  //  N
  n[0][0]=-1./dl;
  n[0][6]= 1./dl;

  //  Qy
  n[1][1] =-2.*gz/dl/(1.+2.*gz);
  n[1][5] =-gz/(1.+2.*gz);
  n[1][7] =-n[1][1];
  n[1][11]= n[1][5];

  //  Qz
  n[2][2] =-2.*gy/dl/(1.+2.*gy);
  n[2][4] = gy/(1.+2.*gy);
  n[2][8] =-n[2][2];
  n[2][10]= n[2][4];

  //  Mx
  n[3][3] =-1./dl;
  n[3][9] = 1./dl;
  
  //  My
  aj1=1./(1.+2*gy);
  n[4][2] = aj1*(6.-12.*s)/ll;
  n[4][4] = aj1*(-2.*(2.+gy)+6.*s)/dl;
  n[4][8] =-n[4][2];
  n[4][10]= aj1*(-2.*(1.-gy)+6.*s)/dl;

  //  Mz
  aj1=1./(1.+2.*gz);
  n[5][1] =-aj1*( 6.-12.*s)/ll;
  n[5][5] = aj1*(-2.*(2.+gz)+6.*s)/dl;
  n[5][7] =-n[5][1];
  n[5][11]= aj1*(-2.*(1.-gz)+6.*s)/dl;
}

/**
   function assembles %matrix of approximation functions
   
   ordering of approximation functions
   u v w - displacement along x,y,z
   fix,fiy,fiz - rotation around x,y,z
   
   @param n - %matrix of approximation functions
   @param s - natural coordinate from segment <0;1>
   @param dl - length of the beam
   @param gy - 6.0*E*I_y/k/G/A/l/l
   @param gz - 6.0*E*I_z/k/G/A/l/l
   
   PF, 20.12.2002, revised by JK 2.9.2006
*/
void beamel3d::bf_matrix (matrix &n,double s,double dl,double gy,double gz)
{
  double ll,aj1;
  ll=dl*dl;
  
  //  u
  n[0][0] = 1.-s;
  n[0][6] = s;

  //  v
  aj1=1./(1.+2.*gz);
  n[1][1] = aj1*(1.+2.*gz-2.*gz*s-3.*s*s+2.*s*s*s);
  n[1][5] =-aj1*(-(1.+gz)*s+(2.+gz)*s*s-s*s*s)*dl;
  n[1][7] = aj1*(2.*gz*s+3.*s*s-2.*s*s*s);
  n[1][11]=-aj1*(gz*s+(1.-gz)*s*s-s*s*s)*dl;

  //  w
  aj1=1./(1.+2.*gy);
  n[2][2] =aj1*(1.+2.*gy-2.*gy*s-3.*s*s+2*s*s*s);
  n[2][4] =aj1*(-(1.+gy)*s+(2.+gy)*s*s-s*s*s)*dl;
  n[2][8] =aj1*(2.*gy*s+3.*s*s-2.*s*s*s);
  n[2][10]=aj1*(gy*s+(1.-gy)*s*s-s*s*s)*dl;

  // fx
  n[3][3] = 1.-s;
  n[3][9] = s;

  // fy
  n[4][2] = aj1*(6.*s-6.*s*s)/dl;
  n[4][4] = aj1*(1.+2.*gy-2.*(2.+gy)*s+3.*s*s);
  n[4][8] =-n[4][2];
  n[4][10]= aj1*(-2.*(1.-gy)*s+3.*s*s);
  
  // fz
  n[5][1] =-aj1*(6.*s-6.*s*s)/dl;
  n[5][5] = aj1*(1.+2.*gz-2.*(2.+gz)*s+3.*s*s);
  n[5][7] =-n[5][1];
  n[5][11]= aj1*(-2.*(1.-gz)*s+3.*s*s);
}


/**
   function assembles %matrix of approximation functions for translational DOFs
   
   ordering of approximation functions
   u v w - displacement along x,y,z
   fix,fiy,fiz - rotation around x,y,z are neglected
   
   ordering of nodal values
   u_1, v_1, w_1, phi_x1, phi_y1, phi_z1, u_2, v_2, w_2, phi_x2, phi_y2, phi_z2

   @param n - %matrix of approximation functions
   @param xi - natural coordinate from segment <0;1>
   @param l - length of the beam
   @param gy - 6.0*E*I_y/k/G/A/l/l
   @param gz - 6.0*E*I_z/k/G/A/l/l
   
   JK, 16.11.2006
*/
void beamel3d::bf_matrix_transl (matrix &n,double xi,double l,double gy,double gz)
{
  fillm (0.0,n);

  //  displacement u (along x)
  n[0][0] = 1.0-xi;
  n[0][6] = xi;
  
  //  displacement v (along y)
  n[1][1] = (1.0 + 2.0*gz - 2.0*gz*xi - 3.0*xi*xi + 2.0*xi*xi*xi)/(1.0 + 2.0*gz);
  n[1][5] = ((1.0+gz)*xi - (2.0+gz)*xi*xi + xi*xi*xi)*l/(1.0+2.0*gz);
  n[1][7] = (2.0*gz*xi + 3.0*xi*xi - 2.0*xi*xi*xi)/(1.0+2.0*gz);
  n[1][11]= (-1.0*gz*xi - (1.0-gz)*xi*xi - xi*xi*xi)*l/(1.0+2.0*gz);
  
  //  displacement w (along z)
  n[2][2] = (1.0+2.0*gy - 2.0*gy*xi - 3.0*xi*xi + 2.0*xi*xi*xi)/(1.0+2.0*gy);
  n[2][4] = (-1.0*(1.0+gy)*xi + (2.0+gy)*xi*xi - xi*xi*xi)*l/(1.0+2.0*gy);
  n[2][8] = (2.0*gy*xi + 3.0*xi*xi - 2.0*xi*xi*xi)/(1.0+2.0*gy);
  n[2][10]= (gy*xi + (1.0-gy)*xi*xi - xi*xi*xi)*l/(1.0+2.0*gy);
}



/**
   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
   
   PF, 20.12.2002, revised by JK, 28.11.2006
*/
void beamel3d::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*6+0][i*6] = Mt->nodes[nodes[i]].e1[0];
      tmat[i*6+1][i*6] = Mt->nodes[nodes[i]].e1[1];
      tmat[i*6+2][i*6] = Mt->nodes[nodes[i]].e1[2];

      tmat[i*6+0][i*6+1] = Mt->nodes[nodes[i]].e2[0];
      tmat[i*6+1][i*6+1] = Mt->nodes[nodes[i]].e2[1];
      tmat[i*6+2][i*6+1] = Mt->nodes[nodes[i]].e2[2];

      tmat[i*6+0][i*6+2] = Mt->nodes[nodes[i]].e3[0];
      tmat[i*6+1][i*6+2] = Mt->nodes[nodes[i]].e3[1];
      tmat[i*6+2][i*6+2] = Mt->nodes[nodes[i]].e3[2];
      
      
      tmat[i*6+3][i*6+3] = Mt->nodes[nodes[i]].e1[0];
      tmat[i*6+4][i*6+3] = Mt->nodes[nodes[i]].e1[1];
      tmat[i*6+5][i*6+3] = Mt->nodes[nodes[i]].e1[2];
      
      tmat[i*6+3][i*6+4] = Mt->nodes[nodes[i]].e2[0];
      tmat[i*6+4][i*6+4] = Mt->nodes[nodes[i]].e2[1];
      tmat[i*6+5][i*6+4] = Mt->nodes[nodes[i]].e2[2];

      tmat[i*6+3][i*6+5] = Mt->nodes[nodes[i]].e3[0];
      tmat[i*6+4][i*6+5] = Mt->nodes[nodes[i]].e3[1];
      tmat[i*6+5][i*6+5] = Mt->nodes[nodes[i]].e3[2];
    }
  }
}


/**
   function assembles transformation %matrix from local to global system x_g = T x_l
   
   columns of the transformation %matrix are created by coordinates
   of local base vectors expressed in global system
   
   @param tmat - transformation %matrix
   @param dl - lenght of the beam
   @param vec - direction %vector of local axis z (in global coordinates)
   @param x,y,z - vectors of nodal coordinates in global system
   @param eid - element id
   
   JK, 3.2.2002, revised 1.9.2006
*/
void beamel3d::beam_transf_matrix (matrix &tmat,double &dl,vector &vec,vector &x,vector &y,vector &z,long eid)
{
  double c;

  fillm (0.0,tmat);
  
  //  local vector x_l
  tmat[0][0]=x[1]-x[0];
  tmat[1][0]=y[1]-y[0];
  tmat[2][0]=z[1]-z[0];
  
  //  length of the beam, it is equal to the norm of x_l
  dl=sqrt((tmat[0][0]*tmat[0][0]+tmat[1][0]*tmat[1][0]+tmat[2][0]*tmat[2][0]));
  
  if (dl<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel3d element",eid);
    fprintf (stderr,"\n in function beamel3d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  //  normed local vector x_l
  tmat[0][0]=tmat[0][0]/dl;
  tmat[1][0]=tmat[1][0]/dl;
  tmat[2][0]=tmat[2][0]/dl;
  
  //  local vector y_l
  tmat[0][1]=vec[1]*tmat[2][0]-vec[2]*tmat[1][0];
  tmat[1][1]=vec[2]*tmat[0][0]-vec[0]*tmat[2][0];
  tmat[2][1]=vec[0]*tmat[1][0]-vec[1]*tmat[0][0];
  
  //  norm of the local vector y_l
  c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));

  if (c<Mp->zero){
    vec[0]=0.0;
    vec[1]=1.0;
    vec[2]=0.0;
    
    //  local vector z_l (defined by vector product)
    tmat[0][2]=tmat[1][0]*vec[2]-tmat[2][0]*vec[1];
    tmat[1][2]=tmat[2][0]*vec[0]-tmat[0][0]*vec[2];
    tmat[2][2]=tmat[0][0]*vec[1]-tmat[1][0]*vec[0];

    //  norm of the local vector z_l
    c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));

    //  normed local vector z_l
    tmat[0][2]=tmat[0][2]/c;
    tmat[1][2]=tmat[1][2]/c;
    tmat[2][2]=tmat[2][2]/c;
    
    if (c<Mp->zero){
      fprintf (stderr,"\n\n zero z base vector of the %ld beamel3d element",eid);
      fprintf (stderr,"\n in function beamel3d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
    }
    
    //  local vector y_l
    tmat[0][1]=tmat[1][2]*tmat[2][0]-tmat[2][2]*tmat[1][0];
    tmat[1][1]=tmat[2][2]*tmat[0][0]-tmat[0][2]*tmat[2][0];
    tmat[2][1]=tmat[0][2]*tmat[1][0]-tmat[1][2]*tmat[0][0];

    //  norm of the local vector y_l
    c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));

    if (c<Mp->zero){
      fprintf (stderr,"\n\n zero y base vector of the %ld beamel3d element",eid);
      fprintf (stderr,"\n in function beamel3d::beam_transf_matrix (file %s, line %d).\n",__FILE__,__LINE__);
    }
    
    //  normed local vector y_l
    tmat[0][1]=tmat[0][1]/c;
    tmat[1][1]=tmat[1][1]/c;
    tmat[2][1]=tmat[2][1]/c;
    
  }
  else{
    
    //  normed local vector y_l
    tmat[0][1]=tmat[0][1]/c;
    tmat[1][1]=tmat[1][1]/c;
    tmat[2][1]=tmat[2][1]/c;
    
    //  local vector z_l (defined by vector product)
    tmat[0][2]=tmat[1][0]*tmat[2][1]-tmat[2][0]*tmat[1][1];
    tmat[1][2]=tmat[2][0]*tmat[0][1]-tmat[0][0]*tmat[2][1];
    tmat[2][2]=tmat[0][0]*tmat[1][1]-tmat[1][0]*tmat[0][1];
    
    //  norm of the local vector z_l
    c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));
    
    if (c<Mp->zero){
      fprintf (stderr,"\n\n zero z base vector of the %ld beamel3d element",eid);
      fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
    }
    
    //  normed local vector z_l
    tmat[0][2]=tmat[0][2]/c;
    tmat[1][2]=tmat[1][2]/c;
    tmat[2][2]=tmat[2][2]/c;
  }
  
  tmat[3][3]  = tmat[0][0];  tmat[3][4]   = tmat[0][1];  tmat[3][5]   = tmat[0][2];
  tmat[4][3]  = tmat[1][0];  tmat[4][4]   = tmat[1][1];  tmat[4][5]   = tmat[1][2];
  tmat[5][3]  = tmat[2][0];  tmat[5][4]   = tmat[2][1];  tmat[5][5]   = tmat[2][2];

  tmat[6][6]  = tmat[0][0];  tmat[6][7]   = tmat[0][1];  tmat[6][8]   = tmat[0][2];
  tmat[7][6]  = tmat[1][0];  tmat[7][7]   = tmat[1][1];  tmat[7][8]   = tmat[1][2];
  tmat[8][6]  = tmat[2][0];  tmat[8][7]   = tmat[2][1];  tmat[8][8]   = tmat[2][2];

  tmat[9][9]  = tmat[0][0];  tmat[9][10]  = tmat[0][1];  tmat[9][11]  = tmat[0][2];
  tmat[10][9] = tmat[1][0];  tmat[10][10] = tmat[1][1];  tmat[10][11] = tmat[1][2];
  tmat[11][9] = tmat[2][0];  tmat[11][10] = tmat[2][1];  tmat[11][11] = tmat[2][2];
}




/**
   function computes stiffness %matrix of 3D beam element
   
   @param eid - number of element
   @param sm - stiffness %matrix
   
   JK, 3.2.2002
*/
void beamel3d::res_stiffness_matrix (long eid,matrix &sm)
{
  stiffness_matrix (eid,0,0,sm);
  
}

/**
   function computes stiffness %matrix of 3D beam element
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param sm - stiffness %matrix
   
   PF, 20.12.2002
*/
void beamel3d::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
{
  long ii,i,i1,transf;
  double e,g,a,ixyz[3],beta[2],dl,ll,gy,gz,aj1;
  ivector nodes(nne);
  vector vec(3),x(nne),y(nne),z(nne);
  matrix d(tncomp,tncomp),tran(ndofe,ndofe);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
  Mc->give_vectorlcs (eid,vec);
  beam_transf_matrix (tran,dl,vec,x,y,z,eid);
  ll=dl*dl;
  
  fillm (0.0,sm);
  ii=Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (d,ii);
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  e=d[0][0];  g=d[1][1];
  
  //  stiffness matrix in local coordinate system
  //  direct s Iy
  if(beta[0]<Mp->zero) gy=0.0;
  else                 gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
  //  direct s Iz
  if(beta[1]<Mp->zero) gz=0.;
  else                 gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
  // axial forces and torsion
  sm[0][0]=  e*a/dl;
  sm[0][6]= -sm[0][0];
  sm[3][3]=  g*ixyz[0]/dl;
  sm[3][9]= -sm[3][3];
  //  direct s Iy
  aj1=ixyz[1]*e/(1.+2*gy);
  sm[2][2] = 12.*aj1/ll/dl;
  sm[2][4] =-6*aj1/ll;
  sm[2][8] =-sm[2][2];
  sm[2][10]= sm[2][4];
  sm[4][4] = 2.*(2.+gy)*aj1/dl;
  sm[4][8] = 6*aj1/ll;
  sm[4][10]= 2.*(1.-gy)*aj1/dl;
  //  direct s Iz
  aj1=ixyz[2]*e/(1.+2.*gz);
  sm[1][1] = 12.*aj1/ll/dl;
  sm[1][5] = 6.*aj1/ll;
  sm[1][7] =-sm[1][1];
  sm[1][11]= sm[1][5];
  sm[5][5] = 2.*(2.+gz)*aj1/dl;
  sm[5][7] =-6.*aj1/ll;
  sm[5][11]= 2.*(1.-gz)*aj1/dl;
  
  for (i=0;i<6;i++){
    i1=i+6;
    sm[i1][i1]=sm[i][i];
    sm[i1][i]=sm[i][i1];
  }
  
  sm[4][2] = sm[2][4];
  sm[5][1] = sm[1][5];
  sm[7][5] = sm[5][7];
  sm[8][4] = sm[4][8];
  sm[11][1]= sm[1][11];
  sm[10][2]= sm[2][10];
  sm[7][11]=-sm[1][5];
  sm[11][7]=-sm[1][5];
  sm[8][10]=-sm[2][4];
  sm[10][8]=-sm[2][4];
  
  //  transformation of stiffness matrix from local element coordinate