#include <stdio.

h>
#include <math.h>
#include<stdlib.h>
#define ABS(x) (x < 0 ? -(x) : (x))
#define EPS 0.00001
/* Begin-Function definitions */
integ(float xni[][2]);
dmat(int,float*,int*,int,float*,float*,int,float d[][3]);
dbmat(int,float*,int*,float*,float*,float d[][3],float b[][8],float db[][8],floa
t*,float,float);
elemstif(int,int,float se[][8],float tl[8],float xni[][2],float d[][3],float*,fl
oat*,float*,float,float,int*);
dbmat2(int,float*,int*,float*,float*,float d[][3],float b[][8],float db[][8],flo
at*,float,float);
elemstif2(int,int,float se[][8],float tl[8],float xni[][2],float d[][3],float*,f
loat*,float*,float,float,int*);
GSolve(float **,int ,float *);
/* End-Function definitions */

/* Begin-Main Code */
main()
{
// Begin-Declaring variables in main code
FILE *fptr; //pointer to a file
int n,i,j,k,ii,jj; //dummy integers for storing values temporarily and f
or using in loops
int nmin,nmax,nrt1,nct1,nrt,nct,it,jt;
int ip,nr1,nc1,nr,nc,in;
char dummy[81], title[81], file1[81], file3[81],file4[81],file5[81];
int ne,nn,nq,nm,nd,nl,nen,ndn,ndim,npr,nbw,opt,nbw2,nq1,nq2,nbw3,nbw4;
// Numbers of---> ne: elements; nn: nodes; nq: total d.o.f; nm: material
s
// nd: disp. BCs; nl: load BCs; nen:el. nodal conn.=4; ndn: dof per node
=2
// ndim: 2-dim; npr: mat. properties=3 (E,nu,alpha); nbw, nbw3: band wid
th of specimen
// nbw2, nbw4: band width for punch; opt: option for plane stress/strain
// nq1,nq2: dof for specimen and punch
int ne1=857,ne2=510,nn1=896,nn2=544;
// ne1: no. of el. in spec.; ne2: no. of el. in punch
 // nn1: no. of nodes in spec.; nn2: no. of nodes in punch
int *noc, *nu, *mat, *noc1, *noc2; //1-d arrays
// noc: element connectivity for full mesh; nu: dof numbering for disp B
Cs
// mat: material array; noc1, noc2: elem. conn. for spec and punch
float *x, *thick, *pm, *u, *tempr, *s1,*s2, *f,*x1,*xx1,*xx2, *g;
// x: x and y coords stored alternately; pm: material properties stored
// u: disp BC values; s1, s2: banded stiffness matrices for spec and pun
ch
//f: load bcs; x1: disps after loading; xx1, xx2: for spec and punch
// g is x,y coordinates after deformation
float **a,**k1,**k2,**kk;
// a:global stiffness matrix including forces after application of BCS (
nq+1)X(nq)
// k1, k2: global stiff. mats for spec. and punch
// kk: full stiffness matrix including spring which goes into a
float xi,eta,th,q[8],c1,sv;
// xi, eta: natural coordinates; th: thickness; q: dof per element
// c1: ; sv:
float c, dj, al, pnu, cnst, reaction;
// dj: det. of Jacobian, al: temp. coeff; pnu: poisson's ratio
// cnst: constant for penalty method used in stiffness matrix
// reaction: Reactions
float b[3][8],d[3][3],db[3][8],se[8][8],tl[8],xni[4][2];
// se: elemental stiffness matrix; tl: temp. load; xni: numerical interg
ration mat.
int *mat1;
float *thick1, *tempr1;
int inc; // increment
float *np; // array used for storing intermediate coordinates during ite
rations
int spc[20]={40,58,57,56,55,54,53,52,51,50,49,48,47,46,45,44,43,42,41,21
}; // Contacting nodes in spec
int puc[20]={ 1080 , 1099, 1100,1101,1102,1103,1104,1105,1106,1107,1108,
1109,1110,1111,1112,1113,1114,1115,1116,928}; // Contacting nodes in pad
int conv; // to verify for convergence
float reldisp[11]; // relative displacement for the contacting nodes
//end-declaring variables in main code
/*-------------------------------------------------------*/

