I am working in a problem of fluid dynamics (diffusion flames, exactly) and I am struggling to prove the continuity of a variable. The function $Z:R^n\to R$ is continuous, but not necessarily smooth. Being $H(Z-Z_0)$ the Heaviside function, the equation for $Z$ is $$ H(Z-Z_0) \nabla \cdot (\mathbf{u} Z) = \nabla^2 Z, $$ with $\nabla \cdot \mathbf{u}=\sum_i \delta(\mathbf{r}-\mathbf{r}_i)$, it is, $\mathbf{u}$ is divergence-free almost everywhere, except for points inside the surface $Z=Z_0$ in which there are sources of $\mathbf{u}$.
What can be said about the continuity of $\nabla Z$ across the surface defined by $Z=Z_0$?
(If it's relevant, the boundary condition to $Z$ is piecewise continuous).
That is my attempt, assuming that $H(Z-Z_0)\nabla Z = \nabla [H(Z-Z_0) Z]$: in that case, the equation can be written as $$ \nabla \cdot (\mathbf{u} H(Z-Z_0) (Z-Z_1)) = \nabla^2 Z, $$ in which $Z_1$ is an arbitrary constant and, supposedly, does not afect the differential operator. Integrating the equation in the volume enclosed by $Z=C$, $$ \oint_{Z=C} H(Z-Z_0) (Z-Z_1) \mathbf{u}\cdot \mathbf{n} dS = \oint_{Z=C} \nabla Z \cdot \mathbf{n} dS. $$ Using $C=Z_0\pm \epsilon$ and taking the limit $\epsilon \to 0$ leads to $$ \oint_{Z=Z_0^+} (Z_0-Z_1) \mathbf{u} \cdot \mathbf{n} dS = \oint_{Z=Z_0^+} \nabla Z \cdot \mathbf{n} dS, $$ $$ 0 = \oint_{Z=Z_0^-} \nabla Z \cdot \mathbf{n} dS. $$ If $Z_1=Z_0$, $\nabla Z$ is continuous across the surface $Z=Z_0$; otherwise, it's discontinuous. However, if $Z_1$ could be chosen arbitrarily, it would make no difference to the behaviour of $\nabla Z$; if $Z_1=Z_0$ cannot be chosen, it's because $\nabla(Z-Z_0) \neq \nabla Z$ and $\nabla Z$ is not continuous at $Z=Z_0$.