CST.rar_cst_finite element c++

/***************************************/ /* program cst */ /* 2-d constant strain triangle */ /***************************************/ #include <stdio.h> #include <math.h> main() { FILE *fptr1, *fptr2; int n,i,j,k,m,i1,i2,i3,ii,jj,m1,nmin,nmax,nrt,nct,it,jt; int nr,nc; char dummy[81], title[81], file1[81], file2[81], file3[81]; int ne,nn,nq,nm,nd,nl,nen,ndn,ndim,npr,nbw,nmpc,lc,ipl; int *noc, *nu, *mat, *mpc; float *x, *thick, *pm, *u, *tempr, *s, *f, *beta; float c, dj, al, pnu, tld, cnst, reaction, s1, s2, s3, ang, r; float b[3][6], d[3][3], db[3][6], se[6][6], str[3], tl[6]; /*-------------------------------------------------------*/ printf("\n"); puts("Input file name < dr:fn.ext >: "); gets(file1); puts("Output file name < dr:fn.ext >: "); gets(file2); printf("\n"); printf(" 1) plane stress analysis\n"); printf(" 2) plane strain analysis\n"); printf(" choose <1 or 2> "); scanf("%d", &lc); if (lc < 1 || lc > 2) lc = 1; 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); fgets(dummy, 80, fptr1); fscanf(fptr1,"%d %d %d %d %d\n", &nd, &nl, &nmpc); npr = 3; /* Material Properties */ /* ----- memory allocation ----- */ x = (float *) calloc(nn*ndim, sizeof(float)); noc = (int *) calloc(ne*nen, sizeof(int)); u = (float *) calloc(nd, sizeof(float)); nu = (int *) calloc(nd, sizeof(int)); mat = (int *) calloc(ne,sizeof(int)); thick = (float *) calloc(ne, sizeof(float)); f = (float *) calloc(nn*ndn, sizeof(float)); tempr = (float *) calloc(ne, sizeof(float)); pm = (float *) calloc(nm*npr, sizeof(float)); mpc = (int *) calloc(2*nmpc, sizeof(int)); beta = (float *) calloc(3*nmpc, sizeof(float)); printf("\n\n PLOT CHOICE\n"); printf(" 1) no plot data\n"); printf(" 2) create data file for in-plane shear stress\n"); printf(" 3) create data file for von mises stress\n"); printf(" choose <1 or 2 or 3> "); scanf("%d%*c", &ipl); if(ipl < 1 || ipl > 3) ipl = 1; /* --- default is no data ---*/ if(ipl > 1){ printf("Output file name < dr:fn.ext >:\n"); gets(file3); } /* ----- 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, "%f\n", &c); x[ndim*(n-1)+j] = c; } } /* ----- connectivity, material, thickness, 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,"%f",&c); thick[n-1] = c; fscanf(fptr1,"%f\n",&c); tempr[n-1] = c; } /* ----- displacement bc ----- */ fgets(dummy,80,fptr1); for (i = 0; i < nd; i++) { fscanf(fptr1, "%d %f\n", &k, &c); nu[i] = k; u[i] = c; } /* ----- component loads ----- */ fgets(dummy,80,fptr1); for (i = 0; i < nl; i++) { fscanf(fptr1, "%d %f\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, "%f\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,"%f",&c); beta[3*j]=c; fscanf(fptr1,"%d",&k); mpc[2*j]=k; fscanf(fptr1,"%f",&c); beta[3*j+1]=c; fscanf(fptr1,"%d",&k); mpc[2*j+1]=k; fscanf(fptr1,"%f",&c); beta[3*j+2]=c; } } fclose (fptr1); /* ----- bandwidth nbw from connectivity noc() and mpc ----- */ nbw = 0; for (i = 0; i < ne; i++) { nmin = noc[nen*i]; nmax = nmin; for (j = 1; j < 3;j++) { n =noc[nen*i+j]; if (nmin > n) nmin = n; if (nmax < n) nmax = n; } n= ndn * (nmax - nmin + 1); if (nbw < n) nbw = n; } for (i = 0; i < nmpc; i++) { n = abs(mpc[2*i] - mpc[2*i+1]) + 1; if (nbw < n) nbw = n; } printf ("the bandwidth is %d\n", nbw); /* ----- allocate memory for stiffness ----- */ s = (float *) calloc(nq*nbw, sizeof(float)); /* ----- global stiffness matrix -----*/ for (n = 0; n < ne; n++) { printf("forming stiffness matrix of element %d\n", n+1); dbmat(&lc,n,b,d,db,mat,pm,npr,x,noc,&dj,&al,&pnu,&i1,&i2,&i3); /* --- element stiffness --- */ for (i = 0; i < 6; i++) { for (j = 0; j < 6; j++) { c = 0.; for (k = 0; k < 3; k++) { c = c + .5 * fabs(dj) * b[k][i] * db[k][j] * thick[n]; } se[i][j] = c; } } /* --- temperature load vector --- */ c = al * tempr[n]; if (lc == 2) c = c * (1 + pnu); for (i = 0; i < 6; i++) { tl[i] = .5 * c * thick[n] * fabs(dj) * (db[0][i] + db[1][i]); } printf (".... placing in global locations\n"); for (ii = 0; ii < nen; ii++) { nrt = ndn * (noc[nen*n + ii] - 1); for (it = 0; it < ndn; it++) { nr = nrt + it; i = ndn * ii + it; for (jj = 0; jj < nen; jj++) { nct = ndn * (noc[nen*n+jj] - 1); for (jt = 0; jt < ndn; jt++) { j = ndn * jj + jt; nc = nct + jt - nr; if (nc >= 0) s[nbw*nr+nc] = s[nbw*nr+nc] + se[i][j]; } } f[nr] = f[nr] + tl[i]; } } } /* ----- decide penalty parameter cnst ----- */ cnst = 0.; for (i = 0; i < nq; i++) { if (cnst < s[i*nbw]) cnst = s[i*nbw]; } cnst = cnst * 10000.; /* ----- modify for displacement boundary conditions ----- */ for (i = 0; i < nd; i++) { k = nu[i]; s[(k-1)*nbw] = s[(k-1)*nbw] + cnst; f[k-1] = f[k-1] + cnst * u[i]; } /* ----- modify for multipoint constraints ----- */ for (i = 0; i < nmpc; i++){ i1 = mpc[2*i]-1; i2 = mpc[2*i+1]-1; s[i1*nbw] = s[i1*nbw] + cnst*beta[3*i]*beta[3*i]; s[i2*nbw] = s[i2*nbw] + cnst*beta[3*i+1]*beta[3*i+1]; n=i1; if (n > i2) n = i2; m = abs(i2-i1); s[n*nbw+m] = s[n*nbw+m]+cnst*beta[3*i]*beta[3*i+1]; f[i1] = f[i1] + cnst*beta[3*i]*beta[3*i+2]; f[i2] = f[i2] + cnst*beta[3*i+1]*beta[3*i+2]; } /* ----- solution of equations using band solver ----- */ bansol(s,f,nq,nbw); /* ----- printing displacements ----- */ fptr1 = fopen(file2, "w"); printf("\n%s\n", title); fprintf(fptr1, "\n%s\n", title); fprintf(fptr1, "bandwidth = %d\n",nbw); if (lc == 1) fprintf(fptr1, "plane stress analysis\n"); if (lc == 2) fprintf(fptr1, "plane strain analysis\n"); fprintf(fptr1, "node# x-displ y-displ\n"); printf ("node# x-displ y-displ\n"); for (i = 0; i < nn; i++) { printf(" %4d %11.4e %11.4e\n",i+1,f[2*i],f[2*i+1]); fprintf(fptr1," %4d %11.4e %11.4e\n",i+1,f[2*i],f[2*i+1]); } /* ----- 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, reactio