⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 plane8nodeele.cpp

📁 平面8节点单元类
💻 CPP
字号:
🔍 
// Plane8NodeEle.cpp: implementation of the CPlane8NodeEle class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "afxtempl.h"
#include "math.h"
#include "fstream.h"
#include "Matrix.h"
#include "SparseMatrix.h"
#include "Node.h"
#include "BaseMaterial.h"
#include "Plane8NodeEle.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

#define PLANE_8NODE_ELE 3
#define PI 3.1415926535897932384626433832795

CMatrix CPlane8NodeEle::m_matKe=CMatrix(16,16);
CMatrix CPlane8NodeEle::m_mat0301=CMatrix(3,1);
CMatrix CPlane8NodeEle::m_mat0303=CMatrix(3,3);
CMatrix CPlane8NodeEle::m_mat0316=CMatrix(3,16);
CMatrix CPlane8NodeEle::m_mat1603=CMatrix(16,3);
CMatrix CPlane8NodeEle::m_mat1602=CMatrix(16,2);
CMatrix CPlane8NodeEle::m_mat0216=CMatrix(2,16);
CMatrix CPlane8NodeEle::m_matNodeDisp=CMatrix(16,1);

double CPlane8NodeEle::m_adGaussKxi[9]={-0.7745967,0.0,0.7745967,
				-0.7745967,0.0,0.7745967,-0.7745967,0.0,0.7745967};
double CPlane8NodeEle::m_adGaussEta[9]={-0.7745967,-0.7745967,
				-0.7745967,0.0,0.0,0.0,0.7745967,0.7745967,0.7745967};


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CPlane8NodeEle::CPlane8NodeEle()
{

}

CPlane8NodeEle::~CPlane8NodeEle()
{

}

void CPlane8NodeEle::ReadNode(int& iCurCharPos, CString &sData)
{
	int iBuf;
	for(int loop=0;loop<8;loop++)
		iBuf=m_aiNode[loop]=ReadInt(iCurCharPos,sData);
}

void CPlane8NodeEle::GetStiffness()
{
	int loop;
	double dJacobi,dWeight;

	m_matKe=0.0;
	for(loop=0;loop<9;loop++)
	{
		GetB(loop,m_mat0316,dJacobi,dWeight);
		m_mat1603.Trans(m_mat0316);
		//m_mat0303=[D];
		GetD(loop,m_mat0303);
		m_mat1603*=m_mat0303;
		m_mat0316*=dJacobi*dWeight;
		m_matKe+=m_mat1603*m_mat0316;
	}
	m_matKe*=m_dThickness;
}

void CPlane8NodeEle::GetB(int iGaussPoint, CMatrix &matB,double& dJacobi,double& dWeight)
{
	int loop;
	double dEta,dKxi,dBuf1,dBuf2,dBuf3,dBuf4,dBuf5,dBuf6,dBuf7,dBuf;
	double m,n;
	double adX[8],adY[8];

	for(loop=0;loop<8;loop++){
		adX[loop]=m_pNode->GetX(m_aiNode[loop]);	
		adY[loop]=m_pNode->GetY(m_aiNode[loop]);
	}
		
	matB=0.0;
	dEta=m_adGaussEta[iGaussPoint];
	dKxi=m_adGaussKxi[iGaussPoint];

	m=1+dEta;n=1-dEta;
	//dN/ddKxi*Xi++;
	dBuf1=(-n/4.0+dKxi*n/2.0+m*n/4.0)*adX[0]
		-dKxi*n*adX[1]
		+(n/4.0+dKxi*n/2.0-m*n/4.0)*adX[2]
		+m*n/2.0*adX[3]
		+(m/4.0+dKxi*m/2.0-m*n/4.0)*adX[4]
		-dKxi*m*adX[5]
		+(-m/4.0+dKxi*m/2.0+m*n/4.0)*adX[6]
		-m*n/2.0*adX[7];
	
	//dN/ddKxi*Yi++;
	dBuf2=(-n/4.0+dKxi*n/2.0+m*n/4.0)*adY[0]
		-dKxi*n*adY[1]
		+(n/4.0+dKxi*n/2.0-m*n/4.0)*adY[2]
		+m*n/2.0*adY[3]
		+(m/4.0+dKxi*m/2.0-m*n/4.0)*adY[4]
		-dKxi*m*adY[5]
		+(-m/4.0+dKxi*m/2.0+m*n/4.0)*adY[6]
		-m*n/2.0*adY[7];
	
	m=1+dKxi;n=1-dKxi;	
	//dN/ddEta*Xi++;
	dBuf3=(-n/4.0+m*n/4.0+dEta*n/2.0)*adX[0]
		-m*n/2.0*adX[1]
		+(-m/4.0+m*n/4.0+dEta*m/2.0)*adX[2]
		-dEta*m*adX[3]
		+(m/4.0-m*n/4.0+dEta*m/2.0)*adX[4]
		+m*n/2.0*adX[5]
		+(n/4.0-m*n/4.0+dEta*n/2.0)*adX[6]
		-dEta*n*adX[7];
	
	//dN/ddEta*Yi++;
	dBuf4=(-n/4.0+m*n/4.0+dEta*n/2.0)*adY[0]
		-m*n/2.0*adY[1]
		+(-m/4.0+m*n/4.0+dEta*m/2.0)*adY[2]
		-dEta*m*adY[3]
		+(m/4.0-m*n/4.0+dEta*m/2.0)*adY[4]
		+m*n/2.0*adY[5]
		+(n/4.0-m*n/4.0+dEta*n/2.0)*adY[6]
		-dEta*n*adY[7];
	
	// |1   2| 
	// |3   4| 

	dBuf5=dBuf1*dBuf4-dBuf3*dBuf2;
	//date|J|
	dJacobi=dBuf5;
	
	//Inv[J]=
	// |4   -2|
	// |-3   1|/date|J|

	dBuf6=dBuf1;
	dBuf2*=-1;
	dBuf3*=-1;
	dBuf1=dBuf4;
	dBuf4=dBuf6;

	dBuf1/=dBuf5;
	dBuf2/=dBuf5;
	dBuf3/=dBuf5;
	dBuf4/=dBuf5;

	//dN1/ddKxi;
	dBuf6=-(1-dEta)/4.0+dKxi*(1-dEta)/2.0+(1-dEta*dEta)/4.0;
	//dN1/ddEta;
	dBuf7=-(1-dKxi)/4.0+(1-dKxi*dKxi)/4.0+dEta*(1-dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,0)=dBuf;
	matB(2,1)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,1)=dBuf;
	matB(2,0)=dBuf;

	//dN2/ddKxi;
	dBuf6=-dKxi*(1-dEta);
	//dN2/ddEta;
	dBuf7=-(1-dKxi*dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,2)=dBuf;
	matB(2,3)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,3)=dBuf;
	matB(2,2)=dBuf;

	//dN3/ddKxi;
	dBuf6=(1-dEta)/4.0+dKxi*(1-dEta)/2.0-(1-dEta*dEta)/4.0;
	//dN3/ddEta;
	dBuf7=-(1+dKxi)/4.0+(1-dKxi*dKxi)/4.0+dEta*(1+dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,4)=dBuf;
	matB(2,5)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,5)=dBuf;
	matB(2,4)=dBuf;

	//dN4/ddKxi;
	dBuf6=(1-dEta*dEta)/2.0;
	//dN4/ddEta;
	dBuf7=-dEta*(1+dKxi);

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,6)=dBuf;
	matB(2,7)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,7)=dBuf;
	matB(2,6)=dBuf;

	//dN5/ddKxi;
	dBuf6=(1+dEta)/4.0+dKxi*(1+dEta)/2.0-(1-dEta*dEta)/4.0;
	//dN5/ddEta;
	dBuf7=(1+dKxi)/4.0-(1-dKxi*dKxi)/4.0+dEta*(1+dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,8)=dBuf;
	matB(2,9)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,9)=dBuf;
	matB(2,8)=dBuf;

	//dN6/ddKxi;
	dBuf6=-dKxi*(1+dEta);
	//dN6/ddEta;
	dBuf7=(1-dKxi*dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,10)=dBuf;
	matB(2,11)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,11)=dBuf;
	matB(2,10)=dBuf;

	//dN7/ddKxi;
	dBuf6=-(1+dEta)/4.0+dKxi*(1+dEta)/2.0+(1-dEta*dEta)/4.0;
	//dN7/ddEta;
	dBuf7=(1-dKxi)/4.0-(1-dKxi*dKxi)/4.0+dEta*(1-dKxi)/2.0;

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,12)=dBuf;
	matB(2,13)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,13)=dBuf;
	matB(2,12)=dBuf;

	//dN8/ddKxi;
	dBuf6=-(1-dEta*dEta)/2.0;
	//dN8/ddEta;
	dBuf7=-dEta*(1-dKxi);

	dBuf=dBuf1*dBuf6+dBuf2*dBuf7;
	matB(0,14)=dBuf;
	matB(2,15)=dBuf;

	dBuf=dBuf3*dBuf6+dBuf4*dBuf7;
	matB(1,15)=dBuf;
	matB(2,14)=dBuf;

	if(dEta==0.0||dKxi==0.0){
		if(dKxi!=0.0||dEta!=0.0){
		dWeight=0.55555556*0.88888889;
		}
		else{
		dWeight=0.88888889*0.88888889;
		}
	}
	else{
		dWeight=0.55555556*0.55555556;
	}
}

