Spectral Numerical Exterior Calculus Methods for Differential Equations on Radial Manifolds

TL;DR

This paper introduces spectral exterior calculus methods for solving PDEs on radial manifolds, achieving high-precision discretizations of differential operators using spherical harmonics and Lebedev quadrature.

Contribution

It develops a spectral exterior calculus framework with hyperinterpolation techniques for accurate numerical solutions on radial manifolds, including convergence analysis and applications to Laplace-Beltrami equations.

Findings

Spectral convergence of numerical exterior derivative and Hodge star operators.

High-precision discretizations achieved with spherical harmonics and Lebedev quadrature.

Successful application to Laplace-Beltrami equations on radial manifolds.

Abstract

We develop exterior calculus approaches for partial differential equations on radial manifolds. We introduce numerical methods that approximate with spectral accuracy the exterior derivative $\mathbf{d}$, Hodge star $\star$, and their compositions. To achieve discretizations with high precision and symmetry, we develop hyperinterpolation methods based on spherical harmonics and Lebedev quadrature. We perform convergence studies of our numerical exterior derivative operator $\overline{\mathbf{d}}$ and Hodge star operator $\overline{\star}$ showing each converge spectrally to $\mathbf{d}$ and $\star$. We show how the numerical operators can be naturally composed to formulate general numerical approximations for solving differential equations on manifolds. We present results for the Laplace-Beltrami equations demonstrating our approach.

Figures

Equations

$$♭:v^j{\partial }_{{\mathbf{x}}^j}\to v_{i}d{\mathbf{x}}^i,♯:v_{i}d{\mathbf{x}}^i\to v^j{\partial }_{{\mathbf{x}}^j}.$$

$$\mathbf{d}\alpha =\frac{1}{k!}\frac{\partial {\alpha }_{i_1\dots i_k}}{\partial x^j}\mathbf{d}{\mathbf{x}}^j\wedge \mathbf{d}{\mathbf{x}}^{i_1}\wedge \cdots \wedge \mathbf{d}{\mathbf{x}}^{i_k}.$$

$$\star \alpha =\frac{\sqrt{|g|}}{(n-k)!k!}{\alpha }^{i_1\dots i_k}{ϵ}_{i_1\dots i_k{j}_1\dots j_{n-k}}\cdot \mathbf{d}{\mathbf{x}}^{j_1}\wedge \cdots \wedge \mathbf{d}{\mathbf{x}}^{j_{n-k}},$$

$$\text{mbox}grad(f)=[\mathbf{d}f{]}^{♯},\text{mbox}div(\mathbf{F})=-(\star \mathbf{d}\star {\mathbf{F}}^{♭})=-\mathbf{\delta }{\mathbf{F}}^{♭},\text{mbox}curl(\mathbf{F})={\left[\star (\mathbf{d}{\mathbf{F}}^{♭})\right]}^{♯},$$

$${\Delta }^H(\mathbf{F})=-{\left[\left(\mathbf{\delta }\mathbf{d}+\mathbf{d}\mathbf{\delta }\right){\mathbf{F}}^{♭}\right]}^{♯},{\Delta }^S(\mathbf{F})=-{\left[\mathbf{\delta }\mathbf{d}{\mathbf{F}}^{♭}\right]}^{♯},{\Delta }^{H}f={\Delta }^{R}f=-(\star \mathbf{d}\star )\mathbf{d}f=-\mathbf{\delta }\mathbf{d}f.$$

$${\int }_{\partial \Omega }\omega ={\int }_{\Omega }\mathbf{d}\omega ,{\int }_{\partial \Omega }\star \omega ={\int }_{\Omega }\mathbf{d}\star \omega .$$

$$\frac{\partial }{\partial t}{\int }_{\Omega }\star u={\int }_{\partial \Omega }\star \omega ={\int }_{\Omega }\mathbf{d}\star \omega .$$

$$\frac{\partial u}{\partial t}=-\star \mathbf{d}\star \omega =-\mathbf{\delta }\omega .$$

$$\omega =\mathbf{d}u.$$

$$\frac{\partial u}{\partial t}=-\mathbf{\delta }\mathbf{d}u.$$

$$\mathcal{P}\left[f\right]=\bar{f}(\theta ,\varphi )=\sum _i{\hat{f}}_i{Y}_i(\theta ,\varphi ),$$

$$⟨u,v{⟩}_Q=\sum _{\ell }{w}_{\ell }u({\mathbf{x}}_{\ell })v({\mathbf{x}}_{\ell }).$$

$${\bar{\mathbf{v}}}^{♯}(\theta ,\varphi )=[\mathcal{P}{\bar{v}}^x,\mathcal{P}{\bar{v}}^y,\mathcal{P}{\bar{v}}^z].$$

$${\partial }_s{v}_i=({\partial }_s{g}_{ij})a_k^j{\bar{v}}^k+g_{ij}({\partial }_s{a}_k^j){\bar{v}}^k+g_{ij}{a}_k^j({\partial }_s{\bar{v}}^k).$$

$$\mathcal{L}u=-g.$$

$$⟨\tilde{\mathcal{L}}\bar{u},Y_i{⟩}_Q=-⟨\bar{g},Y_i{⟩}_Q.$$

$$K\hat{\mathbf{u}}$$

$$\mathcal{L}u=-\mathbf{\delta }\mathbf{d}u=-g.$$

$$Y_n^m(\theta ,\varphi )=\sqrt{\frac{(2n+1)(n-m)!}{4\pi (n+m)!}}{P}_n^m\left(\cos (\varphi )\right)\exp \left(im\theta \right)$$

$$Y_n^m(\theta ,\varphi )=X_n^m(\theta ,\varphi )+i{Z}_n^m(\theta ,\varphi )$$

$${\tilde{Y}}_1=Y_0^0,{\tilde{Y}}_2=Z_1^1,{\tilde{Y}}_3=Y_1^0,{\tilde{Y}}_4=X_1^1,{\tilde{Y}}_5=Z_2^2,{\tilde{Y}}_6=Z_2^1,{\tilde{Y}}_7=Y_2^0,{\tilde{Y}}_8=X_2^1,{\tilde{Y}}_9=X_2^2.$$

$${\partial }_{\theta }{Y}_n^m(\theta ,\varphi )={\partial }_{\theta }\sqrt{\frac{(2n+1)(n-m)!}{4\pi (n+m)!}}{P}_n^m(\cos (\varphi ))\exp \left(im\theta \right)=im{Y}_n^m\left(\theta ,\varphi \right).$$

