#pragma warning(disable:4996);
//计算平面钢架的程序
//程序开始
#include<stdio.h>
#include<math.h>
#define NE 2
#define NJ 3 //定义并输入基本参数
#define NZ 6
#define NPJ 3
#define NPF 2
#define NJ3 9
#define DD 9
#define E0 3.0000E8 // 定义并输入常数
#define A0 0.5
#define I0 4.1667E-2
#define PI 3.141592654
// 这是输入参数的初始化和定义全局变量
int jm[NE+1][3]={{0,0,0},{0,1,2},{0,3,1}};
double gc[NE+1]={0.0,5.0,5.0};
double gj[NE+1]={0.0,0.0,90.0};
double mj[NE+1]={0.0,A0,A0};
double gx[NE+1]={0.0,I0,I0};
int zc[NZ+1]={0,4,5,6,7,8,9};
double pj[NPJ+1][3]={{0.0,0.0,0.0},{0.0,6.0,1.0},{0.0,2.0,2.0},{0.0,-5.0,3.0}};
double pf[NPF+1][5]={{0,0,0,0,0},{0.0,-4.8,5.0,1.0,1.0},{0.0,-8.0,2.5,2.0,2.0}};
double kz[NJ3+1][NJ3+1],p[NJ3+1];
double pe[7],f[7],f0[7],t[7][7];
double ke[7][7],kd[7][7];
// kz[][]---整体刚度矩阵
// ke[][]---整体坐标系下的单元刚度矩阵
// kd[][]---局部坐标系下的单元刚度矩阵
// t[][] ---坐标板换矩阵
//*****函数声明*****
void jdugd(int);
void zb(int);
void gdnl(int);
void dugd(int);
//***主函数开始*****
void main(void)
{
FILE *fp;
int i,j,k,e,dh,h,ii,jj,hz,a1,b1,m,l,dl,zl,z,j0;
double cl,wy[7];
int IM,IN,jn;
//<功能:形成矩阵p>
if(NPJ>0)
{
for(i=1;i<=NPJ;i++)
{
j=(int)pj[i][2];
p[j]=pj[i][1];
}
}
if(NPF>0)
{
for(i=1;i<=NPF;i++)
{
hz=i; //可否省略此行代码,下一行直接写gdnl(i),这样也是将实参i传递给形参hz
gdnl(hz);
e=(int)pf[hz][3];
zb(e);
for(j=1;j<=6;j++)
{
pe[j]=0.0;
for(k=1;k<=6;k++)
{
pe[j]=pe[j]-t[k][j]*f0[k]; //用的是转秩矩阵 zb()函数出来的是直接矩阵 在这先让列不动行动 相当于乘的转秩
}
}
a1=jm[e][1];
b1=jm[e][2];
p[3*a1-2]=p[3*a1-2]+pe[1];
p[3*a1-1]=p[3*a1-1]+pe[2];
p[3*a1]=p[3*a1]+pe[3];
p[3*b1-2]=p[3*b1-2]+pe[4];
p[3*b1-1]=p[3*b1-1]+pe[5];
p[3*b1]=p[3*b1]+pe[6];
}
}
//*********************************************
//<功能:生成整体刚度矩阵ke[][]>
for(e=1;e<=NE;e++) //这里的总刚KE实际是已经处理过的带状矩阵KZ*了!
{
dugd(e);
for(i=1;i<=2;i++)
{
for(ii=1;ii<=3;ii++)
{
h=3*(i-1)+ii;
dh=3*(jm[e][i]-1)+ii;
for(j=1;j<=2;j++)
{
for(jj=1;jj<=3;jj++)
{
l=3*(j-1)+jj;
zl=3*(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<=NZ;i++) //按支撑循环,即约束总数(假设共一个固定端支座和一个固定铰支座,则总数仍为6?)
{
z=zc[i]; //取出支撑对应的位移数(如固定铰支座对应节点码为3,则z=7 8 9)放入变量z中
kz[z][1]=1.0; //Z行第一个元素置为1.0,注意这里的KZ已是处理后的带状矩阵了,即对应于原总刚中的Z行中对角线处的元素置一了!
for(j=2;j<=DD;j++) //这个循环是将带状矩阵中Z行的其他元素全部置为零,由于KZ与KZ*中行元素是不变的,即KZ中的r行元素在KZ*中
kz[z][j]=0.0; //仍在r行。即这个循环的本质是将KZ总刚中Z行除对角线位置元素外的所有元素置为零!
if((z!=1))
{
if(z>DD)j0=DD; //先通过一个判断语句,得出正确的列数J0,要用这个列数来实现Z行斜向上45°所有元素置零,这个操作等同于
else if(z<=DD)j0=z; //将总刚KZ中的Z列除对角线元素外的所有元素置零!
for(j=2;j<=j0;j++) //j从2到j0代表了斜向上45°各元素的列变化规律
kz[z-j+1][j]=0.0; //代表了斜向上45°各元素的行变化规律
}
p[z]=0.0; //将位移约束代码所对应的荷载列阵处的元素置为零
}
/*高斯消元法解方程组*/
///
//消元
for(k=1;k<=NJ3-1;k++) //消元轮次码K取值从1到(n-1)
{
if(NJ3>k+DD-1) //判断消元行码i的取值范围(k+1)~im的条件,因为是在带状矩阵中了,已经不用从k行依次向下消到n行了,最大到im行即可!
IM=k+DD-1; //即n>k+(d-1),取im=k+(d-1),若n<k+(d-1),则取im=n;
else if(NJ3<=k+DD-1)IM=NJ3;
IN=k+1; //这是第K次消元过程,即从k+1行开始向下进行逐行消元,但已经不用消到n行了,在带状矩阵中只需消到im行即可
for(i=IN;i<=IM;i++)
{
l=i-k+1; //确定了Kki在带状矩阵中的列的数目 即总刚中的i列换为带状矩阵中的L列
cl=kz[k][l]/kz[k][1]; //这是小系数(区分开L和数1) 是经历了三次变换后的形式表达了,即由满阵存储的高斯消元法变为半阵存储的高斯消元法,又变为与他相对应的带状矩阵形式
jn=DD-l+1;
for(j=1;j<=jn;j++) //注意,这里的j实际是书上的J,即它是带状矩阵中的循环变量!
{
m=j+i-k; //满阵时是原始的小系数CL×KZ[k][j],现在换为带状矩阵KZ[k][m]了,故k行j列应换位k行m列,m=j-k+1,而带状矩阵中,遇到j就要用J替换,又J=j-i+1 故m=j-k+1=J+i-k 此处的j即是J
kz[i][j]=kz[i][j]-cl*kz[k][m];
}
p[i]=p[i]-cl*p[k];
}
}
//******回代***********
p[NJ3]=(p[NJ3]/kz[NJ3][1]);
for(i=NJ3-1;i>=1;i--)
{
if(DD>NJ3-i+1)j0=NJ3-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]);
}
fp=fopen("D:\\RESULT.txt","wt");
fprintf(fp,"\t\t%s\n","有限元法(李景湧)第二章程序 ");
printf("\n");
printf("\t\t有限元法(李景湧)第二章程序 \n");
printf("__________________________________________________________\n");
fprintf(fp,"________________________________________________________\n");
printf("NJ= U= V= CETA= \n");
fprintf(fp,"NJ= U= V= CETA= \n");//fprintf第二个参数是格式字符串,但是也可以省略,直接写第三个参数即要输出的内容!
for(i=1;i<=NJ;i++)
{
fprintf(fp,"%-9d %-12.11f %-12.11f %-12.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);//负号表示左对齐,12代表数据宽度(不包括小数点),.11表示小数点后有11位数!
printf("%-9d %-12.11f %-12.11f %-12.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]);
}
printf("__________________________________________________________\n");
fprintf(fp,"_________________________________________________________\n");
//*根据E的值输出相应E单元的N,Q,M(A,B)的结果**
printf("E= \t N= Q= M= \n");
fprintf(fp,"E= \t N= Q= M= \n");
//计算轴力N‘,剪力Q,弯矩M*
for(e=1;e<=NE;e++)
{
jdugd(e);
zb(e);
for(i=1;i<=2;i++)
{
for(ii=1;ii<=3;ii++)
{
h=3*(i-1)+ii;
dh=3*(jm[e][i]-1)+ii;
wy[h]=p[dh];
}
}
for(i=1;i<=6;i++)
{
f[i]=0.0;
for(j=1;j<=6;j++)
{
for(k=1;k<=6;k++)
{
f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];
}
}
}
if(NPF>0)
{
for(i=1;i<=NPF;i++)
if(pf[i][3]==e)
{
hz=i;
gdnl(hz);
for(j=1;j<=6;j++)
{
f[j]=f[j]+f0[j];
}
}
}
printf("%-4d(A) %-9.5f %-9.5f %-9.5f\n",e,f[1],f[2],f[3]);
printf(" (B) %-9.5f %-9.5f %-9.5f \n",f[4],f[5],f[6]);
fprintf(fp,"%-4d(A) %-9.5f %-9.5f %-9.5f\n",e,f[1],f[2],f[3]);
fprintf(fp," (B) %-9.5f %-9.5f %-9.5f \n",f[4],f[5],f[6]);
}
return;
}
//**************主程序结束**********************************
//gdnl()函数:<功能:将非节点荷载下的杆端力计算出来存入f0[]> 疑问是 为什么ind=4即力偶荷载的时候没有判断计算?
//***********************************************************
void gdnl(int hz) //主程序中的实参hz即i赋值给此处的形参hz,实现行号变换。
{
int ind,e;
double g,c,l0,d;
g=pf[hz][1]; //取出非节点荷载数组中的每一行的第一列元素,即荷载值 放入变量g中。当然,这个数组中的第零行和第零列已被置零。
c=pf[hz][2]; //取出非节点荷载数组中的每一行的第二列元素,即作用位置 放入变量c中。
e=(int)pf[hz][3]; //由于前面将非节点荷载数组定义为double型,这里int的意思是强制返回整形,取出非节点荷载数组中的每一行的第三列元素 即作用单元 放入变量e中!
ind=(int)pf[hz][4]; //取出非节点荷载数组中的每一行的第四列元素,即荷载类型 放入变量ind 中。
l0=gc[e]; //取出非节点荷载所在单元的杆长 放入变量l0中,注意,不是所有的杆长都放入。因为数组PF中,就是按一个非节点荷载存一行。
d=l0-c; //杆长值减去作用位置 放入变量d中
if(ind==1) //求均布荷载固端反力
{
f0[1]=0.0; //均布荷载作用下,i端1号位移所对应的轴力为零
f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c