I think I may have some rough idea about this type of problem, probably I am wrong though, should be a comment really.
If in $\mathbb{R}^2$, $u_{xy}\neq u_{yx}$, while the first derivatives exists, in weak sense, $u\in H^1$ but $u\notin H^2$, and the problem can be rewritten as: $ \tag{1}\mathrm{div}(R \nabla u) = 0 $ where $R = \begin{pmatrix}0&1 \\1&0\end{pmatrix}$ is the reflection transformation with respect to the line $y=x$, multiply both sides of (1) by a smooth test function $v\in C^{\infty}_0$, suppose $u$ is a piecewisely defined function:
$ \int_{\mathbb{R}^2} \mathrm{div}(R \nabla u) v = -\int_{\mathbb{R}^2} R \nabla u\cdot \nabla v + \sum_{\Omega_i}\int_{\partial \Omega_i} vR\nabla u\cdot \nu $ The first term is gone if $u$ is a weak solution, for second term, if the reflection transformation would make the normal component across the boundary still continuous, ie $(vR\nabla u\cdot \nu|_{\partial \Omega^+} + vR\nabla u\cdot \nu|_{\partial \Omega^-})$ be zero.
A possible candidate $u$ would a piecewise linear function that is $1$ on $y=x$ and $0$ on $y = x\pm 1$, has linear decay within the stripe and $0$ elsewhere, $\mathbb{R}^2$ is being divided into 2 parts, where a.e. $u_{xy} = 0 = u_{yx}$, while the normal component of $R\nabla u$ is continuous across $y=x$, hence we have $\mathrm{div}(R\nabla u)=0$.
I tried for an hour and can't think of an explicit example of $u_{xy} = -u_{yx}$ on a positive Lebesgue measure set, maybe it is as LVK suggested.