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