soilbeam.cpp

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

  //  transformation of stiffness matrix to the global system
  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 soil 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][i]);
//   }
//  }
//  fprintf (Out,"\n\n MT v MT Soilbeam cislo %ld \n",eid);
//  for (i=0;i<3;i++){fprintf (Out," %15.5e",c1[i]);}
//  for (i=0;i<3;i++){fprintf (Out," %15.5e",c2[i]);}
}



void soilbeam::strains (long lcid,long eid,long ri,long ci)
{
  long ipp,i,j,ij,transf;
  double dl;
  ivector cn(ndofe),nodes(nne);
  vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),eps(6);
  matrix tran(3,3);
  
// r  .... end forces last time
// m  .... dl(i)/dl(i-1)
  ipp=Mt->elements[eid].ipp[ri][ci];
  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);
  
  //  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);   

  ij=0;
  for (i=0;i<intordsm[0][0];i++){
    for (j=0;j<6;j++){  
      eps[j]=rl[ij];
	  ij++;
	}
//  fprintf (Out,"\n\n kontrola deformaci v soilbeam zapis,%ld \n",ipp);
//  for (long j=0;j<eps.n;j++){fprintf (Out," %15.5e",eps[j]);}
    Mm->storestrain(lcid,ipp,eps);
	ipp++;
  }
}
	  
void soilbeam::res_internal_forces (long lcid,long eid,vector &ifor)
{
   internal_forces (lcid, eid, 0, 0, ifor);
}

