quadhex.cpp

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

/**
   function returns coordinates of integration points
   
   @param eid - element id
   @param ipp - integration point pointer
   @param ri - row index
   @param ci - column index
   @param coord - vector of coordinates

   10.1.2002
*/
void quadhex::ipcoord (long eid,long ipp,long ri,long ci,vector &coord)
{
  long i,j,k,ii;
  double xi,eta,zeta;
  vector x(nne),y(nne),z(nne),w(intordsm[ri][ci]),gp(intordsm[ri][ci]);
  
  gauss_points (gp.a,w.a,intordsm[ri][ci]);
  Mt->give_node_coord3d (x,y,z,eid);
  ii=Mt->elements[eid].ipp[ri][ci];
  
  for (i=0;i<intordsm[ri][ci];i++){
    xi=gp[i];
    for (j=0;j<intordsm[ri][ci];j++){
      eta=gp[j];
      for (k=0;k<intordsm[ri][ci];k++){
	zeta=gp[k];
	if (ii==ipp){
	  coord[0]=approx (xi,eta,zeta,x);
	  coord[1]=approx (xi,eta,zeta,y);
	  coord[2]=approx (xi,eta,zeta,z);
	}
	ii++;
      }
    }
  }
}



void quadhex::inicipval(long eid, long ri, long ci, matrix &nodval, inictype *ictn)
{
  long i, j, k, l, m, ipp;
  long ii, jj, nv = nodval.n;
  long nstra;
  double xi, eta, zeta, ipval;
  vector w, gp, anv(nne);

  nstra = 0;
  for (j = 0; j < nv; j++) // for all initial values
  {
    for(i = 0; i < nne; i++) // for all nodes on element
      anv[i] = nodval[i][j];
    for (ii = 0; ii < nb; ii++)
    {
      for (jj = 0; jj < nb; jj++)
      {
        ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
        if (intordsm[ii][jj] == 0)
          continue;
        allocv (intordsm[ii][jj],gp);
        allocv (intordsm[ii][jj],w);
        gauss_points (gp.a,w.a,intordsm[ii][jj]);
        for (k = 0; k < intordsm[ii][jj]; k++)
        {
          xi=gp[k];
          for (l = 0; l < intordsm[ii][jj]; l++)
          {
            eta=gp[l];
            for (m = 0; m < intordsm[ii][jj]; m++)
            {
              zeta=gp[m];
              //  value in integration point
              ipval = approx (xi,eta,zeta,anv);
              if ((ictn[i] & inistrain) && (j < Mm->ip[ipp].ncompstr))
              {
                Mm->ip[ipp].strain[j] += ipval;
                ipp++;
                continue;
              }
              if ((ictn[i] & inistress) && (j < nstra + Mm->ip[ipp].ncompstr))
              {
                Mm->ip[ipp].stress[j] += ipval;
                ipp++;
                continue;
              }
              if ((ictn[i] & iniother) && (j < nv))
              {
                Mm->ip[ipp].other[j] += ipval;
                ipp++;
                continue;
              }
              ipp++;
            }
          }
        }
        destrv (gp);  destrv (w);
      }
    }
    if (ictn[i] & inistrain) nstra++;
  }
}

/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void quadhex::ipvolume (long eid,long ri,long ci)
{
  long i,j,k,ii,jj,ipp;
  double xi,eta,zeta,jac;
  vector x(nne),y(nne),z(nne),w,gp;
  
  Mt->give_node_coord3d (x,y,z,eid);
  
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;

      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);
      
      gauss_points (gp.a,w.a,intordsm[ii][jj]);
      
      ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
      
      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];
	  for (k=0;k<intordsm[ii][jj];k++){
	    zeta=gp[k];
	    
	    jac_3d (jac,x,y,z,xi,eta,zeta);
	    jac=fabs(jac);
	    
	    jac*=w[i]*w[j]*w[k];
	    
	    Mm->storeipvol (ipp,jac);
	    ipp++;
	  }
	}
      }
      destrv (gp);  destrv (w);
    }
  }
  
}