void CPlane8NodeEle::GetD(int iGaussPoint, CMatrix &D)
{
	double dBuf,dE,dMiu;
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;
	dE=apMaterial[m_iMaterialIndex]->GetE();
	dMiu=apMaterial[m_iMaterialIndex]->GetMiu();
	dBuf=dE/(1.0-dMiu*dMiu);
	
	D(0,0)=dBuf;	D(0,1)=dMiu*dBuf; D(0,2)=0.0;
	D(1,0)=D(0,1);	D(1,1)=dBuf;		D(1,2)=0.0;
	D(2,0)=0.0;		D(2,1)=0.0;			D(2,2)=dBuf*(1-dMiu)/2.0;
}

void CPlane8NodeEle::OutputParameter(ofstream &fout)
{
	int loop;
	fout<<"Plane8NodeElement  单元厚度: "<<m_dThickness;
	fout<<"  材料类型编号:"<<m_iMaterialIndex;
	fout<<"    结点号:";
	for(loop=0;loop<8;loop++){
		fout<<m_aiNode[loop]<<"  ";
	}
	fout<<endl;
}

void CPlane8NodeEle::StiffAssemble(CSparseMatrix& smatGK)
{
	int loop,loop1,iBuf,aiDOFIndex[16];
	//int nFreeDOF;
	GetStiffness();
	for(loop=0;loop<8;loop++){
		aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<16;loop++){
		//if(aiDOFIndex[loop]<nFreeDOF){
			iBuf=aiDOFIndex[loop];
			smatGK(iBuf,iBuf)+=m_matKe(loop,loop);
			for(loop1=0;loop1<loop;loop1++){
				//if(aiDOFIndex[loop1]<nFreeDOF){
					smatGK(iBuf,aiDOFIndex[loop1])+=m_matKe(loop,loop1);
				//}
			}
		//}
	}
}

void CPlane8NodeEle::OutputInternalForce(ofstream &fout)
{
	int loop;
	
	fout<<"高斯积分点:  0   1   2   3   4   5   6   7   8"<<endl;

	fout<<"X向应变: ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adStrainX[loop]<<"  ";
	}
	fout<<endl;
	fout<<"Y向应变: ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adStrainY[loop]<<"  ";
	}
	fout<<endl;
	fout<<"剪应变 ";
	for(loop=0;loop<9;loop++){
		fout<<m_adStrainXY[loop]<<"  ";
	}
	fout<<endl;

	fout<<"X向应力 ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adStressX[loop]<<"  ";
	}
	fout<<endl;
	fout<<"Y向应力: ";
	for(loop=0;loop<9;loop++){
		fout<<m_adStressY[loop]<<"  ";
	}
	fout<<endl;
	fout<<"剪应力";	
	for(loop=0;loop<9;loop++){
		fout<<m_adStressXY[loop]<<"  ";
	}
	fout<<endl;

	fout<<"主应力1: ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adMainStress[loop]<<"  ";
	}
	fout<<endl;
	fout<<"主应力2: ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adMainStress1[loop]<<"  ";
	}
	fout<<endl;
	fout<<"主应力1方向: ";	
	for(loop=0;loop<9;loop++){
		fout<<m_adMainStressAngle[loop]*180.0/PI<<"  ";
	}
	fout<<endl;
}

int CPlane8NodeEle::GetElementType()
{
	return PLANE_8NODE_ELE;
}

void CPlane8NodeEle::ReadParameter(int &iCurCharPos, CString &sData)
{
	m_dThickness=ReadDouble(iCurCharPos,sData);
	m_iMaterialIndex=ReadInt(iCurCharPos,sData);
}

void CPlane8NodeEle::InternalForceInitial()
{
	for(int loop=0;loop<9;loop++){
		m_adStressX[loop]=0.0;
		m_adStressY[loop]=0.0;
		m_adStressXY[loop]=0.0;
	}
}

void CPlane8NodeEle::CalcuInternalForce(const double* adNodeDisp)
{
	int loop,iBuf;
	double dBuf,dBuf1;

	for(loop=0;loop<8;loop++){
		iBuf=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(2*loop,0)=adNodeDisp[iBuf];
		iBuf=m_pNode->GetYDOFIndex(m_aiNode[loop]);
		m_matNodeDisp(2*loop+1,0)=adNodeDisp[iBuf];
	}

	for(loop=0;loop<9;loop++){
		GetB(loop,m_mat0316,dBuf,dBuf);
		m_mat0301=m_mat0316*m_matNodeDisp;
		m_adStrainX[loop]=m_mat0301(0,0);	
		m_adStrainY[loop]=m_mat0301(1,0);
		m_adStrainXY[loop]=m_mat0301(2,0);
	}

	for(loop=0;loop<9;loop++){
		m_mat0301(0,0)=m_adStrainX[loop];
		m_mat0301(1,0)=m_adStrainY[loop];
		m_mat0301(2,0)=m_adStrainXY[loop];
		GetD(loop,m_mat0303);
		m_mat0301=m_mat0303*m_mat0301;
		m_adStressX[loop]=m_mat0301(0,0);
		m_adStressY[loop]=m_mat0301(1,0);
		m_adStressXY[loop]=m_mat0301(2,0);
	}

	for(loop=0;loop<9;loop++){
		m_adMainStressAngle[loop]
			=0.5*atan2(2.0*m_adStressXY[loop],m_adStressX[loop]-m_adStressY[loop]);
		dBuf=m_adStressX[loop]+m_adStressY[loop];
		dBuf1=m_adStressX[loop]-m_adStressY[loop];
		dBuf1=sqrt(dBuf1*dBuf1+4.*m_adStressXY[loop]*m_adStressXY[loop]);
		m_adMainStress[loop]=0.5*dBuf+0.5*dBuf1;
		m_adMainStress1[loop]=0.5*dBuf-0.5*dBuf1;
	}
}

void CPlane8NodeEle::CalcuGKBandWidth(unsigned long* aiGKDiagAdd)
{
	int loop,iBuf,iMinDOFIndex,aiDOFIndex[16];

	for(loop=0;loop<8;loop++)
	{	aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}

	iMinDOFIndex=aiDOFIndex[0];
	for(loop=1;loop<16;loop++){
		if(iMinDOFIndex>aiDOFIndex[loop]) 
			iMinDOFIndex=aiDOFIndex[loop];
	}

	for(loop=0;loop<16;loop++){
		iBuf=aiDOFIndex[loop]-iMinDOFIndex+1;
		if(aiGKDiagAdd[aiDOFIndex[loop]]<iBuf)
			aiGKDiagAdd[aiDOFIndex[loop]]=iBuf;
	}
}

