I know that for harmonic function $0\leq u\in C^2(D)$, there holds $$ \left| \nabla u\left( x_0 \right) \right|\le \frac{C}{r^{d+1}}\int_{B_r\left( x_0 \right)}{ u\left( x \right)}dx,\quad B_r\left( x_0 \right) \subset D . $$
So I want to know that whether this estimation also holds for the solution of general elliptic equations such as $$ Lu=D_j(a_{ij}D_iu)=0,$$ if so, what assumptions maybe need for $(a_{ij})$?
I know the estimation holds when $u$ is harmonic since $u$ satisfies mean value property, but when $u$ is a solution of $Lu=0$, what formula can be found to replace the mean-value-property.