http://www-math.mit.edu/~persson/mesh/（已删除）

Distmesh 是一个matlab语言写的网格划分软件。
源文件可以从上面的网址获取。
这里按行讲解各个算例。
p01_demo:
概算例是一个单位圆（半径为1）的网格划分，划分后的网格为：
http://s1/bmiddle/618af195g84d42495ed20&690
以下逐行讲解该算例：
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<h<1
% 默认参数值为： 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;
效果如下：
http://s5/bmiddle/618af195g84d481e2bf64&690

http://s12/bmiddle/618af195g84d4864117eb&690
第二步：
根据fd的函数定义，移除边界外的点。
 p = p( feval_r( fd, p, varargin{:} ) <= geps, : );
varagin为fd，fh的附加参数，这里为空。
geps = 0.001 * h0;
也就是保留了到边界的距离以外0.001 * h0以内的点。
http://s5/bmiddle/618af195g84d4b07c3654&690  
根据网格密度函数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 );
这个时候产生的网格如下：
http://s7/bmiddle/618af195g84d4e7bcfb76&690

第三步：迭代
  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.
%
% 生成网格的边的集合，也就是相邻点之间连接的线段  
   bars = [ t(:,[1,2]); t(:,[1,3]); t(:,[2,3]) ];
      bars = unique ( sort ( bars, 2 ), 'rows' );

  end


 %
%  6. Move mesh points based on bar lengths L and forces F
%
%  Make a list of the bar vectors and lengths.
%  Set L0 to the desired lengths, F to the scalar bar forces,
%  and FVEC to the x, y components of the bar forces.
%
%  At the fixed positions, reset the force to 0.
%
    barvec = p(bars(:,1),:) - p(bars(:,2),:);   % 生成bar的矢量
    L = sqrt ( sum ( barvec.^2, 2 ) );  %计算bar的长度

  %根据每个bar的中点坐标，计算需要的三角形边的边长（这个在fh函数里控制）
    hbars = feval_r( fh, (p(bars(:,1),:)+p(bars(:,2),:))/2, varargin{:} );

% 计算 需要的bar的长度，已经乘上了两个scale参数 Fscale, sqrt ( sum(L.^2) / sum(hbars.^2) );
% 具体可参考他们的paper
    L0 = hbars * Fscale * sqrt ( sum(L.^2) / sum(hbars.^2) );

% 计算每个bar上力
    F = max ( L0 - L, 0 );

�r上力的分量，x,y方向                              
    Fvec = F ./ L * [1,1] .* barvec; 

% 计算Ftot， 每个节点上力的残量                  
    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;

% 根据每个节点上的受力，移动该点                        
    p = p + deltat * Ftot; 

%  7. Bring outside points back to the boundary
%
%  Use the numerical gradient of FD to project points back to the boundary.
%

  
    d = feval_r( fd, p, varargin{:} ); %计算点到边界距离
    ix = d > 0;
   
    %计算移动梯度，相对边界
    dgradx = ( feval_r(fd,[p(ix,1)+deps,p(ix,2)],varargin{:}) - d(ix) ) / deps;
    dgrady = ( feval_r(fd,[p(ix,1),p(ix,2)+deps],varargin{:}) - d(ix) ) / deps;

   %将这些移动到边界外的投射回边界上
    p(ix,:) = p(ix,:) - [ d(ix) .* dgradx, d(ix) .* dgrady ];
%
%  I needed the following modification to force the fixed points to stay.
%  Otherwise, they can drift outside the region and be lost.
%  JVB, 21 August 2008.
%
    p(1:nfix,1:2) = pfix;

    N = size ( p, 1 );
%
%  8. Termination criterion: All interior nodes move less than dptol (scaled)
%
    if ( max ( sqrt ( sum ( deltat * Ftot ( d < -geps,:).^2, 2 ) ) / h0 ) < dptol )
      break;
    end

  end

转自：http://blog.sina.com.cn/s/blog_6163bdeb0102dvay.html

源代码的地址已经改变，如下http://people.sc.fsu.edu/~jburkardt/m_src/distmesh/distmesh.html。程序给出了20个划分的例子，文件名为p1_demo.m~p20_demo.m，直接运行程序可以看到划分的动画效果，每个例子基本都是设置一些参数，然后调用distmesh_2d函数进行网格划分，最后得到划分的节点p和各三角形t，最后将划分的三角形绘制出来。划分结果如下

http://s10/middle/6163bdeb4afb1930a4009&amp;690

                  p1_demo                                    p14_demo

http://s4/middle/6163bdeb4afb1931f1743&amp;690

                p5_demo                                      p18_demo

进一步加注释的程序（需要学习的地方用颜色标记）：

