Solution (source code)

= Solution
{c}

Choose $x$ at which $Du(x)$ is finite and the $p$-mean version of the <Lebesgue differentiation theorem> holds:
$$
\varepsilon_x(r):=\left(\frac1{|Q_r(x)|}\int_{Q_r(x)}|Du(z)-Du(x)|^p dz\right)^{1/p}\longrightarrow0.
$$
These conditions hold <almost everywhere>. On a cube about $x$, subtract the affine function and put $v(y)=u^*(y)-u^*(x)-Du(x)\cdot(y-x)$. Its <weak derivative> is $Dv=Du-Du(x)$ and $v(x)=0$. Approximate $u$ as in part (b) and subtract the same affine function from the approximants. <Uniform convergence> and convergence of their <weak derivatives> extend part (a)'s estimate to $v$, giving
$$
|v(y)|\leq Cr^{1-n/p}\|Du-Du(x)\|_{L^p(Q_r(x))}
=Cr\varepsilon_x(r)\qquad(y\in Q_r(x)).
$$
For $h\ne0$, choose $r=4|h|$ and $y=x+h$. It follows that
$$
\frac{|u^*(x+h)-u^*(x)-Du(x)\cdot h|}{|h|}
\leq4C\varepsilon_x(4|h|)\longrightarrow0.
$$
This is precisely classical <Fréchet derivative> differentiability, not merely existence of coordinate derivatives. Therefore the <differentiability almost everywhere of supercritical Sobolev functions> gives
$$
\boxed{D_{\mathrm{classical}}u^*(x)=D_{\mathrm{weak}}u(x)\quad\text{for almost every }x.}
$$