$${\partial }_{\varphi }{Y}_n^m(\theta ,\varphi )=m\cot (\varphi )Y_n^m(\theta ,\varphi )+\sqrt{(n-m)(n+m+1)}\exp \left(-i\theta \right)Y_n^{m+1}(\theta ,\varphi ).$$

$$\mathbf{x}(\theta ,\varphi )=\mathbf{\sigma }(\theta ,\varphi )=r(\theta ,\varphi )\mathbf{r}(\theta ,\varphi )$$

$$\mathbf{x}(\hat{\theta },\hat{\varphi })=r(\hat{\theta },\hat{\varphi })\mathbf{r}(\hat{\theta },\hat{\varphi }),\mathbf{r}(\hat{\theta },\hat{\varphi })=\left[\begin{array}{ccc}\sin (\hat{\varphi })\cos (\hat{\theta }),& \sin (\hat{\varphi })\sin (\hat{\theta }),& \cos (\hat{\varphi })\end{array}\right]$$

$$\mathbf{x}(\bar{\theta },\bar{\varphi })=r(\bar{\theta },\bar{\varphi })\mathbf{r}(\bar{\theta },\bar{\varphi }),\bar{\mathbf{r}}(\bar{\theta },\bar{\varphi })=\left[\begin{array}{ccc}\cos (\bar{\varphi }),& \sin (\bar{\varphi })\sin (\bar{\theta }),& -\sin (\bar{\varphi })\cos (\bar{\theta })\end{array}\right].$$

$${\mathbf{\sigma }}_{\theta }(\theta ,\varphi )$$

$${\mathbf{\sigma }}_{\varphi }(\theta ,\varphi )$$

$$\mathbf{I}=\left[\begin{array}{cc}E& F\\ F& G\end{array}\right]=\left[\begin{array}{cc}{\mathbf{\sigma }}_{\theta }\cdot {\mathbf{\sigma }}_{\theta }& {\mathbf{\sigma }}_{\theta }\cdot {\mathbf{\sigma }}_{\varphi }\\ {\mathbf{\sigma }}_{\varphi }\cdot {\mathbf{\sigma }}_{\theta }& {\mathbf{\sigma }}_{\varphi }\cdot {\mathbf{\sigma }}_{\varphi }\end{array}\right]=\left[\begin{array}{cc}{r}_{\theta }^2+r^2\sin (\varphi {)}^2& r_{\theta }{r}_{\varphi }\\ r_{\theta }{r}_{\varphi }& r_{\varphi }^2+r^2\end{array}\right].$$

$$\mathbf{II}=\left[\begin{array}{cc}L& M\\ N& N\end{array}\right]=\left[\begin{array}{cc}{\mathbf{\sigma }}_{\theta \theta }\cdot \mathbf{n}& {\mathbf{\sigma }}_{\theta \varphi }\cdot \mathbf{n}\\ {\mathbf{\sigma }}_{\varphi \theta }\cdot \mathbf{n}& {\mathbf{\sigma }}_{\varphi \varphi }\cdot \mathbf{n}\end{array}\right].$$

$$\mathbf{n}(\theta ,\varphi )=\frac{{\mathbf{\sigma }}_{\theta }(\theta ,\varphi )\times {\mathbf{\sigma }}_{\varphi }(\theta ,\varphi )}{\Vert {\mathbf{\sigma }}_{\theta }(\theta ,\varphi )\times {\mathbf{\sigma }}_{\varphi }(\theta ,\varphi )\Vert }.$$

Full text

Spectral Numerical Exterior Calculus Methods for Differential Equations on Radial Manifolds

B. Gross and P. J. Atzberger*∗*

W****e develop exterior calculus approaches for partial differential equations on radial manifolds. We introduce numerical methods that approximate with spectral accuracy the exterior derivative $\mathbf{d}$, Hodge star $\star$, and their compositions. To achieve discretizations with high precision and symmetry, we develop hyperinterpolation methods based on spherical harmonics and Lebedev quadrature. We perform convergence studies of our numerical exterior derivative operator $\overline{\mathbf{d}}$ and Hodge star operator $\overline{\star}$ showing each converge spectrally to $\mathbf{d}$ and $\star$. We show how the numerical operators can be naturally composed to formulate general numerical approximations for solving differential equations on manifolds. We present results for the Laplace-Beltrami equations demonstrating our approach.

1 Introduction

There has been a lot of recent interest in numerical methods related to exterior calculus [13, 11, 8, 3]. Application areas include hydrodynamics within fluid interfaces [28, 5], electrodynamics [32], stable methods for finite elements [4, 8, 3], and geometric processing in computer graphics [21, 11, 7, 34]. The exterior calculus of differential geometry provides a cooordinate invariant way to formulate equations on manifolds with close connections to topological and geometric structures inherent in mechanics [20, 16]. The exterior calculus provides less coordinate-centric expressions for analysis and numerical approximation that often can be interpreted more readily in terms of the geometry than alternative approaches such as the tensor calculus [10, 35, 16]. Many discrete exterior calculus approaches have been developed for efficient low-order approximations on triangulated meshes [13, 21, 11]. There has also been recent work on higher-order methods based on collocation for product bases [26], finite element differential forms [4], and more recently subdivision surfaces [7]. In these discrete exterior calculus approaches an effort is made to introduce approximations for fundamental operators such as the exterior derivative $\mathbf{d}$ and Hodge star $\star$ operators that on the discrete level preserve inherent geometric relations [11, 4, 3]. These operators are then used through composition to perform geometric processing tasks or approximate differential equations. We show here how related approaches can be developed using hyperinterpolation to obtain spectrally accurate methods for exterior calculus on radial manifolds.

We take an isogeometric approach based on the hyperinterpolation of spherical harmonics to represent the manifold geometry, differential forms, and related scalar and vector fields [15, 30, 6]. Hyperinterpolation methods use an oversampling of functions to grapple with some of the inherent challenges in designing optimal nodes for interpolation on general domains [30, 36, 29, 25]. This allows for the treatment of approximation instead using approaches such as $L^{2}$-orthogonal projection based on exact quadratures up to a desired order [30, 2]. To achieve discretizations on spherical toplogies with favorable symmetry, we use the nodes of Lebedev quadrature [17, 18]. While the more common approach of using samping points based on lattitude and longitude does provide some computational advantage through fast transforms, the sampling points have poor symmetry and exhibit a significant inhomogeneous distribution over the surface with many sample points clustering near the poles [27, 12, 9]. The Lebedev quadrature points are more regularly distributed on the surface and for a comperable number of points provide over the surface a more uniform resolution of functions. Furthermore, the Lebedev quadrature points are invariant under rotations corresponding to octohedral symmetry [17, 18]. We use an $L^{2}$-projection to spherical harmonics to approximate the exterior derivative operator $\mathbf{d}$ and Hodge star operator $\star$ on the surface. We show that our methods provide spectrally accurate approximations for these operators and their compositions. We show how the methods can be used to develop numerical solvers for differential equations on manifolds. We present results for the Laplace-Beltrami equations demonstrating our approach. The introduced numerical methods can be applied quite generally for approximating with spectral accuracy differential equations or other exterior calculus operations on radial manifolds.

2 Differential geometry and conservation laws on manifolds

We briefly discuss the formulation of conservation laws on manifolds in covariant form and discuss related concepts in differential geometry. More detailed discussions of the associated differential geometry can be found in [1, 31, 24]. Differential forms arise as a natural approach in formulating conservation laws and relations in continuum mechanics [10, 16]. The exterior calculus provides a convenient means to generalize the Stokes Theorem, Divergence Theorem, and Green’s Identities to curved spaces [20]. We introduce numerical methods for approximating exterior calculus operations on radial manifolds. A radial manifold is defined as a surface where every point can be connected by a line segment to the origin without intersecting the surface. We consider compact radial manifolds without boundaries which correspond to shapes with spherical topology, see Figure 1.

2.1 Covariant Form and Exterior Calculus

We briefly discuss the covariant formulation of tensors and exterior calculus. More detailed discussions can be found in [1, 31]. The covariant formulation relies on casting relationships using the cotangent bundle of the manifold. The cotangent space $\mathcal{V}^{*}$ at location $\mathbf{x}$ is a vector space that consists of all of the linear functionals that act on the tangent space $\mathcal{V}$. The dual vector $\mathbf{v}^{*}\in\mathcal{V}^{*}$ acts on a vector $\mathbf{u}\in\mathcal{V}$ as $\mathbf{v}^{*}[\mathbf{u}]=\langle\mathbf{v},\mathbf{u}\rangle=v^{i}g_{ij}u^{j}=v_{j}u^{j}$ where $g_{ij}$ denotes the entries of the metric tensor of the manifold [24, 1]. We use the lower index notational conventions for $v_{j}=g_{ij}v^{i}$ to obtain from the tangent vector components $v^{i}$ the associated covariant vector components $v_{j}$. The Einstein summation conventions are used throughout. We adopt the usual convention of indexing contravariant tensors with superscripts and covariant tensors with subscripts [1, 31]. Since each dual element $\mathbf{v}^{*}$ can be expressed through an inner-product with a representative tangent vector $\mathbf{v}$, the dual space $\mathcal{V}^{*}$ is isomorphic to the tangent space $\mathcal{V}$. We adopt the musical notation $\sharp$ and $\flat$ for raising and lowering indices. This corresponds to the isomorphisms between the tangent and co-tangent spaces of the surface given by

[TABLE]

The $\partial_{\mathbf{x}^{j}}$ denotes for a given choice of coordinates ${x}^{i}$ the associated basis vectors for the tangent space $\mathcal{V}$. The $d{\mathbf{x}^{j}}$ denotes for a given choice of coordinates the $1$-form basis elements for the co-tangent space $\mathcal{V}^{*}$ [1]. The isomorphisms can also be expressed directly in terms of the components as $v_{i}=g_{ij}v^{j}$ and $v^{i}=g^{ij}v_{j}$, where we denote the metric tensor as $g_{ij}$ and its inverse as $g^{ij}$. These conventions extend naturally to higher rank tensors [1].

We consider the exterior derivative $\mathbf{d}$ and the Hodge star $\star$ on the manifold. The exterior derivative $\mathbf{d}$ acts on a $k$-form $\alpha$ to yield the $k+1$-form

[TABLE]

Here we have taken that the $k$-form to be $\alpha=\left(1/k!\right)\alpha_{i_{1}\ldots i_{k}}\mathbf{d}\mathbf{x}^{i_{1}}\wedge\cdots\wedge\mathbf{d}\mathbf{x}^{i_{k}}$, where $\wedge$ denotes the wedge product [31, 1]. Since we are using the Einstein summation conventions, all permutations of the indices can arise and hence the factor of $1/k!$. For notational convenience and to make expressions more compact we will sometimes adsorb this factor implicity into $\alpha_{i_{1}\ldots i_{k}}$. The Hodge star $\star$ acts on a $k$-form and yields the $(n-k)$-form

[TABLE]

where $\alpha^{i_{1}\ldots i_{k}}=g^{i_{1}\ell_{1}}\cdots g^{i_{k}\ell_{k}}\alpha_{\ell_{1}\ldots\ell_{k}}$, $\sqrt{|g|}$ is the square-root of the determinant of the metric tensor, and $\epsilon_{i_{1}\ldots i_{k}j_{1}\ldots j_{n-k}}$ denotes the Levi-Civita tensor [20]. The Hodge star $\star$ can be thought of intuitively as yielding the orthogonal compliment to the $k$-form in the sense that the wedge product recovers a scaled volume form [13, 11]. This interpretation of the Hodge star is easiet to see when the local coordinates are chosen to be orthonormal at location $\mathbf{x}$ [31, 1].

Many of the operations in vector calculus can be generalized to manfolds using the exterior calculus as

[TABLE]

where $f$ is a smooth scalar function and the $\mathbf{F}$ is a smooth vector field. We use the notation $\boldsymbol{\delta}$ to denote the co-differential operator given by $\boldsymbol{\delta}=\star\mathbf{d}\star$. As can be seen above, the $\boldsymbol{\delta}$ operator is closely related to the divergence operation in vector calculus. It is also natural to consider common vector calculus differential operators such as the Laplacian. It is important to note that on manifolds there are a few different operators that share many of the features with the Laplacian of vector calculus. This requires care when formulating conservation laws or considering constitutive models. A few generalizations of the Laplacian include

[TABLE]

The $\Delta^{R}={\mbox{div}}({\mbox{grad}}(\cdot))$ denotes the rough-Laplacian given by the usual divergence of the gradient. For vector fields, $\Delta^{H}(\mathbf{F})$ denotes the Hodge-de Rham Laplacian, which has similarities to taking the curl of the curl [1]. In fact, in the case that ${\mbox{div}}(\mathbf{F})=-\boldsymbol{\delta}\mathbf{F}^{\flat}=0$ we have $\Delta^{H}(\mathbf{F})=\Delta^{S}(\mathbf{F})=-\left[\boldsymbol{\delta}\mathbf{d}\mathbf{F}^{\flat}\right]^{\sharp}$.

The exterior calculus provides generalizations of the Stokes Theorem and the Divergence Theorem to manifolds. These take the form respectively

[TABLE]

The $\omega$ is a $k$-form and $\Omega$ denotes a general smooth domain with smooth boundary $\partial\Omega$ within the manifold. When $\Omega$ is in two dimensions and $\omega$ is a $1$-form, the integrals on the left-hand-side perform a type of line integral over the boundary contour $\partial\Omega$. For the Stokes Theorem the component of the vector field that is tangent to the contour is integrated. For the Divergence Theorem the Hodge star $\star$ ensures that only the component of the vector field that is normal to the contour is integrated. The right-hand side then integrates over the interior region of $\Omega$. For the Stokes Theorem $\mathbf{d}\omega$ corresponds to a generalized representation of the curl operation. For the Divergence Theorem $\mathbf{d}\star\omega$ provides a representation of the divergence. Additional applications of the Hodge star $\star$ and isomorphisms $\flat$, $\sharp$ can also be useful to make conversions that bring these expressions into closer agreement with the standard vector calculus interpretations and intuition. Notice the relation of these expressions when $\omega=\mathbf{F}^{\flat}$ to equation 4. A useful feature of equation 6 is that it also generalizes readily to higher dimensions and $k$-forms.

2.2 Conservation Laws on Manifolds

For conservation laws, the exterior caclulus provides convenient ways to formulate relations and constitutive laws without relying on cumbersome coordinate expressions. One application of this idea can be seen by building on the generalized Stokes Theorem and the generalized Divergence Theorem in equation 6. For a conserved scalar quantity $u$, let $\omega$ denote a differential $1$-form that represents the local flux over a boundary, such as mass or energy. Conservation of $u$ and application of the generalized Divergence Theorem gives

[TABLE]

Since $u$ is a scalar field on the manifold ([math]-form), we represent its integral over the surface as the $2$-form $\star u=\tilde{u}_{ij}d\mathbf{x}^{i}\wedge d\mathbf{x}^{j}$. The flux is represented by the $1$-form $\omega=\omega_{i}d\mathbf{x}^{i}$. By the generalized Divergence Theorem, we have that the exterior derivative $\mathbf{d}$ gives the $2$-form $\mathbf{d}\star\omega=\tilde{\omega}_{ij}d\mathbf{x}^{i}\wedge d\mathbf{x}^{j}$ that integrates over $\Omega$ to give the same value as the flux $\star\omega$ integrated over the boundary $\partial\Omega$. Since $\Omega$ is arbitrary, this provides for smooth fields the local expression for the conservation law in covariant form

[TABLE]

To obtain this result, we applied the Hodge star $\star$ to both sides of equation 7 and used that $\star\star=-1^{m}$, where $m=k\cdot(n-k)$ for a $k$-form in $n$ dimensional space [1, 31]. The co-differential is given by $\boldsymbol{\delta}=\star\mathbf{d}\star$ and we see it plays here a role very similar to the vector calculus operation of divergence.

In many cases the flux itself follows a law depending on the conserved quantity $u$. In the case of Fourier’s Law for heat conduction or Fick’s Law for mass diffusion, the flux depends on the gradient of $u$. This generalizes on the manifold to

[TABLE]

The local conservation law in covariant form is

[TABLE]

When $u$ is heat or mass, this gives the generalization of the heat equation or the diffusion equation to the manifold.

We have given a brief discussion of the utility of exterior calculus methods in formulating scalar conservation laws on manifolds. Similar approaches can also be taken for vector conservation laws on manifolds and in the formulation of further constituitive laws and relations in continuum mechanics [28, 16]. A particular advantage of the exterior calculus is that we can formulate models and constuitive relations following closely intuition developed in the context of vector calculus without the need to deal with cumbersome coordinate expressions. We see that formulations rely on a composition of the exterior derivative $\mathbf{d}$ and Hodge star $\star$ operations. To utilize this approach in numerical calculations we need accurate approximations of the action of these operators on differential forms.

3 Numerical Methods for Exterior Calculus

We take an isogeometric approach based on the hyperinterpolation of spherical harmonics to represent the manifold geometry, differential forms, and related scalar and vector fields [15, 30, 6]. In hyperinterpolation, functions are oversampled to avoid many of the inherent issues associated with trying to design an optimal collection of nodes for Lagrange interpolation [30, 36, 29, 25]. This allows for functions to be approximated through $L^{2}$-orthogonal projections using exact quadratures up to a desired order [30, 2]. To achieve discretizations with favorable symmetry on the sphere, we use the nodes of Lebedev quadrature [17, 18]. This is in contrast to the more common approach of using samping points based on lattitude and longitude coordinates. While lattitude-longitude samplings have a computational advantage through fast transforms, the sampling points have poor symmetry and inhomogeneous distribution over the surface with many points clustering near the poles [27, 12, 9]. The Lebedev quadrature points provide a more regular distribution on the surface. For a comperable number of points, the Lebedev sampling provides a more uniform resolution of functions. The Lebedev quadrature points also have the feature of being invariant under rotations corresponding to octohedral symmetry [17, 18]. We show Lebedev nodes on example radial manifolds in Figure 3. For additional discussion of the differential geometry of radial manifolds see Appendix B.

3.1 Hyperinterpolation and $L^{2}$-Projection

We use hyperinterpolation to obtain a continuum approximation to fields on the manifold surface [30, 2]. To obtain a continuum representation of a function $f$ on the surface, we perform an $L^{2}$-orthogonal projection $\mathcal{P}$ to the space spanned by spherical harmonics up to order $\lfloor L/2\rfloor$,

[TABLE]

where $\hat{f}_{i}=\langle f,Y_{i}\rangle_{Q}$. We take the spherical harmonics $Y_{i}$ to be normalized with $\langle Y_{j},Y_{i}\rangle_{Q}=\delta_{ij}$ see Appendix A. We use the discrete inner-product defined by

[TABLE]

where $w_{\ell},\mathbf{x}_{\ell}$ are the Lebedev quadrature weights and nodes. When the quadrature is of order $L$ and the functions $u,v$ are each band-limited with respect to spherical harmonics up to order $\lfloor L/2\rfloor$, the the inner-product is the same as the $L^{2}$-inner-product $\langle u,v\rangle_{Q}=\langle u,v\rangle_{L^{2}}$. This yields the projection property $\mathcal{P}^{2}=\mathcal{P}$, see equation 11. For additional discussion of how we employ the spherical harmonics see Appendix A.

