soilbeam.cpp

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

#include "soilbeam.h"
#include "global.h"
#include "globmat.h"
#include "genfile.h"
#include "node.h"
#include "element.h"
#include "beamel3d.h"
#include "intpoints.h"
#include "elastisomat.h"
#include <stdlib.h>
#include <math.h>

soilbeam::soilbeam (void)
{
  long i;
  
  nne=2;  ndofe=12;  tncomp=6;  napfun=3;
  nb=1;  intordmm=4;  intordism=2;  ssst=spacebeam;

  ncomp = new long [nb];
  ncomp[0]=1;
  
  cncomp = new long [nb];
  cncomp[0]=0;

  nip = new long* [nb];
  intordsm = new long* [nb];
  for (i=0;i<nb;i++){
    nip[i] = new long [nb];
    intordsm[i] = new long [nb];
  }
  nip[0][0]=2;
  intordsm[0][0]=2;

  c1= new double [3];
  c2= new double [3];

//  c1[0]=1000.;c1[1]=1000.;c1[2]=1000.;
//  c2[0]=0.001;c2[1]=0.001;c2[2]=0.001;
  
  bPod=1.;
  tnip=2;

}

soilbeam::~soilbeam (void)
{
  long i;
  
  for (i=0;i<nb;i++){
    delete [] nip[i];
  }
  delete [] nip;
}

