基于MATLAB的圆形等边三角形网格与正方形均等划分实现

一、圆形等边三角形网格划分

核心思路：通过调整初始节点分布与Delaunay三角剖分，生成近似等边三角形的圆形网格。关键步骤包括非均匀初始节点生成、边界约束和迭代优化。

1. 实现代码

%% 参数设置
h0 = 0.1;          % 初始网格间距
iteration_max = 200; % 最大迭代次数
geps = 0.001*h0;   % 边界容差

%% 定义圆形区域
fd = @(p) sqrt(p(:,1).^2 + p(:,2).^2) - 1; % 圆心(0,0)半径1
fh = @(p) 0.1;    % 均匀网格尺寸函数

%% 生成初始节点（等边三角形布局）
[x, y] = meshgrid(linspace(-1,1,20), linspace(-1,1,20));
x(2:2:end,:) = x(2:2:end,:) + h0/2; % 偶数行右移半格
p = [x(:), y(:)];

%% 过滤边界外节点
p = p(feval(fd, p) <= geps, :);

%% 添加固定边界点（可选）
pfix = [0,1; 1,0; 0,-1; -1,0]; % 四个方向的控制点
p = [p; pfix];

%% Delaunay三角剖分
t = delaunayn(p);

%% 迭代优化节点位置
for iter = 1:iteration_max
    % 计算三角形重心
    pmid = (p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:)) / 3;
    
    % 移除边界外的三角形
    t = t(feval(fd, pmid) <= -geps, :);
    
    % 生成边集合
    bars = unique(sort([t(:,[1,2]); t(:,[1,3]); t(:,[2,3])], 2), 'rows');
    
    % 计算节点受力并更新位置
    barvec = p(bars(:,1),:) - p(bars(:,2),:);
    L = sqrt(sum(barvec.^2, 2));
    F = max(0.1*L - L, 0); % 弹性力模型
    Ftot = sparse(bars(:,[1,2]), bars(:,[2,1]), F, size(p,1), 2);
    p = p + 0.2 * full(Ftot); % 更新节点位置
    
    % 将固定点拉回边界
    p(pfix(:,1), :) = pfix(:,2);
end

%% 可视化
figure;
trisurf(t, p(:,1), p(:,2), zeros(size(t,1),1));
axis equal;
title('圆形等边三角形网格');

2. 关键参数说明

• h0：初始网格间距，控制整体密度。
• fh：尺寸函数，fh=常数表示均匀网格，fh=距离函数可实现非均匀加密。
• pfix：固定点坐标，用于约束边界节点位置。

---

二、正方形均等划分网格

核心思路：通过规则网格划分生成正方形网格，结合Delaunay三角剖分形成四边形或三角形网格。

1. 实现代码

%% 参数设置
h0 = 0.1;          % 初始网格间距
iteration_max = 50;% 最大迭代次数

%% 定义正方形区域
fd = @(p) drectangle(p, -1, 1, -1, 1); % 正方形边界
fh = @(p) 0.1;    % 均匀网格尺寸函数

%% 生成初始节点
[x, y] = meshgrid(linspace(-1,1,20), linspace(-1,1,20));
p = [x(:), y(:)];

%% Delaunay三角剖分
t = delaunayn(p);

%% 迭代优化（可选）
for iter = 1:iteration_max
    pmid = (p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:)) / 3;
    t = t(feval(fd, pmid) <= -geps, :);
    bars = unique(sort([t(:,[1,2]); t(:,[1,3]); t(:,[2,3])], 2), 'rows');
    barvec = p(bars(:,1),:) - p(bars(:,2),:);
    L = sqrt(sum(barvec.^2, 2));
    F = max(0.1*L - L, 0);
    Ftot = sparse(bars(:,[1,2]), bars(:,[2,1]), F, size(p,1), 2);
    p = p + 0.1 * full(Ftot);
end

%% 可视化（四边形网格）
figure;
quads = [t(:,[1,2,3,4])]; % 提取四边形单元（需补充第四点）
quads(:,4) = t(:,1);      % 闭合四边形
patch('Faces', quads, 'Vertices', p, 'FaceColor', 'cyan', 'EdgeColor', 'k');
axis equal;
title('正方形均等四边形网格');

2. 网格类型选择

• 三角形网格：直接使用Delaunay剖分结果。
• 四边形网格：需补充第四点形成闭合单元（如quads(:,4) = t(:,1)）。

---

三、网格质量优化技巧

1. 尺寸函数控制

   通过fh函数实现非均匀加密，例如在圆形边界附近加密：

   fh = @(p) 0.05 + 0.3*sqrt(p(:,1).^2 + p(:,2).^2); % 边界处网格更密
   

2. 节点投影修正

   迭代中将边界外节点拉回几何边界：

   d = feval(fd, p);
   p(d > 0, :) = p(d > 0, :) - d(d > 0, :) .* [1,1] * 0.1; % 梯度投影
   

3. 固定点约束

   在关键位置（如孔洞中心）添加固定点：

   pfix = [0,0]; % 圆心固定点
   

---

四、应用案例对比

网格类型	圆形等边三角形	正方形均等四边形
适用场景	圆形边界问题（如应力分析）	正方形区域模拟（如电磁场计算）
节点分布特点	边缘密集，中心稀疏	均匀分布
计算效率	边界单元更小，收敛更快	计算均匀，适合稳态问题
MATLAB函数	delaunayn+ 迭代优化	delaunayn+ 四边形单元重构

---

五、常见问题解决

1. 网格畸变
   • 原因：初始节点分布不均匀或尺寸函数不合理。
   • 解决：调整h0和fh，增加迭代次数。
2. 边界节点穿透
   • 原因：投影步长过大。
   • 解决：减小投影系数（如0.1改为0.05）。
3. 固定点失效
   • 原因：固定点未正确添加到pfix中。
   • 解决：检查pfix坐标是否在几何域内。

参考代码 网格划分程序 www.youwenfan.com/contentcnp/96340.html

六、扩展应用

• 多孔洞网格：通过dpoly定义多边形区域，结合ddiff实现孔洞。

  fd = @(p) ddiff(drectangle(p,-1,1,-1,1), dcircle(p,0,0,0.3));
  

• 自适应网格：根据场量（如温度梯度）动态调整fh函数。

---

参考文献

[1] Per-Olof Persson. DistMesh: A Simple Mesh Generator in MATLAB. persson.berkeley.edu/distmesh/

[2] MATLAB官方文档：delaunayTriangulation类使用指南