In practice, computing the inner-product $\langle\cdot,\cdot\rangle_{Q}$ only requires we know values of $f$ at the Lebedev nodes $\{\mathbf{x}_{\ell}\}$. This is utilized to represent functions on the surface in numerical calculations. We use this property to represent differential forms on the surface by an equivalent vector field at the Lebedev nodes. We perform in calculations conversions as needed by using the isomorphisms $\flat$, $\sharp$ in equations 1, see Figure 2. In two dimensions the [math]-forms and $2$-forms on the surface are equivalent to scalar fields (a vector field with only one component). For the more interesting case of $1$-forms on the surface, we use as our numerical represention an equivalent vector field with values specified at each of the Lebedev quadrature nodes. To simplify our discussion of our numerical methods we use the terminiology vector field and scalar field interchangably throughout.

The $1$-form $\mathbf{v}^{\flat}$ is equivalent through the isomorphism $\sharp$ to the vector field $\mathbf{v}^{\sharp}$. We represent $1$-forms by the values of $\mathbf{v}^{\sharp}$ stored at the Lebedev quadrature nodes $\{\mathbf{x}_{\ell}\}$. In numerical calculations we avoid issues with charts and coordinate singularities by representing the form as an expansion of spherical harmonics using the coordinates of the embedding space. This is done by representing the components of the associated vector field $\mathbf{v}^{\sharp}$ using the embedding space basis $\boldsymbol{\iota}_{1},\boldsymbol{\iota}_{2},\boldsymbol{\iota}_{3}$ as $\mathbf{v}^{\sharp}(\mathbf{x}_{\ell})=\bar{v}^{x}\boldsymbol{\iota}_{1}+\bar{v}^{y}\boldsymbol{\iota}_{2}+\bar{v}^{z}\boldsymbol{\iota}_{2}=\left[\bar{v}^{x},\bar{v}^{y},\bar{v}^{z}\right]_{\boldsymbol{\iota}_{1},\boldsymbol{\iota}_{2},\boldsymbol{\iota}_{3}}$. Storing the values $\left[\bar{v}^{x},\bar{v}^{y},\bar{v}^{z}\right]$ at the Lebedev nodes provides a convenient numerical representation of the differential form. To simplify the notation, we will often drop the subscript on $[\cdot,\cdot,\cdot]$ for the basis when it can be understood by context. When a continuum representation of the vector field is needed in our numerical calculations, we use the hyperinterpolation in equation 11 to obtain the associated smooth vector field

[TABLE]

We take a similar approach for [math]-forms and $2$-forms which are much easier to handle and are represented by the scalar field $\bar{v}$ to yield $\mathbf{\bar{v}}^{\sharp}(\theta,\phi)=\mathcal{P}\bar{v}$.

3.2 Numerical Exterior Calculus Operators $\overline{\mathbf{d}}$ and $\overline{\star}$

To approximate the exterior derivative $\mathbf{d}$, we need to approximate derivatives of our numerical representation at the Lebdev nodes for differential forms. For this purpose, we make use of the hyperinterpolation provided by equation 13. We remark that our approach in our numerical representation making use of the embedding space basis provides a global description of the differential form over the entire surface and a consistent way to obtain derivatives between different coordinate charts. For a given chart, a differential form has coordinate components for a [math]-form given by $\mathbf{v}^{\flat}=v$, a $1$-form by $\mathbf{v}^{\flat}=v_{i}d\mathbf{x}^{i}$, and a $2$-form by $\mathbf{v}^{\flat}=v_{ij}d\mathbf{x}^{i}\wedge d\mathbf{x}^{j}$. To numerically compute derivatives based on these expressions we perform a conversion from the vector field representation in the embedding space to the local coordinate representation of the components.

The $1$-form presents the most interesting case with the [math]-form and $2$-form handled similarly. For $1$-forms the components $v_{i}$ are related to the components $\bar{v}^{k}$ of the vector field representation by $v_{i}=g_{ij}v^{j}=g_{ij}a^{j}_{k}\bar{v}^{k}$. The term $a^{j}_{k}$ converts between the components $\bar{v}^{k}$ given in the coordinates of the embedding space $\boldsymbol{\iota}_{1},\boldsymbol{\iota}_{2},\boldsymbol{\iota}_{3}$ to the components ${v}^{j}$ given in the local coordinates on the surface $\partial_{\theta},\partial_{\phi}$. The exterior derivative of a $1$-form can be expressed in coordinate components as $\mathbf{d}\mathbf{v}^{\flat}=\partial_{s}v_{i}\hskip 2.84544ptd\mathbf{x}^{s}\wedge d\mathbf{x}^{i}$. For numerical calculations at a given location $\mathbf{x}_{\ell}$ (Lebedev node), we choose an appropriate coordinate chart that is locally non-degenerate. We compute the component derivative as

[TABLE]

The first two terms only depend on the geometry of the manifold and only the values of the differential form at location $\mathbf{x}_{\ell}$ (Lebedev node). This can be obtained readily from the spherical harmonics representation of the geometry of the manifold see Appendix A and Appendix B. In contrast, the last term depends on the derivatives of the coordinate components and requires use of the continuum representation from the hyperinterpolation obtained in equation 13. Putting this together with equation 14 and the coordinate expression for the exterior derivative, we obtain a numerical exterior derivative operator $\overline{\mathbf{d}}$ for $1$-forms. The case of [math]-form and $2$-form can be handled similarly. We should mention that the case of a $2$-form in two dimensions has exterior derivative zero which we also impose in our numerical calculations. In this manner, we obtain a numerical operator $\overline{\mathbf{d}}$ that maps a $k$-form defined at the Lebedev nodes to a $(k+1$)-form defined at the Lebedev nodes. This results in a convenient map between our representations useful in compositions for further application of other numerical exterior calculus operations.

We remark that in practical implimentations of the numerical exterior derivative operator $\overline{\mathbf{d}}$ it is convenient to represent the coordinate conversion in matrix-vector notation as $\mathbf{v}=GA^{-1}\mathbf{\bar{v}}$. The matrix entries $[A^{-1}]_{jk}=a^{j}_{k}$ correspond to the change from the coordinates of the embedding space to the local coordinates of the tangent space. With this notation, the derivatives in local coordinates can be expressed as ${\partial_{s}}\mathbf{v}=({\partial_{s}}G)A^{-1}\mathbf{\bar{v}}+G({\partial_{s}}A^{-1})\mathbf{\bar{v}}+GA^{-1}({\partial_{s}}\mathbf{\bar{v}})$. To avoid differentiating components of the inverse matrix $A^{-1}$, we use the identity ${\partial_{s}}A^{-1}=-A^{-1}({\partial_{s}}A)A^{-1}$ and use a linear algebra solver to compute the action of $A^{-1}$. This is done at each Lebedev node $\mathbf{x}_{\ell}$ with an appropriate choice made for the coordinate chart that is locally non-degenerate. For additional information on how we define the coordinate charts and how we perform practical computations from our representations of the manifold geometry see Appendix B.

