beamgen3d.cpp

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

  sm[2][4] = -integn1[1];
  sm[2][9] = integn1[2];
  sm[2][11]= -integn1[3];
  sm[4][4] = integn2[0];
  sm[4][9] = integn2[1];
  sm[4][11]= -integn2[2];
  sm[9][9] = integn3[0];
  sm[9][11]= -integn3[1];
  sm[11][11]= integn4;
//  direct s Iz   v, fz
  aj1=ixyz[2]*e/(1.+2.*gz);
  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.+gz)/dl*aj1;  integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gz)/dl*aj1;
  integn3[0]=12./ll/dl*aj1;      integn3[1]=-6./ll*aj1;
  integn4   =2.*(2.+gz)/dl*aj1;
  sm[1][1] = integn1[0];
  sm[1][5] = integn1[1];
  sm[1][8] = integn1[2];
  sm[1][12]= integn1[3];
  sm[5][5] = integn2[0];
  sm[5][8] = integn2[1];
  sm[5][12]= integn2[2];
  sm[8][8] = integn3[0];
  sm[8][12]= integn3[1];
  sm[12][12]= integn4;
  
  stiffness_matrixtor (sm, dl,e,g,gy,gz,ixyz,ioyz);

  for (i=0;i<ndofe;i++){for (int j=i+1;j<ndofe;j++){sm[j][i]=sm[i][j];}}

//  transformation of stiffness matrix to the global system
// !!!zmenit   lgmatrixtransf3dblock (sm,tran);

  //  transformation of stiffness matrix to the nodal system
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
  
  fprintf (Out,"\n\n MT prvek beam cislo %ld",eid);
  for (i=0;i<ndofe;i++){
     fprintf (Out,"\n %4ld",i);
     for (int j=0;j<ndofe;j++){
		fprintf (Out," %15.5e",sm[i][j]);
     }
  }

  delete [] ixyz;
  delete [] ioyz;
  delete [] beta;
}

/**
   function computes consistent mass %matrix of 3D 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
   
   PF, 20.9.2006
*/
// *******************************************************
//    TORSION -   exact solver

  void beamgen3d::stiffness_matrixtor (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz)
{
  int i,i1,i2,j;
  double  eh,eh1,s,c,c1,s1,s2,c2,cs,pom1,pom2,pom3,xc,xs;
  vector sioyz(2),sm1(16),sm2(16);
  matrix ck(4,4);
  

//     ioyz[0] - I omega - vys. moment setrvac,   ioyz[1] = omegay,  ioyz[2] = omegaz
//     ixyz[0] - I k  - moment tuhosti ve volnem krouceni

   eh=pow(g*ixyz[0]/e/ioyz[0],0.5); 
//   eh=1.; 
   eh1=eh*dl;
   s=sinh(eh1); c=cosh(eh1); 
   ck_matrix (ck,  s,c,eh,dl);

//    Integraly  cosh, sinh, x.cosh, x.sinh, cosh^2, cosh*sinh, sinh^2
   c1=s/eh; s1=(c-1.)/eh; xc = (dl*s-s1)/eh; xs = (dl*c-c1)/eh; c2=s*c/2/eh+dl/2; cs=s*s/2/eh;  s2=s*c/2/eh-dl/2;
   pom1=e*ioyz[0]*eh*eh*eh*eh;
   pom2=g*ixyz[0];
//   pom2=pom1/eh/eh;
   pom3=pom1/eh/eh;
   for  (i=0;i<4;i++){
       i1=4*i-i*(i+1)/2; //1=4*(i-1)-i*(i-1)/2;
       for  (j=i;j<4;j++){
         i2=j+i1;
         sm1[i2]=(ck[i][2]*ck[j][2]*c2+ck[i][3]*ck[j][3]*s2+
              (ck[i][3]*ck[j][2]+ck[i][2]*ck[j][3])*cs)*pom1;
         sm1[i2]=sm1[i2]+(ck[i][1]*ck[j][1]*dl+eh*(ck[i][1]*ck[j][2]+
               ck[i][2]*ck[j][1])*s1+eh*(ck[i][1]*ck[j][3]+ck[i][3]*ck[j][1])*c1)*pom2;
         sm1[i2]=sm1[i2]+((ck[i][2]*ck[j][3]+ck[i][3]*ck[j][2])*cs+
			   ck[i][2]*ck[j][2]*s2+ck[i][3]*ck[j][3]*c2)*eh*eh*pom2;
		}
   }

   sm[3][3]  = sm1[0] ;
   sm[3][6]  = sm1[1] ;
   sm[3][10] = sm1[2] ;
   sm[3][13] = sm1[3] ;
   sm[6][6]  = sm1[4] ;
   sm[6][10] = sm1[5] ;
   sm[6][13] = sm1[6] ;
   sm[10][10]= sm1[7] ;
   sm[10][13]= sm1[8] ;
   sm[13][13]= sm1[9] ;

//   gy=0.;    gz=0.;
//   druhe der w
//   n[0] = (6.-12.*s)/ll; n[1]= (-2.*(2.+gy)+6.*s)/dl; n[2]=-(6.-12.*s)/ll; n[3]= (-2.*(1.-gy)+6.*s)/dl;
   
   for  (i=0;i<4;i++){
        sm1[i]   = 6.*eh*eh/dl/dl*(-(ck[i][2]*c1+ck[i][3]*s1)+
                 2./dl*(ck[i][2]*xc+ck[i][3]*xs) );
        sm1[i+4] =-2.*eh*eh/dl*(-(2.+gy)*(ck[i][2]*c1+ck[i][3]*s1)+
                 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
        sm2[i+4] =-2.*eh*eh/dl*(-(2.+gz)*(ck[i][2]*c1+ck[i][3]*s1)+
                 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
        sm1[i+8] =-sm1[i];
        sm1[i+12]=-2.*eh*eh/dl*(-(1.-gy)*(ck[i][2]*c1+ck[i][3]*s1)+
                 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
        sm2[i+12]=-2.*eh*eh/dl*(-(1.-gz)*(ck[i][2]*c1+ck[i][3]*s1)+
                 3./dl*(ck[i][2]*xc+ck[i][3]*xs) );
   }

   sioyz[0]=ioyz[1]*e/(1.+2.*gy);
   sioyz[1]=ioyz[2]*e/(1.+2.*gz);
   for  (i=0;i<2;i++){
        i1=1-i;
        i2=i+1;
        sm[i2][3] = sm1[0]*sioyz[i1];
        sm[i2][6] = sm1[1]*sioyz[i1];
        sm[i2][10]= sm1[2]*sioyz[i1];
        sm[i2][13]= sm1[3]*sioyz[i1];
        sm[3][4]  = sm1[4]*sioyz[0];
        sm[3][5]  =-sm2[4]*sioyz[1];
        sm[3][i+8]= sm1[8]*sioyz[i1];
        sm[3][11] = sm1[12]*sioyz[0];
        sm[3][12] =-sm2[12]*sioyz[1];

        sm[4][6]  = sm1[5]*sioyz[0];
        sm[5][6]  =-sm2[5]*sioyz[1];
        sm[4][10] = sm1[6]*sioyz[0];
        sm[5][10] =-sm2[6]*sioyz[1];
        sm[4][13] = sm1[7]*sioyz[0];
        sm[5][13] =-sm2[7]*sioyz[1];
        sm[6][i+8]= sm1[9]*sioyz[i1];
        sm[6][11] = sm1[13]*sioyz[0];
        sm[6][12] =-sm2[13]*sioyz[1];
        i2=i+8;
        sm[i2][10]= sm1[10]*sioyz[i1];
        sm[i2][13]= sm1[11]*sioyz[i1];
        sm[10][11]= sm1[14]*sioyz[0];
        sm[10][12]=-sm2[14]*sioyz[1];

        sm[11][13]= sm1[15]*sioyz[0];
        sm[12][13]=-sm2[15]*sioyz[1];
   }

}

// *******************************************************
//     TORSION      1        cubic function

void beamgen3d::stiffness_matrixtor1 (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz)
{
  double aj1,ll,integn4;
  vector sioyz(2),integn1(4),integn2(3),integn3(2);

//     ioyz[0] - I omega - vys. moment setrvac,   ioyz[1] = omegay,  ioyz[2] = omegaz
//     ixyz[0] - I k  - moment tuhosti ve volnem krouceni

//  gy=0.; gz=0.;
  ll=dl*dl; aj1=e*ioyz[0];
  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.+0.)/dl*aj1;  integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-0.)/dl*aj1;
  integn3[0]=12./ll/dl*aj1;      integn3[1]=-6./ll*aj1;
  integn4=2.*(2.+0.)/dl*aj1;
//first deriv. +, +, -, +      +,- ,-,   +,-,  +, 
  sm[3][3]  = integn1[0]+g*1.2*ixyz[0]/dl;
  sm[3][6]  = integn1[1]+g*0.1*ixyz[0];
  sm[3][10] = integn1[2]-g*1.2*ixyz[0]/dl;
  sm[3][13] = integn1[3]+g*0.1*ixyz[0];
  sm[6][6]  = integn2[0]+g*4./30.*ixyz[0]*dl;
  sm[6][10] = integn2[1]-g*0.1*ixyz[0];
  sm[6][13] = integn2[2]-g*1./30.*ixyz[0]*dl;
  sm[10][10]= integn3[0]+g*1.2*ixyz[0]/dl;
  sm[10][13]= integn3[1]-g*0.1*ixyz[0];
  sm[13][13]= integn4   +g*4./30.*ixyz[0]*dl;

//  direct s Ioy - w (index 2,9), fy (index 4,11)
  aj1=ioyz[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][3] = integn1[0];
   sm[2][6] = integn1[1];
   sm[2][10]= integn1[2];
   sm[2][13]= integn1[3];
    sm[3][4] = -integn1[1];
   sm[3][8]=integn1[2];
    sm[3][11] = -integn1[3];
    sm[4][6]  = -integn2[0];
    sm[4][10] = -integn2[1];
    sm[4][13] = -integn2[2];
   sm[6][9]= integn2[1];
    sm[6][11] = -integn2[2];
   sm[9][10]= integn3[0];
   sm[9][13]= integn3[1];
    sm[10][11]= -integn3[1];
    sm[11][13]= -integn4;

//  direct s Ioz - v (index 1,8), fz (index 5,12)
   aj1=ioyz[2]*e/(1.+2.*gz);
   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.+gz)/dl*aj1;  integn2[1]=-6./ll*aj1; integn2[2]=2.*(1.-gz)/dl*aj1;
   integn3[0]=12./ll/dl*aj1;      integn3[1]=-6./ll*aj1;
   integn4=2.*(2.+gz)/dl*aj1;
   sm[1][3] = integn1[0];
   sm[1][6] = integn1[1];
   sm[1][10]= integn1[2];
   sm[1][13]= integn1[3];
    sm[3][5] = integn1[1];
   sm[3][8]=integn1[2];
    sm[3][12] = integn1[3];
    sm[5][6]  = integn2[0];
    sm[5][10] = integn2[1];
    sm[5][13] = integn2[2];
   sm[6][8]= integn2[1];
    sm[6][12] = integn2[2];
   sm[8][10]= integn3[0];
   sm[8][13]= integn3[1];
    sm[10][12]= integn2[1];
    sm[12][13]= integn4;

 /* 
  for (i=0;i<2;i++){
      i1=1-i;
      i2=i+1;
      sm[i2][3] = 12.*sioyz[i1]/dl;
      sm[i2][6] = 6.*sioyz[i1];
      sm[i2][10]= -sm[i2][3];
      sm[i2][13]=  sm[i2][6];
      sm[3][4]  =-6.*sioyz[0];
      sm[3][5]  = 6.*sioyz[1];
      sm[3][i+8]=-12.*sioyz[i1]/dl;
      sm[3][11] =-6.*sioyz[0];
      sm[3][12] = 6.*sioyz[1];
      sm[4][6]  =-4.*sioyz[0]*dl;
      sm[5][6]  = 4.*sioyz[1]*dl;
      sm[4][10] = 6.*sioyz[0];
      sm[5][10] =-6.*sioyz[1];
      sm[4][13] =-2.*sioyz[0]*dl;
      sm[5][13] = 2.*sioyz[1]*dl;
      sm[6][i+8]=-6.*sioyz[i1];
      sm[6][11] =-2.*sioyz[0]*dl;
      sm[6][12] = 2.*sioyz[1]*dl;
      i2=i+8;
      sm[i2][10]= 12.*sioyz[i1]/dl;
      sm[i2][13]= -6.*sioyz[i1];
      sm[10][11]= 6.*sioyz[0];
      sm[10][12]=-6.*sioyz[1];
      sm[11][13]=-4.*sioyz[0]*dl;
      sm[12][13]= 4.*sioyz[1]*dl;
  }
  */
}

// *******************************************************
//     TORSION      2        kubic function Mindlin element
//
void beamgen3d::stiffness_matrixtor2 (matrix &sm, double dl,double e,double g,double gy,double gz,double *ixyz,double *ioyz,double *iro)
{
  double  kakr,ajy,ajz,ajk;

//  Razeni v matici tuhosti

//  N1, Vy1, Vz1, Mx1, My1, Mz1, B1
//  N2, Vy2, Vz2, Mx2, My2, Mz2, B2

//  K R U T   -  Mx, Qz, My
//     ioyz[0] - I omega - vys. moment setrvac,   ioyz[1] = omegay,  ioyz[2] = omegaz
//     ixyz[0] - I k  - moment tuhosti ve volnem krouceni
//     iro[0] -  Iro   Irosin  Irocos
    kakr  =6.*ioyz[0]*e/iro[0]/g/dl/dl;
    
	ajy= (1.+2.*gy);
    ajk= (1.+2.*kakr);
    ajy=1.;
    sm[2][3] = (12.*ioyz[1]*e/dl/dl/dl+4.*kakr*gy*iro[1]*g/dl)/ajk/ajy;
    sm[3][3] = (12.*ioyz[0]*e/dl/dl/dl+(4.*kakr*kakr*iro[0] 
		      +ixyz[0]*(4.*kakr*(kakr+1.)+1.2) )*g/dl )/ajk/ajk;
    sm[3][4] = (-6.*ioyz[1]*e/dl/dl-2*kakr*gy*iro[1]*g)/ajk/ajy;
    sm[2][6] = sm[4][5];
    sm[3][6] = (-6.*ioyz[0]*e/dl/dl-(2.*kakr*kakr*iro[0]+iro[0]/10.)*g )/ajk/ajk;
    sm[4][6] = ((2.+kakr+gy+2.*kakr*gy)*2.*ioyz[1]*e/dl
		       +kakr*gy*iro[1]*g*dl)/ajk/ajy;
    sm[6][6] = ((1.+(1.+kakr)*kakr)*4.*ioyz[0]*e/dl+(kakr*kakr*iro[0]
		       +(kakr*(kakr+1)+0.4)/3.*ixyz[0])*g*dl)/ajk/ajk;

    sm[3][9] = -sm[2][3]; sm[6][9] = -sm[3][4];
    sm[2][10] = -sm[2][3]; sm[3][10] = -sm[3][3];
    sm[4][10] = -sm[3][4]; sm[6][10] = -sm[3][6];
    sm[9][10] = sm[2][3]; sm[10][10]=  sm[3][3];

    sm[3][11] =  sm[3][4];
    sm[6][11] = ((1.-kakr-gy-2.*kakr*gy)*2.*ioyz[1]*e/dl
                 -kakr*gy*iro[1]*g*dl)/ajk/ajy;
    sm[10][11]= -sm[3][4]; sm[2][13] =  sm[2][6];
    sm[3][13] =  sm[3][6];
    sm[4][13] = ((1.-kakr-gy-2.*kakr*gy)*2.*ioyz[1]*e/dl
                +kakr*gy*iro[1]*g*dl)/ajk/ajy;
    sm[6][13] = ((2.-(4.+4.*kakr)*kakr)*ioyz[0]*e/dl+(kakr*kakr*iro[0]
                +(kakr*kakr-0.1)/3.*ixyz[0] )*g*dl)/ajk/ajk;
    sm[9][13] = -sm[2][6];sm[10][13] = -sm[3][6];
    sm[11][13] =  sm[4][6];sm[13][13] =  sm[6][6];

//  K R U T   -  Mx, Qy, Mz
      ajz= (1.+2.*gz);
      ajz=1.;
      sm[1][3] = (12.*ioyz[2]*e/dl/dl/dl+4.*kakr*gz*iro[2]*g/dl)/ajk/ajz;
      sm[3][5] = (-6.*ioyz[2]*e/dl/dl-2.*kakr*gz*iro[2]*g)/ajk/ajz;
      sm[1][6] = sm[3][5];
      sm[5][6] = ((2.+kakr+gz+2.*kakr*gz)*2.*ioyz[2]*e/dl
                +kakr*gz*iro[2]*g*dl)/ajk/ajz;

      sm[3][8]  = -sm[1][3];sm[6][8]  = -sm[3][5];
      sm[1][10] = -sm[1][3];sm[5][10] = -sm[3][5];
      sm[8][10] =  sm[1][3];

      sm[3][12] =  sm[3][5];
      sm[6][12] = ((1.-kakr-gz-2.*kakr*gz)*2.*ioyz[2]*e/dl
                -kakr*gz*iro[2]*g*dl)/ajk/ajz;
      sm[10][12]= -sm[4][6];sm[2][14] =  sm[2][7];
      sm[5][13] = ((1.-kakr-gz-2.*kakr*gz)*2.*ioyz[2]*e/dl
                +kakr*gz*iro[2]*g*dl)/ajk/ajz;
      sm[8][13] = -sm[1][6];
      sm[12][13]=  sm[5][6];

}


/**
   function serves for computation of nodal forces and moments caused by body load
   
   @param eid - element id
   @param lm - load %matrix
   
   PF, 1.9.2006
*/
void beamgen3d::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(3,3);
  
  //  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
//   !!!zmenit   lgmatrixtransf3dblock (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
   
   PF, 20.9.2006
*/
void beamgen3d::nodal_displ (long eid,long lcid)
{
  long ii,transf;
  double l;
  ivector nodes(nne),cn(ndofe);