基于Matlab的三角形网格自动生成小工具 | Distmesh基础

上几期推文已经给大家介绍了二维三角形、四边形单元网格自动生成算法，不过都是局限于矩形区域，那么遇到弧形、不规则几何区域该怎么构建网格呢？

本次分享的是基于Matlab的三角形网格自动生成小工具——Distmesh^[1]。该小工具可有效的在不规则区域内构建三角形网格，仅以三角形网格为例，其他单元类型的网格可参考别的工具（兴趣研究）。

Distmesh介绍

该程序^[2]是由美国加州大学伯克利分校数学系的Per-Olof Persson教授和麻省理工学院数学系的Gilbert教授合力编制，用于生成非结构化三角形和四面体网格。

采用符号距离函数(Signed Distance Functions) 构造几何区域，空间中任一点到域边界的最短距离。区域内符号为负，区域外为正，以此来构建各种复杂边界。我只学会了一些皮毛，抛砖引玉，将这个小工具介绍一下~

实际的网格生成过程中，DistMesh 使用 Delaunay MATLAB 中的三角剖分例程并尝试优化节点位置，通过基于力的平滑程序，拓扑定期更新。

边界点只允许切向移动，使用距离函数的投影边界。这个迭代程序通常会产生形状非常好的网格，该程序使用起来比较容易上手，可根据个人需要进行特殊修改。

Distmesh函数调用格式如下，先介绍一下形参的含义：

function [p,t]=distmesh2d(fd,fh,h0,bbox,pfix,varargin)

形参解读

1. fd 描述几何距离函数，例：fd = @(p) sqrt(sum(p.^2,2)) - 1描述单位圆；
2. fh描述单元尺寸h的函数；
3. h0控制单元尺寸，对于均匀网格，h0=常数，单元实际尺寸会比该值大一点；
4. bbox定义生成网格区域；
5. pfix指定在生成的网格中应始终存在（固定）的点数，用于几何区域角点（矩形角点）；
6. varargin指定fd和fh控制参数。

输出参数

1. p节点坐标矩阵；
2. t单元节点信息；

圆形区域网格生成

看到以上调用函数的介绍，相信大家，还是比较迷，因为我刚接触的时候也很迷，接下来拿几个简单的案例练练手，就会用了。

圆形区域描述：几何边界点到原点的距离：    ，均匀网格分布，初始单元尺寸    控制在0.2左右，网格生成区域定义在以    和    为顶点的矩形区域内，无固定点，所以pfix为[]。

围绕圆形区域，我将根据自己对相关文献的理解给出三种构建方法：

   

   

方式一

fd=@(p) sqrt(sum(p.^2,2))-1;
% fd=@(p) sqrt((p(:,1)-1).^2+(p(:,2)-1).^2)-1;
[p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);

代码解读

1. 采用匿名函数构建几何距离函数fd；
2. p(:,1)、p(:,2)可理解为坐标      和      ；
3. fh=@huniform说明是均匀网格，有关网格密度尺寸控制的讲解将会在下次推文中说明；
4.      ，关于这个h0，文献中是这样解释的：    网格刚开始形成时，是将区域划分为很多规则的矩形小单元，      就是这些小矩形的单元尺寸；
5. [-1,-1;1,1]生成的网格区域定义在该矩形区域内（两个顶点描述）；
6. pfix=[]说明没有要约束的结点。

方式二

[xx,yy]=meshgrid(-1.1:0.1:1.1,-1.1:0.1:1.1); 
dd=sqrt(xx.^2+yy.^2)-1;
[p,t]=distmesh2d(@dmatrix,@huniform,0.2,[-1,-1;1,1],[],xx,yy,dd);

代码解读

1. 采用离散化的方式建立几何区域，首先使用meshgrid建立矩阵点集；
2. 距离函数显式表达：      ，在形参中以@dmatrix表示离散的点集；
3. varargin为xx,yy,dd，该三个参数是@dmatrix的控制参数，通过线性插值的方法描述几何区域；
4. 该方法可以描述任意复杂区域，但是我还用的不到家，望感兴趣的伙伴可以进行更深的理解。

方式三

fd=@(p) dcircle(p,0,0,1);
[p,t]=distmesh2d(fd,@huniform,0.1,[-1,-1;1,1],[]);

代码解读

1. 采用作者编写的dcircle函数构建圆形区域，形参分别为p，圆心坐标(0,0)，半径 1。

以下三图分别将h0设置为0.05，0.1，0.2，0.3：

   

   

   

   

多边形区域网格生成

pv = [-0.4 -0.5;0.4 -0.2;0.4 -0.7;1.5 -0.4;0.9 0.1;1.6 0.8;...
0.5 0.5;0.2 1;0.1 0.4;-0.7 0.7;-0.4 -0.5];
fd = @(p) dpoly(p,pv);
[p,t] = distmesh2d(fd,@huniform,0.1,[-1,-1; 2,1],pv);

代码解读

1. pv数组定义多边形角点坐标；
2. 调用作者编写的dpoly函数构建多边形区域，形参为坐标点p，角点坐标数组pv；
3. 使用多边形时就需要定义固定点的位置了，因为在角点处符号距离函数描述的不准确，导致程序迭代次数增多，不稳定，一般取为多边形角点坐标，之前的pfix=[]，现在pfix=pv；

 

矩形区域

矩形区域应用的较为广泛，所以作者编写了相应的调用函数drectangle(p,x1,x2,y1,y2)，输入矩形对角坐标即可：

fd = @(p) drectangle(p,-1,1,-1,1);
pfix = [-1,-1;-1,1;1,-1;1,1];
[p,t]=distmesh2d(fd,@huniform,0.1,[-1,-1;1,1],pfix);

 

椭圆区域网格生成

椭圆方程：

其中：    .

fd=@(p) p(:,1).^2/2^2+p(:,2).^2/1^2-1;
[p,t]=distmesh2d(fd,@huniform,0.12,[-2,-1;2,1],[]);

 

也可以在标准椭圆的基础上进行坐标偏移：

fd=@(p) (p(:,1)-1).^2/2^2+(p(:,2)-1).^2/1^2-1;
[p,t]=distmesh2d(fd,@huniform,0.12,[-3,-3;3,3],[]);

 

参考资料