soilplatetr.cpp

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

#include "soilplatetr.h"
#include "dkt.h"
#include "global.h"
#include "globmat.h"
#include "genfile.h"
#include "node.h"
#include "element.h"
#include <stdlib.h>
#include <math.h>

soilplatetr::soilplatetr (void)
{
  long i,j;

  nne=3;  ndofe=9;  tncomp=5;  napfun=3; ned=3;  nned=2;  
  intordmm=3; intordb=2; ssst=plate;
  nb=1; 
  
  ncomp = new long [nb];
  ncomp[0]=3;
  
  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]=7;
  
  tnip=0;
  for (i=0;i<nb;i++){
    for (j=0;j<nb;j++){
      tnip+=nip[i][j];
    }
  }
  
  intordsm[0][0]=7;
}

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

  delete [] cncomp;
  delete [] ncomp;
}

void soilplatetr::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];
    }
  }
}

double soilplatetr::approx (double xi,double eta,vector &nodval)
{
  double f;
  vector bf(nne);
  
  bf_lin_4_2d (bf.a,xi,eta);
  
  scprd (bf,nodval,f);
  
  return f;
}
/**
   function extracts components of the soil
   part of stiffness matrix of the material to the matrix ds
   20.3.2002
*/
void soilplatetr::dbmat (matrix &db,double c1, double c2)
{

  db[0][0] = c1;   db[0][1] = 0.0;  db[0][2] = 0.0;
  db[1][0] = 0.0;  db[1][1] = c2;   db[1][2] = 0.0;
  db[2][0] = 0.0;  db[2][1] = 0.0;  db[2][2] = c2;
}


void soilplatetr::geom_matrix (matrix &gm,matrix &cn,vector &ax,vector &ay,vector &dl, vector &l)
{
  int i,i1,i2,i4;
  double lll;
//  vector ll(3);

//  ll[0]=l(0)*l(1)/2.;
//  ll[1]=l(1)*l(2)/2.;
//  ll[2]=l(2)*l(0)/2.;
  lll=l(0)*l(1)*l(2)/2.;
  for (i=0;i<nne;i++){
	 i1=i+1;
	 if (i1>=nne) i1=i1-nne;
     i2=i1+1;     
	 if (i2>=nne) i2=i2-nne;
     i4=i*3;     
     gm[0][i4]  =  l[i];
     gm[0][i4+1]=(l[i]*l[i]*l[i1]+lll)*dl[i]*cn[i][0]-(l[i]*l[i]*l[i2]+lll)*dl[i2]*cn[i2][0];
     gm[0][i4+2]=(l[i]*l[i]*l[i1]+lll)*dl[i]*cn[i][1]-(l[i]*l[i]*l[i2]+lll)*dl[i2]*cn[i2][1];
     gm[1][i4]  = ay[i];
//     gm[1][i4+1]=-ay[i] *( (2.*l[i]*l[i2]+ll[i1])*dl[i2]*cn[i2][1]+(2.*l[i]*l[i1]+ll[i1])*dl[i]*cn[i][1] )
//                 -ay[i1]*( (ll[i2])*dl[i2]*cn[i2][1]+(l[i]*l[1]+ll[i2])*dl[i]*cn[i][1] )
//                 -ay[i2]*( (l[i]*l[i]+ll[i])*dl[i2]*cn[i2][1]+(ll[i])*dl[i]*cn[i][1] )

//     gm[1][i4+2]=-ay[i] *( (2.*l[i]*l[i2]+ll[i1])*dl[i2]*cn[i2][0]+(2.*l[i]*l[i1]+ll[i1])*dl[i]*cn[i][0] )
//                 -ay[i1]*( (ll[i2])*dl[i2]*cn[i2][0]+(l[i]*l[1]+ll[i2])*dl[i]*cn[i][0] )
//                 -ay[i2]*( (l[i]*l[i]+ll[i])*dl[i2]*cn[i2][0]+(ll[i])*dl[i]*cn[i][0] )

     gm[2][i4]  = ax[i];
//     gm[2][i4+1]=-ax[i] *( (2.*l[i]*l[i2]+ll[i1])*dl[i2]*cn[i2][1]+(2.*l[i]*l[i1]+ll[i1])*dl[i]*cn[i][1] )
//                 -ax[i1]*( (ll[i2])*dl[i2]*cn[i2][1]+(l[i]*l[1]+ll[i2])*dl[i]*cn[i][1] )
//                 -ax[i2]*( (l[i]*l[i]+ll[i])*dl[i2]*cn[i2][1]+(ll[i])*dl[i]*cn[i][1] )

//     gm[2][i4+2]=-ax[i] *( (2.*l[i]*l[i2]+ll[i1])*dl[i2]*cn[i2][0]+(2.*l[i]*l[i1]+ll[i1])*dl[i]*cn[i][0] )
//                 -ax[i1]*( (ll[i2])*dl[i2]*cn[i2][0]+(l[i]*l[1]+ll[i2])*dl[i]*cn[i][0] )
//                 -ax[i2]*( (l[i]*l[i]+ll[i])*dl[i2]*cn[i2][0]+(ll[i])*dl[i]*cn[i][0] )

  }
}

