# Why is my parallel version of DensityPlot slow?

I am exploiting the possibility to parallelize DensityPlot.

This is my naive implementation

Clear[parallelDensityPlotData];
parallelDensityPlotData[func_, xrange_, yrange_, partition_,
opts : OptionsPattern[]] := Module[{x, y},
extractPlot3Ddata[graphics_] := graphics[[1, 1]];
newXrangeList =
Partition[Subdivide[Sequence @@ xrange, partition[]], 2, 1];
newYrangeList =
Partition[Subdivide[Sequence @@ yrange, partition[]], 2, 1];
newRangeList = Tuples[{newXrangeList, newYrangeList}];
rangeToParaRange = Flatten /@ {{x, #[]}, {y, #[]}} &;
paraRangeList = rangeToParaRange /@ newRangeList;
graphicsList =
ParallelTable[
Plot3D[func[x, y], Evaluate[Sequence @@ i], Mesh -> None,
PlotRange -> All,
Evaluate@FilterRules[{opts}, Options[Plot3D]]], {i,
paraRangeList}];
ListDensityPlot[
Partition[Flatten[extractPlot3Ddata /@ graphicsList], 3],
Evaluate@FilterRules[{opts}, Options[ListDensityPlot]]]]


It is a simple implementation that partition the whole rectangle into smaller rectangle and Plot3D on each rectangle, then extract the each data, do a final ListDensityPlot.

However, I found my implementation just can't compete with serial running DensityPlot.

For example, take function from this post

f[x_, y_] := (x^2 + y^2) Exp[-x^2 - y^2] Sin[10 Sqrt[x^2 + y^2] + 10 ArcTan[x, y]]^4;


then

DensityPlot[f[x, y], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 24,
PlotRange -> All, MaxRecursion -> 3] // AbsoluteTiming While

parallelDensityPlotData[f, {-3, 3}, {-3, 3}, {3, 3}, PlotPoints -> 8,
MaxRecursion -> 3, PlotRange -> All] // AbsoluteTiming Though, I only have two cores, but this is kind of slow.

My guess is that Plot3D takes too much time calculating other things?

So I am wondering what is the core adaptive part function inside Plot3D so that we can directly use? I've tried to use PrintDefinitions, but got lost in the output.