/**
   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 quadhex::node_forces_surf (long lcid,long eid,long *is,double *nv,vector &nf)
{
  long i,j,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,  12,16,20,13
      if (ii ==0){
	xi=1.0;
        x[0]=gx[0];x[1]=gx[3];x[2]=gx[7];x[3]=gx[4]; x[4]=gx[11];x[5]=gx[15];x[6]=gx[18];x[7]=gx[12];
        y[0]=gy[0];y[1]=gy[3];y[2]=gy[7];y[3]=gy[4]; y[4]=gy[11];y[5]=gy[15];y[6]=gy[18];y[7]=gy[12];
        z[0]=gz[0];z[1]=gz[3];z[2]=gz[7];z[3]=gz[4]; z[4]=gz[11];z[5]=gz[15];z[6]=gz[18];z[7]=gz[12];
	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[15+i]=nv[9+i];
	  av[33+i]=nv[12+i];  av[45+i]=nv[15+i];  av[57+i]=nv[18+i];  av[36+i]=nv[21+i];
	}
	// Load in GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=2 surface node 2,1,5,6,  9,13,17,14
      else if (ii ==1){
	eta=1.0;
        x[0]=gx[1];x[1]=gx[0];x[2]=gx[4];x[3]=gx[5]; x[4]=gx[8];x[5]=gx[12];x[6]=gx[16];x[7]=gx[13];
        y[0]=gy[1];y[1]=gy[0];y[2]=gy[4];y[3]=gy[5]; y[4]=gy[8];y[5]=gy[12];y[6]=gy[16];y[7]=gy[13];
        z[0]=gz[1];z[1]=gz[0];z[2]=gz[4];z[3]=gz[5]; z[4]=gz[8];z[5]=gz[12];z[6]=gz[16];z[7]=gz[13];
	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[24+i];  av[i]=nv[27+i];  av[12+i]=nv[30+i];  av[15+i]=nv[33+i];
	  av[24+i]=nv[36+i];  av[36+i]=nv[39+i];  av[48+i]=nv[42+i];  av[39+i]=nv[45+i];
	}
	// Load in GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=3 surface node 3,2,6,7, 10, 14,18,15
      else if (ii ==2){
	xi=-1.0;
        x[0]=gx[2];x[1]=gx[1];x[2]=gx[5];x[3]=gx[6]; x[4]=gx[9];x[5]=gx[13];x[6]=gx[17];x[7]=gx[14];
        y[0]=gy[2];y[1]=gy[1];y[2]=gy[5];y[3]=gy[6]; y[4]=gy[9];y[5]=gy[13];y[6]=gy[17];y[7]=gy[14];
        z[0]=gz[2];z[1]=gz[1];z[2]=gz[5];z[3]=gz[6]; z[4]=gz[9];z[5]=gz[13];z[6]=gz[17];z[7]=gz[14];
	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[48+i];  av[3+i]=nv[51+i];  av[15+i]=nv[54+i];  av[18+i]=nv[57+i];
	  av[27+i]=nv[60+i];  av[39+i]=nv[63+i];  av[51+i]=nv[66+i];  av[42+i]=nv[69+i];
	}
	// Load in GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=4 surface node 4,3,7,8,  11,15,19,16
      else if (ii ==3){
	eta=-1.0;
        x[0]=gx[3];x[1]=gx[2];x[2]=gx[6];x[3]=gx[7]; x[4]=gx[10];x[5]=gx[14];x[6]=gx[18];x[7]=gx[15];
        y[0]=gx[3];y[1]=gy[2];y[2]=gy[6];y[3]=gy[7]; y[4]=gy[10];y[5]=gy[14];y[6]=gy[18];y[7]=gy[15];
        z[0]=gx[3];z[1]=gz[2];z[2]=gz[6];z[3]=gz[7]; z[4]=gz[10];z[5]=gz[14];z[6]=gz[18];z[7]=gz[15];
	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[72+i];  av[6+i]=nv[75+i];  av[18+i]=nv[78+i];  av[21+i]=nv[81+i];
	  av[30+i]=nv[84+i];  av[42+i]=nv[87+i];  av[54+i]=nv[90+i];  av[45+i]=nv[93+i];
	}
	// Load in GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=5  surface node 1,2,3,4,  9,10,11,12
      else if (ii ==4) {
	zeta=1.0;
        x[0]=gx[0];x[1]=gx[1];x[2]=gx[2];x[3]=gx[3]; x[4]=gx[8];x[5]=gx[9];x[6]=gx[10];x[7]=gx[11];
        y[0]=gy[0];y[1]=gy[1];y[2]=gy[2];y[3]=gy[3]; y[4]=gy[8];y[5]=gy[9];y[6]=gy[10];y[7]=gy[11];
        z[0]=gz[0];z[1]=gz[1];z[2]=gz[2];z[3]=gz[3]; z[4]=gz[8];z[5]=gz[9];z[6]=gz[10];z[7]=gz[11];
	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[96+i];  av[3+i]=nv[99+i];  av[6+i]=nv[102+i];  av[9+i]=nv[105+i];
	  av[24+i]=nv[108+i];  av[27+i]=nv[111+i];  av[30+i]=nv[114+i];  av[33+i]=nv[117+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 GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);
      }
      // is=6  surface node 5,6,7,8, 17,18,19,20
      else if (ii ==5) {
	zeta=-1.0;
        x[0]=gx[4];x[1]=gx[5];x[2]=gx[6];x[3]=gx[7]; x[4]=gx[16];x[5]=gx[17];x[6]=gx[18];x[7]=gx[19];
        y[0]=gy[4];y[1]=gy[5];y[2]=gy[6];y[3]=gy[7]; y[4]=gy[16];y[5]=gy[17];y[6]=gy[18];y[7]=gy[19];
        z[0]=gz[4];z[1]=gz[5];z[2]=gz[6];z[3]=gz[7]; z[4]=gz[16];z[5]=gz[17];z[6]=gz[18];z[7]=gz[19];
	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[120+i];  av[15+i]=nv[123+i];  av[18+i]=nv[126+i];  av[21+i]=nv[129+i];
	  av[48+i]=nv[132+i];  av[51+i]=nv[135+i];  av[54+i]=nv[138+i];  av[57+i]=nv[141+i];
	}
	// Load in GCS
	if (is[ii] ==1){
	  glvectortransf3dblock (av, tran);
	}
        mxv (am,av,v);  
      }
      
      lgvectortransf3dblock(v, tran);
      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 quadhex::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 5,6,7,8
  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
   quadratic approximation functions are used
   
   @param eid - element id
   @param nodval - nodal values
   @param ipval - value at integration points
   
   17.8.2004, JK
*/
void quadhex::intpointval (long eid,vector &nodval,vector &ipval)
{
  long i,j,ii,jj,k,l;
  double xi,eta,zeta;
  vector w,gp;
  
  l=0;
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;
      
      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);

      gauss_points (gp.a,w.a,intordsm[ii][jj]);

      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];
	  for (k=0;k<intordsm[ii][jj];k++){
	    zeta=gp[k];
	    
	    ipval[l]=approx (xi,eta,zeta,nodval);
	    l++;
	  }
	}
      }
      
      destrv (w);
      destrv (gp);
    }
  }
}

/**
   function interpolates the nodal values to the integration points on the element
   linear approximation functions are used
   
   @param eid - element id
   @param nodval - nodal values
   @param ipval - value at integration points
   
   17.8.2004, JK
*/
void quadhex::intpointval2 (long eid,vector &nodval,vector &ipval)
{
  long i,j,ii,jj,k,l;
  double xi,eta,zeta;
  vector w,gp;
  vector modnodval(Lhex->nne);
  
  for (i=0;i<Lhex->nne;i++){
    modnodval[i]=nodval[i];
  }
  
  l=0;
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;
      
      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);
      
      gauss_points (gp.a,w.a,intordsm[ii][jj]);
      
      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];
	  for (k=0;k<intordsm[ii][jj];k++){
	    zeta=gp[k];
	    
	    ipval[l]=Lhex->approx (xi,eta,zeta,modnodval);
	    l++;
	  }
	}
      }
      
      destrv (w);
      destrv (gp);
    }
  }
}



void quadhex::aver_strains (long lcid,long eid,long ri,long ci,vector &averstra,double &volume)
{
  long i,j,k,l,ii,ipp;
  double xi,eta,zeta,jac;
  vector x(nne),y(nne),z(nne),gp,w,eps;

  Mt->give_node_coord3d (x,y,z,eid);
  
  for (ii=0;ii<nb;ii++){
    if (intordsm[ii][ii]==0)  continue;
    
    allocv (intordsm[ii][ii],gp);
    allocv (intordsm[ii][ii],w);
    allocv (ncomp[ii],eps);
    
    gauss_points (gp.a,w.a,intordsm[ii][ii]);
    
    ipp=Mt->elements[eid].ipp[ri+ii][ci+ii];
    for (i=0;i<intordsm[ii][ii];i++){
      xi=gp[i];
      for (j=0;j<intordsm[ii][ii];j++){
	eta=gp[j];
	for (k=0;k<intordsm[ii][ii];k++){
	  zeta=gp[k];
	  
	  //Mm->givestrain (lcid,ipp,eps);
	  Mm->givestress (lcid,ipp,eps);
	  
	  
	  volume+=w[i]*w[j]*w[k]*jac;
	  
	  for (l=0;l<averstra.n;l++){
	    averstra[l]+=eps[l]*w[i]*w[j]*w[k]*jac;
	  }
	  
	  ipp++;
	}
      }
    }
        

    destrv (w);  destrv (gp);  destrv (eps);
  }
  
}