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



			Finite Element Program for

			  Poison’s Equation using

		      First Order Triangular Elements

	      (Modification of Silvester and Ferrari Program)



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

#include	<stdio.h>

#include	<math.h>



#define	NVE	4	/* number of nodal variables per element */ 

#define	MAXNOD	410	/* maximum node number (free + fixed) */

#define	MAXELEM	730	/* maximum number of elements */

#define	BANDMAX	30	/* maximum free node number */



#define	MAX(num1,num2) ((num1)>(num2)) ? num1:num2;

#define	MIN(num1,num2) ((num1)<(num2)) ? num1:num2;



FILE 	*fptrout;		/* output file pointer */

FILE 	*fptrin;		/* input file pointer */



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

external variables

Note: all external variables are automatically initialized with zero

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



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 	nband;               	/* halfbandwidth  */

int 	nvtx[NVE][MAXELEM];  	/* list of nodes for each element */



float 	source[MAXELEM];   			/* source density in each element */

int 	constr[MAXNOD];      		/* set to "1" for fixed potential  */



float 	selm[NVE][NVE];    			/* s for one triangular element  */

float 	telm[NVE][NVE];    			/* t for one triangular element */



int 	intg[NVE];           		/* working array  */

char 	outfile[10];         	/* output file name */



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

main program

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



main()

{

	int	i, j, k;

/*

	input data from file

*/

	meshin();

/*

	calculate halfbandwidth

*/

	detbw();

/*

	calculate each elemental contribution to the matrix via elmatr( i)

	assemble into global matrix via asmbly( i)

*/

	for( i=1; i<nelmts; i++){

		elmatr( i);           

		asmbly( i);     

	}

/*

	print out the matrix equation

*/

	fputs( "\n\n\t Matrix s and rhs  \n\n", fptrout);

	for( i=1; i<=nodes; i++){

		for( j=1; j<=nodes; j++)

			fprintf( fptrout,"%8.4f ", s[i][j]);

		fprintf( fptrout,"\t%8.4f\n", rhs[i]);

	}

/*

	solve set of equations

*/

	eqsolv();



	fputs( "\n\n\t Final solution   ", fptrout);

	fputs( "\n\n\tNode  x        y     Potential \n", fptrout);



	for( i=1; i<=nodes; i++)

	   fprintf( fptrout,"%9d %8.4f %8.4f %10.4e \n",

			      i, x[i], y[i], potent[i]);



	i=0;

	fprintf( fptrout,"%9d %8.4f %8.4f %10.4e \n",

				i, x[i], y[i], potent[i]);

/*

	close file

*/

	fclose( fptrout);



}



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

function meshin(): reads the user specified input file of data

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



meshin() 

{

	char	infile[10];

	int	i;

	float	sourci, potenti, xx, yy;

/*

	format of input data

	--------------------

	nodes:			number of node, x ,  y

				0 0 0

	elements:		triangle node numbers, source

				0 0 0 0

	boundary conditions:	node number, fixed value

				0 0

*/



	printf( "Enter input file name ?");

	scanf( "%s", infile);

	fptrin = fopen( infile,"r");



	printf( "Enter output file name ? ");

	scanf( "%s", outfile);

	fptrout = fopen( outfile,"w");

/*

	read in the node list

	set:	x[i] to x coordinate value of node i

		y[i] to y coordinate value of node i

	increment number of nodes "nodes"

	print:	node number and coordinates to output file

*/

	nodes = -1;

	fputs( "\n\n\t   Input node list \n", fptrout);

	fputs( "\n\tn      x       y\n", fptrout);



	do{

		fscanf( fptrin, "%d %f %f", &i, &xx , &yy);

		x[i] = xx;

		y[i] = yy;

		++nodes;

		if( i!=0) 

			fprintf( fptrout,"%9.0d %8.3f %8.3f \n",

							i, xx, yy);

	} while( i!=0);

/*

	check that total number of "nodes" is acceptable ?

*/

	if( nodes>MAXNOD)

		printf( "Node number %d  exceeds maximum  %d\n",

							nodes, MAXNOD);

/*

	read in the element list

	set:	source[nelmts] to constant source in element

		nvtx[i][j] to node number for ith node of triangle j

	print:	node numbers and source value of triangle to output file

*/

	fputs( "\n\n\t  Input element list \n", fptrout);

	fputs( "\n\ti        j        k      source\n", fptrout);

	nelmts = 0;



	do{

		fscanf( fptrin,"%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( fptrout,"%9d %8d %8d %12.4f\n",

				  intg[1], intg[2], intg[3], sourci);



	} while( intg[1]!=0);



	nelmts--;

/*

	check that total number of elements "nelmts" is acceptable ?

*/

	if( nelmts>MAXELEM)

	   printf( "Element number %d exceeds maximum %d\n",

						nelmts, MAXELEM);

/*

	read Dirichlet boundary conditions from input file

	set:	potent[i] to boundary value

		constr[i] to 1 implying fixed node

	print:	node number and value to output file

*/

	fputs( "\n\n\t Input fixed potentials \n", fptrout);

	fputs( "\n\tnode value\n", fptrout);



	do{

		fscanf( fptrin,"%d %f", &i, &xx);

		potent[i] = xx;

		constr[i] = 1;

		fprintf( fptrout, "%9d %8.3f\n", i, xx);



	} while( i!=0);

/*

	close input file

*/

	fclose( fptrin);

}



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

function detbw(): determines the bandwidth of the global matrix "s"

		  in order to make matrix solution more efficient

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



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( fptrout,"\n\n\n  halfbandwidth = %d", nband);

}



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

function elmatr( ie): determines the elemental matrices "telm" and "selm" for

		      triangular element "ie"

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



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];

/*

	calculate area of elemental triangle

*/

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

/*

	determine elemental "telm" matrix

*/

	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;

/*

	set "selm" to zero

*/

	for( l=1; l<=3; l++)

		for( m=1; m<=3; m++)

			selm[l][m] = 0.0;

/*

	determine "selm"

*/

	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;

	}

}



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

function asmbly( ie): inserts elemental matrix for triangle

		      ie into global matrix s and column matrix rhs

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



asmbly( ie)

int	ie;

{

	int	i, irow, j, icol;

/*

	loop over each node of triangle

*/

	for( i=1; i<=3; i++){

/*

		set global row number "irow" to "ith"

		node number of triangle "ie"

*/

		irow = nvtx[i][ie];

/*

		if irow is a constrained potential (i.e. boundary row)

		then put a "1" in the first column of "s" and set "rhs"

		to value of constrained potential

*/

		if( constr[irow] == 1 ){

			s[irow][1] = 1.0;

			rhs[irow] = potent[irow];

		}

/*

		else potential is free to vary

*/

		else{

/*

			loop over each node of triangle

*/

			for( j=1; j<=3; j++){

/*

				set global column number "icol" to "jth"

				node number of triangle "ie"

*/

				icol = nvtx[j][ie];

/*

				if icol is a constrained potential bring

				selm*potent over to rhs

*/

				if( constr[icol] == 1) 

					rhs[irow] += telm[i][j]*source[ie]

						    -selm[i][j]*potent[icol];

/*

				else insert "selm" into lower triangular

				part of "s" only since "s" is symmetric and

				determine "rhs" from "telm"

*/

				else{

					rhs[irow] += telm[i][j]*source[ie];

					if( icol >= irow )

					   s[irow][icol-irow+1]+=selm[i][j];

				}

			}

		}

	}

}



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

function eqsolv(): solves the matrix equation "s[][]*potent[] = rhs[]"

		   for the potentials potent[]

		   "s" is a banded upper triangular matrix representing

		   the symmetric stiffness matrix of finite element

		   calculations

		   Gaussian ellimination is used:

		   1) Forward Ellimination

		   2) Back Substitution

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



eqsolv()

{

	int	n, i, j, k, l;

	float	q;

/*

	Forward Ellimination (banded symmetric matrix):

	half bandwidth is "nband" and only upper triangular

	portion of "s" is stored i.e. "s[n][1] is diagonal

*/

	for( n=1; n<=nodes; n++){

		i = n;

/*

		proceed across the row --> down the column

		due to symmetry of "s" and make entries zero

*/

		for( l=2; l<=nband; l++){

			i++;

/*

			if entry is not zero then proceed to make it zero

*/

			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];

			}

		}

/*

		after each row in column "n" is cleared divide

		"rhs" by diagonal (really part of back substitution)

		effectively makes diagonal equal to 1

*/

		rhs[n] /= s[n][1];

	}

/*

	Back Substitution

*/

	for( n=nodes-1; n>0; n--){

		l = n;

/*

		subtract all components already found

*/

		for( k=2; k<=nband; k++){

			l++;

			rhs[n] -= s[n][k]*rhs[l];

		}

	}

/*

	potentials "potent[]" are the new right hand side "rhs"

*/

	for( n=1; n<=nodes; n++)

		potent[n] = rhs[n];

}