有限元程序设计002.doc

有限元程序设计作业
有
限
元
程
序
设
计
作
业

程序功能说明

本程序包含四个子程序，其主要功能分别是：
void matmat void matvec
调用结点及单元信息判断并循环形成总刚
double matvec double gs
调用结点位移及结点信息和对应单元信息，回代求出并输出单元杆端力

有限元程序流程图
输入单元数

输入结点数

根据单元数形成动态数组存储单元及结点信息以及单元截面信息

输入非结点
荷载数

有

形成动态数组存储非结点荷载信息

否

根据非结点荷载信息判断并形成列向量

输入结点荷载数

否
调用结点及单元信息判断并循环形成总刚
形成总荷载列向量
输入结点荷载信息并形成结点荷载列向量
有

形成动态数组存储结构总刚信息

引入支承条件，对总刚及荷载列向量进行处理

输入铰结点数及其编号

调用结点位移及结点信息和对应单元信息，回代求出并输出单元杆端力
解线性方程组，求出结点位移
铰结点的处理

程序变量说明

NE---------------------存储单元号，单元结点号，两端约束情况
nelem--------------------单元数
nnode---------------------结点数
NM--------------------存储单元弹模，面积，惯性矩，长度，与总体X轴夹角
nf-----------------------单元自由度
nf1------------------------结点自由度
K-----------------------结构的总刚度矩阵
K0---------------------存储各单元单刚矩阵
nelem1-----------------非结点荷载单元数
nnode1-------------------有结点荷载结点数
k-----------------------单元单刚矩阵
nnode2-----------------有约束的结点数
nnode3-----------------铰结点数
T-----------------------单元坐标转换矩阵
Po---------------------非结点荷载列向列量
F----------------------单元非结点荷载的固端力
nbg-------------------起始结点号
nend------------------终止结点号
NL--------------------单元的非结点荷载信息
LD--------------------单元非结点荷载的大小
Pa---------------------所有单元结点荷载的列向量
LD1------------------存储各单元荷载信息
P-----------------------所有荷载之和列向量
NR---------------------存储支座约束的结点号
F1----------------------存储各单元的杆端力
NLD-----------------存储有结点荷载的结点号
JJ-----------------------存储铰结的结点编号
V-----------------------返回并存储结点位移列向量
V1---------------------调用各单元的结点位移

源程序
// #include "stdafx.h"
#include <iostream.h>
#include <math.h>
#define PI 3.14159265
double mat3[6][6],X[6];
void matmat(int row,int col,double B1[6][6],double B2[6][6]);
void matvec(int row, int col,double mat1[6][6],double Y1[6]);
double matvec(int row, int col,double ** mat1,double * Y1,double * V);
double gs(int row,int col,double ** A, double * B,double * V);

