hexafron.c

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

/*****        program hexafron         *****/
/*   3-d stress analysis using  8-node     */
/*    isoparametric hexahedral element     */
/*          using frontal solver           */
/*   t.r.chandrupatla and a.d.belegundu    */
/*******************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
struct data
  {
  int variable;
  double coefft;
  };
  FILE *fptr;
main()
{
   FILE *fptr1;
   long int icount;
   int n,i,j,k,m,nfron,ntogo,ndcnt,in,ii;
   char dummy[81], title[81], file1[81], file2[81];
   int ne,nn,nq,nm,nd,nl,nen,ndn,ndim,npr,nmpc,ibl;
   int mtn,mtn1,ip;
   int *noc, *nu, *mat, *mpc, *isbl, *iebl, *indx;
   double *x, *pm, *u, *tempr, *s, *f, *beta;
   double c,dj,al,tld,cnst,reaction,s1,s2,s3,pi,cal,siv1,siv2,vm;
   double se[24][24],xi[3][8],xni[3][8],d[6][6];
   double b[6][24],db[6][24],qt[24],str[6];
/*-------------------------------------------------------*/
   printf("\n");
   puts("Input file name < dr:fn.ext >: ");
   gets(file1);
   puts("Output file name < dr:fn.ext >: ");
   gets(file2);
   printf("\n");
   fptr1 = fopen(file1, "r");
   fgets(dummy,80,fptr1);
   fgets(title,80,fptr1);
   fgets(dummy,80,fptr1);
   fscanf(fptr1,"%d %d %d %d %d %d\n", &nn, &ne, &nm, &ndim, &nen, &ndn);
   npr = 3;    /* Material properties E, Nu, Alpha */
   nq = nn * ndn;
   fgets(dummy, 80, fptr1);
   fscanf(fptr1,"%d %d %d %d %d\n", &nd, &nl, &nmpc);
/* ----- memory allocation ----- */
   x = (double *) calloc(nn*ndim, sizeof(double));
   noc = (int *) calloc(ne*nen, sizeof(int));
   u = (double *) calloc(nd, sizeof(double));
   nu = (int *) calloc(nd, sizeof(int));
   mat = (int *) calloc(ne,sizeof(int));
   f = (double *) calloc(nn*ndn, sizeof(double));
   tempr = (double *) calloc(ne, sizeof(double));
   pm = (double *) calloc(nm*npr, sizeof(double));
   mpc = (int *) calloc(2*nmpc, sizeof(int));
   beta = (double *) calloc(3*nmpc, sizeof(double));
/* ----- total dof is  nq ----- */
     nq = ndn * nn;
/* ===============  read data  ==================== */
/* ----- coordinates ----- */
     fgets(dummy,80,fptr1);
     for (i = 0; i < nn; i++){
        fscanf(fptr1, "%d", &n);
        for (j = 0; j < ndim; j++){
		fscanf(fptr1, "%lf\n", &c);
                x[ndim*(n-1)+j] = c;
            }
         }
/* ----- connectivity, material, temp-change ----- */
   fgets(dummy,80,fptr1);
   for (i = 0; i < ne; i++) {
       fscanf(fptr1,"%d", &n);
       for (j = 0; j < nen; j++) {
           fscanf(fptr1,"%d", &k);
           noc[(n-1)*nen+j]=k;
       }
       fscanf(fptr1,"%d", &k);
       mat[n-1] = k;
       fscanf(fptr1,"%lf\n",&c);
       tempr[n-1] = c;
   }
/* ----- displacement bc  ----- */
   fgets(dummy,80,fptr1);
   for (i = 0; i < nd; i++) {
      fscanf(fptr1, "%d %lf\n", &k, &c);
      nu[i] = k;
      u[i] = c;
   }
/* ----- component loads ----- */
   fgets(dummy,80,fptr1);
   for (i = 0; i < nl; i++) {
      fscanf(fptr1, "%d %lf\n", &k, &c);
      f[k-1] = c;
   }
/* ----- material properties ----- */
   fgets(dummy,80,fptr1);
   for (i = 0; i < nm; i++){
      fscanf(fptr1, "%d", &k);
      for (j = 0; j < npr; j++) {
         fscanf(fptr1, "%lf\n", &c);
	 pm[(k-1)*npr+j] = c;
      }
   }
/* ----- multipoint constraints ----- */
   if (nmpc > 0)
      { fgets(dummy,80,fptr1);
	for(j=0;j<nmpc;j++){
	   fscanf(fptr1,"%lf",&c);
	   beta[3*j]=c;
	   fscanf(fptr1,"%d",&k);
	   mpc[2*j]=k;
           fscanf(fptr1,"%lf",&c);
	   beta[3*j+1]=c;
           fscanf(fptr1,"%d",&k);
           mpc[2*j+1]=k;
           fscanf(fptr1,"%lf",&c);
	   beta[3*j+2]=c;
	   }
	}
   fclose (fptr1);
   prefront(nn,ne,nen,ndn,nq,nmpc,noc,mpc,&ibl);
/* ------------------------------------------------- */
   isbl = (int *) calloc(ibl+1, sizeof(int));
   indx = (int *) calloc(ibl+1, sizeof(int));
   s = (double *) calloc(ibl*ibl, sizeof(double));
/* ------------------------------------------------- */
     nfron = 0;
	 ntogo = 0;
	 ndcnt = 0;
     for (i = 1; i <= ibl; i++) {
	    indx[i] = i;
	    }
     icount = 0;
/* =====  frontal assembly & eliminaton etc.  ===== */
/* ----- corner nodes and integration points */
     integ(xi,xni);
/* ----- open scratch file for writing ----- */
     if ((fptr = fopen("scratch.dat","wb"))==NULL)
       {printf("Can't open file scratch.dat"); exit(1);}
/* ----- element loop ----- */
     mtn1 = 0;
     for (n = 0; n < ne; n++) {
        printf("... forming stiffness matrix of element %d\n", n+1);
        mtn = mat[n];
	    if (mtn != mtn1)
	       dmat(mtn,&al,pm,d);
	       mtn1 = mtn;
	elstif(n,se,qt,xi,xni,d,b,db,tempr,x,al,noc);
        if (n == 0) {
           cnst = 0;
	   for (i = 0; i < 24; i++) {
	           cnst = cnst + se[i][i];
	           }
           cnst = 1e+11 * cnst;
           mpcfron(indx,isbl,mpc,nmpc,&nfron,s,f,ibl,beta,cnst);
	  }
     /* ----- add temperature load to the force vector f[] ----- */
       for (i = 0; i < 8; i++) {
	  ii = 3*(abs(noc[nen*n])-1);
	  for (j = 0; j < 3; j++) {
	     f[ii+j] = f[ii+j] + qt[3*i+j];
	     }
	  }
 /* frontal assembly  and forward elimination */
front(n,noc,nen,ndn,nd,&icount,indx,isbl,ibl,s,f,&nfron,&ntogo,&ndcnt,se,nu,cnst,u);
      }
fclose(fptr);
/* ----- assembly and reduction are complete */
/* ----- now backsubstitute */
     if ((fptr = fopen("scratch.dat","rb"))==NULL)
       {printf("Can't open file scratch.dat"); exit(1);}
     backsub(icount,f);
     fclose(fptr);
/* ----- printing displacements ----- */
   fptr1 = fopen(file2, "w");
   printf("\n%s\n", title);
   fprintf(fptr1, "\n%s\n", title);
 fprintf(fptr1, "node#     x-displ      y-displ      z-displ\n");
 printf ("node#     x-displ      y-displ      z-displ\n");
