lianxi.rar_Mesh_区域划分_网格划分

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