非结构化网格内等值线绘制

参考自 Implicit Surface Intersections--Mike Garrity

示例

绘制函数 (y-xcdot mathrm{tan}(z) = 0) 与圆柱 (x^2 + y^2 = 1) 的交线。

一般绘制交线的方法有以下几种：

1. 寻找解析解，
2. 寻找交点，并且使用追踪方法构造曲线，
3. 计算网格中交点，绘制曲面。

本文使用的就是第三种方法，主要过程包括：

1. 寻找等值线所在单元，
2. 寻找单元插值节点，
3. 每个单元内将插值点逐段连接。

首先构造曲面 (y-xcdot mathrm{tan}(z) = 0) 并且找出所有满足 (x^2 + y^2 > 1) 的节点，将结果保存在变量mask内。

[x3, y3,z3] = meshgrid(linspace(-1.25, 1.25, 150), ...
                       linspace(-1.25, 1.25, 150), ...
                       linspace(0, 2*pi, 150));
f2 = y3-x3.*tan(z3);
hel = isosurface(x3, y3, z3, f2, 0);

%% Draw faces
patch(hel,'FaceColor',[1 .5 0],'EdgeColor','none');
view(3)
camlight
ax = gca;
set(ax,'XLim',[-inf inf],'YLim',[-inf inf],'ZLim',[-inf inf],'DataAspectRatio',[1 1 1])
box on
xlabel('x(mm)')
ylabel('y(mm)')
zlabel('Lambda[rad]')

%% Find vertex
r = sqrt(hel.vertices(:,1).^2 + hel.vertices(:,2).^2);
cylrad = 1;
mask = r>cylrad;

下面通过统计每个三角形内满足要求的节点个数，得到插值线所在单元序号cross及其三个节点编号crossing_tris

outcount = sum(mask(hel.faces),2);
cross = (outcount == 2) | (outcount == 1);
crossing_tris = hel.faces(cross,:);

下一步需要把这些单元转化为线段以表示等值线。首先把这些把每个单元内节点进行整理，将插值点所在的两个边交点排在三个节点的第一个位置，具体过程如下:

1. 将只有一个节点满足 (x^2 + y^2 > 1) 的单元内，另外两个节点标记为1，此节点标记为0。这样所有单元内都只有唯一一个标记为0的节点，此节点就是插值点所在两个边的交点。

out_vert = mask(crossing_tris);
flip = sum(out_vert,2) == 1;
out_vert(flip,:) = 1-out_vert(flip,:);

2. 将三个节点交换顺序并将序号储存在变量overt中，交点在前，后面两个仍按顺序（逆时针或顺时针）排列。

ntri = size(out_vert,1);
overt = zeros(ntri,3);
for i=1:ntri
    v1i = find(~out_vert(i,:));
    v2i = 1 + mod(v1i,3);
    v3i = 1 + mod(v1i+1,3);
    overt(i,:) = crossing_tris(i,[v1i v2i v3i]);
end

3. 计算每个单元内插值节点坐标，其中u，v分别为插值点在两边上所占比率。

u = (cylrad - r(overt(:,1))) ./ (r(overt(:,2)) - r(overt(:,1)));
v = (cylrad - r(overt(:,1))) ./ (r(overt(:,3)) - r(overt(:,1)));

uverts = repmat((1-u),[1 3]).*hel.vertices(overt(:,1),:) + repmat(u,[1 3]).*hel.vertices(overt(:,2),:);
vverts = repmat((1-v),[1 3]).*hel.vertices(overt(:,1),:) + repmat(v,[1 3]).*hel.vertices(overt(:,3),:);

最终将获得的插值节点坐标储存在变量x，y与z中，便可得到最终两曲面交线。

x = nan(2,ntri);
x(1,:) = uverts(:,1)';
x(2,:) = vverts(:,1)';
y = nan(3,ntri);
y(1,:) = uverts(:,2)';
y(2,:) = vverts(:,2)';
z = nan(3,ntri);
z(1,:) = uverts(:,3)';
z(2,:) = vverts(:,3)';

h = line(x(:),y(:),z(:));