平面刚架静力分析主程序

1、平面刚架静力分析主程序/     STATIC ANALYSIS OF PLANE FRAMES/#include <iostream.h>#include <fstream.h>#include <stdlib.h>#include <math.h>#include <iomanip.h>/ Functions declaration ( for C+ Only )void INPUT(double X, double Y, int NCO, double PROP,   

2、                    double AL, int IB, double REAC);void ASSEM(double X, double Y, int NCO, double PROP,                  

3、60;    double TK30, double ELST6, double AL);void STIFF(int NEL, double X, double Y, double PROP,                       int NCO, double ELST6, double AL);void ELASS(int NEL, int NCO, double

4、TM30, double ELMAT6);void BOUND(double TK30, double AL, double REAC, int IB);void SLBSI(double A30, double B, double D, int N, int MS,                      int NRMX, int NCMX);void FORCE(int NCO, doub

5、le PROP, double FORC, double REAC,                      double X, double Y, double AL);void OUTPT(int NCO, double AL, double FORC, double REAC);void BTAB3(double A6, double B6, double V, int N, int NX

6、);/  INITIALIZATION OF GLOBAL VARIABLESint NN,NE,NLN,NBN,N,MS;double E,G;/ ASSIGN DATA SET NUMBERS TO IN,FOR INPUT,AND IO FOR OUTPUTifstream READ_IN;ofstream WRITE_IO;/ INITIALIZATION OF PROGRAM PARAMETERSint NRMX=300;int NCMX=30;int NDF=3;int NNE=2;int NDFEL=NDF*NNE; /    

7、/     MAIN PROGRAM/int main()     double X100,Y100,PROP200,TK30030,            AL300,FORC600,REAC300,ELST66,V30;          int NCO200, IB80;/ OPEN ALL FILES 

8、;    READ_IN.open("DATAI.TXT");     WRITE_IO.open("DATAO.TXT");/ DATA INPUT     INPUT(X,Y,NCO,PROP,AL,IB,REAC);/ ASSEMBLING OF TOTAL STIFFNESS MATRIX     ASSEM(X,Y,NCO,PROP,TK,ELST,AL);/ INTRODUCTION OF BO

9、UNDARY CONDITIONS     BOUND(TK,AL,REAC,IB);/      / SOLUTION OF THE SYSTEM OF EQUATIONS     SLBSI(TK,AL,V,N,MS,NRMX,NCMX);/     / COMPUTATION OF MEMBER FORCES     FORCE(NCO,PROP,FORC,REAC,X,Y,AL)

10、;/ OUTPUT     OUTPT(NCO,AL,FORC,REAC);/ CLOSE ALL FILES          READ_IN.close();     WRITE_IO.close();          return 0;子程序INPUT void INPUT(double X, double Y, int NCO,

11、 double PROP,                       double AL, int IB, double REAC)/ INPUT PROGRAM/          int I, NUM, N1, IC3, K, L, L1, L2, N2;    

12、      double W3;          WRITE_IO.setf(ios:fixed);          WRITE_IO.setf(ios:showpoint);                

13、0;   WRITE_IO << " " <<      "*"                    << endl;/ READ BASIC PARAMETERS     READ_IN >> NN >> NE >>

14、; NLN >> NBN >> E;     WRITE_IO << "nn INTERNAL DATA nn" << " NUMBER OF NODES         :"                   

15、;       << setw(5) << NN << "n" << " NUMBER OF ELEMENTS      :" << setw(5)                      

16、;         << NE << "n" << " NUMBER OF LOADED NODES  :" << setw(5) << NLN                        

17、;       << "n" << " NUMBER OF SUPPORT NODES :" << setw(5) << NBN << "n"                         

18、;      << " MODULUS OF ELASTICITY :" << setw(15) << setprecision(0) << E                               <

19、;< "nn" << " NODAL COORDINATESn" << setw(11) << "NODE" << setw(7)                               <<

20、; "X" << setw(10) << "Yn"/ READ NODAL COORDINATES IN ARRAY X AND Y          for (I=1; I<=NN; I+)                     

