sys11.c

an introduction to boundary element methods一书源码

#include "cbox11.h"

void Sys11(X, Y, Xm, Ym, G, H, Bc, F, Code, Dim) 
float X[52], Y[52], Xm[51], Ym[51];
float G[101][101], H[101][101], F[101], Bc[101];
int   Code[101], Dim;
{    
	float temp;
	int i,j,k,NN,kk,found;
	found = 0;
	
	/*[Compute coordinates of the mid-nodes  ]*/

	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;
        }
	}
   
	                
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++)
		{
			if((M-1) > 0.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)
			{
				Quad11(Xm[i],Ym[i],X[j],Y[j],X[kk],Y[kk],&H[2*i-1][2*j-1],
					&H[2*i-1][2*j],&H[2*i][2*j-1],&H[2*i][2*j],
					&G[2*i-1][2*j-1],&G[2*i-1][2*j],&G[2*i][2*j]);
				G[2*i][2*j-1] = G[2*i-1][2*j];
			}          
			else
			{
				Diag11(X[j],Y[j],X[kk],Y[kk],&G[2*i-1][2*j-1],
					&G[2*i-1][2*j],&G[2*i][2*j]);
				H[(2*i-1)][(2*j-1)] = 0.5;
				H[(2*i)][(2*j)] = 0.5;
				H[(2*i-1)][(2*j)] = 0.0;
				H[(2*i)][(2*j-1)] = 0.0;
				G[(2*i)][(2*j-1)] = G[(2*i-1)][(2*j)];
			}
		}
	}

	/*[Reorder the columns of equation as in (5.28)  
		and form system matrix A which is stored in G]*/

	NN = 2*N;  
	for(j=1;j<=NN;j++)			
	{				
		if(Code[j] > 0)     		
		{
			for(i=1;i<=NN;i++)			
			{
				temp =  G[i][j];
				G[i][j] = -H[i][j];
				H[i][j] = - temp;
			}
		}
		else
		{
			for(i=1;i<=NN;i++)			
				G[i][j] *= mu;
		}
	}

	/*[ Form the right-side vector F which is stored in F]*/

	for(i = 1;i<=NN;i++)			
	{
		F[i] = 0.0;
		for(j=1;j<=NN;j++)				
            F[i] +=  H[i][j] * Bc[j];
	}
}