Creating surface plot over non-rectangular domain - getting jaggedness at boundaries (and not resolving boundary layers)

I have a set of x, y, and z values, where x and y are on a non-rectangular domain. These are the output of the solution of a problem (that I solved using COMSOL). In this case the z-values are only defined on a subset of the domain of x and y, and the array of z-values contains NaN entries for those particular x and y values.

I can't really use [X,Y] = meshgrid(x,y), because (a) the x and y values are vertices from a non-rectangular mesh, (b) the vertices are not necessarily ordered, and (c) the operation is very slow due to taking up a lot of memory.

I have found a workaround which involves using scatteredInterpolant:

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

F = scatteredInterpolant(x,y,z)

xR = linspace(min(x),max(x),101);

yR = linspace(min(y),max(y),101);

[XR,YR] = meshgrid(xR,yR);

ZR = F(XR,YR);

surf(XR,YR,ZR)

This gives the following output for a surface plotted over a T-shaped domain:

or, in a three-dimensional projection:

This is very close to the result I want, but not quite there yet. This is because (a) there are chunks missing near the corners, and (b) there is still some jaggedness on the boundaries. I think this is probably because the points on the boundary of the domain are not necessarily included in xR and yR.

As an attempt at a workaround, I have tried repeating my code with the following modification:

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

% xB and yB are coordinates of points that lie on the boundary of the domain

xR = union(linspace(min(x),max(x),101),xB);

yR = union(linspace(min(y),max(y),101),yB);

However, this still doesn't remove the jaggedness at the edges. In fact, I also get the error:

% Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

I have also tried changing the method that scatteredInterpolant uses, and it looks even worse. Furthermore, the value of the function being plotted is supposed to be zero at the boundaries - the solution contains so-called "boundary layers" near the boundary, where the mesh of the solution is very finely resolved. I gather that the linspace command is not taking this into account.

Could somebody suggest a workaround? Is there any analogue to linspace which uses variable spacing?

Answer 1

[SANKALP DEV]

Hi Oliver,

I understand that you are facing challenges with missing chunks near the corners and jaggedness on the boundaries when using the ‘scatteredInterpolant’ function.

To address these issues, I recommend you incorporate variable spacing in your ‘xR’ and ‘yR’ arrays, instead of using ‘linspace’.

This can help to capture the boundary layers and improve the interpolation near the edges of the domain. One effective approach is to utilize a clustering algorithm, such as ‘k-means’ clustering, to determine the variable spacing based on the density of points.

You can modify your code to incorporate variable spacing using ‘k-means’ clustering as follows:

% x, y, and z are pre-defined 1-by-P arrays, where P is the number of points in the domain. 
F = scatteredInterpolant(x, y, z); 
% Determine the number of points you want in the xR and yR arrays 
numPoints = 101; 
% Create a matrix of [x, y] coordinates 
XY = [x', y']; 
% Perform k-means clustering to determine the variable spacing 
[idx, ~] = kmeans(XY, numPoints); 
% Extract the unique cluster centers as the variable spacing 
clusterCenters = unique(XY(idx, :), 'rows'); 
xR = clusterCenters(:, 1)'; 
yR = clusterCenters(:, 2)'; 
% Evaluate the scatteredInterpolant on the variable spacing grid 
ZR = F(xR, yR); 
% Plot the surface 
surf(xR, yR, ZR); 

In the above modified code, I have used the ‘kmeans’ function to cluster the [x, y] coordinates into ‘numPoints’ clusters. The unique cluster centres are then extracted and used as the variable spacing for ‘xR’ and ‘yR’. This approach should provide better resolution near the boundaries and accurately capture the boundary layers.

To know more about ‘kmeans’ you can refer to the following documentation:

k-means clustering - MATLAB kmeans (mathworks.com)

Hope this helps.