be11.c

an introduction to boundary element methods一书源码

#include "cbox11.h"

main()
{
	/*[The lengths of arrays has been increased by one, first 
					one is not used ]*/
	int Code[101],Dim,NN;
	float X[52], Y[52], Xm[51], Ym[51],Xi[21], Yi[21],Bc[101],F[101];
	float G[101][101],H[101][101],D,stress[61],displ[41];
       
	Dim = 100;	/*[Dim = Max dimension of the system AX = F]*/

	/*[ Read Data ]*/
	
	Input11(Xi,Yi,X,Y,Code,Bc);
	
	
	/*[ Compute matrices G and H and form the system AX = F  ]*/
	
	Sys11(X,Y,Xm,Ym,G,H,Bc,F,Code,Dim);
	
	/*[Solve the system AX = F ]*/
	
	NN = 2*N;
	
	Solve(G,F,&D,NN,Dim);
	
	/*[ Compute stress and displacement at interior points ]*/
	
	Inter11(Bc,F,Code,Xi,Yi,X,Y,stress,displ);
	
	/*[Output solution at boundary nodes and interior Points ]*/
	
	Out11(Xm,Ym,Bc,F,Xi,Yi,stress,displ);
}



void Input11(Xi,Yi,X,Y,Code,Bc)
float Xi[21],Yi[21],X[52],Y[52],Bc[101];
int Code[101];
{
	int i,j,k,Dim,NN,lnsize;
	FILE *infile,*outfile;
	char Title[80],line1[100];
	lnsize = 120;
	printf("\nNOTE: FIRST LINE IN THE INPUT FILE SHOULD BE EITHER BLANK OR THE TITLE NOT DATA \n");
	printf("\nEnter the name of the input file:");
    scanf("%12s", inname);
	printf("\nEnter the name of the output file:");
	scanf("%12s", outname);
       
       
	infile = fopen(inname, "r");
	outfile = fopen(outname,"w");
	
	fgets(line1,lnsize,infile);
	fprintf(outfile,"%s \n",line1);
	fprintf(outfile,"Input \n");
	
	
	fscanf(infile,"%d %d %d", &N, &L, &M); 
	fprintf(outfile,"\nNumber of Boundary Elements = %2d \n",N);
	fprintf(outfile,"Number of Interior Points = %2d \n",L);
	
	for(i = 1; i<=5; i++) 
	   fscanf(infile,"%d",&Last[i]); 

	fscanf(infile,"%f %f", &mu, &nu);
	fprintf(outfile, "\nShear Modulus = %10.4f \n",mu);
	fprintf(outfile,"Poisson Ratio = %10.4f \n",nu);

	if(M > 0)
		{
			fprintf(outfile,"Number of different Boundaries = %d\n\n",M);
			for(i=1;i<=M;i++)
				fprintf(outfile,"Last node on  boundary%2d = %2d\n\n",i,Last[i]);
		}


	/*[ Read coordinates of extreme points ]*/

  fprintf(outfile,"\nCOORDINATES OF EXTREME POINTS OF THE BOUNDARY ELEMENTS\n");
		fprintf(outfile,"\n\nPoint      X 	 	    Y\n");
        for(i = 1; i<=N; i++)
        {
        	fscanf(infile,"%f %f", &X[i],&Y[i]);
	  		fprintf(outfile,"%2d  %10.4f \t %10.4f \n",i,X[i],Y[i]);
        }

	/*[Read boundary conditions in Bc vector. If Code[i] = 0, the Bc[] 	value is a prescribed displacement; if Code[] = 1, the Bc[] value is a prescribed traction. ]*/ 

	fprintf(outfile,"\nBoundary Conditions\n\n");
	fprintf(outfile,"	  Prescribed Value			  Prescribed Value\n");
	fprintf(outfile,"Node   X-direction		Code	   Y-direction		Code\n");
	
	for( i = 1; i<= N; i++)
	{ 
		fscanf(infile,"%d %f %d %f",&Code[2*i-1],&Bc[2*i-1],										
				&Code[2*i],&Bc[2*i]);   
	 
	   fprintf(outfile,"%2d \t %10.4f \t\t %2d \t %10.4f \t\t %2d\n",i,
	   Bc[2*i-1],Code[2*i-1],Bc[2*i-1],Code[2*i]);										 
	 
	}

	if( L)
	{
		fprintf(outfile,"\nInterior Point Coordinates\n");
		fprintf(outfile,"\nPoint 	   Xi 	  	 	   Yi\n"); 
		for(i=1;i<=L;i++)
		{
			fscanf(infile,"%f %f",&Xi[i],&Yi[i]);
			fprintf(outfile,"%2d \t %10.4f \t %10.4f\n",i,Xi[i],Yi[i]);
		    }
		  }
  
  		for(i=1;i<=80;i++)
			fprintf(outfile,"*");
		fprintf(outfile,"\n");
  		fclose(infile);
		fclose(outfile);
}


void Inter11( Bc,F,Code,Xi,Yi,X,Y,stress,displ)
float Bc[101],F[101],Xi[21],Yi[21],X[52],Y[52],stress[61],displ[41];
int Code[101]; 
{
	int NN,i,j,k,kk,lk,found;
	float temp,dx11,dy11,dx12,dy12,dx22,dy22,sx11,sy11,sx12,sy12,sx22,sy22;
    float H11,H12,H21,H22,G11,G12,G22;

	found = 0;

	/*[ Enter all displacements in  Bc and all tractions in F ]*/

	NN = 2*N;
	for(i=1;i<=NN;i++)				
	{		
		if( Code[i] > 0 ) 
		{
			temp = Bc[i];
            Bc[i] = F[i];
            F[i] = temp;
		}
		else
			F[i] *= mu;
	}
     
	/*[Compute stress and displacement at interior points]*/

	if(L)
	{
		for(k=1;k<=L;k++)       
		{
			displ[2*k-1] = 0.0;
			displ[2*k] = 0.0;
			stress[3*k-2] = 0.0;
			stress[3*k-1] = 0.0;
			stress[3*k] = 0.0;

			for(j=1;j<=N;j++)				
			{
				if((M-1) > 0)           
				{
					if (!(j-Last[1]))   
					kk = 1;
					else
					{
			  			found = 0;
						for(lk=2;lk<=M;lk++)		   
			    		{
			     			if(!(j-Last[lk])) 
							{
			       				kk = Last[lk-1]+1;
								found = 1;
								break;
							}
			    		}
						if(!found)
							kk = j+1;
			  		}
				}
				else
				kk = j+1;
				Quad11(Xi[k],Yi[k],X[j],Y[j],X[kk],Y[kk],&H11,&H12,&H21,&H22,
						&G11,&G12,&G22);

          		displ[2*k-1]+=F[2*j-1]*G11+F[2*j]*G12-
							Bc[2*j-1]*H11-Bc[2*j]*H12;

          		displ[2*k] += F[2*j-1]*G12+F[2*j]*G22-
                           		Bc[2*j-1]*H21-Bc[2*j]*H22;
		 		Stress(Xi[k],Yi[k],X[j],Y[j],X[kk],Y[kk],&dx11,&dy11,
							&dx12,&dy12,&dx22,&dy22,&sx11,&sy11,
							&sx12,&sy12,&sx22,&sy22);
				stress[3*k-2] += F[2*j-1]*dx11+F[2*j]*dy11-
	                          Bc[2*j-1]*sx11-Bc[2*j]*sy11;
				stress[3*k-1] += F[2*j-1]*dx12+F[2*j]*dy12-
	                          	Bc[2*j-1]*sx12-Bc[2*j]*sy12;
				stress[3*k] += F[2*j-1]*dx22+F[2*j]*dy22-
	                          Bc[2*j-1]*sx22-Bc[2*j]*sy22;
			}
		}
	}
}


void Out11(Xm,Ym,Bc,F,Xi,Yi,stress,displ)
float Xm[51],Ym[51],Bc[101],F[101],Xi[21],Yi[21],stress[61],displ[41];
{
	
	FILE *outfile;
	int i,j,k;
	
	outfile = fopen(outname,"a");
		
	fprintf(outfile,"\nResults:\n\n Boundary Nodes\n\n");
	fprintf(outfile,"	   X		  Y		     Displ X	Displ Y  Traction X  Traction Y\n");
	for(i=1;i<=N;i++)
		fprintf(outfile,"(%10.4f, %10.4f) %10.4f %10.4f %10.4f %10.4f\n",
						Xm[i],Ym[i],Bc[2*i-1],Bc[2*i],F[2*i-1],F[2*i]);
		
	if(L)
	{	
		fprintf(outfile,"\nInterior point displacements\n\n");
		fprintf(outfile,"	   Xi		  Yi		 Displacement X		Displacement Y\n");
		for(k=1;k<=L;k++)
			fprintf(outfile,"(%10.4f,%10.4f) \t %10.4f \t %10.4f\n",
				Xi[k],Yi[k],displ[2*k-1],displ[2*k]); 
		
		fprintf(outfile,"\nInterior point stresses\n\n");
		fprintf(outfile,"	   Xi         Yi        Sigma X           Tau XY       Sigma Y\n");
	
		for(k=1;k<=L;k++)
			fprintf(outfile,"(%10.4f,%10.4f)\t %10.4f \t  %10.4f \t%10.4f\n",
					Xi[k],Yi[k],stress[3*k-2],stress[3*k-1],stress[3*k]);
					
	}
	fclose(outfile);
}



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