  We present a new meshless method for scalar diffusion equations which is
motivated by their compatible discretizations on primal-dual grids. Unlike the
latter though, our approach is truly meshless because it only requires the
graph of nearby neighbor connectivity of the discretization points $\bm{x}_i$.
This graph defines a local primal-dual grid complex with a \emph{virtual} dual
grid, in the sense that specification of the dual metric attributes is implicit
in the method's construction. Our method combines a topological gradient
operator on the local primal grid with a Generalized Moving Least Squares
approximation of the divergence on the local dual grid. We show that the
resulting approximation of the div-grad operator maintains polynomial
reproduction to arbitrary orders and yields a meshless method, which attains
$O(h^{m})$ convergence in both $L^2$ and $H^1$ norms, similar to mixed finite
element methods. We demonstrate this convergence on curvilinear domains using
manufactured solutions. Application of the new method to problems with
discontinuous coefficients reveals solutions that are qualitatively similar to
those of compatible mesh-based discretizations.
