有限元程序设计作业21p
有限元程序设计作业
姓名：陆文林王子昱
学号：071053205 071054236
班级：07级研究生
专业：结构工程
院系：土木与建筑学院
有限元程序流程图
程序变量说明
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
#include
#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 { NE[i1]=new int [5]; } cout<<"输入:单元编号，起节点号，止节点号，起点是否刚结(是0,否1)，终点是否刚结(是0,否1)"; for(i1=0;i1 { for(i2=0;i2<5;i2++) { cin>>NE[i1][i2]; } } double * *NM; NM=new double * [5]; for(i1=0;i1 { NM[i1]=new double [5]; } cout<<"依次输入单元的:弹模，面积，惯性矩，长度，与x轴夹角"; for(i1=0;i1 { for(i2=0;i2<5;i2++) { cin>>NM[i1][i2]; } } nf1=3*nnode; double ** K; K=new double *[nf1]; for(i1=0;i1 { K[i1]=new double [nf1]; } for(i1=0;i1 { for(i2=0;i2 { K[i1][i2]=0.0; } } nf=6*nelem; double ** K0; K0=new double * [nf]; for(i1=0;i1 { K0[i1]=new double [nf]; } for(i1=0;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 { Po[i1]=0.0; } double * * F; F=new double *[6]; for(i1=0;i1 { F[i1]=new double [6]; } for(i1=0;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 { NL[i1]=new int[3]; } cout<<"依次输入非结点荷载信息:单元编号，非结点荷载类（均布0，三角形荷载1）,荷载在单元上的布设情况(均布0，三角形荷载:起点为零取1,起点为q取2)"; for(i1=0;i1 { or(i2=0;i2<3;i2++) { cin>>NL[i1][i2]; } } double * LD; LD=new double [nelem1]; cout<<"依次输入各单元的均布荷载值或三角形荷载的最大值:"; for(i1=0;i1 { cin>>LD[i1]; } for(i1=0;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 { Pa[i1]=0.0; } cout<<"输入有结点荷载的结点数:"; cin>>nnode1; if(nnode1!=0) { int * NLD; NLD=new int [nnode1]; cout<<"输入有结点荷载的结点序号:"; for(i1=0;i1 { cin>>NLD[i1]; } double ** LD1; LD1=new double *[3]; for(i1=0;i1 { LD1[i1]=new double [3]; } cout<<"按以上结点号顺序依次输入各结点的结点荷载的:x轴方向 力，y轴方向力，弯矩"; for(i1=0;i1 { for(i2=0;i2<3;i2++) { cin>>LD1[i1][i2]; } } for(i1=0;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 delete [] LD1[i1]; delete [] LD1; } double * P; P=new double [nf1]; for(i1=0;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 { cin>>NR[i1]; for(i1=0;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 { cin>>JJ[i1]; } for(i1=0;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"< for(i1=0;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< for(i2=0;i2<6;i2++) { cout< } cout< } return 0; } void matmat(int row,int col,double B1[6][6],double B2[6][6]) { int i1,i2,i3; for(i1=0;i1 { for(i2=0;i2 { mat3[i1][i2]=0.0; for(i3=0;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 { X[i1]=0.0; for(i2=0;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 { V[i1]=0.0; for(i2=0;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 { for(i2=i1;i2 { k1[i2]=A[i1][i2]; } irow=i1; db=fabs(A[i1][i1]); for(i2=i1+1;i2 { if(db { db=fabs(A[i2][i1]); irow=i2; if(db==0.0) { cout<<"对角线元素为零,方程组无解."; return false; } for(i2=i1;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 { A[i1][i2]/=db; } B[i1]/=db; for(i2=i1+1;i2 { db=A[i2][i1]; for(i3=i1;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 { V[i1]=B[i1]; } return * V; delete []k1; delete []B; for(i1=0;i1 delete [] A; delete [] V; } 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）,荷载在单元