function [ p, t ] = distmesh_2d ( fd, fh, h0, box, iteration_max, pfix, varargin )
%% DISTMESH_2D is a 2D mesh generator using distance functions.
%  Example:
%    Uniform Mesh on Unit Circle:
%      fd = inline('sqrt(sum(p.^2,2))-1','p');
%      [p,t] = distmesh_2d(fd,@huniform,0.2,[-1,-1;1,1],100,[]);
%    Rectangle with circular hole, refined at circle boundary:
%      fd = inline('ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5))','p');
%      fh = inline('min(4*sqrt(sum(p.^2,2))-1,2)','p');
%      [p,t] = distmesh_2d(fd,fh,0.05,[-1,-1;1,1],500,[-1,-1;-1,1;1,-1;1,1]);
%  Parameters:
%    Input, function FD, signed distance function d(x,y).
%    fd:d=fd(p)，p=[x y],fd为给定任一点到边界的有符号距离函数，负号表示在区域内，正号为在区域外
%    Input, function FH, scaled edge length function h(x,y).
%    fh:就是网格大小的函数
%    Input, real H0, the initial edge length.
%    h0:也就是h, 网格的初始大小
%    Input, real BOX(2,2), the bounding box [xmin,ymin; xmax,ymax].
%    box:最大外围矩形范围
%    Input, integer ITERATION_MAX, the maximum number of iterations.
%    The iteration might terminate sooner than this limit, if the program decides
%    that the mesh has converged.
%    iteration_max:允许的最大迭代次数
%    Input, real PFIX(NFIX,2), the fixed node positions.
%    pfix：网格中需要固定的点坐标，也就是一定需要出现在网格中的点
%    Input, VARARGIN, aditional parameters passed to FD and FH.%
%    Output, real P(N,2), the node positions.
%    p:网格点的x,y坐标
%    Output, integer T(NT,3), the triangle indices.
%    t:输出网格任意一个三角形的三个顶点
%  Local parameters:
%    Local, real GEPS, a tolerance for determining whether a point is "almost" inside
%    the region.  Setting GEPS = 0 makes this an exact test.  The program currently
%    sets it to 0.001 * H0, that is, a very small multiple of the desired side length
%    of a triangle.  GEPS is also used to determine whether a triangle falls inside
%    or outside the region.  In this case, the test is a little tighter.  The centroid
%    PMID is required to satisfy FD ( PMID ) <= -GEPS.
%    局部变量geps：容忍度，一个点是否属于区域，也在判断三角形是否属于区域内时使用

dptol = 0.001;              % 收敛精度
ttol = 0.1;                 % 三角形划分精度（百分比）
Fscale = 1.2;               % 放大比例
deltat = 0.2;               % 相当于柔度
geps = 0.001 * h0;          % 容忍度
deps = sqrt ( eps ) * h0;   % 微小变化dx
iteration = 0;              % 迭代次数
triangulation_count = 0;    % 三角形划分次数

%  1. Create the initial point distribution by generating a rectangular mesh
%  in the bounding box.
% 根据初始网格的大小h0，先把能涵盖欲划分区域的最大矩形划分为结构网格。
[ x, y ] = meshgrid ( box(1,1) : h0           : box(2,1), ...
    box(1,2) : h0*sqrt(3)/2 : box(2,2) );
%  Shift the even rows of the mesh to create a "perfect" mesh of equilateral triangles.
%  Store the X and Y coordinates together as our first estimate of "P", the mesh points
%  we want.
% 然后把偶数行的点整体向右平移半格，得到正三角形划分
x(2:2:end,:) = x(2:2:end,:) + h0 / 2;
p = [ x(:), y(:) ];
%  2. Remove mesh points that are outside the region,
%  then satisfy the density constraint.
%  Keep only points inside (or almost inside) the region, that is, FD(P) < GEPS.
% 根据fd的函数定义，移除边界外的点
p = p( feval_r( fd, p, varargin{:} ) <= geps, : ); % 1
%  Set R0, the relative probability to keep a point, based on the mesh density function.
r0 = 1 ./ feval_r( fh, p, varargin{:} ).^2;
%  Apply the rejection method to thin out points according to the density.
% 根据网格密度函数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' ); % 2

N = size ( p, 1 );

%  If ITERATION_MAX is 0, we're almost done.
%  For just this case, do the triangulation, then exit.
%  Setting ITERATION_MAX to 0 means that we can see the initial mesh
%  before any of the improvements have been made.
% 如果最大迭代次数为负，则直接结束
if ( iteration_max <= 0 )
    t = delaunayn ( p ); % 3
    triangulation_count = triangulation_count + 1;
    return
