?? be11.c
字號:
#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);
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -