globalelement.cpp

一个计算悬臂梁的有限元vc源码

			for(loop1=0;loop1<nFreeDOF;loop1++){
				adBuf[loop1]=0.0;
				for(loop2=0;loop2<m_nMode;loop2++){
					adBuf[loop1]+=adY[loop2]*m_matEigenVector(loop1,loop2);
				}
			}
			fout.write((char*)adBuf,sizeof(double)*nFreeDOF);

		}
	}
	else{
		fin.open(m_sGroundAccFile,ios::nocreate);
		if(!fin.good()){
			AfxMessageBox("Can't open seismic acceleration file", MB_OK, 0 );
			return;
		}
		nGroundAccData=-1;
		dBuf1=0.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);
		for(loop=0;loop<nGroundAccData;loop++){
			fin>>adGroundAcc[loop];
			adGroundAcc[loop]*=dBuf;
		}
		fin.close();

		GetMassMatrix(smatMass);
		matBuf.Realloc(nFreeDOF,1);
		matBuf1.Realloc(nFreeDOF,1);
		nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
		for(loop=0;loop<nTimeStep;loop++){
			dTime=loop*m_dTimeStep;

			if(loop<nGroundAccData){
				for(loop1=0;loop1<nFreeDOF;loop1++){
					matBuf(loop1,0)=0.0;
				}
				for(loop1=0;loop1<m_nNode;loop1++){
					iBuf=m_Node.GetXDOFIndex(loop1);
					if(iBuf<nFreeDOF){
						matBuf(iBuf,0)=-adGroundAcc[loop];
					}
				}
				//matBuf1=smatMass.MultiplyMat(nFreeDOF,matBuf);
				matBuf1=smatMass*matBuf;
			}
			else{
				for(loop1=0;loop1<nFreeDOF;loop1++)
					matBuf1(loop1,0)=0.0;
			}

			for(loop1=0;loop1<m_nMode;loop1++){
				adP[loop1]=0.0;
				for(loop2=0;loop2<nFreeDOF;loop2++){
					adP[loop1]+=m_matEigenVector(loop2,loop1)
						*matBuf1(loop2,0);
				}
			}

			for(loop1=0;loop1<m_nMode;loop1++){
				dBuf=m_adDampRatio[loop1];
				dBuf*=dBuf;
				dOmigaD=m_adEigen[loop1]*sqrt(1.0-dBuf);
				dBuf=-m_adDampRatio[loop1]*m_adEigen[loop1]*m_dTimeStep;
				dBuf=pow(de,dBuf);
				dAti[loop1]=dBuf*(dAti[loop1]+adP[loop1]*cos(dOmigaD*dTime)*m_dTimeStep);
				dBti[loop1]=dBuf*(dBti[loop1]+adP[loop1]*sin(dOmigaD*dTime)*m_dTimeStep);
				adY[loop1]=(dAti[loop1]*sin(dOmigaD*dTime)-dBti[loop1]*cos(dOmigaD*dTime))/dOmigaD;
			}
			for(loop1=0;loop1<nFreeDOF;loop1++){
				adBuf[loop1]=0.0;
				for(loop2=0;loop2<m_nMode;loop2++){
					adBuf[loop1]+=adY[loop2]*m_matEigenVector(loop1,loop2);
				}
			}
			fout.write((char*)adBuf,sizeof(double)*nFreeDOF);
		}
		delete adGroundAcc;
	}
	fout.close();
	delete adP;
	delete adY;
	delete dAti;
	delete dBti;
	delete adBuf;
}

