be5.c

an introduction to boundary element methods一书源码

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