/*[ Program Be1	]*/  #include"cbox1.h"main(){float X[102],Y[102],Xm[101],Ym[101],Bc[101],F[101],Xi[21],Yi[21],u[21];	float q1[21],q2[21],G[100][100],H[100][100],D;	int Code[101],Dim;/*[ Dim is a maximum dimension of the system of equations AX = F; it must be at most equal the dimension of (Xm, Ym)  ]*/	Dim = 100;/*[ Read input file ]*/	Input1(X,Y,Xi,Yi,Code,Bc);/*[ Compute G and H matrices and form the system A X = F ]*/   	Sys1(X,Y,Xm,Ym,G,H,Bc,F,Code,Dim);  /*[ Solve  system AX = F  ]*/	Solve(G,F,&D,N,Dim);/*[ Compute the potential and fluxes at interior points  ]*/	Inter1(Bc,F,Code,Xi,Yi,X,Y,u,q1,q2);/*[ Output solution at the nodes and interior points  ]*/	Out1(Xm,Ym,Bc,F,Xi,Yi,u,q1,q2);}void Input1(X,Y,Xi,Yi,Code,Bc)	float X[102], Y[102],Xi[21],Yi[21],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 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);        printf("\n");               infile = fopen(inname, "r" );        outfile = fopen(outname,"w" );	fgets(line1,lnsize,infile);	fprintf(outfile,"%s \n",line1);	fprintf(outfile," Input \n");/*[  Read boundary elements and internal points ]*/ 	fscanf(infile,"%d %d",&N,&L);	fprintf(outfile," \nNumber of Boundary Elements = %2d \n",N);	fprintf(outfile,"Number of Interior Points = %2d \n",L);/*[ Read coordinates of extreme points in arrays X and Y  ]*/	fprintf(outfile," \nCordinates of Extreme Points \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 the vector Bc[i].]*/	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  %2d %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     %2.4f %2.4f \n",i,Xi[i],Yi[i]);		    }	}    	for(i=1;i<=80;i++)	fprintf(outfile,"*");	fprintf(outfile," \n");  	fclose(infile);	fclose(outfile);  }void Sys1(X,Y,XM,YM,G,H,Bc,F,Code,Dim)  float X[102],Y[102],XM[101],YM[101],Bc[101],F[101],G[100][100],H[100][100];  int Code[101],Dim;{	float temp,DQ1,DQ2,DU1,DU2;	int  i,j,l,m;   /*[ Compute coordinates of nodes and store in arrays XM and YM ]*/	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;    }/*[ Compute the coeficients of matrices G and H  ]*/	for(i=1;i<=N;i++)	{	 	for(j=1;j<=N;j++)	  	{	    	m = j+1;			if((i-j))		 		Quad1(XM[i],YM[i],X[j],Y[j],X[m],Y[m],&H[i][j],&G[i][j],&DQ1,&DQ2,&DU1,&DU2,0);		 else		 {		  Diag1(X[j],Y[j],X[m],Y[m],&G[i][j]);		  H[i][j] = pi;		  	}		}	  }/*[ Reorder the columns of the system AX = F in accordance with the boundary conditions and form system matrix A which is stored in G ]*/	for(i=1;i<=N;i++)	{	 	if(Code[i] > 0)	  	{	    	for(j=1;j<=N;j++)			{		  	temp = G[j][i];		  	G[j][i] = -H[j][i];		  	H[j][i] = -temp;		 	}	  	} 	 }/*[ Form the right-side Vector F which is Stored in F ]*/	for(i=1;i<=N;i++)	{	 	F[i] = 0;	 	for(j=1;j<=N;j++)	  		F[i] += H[i][j]*Bc[j];	 }  }void Quad1(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;/*[ Quad1 computes the integral of several nonsingular functions along the boundary elements using a four points Gauss Quadrature.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, as defined in section 5.3 ]*/{const float Z[] =  {0.0,0.861136311594053,-0.861136311594053,0.339981043584856,-0.339981043584856};const float W[] = {0.0,0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};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 coefficients ]*/  						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 Diag1(X1,Y1,X2,Y2,G)float X1,Y1,X2,Y2,*G;{  float Ax,Ay,HL;  Ax = (X2 - X1)/2.0;  Ay = (Y2 - Y1)/2.0;  HL = sqrt(SQ(Ax) + SQ(Ay));  (*G) = 2* HL*(1.0 - log(HL)); }/*[function Solve for the solution of a linear system of equations by the Gauss Elimination method providing for interchange of rows when a zero diagonal coeficient is encountered.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[100][100],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 nonzero 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++)  /*[To interchange rows to get nonzero coeff  ]*/ 	{         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 Inter1(Bc,F,Code,Xi,Yi,X,Y,u,q1,q2)	float Bc[101],F[101],Xi[21],Yi[21],X[102],Y[102];	float u[21],q1[21],q2[21];	int Code[101];{	float temp,A,B,qx,qy,ux,uy;	int i,j,k,m;/*[A and B are 2-D arrays; Quad1 takes these arrays when called.The arrays Bc and F are rearranged to store all the values of u in Bc and all the values of q in F. ]*/	for(i=1;i<=N;i++)	{		if(Code[i] > 0)		{			temp = Bc[i];			Bc[i] = F[i];			F[i] = temp;		}	 } /*[ Compute u and q 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++)			{			  m = j + 1;Quad1(Xi[k],Yi[k],X[j],Y[j],X[m],Y[m],&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.0 * pi;			q1[k] /= 2.0 * pi;			q2[k] /= 2.0 * pi;	       }	  }  }void Out1(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,"Solution: \n");		fprintf(outfile," \nBoundary Nodes \n\n");		fprintf(outfile,"    Nodes (X,Y)			u	q\n\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\nInterior Points (Xi,Yi)  	u	  q_x	    q_y\n\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);}