Reconstructing a surface from mesh points and function values (from COMSOL)

Hi all,

After having gotten more to grips with MATLAB/COMSOL LiveLink, I think I have found how to extract the mesh vertex coordinates and values of the dependent variable at those points.

I have stored these in the variables x, y, and z.

When I run plot3(x,y,z) I get the following (note that x and y are from a triangular mesh and not a rectangular grid):

I then run the following code:

x = x'; y = y'; z = z'; % Making sure these are column vectors
F = scatteredInterpolant(x,y,z);
[X,Y] = meshgrid(x,y);
Z = F(X,Y);
surf(X,Y,Z);

This seems to be a common way of constructing a surface plot from vectors that I have found online.

However, the surface plot I get looks awful:

This is the plot that the COMSOL GUI provides.

Does anyone know how I can fix it so that the surface plot produced by MATLAB (a) is much smoother, and (b) doesn't leave jagged edges at the boundary of the domain?

Answer 1

[Bjorn Gustavsson]

Are you sure that you need to use scatteredInterpolant at all? The points in the first plot seem to be on a regular grid. If that is the case then you're already done. If they are on a regular grid but given in some obscure ordering, you should be able to sort them into a correct order, hopefully using reshape and permute.

If your x and y are 2-D arrays you can check that things are OK byt looking at them:

x = x'; y = y'; z = z';
x = reshape(x,[14 14]);
y = reshape(y,[14 14]);
z = reshape(z,[14 14]);
subplot(1,3,1)
imagesc(x)
subplot(1,3,2)
imagesc(y)
subplot(1,3,3)
imagesc(z)

It is definitely your call to meshgrid with repeated values in x and y that generates peculiarly ordered X or Y. Test this snippet just after your meshgrid call.

figure
subplot(1,2,1)
imagesc(X)
subplot(1,2,2)
imagesc(Y)

HTH

[Oliver Bond]

Do you mean "regular grid" as in a rectangular grid with evenly-shaped points? If so, then the points I have plotted do not lie on this sort of grid; they are vertices on a triangular mesh.

Also, x, y and z are each P-by-1 arrays, where P is the number of points in the domain.

[Bjorn Gustavsson]

Then generate X and Y something like this:

[X,Y] = meshgrid(0:0.05:1,0:0.05:1);

Because what you get with your meshgrid call cannot be pretty?

[Bjorn Gustavsson]

Or use trisurf, if you only need to display the surface COMSOL produced.

[Oliver Bond]

I've tried this, but it doesn't seem to work with the triangulation matrix that COMSOL provides.

[Bjorn Gustavsson]

Read the help for trisurf to see the example:

    tri = delaunay(x,y);
    z = peaks(15);
    trisurf(tri,x,y,z)

Then adapt this to your case:

    tri = delaunay(x,y);
    trisurf(tri,x,y,z)

You don't need to use the triangulation that COMSOL returns, it might very well use a completely different structure. You have enough to make do with the x, y and z arrays. My previous correction for how to create X and Y should also work if done in your script.

HTH

[Bjorn Gustavsson]

Did you get to a solution of your problem?

[Oliver Bond]

Sorry for the profusely long delay in replying. I have found something that very nearly does what I want. I discovered the surfir script that somebody created. It almost works in my case - however there are still a couple of pairs of points that are being incorrectly interpolated between. The domain is T-shaped (as I have illustrated with the black lines that I've drawn on myself), but there are two triangular portions that should not be there because there are no z-values that have points inside those domains as their corresponding x- and y-values.

Since it has been a long time since I asked this question, I will create a new one.

[Bjorn Gustavsson]

If you look into surfir you will see that it is an explicit use of the suggestions I've made:

% Below I will mark my comments BG

% h=surfir(x,y,z,s,opt)
%
% Surface plot based on irregularly spaced data points.
% x,y,z are coordinate vectors of the same length.
% s (default 0.5) is the shrink factor in the Matlab function 'boundary' 
%   (type 'help boundary' for details)
% opt is a list of options to be passed to the 'trisurf' function 
%   (type 'help trisurf' for details)
function h=surfir(x,y,z,s,opt)
%% default parameters
if (nargin<4)||isempty(s)                                                   % no shrink factor provided
   s=0.5;                                                                   % default value
end
if nargin<5                                                                 % no options provided
   opt={'FaceColor','interp','edgecolor','none'};                           % default
end
%% Remove duplicate data points
[xy,ind] = unique([x,y],'rows');
z=z(ind);
x=xy(:,1);
y=xy(:,2);
%% triangulate data points
% BG - here the dealaunay-triangulation I suggested is done:
dt = delaunayTriangulation(x, y);                                           % Delaunay triangulation
x=dt.Points(:,1);
y=dt.Points(:,2);
%% find enclosing boundary
k=boundary(x,y,s);                                                          % define boundary enclosing all data points
c=[k(1:end-1),k(2:end)];                                                    % constraints
dt.Constraints=c;                                                           % add constraints
io = dt.isInterior();                                                       % triangles which are in the interior of domain
tri=dt.ConnectivityList;                                                    % read triangles
tri=tri(io,:);                                                              % use only triangles in the interior
%% plot
% BG - Here the function uses trisurf to make the plot
h=trisurf(tri,x,y,z,z,opt{:});                                              % plot triangles and interpolate colors

The triangulation will return a triangulation that covers the convex hull around your points. You would somehow find a way to remove those triangles that lie outside of the boundary you like.

Perhaps you can figure out which triangles that are and what corner-points they have and simply remove those triangles. You could do that programming work yourself, or test different values for the fourth input-argument to surfir.