This is slightly more/different than OP requested, but why not?
g = Plot3D[Sqrt[16 - x^2 - y^2], {x, -4, 4}, {y, -4, 4},
Mesh -> None, PlotStyle -> Opacity[.7]]
Next define a function that gets us all the rectangles we need, wraps them into a Head
of Polygon
and overlays them onto the original plot:
addgrids[g_Graphics3D] :=
Module[{ranges = {Floor@#1, Ceiling@#2} & @@@ PlotRange@g, rects},
rects =
Function[{perm},
Map[Function[{points}, Prepend[points, #]],
Tuples[Rest@#][[{1, 2, 4, 3}]] &@perm] & /@
Range @@ First@perm] /@
(RotateRight[ranges, #] & /@ Range[3]);
rects =
MapIndexed[Function[{list, amount},
Map[RotateLeft[#, First@amount] &, list, {2}]], rects];
Show[{g,
Graphics3D[{Opacity[.2], Polygon@Flatten[rects, 1]}]
}]
]
Then simply post-process the plot with this function:
g // addgrids
While I'm at it, here's a similar approach for adding grid-lines, rather than grid-surfaces:
addlines[g_Graphics3D] :=
Module[{ranges = {Floor@#1, Ceiling@#2} & @@@ PlotRange@g, lines},
lines =
Function[{perm},
Flatten[
Block[{i, j, Table = Table}, ClearAttributes[Table, HoldAll];
(Table[{{First@First@perm, i, j}, {Last@First@perm, i, j}},
{i, Sequence @@ #1}, {j, Sequence @@ #2}] & @@ (Rest@
perm))], 1]] /@
(RotateRight[ranges, #] & /@ Range[3]);
lines =
MapIndexed[Function[{list, amount},
Map[RotateLeft[#, First@amount] &, list, {2}]], lines];
Show[{g,
Graphics3D[Line@Flatten[lines, 1]]
}]
]
g // addlines