📄 rkpmdlg.cpp
字号:
((-1+b[nk])*y[ni*m_hdotnumber+nj]+(1-b[nk])*y[ni*m_hdotnumber+nj+1]+(1+b[nk])*y[(ni+1)*m_hdotnumber+nj+1]+(-1-b[nk])*y[(ni+1)*m_hdotnumber+nj])*
((-1+a[nk])*x[ni*m_hdotnumber+nj]+(-1-a[nk])*x[ni*m_hdotnumber+nj+1]+(1+a[nk])*x[(ni+1)*m_hdotnumber+nj+1]+(1-a[nk])*x[(ni+1)*m_hdotnumber+nj]))/4;
}
ofstream outfile2;
outfile2.open("D:\\Paper\\Program\\RKPM\\Result\\ya.txt");
outfile2<<setiosflags(ios::left);
for(int yi=0;yi<(4*(m_hdotnumber-1)*(m_vdotnumber-1));yi++)
outfile2<<setw(20)<<yi<<setw(20)<<ya[yi]<<endl;
outfile2.close();
}
void CRKPMDlg::setbasicfun(double a,double b,double c,double d)
{
H[0]=1;
H[1]=a-c;
H[2]=b-d;
}
mxArray * CRKPMDlg::setM(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melement[9];//M矩阵元素
mxArray *M=NULL;//申明M矩阵
mlfAssign(&M,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M矩阵,矩阵阶数3x3
for(int j=0;j<NP;j++)//总体坐标循环
{
//计算M矩阵,参照公式推导(4.7)式
setkernal(gx,gy,x[j],y[j]);
Melement[0]=ker;//M矩阵M11
Melement[1]=(gx-x[j])*ker;//M矩阵M21
Melement[2]=(gy-y[j])*ker;//M矩阵M31
Melement[3]=(gx-x[j])*ker;//M矩阵M12
Melement[4]=(gx-x[j])*(gx-x[j])*ker;//M矩阵M22
Melement[5]=(gx-x[j])*(gy-y[j])*ker;//M矩阵M32
Melement[6]=(gy-y[j])*ker;//M矩阵M13
Melement[7]=(gx-x[j])*(gy-y[j])*ker;//M矩阵M23
Melement[8]=(gy-y[j])*(gy-y[j])*ker;//M矩阵M33
mxArray *L=NULL;//申明临时矩阵,以便求M矩阵
mlfAssign(&L,mlfDoubleMatrix(3,3,Melement,NULL));//给L阵赋值
mlfAssign(&M,mlfPlus(M,L));
mxDestroyArray(L);//释放临时矩阵内存
}
return M;
}
mxArray * CRKPMDlg::setMDX(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melementofdx[9];//M矩阵元素对x的导数
mxArray *MDX=NULL;//申明M阵对x的导数阵
mlfAssign(&MDX,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M阵对x的导数,参照公式推导(4.8)式
for(int j=0;j<NP;j++)//总体坐标循环
{
setkernal(gx,gy,x[j],y[j]);
setkernalofdx(gx,gy,x[j],y[j]);
Melementofdx[0]=kerofdx;
Melementofdx[1]=ker+(gx-x[j])*kerofdx;
Melementofdx[2]=(gy-y[j])*kerofdx;
Melementofdx[3]=ker+(gx-x[j])*kerofdx;
Melementofdx[4]=2*(gx-x[j])*ker+(gx-x[j])*(gx-x[j])*kerofdx;
Melementofdx[5]=(gy-y[j])*ker+(gx-x[j])*(gy-y[j])*kerofdx;
Melementofdx[6]=(gy-y[j])*kerofdx;
Melementofdx[7]=(gy-y[j])*ker+(gx-x[j])*(gy-y[j])*kerofdx;
Melementofdx[8]=(gy-y[j])*(gy-y[j])*kerofdx;
mxArray *LX=NULL;
mlfAssign(&LX,mlfDoubleMatrix(3,3,Melementofdx,NULL));
mlfAssign(&MDX,mlfPlus(MDX,LX));
mxDestroyArray(LX);
}
return MDX;
}
mxArray * CRKPMDlg::setMD2X(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melementofd2x[9];//M矩阵元素对x的导数
mxArray *MD2X=NULL;//申明M阵对x的导数阵
mlfAssign(&MD2X,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M阵对x的2阶导数,参照公式推导(4.9)式
for(int j=0;j<NP;j++)//总体坐标循环
{
setkernal(gx,gy,x[j],y[j]);
setkernalofdx(gx,gy,x[j],y[j]);
setkernalofd2x(gx,gy,x[j],y[j]);
Melementofd2x[0]=kerofd2x;
Melementofd2x[1]=2*kerofdx+(gx-x[j])*kerofd2x;
Melementofd2x[2]=(gy-y[j])*kerofd2x;
Melementofd2x[3]=2*kerofdx+(gx-x[j])*kerofd2x;
Melementofd2x[4]=2*ker+4*(gx-x[j])*kerofdx+(gx-x[j])*(gx-x[j])*kerofd2x;
Melementofd2x[5]=2*(gy-y[j])*kerofdx+(gx-x[j])*(gy-y[j])*kerofd2x;
Melementofd2x[6]=(gy-y[j])*kerofd2x;
Melementofd2x[7]=2*(gy-y[j])*kerofdx+(gx-x[j])*(gy-y[j])*kerofd2x;
Melementofd2x[8]=(gy-y[j])*(gy-y[j])*kerofd2x;
mxArray *L2X=NULL;
mlfAssign(&L2X,mlfDoubleMatrix(3,3,Melementofd2x,NULL));
mlfAssign(&MD2X,mlfPlus(MD2X,L2X));
mxDestroyArray(L2X);
}
return MD2X;
}
mxArray * CRKPMDlg::setMDY(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melementofdy[9];//M矩阵元素对x的导数
mxArray *MDY=NULL;//申明M阵对x的导数阵
mlfAssign(&MDY,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M阵对y的导数,参照公式推导(4.10)式
for(int j=0;j<NP;j++)//总体坐标循环
{
setkernal(gx,gy,x[j],y[j]);
setkernalofdy(gx,gy,x[j],y[j]);
Melementofdy[0]=kerofdy;
Melementofdy[1]=(gx-x[j])*kerofdy;
Melementofdy[2]=ker+(gy-y[j])*kerofdy;
Melementofdy[3]=(gx-x[j])*kerofdy;
Melementofdy[4]=(gx-x[j])*(gx-x[j])*kerofdy;
Melementofdy[5]=(gx-x[j])*ker+(gx-x[j])*(gy-y[j])*kerofdy;
Melementofdy[6]=ker+(gy-y[j])*kerofdy;
Melementofdy[7]=(gx-x[j])*ker+(gx-x[j])*(gy-y[j])*kerofdy;
Melementofdy[8]=2*(gy-y[j])*ker+(gy-y[j])*(gy-y[j])*kerofdy;
mxArray *LY=NULL;
mlfAssign(&LY,mlfDoubleMatrix(3,3,Melementofdy,NULL));
mlfAssign(&MDY,mlfPlus(MDY,LY));
mxDestroyArray(LY);
}
return MDY;
}
mxArray * CRKPMDlg::setMD2Y(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melementofd2y[9];//M矩阵元素对x的导数
mxArray *MD2Y=NULL;//申明M阵对x的导数阵
mlfAssign(&MD2Y,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M阵对y的2阶导数,参照公式推导(4.11)式
for(int j=0;j<NP;j++)//总体坐标循环
{
setkernal(gx,gy,x[j],y[j]);
setkernalofdy(gx,gy,x[j],y[j]);
setkernalofd2y(gx,gy,x[j],y[j]);
Melementofd2y[0]=kerofd2y;
Melementofd2y[1]=(gx-x[j])*kerofd2y;
Melementofd2y[2]=2*kerofdy+(gy-y[j])*kerofd2y;
Melementofd2y[3]=(gx-x[j])*kerofd2y;
Melementofd2y[4]=(gx-x[j])*(gx-x[j])*kerofd2y;
Melementofd2y[5]=2*(gx-x[j])*kerofdy+(gx-x[j])*(gy-y[j])*kerofd2y;
Melementofd2y[6]=2*kerofdy+(gy-y[j])*kerofd2y;
Melementofd2y[7]=2*(gx-x[j])*kerofdy+(gx-x[j])*(gy-y[j])*kerofd2y;
Melementofd2y[8]=2*ker+4*(gy-y[j])*kerofdy+(gy-y[j])*(gy-y[j])*kerofd2y;
mxArray *L2Y=NULL;
mlfAssign(&L2Y,mlfDoubleMatrix(3,3,Melementofd2y,NULL));
mlfAssign(&MD2Y,mlfPlus(MD2Y,L2Y));
mxDestroyArray(L2Y);
}
return MD2Y;
}
mxArray * CRKPMDlg::setMDXDY(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
double Melementofdxdy[9];//M矩阵元素对x的导数
mxArray *MDXDY=NULL;//申明M阵对x的导数阵
mlfAssign(&MDXDY,mlfDoubleMatrix(3,3,NULL,NULL));
//以下循环计算M阵对xy的偏导数,参照公式推导(4.12)式
for(int j=0;j<NP;j++)//总体坐标循环
{
setkernal(gx,gy,x[j],y[j]);
setkernalofdx(gx,gy,x[j],y[j]);
setkernalofdy(gx,gy,x[j],y[j]);
setkernalofdxdy(gx,gy,x[j],y[j]);
Melementofdxdy[0]=kerofdxdy;
Melementofdxdy[1]=kerofdy+(gx-x[j])*kerofdxdy;
Melementofdxdy[2]=kerofdx+(gy-y[j])*kerofdxdy;
Melementofdxdy[3]=kerofdy+(gx-x[j])*kerofdxdy;
Melementofdxdy[4]=2*(gx-x[j])*kerofdy+(gx-x[j])*(gx-x[j])*kerofdxdy;
Melementofdxdy[5]=ker+(gy-y[j])*kerofdy+(gx-x[j])*kerofdx+(gx-x[j])*(gy-y[j])*kerofdxdy;
Melementofdxdy[6]=kerofdx+(gy-y[j])*kerofdxdy;
Melementofdxdy[7]=ker+(gy-y[j])*kerofdy+(gx-x[j])*kerofdx+(gx-x[j])*(gy-y[j])*kerofdxdy;
Melementofdxdy[8]=2*(gy-y[j])*kerofdx+(gy-y[j])*(gy-y[j])*kerofdxdy;
mxArray *LXY=NULL;
mlfAssign(&LXY,mlfDoubleMatrix(3,3,Melementofdxdy,NULL));
mlfAssign(&MDXDY,mlfPlus(MDXDY,LXY));
mxDestroyArray(LXY);
}
return MDXDY;
}
mxArray * CRKPMDlg::setB(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
mxArray *M=NULL;//申明M矩阵
mlfAssign(&M,setM(gx,gy));
mxArray *MDX=NULL;//申明M阵对x的导数阵
mlfAssign(&MDX,setMDX(gx,gy));
mxArray *MD2X=NULL;//申明M阵对x的2阶导数阵
mlfAssign(&MD2X,setMD2X(gx,gy));
mxArray *MDY=NULL;//申明M阵对y的导数阵
mlfAssign(&MDY,setMDY(gx,gy));
mxArray *MD2Y=NULL;//申明M阵对y的2阶导数阵
mlfAssign(&MD2Y,setMD2Y(gx,gy));
mxArray *MDXDY=NULL;//申明M阵对xy的偏导数阵
mlfAssign(&MDXDY,setMDXDY(gx,gy));
//计算M阵的逆阵及其逆导数阵
mxArray *MI=NULL;//M阵的逆阵
mlfAssign(&MI,mlfInv(M));//给M阵逆阵赋值
mxArray *MIDX=NULL;//M阵的逆阵对x的导数阵,参照公式推导(4.16)式
mlfAssign(&MIDX,mlfMtimes(mlfMtimes(mlfScalar(-1),MI),mlfMtimes(MDX,MI)));
mxArray *MID2X=NULL;//M阵的逆阵对x的2阶导数阵,参照公式推导(4.16)式
mxArray *LID2X1=NULL,*LID2X2=NULL;//申明两个临时阵,以便计算MID2X
mlfAssign(&LID2X1,mlfMtimes(mlfScalar(2),mlfMtimes(MI,mlfMtimes(mlfMtimes(MDX,MI),mlfMtimes(MDX,MI)))));
mlfAssign(&LID2X2,mlfMtimes(MI,mlfMtimes(MD2X,MI)));
mlfAssign(&MID2X,mlfMinus(LID2X1,LID2X2));
mxDestroyArray(LID2X1);
mxDestroyArray(LID2X2);
mxArray *MIDY=NULL;//M阵的逆阵对y的导数阵
mlfAssign(&MIDY,mlfMtimes(mlfMtimes(mlfScalar(-1),MI),mlfMtimes(MDY,MI)));
mxArray *MID2Y=NULL;//M阵的逆阵对y的2阶导数阵
mxArray *LID2Y1=NULL,*LID2Y2=NULL;//申明两个临时阵,以便计算MID2Y
mlfAssign(&LID2Y1,mlfMtimes(mlfScalar(2),mlfMtimes(MI,mlfMtimes(mlfMtimes(MDY,MI),mlfMtimes(MDY,MI)))));
mlfAssign(&LID2Y2,mlfMtimes(MI,mlfMtimes(MD2Y,MI)));
mlfAssign(&MID2Y,mlfMinus(LID2Y1,LID2Y2));
mxDestroyArray(LID2Y1);
mxDestroyArray(LID2Y2);
mxArray *MIDXDY=NULL;//M阵的逆阵对xy的2阶导数阵
mxArray *LIDXDY1=NULL,*LIDXDY2=NULL;//申明两个临时阵,以便计算MIDXDY
mlfAssign(&LIDXDY1,mlfMtimes(mlfScalar(2),mlfMtimes(MI,mlfMtimes(MDX,mlfMtimes(MI,mlfMtimes(MDY,MI))))));
mlfAssign(&LIDXDY2,mlfMtimes(MI,mlfMtimes(MDXDY,MI)));
mlfAssign(&MIDXDY,mlfMinus(LIDXDY1,LIDXDY2));
mxDestroyArray(LIDXDY1);
mxDestroyArray(LIDXDY2);
//计算几何阵,参照公式推导(4.13),(4.14),(4.15)式
mxArray *B=NULL;//几何阵
mlfAssign(&B,mlfDoubleMatrix(3,NP,NULL,NULL));
double h0[3]={1,0,0};//H[0]
double hdx[3]={0,1,0};//H(X-XI)对x的导数
double hdy[3]={0,0,1};//H(X-XI)对y的导数
for(int k=0;k<NP;k++)//给几何阵赋值
{
setkernal(gx,gy,x[k],y[k]);
setkernalofdx(gx,gy,x[k],y[k]);
setkernalofd2x(gx,gy,x[k],y[k]);
setkernalofdy(gx,gy,x[k],y[k]);
setkernalofd2y(gx,gy,x[k],y[k]);
setkernalofdxdy(gx,gy,x[k],y[k]);
setbasicfun(gx,gy,x[k],y[k]);
mxArray *H0=NULL;
mxArray *HX=NULL;//申明基函数阵
mxArray *HXDX=NULL;//申明基函数阵对x的导数阵
mxArray *HXDY=NULL;//申明基函数阵对y的导数阵
mlfAssign(&H0,mlfDoubleMatrix(1,3,h0,NULL));
mlfAssign(&HX,mlfDoubleMatrix(3,1,H,NULL));
mlfAssign(&HXDX,mlfDoubleMatrix(3,1,hdx,NULL));
mlfAssign(&HXDY,mlfDoubleMatrix(3,1,hdy,NULL));
//计算B11
mxArray *LB11=NULL,*LB12=NULL,*LB13=NULL,*LB14=NULL,*LB15=NULL;//申明5个临时阵,以便计算B11
mlfAssign(&LB11,mlfMtimes(H0,mlfMtimes(MID2X,mlfMtimes(HX,mlfScalar(ker)))));
mlfAssign(&LB12,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MIDX,mlfMtimes(HX,mlfScalar(kerofdx))))));
mlfAssign(&LB13,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MIDX,mlfMtimes(HXDX,mlfScalar(ker))))));
mlfAssign(&LB14,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HXDX,mlfScalar(kerofdx))))));
mlfAssign(&LB15,mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HX,mlfScalar(kerofd2x)))));
mlfIndexAssign(&B,"(?,?)",mlfScalar(1),mlfScalar(k+1),mlfMtimes(mlfScalar(-1),mlfPlus(LB11,mlfPlus(LB12,mlfPlus(LB13,mlfPlus(LB14,LB15))))));
mxDestroyArray(LB11);
mxDestroyArray(LB12);
mxDestroyArray(LB13);
mxDestroyArray(LB14);
mxDestroyArray(LB15);
//计算B21
mxArray *LB21=NULL,*LB22=NULL,*LB23=NULL,*LB24=NULL,*LB25=NULL;//申明5个临时阵,以便计算B21
mlfAssign(&LB21,mlfMtimes(H0,mlfMtimes(MID2Y,mlfMtimes(HX,mlfScalar(ker)))));
mlfAssign(&LB22,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MIDY,mlfMtimes(HX,mlfScalar(kerofdy))))));
mlfAssign(&LB23,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MIDY,mlfMtimes(HXDY,mlfScalar(ker))))));
mlfAssign(&LB24,mlfMtimes(mlfScalar(2),mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HXDY,mlfScalar(kerofdy))))));
mlfAssign(&LB25,mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HX,mlfScalar(kerofd2y)))));
mlfIndexAssign(&B,"(?,?)",mlfScalar(2),mlfScalar(k+1),mlfMtimes(mlfScalar(-1),mlfPlus(LB21,mlfPlus(LB22,mlfPlus(LB23,mlfPlus(LB24,LB25))))));
mxDestroyArray(LB21);
mxDestroyArray(LB22);
mxDestroyArray(LB23);
mxDestroyArray(LB24);
mxDestroyArray(LB25);
//计算B31
mxArray *LB31=NULL,*LB32=NULL,*LB33=NULL,*LB34=NULL,*LB35=NULL,
*LB36=NULL,*LB37=NULL,*LB38=NULL;//申明8个临时阵,以便计算B31
mlfAssign(&LB31,mlfMtimes(H0,mlfMtimes(MIDXDY,mlfMtimes(HX,mlfScalar(ker)))));
mlfAssign(&LB32,mlfMtimes(H0,mlfMtimes(MIDX,mlfMtimes(HXDY,mlfScalar(ker)))));
mlfAssign(&LB33,mlfMtimes(H0,mlfMtimes(MIDX,mlfMtimes(HX,mlfScalar(kerofdy)))));
mlfAssign(&LB34,mlfMtimes(H0,mlfMtimes(MIDY,mlfMtimes(HXDX,mlfScalar(ker)))));
mlfAssign(&LB35,mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HXDX,mlfScalar(kerofdy)))));
mlfAssign(&LB36,mlfMtimes(H0,mlfMtimes(MIDY,mlfMtimes(HX,mlfScalar(kerofdx)))));
mlfAssign(&LB37,mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HXDY,mlfScalar(kerofdx)))));
mlfAssign(&LB38,mlfMtimes(H0,mlfMtimes(MI,mlfMtimes(HX,mlfScalar(kerofdxdy)))));
mlfIndexAssign(&B,"(?,?)",mlfScalar(3),mlfScalar(k+1),mlfMtimes(mlfScalar(-2),mlfPlus(LB31,mlfPlus(LB32,
mlfPlus(LB33,mlfPlus(LB34,mlfPlus(LB35,mlfPlus(LB36,mlfPlus(LB37,LB38)))))))));
mxDestroyArray(LB31);
mxDestroyArray(LB32);
mxDestroyArray(LB33);
mxDestroyArray(LB34);
mxDestroyArray(LB35);
mxDestroyArray(LB36);
mxDestroyArray(LB37);
mxDestroyArray(LB38);
mxDestroyArray(H0);
mxDestroyArray(HX);
mxDestroyArray(HXDX);
mxDestroyArray(HXDY);
}
return B;
}
mxArray * CRKPMDlg::setN(double gx,double gy)
{
UpdateData();
mlfEnterNewContext(0,0);//自动内存管理
//计算型函数矩阵N,参照公式推导(1.17),(5.11)式
mxArray *N=NULL;//声明型函数矩阵N
mlfAssign(&N,mlfDoubleMatrix(1,NP,NULL,NULL));
double h0[3]={1,0,0};//H[0]
mxArray *MI=NULL;//M阵的逆阵
mlfAssign(&MI,mlfInv(setM(gx,gy)));//给M阵逆阵赋值
for(int n=0;n<NP;n++)
{
setkernal(gx,gy,x[n],y[n]);
setbasicfun(gx,gy,x[n],y[n]);
mxArray *H0=NULL;
mxArray *HX=NULL;//申明基函数阵
mlfAssign(&H0,mlfDoubleMatrix(1,3,h0,NULL));
mlfAssign(&HX,mlfDoubleMatrix(3,1,H,NULL));
mlfIndexAssign(&N,"(?,?)",mlfScalar(1),mlfScalar(n+1),mlfMtimes(H0,
mlfMtimes(MI,mlfMtimes(HX,mlfScalar(ker)))));
mxDestroyArray(H0);
mxDestroyArray(HX);
}
return N;
}
void CRKPMDlg::setBgaosicordinate()
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -