LevelPlot3D

By Gianluca Gorni, University of Udine, Italy.

Last revised: June 7, 2004.

Description: Plotting the level sets of a function in 3D, slicing a surface with parallel planes, drawing the intersection between a surface in three dimensions and a half-space.

Instructions: load the package by clicking on the button below, then evaluate the examples.

Button[
      StyleForm["Click here to load package!", 
        FontColor -> RGBColor[1, 0, 0], FontSize -> 14], 
      ButtonFunction :> 
        FrontEndExecute[FrontEndToken["EvaluateInitialization"]], 
      ButtonEvaluator -> None, Active -> True]

Mathematica package (to load into the kernel)

preamble

\(\(BeginPackage["LevelPlot3D`", "Utilities`FilterOptions`"];\)\)

Unprotect["`*"];
ClearAll["`*"];

\(\(LevelPlot3D::usage = 
"LevelPlot3D[f,{x,x0,x1},{y,y0,y1}] draws a 3D plot of the function f with the level sets emphasized. You can control the number (or precise position) of level sets to be drawn with the option Contours. Example: LevelPlot3D[x^2+y^2,{x,-1,1},{y,-1,1},PlotPoints->11,Contours->5}].";\)\)

SliceSurface::usage = "SliceSurface[3DSurface, \
LinearForm,variables] draws on the surface the level sets of the LinearForm. \
You can control how many level sets (or their precise values) with the option \
Slices . Example: SliceSurface[Plot3D[Sin[x*y],{x,0,Pi},{y,0,Pi}],x+y,x,y,z].\
";

LinearClip::usage = "LinearClip[3DSurface, \
LinearInequality,variables] shows the part of the surface that lies in the \
half-plane defined by the linear inequality. Example: \
LinearClip[Plot3D[Sin[x*y],{x,0,Pi},{y,0,Pi}],x+y>=8z,x,y,z].";

PrepareForSlices::usage = "PrepareForSlices is an option for \
LinearClip, specifying whether to keep both the two parts of the surface when \
cutting with a plane. It can take the values True or False.";

Slices::usage = "Slices is an option for SliceSurface, specifying \
how many level sets (or their precise values) are to be drawn. It can take \
either an integer value or a list of real values.";

SliceEdgeStyle::usage = "SliceEdgeStyle is an option for \
LevelPlot3D and SliceSurface, specifying how to draw the edge of the slices, \
i.e., of the level sets.";

ShowProjectedLevelSets::usage = "ShowProjectedLevelSets is an \
option for LevelPlot3D. When set to True, the level sets will be shown both \
on the 3D surface and as projected onto a horizontal plane. The z-level of \
the plane is controlled with ProjectAt. The default is \
ShowProjectedLevelSets->False";

ProjectAt::usage = "ProjectAt is an option for LevelPlot3D. \
ProjectAt->realNumber tells at what z-level to draw the projected level sets. \
The default is ProjectAt->0";

SetAttributes[LinearClip,  HoldRest];

Options[LinearClip] = {PrepareForSlices -> False};

SetAttributes[SliceSurface,  HoldRest];

Options[SliceSurface] = {Slices -> 11, 
          MeshStyle -> {EdgeForm[]}, 
          SliceEdgeStyle -> {}};

SetAttributes[LevelPlot3D, HoldAll];

Options[LevelPlot3D] = {Contours -> 11, 
          MeshStyle -> {EdgeForm[]}, SliceEdgeStyle -> {}, 
          ShowProjectedLevelSets -> False, 
          ProjectAt -> 0};

Begin["`Private`"];

LinearClip

LinearClip[obj_SurfaceGraphics, other___] := 
        LinearClip[Graphics3D[obj], other];

LinearClip[obj_Graphics3D, 
            rel_ /; MemberQ[{Less, LessEqual, Greater, GreaterEqual}, 
                Head[rel]], var__Symbol, 
            opts___Rule] /; (PolynomialQ[Subtract @@ rel, {var}] && 
              Max[(Exponent[Expand[Subtract @@ rel], #1] &) /@ {var}] === 
                1) := Module[{ratObj, strictRel, weakRel, strictNegRel, 
            weakNegRel, eqRel, interiorQ, weakInteriorQ, exteriorQ, 
            weakExteriorQ, onBoundaryQ, interiorQList, weakInteriorQList, 
            exteriorQList, weakExteriorQList, borderTouchingQ, toThrowQ, 
            winnowed, inters, t, clipLine, clipContour, beyondClipLine, 
            beyondClipContour, prepareQ, 
            boundaryLine}, 
ratObj = 
            Rationalize[obj, 1/10^6]; 
strictRel = 
            Rationalize[
              rel /. {LessEqual -> Less, 
                  GreaterEqual -> Greater}, 
              1/10^6]; 
weakRel = 
            strictRel /. {Less -> LessEqual, 
                Greater -> 
                  GreaterEqual}; 
strictNegRel = 
            strictRel /. {Less -> Greater, 
                Greater -> Less}; 
weakNegRel = 
            strictNegRel /. {Less -> LessEqual, 
                Greater -> GreaterEqual}; 
eqRel = 
            strictRel /. {Less -> Equal, 
                Greater -> Equal}; 
interiorQ[
              pt_List] := 
            strictRel /. Thread[{var} -> 
                  pt]; 
weakInteriorQ[pt_List] := 
            weakRel /. Thread[{var} -> 
                  pt]; 
exteriorQ[pt_List] := 
            strictNegRel /. 
              Thread[{var} -> pt]; 
weakExteriorQ[
              pt_List] := 
            weakNegRel /. Thread[{var} -> pt]; 
          onBoundaryQ[pt_List] := 
            eqRel /. Thread[{var} -> 
                  pt]; 
interiorQList[pts_List] := 
            interiorQ /@ pts; 
weakInteriorQList[pts_List] := 
            weakInteriorQ /@ pts; 
exteriorQList[pts_List] := 
            exteriorQ /@ pts; 
weakExteriorQList[pts_List] := 
            weakExteriorQ /@ pts; 
          borderTouchingQ[pts_List] := 
            And[Or @@ weakInteriorQList[pts], Or @@ weakExteriorQList[pts]]; 
          toThrowQ[Polygon[pts_List]] := And @@ exteriorQList[pts]; 
          winnowed = 
            obj /. pol_Polygon?toThrowQ :> 
                Sequence @@ {}; 
          inters[pt1_List, pt2_List] := 
            t pt1 + (1 - t) pt2 /. Solve[
                  eqRel /. Thread[{var} -> 
                        t*pt1 + (1 - t) pt2], t][[1]]; 
          clipLine[pt1_List?exteriorQ, pt2_List?exteriorQ] := Sequence @@ {}; 
          clipLine[pt1_List?interiorQ, pt2_List?weakExteriorQ] := {pt1, 
              inters[pt1, pt2]}; 
          clipLine[pt1_List?weakInteriorQ, pt2_List?exteriorQ] := {pt1, 
              inters[pt1, pt2]}; 
          clipLine[pt1_List?weakExteriorQ, 
              pt2_List?interiorQ] := {inters[pt1, pt2], pt2}; 
          clipLine[pt1_List?exteriorQ, 
              pt2_List?weakInteriorQ] := {inters[pt1, pt2], 
              pt2}; 
clipLine[pt1_List?weakInteriorQ, 
              pt2_List?weakInteriorQ] := {pt1, 
              pt2}; 
clipContour[pts_List] := 
            Flatten[Apply[clipLine, Partition[Append[pts, First[pts]], 2, 1], 
                  1], 1] //. {a___, b_List, b_, c___} :> {a, b, 
                  c}; 
beyondClipLine[pt1_List?interiorQ, 
              pt2_List?interiorQ] := Sequence @@ {}; 
          beyondClipLine[pt1_List?weakInteriorQ, 
              pt2_List?exteriorQ] := {inters[pt1, pt2], 
              pt2}; 
beyondClipLine[pt1_List?interiorQ, 
              pt2_List?weakExteriorQ] := {inters[pt1, pt2], pt2}; 
          beyondClipLine[pt1_List?exteriorQ, pt2_List?weakInteriorQ] := {pt1, 
              inters[pt1, pt2]}; 
          beyondClipLine[pt1_List?weakExteriorQ, pt2_List?interiorQ] := {pt1, 
              inters[pt1, pt2]}; 
          beyondClipLine[pt1_List?weakExteriorQ, 
              pt2_List?weakExteriorQ] := {pt1, 
              pt2}; 
beyondClipContour[pts_List] := 
            Flatten[Apply[beyondClipLine, 
                  Partition[Append[pts, First[pts]], 2, 1], 1], 
                1] //. {a___, b_List, b_, c___} :> {a, b, c}; 
          prepareQ = PrepareForSlices /. {opts} /. Options[LinearClip]; 
          boundaryLine[pts_List] := 
            Cases[Partition[clipContour[pts], 2, 
                1], {p__List?onBoundaryQ} :> 
                Line[{p}]]; 
Show[
            If[prepareQ, 
              ratObj /. Polygon[
                    pts_List?borderTouchingQ] :> 
                  Sequence @@ {Polygon[clipContour[pts]], 
                      Polygon[beyondClipContour[pts]], 
                      Sequence @@ boundaryLine[clipContour[pts]]}, 
              winnowed /. Polygon[
                    pts_List?borderTouchingQ] :> 
                  Polygon[clipContour[pts]]], 
            FilterOptions[Graphics3D, opts, 
              Sequence @@ Rest[obj]]]];

SliceSurface

SliceSurface[obj_SurfaceGraphics, other___] := 
        SliceSurface[Graphics3D[obj], other];

SliceSurface[obj_Graphics3D, linForm_, x_Symbol, y_Symbol, z_Symbol, 
            opts___?OptionQ] /; (PolynomialQ[linForm, {x, y, z}] && 
              Max[Map[Exponent[Expand[linForm], #] &, {x, y, z}]] === 1) := 
        Module[{sl, mshsty, lnsty, ratObj, vertList, extremaOfLinForm, 
            sliceList, 
            singleStep}, 
sl = Slices /. {opts} /. 
              Options[SliceSurface]; 
mshsty = 
            Flatten[{MeshStyle} /. {opts} /. 
                Options[SliceSurface]]; 
lnsty = 
            Flatten[{SliceEdgeStyle} /. {opts} /. 
                Options[SliceSurface]]; 
ratObj = 
            Rationalize[obj[[1]], 1/10^6]; 
vertList = 
            Cases[ratObj, 
              Polygon[pts_List] :> Sequence @@ pts, {1, 
                Infinity}]; 
extremaOfLinForm = {Min[#], 
                  Max[#]} &[
              linForm /. 
                Map[Thread[{x, y, z} -> #] &, 
                  vertList]]; 
If[
            Length[Union[extremaOfLinForm]] === 1, 
            Print["The object is flat in the direction of the slices"]; 
            Abort[]]; 
sliceList = 
            Which[IntegerQ[sl] && sl > 0, 
              Drop[Rest[
                  Range @@ (extremaOfLinForm~
                        Join~{Abs[
                              Subtract @@ extremaOfLinForm]/(sl + 
                                1)})], -1], 
              MatchQ[sl, {val__?((NumericQ[#] && 
                            FreeQ[N[#], Complex]) &)}], 
              sl, 
True, 
              Print["can't handle these slices"]; 
              Abort[]]; 
singleStep[gr_, 
              val_] := (Graphics3D[LinearClip[gr, 
                    linForm ≤ val, x, y, z, 
                    DisplayFunction -> Identity, 
                    PrepareForSlices -> 
                      True][[1]]]); 
Show[
            Graphics3D[
              Join[mshsty, 
                lnsty, Fold[singleStep, Graphics3D[ratObj], 
                    sliceList][[1]]], 
              FilterOptions[Graphics3D, opts, 
                Sequence @@ Rest[List @@ obj]]]]];

LevelPlot3D

LevelPlot3D[f_, {x_Symbol, xMin_, xMax_}, {y_Symbol, yMin_, yMax_}, 
          opts___] := 
        Module[{projQ, zlev, srf, 
            squished}, 
{projQ, 
              zlev} = {ShowProjectedLevelSets, ProjectAt} /. {opts} /. 
              Options[LevelPlot3D]; 
          srf = SliceSurface[
              Plot3D[f, {x, xMin, xMax}, {y, yMin, yMax}, 
                DisplayFunction -> Identity, 
                Evaluate[FilterOptions[Plot3D, opts]]], z, x, y, z, 
              DisplayFunction -> Identity, 
              Evaluate[Sequence @@ ({opts} /. Contours -> Slices)], 
              Evaluate[
                Sequence @@ 
                  Join[Options[LevelPlot3D]]]]; 
squished := 
            Cases[srf, 
              Line[pts_] :> 
                Line[Map[{{1, 0, 0}, {0, 1, 0}, {0, 0, 0}} . # + {0, 0, 
                          zlev} &, pts]], {1, 
                Infinity}]; 
If[
            projQ, 
Show[{srf, 
                Graphics3D[squished]}, 
FilterOptions[
                Graphics3D, opts], 
              DisplayFunction -> $DisplayFunction], 
Show[srf, FilterOptions[Graphics3D, opts], 
              DisplayFunction -> $DisplayFunction]]];

End[];

Protect["`*"];

EndPackage[];

Examples

LevelPlot3D

You can draw level sets of a function (it takes a while for a default mesh):

LevelPlot3D[x^2 - y^2, {x, -1, 1}, {y, -1, 1}, 
        BoxRatios -> {1, 1, 1}];

With the option Contours you can control the level sets to be drawn. If you set it to an integer, the level sets will be equally spaced between the supremum and the infimum of the function on the domain. The default is Contours->11. Contours can also be an explicit list of values.

LevelPlot3D[x^2 - y^2,{x, -1, 1},{y, -1, 1},BoxRatios -> {1, 1, 1},Contours -> {0, 1/6}];

With the option ShowProjectedLevelSets set to True the level lines are shown both on the surface and projected onto a horizontal plane:

LevelPlot3D[x^2 + y^2,{x, -1, 1},{y, -1, 1},ShowProjectedLevelSets -> True,BoxRatios -> {1, 1, 1}];

The horizontal plane is by default at z = 0. You can choose a different value of z with the option ProjectAt:

LevelPlot3D[x^2 - y^2,{x, -1, 1},{y, -1, 1},ShowProjectedLevelSets -> True,ProjectAt -> 1,BoxRatios -> {1, 1, 1}];

SliceSurface

The function SliceSurface is more general than LevelPlot3D, because you can give an arbitrary surface and slice it in arbitrary directions. The following example takes a while to compute:

ball = 
        ParametricPlot3D[{4  Cos[u] Cos[v], 3  Sin[u] Cos[v], 2  Sin[v]}, {u, 
            0, 2  Pi}, {v, -Pi/2, Pi/2}];
SliceSurface[ball, x + 2  y - z, x, y, z, Slices -> 8, 
        Ticks -> None];

LinearClip

The function LinearClip is mainly a component of SliceSurface, but it can be used independently. It draws the intersection of a surface with a half-plane.

Start with a paraboloid:

paraboloid = 
        ParametricPlot3D[
          r {Cos[t], Sin[t], r}, {r, 0, 1}, {t, 0, 2  Pi}];

The paraboloid cut by a plane:

Show[{paraboloid, 
          Graphics3D[{Polygon[{x, 
                    y, -x/4 - y/3 + 1/3} /. {{x -> -1, 
                      y -> -1}, {x -> 1, 
                      y -> -1}, {x -> 1, 
                      y -> 1}, {x -> -1, 
                      y -> 1}}]}]}];

The part of the paraboloid above the plane:

LinearClip[paraboloid,z ≥ -x/4 - y/3 + 1/3,x,y,z];

The part of the paraboloid below the plane:

LinearClip[paraboloid,z ≤  -x/4 - y/3 + 1/3,x,y,z];