How to plot a surface z=f(x,y) with a hole?

I want to make a plot similar to the one at left (the goal is to illustrate an interpolation: the plot at right).

I did this plot with a constrained Delaunay tessellation (with the delaunay package) applied to the points (x,y) given on a grid with a hole and then I constructed a 3D rgl mesh. I'm wondering whether there exists a more straightforward method. I tried rgl::perps3d.deldir but it fills the hole and even if it is possible to remove this part, this way is really too slow.

Here is the grid with the hole and the z values:

x <- y <- seq(-3, 3, by = 0.1)
Grid <- expand.grid(X = x, Y = y)
toremove <- with(Grid, X >= -0.5 & X <= 1 & Y >= -0.5 & Y <= 1)
Grid <- Grid[!toremove, ]
z <- with(Grid, dt(X, df = 1) * dt(Y, df = 1))

I also want to be able to add another plot to the first one (the black interpolated part).

I'd just set the values in the hole to NA, e.g.

library(rgl)
x <- y <- seq(-3, 3, by = 0.1)
Grid <- expand.grid(X = x, Y = y)
toremove <- with(Grid, X >= -0.5 & X <= 1 & Y >= -0.5 & Y <= 1)
z <- with(Grid, dt(X, df = 1) * dt(Y, df = 1))
z[toremove] <- NA
persp3d(x, y, z, col = "red")

^{Created on 2023-03-08 with reprex v2.0.2}

If you have a situation where some x, y pairs actually cause errors in the z calculation, then it's a little trickier: put the real function in a wrapper that returns NA for the hole.