beamgen3d.cpp

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

  vector x(nne),y(nne),z(nne),vec(3),d(ndofe),f(ndofe);
  matrix tm(3,3);
  
  //  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){
    matrix tmat (ndofe,ndofe);
    //  assembling of transformation matrix
    transf_matrix (nodes,tmat);
    //  transformation of nodal displacements and rotations
    mtxv (tmat,d,f);
    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 (tm,l,vec,x,y,z,eid);
  
  //  transformation from global to local coordinate system
//   !!!zmenit   glvectortransf3dblock (d,tm);

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

/**
   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 beamgen3d::nodal_forces (long eid,long lcid)
{
  long ii,j,transf;
  double l;
  ivector nodes(nne),cn(ndofe);
  vector x(nne),y(nne),z(nne),vec(3),d(ndofe),f(ndofe);
  matrix sm(ndofe,ndofe),tm(3,3);
  
  //  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);
 	fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
	fprintf (Out," \n u, v, w, fx, fy, fz, o  \n ");
	for (j=0;j<ndofe;j++){
		fprintf (Out," %15.5e",d[j]);
    }
	fprintf (Out," \n N,  Vy,  Vz, Mx, My, Mz, B,\n ");
	for (j=0;j<ndofe;j++){
		fprintf (Out," %15.5e",f[j]);
    }
 
  //  element node numbers
  Mt->give_elemnodes (eid,nodes);

  //  transformation from nodal to global coordinate system
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    //  assembling of transformation matrix
    transf_matrix (nodes,tmat);
    //  transformation of nodal forces and moments
    mtxv (tmat,f,d);
    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 (tm,l,vec,x,y,z,eid);

  internal_forces (lcid, eid, 0, 0, f);

  //  transformation from global to local coordinate system
  //glvectortransf3dblock (f,tm);
  
  //  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,6,6,f);
  
  
}

/**
   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
   
   PF 20.9.2006
*/
void beamgen3d::res_internal_forces (long lcid,long eid,vector &ifor)
{
   internal_forces1 (lcid, eid, 0, 0, ifor);
}
// *******************************************************
//     TORSION - Forces
void beamgen3d::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long i,j,integr=2,ipp,transf;
  double dl,ll,a,ixyz[3],ioyz[3],beta[2],gy,gz,kakr,e,g;
  double eh,ehx,ajy,ajz,ajk,xs, s,c;
  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);
  vector dderno(4),dderny(4),ddernz(4);
  matrix n(tncomp,ndofe),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
  matrix bb(tncomp,ndofe), ck(4,4);
  
// 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)
  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);
  // 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]=10.0; ixyz[0]=1.; ixyz[1]=10.; ixyz[2]=1.;

  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;
    gy=0.; gz=0.; 
	kakr=0.;

  eh=pow(g*ixyz[0]/e/ioyz[0],0.5); 
