Supose we have $u \in W^{1,p}$ (i.e $u$ has weak partial derivatives, which we denote by $Du$), and that both $u$ and $Du$ are continuous (More precisely, there is a continuous representative in the equivalence class of $u$ and $Du$ that is continuous).
Can we ensure that $u \in C^1$?
I mean, Can we ensure that $u$ (The continuous representative) has classical partial derivatives and that they are continuous?
I am working with the Variational Approach and at a certain step it is required to show that $u$ is $C^2$ and the results I have prove that there is a continuous representative for all the derivatives of $u$ and $u$, which seems to be not enough for me.
Let $\Omega \subset \mathbb{R}^n$ open. If $u \in \mathscr{C}(\Omega)$, and the distributional derivatives of $u$ are given by integration against continuous functions, that is, there exist $v_1,\dotsc, v_n \in \mathscr{C}(\Omega)$ such that for all $\varphi \in \mathscr{D}(\Omega)$ we have
$$- \int_{\Omega} u(x)D_k\varphi(x)\,dx = \int_{\Omega} v_k(x)\varphi(x)\,dx$$
for $1 \leqslant k \leqslant n$, then $u \in \mathscr{C}^1(\Omega)$.
Weak derivatives, when they exist, coincide with distributional derivatives, so this includes your setting.
Proof sketch: Let $\psi \in \mathscr{D}(\mathbb{R}^n)$ with $\psi \geqslant 0$, $\operatorname{supp} \psi \subset \overline{B_1(0)}$ and $\int_{\mathbb{R}^n} \psi(x)\,dx = 1$. For $\varepsilon > 0$, define $\psi_{\varepsilon}(x) := \varepsilon^{-n}\psi(x/\varepsilon)$. For $\delta > 0$ and $0 < \varepsilon < \delta$, let
$$u_{\varepsilon}(x) := (\psi_{\varepsilon} \ast u)(x) = \int_{\Omega} u(x-y) \psi_{\varepsilon}(y)\,dy = \int_{B_1(0)} u(x-\varepsilon y) \psi(y)\,dy$$
for $x \in \Omega_{\delta} = \{x \in \Omega : d(x,\partial \Omega) > \delta\}$.
Then $u_{\varepsilon} \in \mathscr{C}^{\infty}(\Omega_{\delta})$, and
$$D_k u_{\varepsilon} = D_k (\psi_{\varepsilon}\ast u) = (D_k \psi_{\varepsilon}) \ast u = \psi_{\varepsilon} \ast v_k. \tag{1}$$
Since $u, v_1,\dotsc, v_n \in \mathscr{C}(\Omega)$, it follows that $u_{\varepsilon} \to u$, and $D_k u_{\varepsilon} \to v_k$ locally uniformly on $\Omega_{\delta}$ as $\varepsilon \to 0$.
Then fix $x^{(0)} \in \Omega_{\delta}$, and note that for $0 < \eta < d(x^{(0)}, \partial \Omega) - \delta$ we have
\begin{align} u(x^{(0)} + he_k) - u(x^{(0)}) &= \lim_{\varepsilon \to 0} u_{\varepsilon}(x^{(0)} + he_k) - u_{\varepsilon}(x^{(0)})\\ &= \lim_{\varepsilon \to 0} \int_0^h D_ku_{\varepsilon}(x^{(0)} + te_k)\,dt\\ &= \int_0^h v_k(x^{(0)} + te_k)\,dt \end{align}
for $\lvert h\rvert < \eta$, which shows that $u$ is partially differentiable at $x^{(0)}$ with $\dfrac{\partial u}{\partial x_k}(x^{(0)}) = v_k(x^{(0)})$.
Since $\delta > 0$ was arbitrary, it follows that $u$ is partially differentiable in the classical sense on all of $\Omega$, and the partial derivatives are continuous by assumption, so $u \in \mathscr{C}^1(\Omega)$.