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