Given the Poisson's equation
\begin{equation} -\nabla^2 u = f \quad \mathrm{in} \ \Omega \end{equation} \begin{equation} u=0 \quad \mathrm{on} \ \Gamma_D, \quad \frac{\partial u}{\partial n} = 0 \quad \mathrm{on} \ \Gamma_N \end{equation}
Is it possible to prove that if $f(x) \geq 0 \ \forall x \in \Omega$ then $u(x) \geq 0 \ \forall x \in \Omega$?
As long as $\Gamma_D$ is not empty, then yes. To prove this, consider the perturbed problem \begin{equation} -\nabla^2 u_\varepsilon = f + \varepsilon \quad \mathrm{in} \ \Omega \end{equation} \begin{equation} u_\varepsilon=0 \quad \mathrm{on} \ \Gamma_D, \quad \frac{\partial u_\varepsilon}{\partial n} = \varepsilon \quad \mathrm{on} \ \Gamma_N \end{equation} where $\varepsilon > 0$. Here $\nabla^2$ is the Laplacian. Since $f\geq 0$, the Laplacian $\nabla^2 u_\varepsilon$ is strictly negative everywhere, and so $u_\varepsilon$ cannot have an interior minimum (since at any interior minimum, the Hessian would be nonnegative, and so the Laplacian---the trace of the Hessian--would also be nonnegative). Likewise, since the normal derivative is strictly positive, the minimum cannot occur on $\Gamma_N$. Hence, the minimum of $u_\varepsilon$ occurs on $\Gamma_D$, where $u_\varepsilon=0$. So $u_\varepsilon\geq0$. Then send $\varepsilon\to 0$, and use that $u_\varepsilon \to u$ uniformly.