有限元-区域划分

区域划分

这里采用例子$X=[0,8],Y=[0,8],h_x=2,h_y=2$，绘制如下网格划分图像，以这个图像为例，将说明根据小正方形编号确定小正方形的四个顶点$(left\_down,right\_down,right\_up,left\_up)$

在图中节点的序号为黑色数字，小正方形的序号为括号外红色数字，括号内的内容表示小正方形的行和列。绿色线条表示将每个小正方形划分为两个小三角形。

这里需要处理的是根据小正方形的序号找到小正方形的四个节点序号，记节点列数为$lenx$，节点行数为$leny$，显然，小长方形的列数为$lenx-1$，行数为$leny-1$。

如图所示，可知 $lenx=5,leny=5$.

因此可以根据小正方形的序号计算出小正方形所在的行和列：

记小正方形的序号为 $i$ 这里 $i=1,2,3,...24$.

• 行：$row=ceil(i/(lenx-1))$，这里$ceil$函数表示向上取整

• 列：$list=i-(row-1)\ast (lenx-1)$

• 例如：当$i=6$时，

  $row=ceil(6/(5-1))=ceil(1.5)=2，list=6-(2-1)\times (5-1)=2$，因此方块序号$6$对应的位置为$(2,2)$。

进一步，根据方格位置坐标可以确定方格的四个顶点位置，如下：

$left\_down=(row-1)\ast lenx+list$
$right\_down=row\ast lenx+list$
$right\_up=row\ast lenx+(list+1)$
$left\_up=(row-1)\ast lenx+(list+1)$

例如：当$i=6$时，$row=2,list=2$，从而有

$left\_down=(row-1)\ast lenx+list=(2-1)\ast 5+2=7$
$right\_down=row\ast lenx+list=2\ast 5+2=12$
$right\_up=row\ast lenx+(list+1)=2\ast 5+(2+1)=13$
$left\_up=(row-1)\ast lenx+(list+1)=(2-1)\ast 5+(2+1)=8$
这完全可以和我们的小正方形的四个顶点对应起来，因此，在这个小正方形里面划分的两个小三角形的三个顶点描述是容易的（采用逆时针的顶点描述），即

$[left\_down,right\_down,right\_up{]}^T$和$[right\_up,left\_up,left\_down{]}^T$

于是，按照上述节点序号顺序保存节点（共两行，第一行为节点的x坐标，第二行为节点的y坐标，列数对应节点的序号）从而得到节点矩阵$Pos$如下:（这里仅展示部分数据，全部数据运行查看）

节点序号（列）	1	2	3	4	5	6	7	8	9	10
$x$	0	0	0	0	0	2	2	2	2	2
$y$	0	2	4	6	8	0	2	4	6	8

• 存储矩阵Pos时，并没有保存节点序号，所以Pos矩阵大小为$2\times (lenx\times leny)$

以及连接矩阵$T$（表中的数值表示节点序号）如下：（这里仅展示部分数据，全部数据运行查看）

三角形顶点A	1	7	2	8	3	9	4	10
三角形顶点B	6	2	7	3	8	4	9	5
三角形顶点C	7	1	8	2	9	3	10	4

• 连接矩阵$T$保存每个小三角形的节点序号，共有$(lenx-1)\times (leny-1)$个小正方形，每个小正方形下面有两个小三角形，所以$T$的大小为$3\times [2\times (lenx-1)\times (leny-1)]$

最后附上这段实例的代码

clc,clear
nx = 4;
ny = 4;
X =[0,8];
Y =[0,8];
hx = (X(2)-X(1))/nx;                                       
hy = (Y(2)-Y(1))/ny;
x = X(1):hx:X(2);
y = Y(1):hy:Y(2);
lenx = length(x);
leny = length(y);

% 网格划分 pos对应节点， T的每一列对应三角形的三个顶点（逆时针）
pos = zeros(2,lenx*leny);
for i = 1:lenx
    pos(1,(1 + (i-1)*leny):i*(leny)) = x(i);
    pos(2,(1 + (i-1)*leny):i*(leny)) = y;
    plot(x(i)*ones(leny),y,'color',[0.30,0.75,0.93])
    hold on
    plot(x,y(i)*ones(lenx),'color',[0.30,0.75,0.93])
end
for i = 1:length(pos)
    text(pos(1,i),pos(2,i),num2str(i),'HorizontalAlignment','center')
end
T = zeros(3,2*(lenx-1)*(leny-1));
for i = 1:(lenx-1)*(leny-1)
    row = ceil(i/(lenx-1));
    list = i - (row - 1) * (lenx - 1);
    % 第 i 个正方形的四个顶点为 x_row x_(row+1) y_list y_(list+1)
    T(1,2*i-1) = (row - 1) * lenx + list;
    T(2,2*i-1) = row * lenx + list;
    T(3,2*i-1) = row * lenx + (list + 1);
    T(1,2*i) = row * lenx + (list + 1);
    T(2,2*i) = (row - 1) * lenx + (list + 1);
    T(3,2*i) = (row - 1) * lenx + list;
    x_pos = (pos(1,T(1,2*i-1))+pos(1,T(2,2*i-1)))/2;
    y_pos = (pos(2,T(2,2*i-1))+pos(2,T(3,2*i-1)))/2;
    text(x_pos,y_pos,[num2str(i),'(',num2str(row+1),',',num2str(list),')'],...
        'HorizontalAlignment','center',...
        'Color','red')
    hold on
    cen_Line=[pos(:,T(1,2*i-1)),pos(:,T(3,2*i-1))];
    plot(cen_Line(1,:),cen_Line(2,:),'color',[0.73,0.88,0.15])
end

set(gca,'XLim',[X(1)-1,X(2)+1],'YLim',[Y(1)-1,Y(2)+1])
axis equal

• 小结：这个代码中，矩形区域的大小由X，Y控制，你可以随意更改，区间划分采用划分，划分个数分别由$nx,ny$两个参数控制，这两个参数你也可以随意更改。

• ~~需要注意到的是，在产生网格划分示意图时，标记的字体大小是使用系统默认字体大小，所以当网格划分较密集时，所出现带标记的图像并不好看，你可以将标记文本的代码注释掉。~~如下的图像：