linhex.cpp

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

    mxv (am,av,v);
    addv (v,nf,nf);
  }

  //  surface number 5

  if (is[4]>0 ){
    zeta=1.0;
    for (i=0;i<intordb;i++){
      xi=gp[i];
      for (j=0;j<intordb;j++){
	eta=gp[j];
	
	jac2d_3d (jac,x,y,z,xi,eta,4);
	bf_matrix (n,xi,eta,zeta);
	jac = jac*w[i]*w[j];
	nnj (am.a,n.a,jac,n.m,n.n);
      }
    }
    
    if (is[4]==1){
      av[0] = nv[48+0*3+0];
      av[1] = nv[48+0*3+1];
      av[2] = nv[48+0*3+2];
      
      av[3] = nv[48+1*3+0];
      av[4] = nv[48+1*3+1];
      av[5] = nv[48+1*3+2];
      
      av[6] = nv[48+2*3+0];
      av[7] = nv[48+2*3+1];
      av[8] = nv[48+2*3+2];

      av[9] = nv[48+3*3+0];
      av[10] = nv[48+3*3+1];
      av[11] = nv[48+3*3+2];
    }


    if (is[4]==2){
      
      av[0] = nv[48+0*3+0];
      av[1] = nv[48+0*3+1];
      av[2] = nv[48+0*3+2];
      
      av[3] = nv[48+1*3+0];
      av[4] = nv[48+1*3+1];
      av[5] = nv[48+1*3+2];
      
      av[6] = nv[48+2*3+0];
      av[7] = nv[48+2*3+1];
      av[8] = nv[48+2*3+2];

      av[9]  = nv[48+3*3+0];
      av[10] = nv[48+3*3+1];
      av[11] = nv[48+3*3+2];
      
      locglob_nodeval (4,av,tnv,x,y,z);
      
      av[0] = tnv[0*3+0];
      av[1] = tnv[0*3+1];
      av[2] = tnv[0*3+2];
      
      av[3] = tnv[1*3+0];
      av[4] = tnv[1*3+1];
      av[5] = tnv[1*3+2];
      
      av[6] = tnv[2*3+0];
      av[7] = tnv[2*3+1];
      av[8] = tnv[2*3+2];

      av[9] = tnv[3*3+0];
      av[10] = tnv[3*3+1];
      av[11] = tnv[3*3+2];

    }

    mxv (am,av,v);
    addv (v,nf,nf);
  }


  //  surface number 6

  if (is[5]>0 ){
    zeta=-1.0;
    for (i=0;i<intordb;i++){
      xi=gp[i];
      for (j=0;j<intordb;j++){
	eta=gp[j];
	
	jac2d_3d (jac,x,y,z,xi,eta,5);
	bf_matrix (n,xi,eta,zeta);
	jac = jac*w[i]*w[j];
	nnj (am.a,n.a,jac,n.m,n.n);
      }
    }
    
    if (is[5]==1){
      av[12] = nv[60+0*3+0];
      av[13] = nv[60+0*3+1];
      av[14] = nv[60+0*3+2];
      
      av[21] = nv[60+1*3+0];
      av[22] = nv[60+1*3+1];
      av[23] = nv[60+1*3+2];
      
      av[18] = nv[60+2*3+0];
      av[19] = nv[60+2*3+1];
      av[20] = nv[60+2*3+2];

      av[15] = nv[60+3*3+0];
      av[16] = nv[60+3*3+1];
      av[17] = nv[60+3*3+2];
    }


    if (is[5]==2){
      
      av[0] = nv[60+0*3+0];
      av[1] = nv[60+0*3+1];
      av[2] = nv[60+0*3+2];
      
      av[3] = nv[60+1*3+0];
      av[4] = nv[60+1*3+1];
      av[5] = nv[60+1*3+2];
      
      av[6] = nv[60+2*3+0];
      av[7] = nv[60+2*3+1];
      av[8] = nv[60+2*3+2];

      av[9]  = nv[60+3*3+0];
      av[10] = nv[60+3*3+1];
      av[11] = nv[60+3*3+2];
      
      locglob_nodeval (5,av,tnv,x,y,z);
      
      av[12] = tnv[0*3+0];
      av[13] = tnv[0*3+1];
      av[14] = tnv[0*3+2];
      
      av[21] = tnv[1*3+0];
      av[22] = tnv[1*3+1];
      av[23] = tnv[1*3+2];
      
      av[18] = tnv[2*3+0];
      av[19] = tnv[2*3+1];
      av[20] = tnv[2*3+2];

      av[15] = tnv[3*3+0];
      av[16] = tnv[3*3+1];
      av[17] = tnv[3*3+2];

    }

    mxv (am,av,v);
    addv (v,nf,nf);
  }






  destrv (gp);  destrv (w);
  delete [] tnv;
}

