$C^0$ estimate for solutions of elliptic PDE with Neumann BC

I am interested in a reference for (or counterexample to) a particular $C^0$ estimate for solutions of the Laplace equation with Neumann boundary conditions. More precisely, let $(M,g)$ be a $C^\infty$ compact Riemannian manifold with non-empty smooth boundary $\partial M$. Let $\Delta_g=-div_g\nabla$ denote the Laplace operator and $\nu$ the outward pointing unit normal on $\partial M$.

Consider (strong) solutions $u$ to the following system \begin{cases} \Delta_gu=f&\text{ in }M\\ \partial_\nu u=h&\text{ on }\partial M \end{cases} where $f\in C^\infty(M)$ and $h\in C^\infty(\partial M)$ satisfy $\int_M fdvol_g=\int_{\partial M}hd\sigma_g$ (the smoothness condition is for simplicity but the later condition is necessary for the problem to be well-posed). Here $dvol_g$ denotes the Riemannian volume form and $d\sigma_g$ denotes its pullback under the inclusion $\partial M\to M$. To specify $u$ uniquely, let us also suppose that it has zero average value i.e. $\int_Mudvol_g=0$. My question is this: does there exist $C>0$, independent of $f$ and $h$, such that $\sup_M|u|\leq C(\sup_M|f|+\sup_{\partial M}|h|)$?

Most classical sources (Evans, Gilbarg+Trudinger, and Aubin) either only consider Dirichlet conditions (in which case the analogous estimate is essentially just the maximum principle) or only provide Schauder estimates for the Neumann problem -- this is not quite what I want.

Does anyone know of a resource for this kind of estimate or of a counterexample?

• It might be easier to study a first order system for a vector field: Set $w=\nabla u$. Then the system becomes $-\text{div}(w)=f$ in $M$ and $\nu\cdot w=h$ on $\partial M$. Since $u$ must have a zero somewhere (assuming connectedness), $\sup_M|u|\leq\text{diam}(M)\sup_M|w|$, so it suffices to estimate $w$ uniformly. – Joonas Ilmavirta Sep 23 '14 at 19:42
• If f and h are bounded, they are in every $L^p$. You can then get what you want by using $L^p$ estimates and Sobolev embedding. – Michael Renardy Sep 24 '14 at 0:07
• they are both smooth and the domain is compact, so they are certainly bounded. – Matthias Ludewig Sep 27 '14 at 18:40

Yes, I think the estimate you propose is true. As a simple case take $f = 0$ and $M = B_1 \subset \mathbb{R}^n$. By adding a constant assume that $\inf_{B_1} u = 0$ and $\sup_{B_1}u = K,$ and note that these extrema are on the boundary. Note that $u(0)$ is closer to $0$ or $K$, say $u(0) > K/2$ without loss of generality. Then by the Harnack inequality we have $u > cK$ in $B_{1/2}$, so the function $$c(n)K(|x|^{2-n} - 1)$$ is a lower barrier in the annulus $B_1-B_{1/2}$. In particular, at the point where $u$ takes its minimum we have $|h| > c(n)K$ which gives the estimate.
In a more general domain we can argue similarly; take an interior ball touching the point on the boundary where $u = 0$ and argue that $u$ is no larger than $C\sup|h|$ at its center, and similarly take a ball touching the point where $u = K$ and argue that $u$ is larger than $K - C\sup|h|$ at its center. If $K$ is much larger than $\sup|h|$ we contradict the Harnack inequality (with a constant depending on the geometry of $M$).
Finally, if $f$ is nonzero similar arguments should work since we still have the Harnack inequality (now depending on $\sup|f|$) and for a barrier we can take $|x|^{-\alpha}$ for $\alpha$ large depending on $\sup|f|$.