globalelement.cpp

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

	int loop,loop1,nLoad,nMaterial,nLoadGroup;

	fout<<"结点数  "<<m_nNode<<endl;
	fout<<"结点编号  结点类型  结点坐标X  结点坐标Y"<<endl;
	for(loop=0;loop<m_nNode;loop++){
		fout<<loop<<"  "<<m_Node.GetType(loop)<<"  "<<
		m_Node.GetX(loop)<<"  "<<m_Node.GetY(loop)<<endl;
	}

	m_Node.OutputConstrainedNodeDescription(fout);

	fout<<endl<<"单元数  "<<m_nEle<<endl;
	for(loop=0;loop<m_nEle;loop++){
		fout<<"单元:"<<loop<<"  ";
		m_apEle[loop]->OutputParameter(fout);
	}

	nMaterial=m_apMaterial.GetSize();
	fout<<endl<<"材料总数  "<<nMaterial<<endl;
	for(loop=0;loop<nMaterial;loop++){
		fout<<"  材料编号:"<<loop<<"  ";
		m_apMaterial[loop]->OutputParameter(fout);
	}

	nLoadGroup=m_apLoadGroup.GetSize();
	fout<<endl<<"荷载组数: "<<nLoadGroup<<endl;
	for(loop=0;loop<nLoadGroup;loop++){
		fout<<"  荷载组 "<<loop;
		nLoad=m_apLoadGroup[loop]->m_apLoad.GetSize();
		fout<<"  荷载数:"<<nLoad<<endl;
		for(loop1=0;loop1<nLoad;loop1++){
			fout<<"    荷载:"<<loop1<<"  ";
			m_apLoadGroup[loop]->m_apLoad[loop1]->OutputParameter(fout);
		}
	}
	fout<<endl;

	if(m_bFlagForceAnalyzed){
		fout<<"荷载组 "<<m_iCurLoadGroup<<"作用下计算结果:";
		fout<<endl<<"结点位移:"<<endl;
		fout<<"结点编号  X向位移  Y向位移  转角位移"<<endl;
		for(loop=0;loop<m_nNode;loop++){
			fout<<loop<<"  ";
			fout<<m_Node.GetCurX(loop)-m_Node.GetX(loop)<<"  ";
			fout<<m_Node.GetCurY(loop)-m_Node.GetY(loop)<<"  ";
			if(m_Node.GetType(loop)==THREE_DOF)
				fout<<m_Node.GetCurR(loop);
			fout<<endl;
		}

		fout<<endl<<"单元内力:"<<endl;
		for(loop=0;loop<m_nEle;loop++){
			fout<<"单元: "<<loop<<"    "<<endl;
			m_apEle[loop]->OutputInternalForce(fout);
		}

		//appended in chapter5
		fout<<endl<<"支座反力:"<<endl;
		m_Node.OutputSupportReaction(fout,m_adLoadVector);
	}
	//appended in chapter7
	if(m_bFlagModeAnalyzed){
		fout<<endl<<"自振频率(1/s):"<<endl;
		for(loop=0;loop<m_nMode;loop++){
			fout<<"f"<<loop+1<<": "<<m_adEigen[loop]/(2.0*PI);
			fout<<"    ";
		}
		fout<<endl;
	}
}

void CGlobalElement::CalcuSupportReaction()
{
	int loop,loop1;
	int nTotalDOF,nFreeDOF;
	double dBuf;

	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();
	for(loop=nFreeDOF;loop<nTotalDOF;loop++){
		m_adLoadVector[loop]=-m_adLoadVector[loop];
		for(loop1=0;loop1<nTotalDOF;loop1++){
			dBuf=m_smatGK(loop,loop1);
			if(dBuf!=0.0){
				m_adLoadVector[loop]+=dBuf*m_adDisp[loop1];
			}
		}
	}
}

