之前轉載了一篇博客http://blog.sina.com.cn/s/blog_6163bdeb0102dvay.html，講Matlab網格划分程序Distmesh，看了看程序，感覺程序寫得有很多值得學的地方，所以又自己重新看了一看，加了一些注釋，最后再總結一下學到的東西吧。
源代碼的地址已經改變，如下http://people.sc.fsu.edu/~jburkardt/m_src/distmesh/distmesh.html。程序給出了20個划分的例子，文件名為p1_demo.m~p20_demo.m，直接運行程序可以看到划分的動畫效果，每個例子基本都是設置一些參數，然后調用distmesh_2d函數進行網格划分，最后得到划分的節點p和各三角形t，最后將划分的三角形繪制出來。划分結果如下
p1_demo p14_demo
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函數傳入的量是離散的，感覺頂多求差分，突然發現這種想法太簡單了，函數的定義域是連續的，連續的就可以求微分，這里提供了一種很好方法。
這也提供了一種很好的提示，以前一直拿c、c++的編程思想用matlab，現在發現使用matlab的更有效率的方法是使用數學函數的思想，matlab的函數不是c、c++中響應函數的那個函數，而是定義在實數域上的數學函數，用來處理各種數學問題。也可能是之前理解c、c++中的函數就應該用這種思想吧！
此程序中就實現了各種數學意義上的計算，如求點到線的距離、求梯度、算分量等，看了這個程序，才感覺matlab確實是數學的語言工具。