/***************************************/
/* 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
