sys1.c

an introduction to boundary element methods一书源码

#include "cbox1.h"
void Sys1(X,Y,XM,YM,G,H,Bc,F,Code,Dim)
  float X[102],Y[102],XM[101],YM[101],Bc[101],F[101],G[100][100],H[100][100];
  int Code[101],Dim;
{
	float temp,DQ1,DQ2,DU1,DU2;
	int  i,j,l,m; 
  /*[ Compute coordinates of nodes and store in arrays XM and YM ]*/
	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;
    }
/*[ Compute the coeficients of matrices G and H  ]*/
	for(i=1;i<=N;i++)
	{
	 	for(j=1;j<=N;j++)
	  	{
	    	m = j+1;
			if((i-j))
		 		Quad1(XM[i],YM[i],X[j],Y[j],X[m],Y[m],&H[i][j],&G[i][j],&DQ1,&DQ2,&DU1,&DU2,0);
		 else
		 {
		  Diag1(X[j],Y[j],X[m],Y[m],&G[i][j]);
		  H[i][j] = pi;
		  	}
		}
	  }

/*[ Reorder the columns of the system AX = F in accordance with the boundary conditions and form system matrix A which is stored in G ]*/
	for(i=1;i<=N;i++)
	{
	 	if(Code[i] > 0)
	  	{
	    	for(j=1;j<=N;j++)
			{
		  	temp = G[j][i];
		  	G[j][i] = -H[j][i];
		  	H[j][i] = -temp;
		 	}
	  	} 
	 }
/*[ Form the right-side Vector F which is Stored in F ]*/
	for(i=1;i<=N;i++)
	{
	 	F[i] = 0;
	 	for(j=1;j<=N;j++)
	  		F[i] += H[i][j]*Bc[j];
	 }
  }