Plot more points only in specific region of plot

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["NDSolve`FEM`"];
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["NDSolve`FEM`"];
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@Charting`getPlotStyles[
       "DefaultPlotStyle" /. (Method /. 
          Charting`ResolvePlotTheme[$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]}]
    ]