Matlab网格划分程序Distmesh讲解

1、Matlab网格划分程序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

2、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 yfd给定任一点到边界的距离函数，本例中定义为： d = sqrt(

3、x2+y2)-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,1ITERATION_MAX, the maximum number of iterations.real PFIX(NFIX,2), the fixed node positions. 网格中需要固定的点坐标，也就是一定需要出现在网格中的点。输出参数：real P(N,2), the node posit

4、ions. 网格点的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, vara

5、rgin: ) <= 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 includ

6、ed 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 (

7、 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

8、 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 scal

9、ar 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

10、, varargin: );% 计算需要的bar的长度，已经乘上了两个scale参数Fscale, sqrt ( sum(L.2) / sum(hbars.2) );% 具体可参考他们的paperL0 = 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)*

11、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;%计算移动梯度，相对边界dgrad

12、x = ( 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 2019.%p(1:nfix,1:2) = pfix;N =