  Matrix-free finite element implementations of massively parallel geometric
multigrid save memory and are often significantly faster than implementations
using classical sparse matrix techniques. They are especially well suited for
hierarchical hybrid grids on polyhedral domains. In the case of constant
coefficients all fine grid node stencils in the interior of a coarse macro
element are equal. However, for non-polyhedral domains the situation changes.
Then even for the Laplace operator, the non-linear element mapping leads to
fine grid stencils that can vary from grid point to grid point. This
observation motivates a new two-scale approach that exploits a piecewise
polynomial approximation of the fine grid operator with respect to the coarse
mesh size. The low-cost evaluation of these surrogate polynomials results in an
efficient stencil assembly on-the-fly for non-polyhedral domains that can be
significantly more efficient than matrix-free techniques that are based on an
element-wise assembly. The performance analysis and additional hardware-aware
code optimizations are based on the Execution-Cache-Memory model. Several
aspects such as two-scale a priori error bounds and double discretization
techniques are presented. Weak and strong scaling results illustrate the
benefits of the new technique when used within large scale PDE solvers.
