# Plot more points only in specific region of plot

I have the code shown at the bottom of this question which generates a 3D plot. There are locations where the two surfaces you see form cones and touch at infinitesimally small points. These points always occur at z=1 and I have to use many PlotPoints to resolve this touching.

I will be generating thousands of these plots (to form a movie) so need to plot each one as quickly as possible. I would imagine a lot of time would be saved if I used a large number of points only in the regions of interest. Is there a way I can use a large PlotPoints density in a rectangular box that spans the x and y domain but has a small height and is vertically centred around z=1 where I know these points occur (but I don't know their x or y)?

(*The lattice vectors*)
a1 = {Sqrt[3], 0};
a2 = {Sqrt[3]/2, 3/2};

(*Omega/w_0*)
Omega = 0.01;

wp[qx_, qy_, r_] := Module[{},
q = {qx, qy};

(*Nearest neighbour vectors*)
{d1, d2, d3} = # - r & /@ {{0, -1}, {Sqrt[3]/2, 1/2}, {-Sqrt[3]/2, 1/2}};

(*The c_j's*)
{theta, phi} = {0, 0};
{c1, c2, c3} = (1 - 3 Sin[theta]^2 Cos[ArcTan[#[[1]], #[[2]]] - phi + Pi/2]^2)/
Norm[#]^3 & /@ {d1, d2, d3};

modfq =
Sqrt[c1^2 + c2^2 + c3^2 + 2 c1 c2 Cos[q.(d1 - d2)] +
2 c1 c3 Cos[q.(d1 - d3)] +  2 c2 c3 Cos[q.(d2 - d3)]];
{Sqrt[1 + 2 Omega modfq], Sqrt[1 - 2 Omega modfq]}
]

r = {0, 0};

Timing[
With[{plotopts = {Mesh -> None, PlotStyle -> Opacity[0.7],
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}, Automatic},
PlotPoints -> 100, ViewPoint -> {1.43, -2.84, 1.13},
ViewVertical -> {0., 0., 1.}}},
plot1 =
Plot3D[wp[qx, qy, r][[1]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts];
plot2 =
Plot3D[wp[qx, qy, r][[2]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts];
]
]
plot =
Show[plot1, plot2, PlotRange -> {0.96, 1.04},
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}}, LabelStyle -> Medium,
BoxRatios -> {2, 2, 3}, BoxStyle -> Opacity[0.4]]


Note: This question could be a duplicate of here, but from what I gathered it was about feeding it a list of points rather than a region? I could be wrong.

• Related: (10414). Commented Jul 30, 2014 at 19:31

Unfortunately the recursion algorithm fails to increase points near the Dirac points (actually, with a narrow band gap). It seems that it is because $\partial f(x,y)/\partial x,\ \partial f(x,y)/\partial y \ll 1$. However you can scale your function, scale back with a post-processing and obtain a nice plot!

scale = 1000;
postProcess[g_] :=
g /. GraphicsComplex[pts_, opts___] :>
GraphicsComplex[(pts\[Transpose]/{1, 1, scale})\[Transpose],
opts] /. (VertexNormals ->
n_) :> (VertexNormals -> (n\[Transpose] {1, 1,
scale})\[Transpose]);
Timing[With[{plotopts = {Mesh -> None, PlotStyle -> Opacity[0.7],
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}, Automatic},
PlotPoints -> 10, MaxRecursion -> 4,
ViewPoint -> {1.43, -2.84, 1.13},
ViewVertical -> {0., 0., 1.}}},
plot1 =
Plot3D[(wp[qx, qy, r][[1]]) scale, {qx, -Pi, Pi}, {qy, -Pi, Pi},
plotopts] // postProcess;
plot2 =
Plot3D[(wp[qx, qy, r][[2]]) scale, {qx, -Pi, Pi}, {qy, -Pi, Pi},
plotopts] // postProcess;]]
plot = Show[plot1, plot2, PlotRange -> {0.96, 1.04},
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}}, LabelStyle -> Medium,
BoxRatios -> {2, 2, 3}, BoxStyle -> Opacity[0.4]]


You can see a fine mesh near the Dirac points:

Now you can obtain a considerable speedup by tuning PlotPloits and MaxRecursion.

It looks like the plotting may be slow because you're trying to resolve a small feature and upping the number of points you use in the entire plot.

If you instead try upping the number of points only in the region between 0.99 and 1.01, where you know the sharp features are, you cut down on the amount of time spent searching:


Timing[With[{plotopts = {Mesh -> None, PlotStyle -> Opacity[0.7],
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}, Automatic},
ViewPoint -> {1.43, -2.84, 1.13},
ViewVertical -> {0., 0., 1.}}},
plot1 =
Plot3D[wp[qx, qy, r][[1]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts,
PlotRange -> {{-π, π}, {-π, π}, {1.00, 1.01}},
PlotPoints -> 40, ClippingStyle -> None, BoundaryStyle -> None];
plot2 =
Plot3D[wp[qx, qy, r][[2]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts,
PlotRange -> {{-π, π}, {-π, π}, {0.99, 1.00}},
PlotPoints -> 40, ClippingStyle -> None, BoundaryStyle -> None];
plot3 =
Plot3D[wp[qx, qy, r][[1]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts,
PlotRange -> {{-π, π}, {-π, π}, {1.01, 1.04}},
PlotPoints -> 20, ClippingStyle -> None, BoundaryStyle -> None];
plot4 =
Plot3D[wp[qx, qy, r][[2]], {qx, -Pi, Pi}, {qy, -Pi, Pi}, plotopts,
PlotRange -> {{-π, π}, {-π, π}, {0.96, 0.99}},
PlotPoints -> 20, ClippingStyle -> None, BoundaryStyle -> None];
plot = Show[{plot1, plot2, plot3, plot4}, PlotRange -> {0.96, 1.04},
Ticks -> {{-Pi, 0, Pi}, {-Pi, 0, Pi}}, LabelStyle -> Medium,
BoxRatios -> {2, 2, 3}, BoxStyle -> Opacity[0.4]]

]
]



When I run this it looks the same as the plot you've posted but only takes about 4 seconds to plot.

One can use FEM's ToElementMesh, MeshRefinementFunction, and ElementMeshPlot3D to refine the plot. ElementMeshPlot3D does not work like Plot3D. In particular VertexNormals are not computed, so I borrowed the utility function addnormals[] from How to improve this plot? Also there is no PlotStyle option, so some other method needs to be used; it does have a ColorFunction option.

Needs["NDSolveFEM"];
emesh = ToElementMesh[FullRegion[2], {{-Pi, Pi}, {-Pi, Pi}},
MeshRefinementFunction ->
Function[{vertices, area},
Abs[Mean[wp[##, r][[1]] & @@@ vertices] - wp[##, r][[1]] & @@
Mean[vertices]] > 0.0001]
];
Show[emesh["Wireframe"], Frame -> True]


Needs["NDSolveFEM"];
ClearAll[plotwp];

plotwp[wpr_, opts : OptionsPattern[ElementMeshPlot3D]] :=
Module[{mins, wp1C, wp2C, vals1, eIF1, vals2, eIF2, emesh},
wp1C = Compile[{qx, qy}, Evaluate[wpr[qx, qy][[1]] /. Sqrt[s_] :> Sqrt@Abs[s]]];
wp2C = Compile[{qx, qy}, Evaluate[wpr[qx, qy][[2]] /. Sqrt[s_] :> Sqrt@Abs[s]]];
emesh = ToElementMesh[FullRegion[2], {{-Pi, Pi}, {-Pi, Pi}},
MeshRefinementFunction ->
Function[{vertices, area},
Abs[Mean[wp1C @@@ vertices] - wp1C @@ Mean[vertices]] > 0.0001]
];
vals1 = wp1C @@@ emesh["Coordinates"];
eIF1 = ElementMeshInterpolation[{emesh}, vals1];
vals2 = wp2C @@@ emesh["Coordinates"];
eIF2 = ElementMeshInterpolation[{emesh}, vals2];
Show[
addnormals[eIF1][ElementMeshPlot3D[eIF1, BoxRatios -> {1, 1, 0.6}, opts]],
addnormals[eIF2][ElementMeshPlot3D[eIF2, BoxRatios -> {1, 1, 0.6}, opts]],
PlotRange -> All, Options[plot]]
]

pl = plotwp[wp[##, r] &, ColorFunction -> None]; // AbsoluteTiming // First
pl = DeleteCases[pl, HoldPattern[VertexColors -> _], Infinity] /.
GraphicsComplex[p_, g_, r___] :>
GraphicsComplex[
p, {First@ChartinggetPlotStyles[
"DefaultPlotStyle" /. (Method /.
ChartingResolvePlotTheme[\$PlotTheme, Plot3D])][1, Opacity[0.7]],
g}, r]
(* 0.529355  <-- timing  *)


Utility:

(*add VertexNormals to plot assumed to be the plot of the function u*)
addnormals[u_][plot_Graphics3D /; ! FreeQ[plot, GraphicsComplex]] :=
plot /. GraphicsComplex[p_, rest___] :>
GraphicsComplex[p, rest,
VertexNormals ->
With[{dx = Derivative[1, 0][u], dy = Derivative[0, 1][u],
xy = p[[All, {1, 2}]]},
Transpose@{-dx @@@ xy, -dy @@@ xy, ConstantArray[1., Length@p]}]
]