shelltr.cpp

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

#include "shelltr.h"
#include "global.h"
#include "globmat.h"
#include "genfile.h"
#include "node.h"
#include "element.h"
#include "plelemrotlt.h"
#include "dkt.h"
#include "intpoints.h"
#include <math.h>

shelltr::shelltr (void)
{
  long i,j;
  
  if (Perlt==NULL)   Perlt = new planeelemrotlt;
  if (Dkt==NULL)     Dkt = new dktelem;

  nne=3;  ndofe=18;  tncomp=6;  napfun=6;
  ned=3;  nned=2;    ssst=shell;
  
//  nb = Perlt->nb + Cct->nb;
  nb = Perlt->nb + Dkt->nb;

  ncomps = new long [Perlt->nb];
  ncomps[0]=2;
  ncomps[1]=1;
 
  ncompd = new long [Dkt->nb];
  ncompd[0]=3;

  nip = new long* [nb];
  intordsm = new long* [nb];
  for (i=0;i<nb;i++){
    nip[i]=new long [nb];
    intordsm[i] = new long [nb];
  }
  tncomps=Perlt->tncomp;
  tncompd=Dkt->tncomp;
  
  long a00,a01,a10,a11,b00;
  //long b01,b10,b11;

  a00=Perlt->nip[0][0];  a01=Perlt->nip[0][1];
  a10=Perlt->nip[1][0];  a11=Perlt->nip[1][1];
  b00=Dkt->nip[0][0];

  nip[0][0]=a00;  nip[0][1]=a01;  nip[0][2]=0;    
  nip[1][0]=a10;  nip[1][1]=a11;  nip[1][2]=0; 
  nip[2][0]=0;    nip[2][1]=0;    nip[2][2]=b00;
  
  a00=Perlt->intordsm[0][0];  a01=Perlt->intordsm[0][1];
  a10=Perlt->intordsm[1][0];  a11=Perlt->intordsm[1][1];
  b00=Dkt->intordsm[0][0];
  
  intordsm[0][0]=a00;  intordsm[0][1]=a01; intordsm[0][2]=0;
  intordsm[1][0]=a10;  intordsm[1][1]=a11; intordsm[1][2]=0;
  intordsm[2][0]=0;    intordsm[2][1]=0;   intordsm[2][2]=b00;

  dofe = new long [2];
  dofe[0]=Perlt->ndofe;
  dofe[1]=Dkt->ndofe;
  
  tnip=0;
  for (i=0;i<nb;i++){
    for (j=0;j<nb;j++){
      tnip+=nip[i][j];
    }
  }
  
}

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

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

}

/**
   function assembles transformation matrix between coordinate system
   defined on element and global coordinate system
   transformation deals with coordinates
   
   transformation matrix from local to global systems x,y,z
   In vectors   s  are global coordinates of Nodes 1 , 2 , 3
   Vector lokal s. s.( edge 1-2 = x' )  in glob. s. s.
   
   v(g)  = tran * v(l)
               t                t
   v(l)  = tran  * v(g)  =  v(g)  * tran
   
                               t
   sm(g) = tran  * sm(l) * tran
               t
   sm(l) = tran  * sm(g) * tran
   
*/
void shelltr::tran_mat(vector &x, vector &y, matrix &tran, vector &gx, vector &gy, vector &gz)
{
  long i,i1;
  double dl;
  matrix a(3,3);
 
  for (i=0; i<3; i++) {
    i1=i+1; if(i1>2)i1=i1-3;
//    a[i][0]=s2[i]-s1[i];
    a[0][i]=gx[i1]-gx[i];
//    a[i][1]=s3[i]-s2[i];
    a[1][i]=gy[i1]-gy[i];
//    a[i][2]=s1[i]-s3[i];
    a[2][i]=gz[i1]-gz[i];
  }
  
  dl=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);

  tran[0][0]=a[0][0]/dl;
  tran[1][0]=a[1][0]/dl;
  tran[2][0]=a[2][0]/dl;
  
  tran[0][2]=a[1][0]*a[2][1]-a[2][0]*a[1][1];
  tran[1][2]=a[2][0]*a[0][1]-a[0][0]*a[2][1];
  tran[2][2]=a[0][0]*a[1][1]-a[1][0]*a[0][1];

  dl=sqrt(tran[0][2]*tran[0][2]+tran[1][2]*tran[1][2]+tran[2][2]*tran[2][2]);

  tran[0][2]=tran[0][2]/dl;
  tran[1][2]=tran[1][2]/dl;
  tran[2][2]=tran[2][2]/dl;
  
  tran[0][1]=tran[1][2]*tran[2][0]-tran[2][2]*tran[1][0];
  tran[1][1]=tran[2][2]*tran[0][0]-tran[0][2]*tran[2][0];
  tran[2][1]=tran[0][2]*tran[1][0]-tran[1][2]*tran[0][0];

  dl=sqrt(tran[0][1]*tran[0][1]+tran[1][1]*tran[1][1]+tran[2][1]*tran[2][1]);

  tran[0][1]=tran[0][1]/dl;
  tran[1][1]=tran[1][1]/dl;
  tran[2][1]=tran[2][1]/dl;
