beamel3d.cpp

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

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

/**
   function computes consistent mass %matrix of 3D beam element
   rotational DOFs are neglected
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param mm - mass %matrix of one element
   
   JK, 3.2.200
*/
void beamel3d::mass_matrix (long eid,long ri,long ci,matrix &mm)
{
  long i,transf;
  double a,l,gy,gz,xi,jac,rho,beta[2];
  ivector nodes(nne);
  vector vec(3),w(intordmm),gp(intordmm),x(nne),y(nne),z(nne);
  matrix n(3,ndofe),tran(12,12);
  
  fillm (0.0,mm);
  
  // element nodes
  Mt->give_elemnodes (eid,nodes);
  //  area of cross section
  Mc->give_area (eid,a);
  //  shear coefficients
  Mc->give_shearcoeff (eid,beta);
  //  density of the material
  Mc->give_densitye (eid,rho);
  //  vector defining local coordinate system
  Mc->give_vectorlcs (eid, vec);

  gy=beta[0];
  gz=beta[1];

  Mt->give_node_coord3d (x,y,z,eid);

  beam_transf_matrix (tran,l,vec,x,y,z,eid);
  
  //  coordinates and weights of integration points
  gauss_points (gp.a,w.a,intordmm);
  
  
  for (i=0;i<intordmm;i++){
    //  xi has to be between 0.0 and 1.0 while gp[i] are between -1.0 and 1.0
    xi=0.5*(1.0+gp[i]);
    //  matrix of approximation functions
    bf_matrix_transl (n,xi,l,gy,gz);
    jac=w[i]*l/2.0*rho*a;
    
    //  N^T rho N
    nnj (mm.a,n.a,jac,n.m,n.n);
    
  }
  
  //  transformation of mass matrix from local element coordinate
  //  system to the global coordinate system
  lgmatrixtransf (mm,tran);

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

}

/**
   function computes consistent mass %matrix of 3D beam element
   
   @param eid - element id
   @param mm - mass %matrix
   
   JK, 2.9.2006
*/
void beamel3d::res_mass_matrix (long eid,matrix &mm)
{
  mass_matrix (eid,0,0,mm);
}

/**
   @param eid - element id
   @param ri,ci - row and column indices
   @param ism - initial stress %matrix
   
*/
void beamel3d::initstr_matrix (long eid,long ri,long ci,matrix &ism)
{
  long i,i1,ii,transf;
  double nn,  e,g,a,dl,ll,ixyz[3],beta[2],gy,gz,g2,a4;
  ivector nodes(nne);
  vector x(nne),y(nne),z(nne),vec(3);
  matrix d(tncomp,tncomp),tran(12,12);

// nn....... Axial forces

  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);

//  Mc->give_vecbeam (eid,&vec);
  beam_transf_matrix (tran,dl, vec,x,y,z,eid);
  ll=dl*dl;
  
  fillm (0.0,ism);
  ii=Mt->elements[eid].ipp[ri][ci];
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  Mm->matstiff (d,ii);
  
  e=d[0][0];  g=d[1][1];


//  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;
  
  g2=gy*gy;
  a4=nn/(1.+2*gy)/(1.+2*gy);
  ism[2][2] = a4*4*(g2+gy+0.3)/dl;
  ism[4][4] = a4*(g2+gy+0.4)/3.*dl;
  ism[2][4] = - a4/10.;
  ism[2][8] = - a4*4*(g2+gy+0.3)/dl;
  ism[4][10]= a4*(-g2-gy-0.1)/3.*dl;
  ism[4][8] = a4/10.;
  ism[2][10]=- a4/10.;
//  direct s Iz
  g2=gz*gz;
  a4=nn/(1.+2.*gz)/(1.+2.*gz);
  ism[1][1] = a4*4.*(g2+gz+0.3)/dl;
  ism[5][5] = a4*(g2+gz+0.4)/3.*dl;
  ism[1][5] = a4/10.;
  ism[1][7] =- a4*4.*(g2+gz+0.3)/dl;
  ism[5][11]= a4*(-g2-gz-0.1)/3.*dl;
  ism[5][7] =- a4/10.;
  ism[1][11]= a4/10.;

  for (i=0;i<6;i++){
	  i1=i+6;
      ism[i1][i1]=ism[i][i];
      ism[i1][i]=ism[i][i1];
  }

  ism[4][2] = ism[2][4];
  ism[5][1] = ism[1][5];
  ism[7][5] = ism[5][7];
  ism[8][4] = ism[4][8];
  ism[11][1]= ism[1][11];
  ism[10][2]= ism[2][10];
  ism[7][11]=-ism[1][5];
  ism[11][7]=-ism[1][5];
  ism[8][10]=-ism[2][4];
  ism[10][8]=-ism[2][4];

  //  transformation of initial stress matrix from local element coordinate
  //  system to the global coordinate system
  lgmatrixtransf (ism,tran);

  //  transformation of initial stress matrix from global coordinate system
  //  to the local nodal systems
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tran);
    glmatrixtransf (ism,tran);
  }
}

