#include "cbox5.h"main(){	float X[102],Y[102],Xm[101],Ym[101],Bc[101],F[101],Xi[21],Yi[21];	float u[21],q1[21],q2[21],D;	float G[101][101],H[101][101];	int Code[101],Dim;		/*[ Dim =Max dimension of the system AX = F. This number		must be equal or smaller than the dimension of Xm ]*/			Dim = 100;		/*[ Read input  ]*/	Input5(Xi,Yi,X,Y,Code,Bc);		/*[ Compute matrices G and H, and form the system A X = F ]*/	Sys5(X,Y,Xm,Ym,G,H,Bc,F,Code,Dim);			/*[Solve the system AX = F ]*/	Solve(G,F,&D,N,Dim);		/*[ Compute the potential at interior points ]*/	Inter5(Bc,F,Code,Xi,Yi,X,Y,u,q1,q2);		/*[ Output solution at boundary nodes and interior points ]*/	Out5(Xm,Ym,Bc,F,Xi,Yi,u,q1,q2);}void  Input5(Xi,Yi,X,Y,Code,Bc)		float Xi[21],Yi[21],X[102],Y[102],Bc[101];	int Code[101];	{		FILE  *infile, *outfile;		int i,j,k,lnsize;		char title[80],line1[100];		lnsize = 100; 		printf("\nNOTE:FIRST LINE IN INPUT FILE SHOULD BE THE TITLE OR LEFT BLANK NO DATA\n");        		printf("\n Enter the name of the input file:");        scanf("%s", inname);        printf("\n Enter the name of the output file:");        scanf("%s", outname);        printf("\n");               infile = fopen(inname, "r" );        outfile = fopen(outname,"w" );				fgets(line1,lnsize,infile);		fprintf(outfile,"%s\nInput \n",line1);						fscanf(infile,"%d %d %d",&N,&L,&M);		fprintf(outfile," \nNumber of Boundary Elements = %d \n",N);		fprintf(outfile,"Number of Interior Points = %d \n",L);								if(M > 0)		{			fprintf(outfile,"Number of different Boundaries = %d \n \n",M);			for(i=1;i<=M;i++)			{				fscanf(infile,"%d",&Last[i]);				fprintf(outfile,"Last node on  boundary %2d = %2d \n",i,Last[i]);			}		}		fprintf(outfile," \nCordinates of The Extreme"); 		fprintf(outfile," Points of the Boundary Elements \n"); 		fprintf(outfile,"\nPoint    X          Y \n");		for(i=1;i<=N;i++)		 { 		   fscanf(infile,"%f %f",&X[i],&Y[i]);		   fprintf(outfile,"%2d %10.4f %10.4f \n",i,X[i],Y[i]);		 }		fprintf(outfile," \n \n");		  /*[ Read boundary conditions in Bc[i] vector; if Code[i] = 0,  the Bc[i] value is a known potential; if Code[i] = 1, the Bc[i] value is a known potential derivative (flux).]*/		fprintf(outfile,"Boundary Conditions \n");		fprintf(outfile,"\nNode		Code	Prescribed Value \n");		for(i=1;i<=N;i++)		{		  fscanf(infile,"%d %f",&Code[i],&Bc[i]);		  fprintf(outfile,"%2d \t\t\t%2d \t\t\t%10.4f \n",i,Code[i],Bc[i]);		}				fprintf(outfile," \n \n");/*[  Read coordinates of the interior points   ]*/		if( L)		 {		   fprintf(outfile," \nInterior Point Coordinates \n");		   fprintf(outfile,"\nPoint    Xi         Yi\n"); 		   for(i=1;i<=L;i++)		    {		     fscanf(infile,"%f %f",&Xi[i],&Yi[i]);			 fprintf(outfile,"%2d %10.4f %10.4f\n",i,Xi[i],Yi[i]);		    }		  }    		for(i=1;i<=80;i++)		fprintf(outfile,"*");		fprintf(outfile," \n");  		fclose(infile);		fclose(outfile);  }void Inter5(Bc,F,Code,Xi,Yi,X,Y,u,q1,q2)float Bc[101],F[101],Xi[21],Yi[21];float X[102],Y[102],u[21],q1[21],q2[21];int Code[101];{	int i,j,k,lk,kk,found;	float A,B,qx,qy,ux,uy,temp;		found = 0; 	/*[	Initialization ]*/	/*[ Rearrange the Bc and F arrays to store all values of the		potential in Bc and all potential derivatives in F ]*/				for(i=1;i<=N;i++)	{		if(Code[i] > 0)		{			temp = Bc[i];			Bc[i] = F[i];			F[i] = temp;		}	}	/*[ Compute the potential and fluxes at the interior points ]*/		if(L)	{		for(k=1;k<=L;k++)		{			u[k] = 0.0;			q1[k] = 0.0;			q2[k] = 0.0;			for(j=1;j<=N;j++)			{				if((M-1) > 0)				{					if(!(j-Last[1]))						kk = 1;					else						{							found = 0;							for(lk=2;lk<=M;lk++)							{								if(!(j-Last[lk]))								{									kk = Last[lk-1] + 1;									found = 1;									break;								}							}							if(!found)								kk = j + 1;														}				}				else					kk = j + 1;				Quad5(Xi[k],Yi[k],X[j],Y[j],X[kk],Y[kk],						&A,&B,&qx,&qy,&ux,&uy,1); 				u[k]   += F[j]*B - Bc[j]*A;				q1[k] += F[j]*ux - Bc[j]* qx;				q2[k] += F[j]*uy - Bc[j]* qy;			}			u[k]   /= 2 * pi;			q1[k] /= 2 * pi;			q2[k] /= 2 * pi;		}	}}void Out5(Xm,Ym,Bc,F,Xi,Yi,u,q1,q2)float Xm[101],Ym[101],Bc[101],F[101],Xi[21],Yi[21];float u[21],q1[21],q2[21];{	FILE *outfile;	int i,j,k;		outfile = fopen(outname,"a");			fprintf(outfile," \nResults \n\n");		fprintf(outfile,"\n            (X,Y)              u         q\n");		for(i=1;i<=N;i++)		fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f\n",					Xm[i],Ym[i],Bc[i],F[i]);						if(L)	{		fprintf(outfile,"\n\n\n		Interior Points \n\n");		fprintf(outfile,"          (Xi,Yi)              u         q_x       q_y\n");		for(k=1;k<=L;k++)			fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f %10.4f\n",						Xi[k],Yi[k],u[k],q1[k],q2[k]);	}			 	fclose(outfile);}  void Quad5(Xp,Yp,X1,Y1,X2,Y2,H,G,qx,qy,ux,uy,K)	float Xp,Yp,X1,Y1,X2,Y2,*H,*G;	float *qx,*qy,*ux,*uy;	int K;{    /*[ Program Quad5 :This function computes the integral of several nonsingular functions along the boundary elements using a four points Gauss Quadrature.If K = 0, the off-diagonal coefficients of the H and G matrices are computed; when K = 1, all coefficients needed to compute the potential and fluxes at the interior points are computed.Ra = Radius = Distance from the collocation point to the Gauss Integration points on the boundary element; nx,ny = components of the unit normal to the element; rx,ry,rn = Radius derivatives. See section 5.3.       ]*/const float Z[] = { 0.0,0.86113631,-0.86113631,0.33998104,-0.33998104};const float W[] = {0.0,0.34785485, 0.34785485,0.65214515, 0.65214515};float Xg[5],Yg[5];float Ax,Bx,Ay,By,HL,nx,ny,Ra,rx,ry,rn;int i;	Ax = (X2 - X1)/2;	Bx = (X2 + X1)/2;	Ay = (Y2 - Y1)/2;	By = (Y2 + Y1)/2;	HL = sqrt(SQ(Ax) + SQ(Ay));	nx =  Ay/HL;	ny = -Ax/HL;	(*G) = 0.0;	(*H) = 0.0;	(*ux)  = 0.0;	(*uy)  = 0.0;	(*qx)  = 0.0;	(*qy)  = 0.0;/*[ Compute G[x][y], H[x][y], qx, qy, ux, and uy ]*/  							for(i=1;i<=4;i++) 		{		  Xg[i] = Ax*Z[i] + Bx;		  Yg[i] = Ay*Z[i] + By;		  Ra = sqrt( SQ(Xp - Xg[i]) + SQ(Yp - Yg[i]));		  rx = (Xg[i]-Xp)/Ra;		  ry = (Yg[i]-Yp)/Ra;		  rn = rx*nx + ry*ny;		  if(K > 0)		    {			  (*ux) += rx*W[i]*HL/Ra;			  (*uy) += ry*W[i]*HL/Ra;			  (*qx) -= (( 2.0 * SQ(rx) -1.0)*nx + 2.0 *rx*ry*ny)*			  		  	W[i]*HL/SQ(Ra);			  			  (*qy) -= (( 2.0 * SQ(ry) -1.0)*ny + 2.0 *rx*ry*nx)*			  		  	W[i]*HL/SQ(Ra);			 }			 			  (*G) += log(1/Ra) * W[i]*HL;			  (*H) -= rn*W[i]*HL/Ra;		  }}void Diag5(X1,Y1,X2,Y2,G)float X1,Y1,X2,Y2,*G;{  float Ax,Ay,Li;  Ax = (X2 - X1)/2.0;  Ay = (Y2 - Y1)/2.0;  Li = sqrt(SQ(Ax) + SQ(Ay));  (*G) = 2* Li*(1.0 - log(Li)); }/*[ function  Solve : Solution of the linear systems of equations by the Gauss elimination method providing for interchanging rows when encountering a zero diagonal coeficient.A : System MatrixB:  Originally it contains the independent coefficients.  After solutions it	contains the values of the system unknowns.                            				   				N:  Actual number of unknowns.Dim: Row and Column Dimension of A. ]*/void Solve(A,B,D,N,Dim)	float A[101][101],B[101],*D;  /*[ D is assigned some value here ]*/  	int N,Dim;{	int N1,i,j,k,i1,l,k1,found;	float c;	/*[ found is a flag which is used to check if any non zero coeff is found ]*/				N1 = N - 1;	for(k=1;k<=N1;k++)	{		k1 = k + 1;		c = A[k][k];				if( (fabs(c) - tol) <=0 )		{			found = 0;			for(j=k1;j<=N;j++)  /*[Interchange rows to get Nonzero ]*/		        {            	if( (fabs(A[j][k]) - tol) > 0.0 )                 {				 	for(l=k;l<=N;l++)                    {                     	c = A[k][l];                        A[k][l] = A[j][l];                        A[j][l] = c;                     }				  		c = B[k];				  		B[k] = B[j];				  		B[j] = c;				  		c = A[k][k];						found = 1;		/*[ coeff is found ]*/						break;				}			}		 }		if(!found)		{			printf("Singularity in Row %d \n",k);	   		(*D) = 0.0;          	    	return;   /*[ If no coeff is found the control is transferred to main ]*/		}		/*[ Divide row by diagonal coefficient ]*/		  		c = A[k][k];		for(j = k1;j<=N;j++)		   A[k][j] /= c;		B[k] /= c;   /*[ Eliminate unknown X[k] from row i ]*/		for(i = k1;i<=N; i++)		{		  c = A[i][k];		  for(j = k1;j<= N;j++)		 	A[i][j] -= c*A[k][j];		  B[i] -= c*B[k];		}	}/*[ Compute the last unknown ]*/	if((fabs(A[N][N]) - tol) > 0.0)	{	   B[N] /= A[N][N];/*[Apply back substitution to compute the remaining unknowns  ]*/	   for(l = 1;l<= N1;l++)	   {	      k = N - l;	      k1 = k +1;	      for(j = k1;j<=N;j++)			B[k] -= A[k][j]*B[j];	   }/*[Compute the value of the determinent ]*/	(*D) = 1.0;	for(i=1;i<=N;i++)	  (*D) *= A[i][i];   }  else  {  		printf("Singularity in Row %d \n",k);	   	(*D) = 0.0;   }            return;	}void Sys5(X,Y,Xm,Ym,G,H,Bc,F,Code,Dim)float X[102],Y[102],Xm[101],Ym[101],G[101][101];float H[101][101],Bc[101],F[101];int Code[101],Dim;{	int i,j,k,kk,found;	float qx,qy,ux,uy,temp;		found = 0;	/*[ Initialization ]*/		/*[ This function computes the matrices G and H. and forms the system	 A X = F ]*/				X[N+1] = X[1];		Y[N+1] = Y[1];		for(i=1;i<=N;i++)		{			Xm[i] = (X[i] + X[i+1])/2.0;			Ym[i] = (Y[i] + Y[i+1])/2.0;		}				if((M-1) > 0)		{			Xm[Last[1]] = (X[Last[1]]  + X[1])/2.0;			Ym[Last[1]] = (Y[Last[1]]  + Y[1])/2.0;				for(k=2;k<=M;k++)			{				Xm[Last[k]] = (X[Last[k]]  + X[Last[k-1]+1])/2.0;				Ym[Last[k]] = (Y[Last[k]]  + Y[Last[k-1]+1])/2.0;			}		}			/*[ Compute the coefficients of G and H matrices ]*/			for(i=1;i<=N;i++)		{			for(j=1;j<=N;j++)			{				if((M-1) > 0)				{					if(!(j-Last[1]))						kk = 1;					else						{								found = 0;							for(k=2;k<=M;k++)							{								if(!(j-Last[k]))								{									kk = Last[k-1] + 1;									found = 1;									break;								}							}							if(!found)								kk = j + 1;						}				}				else					kk = j + 1;				if(i-j)				{						Quad5(Xm[i],Ym[i],X[j],Y[j],X[kk],Y[kk],&H[i][j],					&G[i][j],&qx,&qy,&ux,&uy,0);									}				else				{						Diag5(X[j],Y[j],X[kk],Y[kk],&G[i][j]);					H[i][j] = pi;  				}												}		}	/*[ Reorder the columns of the system of equations as in (5.28) 	 and form the system matrix A which is stored G]*/		for(j=1;j<=N;j++)	{		if(Code[j] > 0)		{			for(i=1;i<=N;i++)			{				temp = G[i][j];				G[i][j] = -H[i][j];				H[i][j] = -temp;			}		}	}	/*[ Form the right-side vector F which is stored in F ]*/		for(i=1;i<=N;i++)	{		F[i] = 0.0;		for(j=1;j<=N;j++)			F[i] += H[i][j] * Bc[j];	}}