void CPlane8NodeEle::GetNode(int *aiNode, int &nEleNode)
{
	for(int loop=0;loop<8;loop++)
		aiNode[loop]=m_aiNode[loop];
	nEleNode=8;
}

void CPlane8NodeEle::GetGaussPtGCoordinate(double* adGaussX, double* adGaussY)
{
	int loop;
	double m,n,r,s;
	double adX[8],adY[8];
	for(loop=0;loop<8;loop++){
		adX[loop]=m_pNode->GetX(m_aiNode[loop]);	
		adY[loop]=m_pNode->GetY(m_aiNode[loop]);
	}

	for(loop=0;loop<9;loop++){
		m=1.0+m_adGaussKxi[loop];n=1.0-m_adGaussKxi[loop];
		r=1.0+m_adGaussEta[loop];s=1.0-m_adGaussEta[loop];
		//coordinate_x=Ni*Xi
		adGaussX[loop]=n*s*(1-m-r)/4.0*adX[0]
			+m*n*s/2.0*adX[1]+m*s*(1-n-r)/4.0*adX[2]
			+m*r*s/2.0*adX[3]+m*r*(1-n-s)/4.0*adX[4]
			+m*n*r/2.0*adX[5]+n*r*(1-m-s)/4.0*adX[6]
			+n*r*s/2.0*adX[7];
		//coordinate_y=Ni*Yi
		adGaussY[loop]=n*s*(1-m-r)/4.0*adY[0]
			+m*n*s/2.0*adY[1]+m*s*(1-n-r)/4.0*adY[2]
			+m*r*s/2.0*adY[3]+m*r*(1-n-s)/4.0*adY[4]
			+m*n*r/2.0*adY[5]+n*r*(1-m-s)/4.0*adY[6]
			+n*r*s/2.0*adY[7];	
	}
}

void CPlane8NodeEle::GetMainStress(double *adMainStress, double *adMainStress1,
								   double *adMainStressAngle)
{
	for(int loop=0;loop<9;loop++){
		adMainStress[loop]=m_adMainStress[loop];
		adMainStress1[loop]=m_adMainStress1[loop];
		adMainStressAngle[loop]=m_adMainStressAngle[loop];
	}
}