We approximate next the Hodge star $\star$ operator on differential forms. We remark that in the related area of discrete exterior calculus (DEC) efforts are made to preserve geometric structure in the discrete setting, often on triangulated meshes. Interesting issues arise in DEC from the discrete geometry with which one must grapple and extensive studies have been conducted to formulate good approximations for the Hodge star $\star$ operator [23, 14]. Here we avoid many of these issues since we treat the operator at the continuum level and have more geometric information available to us from our spectral representation of both the manifold and the differential forms.

We approximate the Hodge star $\star$ operator on differential forms by a numerical operator $\overline{\star}$ which makes use of the representation at the Lebedev nodes. The Hodge star $\star$ has the feature that it is a local operation that involves values of the differential form and metric tensor only at an individual Lebedev node $\mathbf{x}_{\ell}$. We obtain a numerical operator $\overline{\star}$ by applying the isomorphisms and metric tensor using equation 1 and equation 3. The main consideration numerically is to choose well the coordinate chart so it is locally non-degenerate. The approximation enters through the fidelity of the metric tensor computed from our representation of the manifold geometry. The geometry for the radial manifold is determined by the radial function $r(\theta,\phi)$ which is represented as an expansion in a finite number of spherical harmonics up to order $L$, $r(\theta,\phi)=\sum_{i}\hat{r}_{i}Y_{i}(\theta,\phi)$. This provides in a given coordinate chart at location $\mathbf{x}_{\ell}$ the metric tensor and curvature tensor along with the local coordinate basis vectors $\partial_{\theta},\partial_{\phi}$ and their derivatives see Appendix B. The numerical Hodge star operator $\overline{\star}$ maps $k$-forms defined at the Lebedev nodes to $(n-k)$-forms defined at the Lebedev nodes. This provides a convenient map between representations for further application of numerical exterior calculus operations. In this manner, the numerical exterior derivative $\overline{\mathbf{d}}$ operator and numerical Hodge star $\overline{\star}$ operator can be used through composition to numerically approximate more complex exterior calculus operations on the manifold.

3.3 Solving Differential Equations on Manifolds

We consider differential equations on the manifold of the form

[TABLE]

The $\mathcal{L}$ denotes a linear differential operator that can be expressed in terms of a composition of the exterior calculus operations $\mathbf{d}$ and $\star$. For example, in the case of the Laplace-Beltrami equation this would correspond to the operator $\mathcal{L}=-\boldsymbol{\delta}\mathbf{d}=-\star\mathbf{d}\star\mathbf{d}$. We discretize the operator $\mathcal{L}$ by using a composition of the numerical exterior calculus operators $\overline{\star}$ and $\overline{\mathbf{d}}$ to obtain $\mathcal{\tilde{L}}$. For the Laplace-Beltrami equation, this corresponds to $\mathcal{\tilde{L}}=-\bar{\boldsymbol{\delta}}\mathbf{\bar{d}}=-\overline{\star}\hskip 2.84544pt\overline{\mathbf{d}}\hskip 2.84544pt\overline{\star}\hskip 2.84544pt\overline{\mathbf{d}}$. We approximate equation 15 on the manifold by

[TABLE]

The $\bar{u}=\sum_{j}\hat{u}_{j}Y_{j}$, $\bar{g}=\sum_{j}\hat{g}_{j}Y_{j}$ denote expansions up to order $\lfloor L/2\rfloor$ and $\langle\cdot,\cdot\rangle_{Q}$ denotes the Lebedev inner-product computed up to order $L$ in equation 12. This provides a Galerkin approximation where some additional sources of approximation arise from the treatment of the differential operator by $\mathcal{\tilde{L}}$. The approximation can be expressed in terms of the solution coefficients $\mathbf{\hat{u}}$ as

[TABLE]

The $\mathbf{\hat{u}}$ and $\mathbf{\hat{g}}$ denote the collection of coefficients $\hat{u}_{j}$ and $\hat{g}_{j}$ in the expansion of $\bar{u}$ and $\bar{g}$. The $K$ denotes the stiffness matrix with entries $K_{ij}=\langle\mathcal{\tilde{L}}Y_{j},Y_{i}\rangle_{Q}$ and $M$ denotes the mass matrix with entries $M_{ij}=\langle Y_{j},Y_{i}\rangle_{Q}$. We again mention that throughout the calculations we use expansions of functions up to order $\lfloor{L/2}\rfloor$ and Lebedev quadratures of order $L$. This is important to yield the needed over-sampling of functions to ensure accurate computation of the inner-product $\langle\cdot,\cdot\rangle_{Q}$.

4 Convergence Results

We discuss how our numerical methods converge in approximating the fundamental exterior calculus operations of the exterior derivative $\mathbf{d}$ and the Hodge star $\star$ when applied to [math]-forms, $1$-forms and $2$-forms. We then discuss the convergence of our methods for compositions of operators and present results for the Laplace-Beltrami equation.

We remark that throughout our convergence studies, we describe test functions using for a point $\mathbf{x}$ on the manifold its location within the embedding space. To perform calculations we make use of the embedding space coordinates $[x,y,z]$ corresponding to $\mathbf{x}=x\boldsymbol{\iota}_{1}+y\boldsymbol{\iota}_{2}+z\boldsymbol{\iota}_{3}$, where $\boldsymbol{\iota}_{1},\boldsymbol{\iota}_{2},\boldsymbol{\iota}_{3}$ is the basis for the embedding space. In this manner our test data is not tied to a specific choice of local coordinates on the manifold. All figures report the relative error $\epsilon_{\mbox{\tiny rel}}=\|\mathbf{\bar{w}}^{\sharp}-\mathbf{w}^{\sharp}\|_{2}/\|\mathbf{w}^{\sharp}\|_{2}$. The $\mathbf{w}$ is the exact result and $\mathbf{\bar{w}}$ is the numerically computed result. For $\|\cdot\|_{2}$, we use the $L^{2}$-norm of the embedding space.

4.1 Convergence of the Hodge Star Operator

We approximate the Hodge star $\star$ by the numerical operator $\overline{\star}$ using the hyperinterpolation approach discussed in Section 3.2. We investigate in practice the accuracy of this approach on a few different geometries and differential forms.

We first consider a [math]-form defined on Manifold B defined in Figure 1. We take the [math]-form to be $f=\exp(z)/(3-y)$. We show the accuracy of our numerical operator $\overline{\star}$ in approximating the Hodge star $\star$ as the number of Lebedev nodes increases. We find that the main limitation in the accuracy of the $\overline{\star}$ is the resolution of the geometry of the manifold. This is seen in our results where once a sufficient number of spherical harmonic modes are reached the relative error rapidly decays in approximating $\star f$. The convergence of $\overline{\star}$ as the number of Lebedev nodes is increased and when the geometry is varied is shown in Figure 4.