void CGlobalElement::ModeAnalysis()
{
	bool bLargerThanError;
	int loop,loop1;
	int nTotalDOF,nFreeDOF,nMode,nLumpMass;
	CUIntArray aiBuf;
	double dBuf,dBuf1,dError,dRatioGKtoMass,dQm;
	CArray <double,double&> adBuf,adEigenBak; 
	CMatrix matBuf,matBuf1;
	CMatrix matV,matU,matA,matB;
	CSparseMatrix smatGKBak;
	CSparseMatrix smatMass;

	if(m_apEle.GetSize()==0) return;
	m_bFlagModeAnalyzed=true;
	dError=1.0e-5;

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

	m_adDampRatio.SetSize(m_nMode);
	m_smatGK.ElementBufRealloc();
	m_smatGK=0.0;

	GetMassMatrix(smatMass);

	for(loop=0;loop<m_nEle;loop++){
		m_apEle[loop]->StiffAssemble(m_smatGK);
	}

	dBuf=0.0;
	for(loop=0;loop<nFreeDOF;loop++)
		dBuf+=m_smatGK(loop,loop);
	dBuf/=nFreeDOF;
	dBuf1=0.0;
	for(loop=0;loop<nFreeDOF;loop++){
		dBuf1+=smatMass(loop,loop);
	}
	if(m_iMassMatrixType==CONSISTENT_MASS){
		dBuf1/=nFreeDOF;
	}
	else{
		nLumpMass=m_LumpMass.m_apLoad.GetSize();
		dBuf1/=nLumpMass*2.0;
	}
	dRatioGKtoMass=dBuf/dBuf1;

	dBuf=1.0/dBuf;
	m_smatGK*=dBuf;

	dQm=1.0/dBuf1;
	smatMass*=dQm;

	aiBuf.SetSize(nFreeDOF);
	adBuf.SetSize(nFreeDOF);
	dBuf=0.0;
	for(loop=0;loop<nFreeDOF;loop++){
		adBuf[loop]=-1.0;
		if(smatMass(loop,loop)>0.0){
			adBuf[loop]=m_smatGK(loop,loop)/smatMass(loop,loop);
		}
		if(dBuf<adBuf[loop]) dBuf=adBuf[loop];	
	}

	for(loop=0;loop<nFreeDOF;loop++){
		if(adBuf[loop]<0.0) adBuf[loop]=dBuf+100.0;	
	}

	aiBuf[0]=0;
	for(loop=0;loop<nFreeDOF;loop++){
		for(loop1=0;loop1<loop;loop1++){
			if(adBuf[aiBuf[loop1]]>adBuf[loop]){
				aiBuf.InsertAt(loop1,loop);
				break;
			}
		}
		if(loop1==loop){
			aiBuf[loop1]=loop;
		}
	}

	nMode=m_nMode+8;
	if(2*m_nMode<nMode) nMode=2*m_nMode;
	if(nMode>nFreeDOF)	nMode=nFreeDOF;
	if(m_iMassMatrixType==LUMP_MASS&&nMode>2*nLumpMass) nMode=2*nLumpMass;
	if(m_nMode>nMode)	m_nMode=nMode;

	matV.Realloc(nFreeDOF,nMode);
	matU.Realloc(nFreeDOF,nMode);
	matA.Realloc(nMode,nMode);
	matB.Realloc(nMode,nMode);

	m_adEigen.SetSize(nMode);
	adEigenBak.SetSize(nMode);
	for(loop=0;loop<nMode;loop++)
		m_adEigen[loop]=0.0;

	matV=0.0;
	for(loop=0;loop<nFreeDOF;loop++)
		matV(loop,0)=1.0;
	for(loop=1;loop<nMode;loop++)
		matV(aiBuf[loop],loop)=1.0;

	smatGKBak=m_smatGK;
	do{
		bLargerThanError=false;
		for(loop=0;loop<nMode;loop++)
			adEigenBak[loop]=m_adEigen[loop];
		matU=smatMass*matV;

		if(!m_smatGK.LdltSolve(nFreeDOF,matU))	return;
		matBuf=smatGKBak*matU;
		matBuf1.Trans(matU);
		matA=matBuf1*matBuf;
		matBuf=smatMass*matU;
		matB=matBuf1*matBuf;

		JacobiEigen(matA,matB);
		matV=matU*matB;
		for(loop=0;loop<nMode;loop++){
			dBuf=fabs((m_adEigen[loop]-adEigenBak[loop])/m_adEigen[loop]);
			if(dBuf>dError){
				bLargerThanError=true;
				break;
			}
		}
	}while(bLargerThanError);

	dBuf=sqrt(dQm);
	m_matEigenVector.Realloc(nTotalDOF,m_nMode);
	m_matEigenVector=0.0;
	for(loop=0;loop<m_nMode;loop++)
	for(loop1=0;loop1<nFreeDOF;loop1++)
		m_matEigenVector(loop1,loop)=matV(loop1,loop)*dBuf;

	for(loop=0;loop<m_nMode;loop++){
		m_adEigen[loop]*=dRatioGKtoMass;
		m_adEigen[loop]=sqrt(m_adEigen[loop]);
	}
	m_smatGK.ElementBufEmpty();
}

