axiquad.c

工程中有限元程序,采用C语言编制,包括所有经典的有限元问题!

	   fprintf(fptr1, " %10.4e ", sv);
	   if (ipl == 2)
	      fprintf(fptr2, " %f ", sv);
       }   
	   fprintf(fptr1, "\n");
	   if (ipl > 1)
	      fprintf(fptr2, "\n");
     }
     fclose(fptr1);
     printf("complete results are in file %s\n", file2);
     printf("view using a text processor\n");
     if (ipl > 1) {
        fclose(fptr2);
        printf("element stress data in file %s\n", file3);
        printf("run bestfitq and then contourA or contourB to plot stresses\n");
        }
     return(0);
   }
integ(xni)
  float xni[][2];
  {
    float c;
   /* ----- integration points xni() ----- */
     c = .57735026919;
     xni[0][0] = -c;
     xni[0][1] = -c;
     xni[1][0] = c;
     xni[1][1] = -c;
     xni[2][0] = c;
     xni[2][1] = c;
     xni[3][0] = -c;
     xni[3][1] = c;
     return(0);
  }

dmatrix(n,pm,mat,npr,al1,d)
   int n,npr,*mat;
   float *pm,*al1,d[][4];
  {
     int m;
     float e,c,c1,c2,c3,pnu,al;
   /* -----  d() matrix  ----- */
     /* --- material properties --- */
     m = mat[n]-1;
     e = pm[npr*m];
     pnu= pm[npr*m+1];
     al = pm[npr*m+2];
     *al1 = al;
     /* --- d() matrix --- */
     c1 = e * (1 - pnu) / ((1 + pnu) * (1 - 2 * pnu));
     c2 = pnu / (1 - pnu);
     d[0][0] = c1;
     d[0][1] = c1*c2;
     d[0][2] = 0;
     d[0][3] = c1*c2;
     d[1][0] = c1*c2;
     d[1][1] = c1;
     d[1][2] = 0;
     d[1][3] = c1*c2;
     d[2][0] = 0;
     d[2][1] = 0;
     d[2][2] = .5*e/(1+pnu);
     d[2][3] = 0;
     d[3][0] = c1*c2;
     d[3][1] = c1*c2;
     d[3][2] = 0;
     d[3][3] = c1;
     return(0);
  }

elstif(n,se,tl,xni,d,tempr,x,al,noc)
    int n,*noc;
    float al;
    float *x,*tempr,d[][4],tl[8],se[][8],xni[][2];
  {
    int i,j,k,ip;
    float dte,c,xi,eta,dj,rad,b[4][8],db[4][8];
    float pi = 3.14159;
    /* -----  element stiffness and temperature load  ----- */
     for (i = 0; i < 8; i++) {
	 for (j = 0; j < 8; j++) {
	     se[i][j] = 0.;
	     }
	 tl[i] = 0.;
	 }
     dte = tempr[n];
     /* --- weight factor is one --- */
     /* --- loop on integration points --- */
     for (ip = 0; ip < 4; ip++) {
        /* ---  get db matrix at integration point ip --- */
        xi = xni[ip][0];
	eta = xni[ip][1];
	dbmat(n,x,noc,d,b,db,&dj,&rad,xi,eta);
	/* --- element stiffness matrix  se --- */
	for (i = 0; i < 8; i++) {
           for (j = 0; j < 8; j++) {
              c = 0;
              for (k = 0; k < 4; k++) {
                 c = c + 2 * pi * rad * b[k][i] * db[k][j] * dj;
                 }
              se[i][j] = se[i][j] + c;
              }
	   }
	/* --- determine temperature load tl --- */
	c = 2 * pi * rad * dj * al * dte;
        for (i = 0; i < 8; i++) {
           tl[i] = tl[i] + c * (db[0][i] + db[1][i]);
           }
	}
     return(0);
   }
dbmat(n,x,noc,d,b,db,dj1,rad1,xi,eta)
  float *x,*dj1,*rad1,xi,eta;
  float d[][4],b[][8],db[][8];
  int n,*noc;
  {
   int n1,n2,n3,n4,i,j,k;
   float x1,y1,x2,y2,x3,y3,x4,y4,tj11,tj12,tj21,tj22,dj,c;
   float rad,sh1,sh2,sh3,sh4;
   float a[3][4],g[4][8];
   /* -----  db()  matrix  ----- */
     /* --- nodal coordinates --- */     
     n1 = noc[4*n];
     n2 = noc[4*n+1];
     n3 = noc[4*n+2];
     n4 = noc[4*n+3];
     x1 = x[2*(n1-1)];
     y1 = x[2*(n1-1)+1];
     x2 = x[2*(n2-1)];
     y2 = x[2*(n2-1)+1];
     x3 = x[2*(n3-1)];
     y3 = x[2*(n3-1)+1];
     x4 = x[2*(n4-1)];
     y4 = x[2*(n4-1)+1];
     /* --- formation of jacobian  tj --- */
     tj11 = ((1 - eta) * (x2 - x1) + (1 + eta) * (x3 - x4)) / 4;
     tj12 = ((1 - eta) * (y2 - y1) + (1 + eta) * (y3 - y4)) / 4;
     tj21 = ((1 - xi) * (x4 - x1) + (1 + xi) * (x3 - x2)) / 4;
     tj22 = ((1 - xi) * (y4 - y1) + (1 + xi) * (y3 - y2)) / 4;
     /* --- determinant of the jacobian --- */
     dj = tj11 * tj22 - tj12 * tj21;
     *dj1 = dj;
     /* --- a[3,4] matrix relates strains to --- */
     /* --- local derivatives of u --- */
     a[0][0] = tj22 / dj;
     a[1][0] = 0;
     a[2][0] = -tj21 / dj;
     a[0][1] = -tj12 / dj;
     a[1][1] = 0;
     a[2][1] = tj11 / dj;
     a[0][2] = 0;
     a[1][2] = -tj21 / dj;
     a[2][2] = tj22 / dj;
     a[0][3] = 0;
     a[1][3] = tj11 / dj;
     a[2][3] = -tj12 / dj;
     /* --- g[4,8] matrix relates local derivatives of u --- */
     /* --- to local nodal displacements q[8] --- */
     for (i = 0; i < 4; i++) {
	     for (j = 0; j < 8; j++) {
             g[i][j] = 0;
	         }
	     }
     g[0][0] = -(1 - eta) / 4;
     g[1][0] = -(1 - xi) / 4;
     g[2][1] = -(1 - eta) / 4;
     g[3][1] = -(1 - xi) / 4;
     g[0][2] = (1 - eta) / 4;
     g[1][2] = -(1 + xi) / 4;
     g[2][3] = (1 - eta) / 4;
     g[3][3] = -(1 + xi) / 4;
     g[0][4] = (1 + eta) / 4;
     g[1][4] = (1 + xi) / 4;
     g[2][5] = (1 + eta) / 4;
     g[3][5] = (1 + xi) / 4;
     g[0][6] = -(1 + eta) / 4;
     g[1][6] = (1 - xi) / 4;
     g[2][7] = -(1 + eta) / 4;
     g[3][7] = (1 - xi) / 4;
     /* --- Shape Function Values --- */
     sh1 = .25 * (1 - xi) * (1 - eta);
     sh2 = .25 * (1 + xi) * (1 - eta);
     sh3 = .25 * (1 + xi) * (1 + eta);
     sh4 = .25 * (1 - xi) * (1 + eta);
     rad = sh1 * x1 + sh2 * x2 + sh3 * x3 + sh4 * x4;
     *rad1 = rad;
     /* --- b[4,8] matrix relates strains to q --- */
     for (i = 0; i < 3; i++) {
        for (j = 0; j < 8; j++) {
           c = 0;
           for (k = 0; k < 4; k++) {
               c = c + a[i][k] * g[k][j];
               }
           b[i][j] = c;
           }
        }
     b[3][0] = sh1 / rad; 
     b[3][1] = 0;
     b[3][2] = sh2 / rad; 
     b[3][3] = 0;
     b[3][4] = sh3 / rad;
     b[3][5] = 0;
     b[3][6] = sh4 / rad; 
     b[3][7] = 0;
     /* --- db[3,8] matrix relates stresses to q[8] --- */
     for (i = 0; i < 4; i++) {
        for (j = 0; j < 8; j++) {
           c = 0;
           for (k = 0; k < 4; k++) {
               c = c + d[i][k] * b[k][j];
               }
	   db[i][j] = c;
	   }
	}
     return(0);
   }

/* ----- band solver ----- */
bansol(s,f,nq,nbw)
  int nq, nbw;
  float *s, *f;
{
 int n1,k,nk,i,i1,j,j1,kk;
  float c1;
  /* ----- band solver ----- */
  n1 = nq - 1;
  /* --- forward elimination --- */
  for (k = 1; k <= n1; k++) {
     nk = nq - k + 1;
     if (nk > nbw)
	nk = nbw;
     for (i = 2; i <= nk; i++) {
       c1 = s[nbw*(k-1)+i-1] / s[nbw*(k-1)];
       i1 = k + i - 1;
       for (j = i; j <= nk; j++) {
	j1 = j - i + 1;
	s[nbw*(i1-1)+j1-1] = s[nbw*(i1-1)+j1-1] - c1 * s[nbw*(k-1)+j-1];
	}
       f[i1-1] = f[i1-1] - c1 * f[k-1];
       }
     }
  /* --- back-substitution --- */
  f[nq-1] = f[nq-1] / s[nbw*(nq-1)];
  for (kk = 1; kk <= n1;kk++) {
     k = nq - kk;
     c1 = 1 / s[nbw*(k-1)];
     f[k-1] = c1 * f[k-1];
     nk = nq - k + 1;
     if (nk > nbw)
       nk = nbw;
       for (j = 2; j <= nk; j++) {
	 f[k-1] = f[k-1] - c1 * s[nbw*(k-1)+j-1] * f[k + j - 2];
	}
     }
    return(0);
}