3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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$. enter image description here

enter image description here

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

$\endgroup$
1
  • $\begingroup$ See my updated answer. I think you're probably much better off looking for the elements that straddle the y=0 plane rather than choosing points and then looking for the elements that contain them. $\endgroup$ – Doug Lipinski Mar 25 '14 at 21:17
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ I didn't know this function existed, but unfortunately it didn't work for my data :/ $\endgroup$ – Justin Dong Mar 25 '14 at 7:18
0
$\begingroup$

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);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.