inter11.c

an introduction to boundary element methods一书源码

#include "cbox11.h"

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