sys2.c

an introduction to boundary element methods一书源码

#include "cbox2.h"


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];
	}  
}