end
pold = inf; % 第一次迭代前设置旧的点的坐标为无穷
while ( iteration < iteration_max )
    iteration = iteration + 1;
    if ( mod ( iteration, 10 ) == 0 )
        fprintf ( 1, '  %d iterations,', iteration );
        fprintf ( 1, '  %d triangulations.\n', triangulation_count );
    end
   
    %  3. Retriangulation by the Delaunay algorithm.
    %  Was there large enough movement to retriangulate?
    %  If so, save the current positions, get the list of
    %  Delaunay triangles, compute the centroids, and keep
    %  the interior triangles (whose centroids are within the region).
    %
    % 先判断上次移动后的点和旧的点之间的移动距离，如果小于某个阀值，停止迭代
    if ( ttol < max ( sqrt ( sum ( ( p - pold ).^2, 2 ) ) / h0 ) )
        pold = p;               % 如果还可以移动，保存当前的节点
        t = delaunayn ( p );    % 利用delauny算法，生成三角形网格
        triangulation_count = triangulation_count + 1;          % 划分次数加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.
        % 生成网格的边的集合，也就是相邻点之间连接的线段
        bars = [ t(:,[1,2]); t(:,[1,3]); t(:,[2,3]) ];
        bars = unique ( sort ( bars, 2 ), 'rows' );
       
        %  5. Graphical output of the current mesh
        trimesh ( t, p(:,1), p(:,2), zeros(N,1) )               % 绘制划分的三角形% 3
        view(2), axis equal, axis off, drawnow
    end
   
    %  6. Move mesh points based on bar lengths L and forces F
    %  Make a list of the bar vectors and lengths.
    %  Set L0 to the desired lengths, F to the scalar bar forces,
    %  and FVEC to the x, y components of the bar forces.
    %  At the fixed positions, reset the force to 0.
    barvec = p(bars(:,1),:) - p(bars(:,2),:);   % 生成bar的矢量
    L = sqrt ( sum ( barvec.^2, 2 ) );          % 计算bar的长度
    % 根据每个bar的中点坐标，计算需要的三角形边的边长（这个在fh函数里控制）
    hbars = feval_r( fh, (p(bars(:,1),:)+p(bars(:,2),:))/2, varargin{:} );
    % 计算需要的bar的长度，已经乘上了两个scale参数 Fscale, sqrt ( sum(L.^2) / sum(hbars.^2) );
    % 具体可参考他们的paper
    L0 = hbars * Fscale * sqrt ( sum(L.^2) / sum(hbars.^2) );
    % 计算每个bar上力
    F = max ( L0 - L, 0 );
    % bar上力的分量，x,y方向
    Fvec = F ./ L * [1,1] .* barvec;
    % 计算Ftot， 每个节点上力的残量
    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;
    % 根据每个节点上的受力，移动该点
    p = p + deltat * Ftot;
   
    %  7. Bring outside points back to the boundary
    %  Use the numerical gradient of FD to project points back to the boundary.
    d = feval_r( fd, p, varargin{:} ); % 计算点到边界距离
    ix = d > 0;
    % 计算移动梯度，相对边界
    dgradx = ( feval_r(fd,[p(ix,1)+deps,p(ix,2)],varargin{:}) - d(ix) ) / deps; % 4
    dgrady = ( feval_r(fd,[p(ix,1),p(ix,2)+deps],varargin{:}) - d(ix) ) / deps;
    % 将这些移动到边界外的投射回边界上
    p(ix,:) = p(ix,:) - [ d(ix) .* dgradx, d(ix) .* dgrady ];
    %  I needed the following modification to force the fixed points to stay.
    %  Otherwise, they can drift outside the region and be lost.
    % 修正，以免一些点移到区域外而丢失
    p(1:nfix,1:2) = pfix;
    N = size ( p, 1 );
   
    %  8. Termination criterion: All interior nodes move less than dptol (scaled)
    if ( max ( sqrt ( sum ( deltat * Ftot ( d < -geps,:).^2, 2 ) ) / h0 ) < dptol )
        break;
    end
end
end

值得学习的地方：

1.feval_r(fd, p, varargin{:})

这相当于是回调函数的用法，将某些完成特定功能的函数作为输入参数传递进来，需要实现某些功能时则调用此函数，这样当想改变特定功能时，直接改变传进来的函数句柄就可以了。而且完成特定功能的函数的额外参数可以由varargin传入。

2.unique ( p, 'rows' );

unique函数属于时间序列分析中的功能，最近借了一本书正好看到，时间序列有一些列处理方法函数，可以很方便的处理作为时间序列的向量。

3.delaunayn 、trimesh

这是涉及三角形划分的matlab功能，感觉挺使用的。一方面，有限元方法需要网格划分，有这些函数，如虎添翼，另外在matlab实现三维的类似OpenGL的显示、操作等，使用patch最方便了，而patch最方便的使用方法是传入点和面，而如何构造点和面呢，这里给出了很好的答案。

4.deps = sqrt ( eps ) * h0;   % 微小变化dx

dgradx = ( feval_r(fd,[p(ix,1)+deps,p(ix,2)],varargin{:}) - d(ix) ) / deps;

这个相当于是求函数微分、梯度，本来matlab函数传入的量是离散的，感觉顶多求差分，突然发现这种想法太简单了，函数的定义域是连续的，连续的就可以求微分，这里提供了一种很好方法。

转自：http://blog.sina.com.cn/s/blog_6163bdeb0102dvbw.html