/*  
// local coordinate x=s12
    sl1[0]=0.0;    sl1[1]=0.0;    sl1[2]=0.0;
    sl2[0]=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
    sl2[1]=0.0;    sl2[2]=0.0;
    sl3[0]=(s3(0)-s1(0))*tran(0,0)+(s3(1)-s1(1))*tran(1,0)+(s3(2)-s1(2))*tran(2,0);
    sl3[1]=(s3(0)-s1(0))*tran(0,1)+(s3(1)-s1(1))*tran(1,1)+(s3(2)-s1(2))*tran(2,1);
    sl3[2]=0.0;
*/
// local coordinate x=s12
    x[0]=0.0;
    y[0]=0.0;
    x[1]=sqrt(a[0][0]*a[0][0]+a[1][0]*a[1][0]+a[2][0]*a[2][0]);
    y[1]=0.0;
    x[2]=(gx(2)-gx(0))*tran(0,0)+(gy(2)-gy(0))*tran(1,0)+(gz(2)-gz(0))*tran(2,0);
    y[2]=(gx(2)-gx(0))*tran(0,1)+(gy(2)-gy(0))*tran(1,1)+(gz(2)-gz(0))*tran(2,1);
}



/**
   function assembles transformation matrix between global coordinate system
   and local nodal coordinate ssytem
   transformation deals with nodal forces

   transformation matrix x_g = T x_l

   9.5.2002
*/
void shelltr::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];
    }
  }
}


/**
   function assembles stiffness matrix of triangular shell element
   
   @param eid - element id
   @param sm - stiffness matrix
   
   4.5.2002
*/
void shelltr::res_stiffness_matrix (long eid,matrix &sm)
{
  long *rcn,*ccn,transf;
  vector gx(nne),gy(nne),gz(nne),x(nne),y(nne);
  ivector nodes(nne);
  matrix lsm,tran(3,3);
  
  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (gx,gy,gz,eid);
  //  funkce transformujici x,y,z do roviny xy 
  tran_mat(x,y, tran, gx,  gy,  gz);
  
  fillm (0.0,sm);
  
  rcn = new long [dofe[0]];
  ccn = new long [dofe[0]];
  allocm (dofe[0],dofe[0],lsm);
  
  //  contribution from triangular plane element with rotational degrees of freedom
  Perlt->stiffness_matrix (eid,0,0,lsm,x,y);
  
  rcn[0]=1;  rcn[1]=2;  rcn[2]=6;  rcn[3]=7;  rcn[4]=8;  rcn[5]=12;  rcn[6]=13;  rcn[7]=14;  rcn[8]=18;
  ccn[0]=1;  ccn[1]=2;  ccn[2]=6;  ccn[3]=7;  ccn[4]=8;  ccn[5]=12;  ccn[6]=13;  ccn[7]=14;  ccn[8]=18;
  
  mat_localize (sm,lsm,rcn,ccn);

  delete [] ccn;  delete [] rcn;  destrm (lsm);
  
  rcn = new long [dofe[1]];
  ccn = new long [dofe[1]];
  allocm (dofe[1],dofe[1],lsm);
  
  //  contribution from constant curve triangular plate element
  //Cct->stiffness_matrix (eid,2,2,lsm,x,y);
  Dkt->stiffness_matrix (eid,2,2,lsm,x,y);
  
  rcn[0]=3;  rcn[1]=4;  rcn[2]=5;  rcn[3]=9;  rcn[4]=10;  rcn[5]=11;  rcn[6]=15;  rcn[7]=16;  rcn[8]=17;
  ccn[0]=3;  ccn[1]=4;  ccn[2]=5;  ccn[3]=9;  ccn[4]=10;  ccn[5]=11;  ccn[6]=15;  ccn[7]=16;  ccn[8]=17;

  mat_localize (sm,lsm,rcn,ccn);
  //  transformation of stiffness matrix to GSS
  lgmatrixtransf3dblock(sm, tran);
  //  transformation of stiffness matrix to nodesystem
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
  /*
  fprintf (Out,"\n\n prvek cislo %ld",eid);
  for (long i=0;i<ndofe;i++){
    fprintf (Out,"\n %10.7e %10.7e %10.7e %10.7e %10.7e %10.7e",sm[i][0],sm[i][1],sm[i][2],sm[i][3],sm[i][4],sm[i][5]);
    fprintf (Out,"\n %10.7e %10.7e %10.7e %10.7e %10.7e %10.7e",sm[i][6],sm[i][7],sm[i][8],sm[i][9],sm[i][10],sm[i][11]);
    fprintf (Out,"\n %10.7e %10.7e %10.7e %10.7e %10.7e %10.7e",sm[i][12],sm[i][13],sm[i][14],sm[i][15],sm[i][16],sm[i][17]);
  }
 */ 
  delete [] ccn;  delete [] rcn;  destrm (lsm);
  
}

