Sum of discrete gradient of a-harmonic functions

Question

Let's say we have a solution to a constant coefficient divergence elliptic equation of the form:
$\nabla.a(\nabla \bar u)=0$ on U with $\bar u=u_e$ on $\partial U$.
(For the sake of this question let's assume U is a smoothened out square with the smooth curves around the corners contained in a unit square around the corner).
Now this $u_e$ is constructed as an extension of u on $\mathbb{Z}_d$ in such a way that the $\|u_e\|_{L_p}=\|u\|_{l_p}$. And by elliptic regularity theory I have bounds on $L_p$ norm of $\nabla \bar u$ using the $L_p$ norm of $\nabla u_e$ or just the discrete sum. I am interested in bounding the sum of the discrete gradient of $\bar u$. For this I basically used the the mean value theorem and then triangle inequality to get a bound on the sum of $\bar u$ on $\mathbb{Z}_d \cap U$ in terms of the $L_p$ norm of $\nabla u$ and $u$:
Fix a ball of size r around z in $\mathbb{Z}_d$ so that r<1. Now $\int_{z+B_r} |f(z)-f(x)|= \int_{z+B_r} |\unicode{x2A0F}_{z+B_r}f(x)-f(x)| \leq \int_{z+B_r} |\nabla f|$. After that I used triangle inequality and then summed over the z in the domain to get $\sum f(z)\leq \int_U |\nabla f|+\int_U|f(x)|$ and to estimate the sum of discrete gradient I just used this inequality and the triangle inequality.
For $L_p$ norm we can get a similar bound using Minkowski's inequality too. But I was wondering if it is possible to get a similar bound just using the gradient(also keep in mind that I will be using this bound for an elliptic regularity iteration so a boundary dependent constant such as from Poincare inequality on the whole domain wouldn't work).
So I tried to estimate the sum discrete gradient itself and not that of $\bar u$.
$\sum_{x\in U} \sum_{x\sim y} |f(x)-f(y)|$. For a fixed x$\in \mathbb{Z}_d \cap U$ I take a ball of radius 1 now which will have all it's neighbors on the boundary, then using the mean value theorem again I get $|f(x)-f(y)|=|\unicode{x2A0F}_{x+B_1}f(z)-f(y)| \leq \unicode{x2A0F}_{x+B_1}|f(z)-f(y)|$. If there was a way to bound this using just the gradient that would work or if there is another way to get the required bound that would be great too. Any kind of help will be appreciated.
Thanks!