DistMesh代码划分网格

对求解域进行网格划分是有限元求解的关键步骤之一，网格划分的好坏直接影响计算结果的精度。在使用Matlab进行有限元编程时，可以采用Per-Olof教授提供的DistMesh代码来生成三角形网格。源代码下载：DistMesh。

代码思路

节点生成

该代码首先通过区域范围bbox和间距h0生成均匀节点，并将偶数行节点进行偏移；之后再通过边界条件，将域外节点舍去；最后通过概率采样的方法来取舍节点实现节点的均匀和非均匀布置，其原理是距离加密区越近的节点被保留的概率越大，而远的节点大概率被舍弃。

[x,y]=meshgrid(bbox(1,1):h0:bbox(2,1),bbox(1,2):h0*sqrt(3)/2:bbox(2,2));     % 根据范围生成节点
x(2:2:end,:)=x(2:2:end,:)+h0/2;                      % 偶数节点偏移
p=[x(:),y(:)];                                     
p=p(feval(fd,p,varargin{:})<geps,:);                 % 舍弃域外节点

r0=1./feval(fh,p,varargin{:}).^2;                    % 节点距离加密位置距离平方的导数 
p=p(rand(size(p,1),1)<r0./max(r0),:);                % 根据节点出现概率与距离加密区距离呈反比筛选节点

网格成型

在上述节点生成之后，采用Delaunay三角剖分得到单元网格，将单元中心在求解域外的单元舍去，之后得到每条单元边的节点编号矩阵bars。

t=delaunayn(p);                                  % 三角形单元节点编号
pmid=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;    % 中心点
t=t(feval(fd,pmid,varargin{:})<-geps,:);         % 舍去中心点在外的单元

bars=[t(:,[1,2]);t(:,[1,3]);t(:,[2,3])];         
bars=unique(sort(bars,2),'rows');                % 单元边节点编号

在得到现有三角形单元之后，在单元期望边界长度的控制下，通过节点力的平衡进行节点位置更新迭代，该部分的详细内容可见Per-Olof教授的论文【1】。

barvec=p(bars(:,1),:)-p(bars(:,2),:);              % 单元边矢量
L=sqrt(sum(barvec.^2,2));                          % 单元边长度
hbars=feval(fh,(p(bars(:,1),:)+p(bars(:,2),:))/2,varargin{:});   % 单元中心点期望边界长度比例
L0=hbars*Fscale*sqrt(sum(L.^2)/sum(hbars.^2));     % 期望长度

F=max(L0-L,0);                                     % 边荷载
Fvec=F./L*[1,1].*barvec;                           % 边荷载分量
Ftot=full(sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2));  % 节点合力
Ftot(1:size(pfix,1),:)=0;                          % 固定点节点合理为0
p=p+deltat*Ftot;                                   % 更新节点位置
 
d=feval(fd,p,varargin{:}); ix=d>0;                 % 寻找外部节点
dgradx=(feval(fd,[p(ix,1)+deps,p(ix,2)],varargin{:})-d(ix))/deps; 
dgrady=(feval(fd,[p(ix,1),p(ix,2)+deps],varargin{:})-d(ix))/deps; 
dgrad2=dgradx.^2+dgrady.^2;
p(ix,:)=p(ix,:)-[d(ix).*dgradx./dgrad2,d(ix).*dgrady./dgrad2];    % 根据梯度移动回域内

代码使用

以上是代码的核心内容和大致思路，对于如何根据自己的需求使用该代码，主要需要理解域定义方式以及单元边长控制方式，这就涉及函数fd和fh如何定义。

函数fd定义为节点到边界的距离，域内结果为负值或者零，这样就可以很容易的判定节点是否为域内节点。工具包内除了提供圆形(dcircle)、矩形(drectangle)以及多边形(dpoly)域代码，还提供dunion 和ddiff 实现求解域范围相加和相减。基于此则可以定义所需求解域，以下给出四分之一圆环和带孔洞圆角矩形两个例子。

% 带孔洞圆角矩形
fd1 = @(p) drectangle(p,-2,2,-2,2);
fd2 = @(p) dunion(ddiff(fd1(p),drectangle(p,1.6,2,1.6,2)),dcircle(p,1.6,1.6,0.4));
fd = @(p) ddiff(fd2(p),dcircle(p,0,0,0.2));
fh = @(p) 0.05-0.4*fd(p)
[p,t]=distmesh2d(fd,fh,0.05,[-2,-2;2,2],[1.6,2;2,1.6;2,-2;2,-2;-2,-2]);

% 四分之一扇形
fd=@(p) max(ddiff(dcircle(p,0,0,1),dcircle(p,0,0,0.5)),-min(p(:,1),p(:,2)));
[p,t]=distmesh2d(fd,@huniform,0.05,[0,0;1,1],[0,0.5;0,1;1,0;0.5,0;0,0.5]);

函数fh定义节点单元边界长度，该函数应当满足加密区附近节点输出值小，远处输出值大。例如上述例子的四分之一圆环在边界处进行网格加密的代码如下：

fd=@(p) max(ddiff(dcircle(p,0,0,1),dcircle(p,0,0,0.5)),-min(p(:,1),p(:,2)));
fh = @(p) 0.05-0.25*fd(p)
[p,t]=distmesh2d(fd,fh,0.05,[0,0;1,1],[0,0.5;0,1;1,0;0.5,0;0,0.5]);

对比四分之一圆环两图可以发现，边缘处节点差不多，但是下面的图内部单元反而更大，这是因为节点初始间距h0是一样的，但是图二由于引入尺寸控制函数fh会删除一些节点。所以需要引入尺寸控制函数加密时，也因降低节点初始间距，增加初始节点。

fd=@(p) max(ddiff(dcircle(p,0,0,1),dcircle(p,0,0,0.5)),-min(p(:,1),p(:,2)));
fh = @(p) 0.05-0.25*fd(p)
[p,t]=distmesh2d(fd,fh,0.025,[0,0;1,1],[0,0.5;0,1;1,0;0.5,0;0,0.5]);

以下再给出一个矩形内部圆形部分区域加密的网格生成例子，源码请查看微信公众号：结构设计札记

通过测试发现要想获得比较理想的网格划分需要不断调试参数，特别是控制单元尺寸函数，一般可以表达成如下形式：
$$fh=@(p)a+b\ast abs(boundary)$$
boundary表示到加密边界函数，a调节边界附近尺寸控制参数，b则控制边界四周加密区的大小，b值大则加密区小。

参考文献

【1】A Simple Mesh Generator in MATLAB