[ 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. %
% ⽣成⽹格的边的集合，也就是相邻点之间连接的线段 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
