进行结构响应分析的Wilson法

风花雪月

进行结构响应分析的Wilson法，是直接积分的基本算法之一，相信进行有限元动力分析编程的诸位都能用上（这是vc代码）

1. void CGlobalElement::WilsonMethod()
   
2. {
   
3.         int loop,loop1,iBuf;
   
4.         int nTimeStep,nLoadTimeStep,nGroundAccData,nTotalDOF,nFreeDOF;
   
5.         double *adGroundAcc;
   
6.         double dBuf,dBuf1,dAlfa,dPsi,da,da1,da2,da3,da4,da5;
   
7.         CSparseMatrix smatMass,smatDK;
   
8.         CMatrix matDisp,matVec,matAcc,matBuf,matBuf1;
   
9.         ofstream fout;
   
10.         ifstream fin;
    
11. 
    
12.         if(m_apEle.GetSize()==0) return;
    
13.         m_bFlagDynamicAnalyzed=true;
    
14. 
    
15.         nTotalDOF=m_Node.GetTotalDOF();
    
16.         nFreeDOF=m_Node.GetFreeDOF();
    
17. 
    
18.         matDisp.Realloc(nFreeDOF,1);
    
19.         matVec.Realloc(nFreeDOF,1);
    
20.         matAcc.Realloc(nFreeDOF,1);
    
21.         matBuf.Realloc(nFreeDOF,1);
    
22.         matBuf1.Realloc(nFreeDOF,1);
    
23. 
    
24.         GetMassMatrix(smatMass);
    
25.         m_smatGK.ElementBufRealloc();
    
26.         m_smatGK=0.0;
    
27.         for(loop=0;loop<m_nEle;loop++){
    
28.                 m_apEle[loop]->StiffAssemble(m_smatGK);
    
29.         }
    
30. 
    
31.         dBuf=(m_adEigen[0]*m_adEigen[0]-m_adEigen[1]*m_adEigen[1]);
    
32.         dAlfa=2.0*m_adEigen[0]*m_adEigen[1]*(m_adDampRatio[1]*m_adEigen[0]
    
33.                 -m_adDampRatio[0]*m_adEigen[1])/dBuf;
    
34.         dPsi=2.0*(m_adDampRatio[0]*m_adEigen[0]-m_adDampRatio[1]*m_adEigen[1])/dBuf;
    
35.         dBuf1=m_dWilsonXita*m_dTimeStep;
    
36.         dBuf=1.0/dBuf1;
    
37.         da=3.0*dBuf*(2.0*dBuf+dAlfa);
    
38.         da1=1.0+3.0*dBuf*dPsi;
    
39.         da2=3.0+dBuf1*dAlfa/2.0;
    
40.         da3=dBuf1*dPsi/2.0;
    
41.         da4=6.0*dBuf+3.0*dAlfa;
    
42.         da5=3.0*dPsi;
    
43.        
    
44.         nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
    
45.         nLoadTimeStep=(int)(m_dLoadDuration/m_dTimeStep);
    
46. 
    
47.         smatDK=m_smatGK;
    
48.         smatDK*=da1;
    
49.         smatMass*=da;
    
50.         smatDK+=smatMass;
    
51.         dBuf=1/da;
    
52.         smatMass*=dBuf;
    
53. 
    
54.         matDisp=0.0;matVec=0.0;matAcc=0.0;
    
55.         fout.open("dydata.tmp",ios::binary);
    
56.         if(m_iDynamicLoadType==DYNAMIC_INSTANT_LOAD){
    
57.                 m_iCurLoadGroup=0;
    
58.                 GetLoadVector();
    
59.                 for(loop=0;loop<nTimeStep;loop++){
    
60.                         matBuf=matVec*da4+matAcc*da2;
    
61.                         matBuf=smatMass*matBuf;
    
62.                         matBuf1=matVec*da5+matAcc*da3;
    
63.                         matBuf1=m_smatGK*matBuf1;
    
64.                         for(loop1=0;loop1<nFreeDOF;loop1++)
    
65.                                 m_adDisp[loop1]=matBuf(loop1,0)+matBuf1(loop1,0);
    
66.                         if(loop==0){
    
67.                                 for(loop1=0;loop1<nFreeDOF;loop1++)
    
68.                                         m_adDisp[loop1]+=m_dWilsonXita*m_adLoadVector[loop1];
    
69.                         }
    
70.                         else if(loop==nLoadTimeStep){
    
71.                                 for(loop1=0;loop1<nFreeDOF;loop1++)
    
72.                                         m_adDisp[loop1]-=m_dWilsonXita*m_adLoadVector[loop1];
    
73.                         }
    
74. 
    
75.                         if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
    
76.                         for(loop1=0;loop1<nFreeDOF;loop1++)
    
77.                                 matBuf(loop1,0)=m_adDisp[loop1];
    
78.                         dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
    
79.                         dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
    
80.                         matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
    
81.                         matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
    
82.                         matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
    
83.                         matAcc=matBuf;
    
84.                         for(loop1=0;loop1<nFreeDOF;loop1++)
    
85.                                 m_adDisp[loop1]=matDisp(loop1,0);
    
86.                         fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
    
87.                 }
    
88.         }
    
89.         else{
    
90.                 fin.open(m_sGroundAccFile,ios::nocreate);
    
91.                 if(!fin.good()){
    
92.                         AfxMessageBox("Can't open seismic acceleration file", MB_OK, 0 );
    
93.                         return;
    
94.                 }
    
95.                 dBuf1=0.0;
    
96.                 nGroundAccData=0;
    
97.                 while(!fin.eof()){
    
98.                         fin>>dBuf;
    
99.                         nGroundAccData++;
    
100.                         if(dBuf1<fabs(dBuf)) dBuf1=fabs(dBuf);
     
101.                 }
     
102.                 fin.close();
     
103.                 dBuf=m_dPeakAcc/dBuf1;
     
104.                 adGroundAcc=new double [nGroundAccData];
     
105.                 fin.open(m_sGroundAccFile,ios::nocreate);
     
106.                 adGroundAcc[0]=0.0;
     
107.                 for(loop=1;loop<nGroundAccData;loop++){
     
108.                         fin>>adGroundAcc[loop];
     
109.                         adGroundAcc[loop]*=dBuf;
     
110.                 }
     
111.                 fin.close();
     
112. 
     
113.                 for(loop=0;loop<nTimeStep;loop++){
     
114.                         matBuf=matVec*da4+matAcc*da2;
     
115.                         if(loop<nGroundAccData){
     
116.                                 for(loop1=0;loop1<m_nNode;loop1++){
     
117.                                         iBuf=m_Node.GetXDOFIndex(loop1);
     
118.                                         if(iBuf<nFreeDOF){
     
119.                                                 matBuf(iBuf,0)+=m_dWilsonXita*(adGroundAcc[loop]-adGroundAcc[loop+1]);
     
120.                                         }
     
121.                                 }
     
122.                         }
     
123.                         matBuf=smatMass*matBuf;
     
124.                         matBuf1=matVec*da5+matAcc*da3;
     
125.                         matBuf1=m_smatGK*matBuf1;
     
126.                         for(loop1=0;loop1<nFreeDOF;loop1++)
     
127.                                 m_adDisp[loop1]=matBuf(loop1,0)+matBuf1(loop1,0);
     
128.                        
     
129.                         if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
     
130.                         for(loop1=0;loop1<nFreeDOF;loop1++)
     
131.                                 matBuf(loop1,0)=m_adDisp[loop1];
     
132.                         dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
     
133.                         dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
     
134.                         matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
     
135.                         matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
     
136.                         matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
     
137.                         matAcc=matBuf;
     
138.                         for(loop1=0;loop1<nFreeDOF;loop1++)
     
139.                                 m_adDisp[loop1]=matDisp(loop1,0);
     
140.                         fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
     
141.                 }
     
142.                 delete adGroundAcc;
     
143.         }
     
144.         fout.close();
     
145. }

复制代码

jxndwl

正需要啊，
谢谢