We next consider on Manifold B the $1$-form $\alpha=\sqrt{|g|}\exp(z)d\theta+\sqrt{|g|}\exp(z)d\phi$. We again find that the main limitation in the accuracy of the $\overline{\star}$ is the resolution of the geometry of the manifold. In this case we find the error rapidly decreases once a sufficient number spherical harmonic modes are used. The convergence of $\overline{\star}$ as the number of Lebedev nodes is increased and when the geometry is varied is shown in Figure 5.

These results indicate that the numerical operator $\overline{\star}$ provides an accurate approximation to the Hodge star $\star$. Given the localized nature of the Hodge star, the main consideration to obtain accurate results with $\overline{\star}$ is to use enough Lebedev nodes to ensure sufficient resolution of the geometry of the manifold.

4.2 Convergence of the Exterior Derivative Operator

We investigate the convergence of the numerical operator $\overline{\mathbf{d}}$ in approximating the exterior derivative $\mathbf{d}$ when applied to [math]-forms and $1$-forms. We do not consider $2$-forms here since in two dimensions the exterior derivative would be zero which we also impose in our numerical calculations [1, 31]. We consider the [math]-form given by $f=\exp(z)$ and the $1$-form given by $\alpha=|g|\exp(z)d\theta+|g|\exp(z)d\phi$. We consider for Manifold B the relative error of $\overline{\mathbf{d}}$ in approximating $\mathbf{d}$ as the number of Lebedev nodes is increased and the manifold geometry is varied. These results are shown in Figure 6 and Figure 7.

We see that the numerical operator $\overline{\mathbf{d}}$ converges spectrally in approximating the exterior derivative $\mathbf{d}$ both for the [math]-forms and $1$-forms. Interestingly, we see that in the case when $r_{0}=0.0$ the approximation converges significantly more rapid than the cases when $r_{0}\neq 0$. This occurs since $r_{0}=0$ corresponds to the case when the shape is a sphere where the geometry is relatively simple and many of the geometric terms to be numerically approximated greatly simplify. The main source of error in this case arises primarily from the hyperinterpolation used for computing the derivatives. In the case when $r_{0}\neq 0$, the isogeometric approach used to compute $\overline{\mathbf{d}}$ approximates the geometry of the manifold using a finite spherical harmonics representation which results in an additional source of approximation error.

We investigate how the geometry contributes to convergence by varying $r_{0}$ over the range $[0.0,0.4]$. The case with $r_{0}=0.4$ corresponds to the final shape of Manifold B shown in Figure 1. The convergence is found to be spectral and comparable in each of these different cases for the geometry. We see that as one might expect the largest errors are incurred in the case with the most pronounced geometry corresponding to $r_{0}=0.4$. Overall, the results indicate that the numerical operator $\overline{\mathbf{d}}$ provides for [math]-forms and $1$-forms an accurate approximation for the exterior derivative $\mathbf{d}$.

4.3 Convergence for the Laplace-Beltrami Equation

We consider the composition of exterior calculus operators to represent partial differential equations on manifolds. Taking the approach we discussed in Section 3.3, we consider the Laplace-Beltrami operator $\mathcal{L}=-\boldsymbol{\delta}\mathbf{d}=-\star\mathbf{d}\star\mathbf{d}$ and consider the Poisson problem on the manifold

[TABLE]

Since there are no boundaries for the manifolds we shall consider, we also impose throughout that the zero mode has $\hat{u}_{0}=0$. The Laplace-Beltrami equation 18 provides a test of using compositions of the numerical exterior calculus operators $\overline{\mathbf{d}}$ and $\overline{\star}$ and the associated convergence and accuracy.

We investigate how our methods solve the Laplace-Beltrami equation 18 on Manifold B and Manifold C of Figure 1. To have a known solution to the Laplace-Beltrami equations with which to compare, we manifacture a source term for the solution on the manifold given by $u=\exp(y)/(3-z)^{4}$, where the $x,y,z$ refer to the coordinates of the embedding space as discussed in the beginning of Section 4. We compute the source term $g=-\mathcal{L}u$ symbolically using the package SymPy and evaluate the expressions numerically when data is needed for the source [22]. In this manner we are able to assess the relative errors of the numerical methods with high-precision.

We show the source term $g$, solution $u$, and relative errors $\epsilon_{rel}$ in Figure 8 and Figure 9. We use the approach discussed for solving PDEs on manifolds in Section 3.3 for the Laplace-Beltrami equations. We show how the composition of the numerical operators $\overline{\mathbf{d}}$ and $\overline{\star}$ perform as the number of Lebedev nodes increases in Figure 10. The results indicate that the numerical methods based on composing $\overline{\mathbf{d}}$ and $\overline{\star}$ perform well in approximating the true compositions of the exterior calculus operators $\star$ and $\mathbf{d}$.

We see that when the geometry has $r_{0}=0.0$ corresponding to a sphere the methods converge most rapidly. In the cases with $r_{0}\neq 0$ we see the convergence is comparable in each case. We see the convergence is slowest for the most pronounced geometry where the representation for the geometry incurs the most approximation error. Overall, we find the compositions of $\overline{\mathbf{d}}$ and $\overline{\star}$ perform well and provide a solver with spectral accuracy for the Laplace-Beltrami equations on Manifold B.

We next consider the more geometrically complex Manifold C, see Figure 1. We again use the manufactured solution on the surface given by $u=\exp(y)/(3-z)^{4}$ and $g=-\mathcal{L}u$ computed symbolically. We mention that given the different geometry of Manifold C relative to Manifold B both the solution and the source data are different on each of the manifold surfaces. We show on the manifold surface the source term $g$, solution $u$, and relative error $\epsilon_{rel}$ in Figure 11 and Figure 12.

For the more geometrically complicated Manifold C, we perform calculations to show how the composition of the numerical operators $\overline{\mathbf{d}}$ and $\overline{\star}$ perform as the number of Lebedev nodes increases and $r_{0}$ is varied, see Figure 13.

We find the numerical methods converge spectrally for Manifold C. As may be expected given its more complicated geometry, we find the convergence is somewhat slower than for Manifold B. We find that the numerical methods perform well for Manifold C once a sufficient number of Lebedev nodes are used to resolve the geometry. We again see the trend that the case with smaller $r_{0}$ convergence more rapidly. Overall, we find that in practice the composition of the numerical operator $\overline{\mathbf{d}}$ and $\overline{\star}$ performs well in approximating differential operators on the manifold. This indicates the approach discussed in Section 3.3 can be applied quite generally in practice for solving partial differential equations on manifolds.