void shelltr::appval (vector &l,vector &eps,double **val)
{
  long j;
  vector epsp(tncompd);
  Perlt->appval (l[0],l[1],0,tncomps,eps,val);
  Dkt->appval (l,2,tncompd,epsp,val);
  for (j=0;j<tncompd;j++){
     eps[j+tncomps] = epsp[j];
  }
}
  
void shelltr::res_mainip_strains (long lcid,long eid)
{
  long *rcn,i;
  vector aux;
  vector gx(nne),gy(nne),gz(nne),x(nne),y(nne),r(ndofe),rs(dofe[0]),rd(dofe[1]);
  ivector cn(ndofe),nodes(nne);
  matrix lsm,tmat,tran(3,3);
  
  Mt->give_node_coord3d (gx,gy,gz,eid);
  Mt->give_elemnodes (eid,nodes);
  Mt->give_code_numbers (eid,cn.a);
    
  rcn = new long [dofe[0]];
  allocm (dofe[0],dofe[0],lsm);
  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);
  }
 //  funkce transformujici x,y,z do roviny xy 
  tran_mat(x,y, tran, gx,  gy,  gz);
  glvectortransf3dblock (r, tran);
    
  rcn[0]=1;  rcn[1]=2;  rcn[2]=6;  rcn[3]=7;  rcn[4]=8;  rcn[5]=12;  rcn[6]=13;  rcn[7]=14;  rcn[8]=18;
  for (i=0;i<Perlt->ndofe;i++){
    if (rcn[i]==0)  rs[i]=0.0;
    else rs[i]=r[rcn[i]-1];
  }
  //  contribution from triangular plane element with rotational degrees of freedom
  Perlt->mainip_strains (lcid,eid,0,0,x,y,rs);
    
  rcn[0]=3;  rcn[1]=4;  rcn[2]=5;  rcn[3]=9;  rcn[4]=10;  rcn[5]=11;  rcn[6]=15;  rcn[7]=16;  rcn[8]=17;
  for (i=0;i<Dkt->ndofe;i++){
    if (rcn[i]==0)  rd[i]=0.0;
    else rd[i]=r[rcn[i]-1];
  }
  //  contribution from triangular plate element
  Dkt->mainip_strains (lcid,eid,2,2,x,y,rd);
  
}
void shelltr::nod_strains (long lcid,long eid)
{
  
  //  contribution from triangular plane element with rotational degrees of freedom
  Perlt->nod_strains (lcid,eid,0,0);
    
  //  contribution from triangular plate element
  Dkt->nod_strains (lcid,eid,2,2);
}
  
void shelltr::elem_strains (double **stra,long lcid,long eid)
{
  long i,j;
  double **strap;
  //  contribution from triangular plane element with rotational degrees of freedom
  strap = new double* [nne];
  for (i=0;i<nne;i++){
    strap[i] = new double [tncomps];
  }
  
  //  contribution from triangular plane element with rotational degrees of freedom
  Perlt->elem_strains (strap,lcid,eid,0,0);
/*  
  fprintf (Out,"\n\n strains na stene (prvek %ld)",eid);
  for (i=0;i<nne;i++){
    fprintf (Out,"\n uzel %ld",i);
    for (j=0;j<tncomps;j++){
      fprintf (Out,"   %f",strap[i][j]);
    }
  }
 */ 

  for (i=0;i<nne;i++){
    for (j=0;j<tncomps;j++){
     stra[i][j] = strap[i][j];
    }
  }
  for (i=0;i<nne;i++){
     delete [] strap[i];