#include "cbox2.h"main(){	float X[82],Y[82],F[81],Bc[161],Xi[21],Yi[21],u[21],D;	int Code[161];	float G[81][162],H[81][81],q1[21],q2[21];	int Twice,Dim;		/*[Dim= Maximum dimension of the system AX = F; it is the same as the	 maximum number of nodes, or maximum number of elements ]*/			Dim = 80;	Twice = 2*Dim;			/*[	Read Data	]*/		Input2(Xi,Yi,X,Y,Code,Bc);		/*[ Compute matrices G and H, and form the system A X = F ]*/			Sys2(X,Y,G,H,F,Bc,Code,Dim,Twice);			/*[Solve the system AX = F ]*/		Solve(H,F,&D,N,Dim);	/*[ Compute the potential and fluxes at interior points ]*/	Inter2(F,Bc,Code,Xi,Yi,X,Y,u,q1,q2); 		/*[ Output the solution at boundary nodes and interior points ]*/		Out2(X,Y,F,Bc,Xi,Yi,u,q1,q2);}	void Input2(Xi,Yi,X,Y,Code,Bc)			float Xi[21],Yi[21],X[82],Y[82],Bc[161];	int Code[161];{	FILE  *infile, *outfile;	int   i,j,k,lnsize;	char Title[80],line1[100]; 	lnsize = 100;		printf("\nNOTE: FIRST LINE IN THE INPUT FILE SHOULD BE EITHER BLANK OR THE TITLE NOT 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);               infile = fopen(inname, "r" );        outfile = fopen(outname,"w" );				fgets(line1,lnsize,infile);		fprintf(outfile,"%s \n",line1);		fprintf(outfile,"Input \n");			fscanf(infile,"%d %d",&N,&L);				fprintf(outfile," \nNumber of Boundary Elements = %d \n",N);		fprintf(outfile," \nNumber of InteriorPoints = %d \n",L);			fprintf(outfile," \nCoordinates of the Extreme Points");		fprintf(outfile,"  of the  Boundary Elements  \n \n"); 		fprintf(outfile,"Point      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]);		}		/*[	Read boundary conditions in Bc[] vector. If Code[i] = 0, the Bc value is a known potential;  if Code[i] = 1, the Bc[] value is a known potential derivative (flux). Two Boundary conditions are read per element.  One node may have two values of the potential derivative but only one value of the potential ]*/	fprintf(outfile," \nBoundary  Conditions \n\n");	fprintf(outfile,"               First Node             Second Node\n");	fprintf(outfile,"               Prescribed             Prescribed\n");	fprintf(outfile,"Element     CODE   Value           CODE    Value\n");		for(i=1;i<=N;i++)	{		fscanf(infile,"%d %f %d %f",&Code[2*i -1],&Bc[2*i -1],					&Code[2*i],&Bc[2*i]);							fprintf(outfile,"%2d\t\t\t%2d %10.4f\t\t\t%2d %10.4f\n",				i,Code[2*i -1],Bc[2*i -1],Code[2*i],Bc[2*i]);	}	/*[ Read coordinates of the interior points ]*/	fprintf(outfile,"\nInterior Point Coordinates\n\n");	fprintf(outfile,"       Xi         Yi\n");			for(i=1;i<=L;i++)		{			fscanf(infile,"%f %f \n",&Xi[i],&Yi[i]);			fprintf(outfile," %10.4f %10.4f\n",Xi[i],Yi[i]);		}	fclose(infile);	fclose(outfile);			}void Inter2(F,Bc,Code,Xi,Yi,X,Y,u,q1,q2)	float F[81],Bc[161],Xi[21],Yi[21],X[82],Y[82],u[21],q1[21],q2[21];				int   Code[161];	{		float A1,B1,A2,B2,q11,q21,q12,q22,u11,u21,u12,u22,temp;		int i,j,k;				/*[ This Subroutine computes the values of the potential and the potential derivatives (fluxes) at the interior points ]*/ 					/*[ Rearrange the F and Bc arrays to store all values of the potential in F and all valus of the derivative in Bc ]*/						for(i=1;i<=N;i++)		{			for(j=1;j<=2;j++)			{				if( Code[2*i - 2+j] <= 0)				{					if(!((i != N) || (j!=2)))					{						if(Code[1] >0)						{							temp       = F[1];							F[1]    = Bc[2*N];							Bc[2*N] = temp;						}						else							Bc[2*N] = Bc[1];					}					else					{						if((i == 1) ||(j == 2)|| (Code[2*i -2] == 1))						{							temp 	= F[i-1+j];							F[i-1+j] = Bc[2*i-2+j];							Bc[2*i-2+j] = temp;						}						else							Bc[2*i-1] = Bc[2*i-2];					}				}			}		}					/*[ Compute the potential and fluxes at 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++)				{					Quad2(Xi[k],Yi[k],X[j],Y[j],X[j+1],Y[j+1],					&A1,&A2,&B1,&B2,&q11,&q21,&q12,&q22,					&u11,&u21,&u12,&u22,1);										if((j - N) < 0)					{						u[k]   += Bc[2*j-1]*B1 + Bc[2*j]*B2 - 									F[j]*A1 - F[j+1]*A2;						q1[k] += Bc[2*j-1]*u11+Bc[2*j]*u12-									F[j]*q11-F[j+1]*q12;							q2[k] += Bc[2*j-1]*u21+Bc[2*j]*u22-									F[j]*q21-F[j+1]*q22;					}					else					{						u[k]   += Bc[2*j-1]*B1+ Bc[2*j]*B2-F[j]									*A1-F[1]*A2;						q1[k] += Bc[2*j-1]*u11+Bc[2*j]*u12									- F[j]*q11-F[1]*q12;						q2[k] += Bc[2*j-1]*u21+Bc[2*j]*u22									- F[j]*q21-F[1]*q22;					}				}				u[k]   /= 2.0 * pi;				q1[k] /= 2.0 * pi;				q2[k] /= 2.0 * pi;  		 	}		} 	}	void Out2(X,Y,F,Bc,Xi,Yi,u,q1,q2)	float 	X[82],Y[82],F[81],Bc[161],Xi[21],Yi[21],u[21],q1[21],q2[21];{				FILE *outfile;		int i,j,k;			outfile = fopen(outname,"a");				fprintf(outfile," \n\n\nSolution: \n\n");		fprintf(outfile,"Boundary Nodes \n");				fprintf(outfile,"          (X,Y)               u       Prenodal q  Postnodal q\n");		fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f %10.4f\n",				X[1],Y[1],F[1],Bc[2*N],Bc[1]); 		for(i=2;i<=N;i++)			fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f %10.4f\n",				X[i],Y[i],F[i],Bc[2*i-2],Bc[2*i-1]); 		if(L)		{			fprintf(outfile,"\n\nInterior Points \n");			fprintf(outfile,"          (Xi,Yi)            u          q_x          q_y\n");			for(i=1;i<=L;i++)				fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f  %10.4f\n",Xi[i],Yi[i],u[i],q1[i],q2[i]);		}		fclose(outfile);}void Quad2(Xp,Yp,X1,Y1,X2,Y2,A1,A2,B1,B2,			q11,q21,q12,q22,u11,u21,u12,u22,K)	float  Xp,Yp,X1,Y1,X2,Y2,*A1,*A2,*B1,*B2;float *q11,*q21,*q12,*q22,*u11,*u21,*u12,*u22;int    K;{		/*[ This function computes the integral of several non-singular functions along the boundary elements using a four points Gauss Quadrature. 	When K = 0, The off-diagonal coefficients of the matrices H and G are computed. When K= 1, all coefficients needed for computation of the potential and fluxes at interior points are computed. Ra = Radius = Distance from the collocation point to the Gauss integration points on the boundary elements; nx,ny = Components of the unit normal to the element; rx,ry,rn = Radius derivatives. See section 5.3 ]*/float Xg[5],Yg[5];float Ax,Ay,Bx,By,HL,nx,ny,Ra,rx,ry,rn,F1,F2,G,H,ux,uy,qx,qy;int   i; float Z[] =  {0.0, 0.86113631, -0.86113631, 0.33998104, -0.33998104};float W[] = {0.0, 0.34785485,  0.34785485, 0.65214515,  0.65214515};	Ax = (X2 - X1)/2.0;	Bx = (X2 + X1)/2.0;	Ay = (Y2 - Y1)/2.0;	By = (Y2 + Y1)/2.0;	HL = sqrt(SQ(Ax) + SQ(Ay));	nx =  Ay/HL;	ny = -Ax/HL;	(*A1) = 0.0;	(*A2) = 0.0;	(*B1) = 0.0;	(*B2) = 0.0;	(*q11) = 0.0;	(*q21) = 0.0;	(*q12) = 0.0;	(*q22) = 0.0;	(*u11) = 0.0;	(*u21) = 0.0;	(*u12) = 0.0;	(*u22) = 0.0;	/*[ Compute the terms to be included in the matrices G and H, or the terms needed to compute the potential and fluxes at interior points ]*/	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;		F1 = ( 1.0 - Z[i])/2.0;		F2 = ( 1.0 + Z[i])/2.0;		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);			(*q11) += qx*F1;			(*q12) += qx*F2;			(*q21) += qy*F1;			(*q22) += qy*F2;			(*u11) += ux*F1;			(*u12) += ux*F2;			(*u21) += uy*F1;			(*u22) += uy*F2;		}		H = rn * W[i]*HL/Ra;		G = log(1.0/Ra)*W[i]*HL;		(*A1) -= F1*H;		(*A2) -= F2*H;		(*B1) += F1*G;		(*B2) += F2*G;	}}	void Diag2(X1,Y1,X2,Y2,B1,B2)float X1,Y1,X2,Y2,*B1,*B2; {	/*[ This function computes the parts of the G matrix coefficients for  integrals along an element that includes the collocation point ]*/	float Li;	 Li  = sqrt(SQ(X2 -X1) + SQ(Y2 - Y1));	(*B1) = Li*(1.5 - log(Li))/2.0;	(*B2) = Li*(0.5 - log(Li))/2.0;}/*[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 : Syestem MatrixB:  Originally it contains the independent coefficients.  After solutions it	contains the values of the syestem unknowns.                            				   				N:  Actual number of unknowns.Dim: Row and Column Dimension of A. ]*/void Solve(A,B,D,N,Dim)	float A[81][81],B[81],*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 1",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 2",k);	   	(*D) = 0.0;   }            return;	}	 void Sys2(X,Y,G,H,F,Bc,Code,Dim,No)	float X[82],Y[82],G[81][162],H[81][81],F[81],Bc[161];	int   Code[161],Dim,No;{	/*[ This function computes the matrices G and H, and forms the system  A X = F; H is a square matrix (N,N); G is a rectangular matrix (N,2*N) ]*/		int NN,i,j,k,NF,NS,jj;	float A1,B1,A2,B2,q11,q21,q12,q22,u11,u21,u12,u22,ch;			NN = 2*N;	for(i=1;i<=N;i++)	{		for(j=1;j<=N;j++)			H[i][j] = 0.0;		for(j=1;j<=NN;j++)			G[i][j] = 0.0;	}	/*[ Compute the coeffitients of G and H  ]*/			X[N+1] = X[1];		Y[N+1] = Y[1];				for(i=1;i<=N;i++)		{			NF = i+1;			NS = i+N-2;			for(jj=NF;jj<=NS;jj++)			{				if((jj - N) > 0)					j = jj-N;				else					j = jj;				Quad2(X[i],Y[i],X[j],Y[j],X[j+1],Y[j+1],&A1,&A2,&B1,&B2,				&q11,&q21,&q12,&q22,&u11,&u21,&u12,&u22,0);							if((j-N) <0)					H[i][j+1] += A2;				else					H[i][1]   += A2;								H[i][j]    += A1;				G[i][2*j-1] = B1;				G[i][2*j]   = B2;				H[i][i]	   -= A1 + A2;			}						NF = i+N-1;			NS = i+N;			for(jj=NF;jj<=NS;jj++)			{				if((jj-N) > 0)					j = jj- N;				else					j = jj;				Diag2(X[j],Y[j],X[j+1],Y[j+1],&B1,&B2);				if((jj - NF) <= 0)				{					ch = B1;					B1 = B2;					B2 = ch;				}				G[i][2*j-1] = B1;				G[i][2*j]	= B2;			}		/*[ Add one to the diagonal coefficients for exterior problems ]*/						if(H[i][i] < 0.0)					H[i][i] += 2. * pi;		}					/*[ Reorder the columns of the system of Equations as in (5.28)	and form the system Matrix A which is stored in H ]*/		for(i=1;i<=N;i++)	{		for(j=1;j<=2;j++)		{			if(Code[2*i -2+j] <= 0)			{				if(!((i != N) || (j != 2)))					{					if(Code[1] > 0)					{						for(k=1;k<=N;k++)						{							ch      = H[k][1];							H[k][1] = -G[k][2*N]; 							G[k][2*N] = -ch;						}					}					else					{						for(k=1;k<=N;k++)						{							H[k][1] -= G[k][2*N];							G[k][2*N] = 0.0;						}					}				}				else				{					if((i == 1) || (j >1) || (Code[2*i-2] == 1))					{						for(k=1;k<=N;k++)						{							ch 				= H[k][i-1+j];							H[k][i-1+j]  	= -G[k][2*i - 2+j];							G[k][2*i - 2+j] = -ch;						}					}					else					{						for(k=1;k<=N;k++)						{							H[k][i] -= G[k][2*i-1];							G[k][2*i-1] = 0.0;						}					}				}			}		}	}																									/*[ 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<=NN;j++)			F[i] += G[i][j] * Bc[j];	}  }