大家好，最近的课程作业要求用有限元方法解微分方程。在制作矩形网格划分程序的过程中我看见有一本书上几个例题里面都有“网格平面旋转”的字样，不是很了解其中的作用。就是标红的部分。

这里面θ=0(我可以理解成矩形网格并没有在XY坐标系上旋转吗，即长平行于X，高平行于Y)

JXY2是节点的X,Y坐标值

R=sqrt((JXY2(i,1)+1)^2+JXY2(i,2)^2);R看起来似乎就是节点到(-1,0)的距离喽？

theta1=atan(JXY2(i,2)/(JXY2(i,1)+1));Θ1节点与(-1，0)连线与X轴的夹角

JXY2(i,1)=R*cos(theta/180*pi+theta1);那么最后新生成的坐标值的意义是什么呢？

Snipaste_2018-06-05_18-04-48.png (168.28 KB, 下载次数: 0)

2018-6-5 18:03 上传

本人实在有点菜，请大家帮帮忙哈~     ：)

clc

clear

%%%输入网格信息

H=0.06;%高

L=0.08;%长

Nx=4;%X向单元数

Ny=4;%Y向单元数

theta=0;%网格平面旋转角度

E=Nx*Ny;%总共节点数

N_sq=(2*Nx+1)*(2*Ny+1);%总二次节点数

N_line=(Nx+1)*(Ny+1);%总线性节点数

dx=L/Nx/2;%X向单元长

dy=H/Nx/2;%Y向单元高

%%%输入网格信息

%%%生成网格节点%%%

mesh_net=zeros(2*Nx+1,2*Ny+1);%初始化网格节点号

k=0;l=0;

for i=1:2:2*Ny+1

l=l+1;

for j=1:2:2*Nx+1

k=k+1;

mesh_net(i,j)=k;

if j+1<2*Nx+1

mesh_net(i,j+1)=k+(Nx+1)*(Ny+1)+(2*Nx)*(l-1);

end

end

end

l=0;

for i=2:2:2*Ny+1

l=l+1;

k=(Nx+1)*(Ny+1)+(l-1)*(Nx+2*Nx+1)+Nx;

for j=1:2*Nx+1

k=k+1;

mesh_net(i,j)=k;

end

end

%%%生成网格节点%%%

%%%生成网格坐标%%%

for i=1:2*Ny+1

for j=1:2*Nx+1

JXY2(mesh_net(i,j),1)=dx*(j-1);

JXY2(mesh_net(i,j),2)=dy*(i-1);

end

end

%%%生成网格坐标%%%

%%%网格平面旋转%%%

for i=1:length(JXY2(:,1))

R=sqrt((JXY2(i,1)+1)^2+JXY2(i,2)^2);

theta1=atan(JXY2(i,2)/(JXY2(i,1)+1));         %%%%%%%%%%%不懂%%%%%%%%%%%

JXY2(i,1)=R*cos(theta/180*pi+theta1);         %%%%投影之后X坐标加了1%%%%

JXY2(i,2)=R*sin(theta/180*pi+theta1);

end

%%%网格平面旋转%%%

%%%各单元二次节点序号%%%

k=0;

for i=1:Ny

for j=1:Nx

k=k+1;

JM2(k,:)=[mesh_net(2*i-1,2*j-1),mesh_net(2*i-1,2*j), ...

mesh_net(2*i-1,2*j+1),mesh_net(2*i,2*j-1),mesh_net(2*i,2*j), ...

mesh_net(2*i,2*j+1),mesh_net(2*i+1,2*j-1),mesh_net(2*i+1,2*j), ...

mesh_net(2*i+1,2*j+1)];

end

end

%%%各二次单元节点序号%%%

%%%各线性单元节点序号%%%

JM1=[JM2(:,1),JM2(:,3),JM2(:,9),JM2(:,7)];

JXY1=JXY2([1:N_line],:);

%%%各线性单元节点序号%%%

%%%边界点%%%

BPd=mesh_net(1,:);

BPu=mesh_net(2*Ny+1,:);

BPl=mesh_net(:,1);

BPl=mesh_net(:,2*Nx+1);

%%%边界点%%%

%%%连续方程BE值%%%

thetaxd=pi/2-theta/180*pi;%下边界外法线方向与x轴夹角

thetayd=pi-theta/180*pi;%下边界外法线方向与y轴夹角

thetaxu=pi/2+theta/180*pi;

thetayu=-theta/180*pi;        %%%%%%%%%%%这里这么算theta的正负%%%%%%%%%%%

thetaxr=theta/180*pi;        %%%%%%%%%%%这里这么算theta的正负%%%%%%%%%%%

thetayr=pi-thetaxr;

thetaxl=(180+theta)/180*pi;

thetayl=pi/2+theta/180*pi;

mesh_netd1=ones(Nx,1)*cos(thetaxd);%与下边界x方向余弦向量

mesh_netd2=ones(Nx,1)*cos(thetayd);%与下边界y方向余弦向量

mesh_netu1=ones(Nx,1)*cos(thetaxu);

mesh_netu2=ones(Nx,1)*cos(thetayu);

mesh_netl1=ones(Ny,1)*cos(thetaxl);

mesh_netl2=ones(Ny,1)*cos(thetayl);

mesh_netr1=ones(Ny,1)*cos(thetaxr);

mesh_netr2=ones(Ny,1)*cos(thetayr);

BEd=[[1:Nx]',ones(size([1:Nx]')),mesh_netd1,mesh_netd2];%[单元号,单元位置(1下,2右,3上,4左),与x方向余弦向量,与y方向余弦向量]

BEr=[[Nx:Nx:E]',2*ones(size([1:Ny]')),mesh_netr1,mesh_netr2];

BEu=[[Nx*(Ny-1)+1:Nx*Ny]',3*ones(size([1:Nx]')),mesh_netu1,mesh_netu2];

BEl=[[1:Nx:(Ny-1)*Nx+1]',4*ones(size([1:Ny]')),mesh_netl1,mesh_netl2];

%%%连续方程BE值%%%

%%%调用网格绘制程序