/**
   function serves for computation of nodal forces and moments caused by body load
   
   @param eid - element id
   @param lm - load %matrix
   
   JK, 1.9.2006
*/
void beamel3d::load_matrix (long eid,matrix &lm)
{
  long i,ii,transf;
  double jac,a,ixyz[3],beta[2],e,g,gy,gz,l;
  ivector nodes (nne);
  vector vec(3),x(nne),y(nne),z(nne),w(intordmm),gp(intordmm),xl(2);
  matrix d(tncomp,tncomp),n(napfun,ndofe),tran(12,12);
  
  //  nodal coordinates
  Mt->give_node_coord3d (x,y,z,eid);
  //  coordinates of base vector z of the beam in global system
  Mc->give_vectorlcs (eid,vec);
  //  transformation matrix from local to global coordinates
  beam_transf_matrix (tran,l,vec,x,y,z,eid);
  //  area of cross section
  Mc->give_area (eid,a);
  //  shear coefficients for particular cross section
  Mc->give_shearcoeff (eid,beta);
  //  moments of inertia (I_x, I_y, I_z)
  Mc->give_mominer (eid,ixyz);
  //  number of integration point
  ii=Mt->elements[eid].ipp[0][0];
  //  stiffness matrix of material
  Mm->matstiff (d,ii);
  //  Young's modulus of leasticity
  e=d[0][0];
  //  shear modulus
  g=d[1][1];
  
  //  generalized shear modulus G_y
  if(beta[0]<Mp->zero) gy=0.0;
  else                 gy=6.0*e*ixyz[1]/beta[0]/g/a/l;
  //  generalized shear modulus G_z
  if(beta[1]<Mp->zero) gz=0.;
  else                 gz=6.0*e*ixyz[2]/beta[1]/g/a/l;
  
  
  //  local coordinates of nodes for evaluation of jacobian
  xl[0]=0.0;
  xl[1]=l;
  
  //  element node numbers
  Mt->give_elemnodes (eid,nodes);
  
  //  integration points assembling
  gauss_points (gp.a,w.a,intordmm);
  
  //  matrix cleaning
  fillm (0.0,lm);
  
  for (i=0;i<intordmm;i++){
    //  computation of jacobian
    jac_1d (jac,xl,gp[i]);
    
    //  assembling of matrix of approximation functions
    bf_matrix (n,(gp[i]+1.0)/2.0,l,gy,gz);
    
    jac*=a*w[i];
    
    //  multiplication N^T N jac
    nnj (lm.a,n.a,jac,n.m,n.n);
  }
  
  //  transformation of load matrix to the global system
  lgmatrixtransf (lm,tran);

  //  transformation of stiffness matrix to the nodal system
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    //  assembling of transformation matrix from global to nodal system
    transf_matrix (nodes,tmat);
    glmatrixtransf (lm,tmat);
  }
}


/**
   function stores nodal displacements and rotations in local coordinate system
   this is equivalent to functions computing strains
   nodal displacements and rotations are better quantities in case of beams
   
   @param eid - element id
   @param lcid - load case id
   
   JK, 20.2.2002
*/
void beamel3d::nodal_displ (long lcid,long eid)
{
  long i,j,ip,transf;
  double l;
  ivector nodes(nne),cn(ndofe);
  vector x(nne),y(nne),z(nne),vec(3),d(ndofe),dd(6),f(ndofe);
  matrix tm(3,3),tmat (ndofe,ndofe);
  
  //  code numbers on element
  Mt->give_code_numbers (eid,cn.a);
  
  //  nodal displacements and rotations
  eldispl (lcid,eid,d.a,cn.a,ndofe);
  
  //  element node numbers
  Mt->give_elemnodes (eid,nodes);

  //  transformation from nodal to global coordinate system
  transf = Mt->locsystems (nodes);
  if (transf>0){
    //  assembling of transformation matrix
    transf_matrix (nodes,tmat);
    //  transformation of nodal displacements and rotations
    //  transformation of nodal displacements and rotations
    lgvectortransf (f,d,tmat);
    copyv (f,d);
  }
  
  //  nodal coordinates
  Mt->give_node_coord3d (x,y,z,eid);
  
  //  vector of z-axis
  Mc->give_vectorlcs (eid,vec);
  
  //  assembling of transformation matrix from global to local coordinate system
  beam_transf_matrix (tmat,l,vec,x,y,z,eid);
  
  //  transformation from global problem coordinate system
  //  to local element coordinate system
  glvectortransf (d,f,tmat);
  copyv (f,d);
  
  j=6;
  for (i=0;i<6;i++){
    dd[i]=d[j];
    j++;
  }

  //  number of integration point
  ip=Mt->elements[eid].ipp[0][0];
  
  
  Mm->storestrain (lcid,ip,0,6,d);
  Mm->storestrain (lcid,ip+1,0,6,dd);
  
}

