beamgen3d.cpp

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

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


beamgen3d::beamgen3d (void)
{
  long i,j;
  
  nne=2;  ndofe=14;  tncomp=8;  napfun=8;
  nb=1;  intordmm=4;  intordism=2;  ssst=spacebeam;

  ncomp = new long [nb];
  ncomp[0]=1;
  
  cncomp = new long [nb];
  cncomp[0]=0;

  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];
    }
  }
}

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

/**
   function initiates data on element
   
   JK
*/
void beamgen3d::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 transformation %matrix from global to nodal system x_n = T x_g
   
   @param nodes - array containing node numbers
   @param tmat - transformation %matrix
   
   PF, 20.12.2002
*/
void beamgen3d::transf_matrix (ivector &nodes,matrix &tmat)
{
  long i,i6,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){
      i6=i*7;
      tmat[i6][i6]   = Mt->nodes[nodes[i]].e1[0];  tmat[i6][i6+1]   = Mt->nodes[nodes[i]].e2[0];  tmat[i6][i6+2]   = Mt->nodes[nodes[i]].e3[0];
      tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1];  tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1];  tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
      tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2];  tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2];  tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
      i6=i*7+3;
      tmat[i6][i6]   = Mt->nodes[nodes[i]].e1[0];  tmat[i6][i6+1]   = Mt->nodes[nodes[i]].e2[0];  tmat[i6][i6+2]   = Mt->nodes[nodes[i]].e3[0];
      tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1];  tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1];  tmat[i6+1][i6+2] = Mt->nodes[nodes[i]].e3[1];
      tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2];  tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2];  tmat[i6+2][i6+2] = Mt->nodes[nodes[i]].e3[2];
      tmat[i6+3][i6] = 1.;      
    }
  }
}


/**
   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 beamgen3d::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 beamgen3d element",eid);
    fprintf (stderr,"\n in function beamgen3d::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 beamgen3d element",eid);
      fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
    }
    
    //  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 beamgen3d element",eid);
      fprintf (stderr,"\n in function beamgen3d::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 beamgen3d 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;
  }
}
void beamgen3d::ck_matrix (matrix &ck,  double s,double c,double eh,double dl)
{
  double eh1,c1,cjm;

  fillm (0.0,ck);
  eh1=eh*dl;
  
// SESTAVENI MATICE Ck konstant
      c1=c-1.;
      cjm=c1*c1-(s-eh1)*s;

      ck[0][3]=s/cjm;
      ck[0][2]=-ck[0][3]*c1/s;
      ck[0][1]=-eh*ck[0][3];
      ck[0][0]=1.-ck[0][2];

      ck[1][3]= (eh1*s-c1)/eh/cjm;
      ck[1][2]= -(1./eh+ck[1][3]*c1)/s;
      ck[1][1]= 1.-eh*ck[1][3];
      ck[1][0]=-ck[1][2];

      ck[2][3]=-ck[0][3];
      ck[2][2]=-ck[2][3]*c1/s;
      ck[2][1]=-eh*ck[2][3];
      ck[2][0]=-ck[2][2];

      ck[3][3]= c1/eh/cjm;
      ck[3][2]=-ck[3][3]*(s-eh1)/c1;
      ck[3][1]=-eh*ck[3][3];
      ck[3][0]=-ck[3][2];
 }

/**
   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.9.2006
*/
void beamgen3d::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[0][0]=-1./dl;
     n[0][6]= 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;
// 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];
// 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;
// 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];
// Mx
     n[3][3] =-1./dl;
     n[3][9] = 1./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.9.2006
*/
void beamgen3d::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;
//  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;
// 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);

//  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;
// 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);
// fx
        n[3][3] = 1.-s;
        n[3][9] = s;
  
}






/**
   function computes stiffness %matrix of 3D beam element
   
   @param eid - number of element
   @param sm - stiffness %matrix
   
   PF 20.9.2006
*/
void beamgen3d::res_stiffness_matrix (long eid,matrix &sm)
{
  stiffness_matrix (eid,0,0,sm);
  
  long i;
  for (i=0;i<sm.m;i++){
    if (sm[i][i]<Mp->zero){
      fprintf (stderr,"\n\n nulovy diagonalni prvek v matici tuhosti jednoho prvku (file %s, line %d)",__FILE__,__LINE__);
    }
  }  
}

/**
   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.9.2006
*/
void beamgen3d::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
{
  long ipp,i,transf;
  double e,g,a,*ixyz,*ioyz,*beta,dl,ll,gy,gz,aj1, integn4;
  ivector nodes(nne);
  vector vec(3),x(nne),y(nne),z(nne),integn1(4),integn2(3),integn3(2);
  matrix d(tncomp,tncomp),tran(3,3);
  
  ixyz = new double [3];
  ioyz = new double [3];
  beta = new double [2];

  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);
  ipp=Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (d,ipp);
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  // vysecovy Io, Ioy, Ioz
//!!!  Mc->give_iomega (eid,ioyz);
   ioyz[0]=4.2692e-6 ; ioyz[1]=-2.6687e-5;  ioyz[2]=0.0; ixyz[0]=2.646e-7; ixyz[1]=2.135e-4; ixyz[2]=3.333e-5; 
 ioyz[0]=1.1924e-6 ; ioyz[1]=0.0;  ioyz[2]=7.5625e-6; ixyz[0]=2.333e-7; ixyz[1]=3.048e-5; ixyz[2]=6.25e-5; 
//  ioyz[0]=1.1924e-6 ; ioyz[1]=7.5625e-6;  ioyz[2]=0.0;  ixyz[0]=2.333e-7; ixyz[1]=6.25e-5; ixyz[2]=3.048e-5; 
//  ioyz[0]=0.2773e-6 ; ioyz[1]=0.0;  ioyz[2]=0.0;
//  ioyz[0]=1. ; ioyz[1]=0.;  ioyz[2]=0.0; ixyz[0]=1.; ixyz[1]=1.; ixyz[2]=1.;

  e=d[0][0];  g=d[1][1];

//  stiffness matrix in local coordinate system
// axial forces and torsion
  sm[0][0]=  e*a/dl; sm[0][7]= -sm[0][0]; sm[7][7]= sm[0][0];
  //sm[3][3]=  g*ixyz[0]/dl; sm[3][10]= -sm[3][3]; sm[10][10]= sm[3][3];
//  direct s Iy   w, fy
  if(beta[0]<Mp->zero) gy=0.0;
  else                 gy=6.0*e*ixyz[1]/beta[0]/g/a/ll;
  if(beta[1]<Mp->zero) gz=0.;
  else                 gz=6.0*e*ixyz[2]/beta[1]/g/a/ll;
  gy=0.;gz=0.;  
  aj1=ixyz[1]*e/(1.+2*gy);
  integn1[0]=12./ll/dl*aj1;      integn1[1]= 6./ll*aj1; integn1[2]=-12./ll/dl*aj1; integn1[3]= 6./ll*aj1;
  integn2[0]=2.*(2.+gy)/dl*aj1;  integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gy)/dl*aj1;
  integn3[0]=12./ll/dl*aj1;      integn3[1]=-6./ll*aj1;
  integn4   =2.*(2.+gy)/dl*aj1;
  sm[2][2] = integn1[0];