之前转载了一篇博客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的矢量
加载中…