void soilbeam::eleminit (long eid)
{
  long ii,jj;
  Mt->elements[eid].nb=nb;
  Mt->elements[eid].intordsm = new long* [nb];
  Mt->elements[eid].nip = new long* [nb];

  for (ii=0;ii<nb;ii++){
    Mt->elements[eid].intordsm[ii] = new long [nb];
    Mt->elements[eid].nip[ii] = new long [nb];
    for (jj=0;jj<nb;jj++){
      Mt->elements[eid].intordsm[ii][jj]=intordsm[ii][jj];
      Mt->elements[eid].nip[ii][jj]=nip[ii][jj];
    }
  }
}
void soilbeam::transf_matrix (ivector &nodes,matrix &tmat)
{
  long i,i6,n,m;

  fillm (0.0,tmat);

  n=nodes.n;
  m=tmat.m;
  for (i=0;i<m;i++){
    tmat[i][i]=1.0;
  }
  
  for (i=0;i<n;i++){
    if (Mt->nodes[nodes[i]].transf>0){
	  i6=i*6;
      tmat[i6][i6]   = Mt->nodes[nodes[i]].e1[0];  tmat[i6][i6+1]   = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2]     = Mt->nodes[nodes[i]].e3[0];
      tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1];  tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2]   = Mt->nodes[nodes[i]].e3[1];
      tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2];  tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2]   = Mt->nodes[nodes[i]].e3[2];
	  i6=i*6+3;
      tmat[i6][i6]   = Mt->nodes[nodes[i]].e1[0];  tmat[i6][i6+1]   = Mt->nodes[nodes[i]].e2[0]; tmat[i6][i6+2]     = Mt->nodes[nodes[i]].e3[0];
      tmat[i6+1][i6] = Mt->nodes[nodes[i]].e1[1];  tmat[i6+1][i6+1] = Mt->nodes[nodes[i]].e2[1]; tmat[i6+1][i6+2]   = Mt->nodes[nodes[i]].e3[1];
      tmat[i6+2][i6] = Mt->nodes[nodes[i]].e1[2];  tmat[i6+2][i6+1] = Mt->nodes[nodes[i]].e2[2]; tmat[i6+2][i6+2]   = Mt->nodes[nodes[i]].e3[2];
      
    }
  }
}
void soilbeam::beam_transf_matrix (matrix &tmat,double &dl,vector &vec,vector &x,vector &y,vector &z,long eid)
{
  double c;
  fillm (0.0,tmat);
  
// tmat[0][0] = c;    tmat[0][1] = -1.0*s;    tmat[0][2] = 0.0;
// tmat[1][0] = s;    tmat[1][1] = c;         tmat[1][2] = 0.0;
// tmat[2][0] = 0.0;  tmat[2][1] = 0.0;       tmat[2][2] = 1.0;

// tmat[3][3] = c;    tmat[3][4] = -1.0*s;    tmat[3][5] = 0.0;
// tmat[4][3] = s;    tmat[4][4] = c;         tmat[4][5] = 0.0;
// tmat[5][3] = 0.0;  tmat[5][4] = 0.0;       tmat[5][5] = 1.0;
   tmat[0][0]=x[1]-x[0];
   tmat[1][0]=y[1]-y[0];
   tmat[2][0]=z[1]-z[0];
   dl=sqrt((tmat[0][0]*tmat[0][0]+tmat[1][0]*tmat[1][0]+tmat[2][0]*tmat[2][0]));
   if (dl<Mp->zero){
    fprintf (stderr,"\n\n zero length of the %ld beamel3d element",eid);
    fprintf (stderr,"\n in function beamel3d::beam_transf_matrix.\n");
   }
   tmat[0][0]=tmat[0][0]/dl;
   tmat[1][0]=tmat[1][0]/dl;
   tmat[2][0]=tmat[2][0]/dl;

   tmat[0][1] =vec[1]*tmat[2][0]-vec[2]*tmat[1][0];
   tmat[1][1] =vec[2]*tmat[0][0]-vec[0]*tmat[2][0];
   tmat[2][1] =vec[0]*tmat[1][0]-vec[1]*tmat[0][0];
   c=sqrt((tmat[0][1]*tmat[0][1]+tmat[1][1]*tmat[1][1]+tmat[2][1]*tmat[2][1]));
   if (c<Mp->zero){
    fprintf (stderr,"\n\n zero vec of the %ld beamel3d element",eid);
    fprintf (stderr,"\n in function beamel3d::beam_transf_matrix.\n");
   }
   tmat[0][1]=tmat[0][1]/c;
   tmat[1][1]=tmat[1][1]/c;
   tmat[2][1]=tmat[2][1]/c;

   tmat[0][2]=tmat[1][0]*tmat[2][1]-tmat[2][0]*tmat[1][1];
   tmat[1][2]=tmat[2][0]*tmat[0][1]-tmat[0][0]*tmat[2][1];
   tmat[2][2]=tmat[0][0]*tmat[1][1]-tmat[1][0]*tmat[0][1];
   c=sqrt((tmat[0][2]*tmat[0][2]+tmat[1][2]*tmat[1][2]+tmat[2][2]*tmat[2][2]));
   if (c<Mp->zero){
    fprintf (stderr,"\n\n zero vec of the %ld beamel3d element",eid);
    fprintf (stderr,"\n in function beamel2d::beam_transf_matrix.\n");
   }
   tmat[0][2]=tmat[0][2]/c;
   tmat[1][2]=tmat[1][2]/c;
   tmat[2][2]=tmat[2][2]/c;
}
void soilbeam::geom_matrix (matrix &n,double s,double dl,double gy,double gz)
{
  double ll,aj1;
  ll=dl*dl;
//  u
        n[0][0] = 1.-s;
        n[0][6] = s;
//  w
        aj1=1./(1.+2.*gy);
        n[2][2] =aj1*(1.+2.*gy-2.*gy*s-3.*s*s+2*s*s*s);
        n[2][4] =aj1*(-(1.+gy)*s+(2.+gy)*s*s-s*s*s)*dl;
        n[2][8] =aj1*(2.*gy*s+3.*s*s-2.*s*s*s);
        n[2][10]=aj1*(gy*s+(1.-gy)*s*s-s*s*s)*dl;
// fy
        n[4][2] = aj1*(6.*s-6.*s*s)/dl;
        n[4][4] = aj1*(1.+2.*gy-2.*(2.+gy)*s+3.*s*s);
        n[4][8] =-n[4][2];
        n[4][10]= aj1*(-2.*(1.-gy)*s+3.*s*s);

//  v
        aj1=1./(1.+2.*gz);
        n[1][1] = aj1*(1.+2.*gz-2.*gz*s-3.*s*s+2.*s*s*s);
        n[1][5] =-aj1*(-(1.+gz)*s+(2.+gz)*s*s-s*s*s)*dl;
        n[1][7] = aj1*(2.*gz*s+3.*s*s-2.*s*s*s);
        n[1][11]=-aj1*(gz*s+(1.-gz)*s*s-s*s*s)*dl;
// fz
        n[5][1] =-aj1*(6.*s-6.*s*s)/dl;
        n[5][5] = aj1*(1.+2.*gz-2.*(2.+gz)*s+3.*s*s);
        n[5][7] =-n[5][1];
        n[5][11]= aj1*(-2.*(1.-gz)*s+3.*s*s);
// fx
        n[3][3] = 1.-s;
        n[3][9] = s;
  
}


