be2.c

an introduction to boundary element methods一书源码

#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 Matrix
B:  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];
	}  
}