void linhex::locglob_nodeval (long is,vector &nv,double *tnv,vector &x,vector &y,vector &z)
{
  double norm;
  vector g1(3),g2(3),g3(3),lv(3),gv(3);
  matrix t(3,3);
  
  if (is==0){
    g1[0]=x[4]-x[7];
    g1[1]=y[4]-y[7];
    g1[2]=z[4]-z[7];
    
    g2[0]=x[3]-x[7];
    g2[1]=y[3]-y[7];
    g2[2]=z[3]-z[7];
  }
  if (is==1){
    g1[0]=x[5]-x[4];
    g1[1]=y[5]-y[4];
    g1[2]=z[5]-z[4];
    
    g2[0]=x[0]-x[4];
    g2[1]=y[0]-y[4];
    g2[2]=z[0]-z[4];
  }
  if (is==2){
    g1[0]=x[6]-x[5];
    g1[1]=y[6]-y[5];
    g1[2]=z[6]-z[5];
    
    g2[0]=x[1]-x[5];
    g2[1]=y[1]-y[5];
    g2[2]=z[1]-z[5];
  }
  if (is==3){
    g1[0]=x[7]-x[6];
    g1[1]=y[7]-y[6];
    g1[2]=z[7]-z[6];
    
    g2[0]=x[2]-x[6];
    g2[1]=y[2]-y[6];
    g2[2]=z[2]-z[6];
  }
  if (is==4){
    g1[0]=x[3]-x[2];
    g1[1]=y[3]-y[2];
    g1[2]=z[3]-z[2];
    
    g2[0]=x[1]-x[2];
    g2[1]=y[1]-y[2];
    g2[2]=z[1]-z[2];
  }
  if (is==5){
    g1[0]=x[5]-x[6];
    g1[1]=y[5]-y[6];
    g1[2]=z[5]-z[6];
    
    g2[0]=x[7]-x[6];
    g2[1]=y[7]-y[6];
    g2[2]=z[7]-z[6];
  }
  
  
  scprd (g1,g1,norm);
  norm=1.0/sqrt(norm);
  cmulv (norm,g1,g1);
  
  scprd (g1,g2,norm);
  g2[0]=g2[0]-norm*g1[0];
  g2[1]=g2[1]-norm*g1[1];
  g2[2]=g2[2]-norm*g1[2];
  
  scprd (g2,g2,norm);
  norm=1.0/sqrt(norm);
  cmulv (norm,g2,g2);
  
  g3[0]=g1[1]*g2[2]-g1[2]*g2[1];
  g3[1]=g1[2]*g2[0]-g1[0]*g2[2];
  g3[2]=g1[0]*g2[1]-g1[1]*g2[0];
  
  t[0][0]=g1[0];
  t[1][0]=g1[1];
  t[2][0]=g1[2];

  t[0][1]=g2[0];
  t[1][1]=g2[1];
  t[2][1]=g2[2];

  t[0][2]=g3[0];
  t[1][2]=g3[1];
  t[2][2]=g3[2];
  
  mxv (t.a,nv.a,tnv,3,3);
  mxv (t.a,nv.a+3,tnv+3,3,3);
  mxv (t.a,nv.a+6,tnv+6,3,3);
  mxv (t.a,nv.a+9,tnv+9,3,3);
  
}