int main(int argc, char* argv[])
{
double m=1e8;
int nelem,nnode,nelem1,nnode1,nnode2,nnode3,nf,nf1;
int i1,i2,i3,nbg,nend,nbg1,nend1,nbg2,nend2,nbg3,nend3,N,N1,N2,n,n1;
double L,a,E,A,I,q,T[6][6],T1[6][6],k[6][6],mat1[6][6];
double P1[6],V1[6],F1[6],X1[6];
cout<<"输入单元总数:";
cin>>nelem;
cout<<"输入结点总数:";
cin>>nnode;
int * *NE;
NE=new int * [5];
for(i1=0;i1<nelem;i1++)
{
NE[i1]=new int [5];
}
cout<<"输入:单元编号，起节点号，止节点号，起点是否刚结(是0,否1)，终点是否刚结(是0,否1)";
for(i1=0;i1<nelem;i1++)
{
for(i2=0;i2<5;i2++)
{
cin>>NE[i1][i2];
}
}
double * *NM;
NM=new double * [5];
for(i1=0;i1<nelem;i1++)
{
NM[i1]=new double [5];
}
cout<<"依次输入单元的:弹模，面积，惯性矩，长度，与x轴夹角";
for(i1=0;i1<nelem;i1++)
{
for(i2=0;i2<5;i2++)
{
cin>>NM[i1][i2];
}
}
nf1=3*nnode;
double ** K;
K=new double *[nf1];
for(i1=0;i1<nf1;i1++)
{
K[i1]=new double [nf1];
}
for(i1=0;i1<nf1;i1++)
{
for(i2=0;i2<nf1;i2++)
{
K[i1][i2]=0.0;
}
}
nf=6*nelem;
double ** K0;
K0=new double * [nf];
for(i1=0;i1<nf;i1++)
{
K0[i1]=new double [nf];
}
for(i1=0;i1<nelem;i1++)
{
N=NE[i1][0],N1=6*(N-1),nbg=NE[i1][1],nend=NE[i1][2];
nbg1=NE[i1][3],nend1=NE[i1][4];
E=NM[i1][0],A=NM[i1][1],I=NM[i1][2];
L=NM[i1][3],a=NM[i1][4]*PI/180.0;
if((nbg1==0)&&(nend1==0))
{
k[0][1]=k[0][2]=k[0][4]=k[0][5]=k[1][0]=k[1][3]=k[2][0]=k[2][3]=0; k[3][1]=k[3][2]=k[3][4]=k[3][5]=k[4][0]=k[4][3]=k[5][0]=k[5][3]=0;
k[0][0]=k[3][3]=(E*A)/L,k[0][3]=k[3][0]=-(E*A)/L;
k[1][1]=k[4][4]=(12*E*I)/pow(L,3);
k[1][4]=k[4][1]=-(12*E*I)/pow(L,3);
k[1][2]=k[2][1]=k[1][5]=k[5][1]=(6*E*I)/pow(L,2);
k[2][4]=k[4][2]=k[4][5]=k[5][4]=-(6*E*I)/pow(L,2);
k[2][2]=k[5][5]=4*E*I/L,k[2][5]=k[5][2]=2*E*I/L;
}
else if((nbg1==0)&&(nend1==1))
{
k[0][1]=k[0][2]=k[0][4]=k[0][5]=0;
k[1][0]= k[1][3]=k[1][5]=k[2][0]=k[2][3]=k[2][5]=0;
k[3][1]=k[3][2]=k[3][4]=k[3][5]=0;
k[4][0]=k[4][3]=k[4][5]=0;
k[5][0]=k[5][1]=k[5][2]=k[5][3]=k[5][4]=k[5][5]=0;
k[0][0]=k[3][3]=(E*A)/L,k[0][3]=k[3][0]=-(E*A)/L;
k[1][1]=k[4][4]=3*E*I/pow(L,3);
k[1][4]=k[4][1]=-(3*E*I)/pow(L,3);
k[1][2]=k[2][1]=3*E*I/pow(L,2),k[2][4]=k[4][2]=-3*E*I/pow(L,2);
k[2][2]=3*E*I/L;
}
else if((nbg1==1)&&(nend1==1))
{
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
k[i2][i3]=0.0;
}
}
k[0][0]=k[3][3]=E*A/L,k[0][3]=k[3][0]=-E*A/L;
}
else
{
k[0][1]=k[0][2]=k[0][4]=k[0][5]=k[1][0]=k[1][2]=k[1][3]=0;
k[2][0]=k[2][1]=k[2][2]=k[2][3]=k[2][4]=k[2][5]=0;
k[3][1]=k[3][2]=k[3][4]=k[3][5]=0;
k[4][0]=k[4][2]=k[4][3]=k[5][0]=k[5][2]=k[5][3]=0;
k[0][0]=k[3][3]=(E*A)/L,k[0][3]=k[3][0]=-(E*A)/L;
k[1][1]=k[4][4]=3*E*I/pow(L,3);
k[1][4]=k[4][1]=-(3*E*I)/pow(L,3); k[1][5]=k[5][1]=3*E*I/pow(L,2);
k[4][5]=k[5][4]=-3*E*I/pow(L,2),k[5][5]=3*E*I/L;
}
for(i2=0;i2<6;i2++)
{
N2=N1+i2;
for(i3=0;i3<6;i3++)
{
K0[N2][i3]=k[i2][i3];
}
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
T[i2][i3]=0.0;
}
}
T[2][2]=T[5][5]=1,T[0][1]=T[3][4]=-sin(a),T[1][0]=T[4][3]=sin(a);
T[0][0]=T[1][1]=T[3][3]=T[4][4]=cos(a);
matmat(6,6,T,k);
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
T1[i2][i3]=T[i3][i2];
}
}
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
T[i2][i3]=T1[i2][i3];
}
}
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
mat1[i2][i3]=mat3[i2][i3];
}
}
matmat(6,6,mat1,T);
nbg3=3*(nbg-1),nend3=3*(nend-1);
for(i2=0;i2<3;i2++)
{
nbg2=nbg3+i2;
for(i3=0;i3<3;i3++)
{
nend2=nbg3+i3;
K[nbg2][nend2]+=mat3[i2][i3];
}
}
for(i2=0;i2<3;i2++)
{
nbg2=nbg3+i2;
for(i3=3;i3<6;i3++)
{
nend2=nend3+i3-3;
K[nbg2][nend2]+=mat3[i2][i3];
}
}
for(i2=3;i2<6;i2++)
{
nbg2=nend3+i2-3;
for(i3=0;i3<3;i3++)
{
nend2=nbg3+i3;
K[nbg2][nend2]+=mat3[i2][i3];
}
}
for(i2=3;i2<6;i2++)
{
nbg2=nend3+i2-3;
for(i3=3;i3<6;i3++)
{
nend2=nend3+i3-3;
K[nbg2][nend2]+=mat3[i2][i3];
}
}
}
double * Po;
Po=new double [nf1];
for(i1=0;i1<nf1;i1++)
{
Po[i1]=0.0;
}
double * * F;
F=new double *[6];
for(i1=0;i1<nelem;i1++)
{
F[i1]=new double [6];
}
for(i1=0;i1<nelem;i1++)
{
for(i2=0;i2<6;i2++)
{
F[i1][i2]=0.0;
}
}
cout<<"输入有非结点荷载的单元数:";
cin>>nelem1;
if(nelem1!=0)
{
int * * NL;
NL=new int *[3];
for(i1=0;i1<nelem1;i1++)
{
NL[i1]=new int[3];
}
cout<<"依次输入非结点荷载信息:单元编号，非结点荷载类（均布0，三角形荷载1）,荷载在单元上的布设情况(均布0，三角形荷载:起点为零取1,起点为q取2)";
for(i1=0;i1<nelem1;i1++)
{
or(i2=0;i2<3;i2++)
{
cin>>NL[i1][i2];
}
}
double * LD;
LD=new double [nelem1];
cout<<"依次输入各单元的均布荷载值或三角形荷载的最大值:";
for(i1=0;i1<nelem1;i1++)
{
cin>>LD[i1];
}
for(i1=0;i1<nelem1;i1++)
{
N=NL[i1][0]-1,n=NL[i1][1],n1=NL[i1][2];
nbg=NE[N][1],nend=NE[N][2],nbg1=NE[N][3];
nend1=NE[N][4];
E=NM[N][0],A=NM[N][1],I=NM[N][2];
L=NM[N][3],a=NM[N][4]*PI/180.0,q=LD[i1];
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
T[i2][i3]=0.0;
}
}
T[2][2]=T[5][5]=1,T[0][1]=T[3][4]=-sin(a),T[1][0]=T[4][3]=sin(a);
T[0][0]=T[1][1]=T[3][3]=T[4][4]=cos(a);
if(n==0)
{
if((nbg1==0)&&(nend1==0))
{
P1[0]=P1[3]=0,P1[1]=P1[4]=q*L/2;
P1[2]=q*pow(L,2)/12,P1[5]=-q*pow(L,2)/12;
}
else if((nbg1==0)&&(nend1==1))
{
P1[0]=P1[3]=P1[5]=0,P1[1]=5*q*L/8;
P1[4]=3*q*L/8,P1[2]=q*pow(L,2)/8;
}
else if ((nbg1==1)&&(nend1==0))
{
P1[0]=P1[2]=P1[3]=0,P1[1]=3*q*L/8;
P1[4]=5*q*L/8,P1[5]=-q*pow(L,2)/8;
}
else
{
P1[0]=P1[2]=P1[3]=P1[5]=0,P1[1]=q*L/2,P1[4]=q*L/2;
}
}
else
{
if(n1==1)
{
if((nbg1==0)&&(nend1==0))
{
P1[0]=P1[3]=0,P1[1]=3*q*L/20,P1[4]=7*q*L/20;
P1[2]=q*pow(L,2)/30,P1[5]=-q*pow(L,2)/20;
}
else if((nbg1==0)&&(nend1==1))
{
P1[0]=P1[3]=P1[5]=0,P1[1]=9*q*L/40;
P1[4]=11*q*L/40,P1[2]=7*q*pow(L,2)/120;
}
else if ((nbg1==1)&&(nend1==0))
{
P1[0]=P1[2]=P1[3]=0,P1[1]=q*L/10;
P1[4]=2*q*L/5,P1[5]=-q*pow(L,2)/15;
}
else
{
P1[0]=P1[2]=P1[3]=P1[5]=0,P1[1]=q*L/6,P1[4]=q*L/3;
}
}
else
{
if((nbg1==0)&&(nend1==0))
{
P1[0]=P1[3]=0,P1[1]=7*q*L/20,P1[4]=3*q*L/20;
P1[2]=q*pow(L,2)/20,P1[5]=-q*pow(L,2)/30;
}
else if((nbg1==0)&&(nend1==1))
{
P1[0]=P1[3]=P1[5]=0,P1[1]=2*q*L/5;
P1[4]=q*L/10,P1[2]=q*pow(L,2)/15;
}
else if ((nbg1==1)&&(nend1==0))
{
P1[0]=P1[2]=P1[3]=0,P1[1]=11*q*L/40;
P1[4]=9*q*L/40,P1[5]=-7*q*pow(L,2)/120;
}
else
{
P1[0]=P1[2]=P1[3]=P1[5]=0,P1[1]=q*L/3,P1[4]=q*L/6;
}
}
}
for(i2=0;i2<6;i2++)
{
F[N][i2]=P1[i2];
}
matvec(6,6,T,P1);
for(i2=0;i2<6;i2++)
{
P1[i2]=-X[i2];
}
nbg3=3*(nbg-1),nend3=3*(nend-1);
for(i2=0;i2<3;i2++)
{
nbg2=nbg3+i2;
Po[nbg2]+=P1[i2];
}
for(i2=3;i2<6;i2++)
{
nbg2=nend3+i2-3;
Po[nbg2]+=P1[i2];
}
}
}
double * Pa;
Pa=new double [nf1];
for(i1=0;i1<nf1;i1++)
{
Pa[i1]=0.0;
}
cout<<"输入有结点荷载的结点数:";
cin>>nnode1;
if(nnode1!=0)
{
int * NLD;
NLD=new int [nnode1];
cout<<"输入有结点荷载的结点序号:";
for(i1=0;i1<nnode1;i1++)
{
cin>>NLD[i1];
}
double ** LD1;
LD1=new double *[3];
for(i1=0;i1<nnode1;i1++)
{
LD1[i1]=new double [3];
}
cout<<"按以上结点号顺序依次输入各结点的结点荷载的:x轴方向力，y轴方向力，弯矩";
for(i1=0;i1<nnode1;i1++)
{
for(i2=0;i2<3;i2++)
{
cin>>LD1[i1][i2];
}
}
for(i1=0;i1<nnode1;i1++)
{
N=NLD[i1]-1,nbg3=3*N;
for(i2=0;i2<3;i2++)
{
nbg2=nbg3+i2;
Pa[nbg2]=LD1[i1][i2];
}
}
delete []NLD;
for(i1=0;i1<nnode1;i1++)
delete [] LD1[i1];
delete [] LD1;
}
double * P;
P=new double [nf1];
for(i1=0;i1<nf1;i1++)
{
P[i1]=Po[i1]+Pa[i1];
}
delete [] Po;
delete [] Pa;
cout<<"输入支座约束的结点数:";
cin>>nnode2;
int * NR;
NR=new int [nnode2];
cout<<"输入有支座约束的结点号:";
for(i1=0;i1<nnode2;i1++)
{
cin>>NR[i1];
for(i1=0;i1<nnode2;i1++)
{
N=NR[i1]-1;
for(i2=0;i2<3;i2++)
{
nbg=3*N+i2;
K[nbg][nbg]=m*K[nbg][nbg],P[nbg]=0;
}
}
cout<<"输入铰结点数:";
cin>>nnode3;
if(nnode3!=0)
{
int * JJ;
JJ=new int [nnode3];
cout<<"输入铰结点的结点号:";
for(i1=0;i1<nnode3;i1++)
{
cin>>JJ[i1];
}
for(i1=0;i1<nnode3;i1++)
{
nbg=3*JJ[i1]-1;
K[nbg][nbg]=m,P[nbg]=0;
}
delete [] JJ;
}
double * V;
V=new double [nf1];
gs(nf1,nf1,K,P,V);
delete [] NR;
cout<<"单元号"<<" "<<"IN"<<" "<<"IQ"<<"
"<<"IM"<<" "<<"JN"<<" "<<"JQ"<<" "<<"JM"<<endl;
for(i1=0;i1<nelem;i1++)
{
N=NE[i1][0],N1=6*(N-1),nbg=NE[i1][1];
nend=NE[i1][2],a=NM[i1][4]*PI/180.0;
nbg1=3*(nbg-1),nend1=3*(nend-1);
for(i2=0;i2<6;i2++)
{
for(i3=0;i3<6;i3++)
{
T[i2][i3]=0.0;
}
}
T[2][2]=T[5][5]=1,T[0][1]=T[3][4]=sin(a),T[1][0]=T[4][3]=-sin(a);
T[0][0]=T[1][1]=T[3][3]=T[4][4]=cos(a);
for(i2=0;i2<3;i2++)
{
nbg2=nbg1+i2;
V1[i2]=V[nbg2];
}
for(i2=3;i2<6;i2++)
{
nend2=nend1+i2-3;
V1[i2]=V[nend2];
}
matvec(6,6,T,V1);
for(i2=0;i2<6;i2++)
{
nbg2=N1+i2;
for(i3=0;i3<6;i3++)
{
k[i2][i3]=K0[nbg2][i3];
}
}
for(i2=0;i2<6;i2++)
{
X1[i2]=X[i2];
}
matvec(6,6,k,X1);
for(i2=0;i2<3;i2++)
{
nbg2=N-1;
F1[i2]=F[nbg2][i2];
}
for(i2=3;i2<6;i2++)
{
nbg2=N-1;
F1[i2]=F[nbg2][i2];
}
for(i2=0;i2<6;i2++)
{
F1[i2]=X[i2]+F1[i2];
}
cout<<N<<" ";
for(i2=0;i2<6;i2++)
{
cout<<F1[i2]<<" ";
}
cout<<endl;
}
return 0;
}

void matmat(int row,int col,double B1[6][6],double B2[6][6])
{
int i1,i2,i3;
for(i1=0;i1<row;i1++)
{
for(i2=0;i2<col;i2++)
{
mat3[i1][i2]=0.0;
for(i3=0;i3<col;i3++)
mat3[i1][i2]+=B1[i1][i3]*B2[i3][i2];
}
}
}

void matvec(int row, int col,double mat1[6][6],double Y1[6])
{
int i1,i2;
for(i1=0;i1<row;i1++)
{
X[i1]=0.0;
for(i2=0;i2<col;i2++)
{
X[i1]+=mat1[i1][i2]*Y1[i2];
}
}
}
double matvec(int row, int col,double ** mat1,double * Y1,double * V)
{
int i1,i2;
for(i1=0;i1<row;i1++)
{
V[i1]=0.0;
for(i2=0;i2<col;i2++)
{
V[i1]+=mat1[i1][i2]*Y1[i2];
}
}
return * V;
}

double gs(int row,int col,double ** A, double * B,double * V)
{
int i1,i2,i3,irow;
double db,db1;
double * k1;
k1=new double [col];
for(i1=0;i1<row;i1++)
{
for(i2=i1;i2<col;i2++)
{
k1[i2]=A[i1][i2];
}
irow=i1;
db=fabs(A[i1][i1]);
for(i2=i1+1;i2<row;i2++)
{
if(db<fabs(A[i2][i1]))
{
db=fabs(A[i2][i1]);
irow=i2;
if(db==0.0)
{
cout<<"对角线元素为零,方程组无解.";
return false;
}
for(i2=i1;i2<col;i2++)
{
A[i1][i2]=A[irow][i2];
A[irow][i2]=k1[i2];
}
db=A[i1][i1];
db1=B[i1];
B[i1]=B[irow];
B[irow]=db1;
for(i2=i1;i2<col;i2++)
{
A[i1][i2]/=db;
}
B[i1]/=db;
for(i2=i1+1;i2<row;i2++)
{
db=A[i2][i1];
for(i3=i1;i3<col;i3++)
{
A[i2][i3]-=A[i1][i3]*db;
}
B[i2]-=B[i1]*db;
}
}
for(i1=row-2;i1>=0;i1--)
{
for(i2=col-1;i2>i1;i2--)
{
B[i1]-=B[i2]*A[i1][i2];
}
}
for(i1=0;i1<col;i1++)
{
V[i1]=B[i1];
}
return * V;
delete []k1;
delete []B;
for(i1=0;i1<row;i1++)
delete [] A[i1];
delete [] A;
delete [] V;
}

输出结果

输入单元总数:5
输入结点总数:6
输入:单元编号，起节点号，止节点号，起点是否刚结(是0, 否1)，终点是否刚结(是0,否1)
1 1 2 1 0
2 3 2 0 0
3 2 4 0 0
4 5 4 0 0
5 4 6 0 1
依次输入单元的:弹模，面积，惯性矩，长度，与x轴夹角
1 1 0.08333 3 0
1 1 0.08333 3 90
1 1 0.08333 3 0
1 1 0.08333 3 90
1 1 0.08333 3 0
输入有非结点荷载的单元数:3
依次输入非结点荷载信息:单元编号，非结点荷载类型（均布0，三角形荷载1）,荷载在单元
上的布设情况(均布0，三角形荷载:起点为零取1,起点为q取2)
1 0 0
3 0 0
5 0 0
依次输入各单元的均布荷载值或三角形荷载的最大值:1 1 1
输入有结点荷载的结点数:0
输入支座约束的结点数:4
输入有支座约束的结点号:1 3 5 6
输入铰结点数:2
输入铰结点的结点号:1 6
单元号 IN IQ IM JN JQ JM
1 0.0073948 1.2274 0 -0.0073948 1.7726 -0.817786
2 3.2726 0.0221844 0.0217736 -3.2726 -0.0221844 0.0447796
3 -0.0147896 1.5 0.773006 0.0147896 1.5 -0.773006
4 3.2726 -0.0221844 -0.0217736 -3.2726 0.0221844 -0.0447796
5 0.00739481 1.7726 0.817786 -0.00739481 1.2274 0

结束语：
学有所疏，请屈老师多多指教！同时感谢屈老师这一年来对我们认真的授课与无私的教诲，在此深情的祝福：
屈老师身体健康，工作顺利！
第 20 页