/**
   function computes nodal forces and moments expressed in local coodinate system
   this is equivalent to functions computing stresses
   nodal forces and moments are better quantities in case of beams

   x-axis goes through the beam
   z-axis is defined in cross section
   y-axis is defined by right hand orientation
   
   @param eid - element id
   @param lcid - load case id
   
   JK, 20.2.2002, revised 2.9.2006
*/
void beamel3d::nodal_forces (long lcid,long eid)
{
  long i,j,ii,transf;
  double l;
  ivector nodes(nne),cn(ndofe);
  vector x(nne),y(nne),z(nne),vec(3),d(ndofe),f(ndofe),ff(6);
  matrix sm(ndofe,ndofe),tmat(ndofe,ndofe);
  
  //  code numbers on element
  Mt->give_code_numbers (eid,cn.a);
  
  //  nodal displacements and rotations
  eldispl (lcid,eid,d.a,cn.a,ndofe);
  
  //  assembling of stiffness matrix
  stiffness_matrix (eid,0,0,sm);
  
  //  computation of nodal forces and moments
  mxv (sm,d,f);


  //  transformation of displacement vector from local nodal systems
  //  to global coordinate system
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    //  assembling of transformation matrix
    transf_matrix (nodes,tmat);
    //  transformation of nodal displacements and rotations
    lgvectortransf (d,f,tmat);
    copyv (d,f);
  }


  //  nodal coordinates
  Mt->give_node_coord3d (x,y,z,eid);
  
  //  vector of z-axis
  Mc->give_vectorlcs (eid,vec);
  
  //  assembling of transformation matrix from global to local coordinate system
  beam_transf_matrix (tmat,l,vec,x,y,z,eid);
  
  //  transformation from global problem coordinate system
  //  to local element coordinate system
  glvectortransf (f,d,tmat);
  copyv (d,f);
  
  j=6;
  for (i=0;i<6;i++){
    ff[i]=f[j];
    j++;
  }
  
  //  number of integration point
  ii=Mt->elements[eid].ipp[0][0];
  
  //  storage of nodal forces and moments
  //  nodal forces and moments are expressed in local coordinate system
  Mm->storestress (lcid,ii,0,6,f);
  Mm->storestress (lcid,ii+1,0,6,ff);

}

/**
   function computes internal forces
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - vector of internal forces
   
   20.12.2002
*/
void beamel3d::res_internal_forces (long lcid,long eid,vector &ifor)
{
  internal_forces (lcid, eid, 0, 0, ifor);
}

