Firstly, turning to the $\mathbb{R}^1$ case, we note that
$$
   -u_{xx}+u^2u_x=-\partial_x\left(u_x-\frac{1}{3}u^3\right)=f(x)
$$
and so we get a Chini equation (see http://mathoverflow.net/questions/104006/inhomogeneous-bernoulli-equation):
$$
   -u_x+\frac{1}{3}u^3=\int_{x_0}^xf(x')dx'+C
$$
and a general solution to this equation is not known. Turning to $\mathbb{R}^3$, you can see that your equation can be written down in the form
$$ 
   \nabla\cdot {\bf v}=-f
$$
being
$$
   {\bf v}=\left(u_x-\frac{1}{3}u^3,u_y-\frac{1}{3}u^3,u_z-\frac{1}{3}u^3\right)
$$
and you have to satisfy the Gauss identity
$$
   \int_{\partial\Omega}{\bf v}\cdot d{\bf S}=-\int_{\Omega}dV f.
$$
This means that one has to find three functions $v_1,\ v_2,\ v_3$ that satisfy three different Chini equations. Again, a general solution misses in this case. The conclusion to draw is that you need to specialize the problem to see if this becomes manageable in some case or resort to some perturbation technique.