5 Conclusion

In summary, we have developed numerical methods for exterior calculus on radial manifolds based on hyperinterpolation and Lebedev quadrature. We have shown our methods are spectrally accurate in approximating the exterior derivative $\mathbf{d}$, Hodge star $\star$, and their compositions. We presented results for the Laplace-Beltrami equations that demonstrate the utility of the approach. The introduced numerical methods can be applied quite generally for approximating exterior calculus operations and differential equations on radial manifolds.

6 Acknowledgments

We acknowledge support to P.J.A. and B.G. from research grants DOE ASCR CM4 DE-SC0009254, NSF CAREER Grant DMS-0956210, and NSF Grant DMS - 1616353.

Appendix

Appendix A Spherical Harmonics

The spherical harmonics are given by

[TABLE]

where $m$ denotes the order and $n$ the degree for $n\geq 0$ and $m\in\{-n,\dots,n\}$. The $P^{m}_{n}$ denote the Associated Legendre Polynomials. In our notation, $\theta$ denotes the azmuthal angle and $\phi$ the polar angle of the spherical coordinates [6].

Since we work throughout only with real-valued functions, we have that the modes are self-conjugate and we use that $Y^{m}_{n}=\overline{Y^{-m}_{n}}$. We have found it convenient to represent the spherical harmonics as

[TABLE]

where $X_{n}^{m}$ and $Z_{n}^{m}$ denote the real and imaginary parts. In our numerical methods we use this splitting to construct a purely real set of basis functions on the unit sphere with maximum degree $N$. We remark that this consists of $(N+1)^{2}$ basis elements. In the case of $N=2$ we have the basis elements

[TABLE]

We use a similar convention for the basis for the other values of $N$. We take our final basis elements $Y_{i}$ to be the normalized as $Y_{i}=\tilde{Y}_{i}/\sqrt{\langle\tilde{Y}_{i},\tilde{Y}_{i}\rangle}$.

We compute derivatives of our finite expansions by evaluating analytic formulas for the spherical harmonics in order to try to minimize approximation error [6]. Approximation errors are incurred when sampling the values of expressions involving these derivatives at the Lebedev nodes and performing quadratures. The derivative in the azmuthal coordinate $\theta$ of the spherical harmonics is given by

[TABLE]

This maps the spherical harmonic of degree $n$ to again a spherical harmonic of degree $n$. In our numerics, this derivative can be represented in our finite basis which allows us to avoid projections. This allows for computing the derivative in $\theta$ without incurring an approximation error. For the derivative in the polar angle $\phi$ we have that

[TABLE]

We remark that the expression can not be represented in terms of a finite expansion of spherical harmonics. We use this expression for $\partial_{\phi}Y^{m}_{n}(\theta,\phi)$ when we compute values at the Lebedev quadrature nodes in equation 14. This provides a convenient way to compute derivatives of differential forms following the approach discussed in Section 3. We remark that it is the subsequent hyperinterpolation of the resulting expressions where the approximation error is incurred. We adopt the notational convention that $Y^{m}_{n}=0$ when $m\geq n+1$. For further discussion of spherical harmonics see [6].

Appendix B Differential Geometry of Radial Manifolds

We consider throughout manifolds of radial shape. A radial manifold is defined as a surface where each point can be connected by a line segment to the origin without intersecting the surface. In spherical coordinates, any point $\mathbf{x}$ on the radial manifold can be expressed as

[TABLE]

where $\mathbf{r}$ is the unit vector from the origin to the point on the sphere corresponding to angle $\theta,\phi$ and $r$ is a positive scalar function.

We take an isogeometric approach to representing the manifold $M$. We sample the scalar function $r$ at the Lebedev nodes and represent the geometry using the finite spherical harmonics expansion $r(\theta,\phi)=\sum_{i}\bar{r}_{i}Y_{i}$ up to the order $\lfloor L/2\rfloor$ where $\bar{r}_{i}=\langle r,Y_{i}\rangle_{Q}$ for a quadrature of order $L$.

We consider two coordinate charts for our calculations. The first is referred to as Chart A and has coordinate singularities at the north and south pole. The second is referred to as Chart B and has coordinate singularities at the east and west pole. For each chart we use spherical coordinates with $(\theta,\phi)\in[0,2\pi)\times[0,\pi]$ but to avoid singularities only use values in the restricted range $\phi\in[\phi_{min},\phi_{max}]$, where $0<\phi_{min}\leq\frac{\pi}{4}$, and $\frac{3\pi}{4}\leq\phi_{max}<\pi$. In practice, one typically takes $\phi_{min}=0.8\times\frac{\pi}{4}$ and $\phi_{max}=0.8\times\pi$. For Chart A, the manifold is parameterized in the embedding space $\mathbb{R}^{3}$ as

[TABLE]

and for Chart B

[TABLE]

With these coordinate representations, we can derive explicit expressions for geometric quantities associated with the manifold such as the metric tensor and shape tensor. The derivatives used as the basis $\partial_{\theta},\partial_{\phi}$ for the tangent space can be expressed as

[TABLE]

We have expressions for $\mathbf{r}_{\theta}$ and $\mathbf{r}_{\phi}$ in the embedding space $\mathbb{R}^{3}$ using equation 24 or equation 25 depending on the chart being used. The first fundamental form $\mathbf{I}$ (metric tensor) and second fundamental form $\mathbf{II}$ (shape tensor) are given by

[TABLE]

and

[TABLE]

The $\mathbf{n}$ denotes the outward normal on the surface and is computed using

[TABLE]

The terms $\boldsymbol{\sigma}_{\theta\theta}$, $\boldsymbol{\sigma}_{\theta\phi}$, and $\boldsymbol{\sigma}_{\phi,\phi}$ are obtained by further differentiation from equation 26 and equation 27. We use the notation for the metric tensor $\mathbf{g}=\mathbf{I}$ interchangably. In practical calculations whenever we need to compute the action of the inverse metric tensor we do so through numerical linear algebra (Gaussian elimination with pivoting) [19, 33]. For notational convenience, we use the tensor notation for the metric tensor $g_{ij}$ and its inverse $g^{ij}$ which has the formal correspondence

[TABLE]

For the metric factor we also have that

[TABLE]

To ensure accurate numerical calculations in each of the above expressions the appropriate coordinates either Chart A or Chart B are used to ensure sufficient distance from coordinate singularities at the poles. To compute quantities associated with curvature of the manifold we use the Weingarten map [24] which can be expressed as

[TABLE]

To compute the Gaussian curvature $K$, we use

[TABLE]

For further discussion of the differential geometry of manifolds see [24, 1, 31].