for (i = 0; i < nn; i++) {
 printf(" %4d  %11.4e  %11.4e  %11.4e\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);
 fprintf(fptr1," %4d  %11.4e  %11.4e  %11.4e\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);
}
/* ----- reaction calculation ----- */
   printf("node#    reaction\n");
   fprintf(fptr1, "node#     reaction\n");
   for (i = 0; i < nd; i++) {
      k = nu[i];
      reaction = cnst * (u[i] - f[k-1]);
      printf(" %4d  %11.4e\n", k, reaction);
      fprintf(fptr1, " %4d  %11.4e\n", k, reaction);
   }
/* -----  stress calculations ----- */
     mtn1 = 0;
     for (n = 0; n < ne; n++) {
     fprintf(fptr1, "von mises stress at 8 int. pts. in elem# %d\n", n+1);
        mtn = mat[n];
	  if (mtn != mtn1)
	     dmat(mtn,&al,pm,d);
	     mtn1 = mtn;
      cal = al * tempr[n];
        for (ip = 0; ip < 8; ip++) {
           /* --- von mises stress at integration points */
	   dbmat(n,ip,x,noc,d,b,db,&dj,xi,xni);
      /* --- element nodal displacements stored in qt() */
           for (i = 0; i < 8; i++) {
              in = 3 * (abs(noc[nen*n+i]) - 1);
	      ii = 3 * i;
              for (j = 0; j < 3; j++) {
                 qt[ii + j] = f[in + j];
                }
             }
           /* --- stress calculation str = db * q */
           for (i = 0; i < 6; i++) {
              str[i] = 0;
              for (j = 0; j < 24; j++) {
                 str[i] = str[i] + db[i][j] * qt[j];
                 }
              str[i] = str[i] - cal * (d[i][1] + d[i][2] + d[i][3]);
              }
           /* --- calculation of von mises stress at ip */
           siv1 = str[0] + str[1] + str[2];
           siv2 = str[0] * str[1] + str[1] * str[2] + str[2] * str[0];
           siv2 = siv2 - str[3]*str[3] - str[4]*str[4] - str[5]*str[5];
           vm = sqrt((double) siv1 * siv1 - 3 * siv2);
           if (ip == 4)
              fprintf(fptr1,"\n");
           fprintf(fptr1, "   %11.4e", vm);
          }
        fprintf(fptr1,"\n");
       }
     printf("the results are saved in the file %s\n", file2);
     fclose(fptr1);
     return(0);
   }
/* end of main function */

integ(xi,xni)
double xi[][8], xni[][8];
{
  int i;
  double c;
/* ------- integration points xni() -------- */
     c = .57735026919;
     xi[0][0] = -1;
     xi[1][0] = -1;
     xi[2][0] = -1;
     xi[0][1] = 1;
     xi[1][1] = -1;
     xi[2][1] = -1;
     xi[0][2] = 1;
     xi[1][2] = 1;
     xi[2][2] = -1;
     xi[0][3] = -1;
     xi[1][3] = 1;
     xi[2][3] = -1;
     xi[0][4] = -1;
     xi[1][4] = -1;
     xi[2][4] = 1;
     xi[0][5] = 1;
     xi[1][5] = -1;
     xi[2][5] = 1;
     xi[0][6] = 1;
     xi[1][6] = 1;
     xi[2][6] = 1;
     xi[0][7] = -1;
     xi[1][7] = 1;
     xi[2][7] = 1;
     for (i = 0; i < 8; i++) {
        xni[0][i] = c * xi[0][i];
        xni[1][i] = c * xi[1][i];
        xni[2][i] = c * xi[2][i];
       }
     return(0);
 }
dmat(mtn,al,pm,d)
   int mtn;
   double *al,*pm,d[][6];
  {
   double e,pnu,c1,c2;
   int i,j;
   /* --- d() matrix relating stresses to strains */
     e = pm[3*(mtn-1)];
	 pnu = pm[3*(mtn-1)+1];
	 *al = pm[3*(mtn-1)+2];
     c1 = e / ((1 + pnu) * (1 - 2 * pnu));
     c2 = .5 * e / (1 + pnu);
     for (i = 0; i < 6; i++) {
	     for (j = 0; j < 6; j++) {
	        d[i][j] = 0;
	       }
	    }
     d[0][0] = c1 * (1 - pnu);
     d[0][1] = c1 * pnu;
     d[0][2] = d[0][1];
     d[1][0] = d[0][1];
     d[1][1] = d[0][0];
     d[1][2] = d[0][1];
     d[2][0] = d[0][2];
     d[2][1] = d[1][2];
     d[2][2] = d[0][0];
     d[3][3] = c2;
     d[4][4] = c2;
     d[5][5] = c2;
     return(0);
 }
elstif(n,se,qt,xi,xni,d,b,db,tempr,x,al,noc)
   int n,*noc;
   double al,*x,*tempr,d[][6],qt[],xi[][8],xni[][8];
   double se[][24],b[][24],db[][24];
{
   int i,j,k,ip;
   double dj,dte,c,dsum;
/* --------  element stiffness and temperature load  ----- */
     for (i = 0; i < 24; i++) {
	    for (j = 0; j < 24; j++) {
           se[i][j] = 0;
	       }
	     qt[i] = 0;
	     }
     dte = tempr[n];
     /* --- weight factor is one */
     /* --- loop on integration points */
     for (ip = 0; ip < 8; ip++) {
        printf("integration point = %d\n", ip);
        dbmat(n,ip,x,noc,d,b,db,&dj,xi,xni);
        /* --- element stiffness matrix  se */
        for (i = 0; i < 24; i++) {
           for (j = 0; j < 24; j++) {
              for (k = 0; k < 6; k++) {
                 se[i][j] = se[i][j] + b[k][i] * db[k][j] * dj;
                }
             }
          }
        /* --- determine temperature load qt() */
        c = al * dte;
        for (i = 0; i < 24; i++) {
	   dsum = db[0][i] + db[1][i] + db[2][i];
	   qt[i] = qt[i] + c * fabs(dj) * dsum / 6;
           }
     }
     return(0);
 }
dbmat(n,ip,x,noc,d,b,db,dj,xi,xni)
  int n,*noc,ip;
  double d[][6],b[][24],db[][24],xi[][8],xni[][8],*x,*dj;
{
 int i,j,k,kn,ir,ic;
 double aj[3][3],tj[3][3],gn[3][8],h[9][24],g[6][9];
 double dj1,dj2,dj3,c;
/* -------  db()  matrix  ------ */
/* --- gradient of shape functions - the gn() matrix */
     for (i = 0; i < 3; i++) {
        for (j = 0; j < 8; j++) {
           c = 1;
           for (k = 0; k < 3; k++) {
              if (k != i)
                 c = c * (1 + xi[k][j] * xni[k][ip]);