/* User inputs */
printf("\n");
puts("Input file name < dr:fn.ext >: ");
gets(file1);
puts("Old coordinate file name < dr:fn.ext >: ");
gets(file3);
puts("new coordinate file name <dc:fn.ext>: ");
gets(file4);
 printf("output file name<dr:fn.ext>:");
gets(file5);
printf("\n");
printf(" 1 plane stress analysis\n");
printf(" 2 plane strain analysis\n");
printf(" enter 1 or 2 ");
scanf("%d", &opt);
if (opt < 1 || opt > 2)
opt = 2;
// Begin-reading data from input file
fptr = fopen(file1, "r");
fgets(dummy,80,fptr); // stores first line in a dummy variable
fgets(title,80,fptr); // stores second line in a dummy variable
fgets(dummy,80,fptr); // stores third line in a dummy variable
fscanf(fptr,"%d %d %d %d %d %d\n", &nn, &ne, &nm, &ndim, &nen, &ndn);
// scanning values from fourth line
fgets(dummy, 80, fptr); // stores fifth line in a dummy variable
fscanf(fptr,"%d %d\n", &nd, &nl); // scanning values from 6th line
npr = 3; /* Material properties E, Nu, Alpha */

/* ----- begin-memory allocation ----- */

x = (float *) calloc(nn*ndim, sizeof(float));
g =(float *)calloc(nn*ndim, sizeof(float));
xx1 = (float *)calloc(nn1*ndim, sizeof(float));
xx2 = (float *)calloc(nn2*ndim, sizeof(float));
noc = (int *) calloc(ne*nen, sizeof(int));
noc1 = (int *) calloc(ne1*nen, sizeof(int));
noc2 = (int *) calloc(ne2*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));
mat1=(int *) calloc(ne, sizeof(int));
thick1=(float *) calloc(ne, sizeof(float));
tempr1=(float *) calloc(ne, sizeof(float));

/* ----- total dof is nq ----- */

nq = ndn * nn;
nq1 = ndn*nn1;
nq2 = ndn*nn2;
x1 = (float *) calloc(nq, sizeof(float));
/* ----- end-memory allocation ----- */
/* =============== read data again ==================== */

/* ----- connectivity, material, thickness, temp-change ----- */

fgets(dummy,80,fptr); // stores eighth line in dummy variable
printf("%s I am at beginning of nodes \n",dummy);

/* begin-loop to read element connectivity data line by line */

for (i = 0; i < ne; i++)
{
fscanf(fptr,"%d", &n); // Element number stored in n
for (j = 0; j < nen; j++)
{
fscanf(fptr,"%d", &k); // nen=4 Nodes corresponding to nth eleme
nt
noc[(n-1)*nen+j]=k; // list of nodes attached to all elements in
groups of 4
}

fscanf(fptr,"%d", &k);
mat[n-1] = k; // storing material number
fscanf(fptr,"%f",&c);
thick[n-1] = c; // storing thickness
fscanf(fptr,"%f\n",&c);
tempr[n-1] = c; // storing temperature
 }
/* end-loop to read element connectivity data line by line */

for (i=0; i<ne1*nen; i++)

{ // storing nodal connectivity for spec
noc1[i]=noc[i];
}

for (i=0;i<ne2*nen;i++)
{// storing nodal connectivity for punch
noc2[i]=noc[ne1*nen+i];
}

for(i=0;i<ne2;i++)
{// separating material definition for punch
mat1[i]=mat[857+i];
thick1[i]=thick[857+i];
tempr1[i]=tempr[857+i];
}

/* ----- displacement bc ----- */

fgets(dummy,80,fptr); // reading the dummy line pointing to dof disp bcs
for (i = 0; i < nd; i++)
{
fscanf(fptr, "%d %f\n", &k, &c);
nu[i] = k; // dof number for disp BC
u[i] = c; // value of disp BC
}

/* ----- component loads ----- */

