list plot over arbitrary domain

Question

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.

Answer 1

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  *)

Answer 2

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.

Answer 3

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]]