sys5.c

an introduction to boundary element methods一书源码

#include "cbox5.h"


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