  We study the mean-value harmonic functions on open subsets of $\mathbb{R}^n$
equipped with weighted Lebesgue measures and norm induced metrics. Our main
result is a necessary condition saying that all such functions solve a certain
homogeneous system of elliptic PDEs. Moreover, a converse result is established
in case of analytic weights. Assuming Sobolev regularity of weight $w \in
W^{l,\infty}$ we show that strongly harmonic functions are as well in
$W^{l,\infty}$ and that they are analytic, whenever the weight is analytic.
  The analysis is illustrated by finding all mean-value harmonic functions in
$\mathbb{R}^2$ for the $l^p$-distance $1 \leq p \leq \infty$. The essential
outcome is a certain discontinuity with respect to $p$, i.e. that for all $p
\ne 2$ there are only finitely many linearly independent mean-value harmonic
functions, while for $p=2$ there are infinitely many of them. We conclude with
a remarkable observation that strongly harmonic functions in $\mathbb{R}^n$
possess the mean value property with respect to infinitely many weight
functions obtained from a given weight.