//  eh=1.; 
  eh=eh*dl; 
  s=sinh(eh); c=cosh(eh);
  eh=eh/dl; 
  ck_matrix (ck,  s,c,eh,dl);
  xs=0.0;
  for (i=0;i<2;i++){

    ehx=eh*xs*dl;
    s=sinh(ehx); c=cosh(ehx);


//    Vnitrni sily   N, My, Mz, B , Mvaz,  Mvol
    bb[0][0]=-e*a/dl;
    bb[0][7]= e*a/dl;
//        BB(i,4/7/11/14)=  e*pru(10+i)*eh*eh*(CK(i,3)*C+CK(i,4)*S)
    ajy= (1.+2.*gy); ajz= (1.+2.*gz); ajk= (1.+2.*kakr);
    dderno[0]=(-6.+12.*xs)/ll;  dderno[1]=(-2.*(2.+kakr)+6.*xs)/dl; dderno[2]=( 6.-12.*xs)/ll;  dderno[3]=(-2.*(1.-kakr)+6.*xs)/dl;
//    dderno[0]=eh*eh*(ck[0][2]*s+ck[0][3]*c);  dderno[1]=eh*eh*(ck[1][2]*s+ck[1][3]*c); dderno[2]=eh*eh*(ck[2][2]*s+ck[2][3]*c);  dderno[3]=eh*eh*(ck[3][2]*s+ck[3][3]*c);
    dderny[0]=(-6.+12.*xs)/ll;  dderny[1]=(-2.*(2.+gy)+6.*xs)/dl;   dderny[2]=( 6.-12.*xs)/ll;  dderny[3]=(-2.*(1.-gy)+6.*xs)/dl;
//  ex= - o O'' - y v'' - z w''
//  My = Int. z E ex = - E Iy w''- E Ioy O''=  E Iy fiy'- E Ioy O''
	bb[1][2]=  -e*ixyz[1]/ajy*dderny[0];
    bb[1][3]=  -e*ioyz[1]/ajk*dderno[0];
    bb[1][4]=   e*ixyz[1]/ajy*dderny[1];
    bb[1][6]=  -e*ioyz[1]/ajk*dderno[1];
    bb[1][9]=  -e*ixyz[1]/ajy*dderny[2];
    bb[1][10]= -e*ioyz[1]/ajk*dderno[2];
    bb[1][11]=  e*ixyz[1]/ajy*dderny[3];
    bb[1][13]= -e*ioyz[1]/ajk*dderno[3];
// Qz
    bb[5][2]=  e*ixyz[1]/ajy*12./ll/dl;
    bb[5][3]=  e*ioyz[1]/ajk*12./ll/dl;
	bb[5][4]= -e*ixyz[1]/ajy*6./ll;
    bb[5][6]=  e*ioyz[1]/ajk*6./ll;
    bb[5][9]=  e*ixyz[1]/ajy*(-12.)/ll/dl;
    bb[5][10]= e*ioyz[1]/ajk*(-12.)/ll/dl;
    bb[5][11]=-e*ixyz[1]/ajy*6./ll;
    bb[5][13]= e*ioyz[1]/ajk*6./ll;

// Mz=  -Int. y E ex = E Iz v''+ E Ioz O''=  E Iz fiz'+ E Ioy O''
    ddernz[0]=(-6.+12.*xs)/ll;  ddernz[1]=(-2.*(2.+gz)+6.*xs)/dl; ddernz[2]=( 6.-12.*xs)/ll;  ddernz[3]=(-2.*(1.-gz)+6.*xs)/dl;
    bb[2][1]=  e*ixyz[2]/ajz*ddernz[0];
    bb[2][3]=  e*ioyz[2]/ajk*dderno[0];
    bb[2][5]=  e*ixyz[2]/ajz*ddernz[1];
    bb[2][6]=  e*ioyz[2]/ajk*dderno[1];
    bb[2][ 8]= e*ixyz[2]/ajz*ddernz[2];
    bb[2][10]= e*ioyz[2]/ajk*dderno[2];
    bb[2][12]= e*ixyz[2]/ajz*ddernz[3];
    bb[2][13]= e*ioyz[2]/ajk*dderno[3];
// Qy
    bb[4][1]=  e*ixyz[2]/ajz*12./ll/dl;
    bb[4][3]=  e*ioyz[2]/ajk*12./ll/dl;
    bb[4][5]=  e*ixyz[2]/ajz*6./ll;
    bb[4][6]=  e*ioyz[2]/ajk*6./ll;
    bb[4][8] = e*ixyz[2]/ajz*(-12.)/ll/dl;
    bb[4][10]= e*ioyz[2]/ajk*(-12.)/ll/dl;
    bb[4][12]= e*ixyz[2]/ajz*6./ll;
    bb[4][13]= e*ioyz[2]/ajk*6./ll;

// B =  Int. o E ex = -E Io O'' - E Ioz v'' - E Ioy w''
    bb[3][1]= -e*ioyz[2]/ajz*ddernz[0];
//    bb[3][1]= -e*ioyz[2]/ajz*dderno[0];
    bb[3][2]= -e*ioyz[1]/ajy*dderny[0];
//    bb[3][2]= -e*ioyz[1]/ajy*dderno[0];
    bb[3][3]= -e*ioyz[0]*eh*eh*(ck[0][2]*c+ck[0][3]*s);
    bb[3][4]=  e*ioyz[1]/ajy*dderny[1];
//    bb[3][4]=  e*ioyz[1]/ajy*dderno[1];
    bb[3][5]= -e*ioyz[2]/ajz*ddernz[1];
//    bb[3][5]= -e*ioyz[2]/ajz*dderno[1];
    bb[3][6]= -e*ioyz[0]*eh*eh*(ck[1][2]*c+ck[1][3]*s);
    bb[3][8]= -e*ioyz[2]/ajz*ddernz[2];
//    bb[3][8]= -e*ioyz[2]/ajz*dderno[2];
    bb[3][9]= -e*ioyz[1]/ajy*dderny[2];
//    bb[3][9]= -e*ioyz[1]/ajy*dderno[2];
    bb[3][10]=-e*ioyz[0]*eh*eh*(ck[2][2]*c+ck[2][3]*s);
    bb[3][11]= e*ioyz[1]/ajy*dderny[3];
//    bb[3][11]= e*ioyz[1]/ajy*dderno[3];
    bb[3][12]=-e*ioyz[2]/ajz*ddernz[3];
//    bb[3][12]=-e*ioyz[2]/ajz*dderno[3];
    bb[3][13]=-e*ioyz[0]*eh*eh*(ck[3][2]*c+ck[3][3]*s);
   
//third deriv. +, +, -, + 
// Mvazane
    bb[6][1]= -e*ioyz[2]/ajz*12./ll/dl;
    bb[6][2]= -e*ioyz[1]/ajy*12./ll/dl;
    bb[6][3]= -e*ioyz[0]*eh*eh*eh*(ck[0][2]*s+ck[0][3]*c);
    bb[6][4]=  e*ioyz[1]/ajy*6./ll;
    bb[6][5]= -e*ioyz[2]/ajz*6./ll;
    bb[6][6]= -e*ioyz[0]*eh*eh*eh*(ck[1][2]*s+ck[1][3]*c);
    bb[6][8]= -e*ioyz[2]/ajz*(-12.)/ll/dl;
    bb[6][9]= -e*ioyz[1]/ajy*(-12.)/ll/dl;
    bb[6][10]=-e*ioyz[0]*eh*eh*eh*(ck[2][2]*s+ck[2][3]*c);
    bb[6][11]= e*ioyz[1]/ajy*6./ll;
    bb[6][12]=-e*ioyz[2]/ajz*6./ll;
    bb[6][13]=-e*ioyz[0]*eh*eh*eh*(ck[3][2]*s+ck[3][3]*c);
// Mvolne
    bb[7][3]=  g*ixyz[0]*(ck[0][1]+eh*(ck[0][2]*s+ck[0][3]*c));
    bb[7][6]=  g*ixyz[0]*(ck[1][1]+eh*(ck[1][2]*s+ck[1][3]*c));
    bb[7][10]= g*ixyz[0]*(ck[2][1]+eh*(ck[2][2]*s+ck[2][3]*c));
    bb[7][13]= g*ixyz[0]*(ck[3][1]+eh*(ck[3][2]*s+ck[3][3]*c));


    mxv (bb,rl,fx);

    xs=1.0;
	fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
	fprintf (Out," \n u, v, w, fx, fy, fz, o  \n ");
	for (j=0;j<ndofe;j++){
		fprintf (Out," %15.5e",rl[j]);
    }
	fprintf (Out," \n N, My, Mz, B, Vy,  Vz, Mo  Mt\n ");
	for (j=0;j<tncomp;j++){
		fprintf (Out," %15.5e",fx[j]);
    }
  }

}
// *******************************************************
//     TORSION - Forces

