📄 truss.cpp
字号:
//
// STATIC ANALYSIS OF PLANE TRUSSES
//
#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[],
double AL[], int IB[], double REAC[]);
void ASSEM(double X[], double Y[], int NCO[], double PROP[],
double TK[][20], double ELST[][5], double AL[]);
void STIFF(int NEL, double X[], double Y[], int NCO[],
double PROP[], double ELST[][5], double AL[]);
void ELASS(int NEL, int NCO[], double TM[][20], double ELMAT[][5]);
void BOUND(double TK[][20], double AL[], double REAC[], int IB[]);
void SLBSI(double A[][20], double B[], double D[], int N, int MS,
int NRMX, int NCMX);
void FORCE(int NCO[], double PROP[], double FORC[], double REAC[],
double X[], double Y[], double AL[]);
void OUTPT(double AL[], double FORC[], double REAC[]);
// INITIALIZATION OF GLOBAL VARIABLES
int NN,NE,NLN,NBN,N,MS;
double E,G;
// ASSIGN DATA SET NUMBERS TO IN,FOR INPUT,AND IO FOR OUTPUT
ifstream READ_IN;
ofstream WRITE_IO;
//
// INITIALIZATION OF PROGRAM PARAMETERS
int NRMX=200;
int NCMX=20;
int NDF=2;
int NNE=2;
int NDFEL=NDF*NNE;
//
// MAIN PROGRAM
//
int main()
{
double X[100],Y[100],PROP[100],TK[200][20],
AL[200],FORC[100],REAC[200],ELST[5][5],V[20];
int NCO[200], IB[60];
//
// OPEN ALL FILES
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 BOUNDARY 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);
//
// OUTPUT
OUTPT(AL,FORC,REAC);
//
// CLOSE ALL FILES
READ_IN.close();
WRITE_IO.close();
return 0;
}
//
//
void INPUT(double X[], double Y[], int NCO[], double PROP[],
double AL[], int IB[], double REAC[])
{
// INPUT PROGRAM
//
int I, NUM, N1, IC[2], K, L, L1, L2, N2;
double W[3];
WRITE_IO.setf(ios::fixed);
WRITE_IO.setf(ios::showpoint);
WRITE_IO << " " <<
"*******************************************************************"
<< endl;
//
// READ BASIC PARAMETERS
READ_IN >> NN >> NE >> NLN >> NBN >> E;
WRITE_IO << "\n\n INTERNAL DATA \n\n" << " NUMBER OF NODES :"
<< setw(5) << NN << "\n" << " NUMBER OF ELEMENTS :" << setw(5)
<< NE << "\n" << " NUMBER OF LOADED NODES :" << setw(5) << NLN
<< "\n" << " NUMBER OF SUPPORT NODES :" << setw(5) << NBN << "\n"
<< " MODULUS OF ELASTICITY :" << setw(15) << setprecision(0) << E
<< "\n\n" << " NODAL COORDINATES\n" << setw(11) << "NODE" << setw(7)
<< "X" << setw(10) << "Y\n";
//
// READ NODAL COORDINATES IN ARRAY X AND Y
for (I=1; I<=NN; I++)
{
READ_IN >> NUM >> X[I] >> Y[I];
WRITE_IO.precision(2);
WRITE_IO << setw(10) << NUM << setw(10) << X[I] << setw(10) << Y[I]
<< "\n";
}
// READ ELEMENT CONNECTIVITY IN ARRAY NCO AND
// ELEMENT PROPERTIES IN ARRAY PROP
WRITE_IO << "\n ELEMENT CONNECTIVITY AND PROPERTIES\n" << setw(11)
<< "ELEMENT" << setw(23) << "START NODE END NODE" << setw(9)
<< "AREA" << endl;
for (I=1; I<=NE; I++)
{
N1=NNE*(I-1);
READ_IN >> NUM >> IC[0] >> IC[1] >> PROP[I];
WRITE_IO.precision(5);
WRITE_IO << setw(10) << NUM << setw(10) << IC[0] << setw(10) << IC[1]
<< setw(15) << PROP[I] << "\n";
NCO[N1+1]=IC[0];
NCO[N1+2]=IC[1];
}
//
// COMPUTE ACTUAL NUMBER OF UNKNOWNS AND CLEAR THE LOAD VECTOR
N=NN*NDF;
for (I=1; I<=N; I++)
{
REAC[I]=0.0;
AL[I]=0.0;
}
//
// READ THE NODAL LOADS AND STORE THEM IN ARRAY AL
WRITE_IO << "\n NODAL LOADS\n" << setw(11) << "NODE" << setw(7) << "PX"
<< setw(10) << "PY" << endl;
for (I=1; I<=NLN; I++)
{
READ_IN >> NUM >> W[0] >> W[1];
WRITE_IO.precision(2);
WRITE_IO << setw(10) << NUM << setw(10) << W[0] << setw(10) << W[1]
<< "\n";
for (K=1; K<=NDF; K++)
{
L=NDF*(NUM-1)+K;
AL[L]=W[K-1];
}
}
// READ BOUNDARY NODES DATA. STORE UNKNOWN STATUS INDICATORS
// IN ARRAY IB, AND PRESCRIBED UNKNOWN VALUES IN ARRAY REAC
WRITE_IO << "\n BOUNDARY CONDITION DATA\n" << setw(29) << "STATUS"
<< setw(31) << "PRESCRIBED VALUES\n" << setw(37)
<< "(0:PRESCRIBED, 1:FREE)\n" << setw(11) << "NODE" << setw(9)
<< "U" << setw(10) << "V" << setw(17) << "U" << setw(10) << "V"
<< endl;
for (I=1; I<=NBN; I++)
{
READ_IN >> NUM >> IC[0] >> IC[1] >> W[0] >> W[1];
WRITE_IO.precision(4);
WRITE_IO << setw(10) << NUM << setw(10) << IC[0] << setw(10) << IC[1]
<< setw(20) << W[0] << setw(10) << W[1] << "\n";
L1=(NDF+1)*(I-1)+1;
L2=NDF*(NUM-1);
IB[L1]=NUM;
for (K=1; K<=NDF; K++)
{
N1=L1+K;
N2=L2+K;
IB[N1]=IC[K-1];
REAC[N2]=W[K-1];
}
}
//
return;
}
//
//
void ASSEM(double X[], double Y[], int NCO[], double PROP[],
double TK[][20], double ELST[][5], double AL[])
{
// ASSEMBLING 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<=NE; I++)
{
L1=NNE*(I-1);
for (J=1; J<=N1; J++)
{
L2=L1+J;
J1=J+1;
for (K=J1; K<=NNE; K++)
{
L3=L1+K;
L=abs(NCO[L2]-NCO[L3]);
if ((MS-L)<=0)
{
MS=L;
}
}
}
}
MS=NDF*(MS+1);
//
//
// CLEAR THE TOTAL STIFFNESS MATRIX
for (I=1; I<=N; I++)
{
for (J=1; J<=MS; J++)
{
TK[I][J]=0.0;
}
}
for (I=1; I<=5; I++)
{
for (J=1; J<=5; J++)
{
ELST[I-1][J-1]=0.0;
}
}
//
for (NEL=1; NEL<=NE; NEL++)
{
// COMPUTE THE STIFFNESS MATRIX FOR ELEMENT NEL
STIFF(NEL,X,Y,NCO,PROP,ELST,AL);
// PLACE THE MATRIX IN THE TOTAL STIFFNESS MATRIX
ELASS(NEL,NCO,TK,ELST);
}
//
return;
}
//
//
void STIFF(int NEL, double X[], double Y[], int NCO[],
double PROP[], double ELST[][5], double AL[])
{
// COMPUTATION OF ELEMENT STIFFNESS MATRIX FOR CURRENT ELEMENT
//
int L, N1, N2, I, J, K1, K2;
double D, CO, SI, COEF;
//
//
L=NNE*(NEL-1);
N1=NCO[L+1];
N2=NCO[L+2];
//
// COMPUTE LENGTH OF ELEMENT, AND SINE AND COSINE OF ITS LOCAL X AXIS
D=sqrt(pow((X[N2]-X[N1]),2)+pow((Y[N2]-Y[N1]),2));
CO=(X[N2]-X[N1])/D;
SI=(Y[N2]-Y[N1])/D;
//
// COMPUTE ELEMENT STIFFNESS MATRIX
COEF= E*PROP[NEL]/D;
ELST[1][1]= COEF*CO*CO;
ELST[1][2]= COEF*CO*SI;
ELST[2][2]= COEF*SI*SI;
for (I=1; I<=2; I++)
{
for (J=I; J<=2; J++)
{
K1=I+NDF;
K2=J+NDF;
ELST[K1][K2]= ELST[I][J];
ELST[I][K2]= -ELST[I][J];
}
}
ELST[2][3]= -COEF*SI*CO;
//
return;
}
//
//
void ELASS(int NEL, int NCO[], double TM[][20], double ELMAT[][5])
{
// 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;
//
//
L1=NNE*(NEL-1);
for (I=1; I<=NNE; I++)
{
L2=L1+I;
N1=NCO[L2];
I1=NDF*(I-1);
J1=NDF*(N1-1);
for (J=I; J<=NNE; J++)
{
L2=L1+J;
N2=NCO[L2];
I2=NDF*(J-1);
J2=NDF*(N2-1);
for (K=1; K<=NDF; K++)
{
KI=1;
if ((N1-N2)==0)
{
//
// STORE A DIAGONAL SUBMATRIX
KI=K;
}
if ((N1-N2)<=0)
{
//
// STORE AN OFF DIAGONAL SUBMATRIX
KR=J1+K;
IC=J2-KR+1;
K1=I1+K;
}
else
{
//
// STORE THE TRANSPOSE OF AN OFF DIAGONAL MATRIX
KR=J2+K;
IC=J1-KR+1;
K2=I2+K;
}
for (L=KI; L<=NDF; L++)
{
KC=IC+L;
if ((N1-N2)<=0)
{
K2=I2+L;
}
else
{
K1=I1+L;
}
TM[KR][KC]=TM[KR][KC]+ELMAT[K1][K2];
}
}
}
}
//
return;
}
void BOUND(double TK[][20], double AL[], double REAC[], int IB[])
{
// INTRODUCTION OF THE BOUNDARY CONDITIONS
//
int L, L1, NO, K1, I, L2, KR, J, KV;
//
//
for (L=1; L<=NBN; L++)
{
L1=(NDF+1)*(L-1)+1;
NO=IB[L1];
K1=NDF*(NO-1);
for (I=1; I<=NDF; I++)
{
L2=L1+I;
if (IB[L2]==0)
{
// PRESCRIBED UNKNOWN TO BE CONSIDERED
KR=K1+I;
for (J=2; J<=MS; J++)
{
KV=KR+J-1;
if ((N-KV)>=0)
{
//
// MODIFY ROW OF TK AND CORRESPONDINF ELEMENTS IN AL
AL[KV]=AL[KV]-TK[KR][J]*REAC[KR];
TK[KR][J]=0.0;
}
KV=KR-J+1;
if (KV>0)
{
//
// MODIFY COLUMN IN TK AND CORRESPONDING ELEMENTS IN AL
AL[KV]=AL[KV]-TK[KV][J]*REAC[KR];
TK[KV][J]=0.0;
}
}
//
// SET DIAGONAL COEFFICIENT OF TK EQUAL TO 1 PLACE PRESCRIBED UNKNOWN
// VALUE IN AL
TK[KR][1]=1.0;
AL[KR]=REAC[KR];
}
}
}
//
return;
}
void SLBSI(double A[][20], double B[], double D[], int N, int MS,
int NRMX, int NCMX)
{
// SOLUTION OF SIMUTANEOUS SYSTEMS OF EQUATIONS BY THE GAUSS
// ELIMINATION METHOD,FOR SYMMETRIC BANDED MATRICES
int N1, K, K1, NI, L, J, K2, I, K3;
double C;
//
//
N1=N-1;
for (K=1; K<=N1; K++)
{
C=A[K][1];
K1=K+1;
if (C<=0.000001 && C>=-0.000001)
{
WRITE_IO << " **** SINGULARITY IN ROW" << setw(5) << K;
return;
}
else
{
//
// DIVIDE ROW BY DIAGONAL COEFFICIENT
NI=K1+MS-2;
if (NI<=N) {L=NI;} else {L=N;}
for (J=2; J<=MS; J++)
{
D[J]=A[K][J];
}
for (J=K1; J<=L; J++)
{
K2=J-K+1;
A[K][K2]=A[K][K2]/C;
}
B[K]=B[K]/C;
//
// ELIMINATE UNKNOWN X(K) FROM ROW I
for (I=K1; I<=L; I++)
{
K2=I-K1+2;
C=D[K2];
for (J=I; J<=L; J++)
{
K2=J-I+1;
K3=J-K+1;
A[I][K2]=A[I][K2]-C*A[K][K3];
}
B[I]=B[I]-C*B[K];
}
}
}
//
// COMPUTE LAST UNKNOWN
if (A[N][1]<=0.000001 && A[N][1]>=0.000001)
{
WRITE_IO << " **** SINGULARITY IN ROW" << setw(5) << K;
return;
}
else
{
B[N]=B[N]/A[N][1];
//
// APPLY BACKSUBSTITUTE PROCESS TO COMPUTE REMAINING UNKNOWNS
for (I=1; I<=N1; I++)
{
K=N-I;
K1=K+1;
NI=K1+MS-2;
if (NI<=N) {L=NI;} else {L=N;}
for (J=K1; J<=L; J++)
{
K2=J-K+1;
B[K]=B[K]-A[K][K2]*B[J];
}
}
}
//
return;
}
void FORCE(int NCO[], double PROP[], double FORC[], double REAC[],
double X[], double Y[], double AL[])
{
// COMPUTATION OF ELEMENT FORCES
int I, NEL, L, N1, N2, K1, K2;
double D, CO, SI, COEF;
//
// CLEAR THE REACTIONS ARRAY
for (I=1; I<=N; I++)
{
REAC[I]=0.0;
}
for (NEL=1; NEL<=NE; NEL++)
{
L=NNE*(NEL-1);
N1=NCO[L+1];
N2=NCO[L+2];
K1=NDF*(N1-1);
K2=NDF*(N2-1);
//
// COMPUTE LENGTH OF ELEMENT, AND SINE/COSINE OF ITS LOCAL X AXIS
D=sqrt(pow((X[N2]-X[N1]),2)+pow((Y[N2]-Y[N1]),2));
CO=(X[N2]-X[N1])/D;
SI=(Y[N2]-Y[N1])/D;
COEF=E*PROP[NEL]/D;
//
// COMPUTE MEMBER AXIAL FORCE AND STORE IN ARRAY FORC
FORC[NEL]=COEF*((AL[K2+1]-AL[K1+1])*CO+(AL[K2+2]-AL[K1+2])*SI);
//
// COMPUTE NODAL RESULTANTS
REAC[K1+1]=REAC[K1+1]-FORC[NEL]*CO;
REAC[K1+2]=REAC[K1+2]-FORC[NEL]*SI;
REAC[K2+1]=REAC[K2+1]+FORC[NEL]*CO;
REAC[K2+2]=REAC[K2+2]+FORC[NEL]*SI;
}
return;
}
//
//
void OUTPT(double AL[], double FORC[], double REAC[])
{
// OUTPUT PROGRAM
int I, K1, K2, J;
//
// WRITE NODAL DISPLACEMENTS
WRITE_IO <<
"\n\n **********************************************************************\n\n"
<< "RESULTS\n\n" << "NODAL DISPLACEMENTS\n" << setw(11) << "NODE"
<< setw(12) << "U" << setw(15) << "V" << endl;
for (I=1; I<=NN; I++)
{
K1=NDF*(I-1)+1;
K2=K1+NDF-1;
WRITE_IO.precision(6);
WRITE_IO << setw(10) << I;
for (J=K1; J<=K2; J++)
{
WRITE_IO << setw(15) << AL[J];
}
WRITE_IO << endl;
}
//
// WRITE NODAL REACTIONS
WRITE_IO << "\nNODAL REACTIONS\n" << setw(11) << "NODE" << setw(12) << "PX"
<< setw(15) << "PY\n";
for (I=1; I<=NN; I++)
{
K1=NDF*(I-1)+1;
K2=K1+NDF-1;
WRITE_IO << setw(10) << I;
for (J=K1; J<=K2; J++)
{
WRITE_IO << setw(15) << REAC[J];
}
WRITE_IO << endl;
}
//
// WRITE MEMBER AXIAL FORCES
WRITE_IO << "\nMEMBER FORCES" << setw(27) << "MEMBER AXIAL FORCE\n";
for (I=1; I<=NE; I++)
{
WRITE_IO << setw(10) << I << setw(15) << FORC[I] << endl;
}
WRITE_IO <<
"\n\n **********************************************************************\n";
return;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -