📄 rkpmdlg.cpp
字号:
UpdateData();
double r=m_banwidth/(m_vdotnumber-1);
static double a[2]={-0.57735,0.57735};
//计算边界高斯积分点坐标gyp
for(int ki=0;ki<m_vdotnumber;ki++)
for(int kt=0;kt<2;kt++)
bgy[2*ki+kt]=((1-a[kt])/2)*ki*r+((1+a[kt])/2)*(ki+1)*r;
}
mxArray * CRKPMDlg::setBLM(double bgy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double BLMelement[9];
mxArray *BLM=NULL;//申明左边界M矩阵
mlfAssign(&BLM,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算BLM矩阵,矩阵阶数3x3
for(int kj=0;kj<NP;kj++)//总体坐标循环
{
//计算BLM矩阵,参照公式推导(4.7)式
setkernal(0,bgy,x[kj],y[kj]);
BLMelement[0]=ker;//BLM矩阵M11
BLMelement[1]=(0-x[kj])*ker;//BLMM矩阵M21
BLMelement[2]=(bgy-y[kj])*ker;//BLM矩阵M31
BLMelement[3]=(0-x[kj])*ker;//BLM矩阵M12
BLMelement[4]=(0-x[kj])*(0-x[kj])*ker;//BLM矩阵M22
BLMelement[5]=(0-x[kj])*(bgy-y[kj])*ker;//BLM矩阵M32
BLMelement[6]=(bgy-y[kj])*ker;//BLM矩阵M13
BLMelement[7]=(0-x[kj])*(bgy-y[kj])*ker;//BLM矩阵M23
BLMelement[8]=(bgy-y[kj])*(bgy-y[kj])*ker;//BLM矩阵M33
mxArray *BLLM=NULL;//申明临时矩阵,以便求M矩阵
mlfAssign(&BLLM,mlfDoubleMatrix(3,3,BLMelement,NULL));//给L阵赋值
mlfAssign(&BLM,mlfPlus(BLM,BLLM));
mxDestroyArray(BLLM);//释放临时矩阵内存
}
return BLM;
}
mxArray * CRKPMDlg::setBRM(double bgy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double BRMelement[9];
mxArray *BRM=NULL;//申明右边界M矩阵
mlfAssign(&BRM,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算BRM矩阵及其导数阵,矩阵阶数3x3
for(int kj=0;kj<NP;kj++)//总体坐标循环
{
//计算BRM矩阵,参照公式推导(4.7)式
setkernal(m_banlength,bgy,x[kj],y[kj]);
BRMelement[0]=ker;//BLM矩阵M11
BRMelement[1]=(m_banlength-x[kj])*ker;//BLMM矩阵M21
BRMelement[2]=(bgy-y[kj])*ker;//BLM矩阵M31
BRMelement[3]=(m_banlength-x[kj])*ker;//BLM矩阵M12
BRMelement[4]=(m_banlength-x[kj])*(m_banlength-x[kj])*ker;//BLM矩阵M22
BRMelement[5]=(m_banlength-x[kj])*(bgy-y[kj])*ker;//BLM矩阵M32
BRMelement[6]=(bgy-y[kj])*ker;//BLM矩阵M13
BRMelement[7]=(m_banlength-x[kj])*(bgy-y[kj])*ker;//BLM矩阵M23
BRMelement[8]=(bgy-y[kj])*(bgy-y[kj])*ker;//BLM矩阵M33
mxArray *BLRM=NULL;//申明临时矩阵,以便求M矩阵
mlfAssign(&BLRM,mlfDoubleMatrix(3,3,BRMelement,NULL));//给L阵赋值
mlfAssign(&BRM,mlfPlus(BRM,BLRM));
mxDestroyArray(BLRM);//释放临时矩阵内存
}
return BRM;
}
mxArray * CRKPMDlg::setBLN(double bgy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *BLMI=NULL;//声明BLM阵的逆阵
mlfAssign(&BLMI,mlfInv(setBLM(bgy)));
mxArray *BLN=NULL;//声明左边界型函数矩阵BN
mlfAssign(&BLN,mlfDoubleMatrix(1,NP,NULL,NULL));
double h0[3]={1,0,0};
for(int kn=0;kn<NP;kn++)
{
mxArray *BH0=NULL;
mlfAssign(&BH0,mlfDoubleMatrix(1,3,h0,NULL));
//计算左边界型函数矩阵BN
setkernal(0,bgy,x[kn],y[kn]);
setbasicfun(0,bgy,x[kn],y[kn]);
mxArray *BLHX=NULL;//申明基函数阵
mlfAssign(&BLHX,mlfDoubleMatrix(3,1,H,NULL));
mlfIndexAssign(&BLN,"(?,?)",mlfScalar(1),mlfScalar(kn+1),mlfMtimes(BH0,
mlfMtimes(BLMI,mlfMtimes(BLHX,mlfScalar(ker)))));
mxDestroyArray(BH0);
mxDestroyArray(BLHX);
}
return BLN;
}
mxArray * CRKPMDlg::setBRN(double bgy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *BRMI=NULL;//声明BLM阵的逆阵
mlfAssign(&BRMI,mlfInv(setBRM(bgy)));
mxArray *BRN=NULL;//声明右边界型函数矩阵BN
mlfAssign(&BRN,mlfDoubleMatrix(1,NP,NULL,NULL));
double h0[3]={1,0,0};
for(int kn=0;kn<NP;kn++)
{
mxArray *BH0=NULL;
mlfAssign(&BH0,mlfDoubleMatrix(1,3,h0,NULL));
//计算右边界型函数矩阵BN
setkernal(m_banlength,bgy,x[kn],y[kn]);
setbasicfun(m_banlength,bgy,x[kn],y[kn]);
mxArray *BRHX=NULL;//申明基函数阵
mlfAssign(&BRHX,mlfDoubleMatrix(3,1,H,NULL));
mlfIndexAssign(&BRN,"(?,?)",mlfScalar(1),mlfScalar(kn+1),mlfMtimes(BH0,
mlfMtimes(BRMI,mlfMtimes(BRHX,mlfScalar(ker)))));
mxDestroyArray(BH0);
mxDestroyArray(BRHX);
}
return BRN;
}
void CRKPMDlg::setBya()
{
UpdateData();
double r=m_banwidth/(m_vdotnumber-1);
Bya=r/2;//边界雅克比常数
}
mxArray *CRKPMDlg::setDe()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
//计算弹性关系矩阵De,参照公式推导(2.6)式
mxArray *De=NULL;
double d[9];
double d0;
d0=(m_e*pow(m_banheigth,3))/(12*(1-pow(m_u,2)));
d[0]=d0;
d[1]=d0*m_u;
d[2]=0;
d[3]=d0*m_u;
d[4]=d0;
d[5]=0;
d[6]=0;
d[7]=0;
d[8]=d0*(1-m_u)/2;
mlfAssign(&De,mlfDoubleMatrix(3,3,d,NULL));
mlfSave(mxCreateString("D:\\Paper\\Program\\RKPM\\Result\\de.mat"),"w","de",De,NULL);
return De;
}
mxArray *CRKPMDlg::setD()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
//计算弹性关系矩阵D,参照公式推导(2.6)式
mxArray *D=NULL;
double d[9];
double d0;
d0=(m_e*m_banheigth/2)/(12*(1-pow(m_u,2)));
d[0]=d0;
d[1]=d0*m_u;
d[2]=0;
d[3]=d0*m_u;
d[4]=d0;
d[5]=0;
d[6]=0;
d[7]=0;
d[8]=d0*(1-m_u)/2;
mlfAssign(&D,mlfDoubleMatrix(3,3,d,NULL));
mlfSave(mxCreateString("D:\\Paper\\Program\\RKPM\\Result\\d.mat"),"w","d",D,NULL);
return D;
}
mxArray *CRKPMDlg::setKe()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *Ke=NULL;
mlfAssign(&Ke,mlfDoubleMatrix(NP,NP,NULL,NULL));
mxArray *De=NULL;
mlfAssign(&De,setDe());
for(int i=0;i<4*(m_hdotnumber-1)*(m_vdotnumber-1);i++)//高斯积分点循环
{
//计算几何阵,参照公式推导(4.13),(4.14),(4.15)式
mxArray *B=NULL;//几何阵
mlfAssign(&B,setB(gx[i],gy[i]));
//计算刚度矩阵Ke,参照公式推导(5.10)式
mxArray *KL=NULL;//声明临时矩阵以便计算刚度矩阵K
mlfAssign(&KL,mlfMtimes(mlfScalar(HI),mlfMtimes(mlfScalar(HJ),
mlfMtimes(mlfTranspose(B),mlfMtimes(De,mlfMtimes(B,mlfScalar(ya[i])))))));
mlfAssign(&Ke,mlfPlus(Ke,KL));
mxDestroyArray(KL);
mxDestroyArray(B);
}
return Ke;
}
mxArray *CRKPMDlg::setF(double m_f)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *F=NULL;
mlfAssign(&F,mlfDoubleMatrix(NP,1,NULL,NULL));
for(int i=0;i<4*(m_hdotnumber-1)*(m_vdotnumber-1);i++)//高斯积分点循环
{
//计算型函数矩阵N,参照公式推导(1.17),(5.11)式
mxArray *N=NULL;//声明型函数矩阵N
mlfAssign(&N,setN(gx[i],gy[i]));
//计算载荷阵F,参照公式推导(5.12)式
mxArray *FL=NULL;//声明临时矩阵以便计算载荷阵F
mlfAssign(&FL,mlfMtimes(mlfScalar(HI),mlfMtimes(mlfScalar(HJ),
mlfMtimes(mlfTranspose(N),mlfMtimes(mlfScalar(m_f),mlfScalar(10000*ya[i]))))));
mlfAssign(&F,mlfPlus(F,FL));
mxDestroyArray(FL);
mxDestroyArray(N);
}
return F;
}
mxArray *CRKPMDlg::setLKu()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *LKu=NULL;
mlfAssign(&LKu,mlfDoubleMatrix(NP,NP,NULL,NULL));
//计算罚函数刚度矩阵LKu,参照公式推导(5.10)式
for(int j=0;j<2*(m_vdotnumber-1);j++)
{
mxArray *BLN=NULL;//声明左边界型函数矩阵BN
mlfAssign(&BLN,setBLN(bgy[j]));
//计算左罚函数阵
mxArray *LLKu=NULL;
mlfAssign(&LLKu,mlfMtimes(mlfScalar(HJ),mlfMtimes(mlfTranspose(BLN),
mlfMtimes(BLN,mlfMtimes(mlfScalar(m_alpha),mlfScalar(Bya))))));
mlfAssign(&LKu,mlfPlus(LKu,LLKu));
mxDestroyArray(LLKu);
mxDestroyArray(BLN);
}
return LKu;
}
mxArray *CRKPMDlg::setRKu()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *RKu=NULL;
mlfAssign(&RKu,mlfDoubleMatrix(NP,NP,NULL,NULL));
//计算罚函数刚度矩阵RKu,参照公式推导(5.10)式
for(int j=0;j<2*(m_vdotnumber-1);j++)
{
mxArray *BRN=NULL;//声明右边界型函数矩阵BN
mlfAssign(&BRN,setBRN(bgy[j]));
//计算右罚函数阵
mxArray *LRKu=NULL;
mlfAssign(&LRKu,mlfMtimes(mlfScalar(HJ),mlfMtimes(mlfTranspose(BRN),
mlfMtimes(BRN,mlfMtimes(mlfScalar(m_alpha),mlfScalar(Bya))))));
mlfAssign(&RKu,mlfPlus(RKu,LRKu));
mxDestroyArray(LRKu);
mxDestroyArray(BRN);
}
return RKu;
}
void CRKPMDlg::calculate()
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
//计算弹性关系矩阵D,参照公式推导(2.6)式
mxArray *De=NULL;
mlfAssign(&De,setDe());
mxArray *D=NULL;
mlfAssign(&D,setD());
mxArray *Ke=NULL;//声明刚度矩阵K
mlfAssign(&Ke,setKe());//给刚度矩阵置初值
mxArray *LKu=NULL;//申明左罚函数刚度矩阵Ku
mlfAssign(&LKu,setLKu());
mxArray *RKu=NULL;//申明左罚函数刚度矩阵Ku
mlfAssign(&RKu,setRKu());
mxArray *Ku=NULL;//声明罚函数刚度矩阵Ku
mlfAssign(&Ku,mlfPlus(LKu,RKu));
mxArray *KG=NULL;//声明总体刚度矩阵KG=Ke+Ku
mlfAssign(&KG,mlfPlus(Ke,Ku));
mxArray *F=NULL;//声明载荷阵F
mlfAssign(&F,setF(m_f));//给载荷阵置初值
mxArray *UI=NULL;//声明位移阵
mlfAssign(&UI,mlfMldivide(KG,F));
//计算应变，应力
mxArray *S=NULL;//应变量
mlfAssign(&S,mlfDoubleMatrix(3*NP,1,NULL,NULL));
mxArray *T=NULL;//应力变量
mlfAssign(&T,mlfDoubleMatrix(3*NP,1,NULL,NULL));
for(int k=0;k<NP;k++)
{
mxArray *B=NULL;
mlfAssign(&B,setB(x[k],y[k]));
mxArray *SI=NULL;//节点应变量
mlfAssign(&SI,mlfMtimes(B,UI));
mlfIndexAssign(&S,"(?,?)",mlfScalar(3*k+1),mlfScalar(1),
mlfIndexRef(SI,"(?,?)",mlfScalar(1),mlfScalar(1),NULL));
mlfIndexAssign(&S,"(?,?)",mlfScalar(3*k+2),mlfScalar(1),
mlfIndexRef(SI,"(?,?)",mlfScalar(2),mlfScalar(1),NULL));
mlfIndexAssign(&S,"(?,?)",mlfScalar(3*k+3),mlfScalar(1),
mlfIndexRef(SI,"(?,?)",mlfScalar(3),mlfScalar(1),NULL));
mxArray *TI=NULL;
mlfAssign(&TI,mlfMtimes(D,SI));
mlfIndexAssign(&T,"(?,?)",mlfScalar(3*k+1),mlfScalar(1),
mlfIndexRef(TI,"(?,?)",mlfScalar(1),mlfScalar(1),NULL));
mlfIndexAssign(&T,"(?,?)",mlfScalar(3*k+2),mlfScalar(1),
mlfIndexRef(TI,"(?,?)",mlfScalar(2),mlfScalar(1),NULL));
mlfIndexAssign(&T,"(?,?)",mlfScalar(3*k+3),mlfScalar(1),
mlfIndexRef(TI,"(?,?)",mlfScalar(3),mlfScalar(1),NULL));
}
mlfSave(mxCreateString("D:\\Paper\\Program\\RKPM\\Result\\result.mat"),
"w","ke",Ke,"ku",Ku,"kg",KG,"f",F,"ui",UI,"s",S,"t",T,NULL);
mlfRestorePreviousContext(0,0);
}
void CRKPMDlg::OnOK()
{
// TODO: Add extra validation here
UpdateData();
setcordinate();
setgaosicordinate();
setBgaosicordinate();
setya();
setBya();
calculate();
outputcordinate();
CDialog::OnOK();
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -