be11.c

an introduction to boundary element methods一书源码

}



/*[function  Solve : Solution of the linear systems of equations by the Gauss elimination method providing for interchanging rows when encountering a zero diagonal coeficient.

A : System Matrix
B:  Originally it contains the Independent coefficients.  After solutions it
	contains the values of the system unknowns.                            				   				
N:  Actual number of unknowns.
Dim: Row and Column Dimension of A.
 
]*/

void Solve(A,B,D,N,Dim)

	float A[101][101],B[101],*D;  /*[ D is assigned some value here ]*/
  	int N,Dim;

{
	int N1,i,j,k,i1,l,k1,found;
	float c;
	
/*[ found is a flag which is used to check if any non zero coeff is found ]*/	
	
	
	N1 = N - 1;

	for(k=1;k<=N1;k++)
	{
		k1 = k + 1;
		c = A[k][k];
		
		if( (fabs(c) - tol) <=0 )
		{
			found = 0;
			for(j=k1;j<=N;j++)  /*[Try to Interchange rows to get Nonzero ]*/
	
	        {
            	if( (fabs(A[j][k]) - tol) > 0.0 ) 
                {
				 	for(l=k;l<=N;l++)
                    {
                     	c = A[k][l];
                        A[k][l] = A[j][l];
                        A[j][l] = c;
                     }
				  		c = B[k];
				  		B[k] = B[j];
				  		B[j] = c;
				  		c = A[k][k];
						found = 1;		/*[ coeff is found ]*/
						break;
				}
			}
		 }

		if(!found)
		{
			printf("Singularity in Row %d 1",k);
	   		(*D) = 0.0;          
	    	return;   /*[ If no coefficient is found the control is 	transferred to main ]*/
		}
	


	/*[ Divide row by diagonal coefficient ]*/		

  		c = A[k][k];
		for(j = k1;j<=N;j++)
		   A[k][j] /= c;
		B[k] /= c;  
 
/*[ Eliminate unknown X[k] from row i ]*/

		for(i = k1;i<=N; i++)
		{
		  c = A[i][k];
		  for(j = k1;j<= N;j++)
		 	A[i][j] -= c*A[k][j];
		  B[i] -= c*B[k];
		}
	}
/*[ Compute the last unknown ]*/

	if((fabs(A[N][N]) - tol) > 0.0)
	{
	   B[N] /= A[N][N];

/*[Apply back substitution to compute the remaining unknowns  ]*/

	   for(l = 1;l<= N1;l++)
	   {
	      k = N - l;
	      k1 = k +1;
	      for(j = k1;j<=N;j++)
			B[k] -= A[k][j]*B[j];
	   }
/*[Compute the value of the determinent ]*/

	(*D) = 1.0;
	for(i=1;i<=N;i++)
	  (*D) *= A[i][i];

   }
  else
  {
  		printf("Singularity in Row %d 2",k);
	   	(*D) = 0.0;
   }          

  return;	
}


void Stress( Xp,Yp,X1,Y1,X2,Y2,dx11,dy11,dx12,dy12,
             dx22,dy22,sx11,sy11,sx12,sy12,sx22,sy22)
float Xp,Yp,X1,Y1,X2,Y2,*dx11,*dy11,*dx12,*dy12;
float *dx22,*dy22,*sx11,*sy11,*sx12,*sy12,*sx22,*sy22; 
{
	float Xg[5],Yg[5],Z[5],W[5]; 	/*[dimension increased by 1]*/
	float Ax,Bx,Ay,By,nx,ny,slope,Perp,sgn,SR,FA,AL,Denom,rx,ry,Ra;
	int i;

	Z[1] = 0.86113631;
	Z[2] = -Z[1];
	Z[3] = 0.33998104;
	Z[4] = -Z[3];
	W[1] = 0.34785485;
	W[2] = W[1];
	W[3] = 0.65214515;
	W[4] = W[3];
	
	Ax = (X2-X1)/2.0;
	Bx = (X2+X1)/2.0;
	Ay = (Y2-Y1)/2.0;
	By = (Y2+Y1)/2.0;
	SR = sqrt(SQ(Ax)+SQ(Ay));
	nx = (Y2-Y1)/(2*SR);
	ny = (X1-X2)/(2*SR);
	
	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;
        
	(*dx11) = 0.0;
	(*dy11) = 0.0;
	(*dx12) = 0.0;
	(*dy12) = 0.0;
	(*dx22) = 0.0;
	(*dy22) = 0.0;
	(*sx11) = 0.0;
	(*sy11) = 0.0;
	(*sx12) = 0.0;
	(*sy12) = 0.0;
	(*sx22) = 0.0;
	(*sy22) = 0.0;

	/*[ Compute displacement and stress coefficients ]*/
	

	FA = 1.0-4.0*nu;
	AL = 1.0-2.0*nu;
	Denom = 4.0*pi*(1.0-nu); 	
  
         
	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;
		(*dx11) += (AL*rx+2*cube(rx))*W[i]*SR/(Denom*Ra);
		(*dy11) += (2*SQ(rx)*ry-AL*ry)*W[i]*SR/(Denom*Ra);
		(*dx12) += (AL*ry+2*(SQ(rx))*ry)/(Denom*Ra)*W[i]*SR;
		(*dy12) += (AL*rx+2*rx*SQ(ry))/(Denom*Ra)*W[i]*SR;
		(*dx22) += (2*rx*SQ(ry)-AL*rx)/(Denom*Ra)*W[i]*SR;
		(*dy22) += (AL*ry+2*cube(ry))/(Denom*Ra)*W[i]*SR;
		(*sx11) += (2*Perp/Ra*(AL*rx+nu*2*rx-4*cube(rx))+
                       4*nu*nx*SQ(rx)+AL*(2*nx*SQ(rx)+2*nx)-
                       FA*nx)*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sy11) += (2*Perp/Ra*(AL*ry-4*SQ(rx)*ry)+
                       4*nu*nx*rx*ry+AL*2*ny*SQ(rx)-
                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sx12) += (2*Perp/Ra*(nu*ry-4*SQ(rx)*ry)+2*nu*
                       (nx*ry*rx+ny*SQ(rx))+AL*(2*nx*rx*
                       ry+ny))*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sy12) += (2*Perp/Ra*(nu*rx-4*rx*SQ(ry))+2*nu*
                       (nx*SQ(ry)+ny*rx*ry)+AL*(2*ny*rx*ry+
                       nx))*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sx22) += (2*Perp/Ra*(AL*rx-4*rx*SQ(ry))+4*nu*
                       ny*rx*ry+AL*2*nx*SQ(ry)-FA*nx)*2*mu/
                       (Denom*SQ(Ra))*W[i]*SR;
		(*sy22) += (2*Perp/Ra*(AL*ry+2*nu*ry-4*cube(ry))+
                       4*nu*ny*SQ(ry)+AL*(2*ny*SQ(ry)+2*ny)-
                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;

     }   


 }  


void Sys11(X, Y, Xm, Ym, G, H, Bc, F, Code, Dim) 
float X[52], Y[52], Xm[51], Ym[51];
float G[101][101], H[101][101], F[101], Bc[101];
int   Code[101], Dim;
{    
	float temp;
	int i,j,k,NN,kk,found;
	found = 0;
	
	/*[Compute coordinates of the mid-nodes  ]*/

	X[N+1] = X[1];
	Y[N+1] = Y[1];

 	for(i=1;i<=N;i++)
    {
		Xm[i] = (X[i] + X[i+1])/2.0;
		Ym[i] = (Y[i] + Y[i+1])/2.0;
	}
 					

	if( (M-1) > 0 )
	{
		Xm[Last[1]] = (X[Last[1]]+X[1])/2.0;  
		Ym[Last[1]] = (Y[Last[1]]+Y[1])/2.0;
		for(k=2;k<=M;k++)
        {
			Xm[Last[k]] = (X[Last[k]] + X[Last[k-1] + 1])/2.0;
			Ym[Last[k]] = (Y[Last[k]] + Y[Last[k-1] + 1])/2.0;
        }
	}
   
	                
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++)
		{
			if((M-1) > 0.0) 	
			{			   
				if(!(j-Last[1] )) 	
					kk=1;
				else
				{
					found = 0;
					for(k=2;k<=M;k++)
					{	
						if(!(j-Last[k])) 
                        {
							kk=Last[k-1]+1;
							found = 1;
							break;
						}
					}
					if(!found)
						kk = j+1;
				}
			}
			else
				kk=j+1; 
     
			if(i-j)
			{
				Quad11(Xm[i],Ym[i],X[j],Y[j],X[kk],Y[kk],&H[2*i-1][2*j-1],
					&H[2*i-1][2*j],&H[2*i][2*j-1],&H[2*i][2*j],
					&G[2*i-1][2*j-1],&G[2*i-1][2*j],&G[2*i][2*j]);
				G[2*i][2*j-1] = G[2*i-1][2*j];
			}          
			else
			{
				Diag11(X[j],Y[j],X[kk],Y[kk],&G[2*i-1][2*j-1],
					&G[2*i-1][2*j],&G[2*i][2*j]);
				H[(2*i-1)][(2*j-1)] = 0.5;
				H[(2*i)][(2*j)] = 0.5;
				H[(2*i-1)][(2*j)] = 0.0;
				H[(2*i)][(2*j-1)] = 0.0;
				G[(2*i)][(2*j-1)] = G[(2*i-1)][(2*j)];
			}
		}
	}

	/*[Reorder the columns of equation as in (5.28)  
		and form system matrix A which is stored in G]*/

	NN = 2*N;  
	for(j=1;j<=NN;j++)			
	{				
		if(Code[j] > 0)     		
		{
			for(i=1;i<=NN;i++)			
			{
				temp =  G[i][j];
				G[i][j] = -H[i][j];
				H[i][j] = - temp;
			}
		}
		else
		{
			for(i=1;i<=NN;i++)			
				G[i][j] *= mu;
		}
	}

	/*[ Form the right-side vector F which is stored in F]*/

	for(i = 1;i<=NN;i++)			
	{
		F[i] = 0.0;
		for(j=1;j<=NN;j++)				
            F[i] +=  H[i][j] * Bc[j];
	}
}