fgets(dummy,80,fptr); // reading the dummy line pointing to dof load bcs
for (i=0; i< nn*ndn; i++)
{
f[i]=0;
}
for (i = 0; i < nl; i++)
{
 fscanf(fptr, "%d %f\n", &k, &c);
f[k-1] = c; // value of force, c, for the corresponding dof, k, replacing
the initialized zero value
}
/* ----- material properties ----- */
fgets(dummy,80,fptr);
for (i = 0; i < nm; i++)
{
fscanf(fptr, "%d", &k);
for (j = 0; j < npr; j++)
{
fscanf(fptr, "%f\n", &c);
pm[(k-1)*npr+j] = c;
}
}

fclose (fptr);
// End-reading data from input file
/* ----- coordinates ----- */
printf("\nenter the step number or increment number\n");
scanf("%d",&inc);
//begin-reading node coordinates
if(inc==1)
{// Reads from old co-ordinates file
fptr=fopen(file3,"r");
fgets(dummy,80,fptr); // Header saying node no, x-coord and y-coord
for (i = 0; i < nn; i++)
{
fscanf(fptr, "%d", &n); // Node number
for (j = 0; j < ndim; j++)
{
fscanf(fptr, "%f\n", &c);
x[ndim*(n-1)+j] = c; // coordinates stored alter
nately
}
 }

for (i=0;i<nq1;i++)
{// coordinates for spec
xx1[i]=x[i];
}

for (i=0;i<nq2;i++)
{//coords for punch
xx2[i]=x[nq1+i];
}
/* for(i=0;i<nq2;i++)
{
printf("%f\n",xx2[i]);
}*/
fclose(fptr);
}
else
{// Reads from new co-ordinates file, this has no header or node number
fptr=fopen(file4,"r");
for(i=0;i<nn;i++)
{
for(j=0;j<ndim;j++)
{
fscanf(fptr,"%f\n",&c);
g[ndim*i+j]=c;
}
}
for(i=0;i<nq1;i++)
{
xx1[i]=g[i];
}
for(i=0;i<nq2;i++)
{
xx2[i]=g[nq1+i];
}
fclose(fptr);
}
//end-reading node coordinates
/* ----- bandwidth nbw for SPECIMEN from connectivity noc() ----- */
nbw = 0;
for (i = 0; i < ne1; i++)
{// Loop to go element by element
 nmin = noc1[nen*i]; // initializing min node no. to first node
nmax = nmin; // initializing max node no. to first node
for (j = 1; j < 4;j++)
{
n =noc1[nen*i+j];
if (nmin > n)
nmin = n;
if (nmax < n)
nmax = n;
}
n= ndn * (nmax - nmin + 1);
if (nbw < n)
nbw = n;
}
printf ("the bandwidth of specimen is %d\n", nbw);

/* ----- bandwidth nbw3 for PUNCH from connectivity noc() ----- */

nbw3 = 0;
for (i = 0; i < ne2; i++)
{// Loop to go element by element
nmin = noc2[nen*i];
nmax = nmin;
for (j = 1; j < 4;j++)
{
n =noc2[nen*i+j];
if (nmin > n)
nmin = n;
if (nmax < n)
nmax = n;
}
n= ndn * (nmax - nmin + 1);
if (nbw3 < n)
nbw3 = n;
 }
printf ("the bandwidth for punch is %d\n", nbw3);

/* ----- allocate memory for banded stiffness array----- */

s1 = (float *) calloc(nq1*nbw, sizeof(float));
s2 = (float *) calloc(nq2*nbw3, sizeof(float));
for (i=0; i< nq1*nbw; i++)
{
s1[i]=0;
}
for (i=0; i< nq2*nbw3; i++)
{
s2[i]=0;
}

/* ----- global stiffness matrix, corner nodes and integrationpoints ----- */

