wilsonmethod

进行结构响应分析的Wilson法

void CGlobalElement::WilsonMethod()
{
	int loop,loop1,iBuf;
	int nTimeStep,nLoadTimeStep,nGroundAccData,nTotalDOF,nFreeDOF;
	double *adGroundAcc;
	double dBuf,dBuf1,dAlfa,dPsi,da,da1,da2,da3,da4,da5;
	CSparseMatrix smatMass,smatDK;
	CMatrix matDisp,matVec,matAcc,matBuf,matBuf1;
	ofstream fout;
	ifstream fin;

	if(m_apEle.GetSize()==0) return;
	m_bFlagDynamicAnalyzed=true;

	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();

	matDisp.Realloc(nFreeDOF,1);
	matVec.Realloc(nFreeDOF,1);
	matAcc.Realloc(nFreeDOF,1);
	matBuf.Realloc(nFreeDOF,1);
	matBuf1.Realloc(nFreeDOF,1);

	GetMassMatrix(smatMass);
	m_smatGK.ElementBufRealloc();
	m_smatGK=0.0;
	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->StiffAssemble(m_smatGK);
	}

	dBuf=(m_adEigen[0]*m_adEigen[0]-m_adEigen[1]*m_adEigen[1]);
	dAlfa=2.0*m_adEigen[0]*m_adEigen[1]*(m_adDampRatio[1]*m_adEigen[0]
		-m_adDampRatio[0]*m_adEigen[1])/dBuf;
	dPsi=2.0*(m_adDampRatio[0]*m_adEigen[0]-m_adDampRatio[1]*m_adEigen[1])/dBuf;
	dBuf1=m_dWilsonXita*m_dTimeStep;
	dBuf=1.0/dBuf1;
	da=3.0*dBuf*(2.0*dBuf+dAlfa);
	da1=1.0+3.0*dBuf*dPsi;
	da2=3.0+dBuf1*dAlfa/2.0;
	da3=dBuf1*dPsi/2.0;
	da4=6.0*dBuf+3.0*dAlfa;
	da5=3.0*dPsi;
	
	nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
	nLoadTimeStep=(int)(m_dLoadDuration/m_dTimeStep);

	smatDK=m_smatGK;
	smatDK*=da1;
	smatMass*=da;
	smatDK+=smatMass;
	dBuf=1/da;
	smatMass*=dBuf;

	matDisp=0.0;matVec=0.0;matAcc=0.0;
	fout.open("dydata.tmp",ios::binary);
	if(m_iDynamicLoadType==DYNAMIC_INSTANT_LOAD){
		m_iCurLoadGroup=0;
		GetLoadVector();
		for(loop=0;loop<nTimeStep;loop++){
			matBuf=matVec*da4+matAcc*da2;
			matBuf=smatMass*matBuf;
			matBuf1=matVec*da5+matAcc*da3;
			matBuf1=m_smatGK*matBuf1;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				m_adDisp[loop1]=matBuf(loop1,0)+matBuf1(loop1,0);
			if(loop==0){
				for(loop1=0;loop1<nFreeDOF;loop1++)
					m_adDisp[loop1]+=m_dWilsonXita*m_adLoadVector[loop1];
			}
			else if(loop==nLoadTimeStep){
				for(loop1=0;loop1<nFreeDOF;loop1++)
					m_adDisp[loop1]-=m_dWilsonXita*m_adLoadVector[loop1];
			}

			if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				matBuf(loop1,0)=m_adDisp[loop1];
			dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
			dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
			matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
			matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
			matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
			matAcc=matBuf;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				m_adDisp[loop1]=matDisp(loop1,0);
			fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
		}
	}
	else{
		fin.open(m_sGroundAccFile,ios::nocreate);
		if(!fin.good()){
			AfxMessageBox("Can't open seismic acceleration file", MB_OK, 0 );
			return;
		}
		dBuf1=0.0;
		nGroundAccData=0;
		while(!fin.eof()){
			fin>>dBuf;
			nGroundAccData++;
			if(dBuf1<fabs(dBuf)) dBuf1=fabs(dBuf);
		}
		fin.close();
		dBuf=m_dPeakAcc/dBuf1;
		adGroundAcc=new double [nGroundAccData];
		fin.open(m_sGroundAccFile,ios::nocreate);
		adGroundAcc[0]=0.0;
		for(loop=1;loop<nGroundAccData;loop++){
			fin>>adGroundAcc[loop];
			adGroundAcc[loop]*=dBuf;
		}
		fin.close();

		for(loop=0;loop<nTimeStep;loop++){
			matBuf=matVec*da4+matAcc*da2;
			if(loop<nGroundAccData){
				for(loop1=0;loop1<m_nNode;loop1++){
					iBuf=m_Node.GetXDOFIndex(loop1);
					if(iBuf<nFreeDOF){
						matBuf(iBuf,0)+=m_dWilsonXita*(adGroundAcc[loop]-adGroundAcc[loop+1]);
					}
				}
			}
			matBuf=smatMass*matBuf;
			matBuf1=matVec*da5+matAcc*da3;
			matBuf1=m_smatGK*matBuf1;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				m_adDisp[loop1]=matBuf(loop1,0)+matBuf1(loop1,0);
			
			if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				matBuf(loop1,0)=m_adDisp[loop1];
			dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
			dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
			matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
			matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
			matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
			matAcc=matBuf;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				m_adDisp[loop1]=matDisp(loop1,0);
			fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
		}
		delete adGroundAcc;
	}
	fout.close();
}