void soilbeam::res_stiffness_matrix (long eid,matrix &sm)
{
  stiffness_matrix (eid,0,0,sm);
}

void soilbeam::stiffness_matrix (long eid,long ri,long ci,matrix &sm)
{
  long i,ipp,i1,transf, imat;
  double a2,a4,e,g,a,ixyz[3],beta[2],dl,ll,gy,gz,g2;
  ivector nodes(nne);
  vector vec(3),x(nne),y(nne),z(nne);
  matrix c(tncomp,tncomp),tran(3,3);
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (x,y,z,eid);
  Mc->give_vectorlcs (eid, vec);
  beam_transf_matrix (tran,dl, vec,x,y,z,eid);
  
  fillm (0.0,sm);
  ipp=Mt->elements[eid].ipp[ri][ci];
  Mm->matstiff (c,ipp);
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  Mc->give_area (eid,a);
  Mc->give_mominer (eid,ixyz);
  Mc->give_shearcoeff (eid,beta);
//  Mc->give_widht (eid,&bPod);
  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];
  ll=dl*dl;
//  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;

// axial forces and torsion
  sm[0][0]= c1[0]*dl/3;
  sm[0][6]= c1[0]*dl/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.)) );
  sm[3][3]=  a2/3.*dl+a4/dl;
  sm[3][9]=  a2/6.*dl-a4/dl;

//  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;
  sm[2][2] =  a2*(4.*g2/3.+1.4*gy+13./35.)+a4*4.*(g2+gy+.3);
  sm[2][4] =-(a2*(g2/6.+11.*gy/60.+11./210.)+ a4/10. )*dl;
  sm[2][8] =  a2*(2./3.*g2+0.6*gy+9./70.)-a4*4*(g2+gy+.3);
  sm[2][10]= (a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl;
  sm[4][4] = (a2*(g2/30.+gy/30.+1./105.)+a4*(g2+gy+0.4)/3.)*ll;
  sm[4][8] =-(a2*(g2+0.9*gy+13./70.)/6.-a4/10. )*dl;
  sm[4][10]=-(a2*(g2/30.+gy/30.+1./140.)+a4*(g2+gy+0.1)/3.)*ll;
//  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;
  sm[1][1] =  a2*(4.*g2/3.+1.4*gz+13./35.)+a4*4.*(g2+gz+.3);
  sm[1][5] = (a2*(g2/6.+11.*gz/60.+11./210.)+a4/10. )*dl;
  sm[1][7] =  a2*(2./3.*g2+0.6*gz+9./70.)-a4*4.*(g2+gz+.3);
  sm[1][11]=-(a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl;
  sm[5][5] = (a2*(g2/30.+gz/30.+1./105.)+a4*(g2+gz+0.4)/3.)*ll;
  sm[5][11]=-(a2*(g2/30.+gz/30+1./140.)+a4*(g2+gz+0.1)/3.)*ll;
  sm[5][7] = (a2*(g2+0.9*gz+13./70.)/6.-a4/10. )*dl;

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

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