/**
   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 soilbeam::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long ipp,imat,i,j,k;
  double dl,ll,a2,a4,ixyz[3],beta[2],a,e,g,gy,gz,g2;
  ivector nodes(nne);
  vector x(nne),y(nne),z(nne),vec(3),rl(ndofe),rg(ndofe),f(ndofe),eps(6);
  matrix d(tncomp,tncomp),c(tncomp,tncomp),tran(3,3),r(2,ndofe);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
  
  ipp=Mt->elements[eid].ipp[ri][ci];
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  Mm->matstiff (c,ipp);
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  Mc->give_vectorlcs (eid, vec);
  beam_transf_matrix (tran,dl, vec,x,y,z,eid);
  e=Mm->eliso[imat].e;  g=e/(1.+2.*Mm->eliso[imat].nu);
//  Mm->elmatstiff (d,ipp);
//  e=d[0][0];  g=d[1][1]; 
  c1[0]=c[0][0];
  c1[1]=c[1][1];
  c1[2]=c[2][2];
  c2[0]=c[3][3];
  c2[1]=c[4][4];
  c2[2]=c[5][5];

  k=0;
  for (i=0;i<intordsm[0][0];i++){  
    Mm->computenlstresses (ipp);
    Mm->givestress (lcid,ipp,eps);
    for (j=0;j<6;j++){  rl[k]=eps[j]; k++;}
	ipp++;
  }
  
  ll=dl*dl;
  f[0]= c1[0]*dl/3*rl[0]+ c1[0]*dl/6*rl[6]; 
  f[6]= c1[0]*dl/6*rl[0]+ c1[0]*dl/3*rl[6]; 
  g2=(c2[1]+c2[2])/2.;
  a2=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt((c1[1]+c1[2])/2.*g2) );
  a4=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt(g2*g2*g2/((c1[1]+c1[2])/2.)) );
  f[3]=  (a2/3.*dl+a4/dl)*rl[3]+(a2/6.*dl-a4/dl)*rl(9);
  f[9]=  (a2/6.*dl-a4/dl)*rl[3]+(a2/3.*dl+a4/dl)*rl(9);
  
//  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;

//  direct s Iy
  g2=gy*gy;
  a2=bPod*(c1[2]+sqrt(c1[2]*c2[2])*2./bPod)*dl/(1.+2.*gy)/(1.+2.*gy);
  a4=bPod*(c2[2]+sqrt(c2[2]*c2[2]*c2[2]/c1[2])/bPod);
//  a4=0.0;
  f[2] =   (a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4*(g2+gy+0.3))*rl[2]
	      +(a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+0.3))*rl[8]
	      -(a2*(g2/6.+11.*gy/60.+11./210.)+a4/10.)*dl*rl[4]
	      +(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[10];

//  f[4] = -a2*( (g2/6.+11.*gy/60.+11./210.)*rl[2]+(g2+0.9*gy+13./70.)/6.*rl[8] )*dl
//	     -a4/10.*(rl[4]-rl[10])*dl;
  f[4] = -(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[2]
         +(a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll*rl[4]
         -(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[8]
         -(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll*rl[10];
  f[8] =   (a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+0.3))*rl[2]
	      +(a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4*(g2+gy+0.3))*rl[8]
	      -(a2*(g2+0.9*gy+13./70.)/6.-a4/10.)*dl*rl[4]
		  +(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[10];

//  f[10]=  a2*( (g2+0.9*gy+13./70.)/6.*rl[2]+(g2/6.+11.*gy/60.+11./210.)*rl[8] )*dl
//	     +a4/10.*(rl[4]-rl[10])*dl;
  f[10]= +(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl*rl[2]
         -(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll*rl[4]
		 +(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl*rl[8]       
         +(a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll*rl[10];

//  direct s Iz
  g2=gz*gz;
  a2=bPod*(c1[1]+sqrt(c1[1]*c2[1])*2./bPod)*dl/(1.+2.*gz)/(1.+2.*gz);
  a4=bPod*(c2[1]+sqrt(c2[1]*c2[1]*c2[1]/c1[1])/bPod);
//  a4=0.0;
//  f[1] =  a2*( (4.*g2/3.+1.4*gz+13./35.)*rl[1]+(2./3.*g2+0.6*gz+9./70.)*rl[7] )
//	     +a4*4.*(g2+gz+.3)*(rl[5]-rl[11]);
  f[1] =   (a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4*(g2+gz+0.3))*rl[1]
	      +(a2*(2./3.*g2+0.6*gz+9./70.)-a4*4*(g2+gz+0.3))*rl[7]
	      -(a2*(g2/6.+11.*gz/60.+11./210.)+a4/10.)*dl*rl[5]
	      +(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[11];
//  f[5] =  a2*( (g2/6.+11.*gz/60.+11./210.)*rl[1]+(g2+0.9*gz+13./70.)/6.*rl[7] )*dl
//	     +a4/10.*(rl[5]-rl[11])*dl;
  f[5] = -(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[1]
         +(a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll*rl[5]
         -(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[7]
         -(a2*(g2/30.+gz/30.+1./140.)+a4*(g2+gz+0.1)/3.)*ll*rl[11];
//  f[7] =  a2*( (2./3.*g2+0.6*gz+9./70.)*rl[1]+(4.*g2/3.+1.4*gz+13./35.)*rl[7] )
//	     -a4*4.*(g2+gz+.3)*(rl[5]-rl[11]);
  f[7] =   (a2*(2./3.*g2+0.6*gz+9./70.)-a4*4*(g2+gz+0.3))*rl[1]
	      +(a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4*(g2+gz+0.3))*rl[7]
	      -(a2*(g2+0.9*gz+13./70.)/6.-a4/10.)*dl*rl[5]
		  +(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[11];

//		  -(a2*(g2/6.+11.*gz/60.+11./210.)+a4/10.)*(rl[5]-rl[11]);
//  f[11]= -a2*( (g2+0.9*gz+13./70.)/6.*rl[1]+(g2/6.+11.*gz/60.+11./210.)*rl[7] )*dl
//	     -a4/10.*(rl[5]-rl[11])*dl;
  f[11]= +(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl*rl[1]
         -(a2*(g2/30.+gz/30.+1./140.)+a4*(g2+gz+0.1)/3.)*ll*rl[5]
		 +(a2*(g2/6.+11.*gz/60.+11./210.)+ a4/10. )*dl*rl[7]       
         +(a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll*rl[11];

// do GCS z LCS  glvectortransf3dblock (f, tran);
  for (k=0;k<ndofe;k++){
    ifor[k]+=f[k];
  }
//  fprintf (Out,"\n\n R prvek Soilbeam cislo %ld \n",eid);
//  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",rl[k]);}
//  fprintf (Out,"\n\n Fint prvek Soilbeam cislo %ld \n",eid);
//  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",f[k]);}
//  fprintf (Out,"\n\n MT v Soilbeam cislo %ld \n",eid);
//  for (k=0;k<3;k++){fprintf (Out," %15.5e",c1[k]);}
//  for (k=0;k<3;k++){fprintf (Out," %15.5e",c2[k]);}
}
void soilbeam::internal_forces1 (long lcid,long eid,long ri,long ci,vector &ifor)
{
  long i,j,k,integr=2,ipp;
  double dl,ll,a,ixyz[3],beta[2],g2,gy,gz,a2x,a2y,a2z,e,g,s;
  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),eps(6);
  matrix n(tncomp,ndofe),c(tncomp,tncomp),d(tncomp,tncomp),tran(3,3),r(tncomp,2);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
  
  ipp=Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (c,ipp);
  Mm->elmatstiff (d,ipp);
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
  Mc->give_vectorlcs (eid, vec);
  beam_transf_matrix (tran,dl, vec,x,y,z,eid);
  e=d[0][0];  g=d[1][1]; ll=dl*dl;
  c1[0]=c[0][0];
  c1[1]=c[1][1];
  c1[2]=c[2][2];
  c2[0]=c[3][3];
  c2[1]=c[4][4];
  c2[2]=c[5][5];

  k=0;
  for (i=0;i<intordsm[0][0];i++){  
    Mm->computenlstresses (ipp);
    Mm->givestress (lcid,ipp,eps);
    for (j=0;j<6;j++){      rl[k]=eps[j];	 k++; }
	ipp++;
  }
  
//  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=(c2[1]+c2[2])/2.;
  a2x=bPod*bPod/2*( (c1[1]+c1[2])/2.*bPod/6.+sqrt((c1[1]+c1[2])/2.*g2) );
  g2=gy*gy;
  a2y=bPod*(c1[2]+sqrt(c1[2]*c2[2])*2./bPod)*dl/(1.+2.*gy)/(1.+2.*gy);
  g2=gz*gz;
  a2z=bPod*(c1[1]+sqrt(c1[1]*c2[1])*2./bPod)*dl/(1.+2.*gz)/(1.+2.*gz);
//	 s=x/dl
     s=0.; 
  for (i=0;i<2;i++){
     geom_matrix (n,s, dl, gy, gz);
     fx[0]=( n[0][0]*rl[0]+n[0][6]*rl[6] )*c1[0];
     fx[3]=( n[3][3]*rl[3]+n[3][9]*rl[9] )*a2x;
     fx[1]=( n[1][1]*rl[1]+n[1][5]*rl[5]+n[1][7]*rl[7]+n[1][11]*rl[11] )*a2y;
     fx[5]=( n[5][1]*rl[1]+n[5][5]*rl[5]+n[5][7]*rl[7]+n[5][11]*rl[11] )*a2y;
     fx[2]=( n[2][2]*rl[2]+n[2][4]*rl[4]+n[2][8]*rl[8]+n[2][10]*rl[10] )*a2z;
     fx[4]=( n[4][2]*rl[2]+n[4][4]*rl[4]+n[4][8]*rl[8]+n[4][10]*rl[10] )*a2z;
	 s=dl;
  }

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