有限元三角形划分程序

#include
#include

#define NE=36 //单元数
#define NJ=28 //节点数
#define NZ=14 //支承数即约束节点数*2，包括U和V
#define NPJ=7 //载荷作用节点数,即载荷作用于多少个节点
#define NJ2=56 //总位移数即节点数*2
#define DD=16 //半带宽即相邻两个节点最大差值+1再整体*2

int LXM=0; //类型判别码（平面应力=0，平面应变=1）
double E0=25; //弹性模量
double MU=0.17; //泊松比
double LOU=24; //容重
double TE=1.5; //厚度

//注意：此程序数组多第一行多第一列，目的为了表示方便，自行添加
//位移表示意为：将节点位移按节点顺序排列，某节点的U或V位移对应于第几个位移

double AJZ[NJ+1][3]=; //节点坐标数组

int JM[NE+1][4]=; //节点编码即每个单元i,j,m顺序节点数组

int NZC[NZ+1]=; //约束节点的位移表示数组
double PJ[NPJ+1][2+1]=; //节点荷载数值大小及荷载对应位移表示数组

double AE,KZ[NJ2+1][DD+1],P[NJ2+1],S[3+1][6+1],KE[6+1][6+1];
int IE,JE,ME;

void DUGD(int,int);

void main()
{
int NJ1,k,IN,IM,jn,m,i,j,z,J0,ii,jj,h,dh,E,l,zl,dl;
double PE,c,SIG1,SIG2,SIG3,PYL,RYL,MAYL,MIYL,CETA;
double WY[6+1],YL[3+1];

if(LXM!=0)
{
E0=E0/(1.0-MU*MU);
MU=MU/(1.0-MU);
}

for(i=0;i<=NJ2;i++)
{
for(j=0;j<=DD;j++)
KZ[i][j]=0.0;
}

for(E=1;E<=NE;E++)
{
DUGD(E,3);
for(i=1;i<=3;i++)
{
for(ii=1;ii<=2;ii++)
{
h=2*(i-1)+ii;
dh=2*(JM[E][i]-1)+ii;
for(j=1;j<=3;j++)
{
for(jj=1;jj<=2;jj++)
{
l=2*(j-1)+jj;
zl=2*(JM[E][j]-1)+jj;
dl=zl-dh+1;
if(dl>0)
KZ[dh][dl]=KZ[dh][dl]+KE[h][l];
}
}
}
}
}

for(i=1;i<=NJ2;i++)
P[i]=0.0;

if(NPJ>0)
{
for(i=1;i<=NPJ;i++)
{
j=(int)PJ[i][2];
P[j]=PJ[i][1];
}
}

if(LOU>0)
{
for(E=1;E<=NE;E++)
{
DUGD(E,1);
PE=-LOU*(AE)*TE/3;
P[2*IE]=P[2*IE]+PE;
P[2*JE]=P[2*JE]+PE;
P[2*ME]=P[2*ME]+PE;
}
}

for(i=1;i<=NZ;i++)
{
z=NZC[i];
KZ[z][1]=1.0;
for(j=2;j<=DD;j++)
KZ[z][j]=0.0;
if(z!=1)
{
if(z>DD)
J0=DD;
else
J0=z;
for(j=2;j<=J0;j++)
KZ[z-j+1][j]=0.0;
}
P[z]=0.0;
}

NJ1=NJ2-1;
for(k=1;k<=NJ1;k++)
{
if(NJ2>k+DD-1)
IM=k+DD-1;
else
IM=NJ2;

IN=k+1;
for(i=IN;i<=IM;i++)
{
l=i-k+1;
c=KZ[k][l]/KZ[k][1];
jn=DD-l+1;
for(j=1;j<=jn;j++)
{
m=j+i-k;
KZ[i][j]=KZ[i][j]-c*KZ[k][m];
}
P[i]=P[i]-c*P[k];
}
}

P[NJ2]=P[NJ2]/KZ[NJ2][1];
for(i=NJ1;i>=1;i--)
{
if(DD>NJ2-i+1)
J0=NJ2-i+1;
else
J0=DD;
for(j=2;j<=J0;j++)
{
h=j+i-1;
P[i]=P[i]-KZ[i][j]*P[h];
}
P[i]=P[i]/KZ[i][1];
}
printf("\n");
printf("JD U

V\n");
for(i=1;i<=NJ;i++)
{
printf("%2d %-21.6f %-21.6f\n",i,P[2*i-1],P[2*i]);
}

for(E=1;E<=NE;E++)
{
DUGD(E,2);

for(i=1;i<=3;i++)
{
for(j=1;j<=2;j++)
{
h=2*(i-1)+j;
dh=2*(JM[E][i]-1)+j;
WY[h]=P[dh];
}
}
for(i=1;i<=3;i++)
{
YL[i]=0;
for(j=1;j<=6;j++)
YL[i]=YL[i]+S[i][j]*WY[j];
}
SIG1=YL[1];
SIG2=YL[2];
SIG3=YL[3];
PYL=(SIG1+SIG2)/2;
RYL=sqrt(pow((SIG1-SIG2)/2.0,2)+pow(SIG3,2));
MAYL=PYL+RYL;
MIYL=PYL-RYL;
if(SIG2==MIYL)
CETA=0;
else
CETA=90-57.29578*atan2(SIG3,(SIG2-MIYL));
printf("\n");
printf("E=%2d\n",E);
printf("sx=%-21.6f sy=%-21.6f tou=%-21.6f\n",SIG1,SIG2,SIG3);
printf("s1=%-21.6f s2=%-21.6f theta=%-21.6f\n",MAYL,MIYL,CETA);
}
}
void DUGD(int E,int ASK)
{
double B[3+1][6+1],D[3+1][3+1];
int i,j,k,CM1,CM2,CM,BM,CJ,BJ;

IE=JM[E][1];
JE=JM[E][2];
ME=JM[E][3];
CM=AJZ[JE][1]-AJZ[IE][1];
BM=AJZ[IE][2]-AJZ[JE][2];
CJ=AJZ[IE][1]-AJZ[ME][1];
BJ=AJZ[ME][2]-AJZ[IE][2];
AE=(BJ*CM-BM*CJ/2.0);

if(ASK>1)
{
for(i=1;i<=3;i++)
for(j=1;j<=6;j++)
B[i][j]=0.0;

B[1][1]=(-BJ-BM)/(2.0*AE);
B[1][3]=BJ/(2*AE);
B[1][5]=BM/(2*AE);
B[2][2]=(-CJ-CM)/(2*AE);
B[2][4]=CJ/(2*AE);
B[2][6]=CM/(2*AE);
B[3][1]=B[2][2];
B[3][2]=B[1][1];
B[3][3]=B[2][4];
B[3][4]=B[1][3];
B[3][5]=B[2][6];
B[3][6]=B[1][5];
D[1][1]=E0/(1-MU*MU);
D[1][2]=E0*MU/(1-MU*MU);
D[1][3]=0;
D[2][1]=D[1][2];
D[2][2]=D[1][1];
D[2][3]=0;
D[3][1]=0;
D[3][2]=0;
D[3][3]=E0/(2*(1+MU));

for(i=1;i<=3;i++)
{
for(j=1;j<=6;j++)
{
S[i][j]=0.0;
for(k=1;k<=3;k++)
S[i][j]=S[i][j]+D[i][k]*B[k][j];
}
}
if(ASK>2)
{
for(i=1;i<=6;i++)
{
for(j=1;j<=6;j++)
{
KE[i][j]=0.0;
for(k=1;k<=3;k++)
KE[i][j]=KE[i][j]+S[k][i]*B[k][j]*AE*TE;
}
}

}
}
}