Dear all
How to write or have ``3D frame`` analysis for C++ ??
I want to look for Stiffness and element forces.
Thank you very much!!! :P
3D frame`` analysis for C++ ???
Moderators: silvia, selimgunay, Moderators
HI!!
I need to solve out C ++ procedure yard of 3D frame.
I have already finished 2D wins some The following:
#include<iostream>
#include<fstream>
#include<cmath>
#include<cstdio>
#include<direct.h>
using namespace std;
int unit;//---單位選擇
int NJ,NS,NMP,NCP,NM,NJL,NR,NDOF,JB,JE;
double BL,CX,CY,CZ,XB,YB,ZB,XE,YE,ZE,E,A;;
int **MSUP,**MPRP,*JP,*NSC,**B;
double **COORD,*CP,*EM,**S,**PJ,*P,*R;
double BK[2][2],GK[6][6],V[6],U[2],T[2][6],Q[2],F[6];
//================(SUBROUTINE 9) STORER() FUNCTION ================
//Storing member global force in support reaction vector
void STORER()
{
int i,I1,N;
for(i=1;i<=6;i++)
{
if(i<=3)
I1=(JB-1)*3+i;
else
I1=(JE-1)*3+(i-3);
N=NSC[I1];
if(N>NDOF)
R[N-NDOF-1]=R[N-NDOF-1]+F[i-1];
}
return;
}
//================(SUBROUTINE
MFORCEG() FUNCTION ================
//Member global force vector
void MFORCEG()
{
int i,j;
for(i=0;i<6;i++)
F=0;
for(i=0;i<6;i++)
for(j=0;j<2;j++)
F=F+T[j]*Q[j];
/*cout.setf(ios::fixed);cout.precision(4);
cout<<"---Global Force---\n";
for(j=0;j<6;j++)
{
cout.width(8);
cout<<F[j]<<endl;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 7) MFORCEL() FUNCTION ================
//Member local force vector
void MFORCEL(double *tmpforce)
{
int i,j;
for(j=0;j<2;j++)
Q[j]=0;
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
Q=Q+BK[j]*U[j];
}
/*cout<<"---Local Force---\n";
cout.setf(ios::fixed);cout.precision(4); */
*tmpforce=Q[0];
/* for(j=0;j<2;j++)
{
cout.width(8);
cout<<Q[j]<<endl;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 6) MSTIFFL() FUNCTION ================
//Member local stiffness matrix for plane trusses
void MSTIFFL()
{
double Z;
Z=E*A/BL;
BK[0][0]= Z;
BK[0][1]=-Z;
BK[1][0]=-Z;
BK[1][1]= Z;
return;
}
//================(SUBROUTINE 5) MDISPL() FUNCTION ================
//Member local displacement vector
void MDISPL()
{
int i,j;
for(i=0;i<2;i++)
U=0;
for(i=0;i<2;i++)
for(j=0;j<6;j++)
U=U+T[i][j]*V[j];
/*cout<<"---Local Disp.---\n";
cout.setf(ios::fixed);cout.precision(6);
for(i=0;i<2;i++)
{
cout.width(10);
cout<<U[i]<<endl;;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 4) MTRANS() FUNCTION ================
//Member transformation matrix for plane trusses
void MTRANS()
{
int i,j;
for(i=0;i<2;i++)
for(j=0;j<6;j++)
T[i][j]=0;
T[0][0]= CX;
T[0][1]= CY;
T[0][2]= CZ;
T[1][3]= CX;
T[1][4]= CY;
T[1][5]= CZ;
/* for(i=0;i<2;i++)
{
for(j=0;j<6;j++)
cout<<T[i][j];
cout<<endl;
}*/
return;
}
//================(SUBROUTINE 3) MDISPG() FUNCTION ================
//Member global displacement vector
void MDISPG()
{
int i,j,N;
for(i=0;i<6;i++)
V[i]=0;
j=(JB-1)*3;
// printf("================\n");
for(i=0;i<3;i++)
{
j++;
N=NSC[j];
if(N<=NDOF)
V[i]=P[N-1];
}
j=(JE-1)*3;
for(i=3;i<6;i++)
{
j++;
N=NSC[j];
if(N<=NDOF)
V[i]=P[N-1];
}
//for(i=0;i<6;i++)
//printf("%10.2E\n",V[i]);
return;
}
//==================================================
void GB()
{
double Z,Z1,Z2,Z3,Z4,Z5,Z6;
Z=E*A/BL;
Z1=Z*CX*CX;
Z2=Z*CY*CY;
Z3=Z*CZ*CZ;
Z4=Z*CX*CY;
Z5=Z*CX*CZ;
Z6=Z*CY*CZ;
GK[0][0]= Z1;GK[0][1]= Z4;GK[0][2]= Z5;GK[0][3]=-Z1;GK[0][4]=-Z4;GK[0][5]=-Z5;
GK[1][0]= Z4;GK[1][1]= Z2;GK[1][2]= Z6;GK[1][3]=-Z4;GK[1][4]=-Z2;GK[1][5]=-Z6;
GK[2][0]= Z5;GK[2][1]= Z6;GK[2][2]= Z3;GK[2][3]=-Z5;GK[2][4]=-Z6;GK[2][5]=-Z3;
GK[3][0]=-Z1;GK[3][1]=-Z4;GK[3][2]=-Z5;GK[3][3]= Z1;GK[3][4]= Z4;GK[3][5]= Z5;
GK[4][0]=-Z4;GK[4][1]=-Z2;GK[4][2]=-Z6;GK[4][3]= Z4;GK[4][4]= Z2;GK[4][5]= Z6;
GK[5][0]=-Z5;GK[5][1]=-Z6;GK[5][2]=-Z3;GK[5][3]= Z5;GK[5][4]= Z6;GK[5][5]= Z3;
/*cout.setf(ios::fixed);cout.precision(3);
for(i=0;i<6;i++)
{
for(j=0;j<6;j++)
{
cout.width(10);
cout<<GK[i][j];
}
cout<<endl;
}*/
return ;
}
//===============================================================================
void nsc()
{
int n=1;
int i,j;
for(i=0;i<NJ;i++)
{
for(j=0;j<NS;j++)
if(MSUP[j][0]==i+1)
{
B[i][0]=MSUP[j][1];
B[i][1]=MSUP[j][2];
B[i][2]=MSUP[j][3];
break;
}
else
{
B[i][0]=0;
B[i][1]=0;
B[i][2]=0;
}
}
/*for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
printf("%d\t",B[i][j]);
printf("\n");
}*/
for(i=0;i<NJ;i++)
for(j=0;j<3;j++)
if(B[i][j]==0)
{
NSC[i*3+j+1]=n;
n++;
}
for(i=0;i<NJ;i++)
for(j=0;j<3;j++)
if(B[i][j]==1)
{
NSC[i*3+j+1]=n;
n++;
}
/*printf("-----------------\n");
printf("NSC[%d]\n",NJ*3);
for(i=1;i<=(NJ*3);i++)
{
printf("%d\t",NSC[i]);
printf("\n");
} */
}
//-------------------------------------
void load()
{
int I1,I2,N,i,j;
for(i=0;i<NJL;i++)
{
I1=JP[i];
I2=(I1-1)*3;
for(j=0;j<3;j++)
{
I2=I2+1;
N=NSC[I2]-1;
if(N<NDOF)
P[N]=P[N]+PJ[i][j];
}
}
}
//------------------------------------
void Sys()
{
int i,j,I1,N1,N2;
for(i=0;i<6;i++)
{
if(i<3)
I1=(JB-1)*3+i;
else
I1=(JE-1)*3+(i-3);
N1=NSC[I1+1]-1;
if(N1<NDOF)
{
for(j=0;j<6;j++)
{
if(j<3)
I1=(JB-1)*3+j;
else
I1=(JE-1)*3+(j-3);
N2=NSC[I1+1]-1;
if(N2<NDOF)
S[N1][N2]+=GK[i][j];
}
}
}
/* cout.setf(ios::fixed);cout.precision(3);
for(i=0;i<NDOF;i++)
{
for(j=0;j<NDOF;j++)
cout<<S[i][j]<<"\t";
cout<<endl;
}*/
return;
}
//==============================================
int main()
{
int i,j,k=0,I,m,index;
double Z,Z1;
double **tmpdisp;
double *tmpforce;
double **tmpreaction;
char name[20];
char *subname=".txt";
char fold[20];
fstream in;
fstream out;
cout<<"(1)US (2)SI\n";
cout<<"Please choose the unit:";
cin>>unit;
if(unit==1)
unit=1;
else
unit=1000;
cout<<"-----------------------------------------------------------------\n";
cout<<"If you choose US unit\nPlease set Length(in) Area(in^2) Elastic-modulus(Ksi) Force(k)\n";
cout<<"If you choose SI unit\nPlease set Length(m) Area(mm^2) Elastic-modulus(Gpa) Force(kN)\n";
fflush(stdin);
cout<<"-----------------------------------------------------------------\n";
cout<<"Please input the data which you want to analysis(without txt):";
gets(name);
strcpy(fold,name);
mkdir(fold);
strcat(name,subname);
in.open(name,ios::in);
if(!in)
{
cout<<"Can't find the file\n";
exit(0);
}
chdir(fold);
out.open("result.txt",ios::out);
in>>NJ;
COORD=new double*[NJ];
for(i=0;i<NJ;i++)
{
COORD[i]=new double[3];
for(j=0;j<3;j++)
in>>COORD[i][j];
}
in>>NS;
MSUP=new int*[NS];
for(i=0;i<NS;i++)
{
MSUP[i]=new int[4];
for(j=0;j<4;j++)
in>>MSUP[i][j];
}
in>>NMP;
EM=new double[NMP];
for(i=0;i<NMP;i++)
in>>EM[i];
in>>NCP;
CP=new double[NCP];
for(i=0;i<NCP;i++)
in>>CP[i];
in>>NM;
MPRP=new int*[NM];
for(i=0;i<NM;i++)
{
MPRP[i]=new int[4];
for(j=0;j<4;j++)
in>>MPRP[i][j];
}
in>>NJL;
JP=new int[NJL];
PJ=new double*[NJL];
for(i=0;i<NJL;i++)
{
in>>JP[i];
PJ[i]=new double[3];
for(j=0;j<3;j++)
in>>PJ[i][j];
}
/*
for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
cout<<COORD[i][j]<<"\t";
cout<<endl;
}
cout<<"---------------------------------\n";
for(i=0;i<NS;i++)
{
for(j=0;j<4;j++)
cout<<MSUP[i][j]<<"\t";
cout<<endl;
}
cout<<"---------------------------------\n";
for(i=0;i<NMP;i++)
cout<<EM[i];
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NCP;i++)
cout<<CP[i];
//GB();
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NM;i++)
{
for(j=0;j<4;j++)
cout<<MPRP[i][j]<<"\t" ;
cout<<endl;
}
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NJL;i++)
{
cout<<JP[i]<<"\t";
for(j=0;j<3;j++)
cout<<PJ[i][j]<<"\t";
}*/
for(i=0;i<NS;i++)
for(j=1;j<4;j++)
k+=MSUP[i][j];
NDOF=NJ*3-k;
NSC=new int[NJ*3+1];
memset(NSC,0,sizeof(int)*(NJ*3+1));
B=new int*[NJ];
for(i=0;i<NJ;i++)
{
B[i]=new int[3];
memset(B[i],0,sizeof(int)*3 );
}
nsc();
P=new double[NDOF];//------------------配置LOAD VECTOR
memset(P,0,sizeof(double)*NDOF);
load();
S=new double*[NDOF];
for(i=0;i<NDOF;i++)
{
S[i]=new double[NDOF];
memset(S[i],0,sizeof(double)*NDOF);
}
for(i=0;i<NM;i++)
{
// cout<<"=========Member"<<i+1 <<"============"<<endl;
JB=MPRP[i][0];
JE=MPRP[i][1];
I=MPRP[i][2] ; E=EM[I-1];
I=MPRP[i][3] ; A=CP[I-1];
XB=COORD[JB-1][0];
YB=COORD[JB-1][1];
ZB=COORD[JB-1][2];
XE=COORD[JE-1][0];
YE=COORD[JE-1][1];
ZE=COORD[JE-1][2];
BL=sqrt((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB)+(ZE-ZB)*(ZE-ZB));
CX=(XE-XB)/BL;
CY=(YE-YB)/BL;
CZ=(ZE-ZB)/BL;
GB();
Sys();
}
/* cout.setf(ios::fixed);cout.precision(3);
cout<<"-----------------\n";
cout<<"P\n";
for(i=0;i<NDOF;i++)
{
cout.width(10);
cout<<P[i];
}
cout.setf(ios::fixed);cout.precision(3);
cout<<"------Big K--------\n";
for(i=0;i<NDOF;i++)
{
for(j=0;j<NDOF;j++)
{
cout.width(10);
cout<<S[i][j]<<"\t";
}
cout<<endl;
}*/
for(i=0;i<NDOF;i++)
{
Z1=S[i][i];
for(j=0;j<NDOF;j++)
S[i][j]=S[i][j]/Z1;
P[i]=P[i]/Z1;
for(k=0;k<NDOF;k++)
{
if(k!=i)
{
Z=S[k][i];
for(m=i;m<NDOF;m++)
S[k][m]=S[k][m]-S[i][m]*Z;
P[k]=P[k]-P[i]*Z;
}
}
}
/* cout<<"---------Joint DIS.-----\n";
cout.setf(ios::fixed);cout.precision(6);
for(i=0;i<NDOF;i++)
{
cout.width(10);
cout<<P[i]<<endl;
}*/
R=new double[NJ*3-NDOF];
memset(R,0,sizeof(double)*(NJ*3-NDOF));
tmpforce=new double[NM];
//======================================================
for(i=0;i<NM;i++)
{
//cout<<"=========Member"<<i+1 <<"============"<<endl;
JB=MPRP[i][0];
JE=MPRP[i][1];
I=MPRP[i][2] ; E=EM[I-1];
I=MPRP[i][3] ; A=CP[I-1];
XB=COORD[JB-1][0];
YB=COORD[JB-1][1];
ZB=COORD[JB-1][2];
XE=COORD[JE-1][0];
YE=COORD[JE-1][1];
ZE=COORD[JE-1][2];
BL=sqrt((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB)+(ZE-ZB)*(ZE-ZB));
CX=(XE-XB)/BL;
CY=(YE-YB)/BL;
CZ=(ZE-ZB)/BL;
MDISPG();
MTRANS();
MDISPL();
MSTIFFL();
MFORCEL(&tmpforce[i]);
MFORCEG();
STORER();
}
//delete memory of useless data
for(i=0;i<NJ;i++)
delete[] COORD[i];
delete[] COORD;
//-----------------------------
tmpdisp=new double*[NJ];
for(i=0;i<NJ;i++)
tmpdisp[i]=new double[3];
for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
{
if(NSC[i*3+1+j]<=NDOF)
tmpdisp[i][j]=P[NSC[i*3+1+j]-1];
else
tmpdisp[i][j]=0;
}
}
tmpreaction=new double*[NS];
for(i=0;i<NS;i++)
tmpreaction[i]=new double[3];
m=0;
for(i=0;i<NS;i++)
{
for(j=0;j<3;j++)
{
if(MSUP[i][j+1]==1)
{
tmpreaction[i][j]=R[m];
m++;
}
else
tmpreaction[i][j]=0;
}
}
/* cout<<"=======Reaction========\n";
for(i=0;i<(NJ*3-NDOF);i++)
cout<<R[i]<<endl;*/
out.setf(ios::scientific,ios::floatfield);
out.precision(4);
//==================================== Show Result in [Result.txt] ================
out<<"====================================================="<<endl;
out<<"============= Result of Analysis ===================="<<endl;
out<<"====================================================="<<endl;
out<<" ==================="<<endl;
out<<" Joint Displacement"<<endl;
out<<" ==================="<<endl;
out<<" Joint No. X Translation Y Translation Z Translation"<<endl;
out<<" --------- ------------- ------------- -------------"<<endl;
for(i=0;i<NJ;i++)
{
out.width(8);
out<<i+1;
for(j=0;j<3;j++)
{
out.width(20);
out<<tmpdisp[i][j]*unit;
}
out<<endl;
}
for(k=0;k<NJ;k++)
delete[] tmpdisp[k];
delete[] tmpdisp;
out<<endl;
out<<" ==================="<<endl;
out<<" Member Axial Forces"<<endl;
out<<" ==================="<<endl;
out<<" Member Axial Force Axial Stress"<<endl;
out<<" ------ -------------- ------------"<<endl;
for(i=0;i<NM;i++)
{
index=MPRP[i][3];
out.width(8);
out<<i+1;
out.width(20);
if(tmpforce[i]>=0)
out<<tmpforce[i]<<"(C) "<<tmpforce[i]/CP[index-1];
else
out<<tmpforce[i]<<"(T) "<<tmpforce[i]/CP[index-1];
out<<endl;
}
out<<endl;
delete[] tmpforce;
out<<" ================="<<endl;
out<<" Support Reactions"<<endl;
out<<" ================="<<endl;
out<<" Joint No. X Force Y Force Z Force"<<endl;
out<<" --------- ------------- ------------- -------------"<<endl;
for(i=0;i<NS;i++)
{
out.width(8);
out<<MSUP[i][0];
for(j=0;j<3;j++)
{
out.width(20);
out<<tmpreaction[i][j];
}
out<<endl;
}
cout<<"Analysis Over!!"<<endl;
return 0;
}
I need to solve out C ++ procedure yard of 3D frame.
I have already finished 2D wins some The following:
#include<iostream>
#include<fstream>
#include<cmath>
#include<cstdio>
#include<direct.h>
using namespace std;
int unit;//---單位選擇
int NJ,NS,NMP,NCP,NM,NJL,NR,NDOF,JB,JE;
double BL,CX,CY,CZ,XB,YB,ZB,XE,YE,ZE,E,A;;
int **MSUP,**MPRP,*JP,*NSC,**B;
double **COORD,*CP,*EM,**S,**PJ,*P,*R;
double BK[2][2],GK[6][6],V[6],U[2],T[2][6],Q[2],F[6];
//================(SUBROUTINE 9) STORER() FUNCTION ================
//Storing member global force in support reaction vector
void STORER()
{
int i,I1,N;
for(i=1;i<=6;i++)
{
if(i<=3)
I1=(JB-1)*3+i;
else
I1=(JE-1)*3+(i-3);
N=NSC[I1];
if(N>NDOF)
R[N-NDOF-1]=R[N-NDOF-1]+F[i-1];
}
return;
}
//================(SUBROUTINE
MFORCEG() FUNCTION ================//Member global force vector
void MFORCEG()
{
int i,j;
for(i=0;i<6;i++)
F=0;
for(i=0;i<6;i++)
for(j=0;j<2;j++)
F=F+T[j]*Q[j];
/*cout.setf(ios::fixed);cout.precision(4);
cout<<"---Global Force---\n";
for(j=0;j<6;j++)
{
cout.width(8);
cout<<F[j]<<endl;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 7) MFORCEL() FUNCTION ================
//Member local force vector
void MFORCEL(double *tmpforce)
{
int i,j;
for(j=0;j<2;j++)
Q[j]=0;
for(i=0;i<2;i++)
{
for(j=0;j<2;j++)
Q=Q+BK[j]*U[j];
}
/*cout<<"---Local Force---\n";
cout.setf(ios::fixed);cout.precision(4); */
*tmpforce=Q[0];
/* for(j=0;j<2;j++)
{
cout.width(8);
cout<<Q[j]<<endl;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 6) MSTIFFL() FUNCTION ================
//Member local stiffness matrix for plane trusses
void MSTIFFL()
{
double Z;
Z=E*A/BL;
BK[0][0]= Z;
BK[0][1]=-Z;
BK[1][0]=-Z;
BK[1][1]= Z;
return;
}
//================(SUBROUTINE 5) MDISPL() FUNCTION ================
//Member local displacement vector
void MDISPL()
{
int i,j;
for(i=0;i<2;i++)
U=0;
for(i=0;i<2;i++)
for(j=0;j<6;j++)
U=U+T[i][j]*V[j];
/*cout<<"---Local Disp.---\n";
cout.setf(ios::fixed);cout.precision(6);
for(i=0;i<2;i++)
{
cout.width(10);
cout<<U[i]<<endl;;
}
cout<<endl;*/
return;
}
//================(SUBROUTINE 4) MTRANS() FUNCTION ================
//Member transformation matrix for plane trusses
void MTRANS()
{
int i,j;
for(i=0;i<2;i++)
for(j=0;j<6;j++)
T[i][j]=0;
T[0][0]= CX;
T[0][1]= CY;
T[0][2]= CZ;
T[1][3]= CX;
T[1][4]= CY;
T[1][5]= CZ;
/* for(i=0;i<2;i++)
{
for(j=0;j<6;j++)
cout<<T[i][j];
cout<<endl;
}*/
return;
}
//================(SUBROUTINE 3) MDISPG() FUNCTION ================
//Member global displacement vector
void MDISPG()
{
int i,j,N;
for(i=0;i<6;i++)
V[i]=0;
j=(JB-1)*3;
// printf("================\n");
for(i=0;i<3;i++)
{
j++;
N=NSC[j];
if(N<=NDOF)
V[i]=P[N-1];
}
j=(JE-1)*3;
for(i=3;i<6;i++)
{
j++;
N=NSC[j];
if(N<=NDOF)
V[i]=P[N-1];
}
//for(i=0;i<6;i++)
//printf("%10.2E\n",V[i]);
return;
}
//==================================================
void GB()
{
double Z,Z1,Z2,Z3,Z4,Z5,Z6;
Z=E*A/BL;
Z1=Z*CX*CX;
Z2=Z*CY*CY;
Z3=Z*CZ*CZ;
Z4=Z*CX*CY;
Z5=Z*CX*CZ;
Z6=Z*CY*CZ;
GK[0][0]= Z1;GK[0][1]= Z4;GK[0][2]= Z5;GK[0][3]=-Z1;GK[0][4]=-Z4;GK[0][5]=-Z5;
GK[1][0]= Z4;GK[1][1]= Z2;GK[1][2]= Z6;GK[1][3]=-Z4;GK[1][4]=-Z2;GK[1][5]=-Z6;
GK[2][0]= Z5;GK[2][1]= Z6;GK[2][2]= Z3;GK[2][3]=-Z5;GK[2][4]=-Z6;GK[2][5]=-Z3;
GK[3][0]=-Z1;GK[3][1]=-Z4;GK[3][2]=-Z5;GK[3][3]= Z1;GK[3][4]= Z4;GK[3][5]= Z5;
GK[4][0]=-Z4;GK[4][1]=-Z2;GK[4][2]=-Z6;GK[4][3]= Z4;GK[4][4]= Z2;GK[4][5]= Z6;
GK[5][0]=-Z5;GK[5][1]=-Z6;GK[5][2]=-Z3;GK[5][3]= Z5;GK[5][4]= Z6;GK[5][5]= Z3;
/*cout.setf(ios::fixed);cout.precision(3);
for(i=0;i<6;i++)
{
for(j=0;j<6;j++)
{
cout.width(10);
cout<<GK[i][j];
}
cout<<endl;
}*/
return ;
}
//===============================================================================
void nsc()
{
int n=1;
int i,j;
for(i=0;i<NJ;i++)
{
for(j=0;j<NS;j++)
if(MSUP[j][0]==i+1)
{
B[i][0]=MSUP[j][1];
B[i][1]=MSUP[j][2];
B[i][2]=MSUP[j][3];
break;
}
else
{
B[i][0]=0;
B[i][1]=0;
B[i][2]=0;
}
}
/*for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
printf("%d\t",B[i][j]);
printf("\n");
}*/
for(i=0;i<NJ;i++)
for(j=0;j<3;j++)
if(B[i][j]==0)
{
NSC[i*3+j+1]=n;
n++;
}
for(i=0;i<NJ;i++)
for(j=0;j<3;j++)
if(B[i][j]==1)
{
NSC[i*3+j+1]=n;
n++;
}
/*printf("-----------------\n");
printf("NSC[%d]\n",NJ*3);
for(i=1;i<=(NJ*3);i++)
{
printf("%d\t",NSC[i]);
printf("\n");
} */
}
//-------------------------------------
void load()
{
int I1,I2,N,i,j;
for(i=0;i<NJL;i++)
{
I1=JP[i];
I2=(I1-1)*3;
for(j=0;j<3;j++)
{
I2=I2+1;
N=NSC[I2]-1;
if(N<NDOF)
P[N]=P[N]+PJ[i][j];
}
}
}
//------------------------------------
void Sys()
{
int i,j,I1,N1,N2;
for(i=0;i<6;i++)
{
if(i<3)
I1=(JB-1)*3+i;
else
I1=(JE-1)*3+(i-3);
N1=NSC[I1+1]-1;
if(N1<NDOF)
{
for(j=0;j<6;j++)
{
if(j<3)
I1=(JB-1)*3+j;
else
I1=(JE-1)*3+(j-3);
N2=NSC[I1+1]-1;
if(N2<NDOF)
S[N1][N2]+=GK[i][j];
}
}
}
/* cout.setf(ios::fixed);cout.precision(3);
for(i=0;i<NDOF;i++)
{
for(j=0;j<NDOF;j++)
cout<<S[i][j]<<"\t";
cout<<endl;
}*/
return;
}
//==============================================
int main()
{
int i,j,k=0,I,m,index;
double Z,Z1;
double **tmpdisp;
double *tmpforce;
double **tmpreaction;
char name[20];
char *subname=".txt";
char fold[20];
fstream in;
fstream out;
cout<<"(1)US (2)SI\n";
cout<<"Please choose the unit:";
cin>>unit;
if(unit==1)
unit=1;
else
unit=1000;
cout<<"-----------------------------------------------------------------\n";
cout<<"If you choose US unit\nPlease set Length(in) Area(in^2) Elastic-modulus(Ksi) Force(k)\n";
cout<<"If you choose SI unit\nPlease set Length(m) Area(mm^2) Elastic-modulus(Gpa) Force(kN)\n";
fflush(stdin);
cout<<"-----------------------------------------------------------------\n";
cout<<"Please input the data which you want to analysis(without txt):";
gets(name);
strcpy(fold,name);
mkdir(fold);
strcat(name,subname);
in.open(name,ios::in);
if(!in)
{
cout<<"Can't find the file\n";
exit(0);
}
chdir(fold);
out.open("result.txt",ios::out);
in>>NJ;
COORD=new double*[NJ];
for(i=0;i<NJ;i++)
{
COORD[i]=new double[3];
for(j=0;j<3;j++)
in>>COORD[i][j];
}
in>>NS;
MSUP=new int*[NS];
for(i=0;i<NS;i++)
{
MSUP[i]=new int[4];
for(j=0;j<4;j++)
in>>MSUP[i][j];
}
in>>NMP;
EM=new double[NMP];
for(i=0;i<NMP;i++)
in>>EM[i];
in>>NCP;
CP=new double[NCP];
for(i=0;i<NCP;i++)
in>>CP[i];
in>>NM;
MPRP=new int*[NM];
for(i=0;i<NM;i++)
{
MPRP[i]=new int[4];
for(j=0;j<4;j++)
in>>MPRP[i][j];
}
in>>NJL;
JP=new int[NJL];
PJ=new double*[NJL];
for(i=0;i<NJL;i++)
{
in>>JP[i];
PJ[i]=new double[3];
for(j=0;j<3;j++)
in>>PJ[i][j];
}
/*
for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
cout<<COORD[i][j]<<"\t";
cout<<endl;
}
cout<<"---------------------------------\n";
for(i=0;i<NS;i++)
{
for(j=0;j<4;j++)
cout<<MSUP[i][j]<<"\t";
cout<<endl;
}
cout<<"---------------------------------\n";
for(i=0;i<NMP;i++)
cout<<EM[i];
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NCP;i++)
cout<<CP[i];
//GB();
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NM;i++)
{
for(j=0;j<4;j++)
cout<<MPRP[i][j]<<"\t" ;
cout<<endl;
}
cout<<endl;
cout<<"---------------------------------\n";
for(i=0;i<NJL;i++)
{
cout<<JP[i]<<"\t";
for(j=0;j<3;j++)
cout<<PJ[i][j]<<"\t";
}*/
for(i=0;i<NS;i++)
for(j=1;j<4;j++)
k+=MSUP[i][j];
NDOF=NJ*3-k;
NSC=new int[NJ*3+1];
memset(NSC,0,sizeof(int)*(NJ*3+1));
B=new int*[NJ];
for(i=0;i<NJ;i++)
{
B[i]=new int[3];
memset(B[i],0,sizeof(int)*3 );
}
nsc();
P=new double[NDOF];//------------------配置LOAD VECTOR
memset(P,0,sizeof(double)*NDOF);
load();
S=new double*[NDOF];
for(i=0;i<NDOF;i++)
{
S[i]=new double[NDOF];
memset(S[i],0,sizeof(double)*NDOF);
}
for(i=0;i<NM;i++)
{
// cout<<"=========Member"<<i+1 <<"============"<<endl;
JB=MPRP[i][0];
JE=MPRP[i][1];
I=MPRP[i][2] ; E=EM[I-1];
I=MPRP[i][3] ; A=CP[I-1];
XB=COORD[JB-1][0];
YB=COORD[JB-1][1];
ZB=COORD[JB-1][2];
XE=COORD[JE-1][0];
YE=COORD[JE-1][1];
ZE=COORD[JE-1][2];
BL=sqrt((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB)+(ZE-ZB)*(ZE-ZB));
CX=(XE-XB)/BL;
CY=(YE-YB)/BL;
CZ=(ZE-ZB)/BL;
GB();
Sys();
}
/* cout.setf(ios::fixed);cout.precision(3);
cout<<"-----------------\n";
cout<<"P\n";
for(i=0;i<NDOF;i++)
{
cout.width(10);
cout<<P[i];
}
cout.setf(ios::fixed);cout.precision(3);
cout<<"------Big K--------\n";
for(i=0;i<NDOF;i++)
{
for(j=0;j<NDOF;j++)
{
cout.width(10);
cout<<S[i][j]<<"\t";
}
cout<<endl;
}*/
for(i=0;i<NDOF;i++)
{
Z1=S[i][i];
for(j=0;j<NDOF;j++)
S[i][j]=S[i][j]/Z1;
P[i]=P[i]/Z1;
for(k=0;k<NDOF;k++)
{
if(k!=i)
{
Z=S[k][i];
for(m=i;m<NDOF;m++)
S[k][m]=S[k][m]-S[i][m]*Z;
P[k]=P[k]-P[i]*Z;
}
}
}
/* cout<<"---------Joint DIS.-----\n";
cout.setf(ios::fixed);cout.precision(6);
for(i=0;i<NDOF;i++)
{
cout.width(10);
cout<<P[i]<<endl;
}*/
R=new double[NJ*3-NDOF];
memset(R,0,sizeof(double)*(NJ*3-NDOF));
tmpforce=new double[NM];
//======================================================
for(i=0;i<NM;i++)
{
//cout<<"=========Member"<<i+1 <<"============"<<endl;
JB=MPRP[i][0];
JE=MPRP[i][1];
I=MPRP[i][2] ; E=EM[I-1];
I=MPRP[i][3] ; A=CP[I-1];
XB=COORD[JB-1][0];
YB=COORD[JB-1][1];
ZB=COORD[JB-1][2];
XE=COORD[JE-1][0];
YE=COORD[JE-1][1];
ZE=COORD[JE-1][2];
BL=sqrt((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB)+(ZE-ZB)*(ZE-ZB));
CX=(XE-XB)/BL;
CY=(YE-YB)/BL;
CZ=(ZE-ZB)/BL;
MDISPG();
MTRANS();
MDISPL();
MSTIFFL();
MFORCEL(&tmpforce[i]);
MFORCEG();
STORER();
}
//delete memory of useless data
for(i=0;i<NJ;i++)
delete[] COORD[i];
delete[] COORD;
//-----------------------------
tmpdisp=new double*[NJ];
for(i=0;i<NJ;i++)
tmpdisp[i]=new double[3];
for(i=0;i<NJ;i++)
{
for(j=0;j<3;j++)
{
if(NSC[i*3+1+j]<=NDOF)
tmpdisp[i][j]=P[NSC[i*3+1+j]-1];
else
tmpdisp[i][j]=0;
}
}
tmpreaction=new double*[NS];
for(i=0;i<NS;i++)
tmpreaction[i]=new double[3];
m=0;
for(i=0;i<NS;i++)
{
for(j=0;j<3;j++)
{
if(MSUP[i][j+1]==1)
{
tmpreaction[i][j]=R[m];
m++;
}
else
tmpreaction[i][j]=0;
}
}
/* cout<<"=======Reaction========\n";
for(i=0;i<(NJ*3-NDOF);i++)
cout<<R[i]<<endl;*/
out.setf(ios::scientific,ios::floatfield);
out.precision(4);
//==================================== Show Result in [Result.txt] ================
out<<"====================================================="<<endl;
out<<"============= Result of Analysis ===================="<<endl;
out<<"====================================================="<<endl;
out<<" ==================="<<endl;
out<<" Joint Displacement"<<endl;
out<<" ==================="<<endl;
out<<" Joint No. X Translation Y Translation Z Translation"<<endl;
out<<" --------- ------------- ------------- -------------"<<endl;
for(i=0;i<NJ;i++)
{
out.width(8);
out<<i+1;
for(j=0;j<3;j++)
{
out.width(20);
out<<tmpdisp[i][j]*unit;
}
out<<endl;
}
for(k=0;k<NJ;k++)
delete[] tmpdisp[k];
delete[] tmpdisp;
out<<endl;
out<<" ==================="<<endl;
out<<" Member Axial Forces"<<endl;
out<<" ==================="<<endl;
out<<" Member Axial Force Axial Stress"<<endl;
out<<" ------ -------------- ------------"<<endl;
for(i=0;i<NM;i++)
{
index=MPRP[i][3];
out.width(8);
out<<i+1;
out.width(20);
if(tmpforce[i]>=0)
out<<tmpforce[i]<<"(C) "<<tmpforce[i]/CP[index-1];
else
out<<tmpforce[i]<<"(T) "<<tmpforce[i]/CP[index-1];
out<<endl;
}
out<<endl;
delete[] tmpforce;
out<<" ================="<<endl;
out<<" Support Reactions"<<endl;
out<<" ================="<<endl;
out<<" Joint No. X Force Y Force Z Force"<<endl;
out<<" --------- ------------- ------------- -------------"<<endl;
for(i=0;i<NS;i++)
{
out.width(8);
out<<MSUP[i][0];
for(j=0;j<3;j++)
{
out.width(20);
out<<tmpreaction[i][j];
}
out<<endl;
}
cout<<"Analysis Over!!"<<endl;
return 0;
}