http://persson.berkeley.edu/distmesh/
Distmesh 是一个matlab语言写的网格划分软件。
源文件可以从上面的网址获取。
这里按行讲解各个算例。
p01_demo:
概算例是一个单位圆(半径为1)的网格划分,划分后的网格为:
以下逐行讲解该算例:
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;
效果如下:
第二步:
根据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
Distmesh 是一个matlab语言写的网格划分软件。
源文件可以从上面的网址获取。
这里按行讲解各个算例。
p01_demo:
概算例是一个单位圆(半径为1)的网格划分,划分后的网格为:
![Matlab网格划分程序Distmesh讲解(一) Matlab网格划分程序Distmesh讲解(一)](https://i-blog.csdnimg.cn/blog_migrate/1be01cdc8f7ee92d171176030215a8d7.jpeg)
以下逐行讲解该算例:
function p01_demo ( iteration_max, h )
%
%
%
%
% 这里有两个输入参数,一个是ITERATION_MAX,迭代的最大次数。
% 另一个是h, 网格划分的大小。 0<h<1
% 默认参数值为: ITERATION=200
[ p, t ]=distmesh_2d( fd, fh, h0, box, iteration_max, fixed );
函数需要至少六个参数。
d = fd ( p ),
fd
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.
integer T(NT,3), the triangle indices. 输出网格任一一个三角形的三个顶点。
第一步:
根据h0,网格的大小,先把能涵盖欲划分区域的最大矩形划分为结构网格。
然后把偶数行的点整体向右平移半格,
效果如下:
![Matlab网格划分程序Distmesh讲解(一) Matlab网格划分程序Distmesh讲解(一)](https://i-blog.csdnimg.cn/blog_migrate/44553c7b0be0e3cedea0d5ca50a20c7d.jpeg)
![Matlab网格划分程序Distmesh讲解(一) Matlab网格划分程序Distmesh讲解(一)](https://i-blog.csdnimg.cn/blog_migrate/63deeb35c4540f83b003170c24b775c9.jpeg)
第二步:
根据fd的函数定义,移除边界外的点。
varagin为fd,fh的附加参数,这里为空。
geps = 0.001 * h0;
也就是保留了到边界的距离以外0.001 * h0以内的点。
![Matlab网格划分程序Distmesh讲解(一) Matlab网格划分程序Distmesh讲解(一)](https://i-blog.csdnimg.cn/blog_migrate/91168f254953456045723cb73754e1fb.jpeg)
根据网格密度函数fh,每个点上产生一个0-1随机数,判断是否小于r0/max(r0)
大于的话,改点被删除。
当指定了某些点要保留的时候,把保留的点加入,删除重复的点。
%
%
%
这个时候产生的网格如下:
![Matlab网格划分程序Distmesh讲解(一) Matlab网格划分程序Distmesh讲解(一)](https://i-blog.csdnimg.cn/blog_migrate/2c8c4f532e3b6b87c65f551b20a4c5d7.jpeg)
第三步:迭代
while ( iteration < iteration_max )
%先判断上次移动后的点和旧的点之间的移动距离,如果小于某个阀值,停止迭代
%
%
% 生成网格的边的集合,也就是相邻点之间连接的线段
%
%
%
%
%
%
%
%
% 计算 需要的bar的长度,已经乘上了两个scale参数 Fscale, sqrt ( sum(L.^2) / sum(hbars.^2) );
% 具体可参考他们的paper
% 计算每个bar上力
�r上力的分量,x,y方向
% 计算Ftot, 每个节点上力的残量
%对于固定点,力的残量为零
% 根据每个节点上的受力,移动该点
%
%
%
%
%
%
%
%
%
%
%
%