void CGlobalElement::JacobiEigen(CMatrix &matA, CMatrix &matB)
{
	CMatrix matBBak,matR,matRT,matV;
	double dAlfa,dBeta,dBuf,dBuf1;
	double dA,dB,dC;
	double dError;
	int nRow;
	int iRow,iCol;
	int loop,loop1;
	CArray<double,double&> adBuf;
	CUIntArray aiBuf;

	nRow=matA.GetRow();
	matBBak.Realloc(nRow,nRow);
	matR.Realloc(nRow,nRow);
	matRT.Realloc(nRow,nRow);
	matV.Realloc(nRow,nRow);
	matBBak=matB;

	matV=0.0;
	for(loop=0;loop<nRow;loop++)
		matV(loop,loop)=1.0;

	dError=1.0e-24;

	dBuf=0.0;
	for(loop=0;loop<nRow;loop++){
		for(loop1=loop+1;loop1<nRow;loop1++){
			dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));
			if(dBuf<fabs(dBuf1)){
				dBuf=fabs(dBuf1);
				iRow=loop;iCol=loop1;
			}
			dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));
			if(dBuf<fabs(dBuf1)){
				dBuf=fabs(dBuf1);
				iRow=loop;iCol=loop1;
			}
		}
	}

	while(dBuf>dError){
		matR=0.0;
		for(loop=0;loop<nRow;loop++) matR(loop,loop)=1.0;
		dA=matA(iRow,iRow)*matB(iRow,iCol)-matB(iRow,iRow)*matA(iRow,iCol);
		dC=-matA(iCol,iCol)*matB(iRow,iCol)+matB(iCol,iCol)*matA(iRow,iCol);
		if(dA==0.0&&dC==0.0){
			dAlfa=0.0;
			dBeta=-matA(iRow,iCol)/matA(iCol,iCol);
		}
		else{
			dB=matA(iRow,iRow)*matB(iCol,iCol)-matA(iCol,iCol)*matB(iRow,iRow);
			if(dB>=0.0){
				dBuf=dB/2.0;
				dAlfa=-dC/(dBuf+sqrt(dBuf*dBuf-dA*dC));
			}
			else{
				dBuf=dB/2.0;
				dAlfa=-dC/(dBuf-sqrt(dBuf*dBuf-dA*dC));
			}
			dBeta=dA/dC*dAlfa;
		}
		matR(iRow,iCol)=dAlfa;
		matR(iCol,iRow)=dBeta;
		matRT.Trans(matR);

		matA=matRT*matA*matR;
		matB=matRT*matB*matR;
		matV*=matR;

		dBuf=0.0;
		for(loop=0;loop<nRow;loop++){
			for(loop1=loop+1;loop1<nRow;loop1++){
				dBuf1=matA(loop,loop1)*matA(loop,loop1)/(matA(loop,loop)*matA(loop1,loop1));
				if(dBuf<fabs(dBuf1)){
					dBuf=fabs(dBuf1);
					iRow=loop;iCol=loop1;
				}
				dBuf1=matB(loop,loop1)*matB(loop,loop1)/(matB(loop,loop)*matB(loop1,loop1));
				if(dBuf<fabs(dBuf1)){
					dBuf=fabs(dBuf1);
					iRow=loop;iCol=loop1;
				}
			}
		}
	}

	matRT.Trans(matV);
	matR=matRT*matBBak;
	for(loop=0;loop<nRow;loop++){
		dBuf=0.0;
		for(loop1=0;loop1<nRow;loop1++){
			dBuf+=matR(loop,loop1)*matV(loop1,loop);
		}
		if(dBuf<0.0) dBuf=-dBuf;
		dBuf=sqrt(dBuf);
		for(loop1=0;loop1<nRow;loop1++)
			matV(loop1,loop)/=dBuf;
	}

	adBuf.SetSize(nRow);
	for(loop=0;loop<nRow;loop++){
		adBuf[loop]=matA(loop,loop)/matB(loop,loop);
	}
	aiBuf.SetSize(nRow);
	aiBuf[0]=0;
	for(loop=0;loop<nRow;loop++){
		for(loop1=0;loop1<loop;loop1++){
			if(fabs(adBuf[aiBuf[loop1]])>fabs(adBuf[loop])){
				aiBuf.InsertAt(loop1,loop);
				break;
			}
		}
		if(loop1==loop){
			aiBuf[loop1]=loop;
		}
	}

	for(loop=0;loop<nRow;loop++){
		m_adEigen[loop]=adBuf[aiBuf[loop]];
	}
	for(loop=0;loop<nRow;loop++){
		for(loop1=0;loop1<nRow;loop1++){
			matB(loop1,loop)=matV(loop1,aiBuf[loop]);
		}
	}
}

