quad5.c

an introduction to boundary element methods一书源码

#include "cbox5.h"

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       ]*/

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