Matlab网格划分程序Distmesh讲解

Matlab网格划分程序Distmesh讲解（一）

/~persson/mesh/

Distmesh 是一个matlab语言写的网格划分软件。

源文件可以从上面的网址获取。

这里按行讲解各个算例。

p01_demo:

概算例是一个单位圆（半径为1）的网格划分，划分后的网格为：

以下逐行讲解该算例：

function p01_demo ( iteration_max, h )

% Parameters:

% Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH % should take. (The program might take fewer iterations if it detects convergence.)

% Input, real h, the mesh spacing parameter.

% 这里有两个输入参数，一个是ITERATION_MAX,迭代的最大次数。

% 另一个是h, 网格划分的大小。0

% 默认参数值为：ITERATION=200 h=0.1

[ p, t ]=distmesh_2d( fd, fh, h0, box, iteration_max, fixed );

函数需要至少六个参数。

d = fd ( p )，p=[x y]

fd 给定任一点到边界的距离函数，本例中定义为：d = sqrt(x^2+y^2)-1;

fh, scaled edge length function h(x,y). 也就是网格大小的函数。

h0 也就是h, 网格的大小

real BOX(2,2), the bounding box [xmin,ymin; xmax,ymax].

最大外围矩形范围。本例中为[0,0;1,1]

ITERATION_MAX, the maximum number of iterations.

real PFIX(NFIX,2), the fixed node positions. 网格中需要固定的点坐标，也就是一定需要出现在网格中的点。

输出参数：

real P(N,2), the node positions. 网格点的x,y坐标

integer T(NT,3), the triangle indices. 输出网格任一一个三角形的三个顶点。

第一步：

[ x, y ] = meshgrid ( box(1,1) : h0 : box(2,1), ...

box(1,2) : h0*sqrt(3)/2 : box(2,2) );

根据h0，网格的大小，先把能涵盖欲划分区域的最大矩形划分为结构网格。

然后把偶数行的点整体向右平移半格，

x(2:2:end,:) = x(2:2:end,:) + h0 / 2;

效果如下：

第二步：

根据fd的函数定义，移除边界外的点。

p = p( feval_r( fd, p, varargin{:} ) <= geps, : ); varagin为fd，fh的附加参数，这里为空。geps = 0.001 * h0;

也就是保留了到边界的距离以外0.001 * h0以内的点。

根据网格密度函数fh，每个点上产生一个0-1随机数，判断是否小于r0/max(r0) 大于的话，改点被删除。

p = [ pfix; p(rand(size(p,1),1) < r0 ./ max ( r0 ),: ) ];

[ nfix, dummy ] = size ( pfix );

当指定了某些点要保留的时候，把保留的点加入，删除重复的点。

% Especially when the user has included fixed points, we may have a few

% duplicates. Get rid of any you spot.

%

p = unique ( p, 'rows' );

N = size ( p, 1 );

这个时候产生的网格如下：

第三步：迭代

pold = inf; %第一次迭代前设置旧的点的坐标为无穷

while ( iteration < iteration_max )

iteration = iteration + 1;

%先判断上次移动后的点和旧的点之间的移动距离，如果小于某个阀值，停止迭代if ( ttol < max ( sqrt ( sum ( ( p - pold ).^2, 2 ) ) / h0 ) )

pold = p; %如果还可以移动，保存当前的节点

t = delaunayn ( p ); %利用delauny算法，生成三角形网格

triangulation_count = triangulation_count + 1;

pmid = ( p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:) ) / 3; %计算三角形的重心。

t = t( feval_r( fd, pmid, varargin{:} ) <= -geps, : ); % 移除重心在边界外部的三角形

% 4. Describe each bar by a unique pair of nodes.

%

% 生成网格的边的集合，也就是相邻点之间连接的线段