From c9ca4212819fb33c4b7c6f156a5e3eee3adcf14f Mon Sep 17 00:00:00 2001 From: ajt60gaibb Date: Tue, 21 Aug 2018 09:27:19 -0400 Subject: [PATCH 1/4] Add Helmholtz decomposition example. --- ball/HelmholtzDecomposition.m | 137 ++++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 ball/HelmholtzDecomposition.m diff --git a/ball/HelmholtzDecomposition.m b/ball/HelmholtzDecomposition.m new file mode 100644 index 00000000..0c739650 --- /dev/null +++ b/ball/HelmholtzDecomposition.m @@ -0,0 +1,137 @@ +%% Helmholtz-Hodge decomposition of a vector field +% Nicolas BoullÈ and Alex Townsend, August 2018 + +%% +% (Chebfun Example ball/HelmholtzDecomposition.m) +% [Tags: #ballfun, #vector-valued, #Helmholtz-Hodge] + +%% 1. The Helmholtz-Hodge decomposition +% According to the Helmholtz-Hodge theorem [1], any smooth vector field in +% the ball can be uniquely decomposed into a sum of a divergence-free component, +% a curl-free component and a harmonic component. +% In this Example, we show how this decomposition can be computed by +% Ballfun and introduce the command |HelmholtzDecomposition| [2]. + +%% +% Let $v$ be any vector field in the ball, we can write it as the sum +% $$ \mathbf{v} = \nabla f + \nabla \times \psi + \nabla \phi, $$ +% where $f$ and $\phi$ are scalar-valued potential functions and $\psi$ is +% a vector field. Moreover, $f|_{\partial B} = 0$, +% $\psi\times\vec{n}|_{\partial{B}} = 0$, $\psi$ is divergence-free and +% $\phi$ is harmonic is $\Delta \phi = 0$. Here $\nabla$ is the gradient +% on the ball, $\nabla \times \psi$ stands for the divergence of $\psi$ and +% $\vec{n}$ is the unit normal vector. + +%% +% To illustrate the decomposition, we take the following +% vector field: +S = [40,40,40]; +vx = ballfun(@(x,y,z)cos(x.*y).*z,'cart',S); +vy = ballfun(@(x,y,z)sin(x.*z),'cart',S); +vz = ballfun(@(x,y,z)y.*z,'cart',S); +v = ballfunv(vx,vy,vz); +quiver( v ), view([-36 8]) + +%% 2. Computing the curl-free component +% Since the divergence of a curl is zero and $\phi$ is harmonic we know that +% $$ \nabla \cdot \mathbf{c} = \nabla \cdot \nabla \psi = \nabla^2 \phi, $$ +% where the last equality holds because the divergence of a gradient is the +% Laplacian. Therefore, we impose homogeneous boundary conditions on $f$ and +% solve the equation as follows: +f = helmholtz(div(v),0,zeros(40,40)); +quiver( grad( f ) ), hold on, +title('Curl-free component of v') +view([-36 8]), hold off + +%% +% We confirm that this component is curl-free: +norm( curl( grad( f ) ) ) + +%% 3. Computing the harmonic component +% We 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 component $\phi$, we solve the Laplace +% equation +% $$\Delta \phi = 0.$$ +% The Neumann boundary condition is 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: +v_1 = v - grad(f); +v_1_Boundary = ComputeNormalBoundary(v_1); +phi = helmholtz_neumann(ballfun(zeros(40,40,40)),0,v_1_Boundary); +quiver( grad( phi ) ), hold on, +title('harmonic component of v') +view([-36 8]), hold off + +%% +% We check the harmonicity of this component: +norm( laplace( grad( phi ) ) ) + +%% 3. 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 +% $$v^{(2)} = \nabla\times\nabla\times(\mathbb{r}P_{v^{(2)}})+ +% \nabla\times(\mathbb{r}T_{v^{(2)}}),$$ +% $$\psi = \nabla\times\nabla\times(\mathbb{r}P_{\psi})+ +% \nabla\times(\mathbb{r}T_{\psi}).$$ +% Moreover, the uniqueness of the PT decomposition and the 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 the first equation is subjected to the Dirichlet boundary condition +% $P_{\psi}|_{\partial B} = 0$ as +% $\vec{n}\cdot\nabla\times(\mathbb{r}T_{\psi})|_{\partial B}=0.$ +% Therefore, we can solve for $\psi$ in the Helmholtz-Hodge decomposition +% as follows: +v_2 = v_1 - grad(phi); +[Pv,Tv] = PTdecomposition(v_2); +Ppsi = helmholtz(-Tv,0,zeros(40,40)); +Tpsi = Pv; +% We then recover the divergence-free vector field $\psi$ since +% $$\psi = \nabla\times\nabla\times(\mathbb{r}P_{\psi})+ +% \nabla\times(\mathbb{r}T_{\psi}).$$ +psi = ballfunv.PT2ballfunv(Ppsi,Tpsi); +quiver( curl( psi ) ), hold on, +title('divergence-free component of v') +view([-36 8]), hold off + +%% +% By vector identities this component is divergence-free: +norm( div( curl( psi ) ) ) + +%% 4. Plotting the decomposition +% Here is a plot of the decomposition. +subplot(1,4,1) +quiver( v ), title('vector field'), view([-36 8]) +subplot(1,4,2) +quiver( grad(f) ), title('Curl-free'), view([-36 8]) +subplot(1,4,3) +quiver( curl(psi) ), title('Divergence-free'), view([-36 8]) +subplot(1,4,4) +quiver( grad(phi) ), title('Harmonic'), view([-36 8]) + +%% +% As a sanity check we confirm that the decomposition has been successful: +w = grad( f ) + curl( psi ) + grad(phi); +norm( v - w ) + +%% 5. 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); + +%% 6. References +% +% [1] Y. Tong, S. Lombeyda, A. Hirani, and M. Desbrun, Discrete Multiscale +% Vector Field Decomposition, _ACM Trans. Graphics_, 22 (2003), pp. 445-452. +% +% [2] N. BoullÈ, J. Slomka, and A. Townsend, Optimal complexity Navier-Stokes +% simulations in the ball, in preparation. From cfc9d90b8b4f5c8da9115d553466a0606f928536 Mon Sep 17 00:00:00 2001 From: Nicolas Boulle Date: Tue, 7 May 2019 13:11:20 +0100 Subject: [PATCH 2/4] Update ballfun examples --- ball/HelmholtzDecomposition.m | 201 ++++++++++++++++++++-------------- ball/PTdecomposition.m | 125 +++++++++++++++++++++ ball/SolidHarmonics.m | 164 +++++++++++++++++++++++++++ 3 files changed, 410 insertions(+), 80 deletions(-) create mode 100644 ball/PTdecomposition.m create mode 100644 ball/SolidHarmonics.m diff --git a/ball/HelmholtzDecomposition.m b/ball/HelmholtzDecomposition.m index 0c739650..431e4534 100644 --- a/ball/HelmholtzDecomposition.m +++ b/ball/HelmholtzDecomposition.m @@ -1,137 +1,178 @@ %% Helmholtz-Hodge decomposition of a vector field -% Nicolas BoullÈ and Alex Townsend, August 2018 +% Nicolas Boulle and Alex Townsend, March 2019 %% % (Chebfun Example ball/HelmholtzDecomposition.m) % [Tags: #ballfun, #vector-valued, #Helmholtz-Hodge] -%% 1. The Helmholtz-Hodge decomposition -% According to the Helmholtz-Hodge theorem [1], any smooth vector field in -% the ball can be uniquely decomposed into a sum of a divergence-free component, -% a curl-free component and a harmonic component. -% In this Example, we show how this decomposition can be computed by -% Ballfun and introduce the command |HelmholtzDecomposition| [2]. +%% 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 any vector field in the ball, we can write it as the sum -% $$ \mathbf{v} = \nabla f + \nabla \times \psi + \nabla \phi, $$ +% 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. Moreover, $f|_{\partial B} = 0$, -% $\psi\times\vec{n}|_{\partial{B}} = 0$, $\psi$ is divergence-free and -% $\phi$ is harmonic is $\Delta \phi = 0$. Here $\nabla$ is the gradient -% on the ball, $\nabla \times \psi$ stands for the divergence of $\psi$ and -% $\vec{n}$ is the unit normal vector. +% 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. %% -% To illustrate the decomposition, we take the following -% vector field: -S = [40,40,40]; -vx = ballfun(@(x,y,z)cos(x.*y).*z,'cart',S); -vy = ballfun(@(x,y,z)sin(x.*z),'cart',S); -vz = ballfun(@(x,y,z)y.*z,'cart',S); -v = ballfunv(vx,vy,vz); -quiver( v ), view([-36 8]) - -%% 2. Computing the curl-free component -% Since the divergence of a curl is zero and $\phi$ is harmonic we know that -% $$ \nabla \cdot \mathbf{c} = \nabla \cdot \nabla \psi = \nabla^2 \phi, $$ +% 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. Therefore, we impose homogeneous boundary conditions on $f$ and -% solve the equation as follows: -f = helmholtz(div(v),0,zeros(40,40)); -quiver( grad( f ) ), hold on, -title('Curl-free component of v') -view([-36 8]), hold off +% 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 ) ) ) -%% 3. Computing the harmonic component -% We define a vector field $v^{(1)}$ as +%% 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 component $\phi$, we solve the Laplace +% +% and want to find the harmonic function $\phi$, we solve the Laplace % equation +% % $$\Delta \phi = 0.$$ -% The Neumann boundary condition is given by +% +% 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: -v_1 = v - grad(f); -v_1_Boundary = ComputeNormalBoundary(v_1); -phi = helmholtz_neumann(ballfun(zeros(40,40,40)),0,v_1_Boundary); -quiver( grad( phi ) ), hold on, -title('harmonic component of v') -view([-36 8]), hold off +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( laplace( grad( phi ) ) ) +norm( laplacian( grad( phi ) ) ) -%% 3. Computing the divergence-free component +%% 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 -% $$v^{(2)} = \nabla\times\nabla\times(\mathbb{r}P_{v^{(2)}})+ -% \nabla\times(\mathbb{r}T_{v^{(2)}}),$$ -% $$\psi = \nabla\times\nabla\times(\mathbb{r}P_{\psi})+ -% \nabla\times(\mathbb{r}T_{\psi}).$$ -% Moreover, the uniqueness of the PT decomposition and the vector -% identities lead us to the following system of equations for $\psi$: +% 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 the first equation is subjected to the Dirichlet boundary condition -% $P_{\psi}|_{\partial B} = 0$ as -% $\vec{n}\cdot\nabla\times(\mathbb{r}T_{\psi})|_{\partial B}=0.$ +% +% 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: -v_2 = v_1 - grad(phi); -[Pv,Tv] = PTdecomposition(v_2); -Ppsi = helmholtz(-Tv,0,zeros(40,40)); +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(\mathbb{r}P_{\psi})+ -% \nabla\times(\mathbb{r}T_{\psi}).$$ +% +% $$\psi = \nabla\times\nabla\times(\mathbf{r}P_{\psi})+ +% \nabla\times(\mathbf{r}T_{\psi}).$$ +% psi = ballfunv.PT2ballfunv(Ppsi,Tpsi); -quiver( curl( psi ) ), hold on, -title('divergence-free component of v') -view([-36 8]), hold off +quiver( curl( psi ) ), title('Divergence-free component of v') %% % By vector identities this component is divergence-free: norm( div( curl( psi ) ) ) -%% 4. Plotting the decomposition -% Here is a plot of the decomposition. +%% Visualizing the decomposition +% Here is a plot of each component of the decomposition. subplot(1,4,1) -quiver( v ), title('vector field'), view([-36 8]) +quiver( v ), title('Vector field') subplot(1,4,2) -quiver( grad(f) ), title('Curl-free'), view([-36 8]) +quiver( grad(f) ), title('Curl-free') subplot(1,4,3) -quiver( curl(psi) ), title('Divergence-free'), view([-36 8]) +quiver( curl(psi) ), title('Divergence-free') subplot(1,4,4) -quiver( grad(phi) ), title('Harmonic'), view([-36 8]) +quiver( grad(phi) ), title('Harmonic') %% % As a sanity check we confirm that the decomposition has been successful: -w = grad( f ) + curl( psi ) + grad(phi); +w = grad( f ) + curl( psi ) + grad( phi ); norm( v - w ) -%% 5. 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); +%% 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); -%% 6. References +%% References: +% +% [1] G. Backus, Poloidal and toroidal fields in geomagnetic field modelling, +% _Reviews of Geophysics_, 24 (1986), pp. 75-109. % -% [1] Y. Tong, S. Lombeyda, A. Hirani, and M. Desbrun, Discrete Multiscale +% [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. -% -% [2] N. BoullÈ, J. Slomka, and A. Townsend, Optimal complexity Navier-Stokes -% simulations in the ball, in preparation. diff --git a/ball/PTdecomposition.m b/ball/PTdecomposition.m new file mode 100644 index 00000000..4146c18d --- /dev/null +++ b/ball/PTdecomposition.m @@ -0,0 +1,125 @@ +%% Poloidal-Toroidal decomposition of a vector field +% Nicolas Boulle and Alex Townsend, March 2019 + +%% +% (Chebfun Example 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. \ No newline at end of file diff --git a/ball/SolidHarmonics.m b/ball/SolidHarmonics.m new file mode 100644 index 00000000..3784de1e --- /dev/null +++ b/ball/SolidHarmonics.m @@ -0,0 +1,164 @@ +%% Solid harmonics +% Nicolas Boulle and Alex Townsend, March 2019 + +%% +% (Chebfun Example ball/SolidHarmonics.m) +% [Tags: #ballfun, #vector-valued, #Solid-Harmonics] + +%% Introduction +% Solid harmonics are eigenfunctions of the Laplace operator in spherical +% coordinates: +% +% $$\nabla^2\phi =\frac{1}{r^2}\left[\frac{\partial}{\partial +% r}\left(r^2\frac{\partial\phi}{\partial r}\right)+\frac{1}{\sin\theta} +% \frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial\phi}{\partial\theta}\right)+ +% \frac{1}{\sin^2\theta}\frac{\partial^2\phi}{\partial\lambda^2}\right] = 0.$$ +% +% This relationship holds because spherical harmonics are eigenfunctions of +% the surface Laplace (Laplace-Beltrami) operator, i.e., +% +% $$ \frac{1}{\sin\theta}\frac{\partial}{\partial\theta} +% \left(\sin\theta\frac{\partial Y^m_l}{\partial\theta}\right)+ \frac{1}{\sin^2\theta} +% \frac{\partial^2 Y^m_l}{\partial\lambda^2} = -l(l+1)Y^m_l,$$ +% +% where $Y^m_l$ stands for the normalized spherical harmonic of degree $l$ and order +% $m$. Substitution of $\phi=F(r)Y^m_l$ into the Laplace equation gives +% +% $$\frac{\partial}{\partial r}\left(r^2\frac{\partial FY^m_l}{\partial +% r}\right)=l(l+1)FY^m_l.$$ +% +% Therefore, general solutions to this equation are of the form $F(r) = Ar^l$ or +% $F(r)=Br^{-l-1}$. The particular solutions of the Laplace equations are +% regular solid harmonics, i.e., +% +% $$R^m_l(r,\lambda,\theta) = a_{lm}r^lY^m_l(\lambda,\theta),$$ +% +% which vanish at the origin, and irregular solid harmonics, i.e., +% +% $$I^m_l(r,\lambda,\theta) = a_{lm}\frac{Y^m_l(\lambda,\theta)}{r^{l+1}}.$$ + +%% +% In this example, we will use solid harmonics to mean the +% regular solid harmonics $R^m_l$ since the irregular solid harmonic possess a +% singularity at the origin. The solid harmonics are normalized so that their 2-norm +% is equal to $1$: +% +% $$\int_B R^m_l R^{m*}_ldV=1.$$ +% +% Thus, we have +% +% $$a_{lm}^2\int_0^1r^{2l}r^2dr\int_{\partial B}Y^m_lY^{m*}_ldS=1,$$ +% +% so that $a_{lm} = \sqrt{2l+3}$. + +%% Solid harmonics in Ballfun +% Solid harmonics can be constructed in Ballfun by calling the command +% |solharm|. This creates a solid harmonic of a given degree and order. For +% example, $R^2_4$ can be constructed and plotted as follows: +R42 = ballfun.solharm(4, 2); +plot( R42 ), axis off + +%% +% We can verify that this function is an eigenfunction of the Laplace operator +% with zero eigenvalue as follows: +norm( laplacian( R42 ) ) + +%% +% The solid harmonics are also orthonormal on the ball with respect to the +% standard $L^2$ inner-product. This can be verified with the |sum3| command: +R40 = ballfun.solharm(4, 0); +sum3( R42.*conj(R42) ) +sum3( R40.*conj(R40) ) +sum3( R42.*conj(R40) ) + +%% +% Here is a plot of the solid harmonics $R^m_l$, with $l=0,...,4$ and +% $0\leq m\leq l$. +N = 3; +for l = 0:N + for m = 0:l + R = ballfun.solharm(l,m); + subplot(N+1,N+1,l*(N+1)+m+1), plot(R) + axis off + end +end + +%% Computing solid harmonics coefficients +% A fast and stable algorithm for computing the solid harmonics is implemented +% in Ballfun. The computation of the solid harmonics $R^m_l$ of degree $l$ +% and order $m$ requires $\mathcal{O}(l\log l)$ operations. +% The solid harmonics $R^m_l$ can be expressed as +% +% $$R^m_l(r,\lambda,\theta) = \sqrt{2l+3}r^lP^m_l(\theta)e^(im\lambda),$$ +% +% where $P^m_l$ stands for the associated Legendre polynomial of degree $l$ +% and order $m$. + +%% +% The main issue in the computation of $R^m_l$ is to find the Fourier +% coefficients of $P^m_l$. +% Thus, the algorithm used in |solharm| to compute the coefficients of this +% polynomial is the Modified Forward Column (MFC) method described in [3]. +% The most popular recursive algorithm [1] (Forward Column method) that +% computes non-sectoral ($l>m$) $P^m_l$ from previously computed $P^m_{l-1}$ +% and $P^m_{l-2}$ is given by +% +% $$P^m_l(\theta) = a_{lm}\cos(\theta)P^m_{l-1}(\theta)+b_{lm}P^m_{l-1},\quad \forall l>m,$$ +% +% where +% +% $$P^m_l(\theta) = a_{lm}\cos(\theta)P^m_{l-1}(\theta)+b_{lm}P^m_{l-1},\quad \forall l>m,$$ +% +% and +% +% $$b_{lm} = \sqrt{\frac{(2l+1)(l+m+1)(n-m-1)}{(l-m)(l+m)(2l-3)}}.$$ +% +% The sectoral ($l=m$) $P^m_l$ can be computed using the initial values +% $P^0_0(\theta)=1$ and $P^1_1(\theta)=\sqrt(3)\sin(\theta)$ and the +% recursion [1] +% +% $$P^m_m(\theta)=\sin(\theta)\sqrt{\frac{2m+1}{2m}}P^{m-1}_{m-1},\quad +% \forall m>1.$$ +% +% Using these recursions, $P^m_l(\theta)$ is evaluated at $2l+1$ points and +% the coefficients are then recovered by an FFT. + +%% +% The recursive algorithm is unstable and will overflow for large +% degrees $l>1900$ [2] because of the factors $\sin(\theta)^m$ in $P^m_l$. +% The idea of the MFC method is to compute +% $P^m_l(\theta)/\sin(\theta)^m$ in the recursion and then multiply by +% $\sin(\theta)^m$ at the end before the FFT. +% The recursion to compute the sectoral values of $P^m_l/\sin(\theta)^m$ +% remains unchanged: +% +% $$\frac{P^m_l(\theta)}{\sin(\theta)^m} = a_{lm}cos(\theta) +% \frac{P^m_{l-1}(\theta)}{\sin(\theta)^m}+b_{lm} +% \frac{P^m_{l-1}}{\sin(\theta)^m},\quad \forall l>m.$$ +% +% Finally, the sectoral values of $P^m_m(\theta)/\sin(\theta)^m$ are +% computed using the initial values $P^0_0(\theta)/\sin(\theta)^0 = 1$ and +% $P^1_1(\theta)/\sin(\theta) = \sqrt(3)$ and the relationship +% +% $$P^m_m(\theta)=\sqrt{\frac{2m+1}{2m}}P^{m-1}_{m-1},\quad \forall m>1.$$ + +%% +% The computation of the solid harmonics is very fast in Ballfun and a +% solid harmonics of degree $150$ can be computed in a few tenths of +% a second: +tic +ballfun.solharm(150, 50); +toc + +%% References: +% +% [1] C. Colombo, Numerical methods for harmonic analysis on the sphere, +% _Department of Geodetic Science and Surveying_, (1981), The Ohio State +% University, Columbus. +% +% [2] D. Gleason, Partial sums of Legendre series via Clenshaw +% summation, _Manuscr Geod_, 10 (1985), pp. 115-130. +% +% [3] S. Holmes, and W. Featherstone, A unified approach to the Clenshaw +% summation and the recursive computation of very high degree and order +% normalised associated Legendre functions, _Journal of Geodesy_, 76 (2002), pp 279-299. \ No newline at end of file From 5ffffc866e1e3357afb356920ad322d64007cab6 Mon Sep 17 00:00:00 2001 From: Nicolas Boulle Date: Mon, 13 May 2019 18:53:28 +0100 Subject: [PATCH 3/4] Update date --- ball/HelmholtzDecomposition.m | 2 +- ball/PTdecomposition.m | 2 +- ball/SolidHarmonics.m | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ball/HelmholtzDecomposition.m b/ball/HelmholtzDecomposition.m index 431e4534..3cdbcaff 100644 --- a/ball/HelmholtzDecomposition.m +++ b/ball/HelmholtzDecomposition.m @@ -1,5 +1,5 @@ %% Helmholtz-Hodge decomposition of a vector field -% Nicolas Boulle and Alex Townsend, March 2019 +% Nicolas Boulle and Alex Townsend, May 2019 %% % (Chebfun Example ball/HelmholtzDecomposition.m) diff --git a/ball/PTdecomposition.m b/ball/PTdecomposition.m index 4146c18d..cd4d4a1f 100644 --- a/ball/PTdecomposition.m +++ b/ball/PTdecomposition.m @@ -1,5 +1,5 @@ %% Poloidal-Toroidal decomposition of a vector field -% Nicolas Boulle and Alex Townsend, March 2019 +% Nicolas Boulle and Alex Townsend, May 2019 %% % (Chebfun Example ball/PTdecomposition.m) diff --git a/ball/SolidHarmonics.m b/ball/SolidHarmonics.m index 3784de1e..db5e8ad1 100644 --- a/ball/SolidHarmonics.m +++ b/ball/SolidHarmonics.m @@ -1,5 +1,5 @@ %% Solid harmonics -% Nicolas Boulle and Alex Townsend, March 2019 +% Nicolas Boulle and Alex Townsend, May 2019 %% % (Chebfun Example ball/SolidHarmonics.m) From 676ee0e21b2ed72b93df9466ff6e5a007b19d18b Mon Sep 17 00:00:00 2001 From: Nicolas Boulle Date: Tue, 14 May 2019 09:39:09 +0100 Subject: [PATCH 4/4] Merge sphere and ball examples --- {sphere => sphere-ball}/AtmosphericData.mat | Bin {sphere => sphere-ball}/AtmosphericTemperature.m | 2 +- {sphere => sphere-ball}/Gravity.m | 2 +- .../HelmholtzDecompositionBall.m | 4 ++-- .../HelmholtzDecompositionSphere.m | 4 ++-- {ball => sphere-ball}/PTdecomposition.m | 2 +- {sphere => sphere-ball}/RayleighQuotientExample.m | 2 +- {ball => sphere-ball}/SolidHarmonics.m | 2 +- {sphere => sphere-ball}/SphereHeatConduction.m | 2 +- {sphere => sphere-ball}/SpherefunPartition.m | 2 +- {sphere => sphere-ball}/SpherefunRotate.m | 2 +- {sphere => sphere-ball}/SphericalHarmonics.m | 2 +- 12 files changed, 13 insertions(+), 13 deletions(-) rename {sphere => sphere-ball}/AtmosphericData.mat (100%) rename {sphere => sphere-ball}/AtmosphericTemperature.m (99%) rename {sphere => sphere-ball}/Gravity.m (97%) rename ball/HelmholtzDecomposition.m => sphere-ball/HelmholtzDecompositionBall.m (98%) rename sphere/HelmholtzDecomposition.m => sphere-ball/HelmholtzDecompositionSphere.m (97%) rename {ball => sphere-ball}/PTdecomposition.m (96%) rename {sphere => sphere-ball}/RayleighQuotientExample.m (99%) rename {ball => sphere-ball}/SolidHarmonics.m (96%) rename {sphere => sphere-ball}/SphereHeatConduction.m (99%) rename {sphere => sphere-ball}/SpherefunPartition.m (97%) rename {sphere => sphere-ball}/SpherefunRotate.m (99%) rename {sphere => sphere-ball}/SphericalHarmonics.m (99%) diff --git a/sphere/AtmosphericData.mat b/sphere-ball/AtmosphericData.mat similarity index 100% rename from sphere/AtmosphericData.mat rename to sphere-ball/AtmosphericData.mat diff --git a/sphere/AtmosphericTemperature.m b/sphere-ball/AtmosphericTemperature.m similarity index 99% rename from sphere/AtmosphericTemperature.m rename to sphere-ball/AtmosphericTemperature.m index b2e5f8f3..916ecbdd 100644 --- a/sphere/AtmosphericTemperature.m +++ b/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 diff --git a/sphere/Gravity.m b/sphere-ball/Gravity.m similarity index 97% rename from sphere/Gravity.m rename to sphere-ball/Gravity.m index a59c6962..224b2909 100644 --- a/sphere/Gravity.m +++ b/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 diff --git a/ball/HelmholtzDecomposition.m b/sphere-ball/HelmholtzDecompositionBall.m similarity index 98% rename from ball/HelmholtzDecomposition.m rename to sphere-ball/HelmholtzDecompositionBall.m index 3cdbcaff..a792a99d 100644 --- a/ball/HelmholtzDecomposition.m +++ b/sphere-ball/HelmholtzDecompositionBall.m @@ -1,8 +1,8 @@ -%% Helmholtz-Hodge decomposition of a vector field +%% Helmholtz-Hodge decomposition on the ball % Nicolas Boulle and Alex Townsend, May 2019 %% -% (Chebfun Example ball/HelmholtzDecomposition.m) +% (Chebfun example sphere-ball/HelmholtzDecompositionBall.m) % [Tags: #ballfun, #vector-valued, #Helmholtz-Hodge] %% The Helmholtz-Hodge decomposition diff --git a/sphere/HelmholtzDecomposition.m b/sphere-ball/HelmholtzDecompositionSphere.m similarity index 97% rename from sphere/HelmholtzDecomposition.m rename to sphere-ball/HelmholtzDecompositionSphere.m index e4fbc04d..d75f40a0 100644 --- a/sphere/HelmholtzDecomposition.m +++ b/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 diff --git a/ball/PTdecomposition.m b/sphere-ball/PTdecomposition.m similarity index 96% rename from ball/PTdecomposition.m rename to sphere-ball/PTdecomposition.m index cd4d4a1f..82a337d3 100644 --- a/ball/PTdecomposition.m +++ b/sphere-ball/PTdecomposition.m @@ -2,7 +2,7 @@ % Nicolas Boulle and Alex Townsend, May 2019 %% -% (Chebfun Example ball/PTdecomposition.m) +% (Chebfun example sphere-ball/PTdecomposition.m) % [Tags: #ballfun, #vector-valued, #Poloidal-Toroidal] %% The Poloidal-Toroidal decomposition diff --git a/sphere/RayleighQuotientExample.m b/sphere-ball/RayleighQuotientExample.m similarity index 99% rename from sphere/RayleighQuotientExample.m rename to sphere-ball/RayleighQuotientExample.m index 9c705519..bb7b9038 100644 --- a/sphere/RayleighQuotientExample.m +++ b/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; diff --git a/ball/SolidHarmonics.m b/sphere-ball/SolidHarmonics.m similarity index 96% rename from ball/SolidHarmonics.m rename to sphere-ball/SolidHarmonics.m index db5e8ad1..633d1933 100644 --- a/ball/SolidHarmonics.m +++ b/sphere-ball/SolidHarmonics.m @@ -2,7 +2,7 @@ % Nicolas Boulle and Alex Townsend, May 2019 %% -% (Chebfun Example ball/SolidHarmonics.m) +% (Chebfun example sphere-ball/SolidHarmonics.m) % [Tags: #ballfun, #vector-valued, #Solid-Harmonics] %% Introduction diff --git a/sphere/SphereHeatConduction.m b/sphere-ball/SphereHeatConduction.m similarity index 99% rename from sphere/SphereHeatConduction.m rename to sphere-ball/SphereHeatConduction.m index 67d21c2a..f07ef1e1 100644 --- a/sphere/SphereHeatConduction.m +++ b/sphere-ball/SphereHeatConduction.m @@ -2,7 +2,7 @@ % Alex Townsend and Grady Wright, May 2016 %% -% (Chebfun example spherefun) +% (Chebfun example sphere-ball/SphereHeatConduction.m) % [Tags: #spherefun, #heat, #diffusion] %% 1. Introduction diff --git a/sphere/SpherefunPartition.m b/sphere-ball/SpherefunPartition.m similarity index 97% rename from sphere/SpherefunPartition.m rename to sphere-ball/SpherefunPartition.m index f8a95799..356e66c9 100644 --- a/sphere/SpherefunPartition.m +++ b/sphere-ball/SpherefunPartition.m @@ -2,7 +2,7 @@ % Behnam Hashemi, November 2016 %% -% (Chebfun Example sphere/SpherefunPartition.m) +% (Chebfun example sphere-ball/SpherefunPartition.m) % [Tags: #Spherefun] %% diff --git a/sphere/SpherefunRotate.m b/sphere-ball/SpherefunRotate.m similarity index 99% rename from sphere/SpherefunRotate.m rename to sphere-ball/SpherefunRotate.m index 506d7905..84496a81 100644 --- a/sphere/SpherefunRotate.m +++ b/sphere-ball/SpherefunRotate.m @@ -2,7 +2,7 @@ % Alex Townsend and Grady Wright, May 2017 %% -% (Chebfun example spherefun) +% (Chebfun example sphere-ball/SpherefunRotate.m) % [Tags: #spherefun, #rotate, #fourier] %% Introduction diff --git a/sphere/SphericalHarmonics.m b/sphere-ball/SphericalHarmonics.m similarity index 99% rename from sphere/SphericalHarmonics.m rename to sphere-ball/SphericalHarmonics.m index e5e01506..2df7bbda 100644 --- a/sphere/SphericalHarmonics.m +++ b/sphere-ball/SphericalHarmonics.m @@ -2,7 +2,7 @@ % Alex Townsend and Grady Wright, May 2016 %% -% (Chebfun example sphere/SphericalHarmonics.m) +% (Chebfun example sphere-ball/SphericalHarmonics.m) % [Tags: #spherefun, #eigenfunctions, #fourier] %% 1. Introduction