/**
   function computes nodal forces caused by presure on surface
   
   @param eid - element id
   @param ri,ci - row and column indices
   @param nfor - vector of presure 
   @param eis - surface id 
   
   4.2002, PF
*/
void linhex::node_forces_surf_old (long lcid,long eid,long *is,double *nv,vector &nf)
{
  long i,j,i1,ii;
  double xi=0.0,eta=0.0,zeta=0.0,jac, w1,w2;
  ivector nodes(nne);
  vector gx(nne),gy(nne),gz(nne),x(nnsurf),y(nnsurf),z(nnsurf),gp(intordb),w(intordb),av(ndofe),v(ndofe),pom(3);
  matrix n(napfun,ndofe),an(napfun,ndofe),am(ndofe,ndofe), tran(3,3);

  Mt->give_elemnodes (eid,nodes);
  Mt->give_node_coord3d (gx,gy,gz,eid);
  gauss_points (gp.a,w.a,intordb);
  fillv (0.0,nf);
  fillm (0.0,an);
  for (ii=0;ii<6;ii++){
    fillv (0.0,av);
    // is=0 not loading
    if (is[ii] ==0){}
    else {
      tran_mat( tran, gx, gy, gz, ii);
      // is=1 surface node 1,4,8,5
      if (ii ==0){
	xi=1.0;
        x[0]=gx[0];x[1]=gx[3];x[2]=gx[7];x[3]=gx[4]; 
        y[0]=gy[0];y[1]=gy[3];y[2]=gy[7];y[3]=gy[4]; 
        z[0]=gz[0];z[1]=gz[3];z[2]=gz[7];z[3]=gz[4]; 
	for (i=0;i<intordb;i++){
	  eta=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    zeta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,eta,zeta);
	    jac = jac*w1*w2;     
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[i]=nv[i];  av[9+i]=nv[3+i];  av[21+i]=nv[6+i];  av[12+i]=nv[9+i];
	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=2 surface node 2,1,5,6
      else if (ii ==1){
	eta=1.0;
        x[0]=gx[1];x[1]=gx[0];x[2]=gx[4];x[3]=gx[5];
        y[0]=gy[1];y[1]=gy[0];y[2]=gy[4];y[3]=gy[5]; 
        z[0]=gz[1];z[1]=gz[0];z[2]=gz[4];z[3]=gz[5]; 
	for (i=0;i<intordb;i++){
	  xi=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    zeta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,xi,zeta);
	    jac = jac*w1*w2;      
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[3+i]=nv[12+i];  av[i]=nv[15+i];  av[12+i]=nv[18+i];  av[15+i]=nv[21+i];
	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=3 surface node 3,2,6,7
      else if (ii ==2){
	xi=-1.0;
        x[0]=gx[2];x[1]=gx[1];x[2]=gx[5];x[3]=gx[6]; 
        y[0]=gy[2];y[1]=gy[1];y[2]=gy[5];y[3]=gy[6]; 
        z[0]=gz[2];z[1]=gz[1];z[2]=gz[5];z[3]=gz[6]; 
	for (i=0;i<intordb;i++){
	  eta=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    zeta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,eta,zeta);
	    jac = jac*w1*w2;      
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[6+i]=nv[24+i];  av[3+i]=nv[27+i];  av[15+i]=nv[30+i];  av[18+i]=nv[33+i];
	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=4 surface node 4,3,7,8
      else if (ii ==3){
	eta=-1.0;
        x[0]=gx[3];x[1]=gx[2];x[2]=gx[6];x[3]=gx[7]; 
        y[0]=gy[3];y[1]=gy[2];y[2]=gy[6];y[3]=gy[7];
        z[0]=gz[3];z[1]=gz[2];z[2]=gz[6];z[3]=gz[7]; 
	for (i=0;i<intordb;i++){
	  xi=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    zeta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,xi,zeta);
	    jac = jac*w1*w2;
	    
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[9+i]=nv[36+i];  av[6+i]=nv[39+i];  av[18+i]=nv[42+i];  av[21+i]=nv[45+i];
	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=5  surface node 1,2,3,4
      else if (ii ==4) {
	zeta=1.0;
        x[0]=gx[0];x[1]=gx[1];x[2]=gx[2];x[3]=gx[3]; 
        y[0]=gy[0];y[1]=gy[1];y[2]=gy[2];y[3]=gy[3]; 
        z[0]=gz[0];z[1]=gz[1];z[2]=gz[2];z[3]=gz[3]; 
	for (i=0;i<intordb;i++){
	  xi=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    eta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,xi,eta);
	    jac = jac*w1*w2;      
	    //	for constant			for (i1=0;i1<ndofe;i1++){
	    //								for (j1=0;j1<3;j1++){
	    //									an[j1][i1]=an[j1][i1]+n[j1][i1]*jac;}	}
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[i]=nv[48+i];  av[3+i]=nv[51+i];  av[6+i]=nv[54+i];  av[9+i]=nv[57+i];
	}
	//		fprintf (Out,"\n\n nn");
	//		for (i1=0;i1<ndofe;i1++){
	//			fprintf (Out,"\n %4ld   %20.10e %20.10e %20.10e",i1,an[0][i1],an[1][i1],an[2][i1]);}       
	
	// for constant		for (i=0;i<ndofe;i++){
	//						for (j=0;j<pom.n;j++){
	//							v[i]=an[j][i]*pom[j];}	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=6  surface node 6,5,8,7
      else if (ii ==5) {
	zeta=-1.0;
        x[0]=gx[5];x[1]=gx[4];x[2]=gx[7];x[3]=gx[6]; 
        y[0]=gy[5];y[1]=gy[4];y[2]=gy[7];y[3]=gy[6]; 
        z[0]=gz[5];z[1]=gz[4];z[2]=gz[7];z[3]=gz[6]; 
	for (i=0;i<intordb;i++){
	  xi=gp[i];  w1=w[i];      
	  for (j=0;j<intordb;j++){
	    eta=gp[j];  w2=w[j];
	    bf_matrix (n,xi,eta,zeta);
	    jac2d3d (jac,x,y,z,xi,eta);
	    jac = jac*w1*w2;      
	    nnj (am.a,n.a,jac,n.m,n.n);
	  }
	}
	for (i=0;i<3;i++){
	  av[12+i]=nv[60+i];  av[15+i]=nv[63+i];  av[18+i]=nv[66+i];  av[21+i]=nv[69+i];
	}
	// Load in LCS and transformation to GCS
	if (is[ii] ==2){
	  lgvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);  
      }
      
			fprintf (Out,"\n\n zatiz");
			for (i=0;i<nne;i++){
                i1=3*i;
				fprintf (Out,"\n %4ld   %20.10e %20.10e %20.10e",i,av[i1],av[i1+1],av[i1+2]);}       
			fprintf (Out,"\n\n vloc");
			for (i=0;i<nne;i++){
                i1=3*i;
				fprintf (Out,"\n %4ld   %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}       
//      lgvectortransf3dblock(v, tran);
			fprintf (Out,"\n\n vglob");
			for (i=0;i<nne;i++){
                i1=3*i;
				fprintf (Out,"\n %4ld   %20.10e %20.10e %20.10e",i,v[i1],v[i1+1],v[i1+2]);}       
      addv (nf,v,nf);      
    }
  }

  //  transformation of stiffness matrix to nodesystem
  //  transf = Mt->locsystems (nodes);
  //  if (transf>0){
  //    matrix tmat (ndofe,ndofe);
  //    transf_matrix (nodes,tmat);
  //    lgvectortransf3dblock (v,tmat);
  //  }

  //  fprintf (Out,"\n\n zatizeni na prvku cislo %ld",eid);
  //  for (i=0;i<ndofe;i++){
  //	  fprintf (Out,"\n %4ld   %20.10e",i,nf[i]);}
}

/**
   function computes transformation matrix on surface
   
   @param x, y - local coordinate  xL=adge12
   @param tran - tranformation metrix to GCS
   @param gx,gy,gz - vector of node global coordinate 
   @param is - surface id 
   
   4.2002, PF
*/
void linhex::tran_mat(matrix &tran, vector &gx, vector &gy, vector &gz, long is)
{
  long i,i1;
  double dl;
  matrix a(3,3);
// is=5  surface node 1,2,3,4
  if (is ==4) {
    for (i=0; i<3; i++) {
      i1=i+1; if(i1>2)i1=i1-3;
      a[0][i]=gx[i1]-gx[i];
      a[1][i]=gy[i1]-gy[i];
      a[2][i]=gz[i1]-gz[i];
    }
  }
// is=6  surface node 6,5,8,7
  else if (is ==5) {
    for (i=4; i<7; i++) {
      i1=i+1; if(i1>6)i1=i1-3;
      a[0][i-4]=gx[i1]-gx[i];
      a[1][i-4]=gy[i1]-gy[i];
      a[2][i-4]=gz[i1]-gz[i];
    }
  }
// is=1 surface node 1,4,8,5
// is=2 surface node 2,1,5,6
// is=3 surface node 3,2,6,7
// is=4 surface node 4,3,7,8
  else {
      if(is>0)i1=is-4;else i1=is;
      a[0][0]=gx[i1+3]-gx[is];
      a[1][0]=gy[i1+3]-gy[is];
      a[2][0]=gz[i1+3]-gz[is];
      a[0][1]=gx[i1+7]-gx[i1+3];
      a[1][1]=gy[i1+7]-gy[i1+3];
      a[2][1]=gz[i1+7]-gz[i1+3];
      a[0][2]=gx[is]-gx[i1+7];
      a[1][2]=gy[is]-gy[i1+7];
      a[2][2]=gz[is]-gz[i1+7];
  }
  
  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);
//    x[3]=(gx(3)-gx(0))*tran(0,0)+(gy(3)-gy(0))*tran(1,0)+(gz(3)-gz(0))*tran(2,0);
//    y[3]=(gx(3)-gx(0))*tran(0,1)+(gy(3)-gy(0))*tran(1,1)+(gz(3)-gz(0))*tran(2,1);
}


/**
   function interpolates the nodal values to the integration points on the element
   
   @param eid - element id
   
   17.8.2004, JK
*/
void linhex::intpointval (long eid,vector &nodval,vector &ipval)
{
  long i,j,ii,jj,k,l;
  double xi,eta,zeta;
  vector w,gp;
  
  l=0;