Does weakly harmonic function satisfy mean value property or generalized mean value property?

I am taking a higher level PDE class right now and just reached Sobolev space. I got the idea that generalized derivative (or weak derivative) can be used to characterize not so smooth functions. Knowing that a $C^2$ harmonic function satisfies the mean value property, my question is: if $u\in W^{2,p}$, and $\Delta u = 0$ in the sense of test functions, does this $u$ satisfies mean value property as well?

If not, is there some generalized mean value property for Sobolev functions? for example, mean value property on a square or polygon or something.

This example might not illustrate your doubts that well. According to Weyl's lemma, any weakly harmonic function $u \in L^1_{loc}$ is in fact smooth (up to redefining on a set of measure zero), therefore has the mean value property. Note that the assumption $u \in L^1_{loc}$ is essentially weaker than $u \in W^{2,p}$.