function [ p, t ] = distmesh_2d ( fd, fh, h0, box, iteration_max, pfix, varargin )
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; % 三角形划分次数
[ x, y ] = meshgrid ( box(1,1) : h0 : box(2,1), ...
box(1,2) : h0*sqrt(3)/2 : box(2,2) );
x(2:2:end,:) = x(2:2:end,:) + h0 / 2;
p = [ x(:), y(:) ];
r0 = 1 ./ feval_r( fh, p, varargin{:} ).^2;
p = [ pfix; p(rand(size(p,1),1) < r0 ./ max ( r0 ),: ) ];
[ nfix, dummy ] = size ( pfix );
p = unique ( p, 'rows' );
N = size ( p, 1 );
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
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, : );
bars = [ t(:,[1,2]); t(:,[1,3]); t(:,[2,3]) ];
bars = unique ( sort ( bars, 2 ), 'rows' );
trimesh ( t, p(:,1), p(:,2), zeros(N,1) ) % 绘制划分的三角形% 3
view(2), axis equal, axis off, drawnow
end
barvec = p(bars(:,1),:) - p(bars(:,2),:); % 生成bar的矢量
L = sqrt ( sum ( barvec.^2, 2 ) ); % 计算bar的长度
hbars = feval_r( 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;
% 根据每个节点上的受力，移动该点
p = p + deltat * Ftot;
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
