$\Delta u = \operatorname{div}f \ \ \mbox{in} \ \ B_1, f \in L^2 \Rightarrow \nabla u \in L^2$ I'm looking for results like, If $f \in L^p$ 
and 
$$ 
\begin{array}{rclcl}
\Delta u & = & \operatorname{div}f & \mbox{in} &  B_1\\
u&=&0& \mbox{on}& \partial B_1
\end{array}
$$
then 
$$ \int_{B_1} \!\left| \nabla u \right|^2 \le \int_{B_1}\!\left|\,f\right|^2 .$$
This is, $f \in L^2 \Rightarrow \nabla u \in L^2$. If  you know more generally, $f \in L^q \Rightarrow \nabla u \in L^q$. Is great if you can write the solution. If no, reference is good, preferable for non general cases of second order elliptic equations. But if you know only references for general case, I'm very grateful too.
 A: Define $F:H_0^1(B_1)\to\mathbb{R}$ by $$F(u)=\frac{1}{2}\int_{B_1}|\nabla u|^2-\int_{B_1}f\cdot\nabla u$$
I - $F$ is weakly lower semi continuous, i.e. if $u_n\to u$ weakly in $H_0^1(B_1)$, then $$F(u)\leq\liminf F(u_n)$$ 
To prove this, just note that the norm function is weakly lower semicontinous (you can find this result in any good book of functional analysis). On the other hand, the expression $\int_{B_1} f\cdot\nabla v$ for $v\in H_0^1$ is a linear functional on $H_0^1$, then, by definition of weak convergenge, we have that $$\int_{B_1} f\cdot\nabla u_n\to \int_{B_1} f\cdot\nabla u$$
The last convergence implies I.
II - $F$ is coercive.
Note that $-\int_{B_1} f\cdot\nabla v\geq -\|f\|_2\|\nabla u\|_2$, which implies coerciveness.
With I and II, we know that there exist $u\in H_0^1(\Omega)$, which minimizes $F$, then $\langle F'(u),v\rangle=0$ for all $v\in H_0^1$, or equivalently $$\tag{1}\int_{B_1}\nabla u\nabla v=\int_{B_1}f\cdot\nabla v,\ \forall\ v\in H_0^1$$
By taking $v=u$ in the last equality, we conclude by using Holder inequality that $$\int_{B_1}|\nabla u|^2\leq\int_{B_1}|f|^2$$
For the result $f\in L^q$ implies $\nabla u\in L^q$, I suggest you to take a look in Calderón Zigmund theory as suggested in another post.
Remark: Note that $F$ is strictly convex, so the solution is unique.
Remark 2: Note that if $f\in L^2(B_1)^N$, then $\operatorname{div}f\in H^{-1}(B_1)^N$, hence, $\langle \operatorname{div}f, v\rangle =-\int_{B_1}f\cdot\nabla v$ for all $v\in H_0^1(B_1)$, whicih implies that $u$ satisfies (in the weak sense) $$\Delta u=\operatorname{div}f$$