void soilplatetr::transf_matrix (ivector &nodes,matrix &tmat)
{
  long i,n,m;
  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){
      tmat[i*3+1][i*3+1]=Mt->nodes[nodes[i]].e1[0];  tmat[i*3+1][i*3+2]=Mt->nodes[nodes[i]].e2[0];
      tmat[i*3+2][i*3+1]=Mt->nodes[nodes[i]].e1[1];  tmat[i*3+2][i*3+2]=Mt->nodes[nodes[i]].e2[1];
    }
  }
}
void soilplatetr::stiffness_matrix (long eid,long ri, long ci, matrix &sm,vector &x, vector &y)
{
  long i,ipp,ij;
  double det,jac;
  ivector nodes(nne);
  vector w,gp1,gp2,c1(nne),c2(nne),l(3),dl(3),ax(3),ay(3);
  matrix gm,cn(3,2),cc(3,3);
//     jakobian
  Mt->give_elemnodes (eid,nodes);
//  Mc->give_soil (eid,nodes,c1,c2);
//  c1[0]=1000; c1[1]=1000;c1[2]=1000;
//  c2[0]=1000;c2[1]=1000;c2[2]=1000;

// triangular  
  det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]));
  dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
  dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
  dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
  // matrix vector of normal of sides {cos,sin} 
  cn[0][0]= (y[1]-y[0])/dl[0];
  cn[0][1]=-(x[1]-x[0])/dl[0];
  cn[1][1]=-(x[2]-x[1])/dl[1];
  cn[1][0]= (y[2]-y[1])/dl[1];
  cn[2][1]=-(x[0]-x[2])/dl[2];
  cn[2][0]= (y[0]-y[2])/dl[2];
  //   der   dL(1)/dy
  ax[0]=(x[2]-x[1])/det;  ax[1]=(x[0]-x[2])/det;  ax[2]=(x[1]-x[0])/det;
  //   der   dL(1)/dx
  ay[0]=(y[1]-y[2])/det;  ay[1]=(y[2]-y[0])/det;  ay[2]=(y[0]-y[1])/det;

  fillm (0.0,sm);
  
//  ii=Mt->elements[eid].ipp[sip][0];
  allocv (intordsm[0][0],w);  allocv (intordsm[0][0],gp1);  allocv (intordsm[0][0],gp2);
  allocm (ncomp[0],ndofe,gm); allocm (ncomp[0],ncomp[0],cc);
  gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);
  ipp=Mt->elements[eid].ipp[ri][ci];
  for (i=0;i<intordsm[0][0];i++){
     l[0]=gp1[i];  l[1]=gp2[i];  l[2]=1.0-l[0]-l[1];
     jac=det*w[i];
     geom_matrix (gm,cn,ax,ay,dl,l);
     Mm->matstiff (cc,ipp);
//     c11 = c1[0]*l[0]+c1[1]*l[1]+c1[2]*l[2];
//     c22 = c2[0]*l[0]+c2[1]*l[1]+c2[2]*l[2];
//    dbmat (cc,c11,c22);
     bdbjac (sm,gm,cc,gm,jac);
     ipp++;  
  }
  		fprintf (Out,"\n");
  fprintf (Out,"\n\n prvek cislo %ld \n",eid);
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[0][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[1][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[2][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[3][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[4][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[5][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[6][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[7][ij]);}fprintf (Out,"\n");
		for (ij=0;ij<ndofe;ij++){fprintf (Out,"%12.3e",sm[8][ij]);}fprintf (Out,"\n");
destrm (cc);  destrm (gm);  destrv (gp1);destrv (gp2);  destrv (w);
}
	  
void soilplatetr::res_stiffness_matrix (long eid,matrix &sm)
{
  vector x(nne),y(nne);
  Mt->give_node_coord2d (x,y,eid);
  stiffness_matrix (eid,0,0,sm,x,y);
}

void soilplatetr::res_mainip_strains (long lcid,long eid)
{
  vector aux,x(nne),y(nne),r(ndofe);
  ivector cn(ndofe),nodes(nne);
  matrix tmat;

  Mt->give_node_coord2d (x,y,eid);
  Mt->give_elemnodes (eid,nodes);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,r.a,cn.a,ndofe);

  //  transformation of displacement vector
  long transf = Mt->locsystems (nodes);
  if (transf>0){
    allocv (ndofe,aux);
    allocm (ndofe,ndofe,tmat);
    transf_matrix (nodes,tmat);
    //locglobtransf (aux,r,tmat);
    lgvectortransf (aux,r,tmat);
    copyv (aux,r);
    destrv (aux);
    destrm (tmat);
  }

  mainip_strains (lcid,eid,0,0,x,y,r);
}


void soilplatetr::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
{
  long i,ipp,ij;
  double det;
  vector ax(3),ay(3),dl(3),l(3),gp1,gp2,w,eps;
  matrix cn(3,2),gm;

  // translation to center it is necesary for w which in local coord. and no in relation
// triangular  
  det = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
  dl[0] = sqrt((x[1]-x[0])*(x[1]-x[0])+(y[1]-y[0])*(y[1]-y[0]));
  dl[1] = sqrt((x[2]-x[1])*(x[2]-x[1])+(y[2]-y[1])*(y[2]-y[1]));
  dl[2] = sqrt((x[0]-x[2])*(x[0]-x[2])+(y[0]-y[2])*(y[0]-y[2]));
  // matrix vector of normal of sides {cos,sin} 
  cn[0][0]= (y[1]-y[0])/dl[0];
  cn[0][1]=-(x[1]-x[0])/dl[0];
  cn[1][1]=-(x[2]-x[1])/dl[1];
  cn[1][0]= (y[2]-y[1])/dl[1];
  cn[2][1]=-(x[0]-x[2])/dl[2];
  cn[2][0]= (y[0]-y[2])/dl[2];
  //   der   dL(1)/dy
  ax[0]=(x[2]-x[1])/2./det;  ax[1]=(x[0]-x[2])/2./det;  ax[2]=(x[1]-x[0])/2./det;
  //   der   dL(1)/dx
  ay[0]=(y[1]-y[2])/2./det;  ay[1]=(y[2]-y[0])/2./det;  ay[2]=(y[0]-y[1])/2./det;
  
  fprintf (Out,"\n\n deformace");
  fprintf (Out,"\n\n prvek cislo %ld",eid);
  ipp=Mt->elements[eid].ipp[ri][ci];
  allocv (intordsm[0][0],w);  allocv (intordsm[0][0],gp1);  allocv (intordsm[0][0],gp2);
  allocm (ncomp[0],ndofe,gm);
  allocv (ncomp[0],eps);
  gauss_points_tr (gp1.a,gp2.a,w.a,intordsm[0][0]);

    for (i=0;i<intordsm[0][0];i++){
     l[0]=gp1[i];  l[1]=gp2[i];  l[2]=1.0-l[0]-l[1];
     geom_matrix (gm,cn,ax,ay,dl,l);
     mxv (gm,r,eps);
     Mm->storestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
     ipp++;
		fprintf (Out,"\n");
		for (ij=0;ij<ncomp[0];ij++){fprintf (Out,"%20.10e",eps[ij]);}
	}
    destrm (gm);  destrv (eps);  destrv (gp1);  destrv (gp2);  destrv (w);
  
}

void soilplatetr::res_allip_stresses (long lcid,long eid)
{
   allip_stresses (lcid, eid, 0, 0);
}

/**
   function computes stresess
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   
   20.9.2006
*/
void soilplatetr::allip_stresses (long lcid,long eid,long ri,long ci)
{
  long ipp,i;
  ivector nodes(nne);
  vector eps, sig;
  matrix cc;
  
  Mt->give_elemnodes (eid,nodes);
    
  // translation to center it is necesary for w which in local coord. and no in relation
//  elem_strains (stra,lcid,eid,ri,ci); 
  
  allocv (ncomp[0],eps);  allocv (ncomp[0],sig);
  allocm (ncomp[0],ncomp[0],cc);
        
  ipp=Mt->elements[eid].ipp[ri][ci];
 
//  fprintf (Out,"\n\n Eps,sig soilprvek cislo %ld\n",eid);
  for (i=0;i<intordsm[0][0];i++){
      Mm->matstiff (cc,ipp);
	  Mm->givestrain (lcid,ipp,cncomp[0],ncomp[0],eps);
	  mxv (cc,eps,sig);	
	  Mm->storestress (lcid,ipp,cncomp[0],sig);      
//  for (k=0;k<3;k++){
//    printf (Out," %15.5e",eps[k]); fprintf (Out," %15.5e",sig[k]);  }
//    fprintf (Out,"\n");
	  ipp++;  
  }
  destrv (eps);  destrv (sig); destrm (cc); 

//  fprintf (Out,"\n\n Fint soil prvek cislo %ld\n",eid);
//  for (k=0;k<ndofe;k++){fprintf (Out," %15.5e",ifor[k]);}
  
}