void CGlobalElement::NewmarkMethod()
{
	int loop,loop1,iBuf;
	int nTotalDOF,nFreeDOF;
	int nTimeStep,nLoadTimeStep,nGroundAccData;
	double dBuf,dBuf1;
	double dAlfa,dPsi;
	double *adGroundAcc;
	double da,da1,da2,da3;
	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;
	da=1.0+dAlfa*m_dNewmarkGama*m_dTimeStep;
	da1=dPsi*m_dNewmarkGama*m_dTimeStep+m_dNewmarkBeta*m_dTimeStep*m_dTimeStep;
	da2=-dAlfa*m_dTimeStep;
	da3=dPsi*m_dTimeStep+m_dTimeStep*m_dTimeStep/2.0;

	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();
		nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
		nLoadTimeStep=(int)(m_dLoadDuration/m_dTimeStep);
		for(loop=0;loop<nTimeStep;loop++){
			matBuf=matAcc*da2;
			matBuf=smatMass*matBuf;
			matBuf1=matVec*m_dTimeStep+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_adLoadVector[loop1];
			}
			else if(loop==nLoadTimeStep){
				for(loop1=0;loop1<nFreeDOF;loop1++)
					m_adDisp[loop1]-=m_adLoadVector[loop1];
			}

			if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				matBuf1(loop1,0)=m_adDisp[loop1];
			matDisp+=matVec*m_dTimeStep+matAcc*(m_dTimeStep*m_dTimeStep/2.0)
				+matBuf1*(m_dNewmarkBeta*m_dTimeStep*m_dTimeStep);
			matVec+=matAcc*m_dTimeStep+matBuf1*(m_dNewmarkGama*m_dTimeStep);
			matAcc+=matBuf1;
			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();

		nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
		for(loop=0;loop<nTimeStep;loop++){
			matBuf=matAcc*da2;
			if(loop<nGroundAccData){
				for(loop1=0;loop1<m_nNode;loop1++){
					iBuf=m_Node.GetXDOFIndex(loop1);
					if(iBuf<nFreeDOF){
						matBuf(iBuf,0)+=adGroundAcc[loop]-adGroundAcc[loop+1];
					}
				}
			}
			matBuf=smatMass*matBuf;
			matBuf1=matVec*m_dTimeStep+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++)
				matBuf1(loop1,0)=m_adDisp[loop1];
			matDisp+=matVec*m_dTimeStep+matAcc*(m_dTimeStep*m_dTimeStep/2.0)
				+matBuf1*(m_dNewmarkBeta*m_dTimeStep*m_dTimeStep);
			matVec+=matAcc*m_dTimeStep+matBuf1*(m_dNewmarkGama*m_dTimeStep);
			matAcc+=matBuf1;
			for(loop1=0;loop1<nFreeDOF;loop1++)
				m_adDisp[loop1]=matDisp(loop1,0);
			fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
		}
		delete adGroundAcc;
	}
	fout.close();
}

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();
}

void CGlobalElement::GetMassMatrix(CSparseMatrix& smatMass)
{
	int loop,iBuf,iNodeDOFIndex,nLumpMass;
	unsigned long* aiBuf;
	if(m_iMassMatrixType==CONSISTENT_MASS){
		smatMass=m_smatGK;
		smatMass.ElementBufRealloc();
		smatMass=0.0;
		for(loop=0;loop<m_nEle;loop++){
			m_apEle[loop]->MassAssemble(smatMass);
		}

		nLumpMass=m_LumpMass.m_apLoad.GetSize();
		for(loop=0;loop<nLumpMass;loop++){
			iBuf=m_LumpMass.m_apLoad[loop]->GetLoadedNode();
			iNodeDOFIndex=m_Node.GetXDOFIndex(iBuf);
			smatMass(iNodeDOFIndex,iNodeDOFIndex)+=
				m_LumpMass.m_apLoad[loop]->GetLoad();
			iNodeDOFIndex=m_Node.GetYDOFIndex(iBuf);
			smatMass(iNodeDOFIndex,iNodeDOFIndex)+=
				m_LumpMass.m_apLoad[loop]->GetLoad();
		}
	}
	else{
		iBuf=m_Node.GetTotalDOF();
		aiBuf=new unsigned long [iBuf];
		for(loop=0;loop<iBuf;loop++)
			aiBuf[loop]=loop;
		smatMass.Realloc(iBuf,aiBuf);
		smatMass.ElementBufRealloc();
		smatMass=0.0;
		nLumpMass=m_LumpMass.m_apLoad.GetSize();
		if(nLumpMass==0&&m_iMassMatrixType==LUMP_MASS){
			AfxMessageBox("No mass !", MB_OK, 0 );
			return;
		}
		for(loop=0;loop<nLumpMass;loop++){
			iBuf=m_LumpMass.m_apLoad[loop]->GetLoadedNode();
			iNodeDOFIndex=m_Node.GetXDOFIndex(iBuf);
			smatMass(iNodeDOFIndex,iNodeDOFIndex)+=
				m_LumpMass.m_apLoad[loop]->GetLoad();
			iNodeDOFIndex=m_Node.GetYDOFIndex(iBuf);
			smatMass(iNodeDOFIndex,iNodeDOFIndex)+=
				m_LumpMass.m_apLoad[loop]->GetLoad();
		}
		delete aiBuf;
	}
}