21、60;        READ_IN >> NUM >> XI-1 >> YI-1;                    WRITE_IO.precision(2);         WRITE_IO << s

22、etw(10) << NUM << setw(10) << XI-1 << setw(10) << YI-1                                   << "n"

23、;            / READ ELEMENT CONNECTIVITY IN ARRAY NCO AND/ ELEMENT PROPERTIES IN ARRAY PROP          WRITE_IO << "n ELEMENT CONNECTIVITY AND PROPERTIESn" << setw(11)    

24、;                      << "ELEMENT" << setw(23) << "START NODE  END NODE" << setw(9)            

25、                   << "AREA" << endl;          for (I=1; I<=NE; I+)              

26、;                N1=NNE*(I-1);                    READ_IN >> NUM >> IC0 >> IC1 >> W0 >> W1; 

27、60;                  WRITE_IO.precision(5);                    WRITE_IO << setw(10) << NUM << setw(10)

28、 << IC0 << setw(10) << IC1                                   << setw(15) << W0 << setw(15) << W1 &

29、lt;< "n"                    PROPN1=W0;                    PROPN1+1=W1;    

30、60;    NCON1=IC0;                    NCON1+1=IC1;               / COMPUTE ACTUAL NUMBER OF UNKNOWNS AND CLEAR THE LOAD VECTOR

31、     N=NN*NDF;          for (I=1; I<=N; I+)                              ALI-1=0.0; 

32、60;        / READ THE NODAL LOADS AND STORE THEM IN ARRAY AL     WRITE_IO << "n NODAL LOADSn" << setw(11) << "NODE" << setw(7) << "PX"         &#

33、160;                << setw(10) << "PY" << setw(10) << "MZ" << endl;          for (I=1; I<=NLN; I+)     

34、60;                        READ_IN >> NUM >> W0 >> W1 >> W2;         WRITE_IO.precision(2);      

35、;   WRITE_IO << setw(10) << NUM << setw(10) << W0 << setw(10) << W1                                 

36、;  << setw(10) << W2 << "n"              for (K=1; K<=NDF; K+)                        

37、60;                        L=NDF*(NUM-1)+K;                       

38、0;     ALL-1=WK-1;                                   / READ BOUNDARY NODES DATA. STORE UNKNOWN STATUS INDICATORS/ IN ARR

39、AY IB, AND PRESCRIBED UNKNOWN VALUES IN ARRAY REAC     WRITE_IO << "n BOUNDARY CONDITION DATAn" << setw(29) << "STATUS"                      &

40、#160;   << setw(31) << "PRESCRIBED VALUESn" << setw(37)                               << "(0:PRESCRIBED, 1:FREE

41、)n" << setw(11) << "NODE" << setw(9)                               << "U" << setw(10) << "V&qu

42、ot; << setw(10) << "RZ" << setw(17) << "U"                               << setw(10) << "V" &l

43、t;< setw(10) << "RZ" << endl;      for (I=1; I<=NBN; I+)                              READ_IN >> NUM >&

44、gt; IC0 >> IC1 >> IC2 >> W0 >> W1 >> W2;                    WRITE_IO.precision(4);               

45、60;    WRITE_IO << setw(10) << NUM << setw(10) << IC0 << setw(10) << IC1                               

46、<< setw(10) << IC2 << setw(20) << W0 << setw(10) << W1                                << setw(10) << W2 <

47、;< "n"              L1=(NDF+1)*(I-1)+1;                    L2=NDF*(NUM-1);         

48、0;          IBL1-1=NUM;                    for (K=1; K<=NDF; K+)               

49、60;                                 N1=L1+K;                 

50、;            N2=L2+K;                             IBN1-1=ICK-1;       

51、                      REACN2-1=WK-1;                           &#

52、160;  /          return;子程序ASSEM /void ASSEM(double X, double Y, int NCO, double PROP,                       double TK20, double ELST5, double AL)/ AS

53、SEMBLING OF THE TOTAL MATRIX FOR THE PROBLEM/    int N1, I, L1, J, L2, J1, K, L3, L, NEL;  / COMPUTE HALF BAND WIDTH AND STORE IN MS    N1=NNE-1;         MS=0;         for (I=1; I<=N

54、E; I+)                      L1=NNE*(I-1);             for (J=1; J<=N1; J+)            &

55、#160;                       L2=L1+J;                 J1=J+1;           

56、0;     for (K=J1; K<=NNE; K+)                                           

57、      L3=L1+K;                     L=abs(NCOL2-NCOL3);                     if (MS-L)&

58、lt;=0)                                                  

59、;       MS=L;                                           

60、0;                                                    

61、MS=NDF*(MS+1);/ CLEAR THE TOTAL STIFFNESS MATRIX             for (I=1; I<=N; I+)                      for (J=1; J<=MS; J+)   

62、0;                           TKIJ=0.0;                for (I=1; I<=5; I+)     

63、0;                for (J=1; J<=5; J+)                               ELSTI-1J-1=0.0;

64、0;                                    /          for (NEL=1; NEL<=NE; NEL+)   

65、;      /      COMPUTE THE STIFFNESS MATRIX FOR ELEMENT NEL        STIFF(NEL,X,Y,NCO,PROP,ELST,AL);/      PLACE THE MATRIX IN THE TOTAL STIFFNESS MATRIX      

66、60; ELASS(NEL,NCO,TK,ELST);                      /         return;子程序STIFF /void STIFF(int NEL, double X, double Y, double PROP,      &

67、#160;                int NCO, double ELST6, double AL)/ COMPUTATION OF ELEMENT STIFFNESS MATRIX FOR CURRENT ELEMENT/    int L, N1, N2, I, J;     double DX, DY, D, CO, SI, AX, YZ, ROT66, V30;  /&#

68、160;   L=NNE*(NEL-1);         N1=NCOL;         N2=NCOL+1;         AX=PROPL;         YZ=PROPL+1;/ COMPUTE LENGTH OF ELEMENT, AND SI

69、NE AND COSINE OF ITS LOCAL X AXIS         DX=XN2-1-XN1-1;         DY=YN2-1-YN1-1;         D=sqrt(DX*DX+DY*DY);         CO=DX/D;  &

70、#160;      SI=DY/D;/ CLEAR THE ELEMENT STIFFNESS AND ROTATION MATRICES    for (I=1;I<=6;I+)                      for (J=1;J<=6;J+)     

71、0;                     ELSTI-1J-1=0.0;          ROTI-1J-1=0.0;                   

72、0; / FORM ELEMENT ROTATION MATRIX    ROT00=CO;         ROT01=SI;         ROT10=-SI;         ROT11=CO;         ROT22=1.0; 

73、;        for (I=1;I<=3;I+)                     for (J=1;J<=3;J+)                   &

74、#160;  ROTI+2J+2=ROTI-1J-1;                     / COMPUTE ELEMENT LOCAL STIFFNESS MATRIX    ELST00=E*AX/D;         ELST03=-ELST00;   &#

75、160;     ELST11=12*E*YZ/(pow(D,3);         ELST12=6*E*YZ/(D*D);         ELST14=-ELST11;         ELST15=ELST12;         E

76、LST21=ELST12;         ELST22=4*E*YZ/D;         ELST24=-ELST12;         ELST25=2*E*YZ/D;         ELST30=ELST03;     

77、    ELST33=ELST00;         ELST41=ELST14;         ELST42=ELST24;         ELST44=ELST11;         ELST45=ELST24;  

78、60;      ELST51=ELST15;         ELST52=ELST25;         ELST54=ELST45;         ELST55=ELST22;/ ROTATE ELEMENT STIFFNESS MATRIX TO GLOBAL COORDINATES &#

79、160;  BTAB3(ELST,ROT,V,6,6);/         return;子程序BTAB3 void BTAB3(double A6, double B6, double V, int N, int NX)/ COMPUTE THE MATRIX OPERATION A=TRANSPOSE(B)*A*B/   int I, J, K;/ COMPUTE A*B AND STORE IN A   for (I=1;I<=N;I+) &#

80、160;            for (J=1;J<=N;J+)                         VJ-1=0.0;              for (K

81、=1;K<=N;K+)                                VJ-1=VJ-1+AI-1K-1*BK-1J-1;             

82、                          for (J=1;J<=N;J+)                    AI-1J-1=VJ-1;     &

83、#160;        / COMPUTE TRANSPOSE(B)*A AND STORE IN A   for (J=1;J<=N;J+)              for (I=1;I<=N;I+)                  

84、60;      VI-1=0.0;              for (K=1;K<=N;K+)                             

85、;   VI-1=VI-1+BK-1I-1*AK-1J-1;                                     for (I=1;I<=N;I+)      &#

86、160;             AI-1J-1=VI-1;              /   return;子程序ELASS SUBROUTINE ELASS(NEL,NCO,TM,ELMAT)!  STORE THE ELEMENT MATRIX FOR ELEMENT NEL IN THE TOTAL MATRIX!   

87、   COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G           DIMENSION NCO(1),TM(200,20),ELMAT(NDFEL,NDFEL)!      L1=NNE*(NEL-1)           DO I=1,NNE  &#

88、160;           L2=L1+I              N1=NCO(L2)              I1=NDF*(I-1)           &#

89、160;  J1=NDF*(N1-1)              DO J=1,NNE                 L2=L1+J                 N

90、2=NCO(L2)            I2=NDF*(J-1)            J2=NDF*(N2-1)                 DO K=1,NDF     

91、               KI=1                    IF(N1-N2).EQ.0) THEN !              

92、   STORE A DIAGONAL SUBMATRIX                  KI=K                    ENDIF       &#

93、160;          IF(N1-N2).LE.0) THEN!                 STORE AN OFF DIAGONAL SUBMATRIX               &#

94、160;  KR=J1+K                  IC=J2-KR+1                       K1=I1+K     

95、0;              ELSE!                 STORE THE TRANSPOSE OF AN OFF DIAGONAL MATRIX             

96、;     KR=J2+K                  IC=J1-KR+1                       K2=I2+K   &

97、#160;                ENDIF               DO L=KI,NDF                  KC=IC

98、+L                       IF(N1-N2).LE.0) THEN                     K2=I2+L   

99、0;              ELSE                       K1=I1+L              &

100、#160;        ENDIF                  TM(KR,KC)=TM(KR,KC)+ELMAT(K1,K2)               ENDDO   

101、0;             ENDDO              ENDDO           ENDDO!      RETURN        

102、;   END子程序BOUND/void ELASS(int NEL, int NCO, double TM20, double ELMAT5)/ STORE THE ELEMENT MATRIX FOR ELEMENT NEL IN THE TOTAL MATRIX/    int L1, I, L2, N1, I1, J1, J, N2, I2, J2, K, KI, KR, IC, K1, K2, L, KC;   /          

103、L1=NNE*(NEL-1);         for (I=1; I<=NNE; I+)                      L2=L1+I;             N1=NCOL2; 

104、0;           I1=NDF*(I-1);             J1=NDF*(N1-1);             for (J=1; J<=NNE; J+)         &

105、#160;                          L2=L1+J;                 N2=NCOL2;       

106、60;    I2=NDF*(J-1);            J2=NDF*(N2-1);                 for (K=1; K<=NDF; K+)           &#

107、160;                                     KI=1;               

108、;      if (N1-N2)=0)                                     /      &

109、#160;           STORE A DIAGONAL SUBMATRIX                    KI=K;              

110、0;                                         if (N1-N2)<=0)         &

111、#160;                            /                  STORE AN OFF DIAGONAL SUBMAT

112、RIX                    KR=J1+K;                    IC=J2-KR+1;        

113、                 K1=I1+K;                                 &#

114、160;                                        else          &#

115、160;                          /                  STORE THE TRANSPOSE OF AN OFF DIAGONAL MA

116、TRIX                    KR=J2+K;                    IC=J1-KR+1;        

117、;                 K2=I2+K;                                 &

118、#160;                   for (L=KI; L<=NDF; L+)                           

119、0;                             KC=IC+L;                      

120、;   if (N1-N2)<=0)                                             

121、0;                          K2=I2+L;                        

122、                                           else        

123、                                                   

124、;            K1=I1+L;                                     &#

125、160;                             TMKRKC=TMKRKC+ELMATK1K2;                   

126、;                                               /      

127、   return;子程序SLBSI SUBROUTINE SLBSI(A,B,D,N,MS,NX,MX)!  SOLUTION OF SIMUTANEOUS SYSTEMS OF EQUATIONS BY THE GAUSS!  ELIMINATION METHOD,FOR SYMMETRIC BANDED MATRICES!      DIMENSION A(NX,MX),B(NX),D(MX)!      N1=N-1  

128、0;        DO K=1,N1              C=A(K,1)              K1=K+1            &#

129、160; IF(ABS(C)-0.000001).LE.0) THEN                 WRITE(6,2) K   2        FORMAT(' * SINGULARITY IN ROW',I5)         

130、60;  STOP               ELSE!           DIVIDE ROW BY DIAGONAL COEFFICIENT            NI=K1+MS-2     

131、;       L=MIN0(NI,N)                 DO J=2,MS               D(J)=A(K,J)        

132、    ENDDO             DO J=K1,L                    K2=J-K+1            

133、0;  A(K,K2)=A(K,K2)/C            ENDDO            B(K)=B(K)/C!           ELIMINATE UNKNOWN X(K) FROM ROW I     

134、;       DO I=K1,L                    K2=I-K1+2                    C=D(K2) &

135、#160;                  DO J=I,L                       K2=J-I+1       &

136、#160;               K3=J-K+1                  A(I,K2)=A(I,K2)-C*A(K,K3)            &#

137、160;  ENDDO               B(I)=B(I)-C*B(K)                 ENDDO         ENDIF    

138、60;      ENDDO!  COMPUTE LAST UNKNOWN      IF(ABS(A(N,1)-0.000001).LE.0) THEN          WRITE(6,7) K   7     FORMAT(' * SINGULARITY IN ROW',I5)    &

139、#160;    STOP           ELSE         B(N)=B(N)/A(N,1)!  APPLY BACKSUBSTITUTE PROCESS TO COMPUTE REMAINING UNKNOWNS         DO I=1,N1   

140、              K=N-I                 K1=K+1                 NI=K1+MS-2 

141、                L=MIN0(NI,N)                 DO J=K1,L               &

142、#160;    K2=J-K+1               B(K)=B(K)-A(K,K2)*B(J)                 ENDDO          

143、;    ENDDO           ENDIF      RETURN           END子程序FORCE /void FORCE(int NCO, double PROP, double FORC, double REAC,      

144、0;               double X, double Y, double AL)/ COMPUTATION OF ELEMENT FORCES    int I, NEL, L, N1, N2, K1, K2, J1, J2, J, I1, I2;     double DX, DY, D, CO, SI, ROT66, U300, UL6, F6, AX, YZ, FG6;/ CLEAR T

145、HE REACTIONS ARRAY    for (I=1; I<=30; I+)                 REACI-1=0.0;/                 FORCI-1=0.0;    &#

146、160;        for (NEL=1; NEL<=NE; NEL+)                      L=NNE*(NEL-1);             N1=NCOL;   

147、0;         N2=NCOL+1;                   AX=PROPL;             YZ=PROPL+1;/      COMPUTE LENGTH

148、OF ELEMENT, AND SINE/COSINE OF ITS LOCAL X AXIS                   DX=XN2-1-XN1-1;             DY=YN2-1-YN1-1;       

149、0;     D=sqrt(DX*DX+DY*DY);             CO=DX/D;             SI=DY/D;/      FORM ELEMENT ROTATION MATRIX        ROT

150、00=CO;             ROT01=SI;             ROT02=0.0;             ROT10=-SI;           

151、  ROT11=CO;             ROT12=0.0;             ROT20=0.0;             ROT21=0.0;         &

152、#160;   ROT22=1.0;/      ROTATE ELEMENT NODAL DISPLACEMENTS TO ELEMENT/      LOCAL REFERENCE FRAME, AND STORE IN ARRAY UL        K1=NDF*(N1-1);           &

153、#160; K2=NDF*(N2-1);             for (I=1;I<=3;I+)                                   J1=K1+I;                J2=K2+I;                UI-1=ALJ1-1;