quad11.c

an introduction to boundary element methods一书源码

#include "cbox11.h"


 void  Quad11(Xp,Yp,X1,Y1,X2,Y2,H11,H12,H21,H22,G11,G12,G22)
	float Xp,Yp,X1,Y1,X2,Y2,*H11,*H12,*H21,*H22,*G11,*G12,*G22;
	
{
	float Ax,Ay,Bx,By,nx,ny,sgn,Denom,Ra,rx,ry,slope,Perp;
	float Z[] =  {0.0, 0.86113631, -0.86113631, 0.33998104, -0.33998104};
	float W[] = {0.0, 0.34785485,  0.34785485, 0.65214515,  0.65214515};
	float Xg[5],Yg[5],HL;
	int i;

	Ax = (X2-X1)/2.0;
	Bx = (X2+X1)/2.0;
	Ay = (Y2-Y1)/2.0;
	By = (Y2+Y1)/2.0;
	nx = (Y2-Y1)/(2*sqrt(SQ(Ax)+SQ(Ay)));
	ny = (X1-X2)/(2*sqrt(SQ(Ax)+SQ(Ay)));
	if(Ax)
	{
		slope = Ay/Ax;
		Perp = fabs(( slope*Xp-Yp+Y1-slope*X1)/sqrt(SQ(slope)+1));
	}
	else
		Perp = fabs(Xp-X1);
       
	/*[ Determine the direction of the outward normal ]*/

	sgn = (X1-Xp)*(Y2-Yp)-(X2-Xp)*(Y1-Yp);
	if (sgn < 0) 
		Perp = -Perp;

	(*H11) = 0.0;
	(*H12) = 0.0;
	(*H21) = 0.0;
	(*H22) = 0.0;
	(*G11) = 0.0;
	(*G12) = 0.0;
	(*G22) = 0.0;

	/*[ Compute coefficients of the matrices G and H ]*/       

	Denom = 4*pi*(1-nu);
	HL = sqrt(SQ(Ax)+SQ(Ay));
	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;

        (*G11) += ((3-4*nu)*log(1.0/Ra)+ SQ(rx))*W[i]*
                            	HL/(2*Denom*mu);
        (*G12)   += rx*ry*W[i]*HL/(2*Denom*mu);
         
	    (*G22)     += ((3-4*nu)*log(1.0/Ra)+SQ(ry))*W[i]*
                       			HL/(2*Denom*mu);
         
		(*H11) -= Perp*((1-2*nu)+2*SQ(rx))/(SQ(Ra)*Denom)*W[i]*
                       			HL;
        (*H12) -= (Perp*2*rx*ry/Ra+(1-2*nu)*(nx*ry-ny*rx))*
                       		W[i]*HL/(Ra*Denom);
        (*H21) -= (Perp*2*rx*ry/Ra+(1-2*nu)*(ny*rx-nx*ry))*
                       		 W[i]*HL/(Ra*Denom);
        (*H22) -= Perp*((1-2*nu)+2*SQ(ry))*W[i]*
                       			HL/(SQ(Ra)*Denom);
        
	}

       
}

void Diag11(X1,Y1,X2,Y2,G11,G12,G22)
float X1,Y1,X2,Y2,*G11,*G12,*G22;
{
	float Ax,Ay,SR,Denom;

	Ax = (X2-X1)/2;
	Ay = (Y2-Y1)/2;
	SR =sqrt(SQ(Ax)+SQ(Ay));;
	Denom = 4*pi*mu*(1-nu);
 	(*G11) = SR*((3-4*nu)*(1-log(SR))+SQ(X2-X1)/(4*SQ(SR)))/Denom;
 	(*G22) = SR*((3-4*nu)*(1-log(SR))+SQ(Y2-Y1)/(4*SQ(SR)))/Denom;
 	(*G12) = (X2-X1)*(Y2-Y1)/(4*SR*Denom);

}