# How can I improve the plotting speed of plot3D with region calculations?

I'm looking to achieve the desired result as shown in the attached image. There may be various methods to obtain such a plot, but my current code seems to be running quite slowly. Unfortunately, my expertise only allows me to write code of this nature. I'm wondering if anyone can help me optimize my code for better performance or suggest alternative approaches to create this plot more efficiently. Thank you

Clear["Global*"]
s0 = PauliMatrix[1];

s1 = PauliMatrix[2];
s2 = PauliMatrix[3];
s3 = PauliMatrix[4];
Ham =  kx s1 + ky s2 ;
DD = Eigenvalues[Ham];
VV = Eigenvectors[Ham];

region = RegionUnion[Disk[{0, 0}, 0.1],
RegionDifference[Disk[{0, 0}, 0.3], Disk[{0, 0}, 0.2]],
RegionDifference[Disk[{0, 0}, 0.5], Disk[{0, 0}, 0.4]],
RegionDifference[Disk[{0, 0}, 0.7], Disk[{0, 0}, 0.6]],
RegionDifference[Disk[{0, 0}, 0.9], Disk[{0, 0}, 0.8]],
RegionDifference[Disk[{0, 0}, 1.1], Disk[{0, 0}, 1]],
RegionDifference[Disk[{0, 0}, 1.3], Disk[{0, 0}, 1.2]],
RegionDifference[Disk[{0, 0}, 1.5], Disk[{0, 0}, 1.4]],
RegionDifference[Disk[{0, 0}, 1.7], Disk[{0, 0}, 1.6]],
RegionDifference[Disk[{0, 0}, 1.9], Disk[{0, 0}, 1.8]]];

region2 = RegionDifference[Disk[{0, 0}, 1], region];

P1 = Plot3D[DD[[1]], {kx, ky} \[Element] region,
PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1},
ClippingStyle -> None,
ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &),
Mesh -> None, Axes -> True, Boxed -> False, MaxRecursion -> 2,
PlotPoints -> 20];
P2 = Plot3D[DD[[2]], {kx, ky} \[Element] region2,
PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1},
ClippingStyle -> None,
ColorFunction -> (ColorData["SunsetColors"][#3] &), Mesh -> None,
Axes -> True, Boxed -> False, MaxRecursion -> 2, PlotPoints -> 20];
P3 = Graphics3D[{Red, Sphere[{0, 0, 0}, 0.1]}];
Show[P1, P2, BoxRatios -> {1, 1, 1.2}]

• Any insights into what might be causing the slow performance would be greatly appreciated. Oct 2 at 9:16

You can get the desired result much faster using the options MeshFunctions, Mesh and MeshShading as follows:

Plot3D[DD[[1]], {kx, -2, 2}, {ky, -2, 2}, PlotRange -> {-1.5, 1.5},
BoxRatios -> {1, 1, 1.1}, ClippingStyle -> None,
ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &),
Axes -> True, Boxed -> False, MaxRecursion -> 2, PlotPoints -> 20,
MeshFunctions -> {#3 &}, Mesh -> {Range[-2, 2, .1]},
MeshShading -> {Automatic, None}]


Do the same for DD[[1]] and DD[[2]] with association color schemes:

Show[
Plot3D[DD[[#]], {kx, -2, 2}, {ky, -2, 2},
PlotRange -> {-1.5, 1.5},
BoxRatios -> {1, 1, 1.1},
ClippingStyle -> None,
ColorFunction -> ({x, y, z} |-> ColorData[#2][z]),
Axes -> True,
Boxed -> False,
MaxRecursion -> 2,
PlotPoints -> 20,
MeshFunctions -> {#3 &},
Mesh -> {Range[-2, 2, .1]},
MeshShading -> {Automatic, None}] & @@@
{{1, {"SunsetColors", "Reverse"}},
{2, "SunsetColors"}},
PlotRange -> All]


Aside:

If you have to work with regions, you might consider Annulus to define your region:

region2 = Apply[RegionUnion] @
Prepend[Disk[{0, 0}, .1]] @
Map[Annulus[{0, 0}, {#, # + .1}] &, Range[.2, 1.9, .2]];

Region @ region2


• thank you very much! Oct 2 at 11:16
• @ZhongfuLi, my pleasure. Welcome to mmase.
– kglr
Oct 2 at 11:19

We can BoundaryDiscretizeRegion the regions region and region2 to speed up the 3D plot.

region = BoundaryDiscretizeRegion[region];
region2 = BoundaryDiscretizeRegion[region2];


That is

Clear["Global*"]
s0 = PauliMatrix[1];

s1 = PauliMatrix[2];
s2 = PauliMatrix[3];
s3 = PauliMatrix[4];
Ham = kx s1 + ky s2;
DD = Eigenvalues[Ham];
VV = Eigenvectors[Ham];

region =
RegionUnion[Disk[{0, 0}, 0.1],
RegionDifference[Disk[{0, 0}, 0.3], Disk[{0, 0}, 0.2]],
RegionDifference[Disk[{0, 0}, 0.5], Disk[{0, 0}, 0.4]],
RegionDifference[Disk[{0, 0}, 0.7], Disk[{0, 0}, 0.6]],
RegionDifference[Disk[{0, 0}, 0.9], Disk[{0, 0}, 0.8]],
RegionDifference[Disk[{0, 0}, 1.1], Disk[{0, 0}, 1]],
RegionDifference[Disk[{0, 0}, 1.3], Disk[{0, 0}, 1.2]],
RegionDifference[Disk[{0, 0}, 1.5], Disk[{0, 0}, 1.4]],
RegionDifference[Disk[{0, 0}, 1.7], Disk[{0, 0}, 1.6]],
RegionDifference[Disk[{0, 0}, 1.9], Disk[{0, 0}, 1.8]]];

region2 = RegionDifference[Disk[{0, 0}, 1], region];

region = BoundaryDiscretizeRegion[region];
region2 = BoundaryDiscretizeRegion[region2];

P1 = Plot3D[DD[[1]], {kx, ky} ∈ region,
PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1},
ClippingStyle -> None,
ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &),
Mesh -> None, Axes -> True, Boxed -> False, MaxRecursion -> 2,
PlotPoints -> 20];
P2 = Plot3D[DD[[2]], {kx, ky} ∈ region2,
PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1},
ClippingStyle -> None,
ColorFunction -> (ColorData["SunsetColors"][#3] &), Mesh -> None,
Axes -> True, Boxed -> False, MaxRecursion -> 2,
PlotPoints -> 20];
P3 = Graphics3D[{Red, Sphere[{0, 0, 0}, 0.1]}];
Show[P1, P2, BoxRatios -> {1, 1, 1.2}]

• thanks! it works! Oct 3 at 9:37