📄 cele.cpp
字号:
#include "iostream.h"
#include "CEle.h"#include "memory.h"#include "math.h"CEle::CEle(){}CEle::~CEle(){}void CEle::GetN(double mx,double my){
//得到型函数的值,是关于x,y的二元函数 m_mN.Realloc(2,16); m_mN=0; m_mN(0,0)=m_mN(1,1)=0.25*(1-mx)*(1-my)-0.25*(1-mx*mx)*(1-my)-0.25*(1-mx)*(1-my*my); m_mN(0,2)=m_mN(1,3)=0.25*(1+mx)*(1-my)-0.25*(1-mx*mx)*(1-my)-0.25*(1+mx)*(1-my*my); m_mN(0,4)=m_mN(1,5)=0.25*(1+mx)*(1+my)-0.25*(1-mx*mx)*(1+my)-0.25*(1+mx)*(1-my*my); m_mN(0,6)=m_mN(1,7)=0.25*(1-mx)*(1+my)-0.25*(1-mx*mx)*(1+my)-0.25*(1-mx)*(1-my*my); m_mN(0,8)=m_mN(1,9)=0.5*(1-mx*mx)*(1-my); m_mN(0,10)=m_mN(1,11)=0.5*(1+mx)*(1-my*my); m_mN(0,12)=m_mN(1,13)=0.5*(1-mx*mx)*(1+my); m_mN(0,14)=m_mN(1,15)=0.5*(1-mx)*(1-my*my);}void CEle::GetJ(int IM,double mx,double my){
//等到IM单元的Jacob矩阵,
//在平面单元,为一2X2矩阵 m_mJ.Realloc(2,2); m_mJ=0; CMatrix mDN(2,8,0); mDN(0,0)=-0.25*(1-my)+0.5*mx*(1-my)+0.25*(1-my*my); mDN(1,0)=-0.25*(1-mx)+0.5*my*(1-mx)+0.25*(1-mx*mx); mDN(0,1)=0.25*(1-my)+0.5*mx*(1-my)-0.25*(1-my*my); mDN(1,1)=-0.25*(1+mx)+0.5*my*(1+mx)+0.25*(1-mx*mx); mDN(0,2)=0.25*(1+my)+0.5*mx*(1+my)-0.25*(1-my*my); mDN(1,2)=0.25*(1+mx)+0.5*my*(1+mx)-0.25*(1-mx*mx); mDN(0,3)=-0.25*(1+my)+0.5*mx*(1+my)+0.25*(1-my*my); mDN(1,3)=0.25*(1-mx)+0.5*my*(1-mx)-0.25*(1-mx*mx); mDN(0,4)=-1*mx*(1-my); mDN(1,4)=-0.5*(1-mx*mx); mDN(0,5)=0.5*(1-my*my); mDN(1,5)=-1*(1+mx)*my; mDN(0,6)=-1*mx*(1+my); mDN(1,6)=0.5*(1-mx*mx); mDN(0,7)=-0.5*(1-my*my); mDN(1,7)=-1*(1-mx)*my; CMatrix mX(8,2,0); for(int i=1;i<=8;i++) { int i1=int(m_mMP(IM-1,i-1)); mX(i-1,0)=m_mCOORD(i1-1,0); mX(i-1,1)=m_mCOORD(i1-1,1); } m_mJ=mDN*mX; m_dDetJ=m_mJ(0,0)*m_mJ(1,1)-m_mJ(1,0)*m_mJ(0,1); m_mJ(0,1)*=-1; m_mJ(1,0)*=-1; double m=m_mJ(0,0); m_mJ(0,0)=m_mJ(1,1); m_mJ(1,1)=m; m=1/m_dDetJ; m_mJ=m_mJ*m;}void CEle::GetB(int IM,double mx,double my){
//等到单元的B矩阵 m_mB.Realloc(3,16); m_mB=0; CMatrix mDN(4,16,0); mDN(0,0)=mDN(2,1)=-0.25*(1-my)+0.5*mx*(1-my)+0.25*(1-my*my); mDN(1,0)=mDN(3,1)=-0.25*(1-mx)+0.5*my*(1-mx)+0.25*(1-mx*mx); mDN(0,2)=mDN(2,3)=0.25*(1-my)+0.5*mx*(1-my)-0.25*(1-my*my); mDN(1,2)=mDN(3,3)=-0.25*(1+mx)+0.5*my*(1+mx)+0.25*(1-mx*mx); mDN(0,4)=mDN(2,5)=0.25*(1+my)+0.5*mx*(1+my)-0.25*(1-my*my); mDN(1,4)=mDN(3,5)=0.25*(1+mx)+0.5*my*(1+mx)-0.25*(1-mx*mx); mDN(0,6)=mDN(2,7)=-0.25*(1+my)+0.5*mx*(1+my)+0.25*(1-my*my); mDN(1,6)=mDN(3,7)=0.25*(1-mx)+0.5*my*(1-mx)-0.25*(1-mx*mx); mDN(0,8)=mDN(2,9)=-1*mx*(1-my); mDN(1,8)=mDN(3,9)=-0.5*(1-mx*mx); mDN(0,10)=mDN(2,11)=0.5*(1-my*my); mDN(1,10)=mDN(3,11)=-1*(1+mx)*my; mDN(0,12)=mDN(2,13)=-1*mx*(1+my); mDN(1,12)=mDN(3,13)=0.5*(1-mx*mx); mDN(0,14)=mDN(2,15)=-0.5*(1-my*my); mDN(1,14)=mDN(3,15)=-1*(1-mx)*my; GetJ(IM,mx,my); CMatrix mA(3,4,0); mA(0,0)=mA(2,2)=m_mJ(0,0); mA(0,1)=mA(2,3)=m_mJ(0,1); mA(1,2)=mA(2,0)=m_mJ(1,0); mA(2,1)=mA(1,3)=m_mJ(1,1); m_mB=mA*mDN;}void CEle::GetD(){
//得到单元的D矩阵,因为单元材料假设
//相同,故所有单元D阵相同 m_mD.Realloc(3,3); m_mD=0; m_mD(0,0)=m_mD(1,1)=1; m_mD(0,1)=m_mD(1,0)=m_dMu; m_mD(2,2)=0.5*(1-m_dMu); double d0=m_dE/(1-m_dMu*m_dMu); m_mD=m_mD*d0;}void CEle::GetKe(int IM){
//得到单元IM的单元刚度矩阵 m_mKe.Realloc(16,16); m_mKe=0; CMatrix mTranB(16,3,0); double x1=1/sqrt(3); double x2=-1/sqrt(3); GetB(IM,x1,x1); mTranB.Trans(m_mB); m_mKe=m_mKe+mTranB*m_mD*m_mB; GetB(IM,x1,x2); mTranB.Trans(m_mB); m_mKe=m_mKe+mTranB*m_mD*m_mB; GetB(IM,x2,x1); mTranB.Trans(m_mB); m_mKe=m_mKe+mTranB*m_mD*m_mB; GetB(IM,x2,x2); mTranB.Trans(m_mB); m_mKe=m_mKe+mTranB*m_mD*m_mB; m_mKe=m_mKe*m_dDetJ*m_db;}void CEle::GetPe(int IM){
//得到单元荷载矩阵,
//在本程序中,为一二维矩阵,第二维元素为1 m_mPe.Realloc(16,1); m_mPe=0; int i; for(i=1;i<=m_nN1;i++) if(IM-int(m_mML(i-1,0))==0) { CMatrix mNT(16,2,0); CMatrix mQ(2,1,0); mQ(1,0)=m_mML(i-1,1); double x1=1/sqrt(3.0); double x2=-1/sqrt(3.0); CMatrix mA1(16,1,0); GetN(x1,1); GetJ(IM,x1,1); mA1=mNT.Trans(m_mN)*mQ; double x3=m_mJ(1,1)/(m_mJ(0,0)*m_mJ(1,1)-m_mJ(1,0)*m_mJ(0,1)); mA1=mA1*x3; CMatrix mA2(16,1,0); GetN(x2,1); GetJ(IM,x2,1); mA2=mNT.Trans(m_mN)*mQ; x3=m_mJ(1,1)/(m_mJ(0,0)*m_mJ(1,1)-m_mJ(1,0)*m_mJ(0,1)); mA2=mA2*x3; m_mPe=mA1+mA2; }}void CEle::GetK(){
//将所有单元的刚度矩阵集成
//到整体刚度矩阵
//同时集成整体荷载向量 m_mK.Realloc(2*m_nNJ,2*m_nNJ); m_mK=0; m_mP.Realloc(2*m_nNJ,1); m_mP=0; int IM; int i,j,i1,j1; for(IM=1;IM<=m_nNM;IM++) { GetPe(IM); GetKe(IM); for(i=1;i<=8;i++) { i1=int(m_mMP(IM-1,i-1));
m_mP(2*i1-2,0)+=m_mPe(2*i-2,0);
m_mP(2*i1-1,0)+=m_mPe(2*i-1,0); for(j=1;j<=8;j++) { j1=m_mMP(IM-1,j-1); m_mK(2*i1-2,2*j1-2)+=m_mKe(2*i-2,2*j-2); m_mK(2*i1-1,2*j1-2)+=m_mKe(2*i-1,2*j-2); m_mK(2*i1-2,2*j1-1)+=m_mKe(2*i-2,2*j-1); m_mK(2*i1-1,2*j1-1)+=m_mKe(2*i-1,2*j-1); } } }}void CEle::Solve(){
//用高斯消元法解方程 m_mDis.Realloc(2*m_nNJ,1); m_mDis=0; int i,j,m;
double z1,z; for(m=0;m<2*m_nNJ-1;m++)
for(i=m+1;i<2*m_nNJ;i++)
{
z=m_mK(m,i)/m_mK(m,m);
m_mP(i,0)-=z*m_mP(m,0);
for(j=i;j<2*m_nNJ;j++)
m_mK(i,j)-=m_mK(m,j)*z;
}
for(i=2*m_nNJ-1;i>=0;i--)
{
z=0;
for(j=2*m_nNJ-1;j>i;j--)
z+=m_mDis(j,0)*m_mK(i,j);
m_mDis(i,0)=(m_mP(i,0)-z)/m_mK(i,i);
}
m_mP=m_mDis;}//LU:solve the equations P=SDvoid CEle::GetDise(int IM){
//得到单元IM的局部位移向量 m_mDise.Realloc(16,1); m_mDise=0; int i,i1; for(i=1;i<=8;i++) { i1=m_mMP(IM-1,i-1); m_mDise(2*i-2,0)=m_mP(2*i1-2,0); m_mDise(2*i-1,0)=m_mP(2*i1-1,0); }}void CEle::GetSe(int IM){ //单元应力磨平,用高斯外推法 m_mSe.Realloc(3,8); m_mSe=0; CMatrix mS[8]; int i,j,i1; for(i=1;i<=8;i++) mS[i-1].Realloc(3,1); double x1=1/sqrt(3); double x2=-1/sqrt(3); GetB(IM,x2,x2); mS[0]=m_mD*m_mB*m_mDise; GetB(IM,x1,x2); mS[1]=m_mD*m_mB*m_mDise; GetB(IM,x1,x1); mS[2]=m_mD*m_mB*m_mDise; GetB(IM,x2,x1); mS[3]=m_mD*m_mB*m_mDise; CMatrix mN(4,4,0); mN(0,0)=mN(1,1)=mN(2,2)=mN(3,3)=1+sqrt(3)/2; mN(0,1)=mN(0,3)=mN(1,0)=mN(1,2)=mN(2,1)=mN(2,3)=mN(3,0)=mN(3,2)=-0.5; mN(0,2)=mN(1,3)=mN(2,0)=mN(3,1)=1-sqrt(3)/2; CMatrix mS1(4,3,0); for(i=1;i<=4;i++) { mS1(i-1,0)=mS[i-1](0,0); mS1(i-1,1)=mS[i-1](1,0); mS1(i-1,2)=mS[i-1](2,0); } CMatrix mS2(4,3,0); mS2=mN*mS1; CMatrix mS3(3,4,0); mS3.Trans(mS2); GetN2(0,-1); mS[4]=mS3*m_mN2; GetN2(1,0); mS[5]=mS3*m_mN2; GetN2(0,1); mS[6]=mS3*m_mN2; GetN2(-1,0); mS[7]=mS3*m_mN2; for(i=1;i<=3;i++) for(j=1;j<=4;j++) { m_mSe(i-1,j-1)=mS3(i-1,j-1); m_mSe(i-1,j-1+4)=mS[j-1+4](i-1,0); }}void CEle::GetN2(double mx,double my){
//用四节点的型函数
//在单元应力磨平时用到 m_mN2.Realloc(4,1); m_mN2(0,0)=0.25*(1-mx)*(1-my); m_mN2(1,0)=0.25*(1+mx)*(1-my); m_mN2(2,0)=0.25*(1+mx)*(1+my); m_mN2(3,0)=0.25*(1-mx)*(1+my);}void CEle::GetS(){
//等到整体应力矩阵
//即,将单元应力集成到整体阵中 m_mS.Realloc(3,m_nNJ); m_mS=0; int i,j,i1; int *I=new int [m_nNJ]; for(i=1;i<=m_nNJ;i++) I[i-1]=0; int IM; for(IM=1;IM<=m_nNM;IM++) { GetDise(IM); GetSe(IM); for(i=1;i<=8;i++) { i1=m_mMP(IM-1,i-1); m_mS(0,i1-1)+=m_mSe(0,i-1); m_mS(1,i1-1)+=m_mSe(1,i-1); m_mS(2,i1-1)+=m_mSe(2,i-1); I[i1-1]++; } } for(i=1;i<=m_nNJ;i++) for(j=1;j<=3;j++) m_mS(j-1,i-1)/=double(I[i-1]);} CEle::CEle(double dL,double dB,double dH,double dQ,int nN1,int nN2,double dE,double dMu){
//构造函数,并得到梁的
//全部编号信息
// m_db=dB; m_dL=dL; m_dH=dH; m_dQ=dQ; m_nN1=nN1; m_nN2=nN2; m_dE=dE; m_dMu=dMu; int i,j,k,l; //建立节点信息 m_nNJ=(m_nN1+1)*(2*m_nN2+1)+m_nN1*(m_nN2+1); m_mCOORD.Realloc(m_nNJ,2); k=1; for(i=1;i<=2*m_nN1+1;i++) { if(i%2==1) for(j=1;j<=2*m_nN2+1;j++) { m_mCOORD(k-1,0)=m_dL*(i-1)/(2*m_nN1); m_mCOORD(k-1,1)=-1*m_dH*(j-1)/(2*m_nN2); k++; } else if(i%2==0) for(j=1;j<=m_nN2+1;j++) { m_mCOORD(k-1,0)=m_dL*(i-1)/(2*m_nN1); m_mCOORD(k-1,1)=-1*m_dH*(j-1)/m_nN2; k++; } }
//建立单元信息 m_nNM=m_nN1*m_nN2; m_mMP.Realloc(m_nNM,8); k=1; for(i=1;i<=m_nN1;i++) { for(j=1;j<=m_nN2;j++) { l=(3*m_nN2+2)*(i-1)+2*(j-1)+1; m_mMP(k-1,0)=l+2; m_mMP(k-1,1)=l+3*m_nN2+4; m_mMP(k-1,2)=l+3*m_nN2+2; m_mMP(k-1,3)=l; m_mMP(k-1,4)=l+2*m_nN2+1-(j-1)+1; m_mMP(k-1,5)=l+3*m_nN2+2+1; m_mMP(k-1,6)=l+2*m_nN2+1-(j-1); m_mMP(k-1,7)=l+1; k++; } } //建立荷载信息 m_mML.Realloc(m_nN1,2); m_mML=0; for(i=1;i<=m_nN1;i++) { m_mML(i-1,0)=m_nN2*(i-1)+1; m_mML(i-1,1)=m_dQ; }}void CEle::restrict(char flag1,char flag2)
{
int i=0;
switch(flag1)
{
case 'g':
case 'G':
for(i=0;i<=2*m_nN2;i++)
m_mK(2*i,2*i)+=1e10,m_mK(2*i+1,2*i+1)+=1e10;
break;
case 'j':
case 'J':
m_mK(2*m_nN2,2*m_nN2)+=1e10;
m_mK(2*m_nN2+1,2*m_nN2+1)+=1e10;
break;
case 'l':
case 'L':
m_mK(2*m_nN2+1,2*m_nN2+1)+=1e10;
break;
case 'h':
case 'H':
for(i=0;i<=2*m_nN2;i++)
m_mK(2*i,2*i)+=1e10;
break;
}
switch(flag2)
{
case 'g':
case 'G':
for(i=m_nNJ-1;i>=m_nNJ-2*m_nN2-1;i--)
{
m_mK(2*i,2*i)+=1e10;
m_mK(2*i+1,2*i+1)+=1e10;
}
break;
case 'j':
case 'J':
m_mK(2*m_nNJ-2,2*m_nNJ-2)+=1e10;
m_mK(2*m_nNJ-1,2*m_nNJ-1)+=1e10;
break;
case 'l':
case 'L':
m_mK(2*m_nNJ-1,2*m_nN2-1)+=1e10;
break;
case 'h':
case 'H':
for(i=m_nNJ-1;i>=m_nNJ-2*m_nN2-1;i--)
m_mK(2*i,2*i)+=1e10;
break;
}
// m_mP(2*m_nN2,2*m_nN2)+=1e10;
// m_mP(2*m_nN2+1,2*m_nN2+1)+=1e10;
// m_mP(2*m_nNJ,2*m_nNJ)+=1e10;
// m_mP(2*m_nNJ+1,2*m_nNJ+1)+=1e10;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -