### The Helmholtz-Hodge decomposition

In the subject of vector calculus, Helmholtz's theorem states that any sufficiently smooth function in the unit ball can be expressed as a sum of a curl-free, a divergence-free, and a harmonic vector field [4]. In this Example, we show how this decomposition is computed by Ballfun and introduce the command HelmholtzDecomposition [3].

Let $v$ be a vector field defined in the unit ball. Helmholtz's theorem says that we can decompose $v$ as follows: $$\mathbf{v} = \nabla f + \nabla \times \psi + \nabla\phi,$$ where $f$ and $\phi$ are scalar-valued potential functions and $\psi$ is a vector field. The first term, $\nabla f$, is a gradient field and hence curl-free, while the second term, $\nabla\times\psi$, is solenoidal. The third term is a harmonic vector field (the vector Laplacian of $\nabla\phi$ is zero). From vector identities, one knows that the scalar field, $\phi$, is itself harmonic, i.e., $\Delta\phi=0$.

The Helmholtz-Hodge decomposition can be made unique by imposing additional constraints on $f$ and $\psi$ [4]. The standard constraints are: (1) $f$ is zero on the boundary of the unit ball, (2) the normal component of $\psi$ on the boundary is zero, and (3) $\psi$ is divergence-free. The HelmholtzDecomposition command in Ballfun computes the decomposition under these additional constraints.

The Helmholtz-Hogde decomposition is an important tool in fluid dynamics as it is used for compressible flow visualization, in CFD simulations (to impose the incompressibility condition), and in topological analysis. A survey of applications is available in [2].

### Calculating the decomposition

To explain the procedure for computing the decomposition, we take the following vector field:

v = ballfunv(@(x,y,z)cos(x.*y).*z,@(x,y,z)sin(x.*z),@(x,y,z)y.*z);
quiver( v )

### Computing the curl-free component

Since the divergence of a curl is zero and $\phi$ is harmonic, we know that the divergence of $\mathbf{v}$ is the Laplacian of $f$, i.e., $$\nabla \cdot \mathbf{v} = \nabla \cdot \nabla f = \nabla^2 f,$$ where the last equality holds because the divergence of the gradient is the Laplacian. Along with this, the zero Dirichlet conditions defines $f$.

f = poisson(div(v), @(lam,th)0, 50);
quiver( grad( f ) ), title('curl-free component of v')

We confirm that this component is curl-free:

norm( curl( grad( f ) ) )
ans =
9.772028777653057e-14


### Computing the harmonic component

Now, one can define a vector field $v^{(1)}$ as $$v^{(1)} = v - \nabla f = \nabla \times \psi + \nabla \phi.$$ Since we have the identity $$\vec{n} \cdot (\nabla \times \psi)|_{\partial{B}} = \psi \times \vec{n}|_{\partial B} = 0$$ and we want to find the harmonic function $\phi$, we solve the Laplace equation $$\Delta \phi = 0.$$ The Neumann boundary conditions are given by $$\vec{n} \cdot \nabla \phi|_{\partial B} = \frac{\partial \phi}{\partial r}|_{\partial B} = \vec{n} \cdot v^{(1)}|_{\partial B}.$$ Therefore, we can solve for $\psi$ in the Helmholtz-Hodge decomposition as follows:

v1 = v - grad(f);
bc = dot(spherefunv.unormal,v1(1,:,:,'spherical'));
phi = helmholtz(ballfun(0), 0, bc, 50, 'neumann');
quiver( grad( phi ) ), title('harmonic component of v')

We check the harmonicity of this component:

norm( laplacian( grad( phi ) ) )
ans =
3.527524381308947e-09


### Computing the divergence-free component

Let $v^{(2)}$ be the following vector field: $$v^{(2)} = v^{(1)} - \nabla \phi = \nabla \times \psi.$$ Since $v^{(2)}$ and $\psi$ are divergence-free, we can write their Poloidal-Toroidal decomposition [1] as $$v^{(2)} = \nabla\times\nabla\times(\mathbf{r}P_{v^{(2)}})+ \nabla\times(\mathbf{r}T_{v^{(2)}}),$$ $$\psi = \nabla\times\nabla\times(\mathbf{r}P_{\psi})+ \nabla\times(\mathbf{r}T_{\psi}),$$ where $\mathbf{r} = r\hat{r}$ and $\hat{r}$ denotes the unit radial vector in spherical coordinates. Moreover, the uniqueness of the PT decomposition and further vector identities lead us to the following system of equations for $\psi$: $$\Delta P_{\psi} = -T_{v^{(2)}}, \quad T_{\psi} = P_{v^{(2)}},$$ where $P_{\psi}$ is subjected to zero Dirichlet conditions because $\vec{n}\cdot\nabla\times(\mathbf{r}T_{\psi})|_{\partial B}=0$. Therefore, we can solve for $\psi$ in the Helmholtz-Hodge decomposition as follows:

v2 = v1 - grad(phi);
[Pv, Tv] = PTdecomposition(v2);
Ppsi = poisson(-Tv, @(lam,th)0, 50);
Tpsi = Pv;

We then recover the divergence-free vector field $\psi$ since $$\psi = \nabla\times\nabla\times(\mathbf{r}P_{\psi})+ \nabla\times(\mathbf{r}T_{\psi}).$$

psi = ballfunv.PT2ballfunv(Ppsi,Tpsi);
quiver( curl( psi ) ), title('divergence-free component of v')

By vector identities this component is divergence-free:

norm( div( curl( psi ) ) )
ans =
1.536895357752616e-10


### Visualizing the decomposition

Here are plots of each component of the decomposition:

subplot(2,2,1)
quiver( v ,'numpts',20), title('vector field')
subplot(2,2,2)
subplot(2,2,3)
quiver( curl(psi) ,'numpts',20), title('divergence-free')
subplot(2,2,4)
quiver( grad(phi) ,'numpts',20), title('harmonic')

As a sanity check we confirm that the decomposition has been successful:

w = grad( f ) + curl( psi ) + grad( phi );
norm( v - w )
ans =
7.815499404544524e-12


### The HelmholtzDecomposition command

Ballfun has a command HelmholtzDecomposition that computes the Helmholtz-Hodge decomposition of a vector field. Therefore, this example can be replicated with the following code:

[f, Ppsi, Tpsi, phi] = HelmholtzDecomposition( v );
psi = ballfunv.PT2ballfunv(Ppsi, Tpsi);

### References

[1] G. Backus, Poloidal and toroidal fields in geomagnetic field modelling, Reviews of Geophysics, 24 (1986), pp. 75-109.

[2] H. Bhatia, G. Norgard, V. Pascucci, and P.-T. Bremer, The Helmholtz-Hodge decomposition--A survey, IEEE Trans. Vis. Comput. Graphics, 19 (2013), pp. 1386-1404.

[3] N. Boullé and A. Townsend, Computing with functions on the ball, in preparation.

[4] Y. Tong, S. Lombeyda, A. Hirani, and M. Desbrun, Discrete multiscale vector field decomposition, ACM Trans. Graphics, 22 (2003), pp. 445-452.