# list plot over arbitrary domain

I have a weird problem. I present here a minimal example. I wish to ListPlot3D over a non-rectangular domain. For that I build a grid:

grid = Flatten[Outer[If[Abs[#2] > #1, {#1, #2}, (## &[])] &, Range[0, 10],
Range[-10, 10]], 1];


To make sure I plot the grid points:

ListPlot[grid]


I now map a function over my grid and plot it:

f = {#[[1]], #[[2]], (#[[1]] + #[[2]])^2} &;
ListPlot3D[ParallelMap[f, grid]]


As far as I understand I should see a plot of my function over the domain which I specified, but instead I see this:

I see the same issue with the contour plot:

Can someone elaborate on what's going on? Why do the plot function continue my domain. Based on which data is the continuation done? I see the same behaviour even if I set InterpolationOrder->0? How can I force MMA to stick to my domain?

## Update

I noticed the following: If I add RegionFunction alone nothing happens:

ListPlot3D[ParallelMap[f, grid], RegionFunction -> Function[{x, y, z}, Abs[y] > x]]


But if I add also InterpolationOrder->0:

ListPlot3D[ParallelMap[f, grid],
RegionFunction -> Function[{x, y, z}, Abs[y] > x],
InterpolationOrder -> 0]


Somehow the interpolation interferes with the region function.

## Update 2

I found a partial solution eventually here: How to disable interpolation of ListPlot3D across points outside of the domain?

As suggested there the interpolation causes the problem, and therefore limiting the amount of used plot points (together with RegionFunction) resolves the issue:

ListPlot3D[ParallelMap[f, grid],
RegionFunction -> Function[{x, y, z}, Abs[y] > x],
MaxPlotPoints -> 110]


However this solution is only partially successful. In my real data, as I want to have a rather good approximation I use more than 2000 points for the dataset. When limiting with MaxPlotPoints -> 2000, My computer just evaluates forever. Thus I currently use MaxPlotPoints -> 200, which takes around two minutes. A second solution is to use MaxPlotPoints -> 2000, with InrerpolationOrder->0, which evaluates quickly, but the points remain disconnected.

Any workaround to have a complete solution? For instance using the RegionFunction option alone would suffice if I could use Plot3D instead of ListPlot3D.

## Update 3

The issue has been resolved. However I add this in case someone else who encounters the issue sees this. The source of the problem is the behaviour of interpolation on concave domains. Apparently when ListPlot3D interpolates it just does interpolation between each two nearest neighbors, regardless of the domain. Thus when plotting over a concave domain, the interpolation fills in parts which aren't part of the domain, even if it has been clearly defined using the RegionFunction option (At least at MMA 9). Personally I think this is a bug. The suggested workaround in the accepted answer works well, as long as you don't use a ColorFunction. If you wish to use such a color function perhaps the best solution would be to define a color function, and use it for both plots instead of letting it be resolved automatically for each graph.

• What is the domain?, because what you have is a grid of points. p.s. take a look at ListPlot3D[ParallelMap[f, grid], Mesh -> All] – Kuba Mar 6 '17 at 10:52
• The domain is the list points. I'm using ListPlot. As you can see my domain is a "square" minus a "triangle", or the list of all points for which Abs[y]>x. Doesn't that suffice as a domain? – Yair M Mar 6 '17 at 10:55
• Possible duplicate of How to plot customly-shaped set of 3D points each connected to its neighbors? – corey979 Mar 6 '17 at 10:56

Perhaps the domain can be subdivided into convex pieces.

There a possible problem with Abs[y] > x on discrete data: Is the region disconnected at {0, 0} or not? I assume it is, otherwise the region would be Abs[y] >= x. I raise the issue because in the OP, the plots look connected and there's no mention whether that is inaccurate or desired. Anyway, the plot can be fixed by changing the inequality.

Clear[f]; (* if f[] is vectorized; otherwise map as in OP *)
f[g_] := With[{x = g[[All, 1]], y = g[[All, 2]]},
Transpose@{x, y, (x + y)^2}];

top = UnitStep[grid[[All, 2]]];
bottom = UnitStep[-grid[[All, 2]]];
surf = f[grid];

Show[
ListPlot3D[Pick[surf, top, 1]],
ListPlot3D[Pick[surf, bottom, 1]],
PlotRange -> All
]


Here is a grid with 6000+ points, which takes 0.3 sec. to plot:

grid = Flatten[
Outer[If[Abs[#2] > #1, {#1, #2}, (## &[])] &,
Range[0., 10., 1./8], Range[-10., 10., 1./8]], 1];
Length@grid
(*  6480  *)


as a workaround, create an interpolation function and plot that:

grid = Flatten[
Outer[If[Abs[#2] > #1, {#1, #2}, (## &[])] &, Range[0, 10],
Range[-10, 10]], 1];
f = {#[[1]], #[[2]], (#[[1]] + #[[2]])^2} &;
int = Interpolation[ParallelMap[f, grid], InterpolationOrder -> 1]
Plot3D[int[x, y], {x, 0, 10}, {y, -10, 10},
RegionFunction -> Function[{x, y, z}, Abs[y] > x]]


FWIW ListPlot3D respects the RegionFunction if you feed it a structured grid.

grid = Flatten[Table[{x, y}, {x, 0, 10}, {y, -10, 10}], 1];
f = {#[[1]], #[[2]],
If[Abs[#[[2]]] >= #[[1]], (#[[1]] + #[[2]])^2, 0, 0]} &;
ListPlot3D[ParallelMap[f, grid],
RegionFunction -> Function[{x, y, z}, Abs[y] > x]]


not so good because of the edge effects, but maybe it gives a clue whats going on.

• How is what you do with the ListPlot different from what I did above? For some reason when I sued the RegionFunction it didn't work for me as I expected, unless I disabled the interpolation via InterpolationOrder->0. You can see that in my question above. Yet your solution seems to do quite the same thing with a slight change in the definition of the grid and function, but in this way ListPlot respects the RegionFunction. What is the source of the difference? – Yair M Mar 6 '17 at 16:15
• What do you mean by a structured grid, and why isn't mine structured? – Yair M Mar 6 '17 at 16:16
• In my second example I created a complete rectangular grid with dummy values (0) at the missing locations. (That's what I meant by structured, maybe not the best word choice). For this simple example you could fix the edge problem by fully populating the grid with correct values, but I assume the real function is undefined or too expensive to compute over the whole domain. – george2079 Mar 6 '17 at 16:32

Another workaround.

The RegionFunction is required because otherwise Mathematica has no way to know which empty areas are outside the domain and which are simply gaps between points.

Points are tested with the RegionFunction before being included in the mesh. In your example every specified point is inside the region, so the mesh joins them all up. If you add a handful of points which are outside the region, this will stop Mathematica from meshing across the gap and force it to subdivide to find the region boundary.

ListPlot3D[Map[f, grid ~Join~ Table[{x, 0}, {x, 10}]],
RegionFunction -> Function[{x, y, z}, Abs[y] > x]]