integ(xni);
for (n = 0; n < ne1; n++)
{//begin- n loop going element by element in spec
printf("forming stiffness matrix of specimen %d\n", n+1);
dmat(n,pm,mat,npr,&pnu,&al,opt,d); // Constitutive matrix

/* --- element stiffness --- */

elemstif(n,opt,se,tl,xni,d,thick,tempr,xx1,al,pnu,noc1);
printf (".... placing in global locations\n");
for (ii = 0; ii < nen; ii++)
{//begin - ii loop to go through four nodes corresonding to each element
nrt1 = ndn * (noc1[nen*n + ii] - 1);
for (it = 0; it < ndn; it++)
{//begin - it loop to go through dimensions
nr1 = nrt1 + it; // gives position of Ux and Uy (row number in ban
ded matrix)
i = ndn * ii + it; // goes from 0 to 7, corresponding to local sti
ffness matrix
for (jj = 0; jj < nen; jj++)
{// begin - jj loop similar to ii
 nct1 = ndn * (noc1[nen*n+jj] - 1);
for (jt = 0; jt < ndn; jt++)
{//begin - jt loop similar to it
nc1 = nct1 + jt - nr1; // corresponds to column number in ba
nded matrix
j = ndn * jj + jt; // goes from 0 to 7, corresponding to loc
al stiffness matrix
if (nc1 >= 0)
s1[nbw*nr1+nc1] = s1[nbw*nr1+nc1] + se[i][j];
}//end - jt loop
}//end - jj loop
}//end - it loop
}//end - ii loop
}//end- n loop going element by element in spec

for (n = 0;n < ne2;n++)

{
printf("forming stiffness matrix of punch %d\n", n+1);
dmat(n,pm,mat1,npr,&pnu,&al,opt,d);

/* --- element stiffness --- */

elemstif2(n,opt,se,tl,xni,d,thick1,tempr1,xx2,al,pnu,noc2);
printf (".... placing in global locations\n");
for (ii = 0; ii < nen; ii++)
{
nrt = ndn*(noc2[nen*n + ii] - nn1-1);
for (it = 0; it < ndn; it++)
{
nr = nrt + it;
i = ndn * ii + it;
for (jj = 0; jj < nen; jj++)
{
nct = ndn * (noc2[nen*n+jj] - nn1-1);
for (jt = 0; jt < ndn; jt++)
{
j = ndn * jj + jt;
 nc = nct + jt - nr;
if (nc >= 0)
s2[nbw3*nr+nc] = s2[nbw3*nr+nc] + se[i][j];
}
}
}
}
}

