beamel2d.cpp

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

   influence of inertial forces from rotations is accounted
   
   @param eid - number of element
   @param mm - mass matrix
   
   JK, 3.2.2006
*/
void beamel2d::res_mass_matrix (long eid,matrix &mm)
{
  //mass_matrix (eid,0,0,mm);
  mass_matrix_expl (eid,0,0,mm);
}

/**
   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_expl (long eid,long ri,long ci,matrix &mm)
{
  long ipid,transf;
  double e,g,shearcoeff,j,iy,a,l,kappa,rho;
  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::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);
  //  density of the material
  Mc->give_densitye (eid,rho);
  //  moment of inertia of element cross section
  Mc->give_mominer (eid,&iy);
  //  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;
  
  j = (1.0+2.0*kappa)*(1.0+2.0*kappa);
  
  
  mm[0][0]=rho*a*l/3.0;
  mm[0][3]=rho*a*l/6.0;
  
  mm[1][1]=rho*a*l*(13.0/35.0+7.0/5.0*kappa+4.0/3.0*kappa*kappa)/j;
  mm[1][2]=rho*a*l*l*(-11.0/210.0-11.0/60*kappa-1.0/6.0*kappa*kappa)/j;
  mm[1][4]=rho*a*l*(9.0/70+3.0/5.0*kappa+2.0/3.0*kappa*kappa)/j;
  mm[1][5]=rho*a*l*l*(13.0/420+3.0/20.0*kappa+1.0/6.0*kappa*kappa)/j;
  
  mm[2][1]=mm[1][2];
  mm[2][2]=rho*a*l*l*l*(1.0/105.0+1.0/30.0*kappa+1.0/30.0*kappa*kappa)/j;
  mm[2][4]=rho*a*l*l*(-13.0/420.0-3.0/20.0*kappa-1.0/6.0*kappa*kappa)/j;
  mm[2][5]=rho*a*l*l*l*(-1.0/140.0-1.0/30.0*kappa-1.0/30.0*kappa*kappa)/j;
  
  mm[3][0]=mm[0][3];
  mm[3][3]=mm[0][0];
  
  mm[4][1]=mm[1][4];
  mm[4][2]=mm[2][4];
  mm[4][4]=rho*a*l*(13.0/35.0+7.0/5.0*kappa+4.0/3.0*kappa*kappa)/j;
  mm[4][5]=rho*a*l*l*(11.0/210.0+11.0/60.0*kappa+1.0/6.0*kappa*kappa)/j;
  
  mm[5][1]=mm[1][5];
  mm[5][2]=mm[2][5];
  mm[5][4]=mm[4][5];
  mm[5][5]=rho*a*l*l*l*(1.0/105+1.0/30.0*kappa+1.0/30.0*kappa*kappa)/j;

  
  /*
  mm[0][0]=rho*a*l/2.0;
  mm[1][1]=rho*a*l/24.0*12.0;
  mm[2][2]=rho*a*l*l*l/24.0;
  mm[3][3]=rho*a*l/2.0;
  mm[4][4]=rho*a*l/24.0*12.0;
  mm[5][5]=rho*a*l*l*l/24.0;
  */
  

  //  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 geometric-stiffness %matrix
   notation based on R. Clough, J. Penzien: Dynamic of Structures,
   McGraw-Hill, 2nd edition, 1993, pages 195-196
   
   notation initial stress %matrix is also used
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param ism - inital stress %matrix
   
   JK
*/
void beamel2d::initstr_matrix (long eid,long ri,long ci,matrix &ism)
{
  long i,ipid,transf;
  double xi,ww,jac,e,g,shearcoeff,iy,a,l,kappa;
  ivector nodes(nne);
  vector x(nne),z(nne),w(intordism),gp(intordism);
  matrix n(1,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,intordism);

  //  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::initstr_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  fillm (0.0,ism);
  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);
  //  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<intordism;i++){
    xi=(1.0+gp[i])/2.0;  ww=w[i];
    dbf_matrix (n,xi,l,kappa);
    jac=a*l/2.0*ww;
    nnj (ism.a,n.a,jac,n.m,n.n);
  }

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

  //  transformation of initial stress 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 (ism,tmat);
  }
}


