#include	<stdio.h>
#include	<math.h>

#define	NVE	4	/* number of nodal variables per element   */ 
#define	MAXNOD	410	/* maximum node number  */
#define	MAXELEM	730	/* maximum number of elements  */
#define	BANDMAX	30
#define	MAX(num1,num2) ((num1)>(num2)) ? num1:num2;
#define	MIN(num1,num2) ((num1)<(num2)) ? num1:num2;

FILE 	*fptr;
FILE 	*fred;
/*
external variables
*/
float 	x[MAXNOD];		/* x coordinates of nodes */
float 	y[MAXNOD];      	/* y coordinates of nodes */
float 	s[MAXNOD][BANDMAX];	/* total matrix of the system */
float 	rhs[MAXNOD];      	/* right hand side vector    */
float 	potent[MAXNOD];    	/* vector of solution */
int 	nodes;               	/* number of nodes used in problem */
int 	nelmts;              	/* number of elements in model */
int 	nvtx[NVE][MAXELEM];  	/* list of nodes for each element */
float 	source[MAXELEM];   	/* source density in each element */
int 	constr[MAXNOD];      	/* 1 - for fixed potential  */
float 	selm[NVE][NVE];    	/* s for one element  */
float 	telm[NVE][NVE];    	/* t for one element */
int 	intg[NVE];           	/* working array  */
int 	nband;               	/* halfbandwidth  */
char 	outfile[10];           	/* output file name */

/*-----------------------------------------------------------*/

main()
{
int i,j,k;

meshin();          /*   input  data */

detbw();          /*   calculate halfbandwidth  */

for (i=1; i<nelmts; i++) {
     elmatr(i);           
     asmbly(i);     
     }

eqsolv();           /* solve set of equations  */

fputs("\n\n\t Final solution   ",fptr);
fputs("\n\n\tNode  x        y     Potential \n",fptr);

for(i=1; i<=nodes; i++)
fprintf(fptr,"%9d %8.4f %8.4f %10.4e \n",i,x[i],y[i],potent[i]);
i=0;
fprintf(fptr,"%9d %8.4f %8.4f %10.4e \n",i,x[i],y[i],potent[i]);

fclose(fptr);

}

/*------------------------------------------------------*/

meshin() 
{
char infile[10];
int i;
float sourci,potenti,xx,yy;

/*  nodes: number of node, x ,  y   */
/*  elements:  node numbers, source */
/*  boundary conditions:  node number, fixed value  */

printf("Enter input file name ?");
scanf("%s", infile);
fred=fopen(infile,"r");

printf("Enter output file name ? ");
scanf("%s", outfile);
fptr=fopen(outfile,"w");

/* read in the node list  */

nodes=-1;
fputs("\n\n\t   Input node list \n",fptr);
fputs("\n\tn      x       y\n",fptr);

do  {
    fscanf(fred, "%d %f %f", &i, &xx , &yy  );
    x[i]=xx;
    y[i]=yy;
    ++nodes;
    if(i!=0) 
          fprintf(fptr,"%9.0d %8.3f %8.3f \n",i,xx,yy);
    }
while(i!=0);
    /* number of nodes acceptable ? */
    if(nodes>MAXNOD)
    printf("Node number %d  exceeds maximum  %d\n",nodes,MAXNOD);

    /*  read in the element list  */

    fputs("\n\n\t  Input element list \n",fptr);
    fputs("\n\ti        j        k      source\n" ,fptr);
nelmts= 0;
do {
     fscanf(fred,"%d %d %d %f", &intg[1],&intg[2],&intg[3],&sourci);
     ++nelmts;
     source[nelmts]=sourci;
     nvtx[1][nelmts]= intg[1];
     nvtx[2][nelmts]= intg[2];
     nvtx[3][nelmts]= intg[3];
     if(intg[1]!=0) 
	fprintf(fptr,"%9d %8d %8d %12.4f\n", intg[1],intg[2],intg[3],sourci);
     }
while(intg[1]!=0);
nelmts--;
          /*  total number of elemts acceptable ? */
if(nelmts>MAXELEM)
     printf("Element number %d exceeds maximum %d\n",nelmts,MAXELEM);

         /*  read boundary conditions */

fputs("\n\n\t Input fixed potentials \n",fptr);
fputs("\n\tnode value\n",fptr);

do {
     fscanf(fred,"%d %f",&i,&xx);
     potent[i]=xx;
     constr[i]=1;
     fprintf(fptr,"%9d %8.3f\n",i,xx);
     }
while(i!=0);

fclose(fred);
}
/*---------------------------------------------------------*/
detbw()
{
int i,bn,sn,bw;

for(i=1; i<=nelmts; ++i) {
     bn=MAX(nvtx[1][i],nvtx[2][i]);
     bn=MAX(bn,nvtx[3][i]);

     sn=MIN(nvtx[1][i],nvtx[2][i]);
     sn=MIN(sn,nvtx[3][i]);

     bw=bn-sn+1;
     if(bw > nband)
               nband=bw;
     }

fprintf(fptr,"\n\n\n  halfbandwidth = %d", nband);
}
/*--------------------------------------------------------------*/
elmatr(ie)
int ie;
{
float xi,xj,xk,yi,yj,yk;
int i,j,k,l,m,i1,i2,i3,nvrtex,i4;
float area,ctng,ctng2;

i=nvtx[1][ie];
j=nvtx[2][ie];
k=nvtx[3][ie];

area=((x[j]-x[i])*(y[k]-y[i])-(x[k]-x[i])*(y[j]-y[i]))/2.0;

for(l=1; l<=3; l++) {
     for(m=1; m<=3; m++)
            telm[l][m]=area/12.0;
     telm[l][l]=area/6.0;
     }

i1=1;
i2=2;
i3=3;

for(l=1; l<=3; l++)
for(m=1; m<=3; m++)
      selm[l][m]=0.0;

for(nvrtex=1; nvrtex<=3; nvrtex++) {
     ctng2=((x[j]-x[i])*(x[k]-x[i])+(y[j]-y[i])*(y[k]-y[i]))/4.0/area;

     selm[i2][i2]+=ctng2;
     selm[i2][i3]-=ctng2;
     selm[i3][i2]-=ctng2;
     selm[i3][i3]+=ctng2;

     i4=i1;
     i1=i2;
     i2=i3;
     i3=i4;

     l=i;
     i=j;
     j=k;
     k=l;
     }
}
/*-----------------------------------------------------------*/
asmbly(ie)
int ie;
{
int i,irow,j,icol;

for(i=1; i<=3; i++) {
     irow=nvtx[i][ie];
     if( constr[irow] == 1 ) {
             s[irow][1]=1.0;
             rhs[irow]=potent[irow];
             }
     else {
             for(j=1; j<=3; j++) {
                      icol=nvtx[j][ie];
                      if( constr[icol]==1 ) 
                         rhs[irow]+=telm[i][j]*source[ie]
                                   -selm[i][j]*potent[icol];
                      else {
                         rhs[irow]+=telm[i][j]*source[ie];
                         if( icol >= irow )
                              s[irow][icol-irow+1]+=selm[i][j];
                         }
                      }
              }
      }
}

/*------------------------------------------------------*/
eqsolv()
{
int n,i,j,k,l;
float q;

for(n=1; n<=nodes; n++) {
     i=n;
     for( l=2; l<=nband; l++) {
          i++;
          if(s[n][l]!=0.0) {
               q=s[n][l]/s[n][1];
               j=0;
               for(k=l; k<=nband; k++) {
                    j++;
                    if(s[n][k]!=0.0)
                         s[i][j]-=q*s[n][k];
                    }
               s[n][l]=q;
               rhs[i]-=q*rhs[n];
               }
           }
           rhs[n]/=s[n][1];
      }
      for(n=nodes-1; n>0; n--) {
          l=n;
          for(k=2; k<=nband; k++) {
               l++;
               rhs[n]-=s[n][k]*rhs[l];
               }
          }
      for(n=1; n<=nodes; n++)
      potent[n]=rhs[n];

}