linhex_nb1.cpp

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

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

      av[18] = nv[24+3*3+0];
      av[19] = nv[24+3*3+1];
      av[20] = nv[24+3*3+2];
    }


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

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

      av[18] = tnv[3*3+0];
      av[19] = tnv[3*3+1];
      av[20] = tnv[3*3+2];

    }

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


  //  surface number 4

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

      av[21] = nv[36+3*3+0];
      av[22] = nv[36+3*3+1];
      av[23] = nv[36+3*3+2];
    }


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

      av[9]  = nv[36+3*3+0];
      av[10] = nv[36+3*3+1];
      av[11] = nv[36+3*3+2];
      
      locglob_nodeval (3,av,tnv,x,y,z);
      
      av[9] = tnv[0*3+0];
      av[10]= tnv[0*3+1];
      av[11]= tnv[0*3+2];
      
      av[6] = tnv[1*3+0];
      av[7] = tnv[1*3+1];
      av[8] = tnv[1*3+2];
      
      av[18] = tnv[2*3+0];
      av[19] = tnv[2*3+1];
      av[20] = tnv[2*3+2];

      av[21] = tnv[3*3+0];
      av[22] = tnv[3*3+1];
      av[23] = tnv[3*3+2];

    }

    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