void beamgen3d::internal_forces1 (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long i,j,integr=2,ipp,transf;
  double dl,ll,a,ixyz[3],ioyz[3],beta[2],gy,gz,e,g;
  double  ajy,ajz,ajk, xs, kakr ;
  ivector cn(ndofe),nodes(nne);
  vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),w(integr),gp(integr), ddern(4),dderno(4);
  vector fx(4);
  matrix n(tncomp,ndofe),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
  matrix bb(4,ndofe);
  
// 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)
  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);
  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]=10.0; ixyz[0]=1.; ixyz[1]=10.; ixyz[2]=1.;

  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;

  gy=0.; gz=0.; kakr=0.;

  xs=0.0;
  for (i=0;i<2;i++){
//    Vnitrni sily   N, My, B, Mz, B
    bb[0][0]=-e*a/dl;
    bb[0][7]= e*a/dl;

    ajy= (1.+2.*gy); ajz= (1.+2.*gz); ajk= (1.+2.*kakr);
    dderno[0]=(-6.+12.*xs)/ll; dderno[1]=(-2.*(2.+kakr)+6.*xs)/dl; dderno[2]=( 6.-12.*xs)/ll;  dderno[3]=(-2.*(1.-kakr)+6.*xs)/dl;
    ddern[0]=(-6.+12.*xs)/ll;  ddern[1]=(-2.*(2.+gy)+6.*xs)/dl;    ddern[2]=( 6.-12.*xs)/ll;   ddern[3]=(-2.*(1.-gy)+6.*xs)/dl;
//My
	bb[1][2]=  -e*ixyz[1]/ajy*ddern[0];
    bb[1][3]=  -e*ioyz[1]/ajk*dderno[0];
    bb[1][4]=   e*ixyz[1]/ajy*ddern[1];
    bb[1][6]=  -e*ioyz[1]/ajk*dderno[1];
    bb[1][9]=  -e*ixyz[1]/ajy*ddern[2];
    bb[1][10]= -e*ioyz[1]/ajk*dderno[2];
    bb[1][11]=  e*ixyz[1]/ajy*ddern[3];
    bb[1][13]= -e*ioyz[1]/ajk*dderno[3];
// B
	bb[3][2]=  -e*ioyz[1]/ajy*ddern[0];
    bb[3][3]=  -e*ioyz[0]/ajk*dderno[0];
    bb[3][4]=   e*ioyz[1]/ajy*ddern[1];
    bb[3][6]=  -e*ioyz[0]/ajk*dderno[1];
    bb[3][9]=  -e*ioyz[1]/ajy*ddern[2];
    bb[3][10]= -e*ioyz[0]/ajk*dderno[2];
    bb[3][11]=  e*ioyz[1]/ajy*ddern[3];
    bb[3][13]= -e*ioyz[0]/ajk*dderno[3];

//Mz
    ddern[0]=(-6.+12.*xs)/ll;  ddern[1]=(-2.*(2.+gz)+6.*xs)/dl; ddern[2]=( 6.-12.*xs)/ll;  ddern[3]=(-2.*(1.-gz)+6.*xs)/dl;
    bb[2][1]=  e*ixyz[2]/ajz*ddern[0];
    bb[2][3]=  e*ioyz[2]/ajk*dderno[0];
    bb[2][5]=  e*ixyz[2]/ajz*ddern[1];
    bb[2][6]=  e*ioyz[2]/ajk*dderno[1];
    bb[2][ 8]= e*ixyz[2]/ajz*ddern[2];
    bb[2][10]= e*ioyz[2]/ajk*dderno[2];
    bb[2][12]= e*ixyz[2]/ajz*ddern[3];
    bb[2][13]= e*ioyz[2]/ajk*dderno[3];
// B
    bb[3][1]=  -e*ioyz[2]/ajz*ddern[0];
    bb[3][5]=  -e*ioyz[2]/ajz*ddern[1];
    bb[3][8]=  -e*ioyz[2]/ajz*ddern[2];
    bb[3][12]= -e*ioyz[2]/ajz*ddern[3];

    mxv (bb,rl,fx);

    xs=1.0;
	fprintf (Out,"\n\n MT prvek beam cislo %ld",lcid);
	fprintf (Out," \n u, v, w, fx, fy, fz, o  \n ");
	for (j=0;j<ndofe;j++){
		fprintf (Out," %15.5e",rl[j]);
    }
	fprintf (Out," \n N, My, Mz, B\n ");
	for (j=0;j<4;j++){
		fprintf (Out," %15.5e",fx[j]);
    }
  }
}