#include "Stdio.h"
#include "Conio.h"
#include "Math.h"

void FLD(int *LD,int **ME,int NN, int N,int NP, int NE,int ND,int NF)
{
int NW,i,j,j1,k,ig;
LD[1]=1;
for(k=1;k<=NP;k++)
{
ig=NP+1;
for(i=1;i<=NE;i++)
for(j=1;j<=ND;j++)
{
if (ME[j][i]==k)
{
for(j1=1;j1<=ND;j1++)
{
if(ME[j1][i]<ig)
ig=ME[j1][i];
}
}
}
}
for(i=1;i<=NF;i++)
{
j=NF*(k-1)+i;
NW=NF*(k-ig)+i;
if(j!=1)
LD[j]=LD[j-1]+NW;
}
NN=LD[N];
printf("\nNN=%d ",NN);
for(i=1;i<=NF;i++)
printf("\n LD[i]=%d \n",LD[i]);
/* return LD[N]; */
}

void FIS(int **ME,int *IS,int NF,int ND,int IE)
{
int i,ii,j;
for(i=1;i<=ND;i++)
for(j=1;j<=NF;j++)
{
ii=NF*(i-1)+j;
IS[ii]=NF*(ME[i][IE]-1)+j;
}
}

void FKE(float *X,float *Y,float *AE,int **ME,float *T,int IE, int NDF,float E0,float **EK)
{
int NI,NJ,i,j;
float XIJ,YIJ,BL,C;
NI=ME[1][IE];
NJ=ME[2][IE];
XIJ=X[NJ]-X[NI];
YIJ=Y[NJ]-Y[NI];
BL=sqrt(XIJ*XIJ+YIJ*YIJ);
T[1]=XIJ/BL;
T[2]=YIJ/BL;
T[3]=-T[1];
T[4]=-T[2];
C=E0*AE[IE]*1.0/BL;
for(i=1;i<=NDF;i++)
for(j=1;j<=NDF;j++)
{
EK[i][j]=C*T[i]*T[j];
}
printf("\nIE=%d ",IE);
for(i=1;i<=NDF;i++)
{
printf("\n");
for(j=1;j<=NDF;j++)
{
printf("%5.2f",EK[i][j]);
}
}
/* return EK[4][4]; */
}

void FASUM(float **EK,int *IS,float *A,int *LD,int NDF)
{
int NI,NJ,i,j,ij;
for(i=1;i<=NDF;i++)
for(j=1;j<=NDF;j++)
{
NI=IS[i];
NJ=IS[j];
if(NJ<=NI)
{
ij=LD[NI]-(NI-NJ);
A[ij]=A[ij]+EK[i][j];
}
}
}

void FR(int *LD,float *A,int *NRR,int NR)
{
int i,N1,NI;
for(i=1;i<=NR;i++)
{
N1=NRR[i];
NI=LD[N1];
A[NI]=1.0E25;
}
}

MAX(int MI,int MJ)
{
int c;
if (MI>MJ)
c=MI;
else
c=MJ;
return c;
}

void BAND3(float *A,float **B,int *LD,int N,int M, int ISW)
{
int i,j,ij,k,LDI,IM1, MI,I0,J0,JM1,MIJ,MJ,IK,JK,ALDI,LDJ;
float s;
for(i=1;i<=N;i++)
{
LDI=LD[i];
if(i!=1)
{
I0=LDI-i;
IM1=i-1;
MI=1-I0+LD[IM1];
if(MI<=IM1)
{
for(j=MI;j<=IM1;j++)
{
J0=LD[j]-j;
ij=I0+j;
if(j!=1)
{
JM1=j-1;
MJ=1-J0+LD[JM1];
MIJ=MAX(MI,MJ);
if(MIJ<=JM1)
{
s=0.0;
for(k=MIJ;k<=JM1;k++)
{
IK=I0+k;
JK=J0+k;
s=s+A[IK]*A[JK];
}
A[ij]=A[ij]-s;
}
}
for(k=1;k<=M;k++)
B[i][k]=B[i][k]-A[ij]*B[j][k];
}
}
}
ALDI=A[LDI];
if(i!=1&&MI<=IM1)
{
for(j=MI;j<=IM1;j++)
{
ij=I0+j;
LDJ=LD[j];
s=A[ij];
A[ij]=s/A[LDJ];
ALDI=ALDI-s*A[ij];
}
}
A[LDI]=ALDI;
if(ALDI==0.0)
{
ISW=0;
}
else
{
for(k=1;k<=M;k++)
B[i][k]=B[i][k]/ALDI;
}
}
for(LDI=2;LDI<=N;LDI++)
{
i=N-LDI+2;
I0=LD[i]-i;
MI=1-I0+LD[i-1];
IM1=i-1;
if(MI<=IM1)
{
for(j=MI;j<=IM1;j++)
{
ij=I0+j;
for(k=1;k<=M;k++)
B[j][k]=B[j][k]-A[ij]*B[i][k];
}
}
}
ISW=1;
}

main()
{
int N,NP,NE,NL,NR,ND,NF,NDF,NN,NI,NJ,ISW,XIJ,YIJ,BL,i,j,LD[40],IS[4],U[4],ME[2][6]={{1,2,3,4,2,1},{2,3,4,1,4,3}},NRR[3]={1,2,7};
float EK[4][4],C,X[4]={0.0,10.0,10.0,0.0},Y[4]={10.0,10.0,0.0,0.0},AE[6]={1.0,1.0,1.0,1.0,1.414,1.414},E0,P[8][1]={0.0,0.0,0.0,0.0,0.0,-1000.0,00,0.0},A[40],S[6],T[4];
NP=4;
NE=6;
NL=1;
NR=3;
ND=2;
NF=2;
/* X[4]={0.0,10.0,10.0,0.0};
Y[4]={10.0,10.0,0.0,0.0};
ME[2][6]={{1,2,3,4,2,1},{2,3,4,1,4,3}};
AE[6]={1.0,1.0,1.0,1.0,1.414,1.414};
NRR[3]={1,2,7};
P[8][1]={0.0,0.0,0.0,0.0,0.0,-1000.0,0.0,0.0}; */
E0=1000000.0;

N=NF*NP;
NDF=ND*NF;
FLD(*LD,**ME,NN,N,NP,NE,ND,NF);

for(i=1;i<=NN;i++)
A[i]=0.0;
for(i=1;i<=NE;i++)
{
FIS(**ME,*IS,NF,ND,i);
FKE(*X,*Y,*AE,**ME,*T,i,NDF,E0,**EK);
FASUM(**EK,*IS,*A,*LD,NDF);
}
FR(*LD,*A,*NRR,NR);

for(i=1;i<=NN;i++)
{
printf("\nA[%d]=%5.2f \n ",i,A[i]);
}

BAND3(*A,**P,*LD,N,NL,ISW);

for(i=1;i<=NP;i++)
{
printf("%d %5.2f %5.2f\n",i,P[2*i-1][1],P[2*i][1]);
}

for(i=1;i<=NE;i++)
{
S[i]=0.0;
NI=ME[1][i];
NJ=ME[2][i];
XIJ=X[NJ]-X[NI];
YIJ=Y[NJ]-Y[NI];
BL=sqrt(XIJ*XIJ+YIJ*YIJ);
T[1]=XIJ/BL;
T[2]=YIJ/BL;
T[3]=-T[1];
T[4]=-T[2];
C=E0*AE[i]/BL;
for(j=1;j<=NF;j++)
{
IS[j]=NF*(NI-1)+j;
IS[j+2]=NF*(NJ-1)+j;
}
for(j=1;j<=NDF;j++)
{
U[j]=P[IS[j]][1];
S[i]=S[i]-C*T[j]*U[j];
printf("%d %5.2f \n",i,S[i]);
}
}
}