printf("\nNOW ASSEMBLING,CACULATING AND THEN DISPLAYING THE RESULTS,PLEASE BE PA

TIENT\n");
/*allocating memory for SPECIMEN stiffness matrix */
k1 = (float **)malloc((nq1)*sizeof(float *));
for (i=0;i<nq1;i++)
k1[i] = (float *)malloc(nq1*sizeof(float));

/*------------converting banded stiffness matrix to general symmetric matrix---

-------*/

/* begin - initializing matrix to 0 */

for(i=0;i<nq1;i++)
{
for(j=0;j<nq1;j++)
{
k1[i][j]=0;
}
}
/* end - initializing matrix to 0 */

nbw2=nbw; // counter which decreases with every step

for(i=0;i<nq1;i++)
{
if(i>=nq1-nbw+1)
 {
nbw2--; // Last few rows have decreasing band width !!
}
for(j=i;j<i+nbw2;j++)
{
k1[i][j]=s1[i*nbw+j-i];
}
}

for(i=1;i<nq1;i++)
{
for(j=0;j<i;j++)
{
k1[i][j]=k1[j][i];
}
}

/*allocating memory for PUNCH stiffness matrix */

k2 = (float **)malloc((nq2)*sizeof(float *));
for (i=0;i<nq2;i++)
k2[i] = (float *)malloc(nq2*sizeof(float));

/*------------converting banded stiffness matrix to general symmetric matrix---

-------*/
for(i=0;i<nq2;i++)
{
for(j=0;j<nq2;j++)
{
k2[i][j]=0;
}
}
nbw4=nbw3;
for(i=0;i<nq2;i++)
{
if(i>=nq2-nbw3+1)
{
nbw4--;
}
for(j=i;j<i+nbw4;j++)
{
k2[i][j]=s2[i*nbw3+j-i];
}
}

for(i=1;i<nq2;i++)
{
for(j=0;j<i;j++)
{
k2[i][j]=k2[j][i];
}
}

/*------------------------------SPRING COUPLING ASSEMBLY------------------------

-----*/

/*-----------stiffness matrix accounting for springs initialised to zeros ksp

r[3110][3110]---------*/
/*-----------combined final stiffness matrix initialised to zeros kk[3110][311
0] ---------*/
kk = (float **)malloc((nq)*sizeof(float *));
for (i=0;i<nq;i++)
 kk[i] = (float *)malloc(nq*sizeof(float));
for(i=0;i<nq;i++)
{
for(j=0;j<nq;j++)
{
kk[i][j]=0;
}
}
/*-----------spring matrix accounting for 1 node set contact 2 & 804---------*/
for(i=1;i<inc+1;i++)
{
kk[2*spc[i-1]-1][2*spc[i-1]-1]=1e12;
kk[2*spc[i-1]-1][2*puc[i-1]-1]=-1e12;
kk[2*puc[i-1]-1][2*spc[i-1]-1]=-1e12;
kk[2*puc[i-1]-1][2*puc[i-1]-1]=1e12;
}
/*just like 1st node contact*/
/*kk[3][3]=1e12;
kk[3][1607]= -1e12;
kk[1607][3]= -1e12;
kk[1607][1607]=1e12;
/*-----------adding k1 ,k2,& kspr to give kk(final)---------*/

for(i=0;i<1792;i++)
{
for(j=0;j<1792;j++)
{
kk[i][j]=kk[i][j]+k1[i][j];
}
}

for(i=1792;i<2880;i++)
{
for(j=1792;j<2880;j++)
{
 kk[i][j]=kk[i][j]+k2[i-1792][j-1792];
}
}

/* ----- decide penalty parameter cnst ----- */

cnst = 0.;
for (i = 0; i < nq; i++) {
if (cnst < kk[i][i])
cnst = kk[i][i];
}
cnst = cnst * 10000.;
/* ----- modify for displacement boundary conditions ----- */
for (i = 0; i < nd; i++) {
k = nu[i];
kk[k-1][k-1] = kk[k-1][k-1]+ cnst;
f[k-1] = f[k-1] + cnst * u[i];
}
/*------------------------------------------------------------------------------
---------*/

a = (float **)malloc((nq+1)*sizeof(float *));

for (i=0;i<nq+1;i++)
a[i] = (float *)malloc(nq*sizeof(float));
for(i=0;i<nq;i++)
{
for(j=0;j<nq;j++)
{
a[i][j]=kk[i][j];
}
}
for(i=0;i<nq;i++)
{
a[nq][i]=f[i];
}

GSolve(a,nq,x1);
/* ----- solution of equations using gauss solver ----- ------------------------
----------------------------------------*/
/* bandsolv(s,f,nq,nbw); */

/* ----- printing displacements ----- */

fptr = fopen(file5, "w");//opening of a file in which disp will be print
ed
printf("\n%s\n", title);
fprintf(fptr, "\n%s\n", title);
fprintf(fptr,"bandwidth for spcimen is %d\n",nbw);
fprintf(fptr,"bandwidth for punch is %d\n",nbw3);
if (opt == 1)
fprintf(fptr, "plane stress analysis\n");
if (opt == 2)
fprintf(fptr, "plane strain analysis\n");
fprintf(fptr, "nodulal 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,x1[2*i],x1[2*i+1]);
fprintf(fptr," %4d %11.4e %11.4e\n",i+1,x1[2*i],x1[2*i+1]);
}
fclose(fptr);//closing of this file
/*NOW FINDING OVERALL DISP OF EACH NODE*/
/*allocating memory for new displacements*/
np=(float *) calloc(nn*ndn, sizeof(float));
if(inc==1)
{
printf("\n now writing final coordinates for increment1 in user choice file 4
");
for(i=0;i<nn;i++)
{
 g[2*i]=x[2*i]+x1[2*i]; //the final x-coordinate after deformation wit
h increment1
g[2*i+1]=x[2*i+1]+x1[2*i+1];//the final y-coordinat after deformation wi
th increment1
}
fptr=fopen(file4,"w");
for(i=0;i<nn;i++)
{
fprintf(fptr,"%.15f\t %.15f\n",g[2*i],g[2*i+1]);
}
fclose(fptr);
}
else
{
printf("\n now writing final coordinates for increment no-%d in user choice fi
le 4",inc);
for(i=0;i<nn;i++)
{
np[2*i]=g[2*i]+x1[2*i];//the final x-coordinate after deformation with in
crement more than one
np[2*i+1]=g[2*i+1]+x1[2*i+1];//the final y-coordinate after deformation w
ith increment more than one
}
fptr=fopen(file4,"w");
for(i=0;i<nn;i++)
{
fprintf(fptr,"%.15f\t %.15f\n",np[2*i],np[2*i+1]);
}
fclose(fptr);
}
/*now calculating coordinate difference*/
if(inc==1)
{
for(i=0;i<20;i++)
{
reldisp[i]=g[2*spc[i]-1]-g[2*puc[i]-1];
}
}
else
{
for(i=0;i<20;i++)
{
reldisp[i]=np[2*spc[i]-1]-np[2*puc[i]-1];
}
}
//finally printing cooordinate difference on screen
for(i=0;i<20;i++)
{
printf("\n%0.15f",reldisp[i]);
}
printf("\n IS THERE ANY CONVERGENCE OBTAINED WITH THIS STEP,IF YES THEN PRESS 1,
IF NO THEN PRESS 2\n");
scanf("%d",&conv);
if(conv==1)
{
fptr=fopen("reldisp.txt","w");
fprintf(fptr,"this is for step number %d\n",inc);
for(i=0;i<20;i++)
{
fprintf(fptr,"%.15f\n",reldisp[i]);
}
fclose(fptr);
}
printf("\nthanks for ur patience..for complete results,check the output files
in this folder");
return(0);
}

integ(xni)
float xni[][2];
{
float c;
/* ----- integration points xni() ----- */
c = 0.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);
}

dmat(n,pm,mat,npr,pnu1,al1,opt,d)
int opt,n,npr,*mat;
float *pm,*pnu1,*al1,d[][3];
{
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];
*pnu1 = pnu;
*al1 = al;
/* --- d() matrix --- */
if (opt == 1) {
/* --- plane stress --- */
c1 = e / (1 - pnu * pnu);
c2 = c1 * pnu;
}
else {
/* --- plane strain --- */
c = e / ((1 + pnu) * (1 - 2 * pnu));
c1 = c * (1 - pnu);
c2 = c * pnu;
}
c3 = .5 * e / (1 + pnu);
d[0][0] = c1;
d[0][1] = c2;
d[0][2] = 0;
d[1][0] = c2;
d[1][1] = c1;
d[1][2] = 0;
d[2][0] = 0;
d[2][1] = 0;
d[2][2] = c3;
return(0);
}
elemstif(n,opt,se,tl,xni,d,thick,tempr,x,al,pnu,noc)
int n,opt,*noc;
float al,pnu;
float *x,*tempr,*thick,d[][3],tl[8],se[][8],xni[][2];
{
int i,j,k,ip;
float dte,c,xi,eta,th,dj,b[3][8],db[3][8];
/* ----- 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,thick,&th,d,b,db,&dj,xi,eta);
/* --- element stiffness matrix se --- */
for (i = 0; i < 8; i++) {
for (j = 0; j < 8; j++) {
c = 0;
for (k = 0; k < 3; k++) {
c = c + b[k][i] * db[k][j] * dj * th;
}
 se[i][j] = se[i][j] + c;
}
}
/* --- determine temperature load tl --- */
c = al * dte;
if (opt == 2)
c = (1 + pnu) * c;
for (i = 0; i < 8; i++) {
tl[i] = tl[i] + th * dj * c * (db[0][i] + db[1][i]);
}
}
return(0);
}
dbmat(n,x,noc,thick,th1,d,b,db,dj1,xi,eta)
float *x,*dj1,*thick,*th1,xi,eta;
float d[][3],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 th,a[3][4],g[4][8];
/* ----- db() matrix ----- */
/* --- nodal coordinates --- */
th = thick[n];
*th1 = th;
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;
/* --- b[3,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;
 }
}
/* --- db[3,8] matrix relates stresses to q[8] --- */
for (i = 0; i < 3; i++) {
for (j = 0; j < 8; j++) {
c = 0;
for (k = 0; k < 3; k++) {
c = c + d[i][k] * b[k][j];
}
db[i][j] = c;
}
}
return(0);
}
/*repeating functions for punch*/
elemstif2(n,opt,se,tl,xni,d,thick,tempr,x,al,pnu,noc)
int n,opt,*noc;
float al,pnu;
float *x,*tempr,*thick,d[][3],tl[8],se[][8],xni[][2];
{
int i,j,k,ip;
float dte,c,xi,eta,th,dj,b[3][8],db[3][8];
/* ----- 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];
dbmat2(n,x,noc,thick,&th,d,b,db,&dj,xi,eta);
/* --- element stiffness matrix se --- */
for (i = 0; i < 8; i++) {
for (j = 0; j < 8; j++) {
c = 0;
for (k = 0; k < 3; k++) {
c = c + b[k][i] * db[k][j] * dj * th;
}
se[i][j] = se[i][j] + c;
}
}
/* --- determine temperature load tl --- */
c = al * dte;
if (opt == 2)
c = (1 + pnu) * c;
for (i = 0; i < 8; i++) {
tl[i] = tl[i] + th * dj * c * (db[0][i] + db[1][i]);
}
}
return(0);
}
dbmat2(n,x,noc,thick,th1,d,b,db,dj1,xi,eta)
float *x,*dj1,*thick,*th1,xi,eta;
float d[][3],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 th,a[3][4],g[4][8];
/* ----- db() matrix ----- */
/* --- nodal coordinates --- */
th = thick[n];
*th1 = th;
n1 = noc[4*n];
n2 = noc[4*n+1];
n3 = noc[4*n+2];
n4 = noc[4*n+3];
x1 = x[2*(n1-1)-1792];
y1 = x[2*(n1-1)+1-1792];
x2 = x[2*(n2-1)-1792];
y2 = x[2*(n2-1)+1-1792];
x3 = x[2*(n3-1)-1792];
y3 = x[2*(n3-1)+1-1792];
x4 = x[2*(n4-1)-1792];
y4 = x[2*(n4-1)+1-1792];
/* --- 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;
/* --- b[3,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;
}
}
/* --- db[3,8] matrix relates stresses to q[8] --- */
for (i = 0; i < 3; i++) {
for (j = 0; j < 8; j++) {
c = 0;
for (k = 0; k < 3; k++) {
c = c + d[i][k] * b[k][j];
}
db[i][j] = c;
}
}
return(0);
}
/* ----- band solver ----- */
/* -----gauss solver ----- */
GSolve(float **a,int nq,float *x1)
{
int i,j,k,maxrow;
float tmp;

for (i=0;i<nq;i++) {

/* Find the row with the largest first value */

maxrow = i;
for (j=i+1;j<nq;j++) {
if (ABS(a[i][j]) > ABS(a[i][maxrow]))
maxrow = j;
}

/* Swap the maxrow and ith row */

for (k=i;k<nq+1;k++) {
tmp = a[k][i];
a[k][i] = a[k][maxrow];
a[k][maxrow] = tmp;
}
/* Eliminate the ith element of the jth row */
for (j=i+1;j<nq;j++) {
for (k=nq;k>=i;k--) {
a[k][j] -= a[k][i] * a[i][j] / a[i][i];
}
}
}

/* Do the back substitution */

for (j=nq-1;j>=0;j--) {
tmp = 0;
 for (k=j+1;k<nq;k++)
tmp += a[k][j] * x1[k];
x1[j] = (a[nq][j] - tmp) / a[j][j];
}
return(0);