好吧，我想出了一个办法。这有点令人费解，但我想我还是会在这里分享。我很肯定有一种更快的方法可以做到这一点，但由于每次模拟我只需要真正做几次，这并不是什么大不了的事。

% the line y = 0
xy = [linspace(0,1,100)', zeros(100,1)];
zz = zeros(100,1);

for i=1:100
    for j=1:nElements
        % check if the point xy(i,:) is in or on the jth element (I used barycentric coordinates to do this)
        if (xy(i,:) in element j)
            % compute zz(j) (in my case, using finite element solution on element j)
        end
    end
end

nElements是网格元素的总数。通过重心坐标确定该点是在元素中还是在其边界上非常简单。有关更多详细信息，请参见此处。

是的，它将检查所有元素，直到找到正确的元素，但它完成了工作。下面是我从沿的曲面图制作的示例。 $y=0$

如果有人有更快/更优雅的解决方案，请分享。

最简单的方法是使用内置griddata函数。一个问题是此函数不将三角形列表作为输入，而是使用 Delaunay 三角剖分进行插值。如果您的网格不是 Delaunay 三角剖分，则结果可能与直接从三角形插值略有不同。用法如下：

xv = linspace(xmin,xmax,nx); %//desired x points
yv = zeros(size(x)); %// desired y points

zv = griddata(x,y,z,xv,yv,'linear');

查看帮助文档以了解除线性插值之外的选项。

如果您需要确保使用网格而不是 Delaunay 三角剖分，则可能需要更复杂的东西。

__EDIT__
除非您有一个有效的搜索树来查找包含给定点的元素，否则我认为任何有效的解决方案都应该基于找到跨越所需平面（y = 0）的元素，然后仅在这些元素上进行插值而不是选择点，然后试图找出它们所在的元素。这将消除你的嵌套循环，这在 MATLAB 中可能非常慢。

一种这样的解决方案是计算三角形边缘与平面 y=0 的交点。如果您需要更高的分辨率并具有非线性元素，您还可以使用该元素内的有限元解决方案在这些三角形内插值。这可以按如下方式完成（并且可能比您发布的大量元素解决方案要快得多）：

%//loop over elements
xsave=[]; zsave=[];
for i=1:n_elements
    [x y z]=getVertices(elements,i);
    if min(y)*max(y)<=0 %//check if the element intersects y=0
        %// Find intersections between edges and y=0 plane
        %// The edges can be parameterized for s in [0, 1] as:
        %//   x(s) = x(1) + s*(x(2)-x(1));
        %//   y(s) = y(1) + s*(y(2)-y(1));
        %//   z(s) = z(1) + s*(z(2)-z(1));
        %// Then just solve for s setting y(s)=0
        s1 = -y(1)/(y(2)-y(1));
        s2 = -y(2)/(y(3)-y(2));
        s3 = -y(3)/(y(1)-y(3));

        %// Save the desired points:
        if s1>=0 && s1<=1
            xsave(end+1) = x(1) + s1*(x(2)-x(1));
            zsave(end+1) = z(1) + s1*(z(2)-z(1));
        end

        if s2>=0 && s2<=1
            xsave(end+1) = x(2) + s2*(x(3)-x(2));
            zsave(end+1) = z(2) + s2*(z(3)-z(2));
        end

        if s3>=0 && s3<=1
            xsave(end+1) = x(3) + s3*(x(1)-x(3));
            zsave(end+1) = z(3) + s3*(z(1)-z(3));
        end
    end
end

%//sort the resulting points, throw out duplicates and plot
[xsave,I,~] = unique(xsave);
zsave = zsave(I);

plot(xsave,zsave)

作为一个公平的警告，我没有测试上面的代码，可以稍微重写它以删除一些重复，但我认为它是最清楚的。

如果要获取特定位置的横截面，可以通过以下方式使用 MATLAB 函数 interp2：

w = interp2(xgrid,ygrid,zgrid,x,y0);
figure()
plot(x,w)

其中、和是用于绘制曲面的矩阵，$xgrid$$ygrid$$zgrid$$x$是轴之一的坐标，并且$y0$是您要截取横截面的其他轴上的坐标。

例如，如果要绘制横截面$y = 0$那么你必须写：

w = interp2(xgrid,ygrid,zgrid,x,0);