quad1.c

边界元程序

#include "cbox1.h"
void Quad1(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;
/*[ Quad1 computes the integral of several nonsingular functions along the boundary elements using a four points Gauss Quadrature.
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, as defined in section 5.3 ]*/
{
const float Z[] =  {0.0,0.861136311594053,-0.861136311594053,0.339981043584856,-0.339981043584856};const float W[] = 
{0.0,0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};
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 coefficients ]*/  					
	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 Diag1(X1,Y1,X2,Y2,G)
float X1,Y1,X2,Y2,*G;
{
  float Ax,Ay,HL;
  Ax = (X2 - X1)/2.0;
  Ay = (Y2 - Y1)/2.0;
  HL = sqrt(SQ(Ax) + SQ(Ay));
  (*G) = 2* HL*(1.0 - log(HL));
 }