void CPlane8NodeEle::GetMassMatrix()
{
	int loop;
	double dEta,dKxi,dBuf1,dBuf2,dBuf3,dBuf4;
	double dDensity,dJacobi,dWeight;
	double m,n,r,s;
	double adX[8],adY[8];
	CTypedPtrArray<CPtrArray,CBaseMaterial*>& apMaterial=*m_papMaterial;

	for(loop=0;loop<8;loop++){
		adX[loop]=m_pNode->GetX(m_aiNode[loop]);	
		adY[loop]=m_pNode->GetY(m_aiNode[loop]);
	}

	dDensity=apMaterial[m_iMaterialIndex]->GetDensity();
	m_matKe=0.0;
	for(loop=0;loop<9;loop++){
		dEta=m_adGaussEta[loop];
		dKxi=m_adGaussKxi[loop];

		m=1+dEta;n=1-dEta;
		//dN/ddKxi*Xi++;
		dBuf1=(-n/4.0+dKxi*n/2.0+m*n/4.0)*adX[0]
			-dKxi*n*adX[1]
			+(n/4.0+dKxi*n/2.0-m*n/4.0)*adX[2]
			+m*n/2.0*adX[3]
			+(m/4.0+dKxi*m/2.0-m*n/4.0)*adX[4]
			-dKxi*m*adX[5]
			+(-m/4.0+dKxi*m/2.0+m*n/4.0)*adX[6]
			-m*n/2.0*adX[7];
		
		//dN/ddKxi*Yi++;
		dBuf2=(-n/4.0+dKxi*n/2.0+m*n/4.0)*adY[0]
			-dKxi*n*adY[1]
			+(n/4.0+dKxi*n/2.0-m*n/4.0)*adY[2]
			+m*n/2.0*adY[3]
			+(m/4.0+dKxi*m/2.0-m*n/4.0)*adY[4]
			-dKxi*m*adY[5]
			+(-m/4.0+dKxi*m/2.0+m*n/4.0)*adY[6]
			-m*n/2.0*adY[7];
		
		m=1+dKxi;n=1-dKxi;	
		//dN/ddEta*Xi++;
		dBuf3=(-n/4.0+m*n/4.0+dEta*n/2.0)*adX[0]
			-m*n/2.0*adX[1]
			+(-m/4.0+m*n/4.0+dEta*m/2.0)*adX[2]
			-dEta*m*adX[3]
			+(m/4.0-m*n/4.0+dEta*m/2.0)*adX[4]
			+m*n/2.0*adX[5]
			+(n/4.0-m*n/4.0+dEta*n/2.0)*adX[6]
			-dEta*n*adX[7];
		
		//dN/ddEta*Yi++;
		dBuf4=(-n/4.0+m*n/4.0+dEta*n/2.0)*adY[0]
			-m*n/2.0*adY[1]
			+(-m/4.0+m*n/4.0+dEta*m/2.0)*adY[2]
			-dEta*m*adY[3]
			+(m/4.0-m*n/4.0+dEta*m/2.0)*adY[4]
			+m*n/2.0*adY[5]
			+(n/4.0-m*n/4.0+dEta*n/2.0)*adY[6]
			-dEta*n*adY[7];
		
		// |1   2| 
		// |3   4| 
		dJacobi=dBuf1*dBuf4-dBuf3*dBuf2;

		
		m=1.0+dKxi;n=1.0-dKxi;
		r=1.0+dEta;s=1.0-dEta;
		m_mat0216=0.0;
		m_mat0216(0,0)=m_mat0216(1,1)=n*s*(1-m-r)/4.0;
		m_mat0216(0,2)=m_mat0216(1,3)=m*n*s/2.0;
		m_mat0216(0,4)=m_mat0216(1,5)=m*s*(1-n-r)/4.0;
		m_mat0216(0,6)=m_mat0216(1,7)=m*r*s/2.0;
		m_mat0216(0,8)=m_mat0216(1,9)=m*r*(1-n-s)/4.0;
		m_mat0216(0,10)=m_mat0216(1,11)=m*n*r/2.0;
		m_mat0216(0,12)=m_mat0216(1,13)=n*r*(1-m-s)/4.0;
		m_mat0216(0,14)=m_mat0216(1,15)=n*r*s/2.0;

		m_mat1602.Trans(m_mat0216);

		if(dEta==0.0||dKxi==0.0){
			if(dKxi!=0.0||dEta!=0.0){
				dWeight=0.55555556*0.88888889;
			}
			else{
				dWeight=0.88888889*0.88888889;
			}
		}
		else{
			dWeight=0.55555556*0.55555556;
		}

		m_mat0216*=dDensity*m_dThickness*dWeight;
		m_matKe+=m_mat1602*m_mat0216;
	}
}

void CPlane8NodeEle::MassAssemble(CSparseMatrix& smatMass)
{
	int loop,loop1,iBuf,aiDOFIndex[16];
	//int nFreeDOF;
	GetMassMatrix();
	for(loop=0;loop<8;loop++){
		aiDOFIndex[2*loop]=m_pNode->GetXDOFIndex(m_aiNode[loop]);
		aiDOFIndex[2*loop+1]=m_pNode->GetYDOFIndex(m_aiNode[loop]);
	}
	//nFreeDOF=m_pNode->GetFreeDOF();
	for(loop=0;loop<16;loop++){
		//if(aiDOFIndex[loop]<nFreeDOF){
			iBuf=aiDOFIndex[loop];
			smatMass(iBuf,iBuf)+=m_matKe(loop,loop);
			for(loop1=0;loop1<loop;loop1++){
				//if(aiDOFIndex[loop1]<nFreeDOF){
					smatMass(iBuf,aiDOFIndex[loop1])+=m_matKe(loop,loop1);
				//}
			}
		//}
	}
}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -