Add ballfun examples

sphere/AtmosphericTemperature.m → sphere-ball/AtmosphericTemperature.m

@@ -2,7 +2,7 @@
 % Alex Townsend and Grady Wright, May 2016 
 
 %%
-% (Chebfun example sphere/AtmosphericTemperature.m)
+% (Chebfun example sphere-ball/AtmosphericTemperature.m)
 % [Tags: #spherefun]
 
 %% 1. Introduction 

sphere/Gravity.m → sphere-ball/Gravity.m

@@ -2,7 +2,7 @@
 % Nick Trefethen, May 2016
 
 %%
-% (Chebfun example sphere/Gravity.m)
+% (Chebfun example sphere-ball/Gravity.m)
 % [Tags: #spherefun]
 
 %% 1. A little geometry

sphere-ball/HelmholtzDecompositionBall.m

@@ -0,0 +1,178 @@
+%% Helmholtz-Hodge decomposition on the ball
+% Nicolas Boulle and Alex Townsend, May 2019 
+
+%% 
+% (Chebfun example sphere-ball/HelmholtzDecompositionBall.m)
+% [Tags: #ballfun, #vector-valued, #Helmholtz-Hodge] 
+
+%% The Helmholtz-Hodge decomposition 
+% In the area of vector calculus, Helmholtz's theorem states that any
+% sufficiently smooth in the ball can be expressed as a sum of a curl-free,
+% a divergence-free, and an 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 ball of radius $1$, 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 an harmonic vector field (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 constrains on $f$ and $\psi$ [4]. The standard constrains 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 constrains. 
+
+%% 
+% The Helmholtz-Hogde decomposition is an important tool in fluid dynamics as it
+% is used for flow visualization (when the fluid is compressible), in CFD
+% simulations (to impose that a fluid remains incompressible), and
+% topological analysis. A survey of applications is available here [2]. 
+
+%% Calculating the Helmholtz-Hodge 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 a gradient is the
+% Laplacian. Along with this, the zero Dirichlet conditions defines $f$. 
+
+f = poisson(div(v), @(lam,th)0, 50, 50, 50);
+quiver( grad( f ) ), title('Curl-free component of v')
+
+%% 
+% We confirm that this component is curl-free: 
+norm( curl( grad( f ) ) ) 
+
+%% 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 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);
+v1c = v1.comp;
+Vx = coeffs3(v1c{1},50); Vy = coeffs3(v1c{2},50); Vz = coeffs3(v1c{3},50);
+Vx = reshape(sum(Vx,1),50,50); Vy = reshape(sum(Vy,1),50,50); Vz = reshape(sum(Vz,1),50,50);
+MsinL = trigspec.multmat(50, [0.5i;0;-0.5i] ); McosL = trigspec.multmat(50, [0.5;0;0.5] );
+MsinT = trigspec.multmat(50, [0.5i;0;-0.5i] ); McosT = trigspec.multmat(50, [0.5;0;0.5] );
+v1_bc = McosL*Vx*MsinT.' + MsinL*Vy*MsinT.' + Vz*McosT.';
+phi = helmholtz(ballfun(0), 0, v1_bc, 50, 'neumann');
+quiver( grad( phi ) ), title('Harmonic component of v')
+
+%% 
+% We check the harmonicity of this component: 
+norm( laplacian( grad( phi ) ) ) 
+
+%% 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\vec{r}$. 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 ) ) ) 
+
+%% Visualizing the decomposition
+% Here is a plot of each component of the decomposition. 
+subplot(1,4,1) 
+quiver( v ), title('Vector field')
+subplot(1,4,2)
+quiver( grad(f) ), title('Curl-free')
+subplot(1,4,3)
+quiver( curl(psi) ), title('Divergence-free')
+subplot(1,4,4)
+quiver( grad(phi) ), title('Harmonic')
+
+%% 
+% As a sanity check we confirm that the decomposition has been successful:
+w = grad( f ) + curl( psi ) + grad( phi );
+norm( v - w ) 
+
+%% The HelmholtzDecomposition command 
+% Ballfun has a command called |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. Boulle, 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.

sphere/HelmholtzDecomposition.m → sphere-ball/HelmholtzDecompositionSphere.m

@@ -1,8 +1,8 @@
-%% Helmholtz-Hodge decomposition of a vector field
+%% Helmholtz-Hodge decomposition on the sphere
 % Alex Townsend and Grady Wright, May 2016 
 
 %% 
-% (Chebfun Example sphere/HelmholtzDecomposition.m)
+% (Chebfun example sphere-ball/HelmholtzDecompositionSphere.m)
 % [Tags: #spherefun, #vector-valued, #Helmholtz-Hodge] 
 
 %% 1. The Helmholtz-Hodge decomposition 

sphere-ball/PTdecomposition.m

@@ -0,0 +1,125 @@
+%% Poloidal-Toroidal decomposition of a vector field
+% Nicolas Boulle and Alex Townsend, May 2019 
+
+%% 
+% (Chebfun example sphere-ball/PTdecomposition.m)
+% [Tags: #ballfun, #vector-valued, #Poloidal-Toroidal] 
+
+%% The Poloidal-Toroidal decomposition
+% In vector calculus, the Poloidal-Toroidal (PT) decomposition [1] is a 
+% restricted form of the Helmholtz-Hodge decomposition [3]. It states that 
+% any sufficiently smooth and divergence-free vector field in the ball can 
+% be expressed as the sum of a poloidal field and a toroidal field. It is 
+% used in the analysis of divergence-free vector fields in geomagnetism, 
+% flow visualization, and incompressible fluid simulations. 
+
+%% 
+% For a chosen unit vector $\mathbf{r}$, a toroidal field, $\mathbf{T}$, is
+% one that is tangential to $\mathbf{r}$ while a poloidal field, $\mathbf{P}$, 
+% is one whose curl is tangential to $\mathbf{r}$, i.e., 
+%
+% $$ \mathbf{r}\cdot \mathbf{T} = 0, \quad \mathbf{r}\cdot (\nabla \times \mathbf{P}) = 0.$$
+
+%% The Poloidal-Toroidal decomposition for vectors defined in the unit ball
+% Let $w$ be any divergence-free vector field defined in the unit ball, 
+% then the PT decomposition says that $w$ can be written as the following 
+% sum:
+%
+% $$ w = \nabla\times\nabla\times(\mathbf{r}P_w) + \nabla\times(\mathbf{r}T_w),$$
+% 
+% where $P_w$ and $T_w$ are scalar-valued potential functions (called the 
+% poloidal and toroidal scalars). In this setting, the natural unit vector, 
+% $\mathbf{r}$, to select is the unit radial vector that points away from 
+% the origin. The scalars $P_w$ and T_w$ are unique up to the addition of 
+% an arbitrary function that only depends on the radial variable, $r$. 
+
+%%
+% In this Example, we show how this decomposition is computed by
+% Ballfun and introduce the command |PTdecomposition| [2]. 
+
+%% 
+% To illustrate the decomposition, we take the following divergence-free
+% vector field:
+Pw = ballfun(@(x,y,z)cos(x.*y));
+Tw = ballfun(@(x,y,z)sin(y.*z));
+w = ballfunv.PT2ballfunv(Pw, Tw);
+quiver( w )
+
+%%
+% We start by checking that $w$ is a divergence-free vector field:
+norm( div( w ) )
+
+%% Computing the Poloidal-Toroidal decomposition
+% The scalars $P_w$ and $T_w$ are unique up to the addition of 
+% an arbitrary function that only depends on the radial variable, $r$.
+% A popular additional constraint is to select them so that their integral
+% along each sphere of radius $1$ is zero.
+%
+% Under this additional restriction, they satisfy the following relations [1]:
+% 
+% $$\nabla_1^2 P_w = -rv_r,$$
+%
+% $$\nabla_1 T_w = -\Lambda_1\cdot w,$$
+%
+% where $\nabla_1^2$ and $\Lambda_1$ denote the dimensionless Laplacian 
+% and surface curl, respectively. These operations are defined in the spherical
+% coordinates $(r,\lambda,\theta)$ by
+%
+% $$\nabla_1^2 = \frac{1}{\sin\theta}\partial_\theta(\sin\theta\partial_\theta)
+% +\frac{1}{\sin^2\theta}\partial_\lambda^2,$$
+% 
+% $$\Lambda_1\cdot w = \frac{1}{\sin\theta}[\partial_\theta(w_\lambda\sin\theta)
+% -\partial_\lambda w_\theta].$$
+%
+% Moreover, these constraints are imposed in the solver by requiring that
+% the zeroth Fourier modes in $\lambda$ and $\theta$ are zero.
+% This actually means that if $P_w$ is represented by a
+% Chebyshev-Fourier-Fourier series, i.e.,
+%
+% $$\sum_{i=0}^n\sum_{j=-n/2}^{n/2}\sum_{k=-n/2}^{n/2}a_{i,j,k}T_i(r)e^{\mathbf{i}j\lambda}e^{\mathbf{i}k\theta},$$
+%
+% then $a_{i,0,0} = 0$ for $0\leq i\leq n$.
+
+%% The PTdecomposition command
+% Ballfun has a command called |PTdecomposition| that computes the 
+% Poloidal-Toroidal decomposition of a divergence-free vector field.
+% Therefore, the decomposition of $v$ can be computed as follows: 
+[Pw, Tw] = PTdecomposition( w );
+subplot(1,2,1)
+plot(Pw), title('Poloidal scalar')
+subplot(1,2,2)
+plot(Tw), title('Toroidal scalar')
+
+%% Visualizing the decomposition
+% Here is a plot of the decomposition.
+[P,T] = ballfunv.PT2ballfunv(Pw,Tw);
+subplot(1,3,1) 
+quiver( w ), title('Divergence-free vector field')
+subplot(1,3,2)
+quiver( P ), title('Poloidal component')
+subplot(1,3,3)
+quiver( T ), title('Toroidal component')
+
+%% Recover the vector field from the PT scalars
+% The original vector field can be recover from the poloidal and toroidal
+% scalars since
+%
+% $$ w = \nabla\times\nabla\times(\mathbf{r}P_w) + \nabla\times(\mathbf{r}T_w).$$
+%
+% This operation is implemented in Ballfun in the |PT2ballfunv| command
+v = ballfunv.PT2ballfunv(Pw, Tw);
+
+%% 
+% As a sanity check, we confirm that the decomposition has been successful:
+norm( v - w ) 
+
+%% References
+%
+% [1] G. Backus, Poloidal and toroidal fields in geomagnetic field modelling,
+% _Reviews of Geophysics_, 24 (1986), pp. 75-109.
+%
+% [2] N. Boulle, and A. Townsend, Computing with functions on the ball, in 
+% preparation.
+%
+% [3] Y. Tong, S. Lombeyda, A. Hirani, and M. Desbrun, Discrete Multiscale
+% Vector Field Decomposition, _ACM Trans. Graphics_, 22 (2003), pp. 445-452.

sphere/RayleighQuotientExample.m → sphere-ball/RayleighQuotientExample.m

@@ -2,7 +2,7 @@
 % Grady Wright, February 2017
 
 %%
-% (Chebfun example/sphere/RayleighQuotient.m )
+% (Chebfun example sphere-ball/RayleighQuotient.m)
 % [Tags: #spherefun, #spherefunv, #eigenvalues]
 MS = 'MarkerSize'; ms = 22; LW = 'LineWidth'; lw = 2;