Obtaining a 2D profile from a 3D surface

Question

I am plotting a saturation field (as the solution of a finite element problem, using trisurf in MATLAB). Ideally, I would like to obtain the two-dimensional profile of this surface along $y=0$. See Figure 4 of this paper for an example. I am not really familiar with MATLAB's plotting routines but I noticed there's a feature to grab a data point on the surface using the cursor, so I'm hoping this is a built-in feature.

Answer 1

Well I've figured out a way to do it. It's a bit convoluted but I figured I'd share it here anyway. I'm positive there is a faster way to do this but since I only need to really do this a couple of times per simulation, it's not a huge deal.

% 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 is the total number of mesh elements. Determining whether the point is in the element or on its boundary is very simple via barycentric coordinates. See here for more details.

Yes, it's going to check over all the elements until it finds the correct one but it gets the job done. Below is an example I made taken from a surface plot along $y=0$.

If anyone has a faster/more elegant solution, please do share.

Answer 2

The simplest method would be to use the built in griddata function. The one issue is that this function does not take a triangle list as input and instead uses a Delaunay triangulation for interpolation. If your mesh is not a Delaunay triangulation the results may be slightly different than directly interpolating from your triangles. Usage is follows:

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

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

Take a look at the help documentation for options other than linear interpolation.

If you need to be sure to use your mesh rather than a Delaunay triangulation, you may need something more complex.

__EDIT__
Unless you have an efficient search tree to find which element contains a given point, I think any efficient solution should probably be based on finding the elements that straddle the desired plane (y=0) and then interpolating only on those elements rather than choosing points and then trying to figure out which element they're in. This will eliminate your nested loop which can be quite slow in MATLAB.

One such solution would be to compute the intersections of your triangle edges with the plane y=0. If you need more resolution and have nonlinear elements you could also interpolate within those triangles using your finite element solution within that element. This can be done as follows (and could be much faster than your posted solution for large numbers of elements):

%//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)

As a fair warning, I haven't tested the above code and it could be re-written slightly to remove some of the duplication, but I think it's most clear as is.

Answer 3

If you want to obtain the cross-section at a specific location you can use the MATLAB function interp2 in this way:

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

Where $xgrid$, $ygrid$ and $zgrid$ are the matrices you used to plot your surface, $x$ is the coordinate of one of the axis and $y0$ is the coordinate on the other axesat which you want to take the cross-section.

For example, if you want to plot the cross-section for $y = 0$ then you have to write:

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