Let $\mathrm{div}\,(A\,\mathrm{grad}\,u) + b u = f$. Is $(A\,\mathrm{grad}\,u)$ weakly differentiable? Let us consider the basic linear elliptic PDE
$$
\mathrm{div} (A\,\mathrm{grad}\,u) + bu = f,
$$
with $f\in L^p,$ $A,b$ uniformly bounded. Do we have, for a weak solution $u\in W^{1,p}(\Omega')$, 
$$
(A\,\mathrm{grad}\,u)\in (W^{1,p})^d(\Omega),\ \mbox{ in }\Omega\subset\subset\Omega'\subset\mathbb{R}^d?
$$
or at least do we have
$$
\mathbf{n}\cdot(A\,\mathrm{grad}\,u)\in L^{p}(\partial\Omega)?
$$
This seems to me true at least if $A$ is piecewise Lipschitz uniformly because across an interface where $A$ has a jump, the transmission condition $A\,\mathrm{grad}\,u$ still has matched traces.
My final goal is to know under what least conditions, we have Green's representation formula in the sense of Lebesgue integrals:
$$
u(\mathbf{y})=\int_{\partial\Omega}\frac{\partial G(\cdot,\mathbf{y})}
{\partial \mathbf{n}_A}u-
\frac{\partial u}{\partial \mathbf{n}_A}G(\cdot,\mathbf{y})~dS(\mathbf{x})
+\int_{\Omega}fG(\cdot,\mathbf{y})~d\mathbf{x},
$$
for $\mathbf{y}\in\Omega$ and $\frac{\partial v}{\partial\mathbf{n}_A}=\mathbf{n}\cdot(A\,\mathrm{grad}\,v).$
Edit: I'm still waiting a counter-example to the above guess when $A$ is bounded only.
 A: As far as I know, there is no such results if you only assume $A\in L^\infty(\Omega)$.
But if you assume that the domain $\Omega$ is partitioned into two subdomains $\Omega_1\cup\Omega_2$ by a smooth interface $\Gamma$ inside $\Omega$, and in each subdomain you assume that $\,A\in W^{1,\infty}(\Omega_1)\cap W^{1,\infty}(\Omega_2)$ ($\,A$ may be discontinuous across the interface), then the solution of
$$
\left\{\begin{array}{ll}
\nabla\cdot(A\nabla u)=f &\mbox{in}~~\Omega,\\
u=0 &\mbox{on}~~\Omega,
\end{array}\right.
$$
with $f\in L^p(\Omega)$ satisfies 
$$
\|u\|_{ W^{2,p}(\Omega_1)}+\|u\|_{W^{2,p}(\Omega_2)}\leq C\|f\|_{L^p(\Omega)} 
$$ 
and the trace of ${\bf n}\cdot A\nabla u$ on the two sides of the interface coincide, with ${\bf n}\cdot A\nabla u\in L^p(\partial\Omega)$.
If $A\in C^{1+\alpha}(\overline\Omega_1)\cap C^{1+\alpha}(\overline\Omega_2)$, then
$$
\|u\|_{ C^{2+\alpha}(\overline\Omega_1)}+\|u\|_{C^{2+\alpha}(\overline\Omega_2)}\leq C(\|f\|_{C^\alpha(\overline\Omega_1)}+\|f\|_{C^\alpha(\overline\Omega_2)}) .
$$ 
These classical results probably can be found in the following papers
"Gradient estimates for the perfect conductivity problem" by Ellen Shiting Bao, Yanyan Li and Biao Yin
"Gradient Estimates for the Perfect and Insulated Conductivity Problems with Multiple Inclusions" by Ellen Shiting Bao, Yanyan Li and Biao Yin
A: Suppose $p<n$. Remark that just from the equation, 
$$
\partial_i (A_{ij} \partial_j u) = A_{ij} \partial_{ij} u + \partial_i A_{ij} \partial_j u,
$$
so if the left-hand side is in $L^p$, and $u\in W^{2,p}$, then 
$\nabla u \in W^{1,p^*}$ with $p^*=np/(n-p)$ from Sobolev, and for the third term to be in $L^p$ it is natural to expect that 
$$
A \in W^{1,q} \cap L^\infty
$$ 
with, using H\"{o}lder's inequality
$$
\frac{1}{p}= \frac{n-p}{np} + \frac{1}{q},
$$
in other words $q = n$. So $A\in W^{1,n}\cap L^\infty$ is the reasonable general assumption (if you don't look at specific geometries as in Bunyang Li's answer).
Next, you can get rid of $f$: assuming that the domain $\Omega$ is sufficiently smooth, solve in $H^{1}_{0}(\Omega)$
$$
-\Delta \psi = f
$$
and you obtain that $\psi\in W^{2,p}(\Omega)$, and the problem writes
$$
\textrm{div}(A\nabla u + \nabla \psi)=0 \mbox{ in } \Omega, \quad(\star)
$$ 
do I have that at least locally, $(a\nabla u)\in W^{1,p}(\Omega)$, when $A \in W^{1,n} and bounded and coercive?

I think the optimal assumption is VMO and bounded (and elliptic), from Gallouet & Qafsaoui's 2002 paper. But I don't think it is written.

An easy case is $d=3$ and $A$ a scalar (not a matrix). Set $E=\nabla u $. From ($\star$) you obtain
$A E +\nabla \psi = \nabla \times H$ , with $\phi \in W^{1,p}$, or in other words
$$
E= A^{-1} \nabla \times H - A^{-1}\nabla \psi
$$. 
The second term is no problem, it is in $W^{1,p}$. Taking the divergence, as $A$ is a scalar, you obtain 
$$
\textrm{div} E = \nabla \times H \cdot \nabla A^{-1} +  \textrm{div}A^{-1}\nabla \psi
$$
Both terms are in $L^p$, so $W^{1,p}$(div). Of course, $\nabla \times E =0$. So you have $E$ in $W^{1,p}$(curl). 
Then you just need $W^{1,p}$(curl)$\cap W^{1,p}$(div)=$W^{1,p}$, which is true (at least locally, you have to be careful on the boundary).