void beamel3d::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long i,integr=2,ipp,k,transf;
  double dl,ll,a,ixyz[3],beta[2],gy,gz,e,g,s,dj,   m;
  ivector cn(ndofe),nodes(nne);
  vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),fx(tncomp),w(integr),gp(integr);
  matrix n(tncomp,ndofe),d(tncomp,tncomp),tran(12,12),r(tncomp,2);
  
  // r  .... end forces last time
  // m  .... dl(i)/dl(i-1)
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,rl.a,cn.a,ndofe);
  Mc->give_vectorlcs (eid, vec);
  beam_transf_matrix (tran,dl, vec,x,y,z,eid);
  ll=dl*dl;
  // r  .... end forces last time
  // m  .... dl(i)/dl(i-1)
  m=dl;
  fillm (0.0,r);
  
  //  transformation of displacement vector
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    //locglobtransf (rg,rl,tmat);
    lgvectortransf (rg,rl,tmat);
  }
  else{
    copyv (rl,rg);
  }
  // *** transf. from gcs to lcs
  //  globloctransf (rg,rl,tran);
  
  //  stiffness_matrix (eid,ri,ci,sm);
  //  mxv (sm,r,contr);
  
  ipp=Mt->elements[eid].ipp[0][0];
  Mm->matstiff (d,ipp);
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  
  e=d[0][0];  g=d[1][1];
  
  //  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;
  gauss_points (gp.a,w.a,integr);
  for (i=0;i<integr;i++){
    //	 s=x/dl
    s=(gp[i]+1.)/2.; // gaus point in (-1,1) function on beam (0,1)
    dj=w[i]/2.*dl;   // weight divide by 2
    geom_matrix (n,s, dl, gy, gz);
    fx[0]=( n[0][0]*rl[0]+n[0][6]*rl[6] )*a*e;
    fx[3]=( n[3][3]*rl[3]+n[3][9]*rl[9] )*ixyz[0]*g;
    fx[1]=( n[1][1]*rl[1]+n[1][5]*rl[5]+n[1][7]*rl[7]+n[1][11]*rl[11] )*beta[1]*a*g;
    fx[5]=( n[5][1]*rl[1]+n[5][5]*rl[5]+n[5][7]*rl[7]+n[5][11]*rl[11] )*ixyz[2]*e;
    fx[2]=( n[2][2]*rl[2]+n[2][4]*rl[4]+n[2][8]*rl[8]+n[2][10]*rl[10] )*beta[0]*a*g;
    fx[4]=( n[4][2]*rl[2]+n[4][4]*rl[4]+n[4][8]*rl[8]+n[4][10]*rl[10] )*ixyz[1]*e;
    
    //  change of the II piola-Kirchhoff to Cauchyho real stress (through scale m)
    f[0]=f[0]+( r[0][0]*(1-s) + r[0][1]*s+fx[0] )*n[0][0]*dj;
    f[6]=f[6]+( r[0][0]*(1-s) + r[0][1]*s-fx[0] )*n[0][6]*dj;
    f[3]=f[3]+( r[3][0]*(1-s) + r[3][1]*s+fx[3] )*n[3][3]*dj;
    f[9]=f[9]+( r[3][0]*(1-s) + r[3][1]*s-fx[3] )*n[3][9]*dj;
    // Vy, Mz
    //     f[1]=f[1]+( (r[1][0]*(1-s)+r[1][1]*s)*n[1][1] + ((r[5][0]*(1-s)+r[5][1]*s)*dl/m)*n[1][5] )*dj;
    f[1]=f[1]  +(fx[1]*n[1][1] + fx[5]*n[1][5])*dj;
    
    //     f[7]=f[7]+( (r[1][0]*(1-s)+r[1][1]*s)*n[1][1] + ((r[5][0]*(1-s)+r[5][1]*s)*dl/m)*n[5][7] )*dj;
    f[7]=f[7]  +(fx[1]*n[1][7] + fx[5]*n[5][7])*dj;
    
    //     f[5]=f[5]+( ((r[5][0]*(1-s)+r[5][1]*s)*dl/m)*n[5][5] + (r[1][0]*(1-s)+r[1][1]*s)*n[1][5] )*dj;
    f[5]=f[5]  +(fx[5]*n[5][5] + fx[1]*n[1][5])*dj;
    
    //     f[11]=f[11]+( ((r[5][0]*(1-s)+r[5][1]*s)*dl/m)*n[11][5] + (r[1][0]*(1-s)+r[1][1]*s)*n[11][1] )*dj;
    f[11]=f[11]+(fx[5]*n[5][11] + fx[1]*n[1][11])*dj;
    // Vz, My
    //     f[2]=f[2]+( (r[2][0]*(1-s)+r[2][1]*s)*n[2][2] + ((r[4][0]*(1-s)+r[4][1]*s)*dl/m)*n[4][2] )*dj;
    f[2]=f[2]  +(fx[2]*n[2][2] + fx[4]*n[4][2])*dj;
    
    //     f[8]=f[8]+( (r[2][0]*(1-s)+r[2][1]*s)*n[2][8] + ((r[4][0]*(1-s)+r[4][1]*s)*dl/m)*n[4][8] )*dj;
    f[8]=f[8]  +(fx[2]*n[2][8] + fx[4]*n[4][8])*dj;
    
    //     f[4]=f[4]+( ((r[4][0]*(1-s)+r[4][1]*s)*dl/m)*n[4][4] + (r[2][0]*(1-s)+r[2][1]*s)*n[2][4] )*dj;
    f[4]=f[4]  +(fx[4]*n[4][4] + fx[2]*n[2][4] )*dj;
    
    //     f[10]=f[10]+( ((r[4][0]*(1-s)+r[4][1]*s)*dl/m)*n[4][10] + (r[2][0]*(1-s)+r[2][1]*s)*n[2][10] )*dj;
    f[10]=f[10]+(fx[4]*n[4][10] + fx[2]*n[2][10])*dj;
  }
  
  for (k=0;k<ndofe;k++){
    ifor[k]+=f[k];
    ifor[k]=0.0;                     //nekonsoliduje 
  }
  //  fprintf (Out,"\n\n R prvek beam cislo %ld\n",eid);
  //  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",rl[k]); }
  //  fprintf (Out,"\n\n Fint prvek cislo %ld\n",eid);
  //  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",f[k]);}
}