void CGlobalElement::ModeSuperposition()
{
	int loop,loop1,loop2,iBuf;
	int nTotalDOF,nFreeDOF;
	int nTimeStep,nLoadTimeStep,nGroundAccData;
	double	*adP,*adY,*adGroundAcc;
	double dAlfa,dPsi,de;
	double dBuf,dBuf1,*adBuf;
	double dOmigaD;
	double dTime;
	double	*dAti,*dBti;
	ofstream fout;
	ifstream fin;
	CMatrix matBuf,matBuf1;
	CSparseMatrix smatMass;

	if(m_apEle.GetSize()==0) return;
	m_bFlagDynamicAnalyzed=true;
	de=2.718282;
	adP=new double [m_nMode];
	adY=new double [m_nMode];

	nTotalDOF=m_Node.GetTotalDOF();
	nFreeDOF=m_Node.GetFreeDOF();
	adBuf=new double [nFreeDOF];

	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;
	
	for(loop=2;loop<m_nMode;loop++)
		m_adDampRatio[loop]=0.5*(dPsi*m_adEigen[loop]+dAlfa/m_adEigen[loop]);

	dAti=new double [m_nMode];
	dBti=new double [m_nMode];
	for(loop=0;loop<m_nMode;loop++){
		dAti[loop]=dBti[loop]=0.0;
	}

	fout.open("dydata.tmp",ios::binary);
	if(m_iDynamicLoadType==DYNAMIC_INSTANT_LOAD){
		m_iCurLoadGroup=0;
		GetLoadVector();

		for(loop=0;loop<m_nMode;loop++){
			adP[loop]=0.0;
			for(loop1=0;loop1<nFreeDOF;loop1++){
				adP[loop]+=m_matEigenVector(loop1,loop)
							*m_adLoadVector[loop1];
			}
		}
		
		nLoadTimeStep=(int)(m_dLoadDuration/m_dTimeStep);
		nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
		for(loop=0;loop<nTimeStep;loop++){
			dTime=loop*m_dTimeStep;
			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);
				if(loop<nLoadTimeStep){
					dAti[loop1]=dBuf*(dAti[loop1]+adP[loop1]*cos(dOmigaD*dTime)*m_dTimeStep);
					dBti[loop1]=dBuf*(dBti[loop1]+adP[loop1]*sin(dOmigaD*dTime)*m_dTimeStep);
				}
				else{
					dAti[loop1]=dBuf*dAti[loop1];
					dBti[loop1]=dBuf*dBti[loop1];
				}
				adY[loop1]=(dAti[loop1]*sin(dOmigaD*dTime)-dBti[loop1]*cos(dOmigaD*dTime))/dOmigaD;
			}