/**
   function computes consistent geometric-stiffness %matrix
   notation based on R. Clough, J. Penzien: Dynamic of Structures,
   McGraw-Hill, 2nd edition, 1993, pages 195-196
   
   notation initial stress %matrix is also used
   
   initial stress %matrix is formulated explicitly

   @param eid - number of element
   @param ri,ci - row and column indices
   @param ism - inital stress %matrix
   
   JK
*/
void beamel2d::initstr_matrix_expl (long eid,long ri,long ci,matrix &ism)
{
  long ipid,transf;
  double e,g,shearcoeff,iy,a,l,kappa,denom;
  ivector nodes(nne);
  vector x(nne),z(nne),w(intordism),gp(intordism);
  matrix 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,intordism);

  //  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::initstr_matrix (file %s, line %d).\n",__FILE__,__LINE__);
  }
  
  fillm (0.0,ism);
  ipid=Mt->elements[eid].ipp[ri][ci];
  //  area of element cross section
  Mc->give_area (ipid,a);
  //  shear coefficient
  Mc->give_shearcoeff (eid,&shearcoeff);
  //  moment of inertia of element cross section
  Mc->give_mominer (eid,&iy);
  //  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;
  
  denom=(1.0+2.0*kappa)*(1.0+2.0*kappa);
  
  fillm (0.0,ism);
  
  ism[1][1]=(6.0/5.0+4.0*kappa+4.0*kappa*kappa)/denom/l;
  ism[1][2]=-1.0/10.0/denom;
  ism[1][4]=(-6.0/5.0-4.0*kappa-4.0*kappa*kappa)/denom/l;
  ism[1][5]=-1.0/10.0/denom;
  
  ism[2][1]=ism[1][2];
  ism[2][2]=(2.0/15.0+kappa/3.0+kappa*kappa/3.0)*l/denom;
  ism[2][4]=1.0/10.0/denom;
  ism[2][5]=(-1.0/30.0-kappa/3.0-kappa*kappa/3.0)*l/denom;
  
  ism[4][1]=ism[1][4];
  ism[4][2]=ism[2][4];
  ism[4][4]=(6.0/5.0+4.0*kappa+4.0*kappa)/denom/l;
  ism[4][5]=1.0/10.0/denom;
  
  ism[5][1]=ism[1][5];
  ism[5][2]=ism[2][5];
  ism[5][4]=ism[4][5];
  ism[5][5]=(2.0/15.0+kappa/3.0+kappa*kappa/3.0)*l/denom;
  
  
  //  transformation of initial stress matrix from local element coordinate
  //  system to the global coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  lgmatrixtransf (ism,tmat);

  //  transformation of initial stress 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 (ism,tmat);
  }
}



/**
   function stores end displacements and rotations to integration points
   
   @param lcid - load case id
   @param eid - element id
   
   JK, 24.5.2006
*/
void beamel2d::nodal_displ (long lcid,long eid)
{
  long i,j,ipid,transf;
  ivector cn(ndofe),nodes(nne);
  vector x(nne),z(nne),d(ndofe),dd(3),f(ndofe);
  matrix tmat(ndofe,ndofe);
  
  //  node coordinates in x-z plane
  Mt->give_node_coord2dxz (x,z,eid);
  //  code numbers on element
  Mt->give_code_numbers (eid,cn.a);
  //  nodal displacements and rotations
  eldispl (lcid,eid,d.a,cn.a,ndofe);
  
  //  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 (f,d,tmat);
    copyv (f,d);
  }
  
  //  assembling of transformation matrix from global problem coordinate system
  //  to local element coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  //  transformation from global problem coordinate system
  //  to local element coordinate system
  glvectortransf (d,f,tmat);
  copyv (f,d);
  
  j=3;
  for (i=0;i<3;i++){
    dd[i]=d[j];
    j++;
  }
  
  //  number of integration point
  ipid=Mt->elements[eid].ipp[0][0];
  
  //  storage of components
  Mm->storestrain (lcid,ipid,0,3,d);
  Mm->storestrain (lcid,ipid+1,0,3,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

   @param lcid - load case id
   @param eid - element id
   
   JK, 20.2.2002, revised 2.9.2006
*/
void beamel2d::nodal_forces (long lcid,long eid)
{
  long i,j,ipid,transf;
  ivector nodes(nne),cn(ndofe);
  vector x(nne),z(nne),d(ndofe),f(ndofe),ff(3),ifor(ndofe);
  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_expl (eid,0,0,sm);
  
  //  computation of nodal forces and moments
  mxv (sm,d,f);
  
  
  

  if (Mp->tprob == forced_dynamics || Mp->tprob == eigen_dynamics){
    
    //  mass matrix
    mass_matrix_expl (eid,0,0,sm);
    //  contribution from inertial forces
    mxv (sm,d,ifor);
    cmulv (Lsrs->w[lcid],ifor,ifor);
    
    f[0]-=ifor[0];
    f[1]-=ifor[1];
    f[2]-=ifor[2];

    f[3]-=ifor[3];
    f[4]-=ifor[4];
    f[5]-=ifor[5];
  }


  //  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);
  }
  
  //  node coordinates in x-z plane
  Mt->give_node_coord2dxz (x,z,eid);
  //  assembling of transformation matrix from global problem coordinate system
  //  to local element coordinate system
  beam_transf_matrix (eid,x,z,tmat);
  //  transformation from global problem coordinate system
  //  to local element coordinate system
  glvectortransf (f,d,tmat);
  copyv (d,f);
  
  
  j=3;
  for (i=0;i<3;i++){
    ff[i]=f[j];
    j++;
  }
  
  //  number of integration point
  ipid=Mt->elements[eid].ipp[0][0];
  
  //  storage of nodal forces and moments
  //  nodal forces and moments are expressed in local element coordinate system
  Mm->storestress (lcid,ipid,0,3,f);
  Mm->storestress (lcid,ipid+1,0,3,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
   
   JK, 12.8.2001
*/
void beamel2d::internal_forces (long lcid,long eid,vector &ifor)
{
  long i,j,ipid;
  ivector cn(ndofe),nodes(nne);
  vector iforl(ndofe),r(ndofe),contr(ndofe),ff(3);
  matrix sm(ndofe,ndofe),tmat(ndofe,ndofe);
  
  fillv (0.0,ifor);

  //  code numbers on element
  Mt->give_code_numbers (eid,cn.a);
  //  nodal displacements and rotations
  eldispl (lcid,eid,r.a,cn.a,ndofe);
  //  assembling of stiffness matrix
  stiffness_matrix_expl (eid,0,0,sm);
  //  computation of nodal forces and moments
  mxv (sm,r,contr);
  
  addv (contr,ifor,ifor);

  j=3;
  for (i=0;i<3;i++){
    ff[i]=ifor[j];
    j++;
  }
  
  //  number of integration point
  ipid=Mt->elements[eid].ipp[0][0];
  
  //  storage of nodal forces and moments
  //  nodal forces and moments are expressed in global problem or local nodal coordinate system
  Mm->storestress (lcid,ipid,0,3,ifor);
  Mm->storestress (lcid,ipid+1,0,3,ff);
}

/**
   function computes internal forces

   @param lcid - load case id
   @param eid - element id
   @param ifor - vector of internal forces
   
   JK, 12.8.2001
*/
void beamel2d::res_internal_forces (long lcid,long eid,vector &ifor)
{
  internal_forces (lcid,eid,ifor);
}

/**
   function defines meaning of DOFs at nodes
   
   @param eid - element id
   
   1.2.2005, JK
*/
void beamel2d::define_meaning (long eid)
{
  ivector cn(ndofe),nod(nne);
  
  Mt->give_elemnodes (eid,nod);
  Mt->give_code_numbers (eid,cn.a);

  //  displacement in x direction
  if (cn[0]>0)  Mt->nodes[nod[0]].meaning[0]=1;
  //  displacement in z direction
  if (cn[1]>0)  Mt->nodes[nod[0]].meaning[1]=3;
  //  displacement in x direction
  if (cn[3]>0)  Mt->nodes[nod[1]].meaning[0]=1;
  //  displacement in z direction
  if (cn[4]>0)  Mt->nodes[nod[1]].meaning[1]=3;

}