quad2.c

边界元程序

#include "cbox2.h"

void Quad2(Xp,Yp,X1,Y1,X2,Y2,A1,A2,B1,B2,
			q11,q21,q12,q22,u11,u21,u12,u22,K)
	
float  Xp,Yp,X1,Y1,X2,Y2,*A1,*A2,*B1,*B2;
float *q11,*q21,*q12,*q22,*u11,*u21,*u12,*u22;
int    K;

{
	
	/*[ This function computes the integral of several non-singular functions along the boundary elements using a four points Gauss Quadrature. 
	When K = 0, The off-diagonal coefficients of the matrices H and G are computed. When K= 1, all coefficients needed for computation of the potential and fluxes at interior points are computed. 
Ra = Radius = Distance from the collocation point to the Gauss integration points on the boundary elements; nx,ny = Components of the unit normal to the element; rx,ry,rn = Radius derivatives. See section 5.3 ]*/



float Xg[5],Yg[5];
float Ax,Ay,Bx,By,HL,nx,ny,Ra,rx,ry,rn,F1,F2,G,H,ux,uy,qx,qy;
int   i; 
float Z[] =  {0.0, 0.86113631, -0.86113631, 0.33998104, -0.33998104};
float W[] = {0.0, 0.34785485,  0.34785485, 0.65214515,  0.65214515};


	Ax = (X2 - X1)/2.0;
	Bx = (X2 + X1)/2.0;
	Ay = (Y2 - Y1)/2.0;
	By = (Y2 + Y1)/2.0;
	HL = sqrt(SQ(Ax) + SQ(Ay));
	nx =  Ay/HL;
	ny = -Ax/HL;
	(*A1) = 0.0;
	(*A2) = 0.0;
	(*B1) = 0.0;
	(*B2) = 0.0;
	(*q11) = 0.0;
	(*q21) = 0.0;
	(*q12) = 0.0;
	(*q22) = 0.0;
	(*u11) = 0.0;
	(*u21) = 0.0;
	(*u12) = 0.0;
	(*u22) = 0.0;
	
/*[ Compute the terms to be included in the matrices G and H, or the terms needed to compute the potential and fluxes at interior points ]*/

	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;
		F1 = ( 1.0 - Z[i])/2.0;
		F2 = ( 1.0 + Z[i])/2.0;
		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);
			(*q11) += qx*F1;
			(*q12) += qx*F2;
			(*q21) += qy*F1;
			(*q22) += qy*F2;
			(*u11) += ux*F1;
			(*u12) += ux*F2;
			(*u21) += uy*F1;
			(*u22) += uy*F2;
		}
		H = rn * W[i]*HL/Ra;
		G = log(1.0/Ra)*W[i]*HL;
		(*A1) -= F1*H;
		(*A2) -= F2*H;
		(*B1) += F1*G;
		(*B2) += F2*G;
	}
}
	

void Diag2(X1,Y1,X2,Y2,B1,B2)

float X1,Y1,X2,Y2,*B1,*B2; 
{
	/*[ This function computes the parts of the G matrix coefficients for  integrals along an element that includes the collocation point ]*/
	float Li;
	 Li  = sqrt(SQ(X2 -X1) + SQ(Y2 - Y1));
	(*B1) = Li*(1.5 - log(Li))/2.0;
	(*B2) = Li*(0.5 - log(Li))/2.0;
}