Spectral Numerical Exterior Calculus Methods for Differential Equations on Radial Manifolds
SSpectral Numerical Exterior Calculus Methods forDifferential Equations on Radial Manifolds
B. Gross and P. J. Atzberger ∗ Department of Mathematics, Department of Mechanical Engineering,University of California Santa Barbara; [email protected]; http://atzberger.org/
March 14, 2017 W e develop exterior calculus approaches for partial differential equations on ra-dial manifolds. We introduce numerical methods that approximate with spec-tral accuracy the exterior derivative d, Hodge star (cid:63) , and their compositions.To achieve discretizations with high precision and symmetry, we develop hyperinter-polation methods based on spherical harmonics and Lebedev quadrature. We performconvergence studies of our numerical exterior derivative operator d and Hodge staroperator (cid:63) showing each converge spectrally to d and (cid:63) . We show how the numericaloperators can be naturally composed to formulate general numerical approximationsfor 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 [3, 8, 11, 13].Application areas include hydrodynamics within fluid interfaces [5, 28], electrodynamics [32], stablemethods for finite elements [3, 4, 8], and geometric processing in computer graphics [7, 11, 21, 34].The exterior calculus of differential geometry provides a cooordinate invariant way to formulateequations on manifolds with close connections to topological and geometric structures inherent inmechanics [16, 20]. The exterior calculus provides less coordinate-centric expressions for analysisand numerical approximation that often can be interpreted more readily in terms of the geometrythan alternative approaches such as the tensor calculus [10, 16, 35]. Many discrete exterior calculusapproaches have been developed for efficient low-order approximations on triangulated meshes [11,13, 21]. There has also been recent work on higher-order methods based on collocation for productbases [26], finite element differential forms [4], and more recently subdivision surfaces [7]. In thesediscrete exterior calculus approaches an effort is made to introduce approximations for fundamentaloperators such as the exterior derivative d and Hodge star (cid:63) operators that on the discrete levelpreserve inherent geometric relations [3, 4, 11]. These operators are then used through compositionto perform geometric processing tasks or approximate differential equations. We show here howrelated approaches can be developed using hyperinterpolation to obtain spectrally accurate methodsfor exterior calculus on radial manifolds. * Work supported by DOE ASCR CM4 DE-SC0009254, NSF CAREER Grant DMS-0956210, andNSF Grant DMS - 1616353. Page 1 of 24 a r X i v : . [ m a t h . NA ] M a r igure 1: Radial Manifolds: A radial manifold is defined as a surface where each point can be connectedby a line segment to the origin without intersecting the surface. Shown are three radial manifolds whichfor discussions we refer to interchangeably as the (i) Sphere / Manifold A, (ii) Dimple / Manifold B, and(iii) Fountain / Manifold C. The manifolds are generated by the radial functions (i) r ( θ, φ ) = 1 . , (ii) r ( θ, φ ) = 1 + r sin(3 φ ) cos( θ ) with r = 0 . , and (iii) r ( θ, φ ) = 1 + r sin(7 φ ) cos( θ ) with r = 0 . . Foradditional discussion of the differential geometry of radial manifolds see Appendix B. We take an isogeometric approach based on the hyperinterpolation of spherical harmonics torepresent the manifold geometry, differential forms, and related scalar and vector fields [6, 15, 30].Hyperinterpolation methods use an oversampling of functions to grapple with some of the inherentchallenges in designing optimal nodes for interpolation on general domains [25, 29, 30, 36]. Thisallows for the treatment of approximation instead using approaches such as L -orthogonal projectionbased on exact quadratures up to a desired order [2, 30]. To achieve discretizations on sphericaltoplogies with favorable symmetry, we use the nodes of Lebedev quadrature [17, 18]. While themore common approach of using samping points based on lattitude and longitude does providesome computational advantage through fast transforms, the sampling points have poor symmetryand exhibit a significant inhomogeneous distribution over the surface with many sample pointsclustering near the poles [9, 12, 27]. The Lebedev quadrature points are more regularly distributedon the surface and for a comperable number of points provide over the surface a more uniformresolution of functions. Furthermore, the Lebedev quadrature points are invariant under rotationscorresponding to octohedral symmetry [17, 18]. We use an L -projection to spherical harmonicsto approximate the exterior derivative operator d and Hodge star operator (cid:63) on the surface. Weshow that our methods provide spectrally accurate approximations for these operators and theircompositions. We show how the methods can be used to develop numerical solvers for differentialequations on manifolds. We present results for the Laplace-Beltrami equations demonstrating ourapproach. The introduced numerical methods can be applied quite generally for approximating withspectral 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 discussrelated concepts in differential geometry. More detailed discussions of the associated differentialgeometry can be found in [1, 24, 31]. Differential forms arise as a natural approach in formulatingconservation laws and relations in continuum mechanics [10, 16]. The exterior calculus provides a
Page 2 of 24 onvenient means to generalize the Stokes Theorem, Divergence Theorem, and Green’s Identities tocurved spaces [20]. We introduce numerical methods for approximating exterior calculus operationson radial manifolds. A radial manifold is defined as a surface where every point can be connected bya line segment to the origin without intersecting the surface. We consider compact radial manifoldswithout boundaries which correspond to shapes with spherical topology, see Figure 1.
Figure 2:
We consider vectors and covectors on the manifold. The isomorphisms between the tangentspace and co-tangent space are given by (cid:91) [ v ] = v j ∂ x j → v i d x i = v ∗ and (cid:93) [ v ∗ ] = v i d x i → v j ∂ x j = v . Themusical notation (cid:93) and (cid:91) is used to indicate how the indices are either being raised or lowered in the tensornotation [1, 31]. We briefly discuss the covariant formulation of tensors and exterior calculus. More detailed discussionscan be found in [1, 31]. The covariant formulation relies on casting relationships using the cotangentbundle of the manifold. The cotangent space V ∗ at location x is a vector space that consists ofall of the linear functionals that act on the tangent space V . The dual vector v ∗ ∈ V ∗ acts ona vector u ∈ V as v ∗ [ u ] = (cid:104) v , u (cid:105) = v i g ij u j = v j u j where g ij denotes the entries of the metrictensor of the manifold [1, 24]. We use the lower index notational conventions for v j = g ij v i toobtain from the tangent vector components v i the associated covariant vector components v j . TheEinstein summation conventions are used throughout. We adopt the usual convention of indexingcontravariant tensors with superscripts and covariant tensors with subscripts [1, 31]. Since eachdual element v ∗ can be expressed through an inner-product with a representative tangent vector v ,the dual space V ∗ is isomorphic to the tangent space V . We adopt the musical notation (cid:93) and (cid:91) for raising and lowering indices. This corresponds to the isomorphisms between the tangent andco-tangent spaces of the surface given by (cid:91) : v j ∂ x j → v i d x i , (cid:93) : v i d x i → v j ∂ x j . (1)The ∂ x j denotes for a given choice of coordinates x i the associated basis vectors for the tangent space V . The d x j denotes for a given choice of coordinates the 1-form basis elements for the co-tangentspace 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 . Theseconventions extend naturally to higher rank tensors [1].We consider the exterior derivative d and the Hodge star (cid:63) on the manifold. The exterior Page 3 of 24 erivative d acts on a k -form α to yield the k + 1-form d α = 1 k ! ∂α i ...i k ∂x j dx j ∧ dx i ∧ · · · ∧ dx i k . (2)Here we have taken that the k -form to be α = (1 /k !) α i ...i k dx i ∧ · · · ∧ dx i k , where ∧ denotes thewedge product [1, 31]. Since we are using the Einstein summation conventions, all permutationsof the indices can arise and hence the factor of 1 /k !. For notational convenience and to makeexpressions more compact we will sometimes adsorb this factor implicity into α i ...i k . The Hodgestar (cid:63) acts on a k -form and yields the ( n − k )-form (cid:63)α = (cid:112) | g | ( n − k )! k ! α i ...i k (cid:15) i ...i k j ...j n − k · dx j ∧ · · · ∧ dx j n − k , (3)where α i ...i k = g i (cid:96) · · · g i k (cid:96) k α (cid:96) ...(cid:96) k , (cid:112) | g | is the square-root of the determinant of the metric tensor,and (cid:15) i ...i k j ...j n − k denotes the Levi-Civita tensor [20]. The Hodge star (cid:63) can be thought of intuitivelyas yielding the orthogonal compliment to the k -form in the sense that the wedge product recovers ascaled volume form [11, 13]. This interpretation of the Hodge star is easiet to see when the localcoordinates are chosen to be orthonormal at location x [1, 31].Many of the operations in vector calculus can be generalized to manfolds using the exteriorcalculus asgrad( f ) = [ d f ] (cid:93) , div( F ) = − ( (cid:63) d (cid:63) F (cid:91) ) = − δ F (cid:91) , curl( F ) = (cid:104) (cid:63) ( dF (cid:91) ) (cid:105) (cid:93) , (4)where f is a smooth scalar function and the F is a smooth vector field. We use the notation δ todenote the co-differential operator given by δ = (cid:63) d (cid:63) . As can be seen above, the δ operator is closelyrelated to the divergence operation in vector calculus. It is also natural to consider common vectorcalculus differential operators such as the Laplacian. It is important to note that on manifolds thereare 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 fewgeneralizations of the Laplacian include∆ H ( F ) = − (cid:104) ( δ d + d δ ) F (cid:91) (cid:105) (cid:93) , ∆ S ( F ) = − (cid:104) δ dF (cid:91) (cid:105) (cid:93) , ∆ H f = ∆ R f = − ( (cid:63) d (cid:63) ) d f = − δ d f. (5)The ∆ R = div(grad( · )) denotes the rough-Laplacian given by the usual divergence of the gradient. Forvector fields, ∆ H ( F ) denotes the Hodge-de Rham Laplacian, which has similarities to taking the curlof the curl [1]. In fact, in the case that div( F ) = − δ F (cid:91) = 0 we have ∆ H ( F ) = ∆ S ( F ) = − (cid:2) δ dF (cid:91) (cid:3) (cid:93) .The exterior calculus provides generalizations of the Stokes Theorem and the DivergenceTheorem to manifolds. These take the form respectively (cid:90) ∂ Ω ω = (cid:90) Ω d ω, (cid:90) ∂ Ω (cid:63)ω = (cid:90) Ω d (cid:63) ω. (6)The ω is a k -form and Ω denotes a general smooth domain with smooth boundary ∂ Ω within themanifold. When Ω is in two dimensions and ω is a 1-form, the integrals on the left-hand-side performa type of line integral over the boundary contour ∂ Ω. For the Stokes Theorem the component
Page 4 of 24 f the vector field that is tangent to the contour is integrated. For the Divergence Theorem theHodge star (cid:63) ensures that only the component of the vector field that is normal to the contouris integrated. The right-hand side then integrates over the interior region of Ω. For the StokesTheorem d ω corresponds to a generalized representation of the curl operation. For the DivergenceTheorem d (cid:63) ω provides a representation of the divergence. Additional applications of the Hodgestar (cid:63) and isomorphisms (cid:91) , (cid:93) can also be useful to make conversions that bring these expressionsinto closer agreement with the standard vector calculus interpretations and intuition. Notice therelation of these expressions when ω = F (cid:91) to equation 4. A useful feature of equation 6 is that italso generalizes readily to higher dimensions and k -forms. For conservation laws, the exterior caclulus provides convenient ways to formulate relations andconstitutive laws without relying on cumbersome coordinate expressions. One application of thisidea can be seen by building on the generalized Stokes Theorem and the generalized DivergenceTheorem in equation 6. For a conserved scalar quantity u , let ω denote a differential 1-form thatrepresents the local flux over a boundary, such as mass or energy. Conservation of u and applicationof the generalized Divergence Theorem gives ∂∂t (cid:90) Ω (cid:63)u = (cid:90) ∂ Ω (cid:63)ω = (cid:90) Ω d (cid:63) ω. (7)Since u is a scalar field on the manifold (0-form), we represent its integral over the surface as the2-form (cid:63)u = ˜ u ij d x i ∧ d x j . The flux is represented by the 1-form ω = ω i d x i . By the generalizedDivergence Theorem, we have that the exterior derivative d gives the 2-form d (cid:63) ω = ˜ ω ij d x i ∧ d x j that integrates over Ω to give the same value as the flux (cid:63)ω integrated over the boundary ∂ Ω.Since Ω is arbitrary, this provides for smooth fields the local expression for the conservation law incovariant form ∂u∂t = − (cid:63) d (cid:63) ω = − δ ω. (8)To obtain this result, we applied the Hodge star (cid:63) to both sides of equation 7 and used that (cid:63)(cid:63) = − m , where m = k · ( n − k ) for a k -form in n dimensional space [1, 31]. The co-differential isgiven by δ = (cid:63) d (cid:63) and we see it plays here a role very similar to the vector calculus operation ofdivergence.In many cases the flux itself follows a law depending on the conserved quantity u . In the caseof Fourier’s Law for heat conduction or Fick’s Law for mass diffusion, the flux depends on thegradient of u . This generalizes on the manifold to ω = d u. (9)The local conservation law in covariant form is ∂u∂t = − δ d u. (10)When u is heat or mass, this gives the generalization of the heat equation or the diffusion equationto the manifold. Page 5 of 24 e have given a brief discussion of the utility of exterior calculus methods in formulatingscalar conservation laws on manifolds. Similar approaches can also be taken for vector conservationlaws on manifolds and in the formulation of further constituitive laws and relations in continuummechanics [16, 28]. A particular advantage of the exterior calculus is that we can formulate modelsand constuitive relations following closely intuition developed in the context of vector calculuswithout the need to deal with cumbersome coordinate expressions. We see that formulations rely ona composition of the exterior derivative d and Hodge star (cid:63) operations. To utilize this approachin numerical calculations we need accurate approximations of the action of these operators ondifferential forms.
3. Numerical Methods for Exterior Calculus
We take an isogeometric approach based on the hyperinterpolation of spherical harmonics torepresent the manifold geometry, differential forms, and related scalar and vector fields [6, 15, 30]. Inhyperinterpolation, functions are oversampled to avoid many of the inherent issues associated withtrying to design an optimal collection of nodes for Lagrange interpolation [25, 29, 30, 36]. This allowsfor functions to be approximated through L -orthogonal projections using exact quadratures up toa desired order [2, 30]. To achieve discretizations with favorable symmetry on the sphere, we use thenodes of Lebedev quadrature [17, 18]. This is in contrast to the more common approach of usingsamping points based on lattitude and longitude coordinates. While lattitude-longitude samplingshave a computational advantage through fast transforms, the sampling points have poor symmetryand inhomogeneous distribution over the surface with many points clustering near the poles [9,12,27].The Lebedev quadrature points provide a more regular distribution on the surface. For a comperablenumber of points, the Lebedev sampling provides a more uniform resolution of functions. TheLebedev quadrature points also have the feature of being invariant under rotations corresponding tooctohedral symmetry [17, 18]. We show Lebedev nodes on example radial manifolds in Figure 3. Foradditional discussion of the differential geometry of radial manifolds see Appendix B. Figure 3:
Lebedev Quadrature. Shown are the sample points of the Lebedev quadrature in the case of 302nodes. The Lebedev nodes distribute nearly uniformly over the surface and are invariant under the rotationscorresponding to octahedral symmetry [17, 18].
Page 6 of 24 .1. Hyperinterpolation and L -Projection We use hyperinterpolation to obtain a continuum approximation to fields on the manifold surface [2,30]. To obtain a continuum representation of a function f on the surface, we perform an L -orthogonalprojection P to the space spanned by spherical harmonics up to order (cid:98) L/ (cid:99) , P [ f ] = ¯ f ( θ, φ ) = (cid:88) i ˆ f i Y i ( θ, φ ) , (11)where ˆ f i = (cid:104) f, Y i (cid:105) Q . We take the spherical harmonics Y i to be normalized with (cid:104) Y j , Y i (cid:105) Q = δ ij seeAppendix A. We use the discrete inner-product defined by (cid:104) u, v (cid:105) Q = (cid:88) (cid:96) w (cid:96) u ( x (cid:96) ) v ( x (cid:96) ) . (12)where w (cid:96) , x (cid:96) 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 (cid:98) L/ (cid:99) , the the inner-product is the same as the L -inner-product (cid:104) u, v (cid:105) Q = (cid:104) u, v (cid:105) L . This yieldsthe projection property P = P , see equation 11. For additional discussion of how we employ thespherical harmonics see Appendix A.In practice, computing the inner-product (cid:104)· , ·(cid:105) Q only requires we know values of f at theLebedev nodes { x (cid:96) } . 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 atthe Lebedev nodes. We perform in calculations conversions as needed by using the isomorphisms (cid:91) , (cid:93) in equations 1, see Figure 2. In two dimensions the 0-forms and 2-forms on the surface areequivalent to scalar fields (a vector field with only one component). For the more interesting case of1-forms on the surface, we use as our numerical represention an equivalent vector field with valuesspecified at each of the Lebedev quadrature nodes. To simplify our discussion of our numericalmethods we use the terminiology vector field and scalar field interchangably throughout.The 1-form v (cid:91) is equivalent through the isomorphism (cid:93) to the vector field v (cid:93) . We represent1-forms by the values of v (cid:93) stored at the Lebedev quadrature nodes { x (cid:96) } . In numerical calculationswe avoid issues with charts and coordinate singularities by representing the form as an expansion ofspherical harmonics using the coordinates of the embedding space. This is done by representingthe components of the associated vector field v (cid:93) using the embedding space basis ι , ι , ι as v (cid:93) ( x (cid:96) ) = ¯ v x ι + ¯ v y ι + ¯ v z ι = [¯ v x , ¯ v y , ¯ v z ] ι , ι , ι . Storing the values [¯ v x , ¯ v y , ¯ v z ] at the Lebedev nodesprovides a convenient numerical representation of the differential form. To simplify the notation, wewill often drop the subscript on [ · , · , · ] for the basis when it can be understood by context. Whena continuum representation of the vector field is needed in our numerical calculations, we use thehyperinterpolation in equation 11 to obtain the associated smooth vector field ¯v (cid:93) ( θ, φ ) = [ P ¯ v x , P ¯ v y , P ¯ v z ] . (13)We take a similar approach for 0-forms and 2-forms which are much easier to handle and arerepresented by the scalar field ¯ v to yield ¯v (cid:93) ( θ, φ ) = P ¯ v . Page 7 of 24 .2. Numerical Exterior Calculus Operators d and (cid:63) To approximate the exterior derivative d , we need to approximate derivatives of our numericalrepresentation at the Lebdev nodes for differential forms. For this purpose, we make use of thehyperinterpolation provided by equation 13. We remark that our approach in our numericalrepresentation making use of the embedding space basis provides a global description of thedifferential form over the entire surface and a consistent way to obtain derivatives between differentcoordinate charts. For a given chart, a differential form has coordinate components for a 0-formgiven by v (cid:91) = v , a 1-form by v (cid:91) = v i d x i , and a 2-form by v (cid:91) = v ij d x i ∧ d x j . To numerically computederivatives based on these expressions we perform a conversion from the vector field representationin the embedding space to the local coordinate representation of the components.The 1-form presents the most interesting case with the 0-form and 2-form handled similarly.For 1-forms the components v i are related to the components ¯ v k of the vector field representationby v i = g ij v j = g ij a jk ¯ v k . The term a jk converts between the components ¯ v k given in the coordinatesof the embedding space ι , ι , ι to the components v j given in the local coordinates on thesurface ∂ θ , ∂ φ . The exterior derivative of a 1-form can be expressed in coordinate components as dv (cid:91) = ∂ s v i d x s ∧ d x i . For numerical calculations at a given location x (cid:96) (Lebedev node), we choose anappropriate coordinate chart that is locally non-degenerate. We compute the component derivativeas ∂ s v i = ( ∂ s g ij ) a jk ¯ v k + g ij ( ∂ s a jk )¯ v k + g ij a jk ( ∂ s ¯ v k ) . (14)The first two terms only depend on the geometry of the manifold and only the values of thedifferential form at location x (cid:96) (Lebedev node). This can be obtained readily from the sphericalharmonics representation of the geometry of the manifold see Appendix A and Appendix B. Incontrast, the last term depends on the derivatives of the coordinate components and requires use ofthe continuum representation from the hyperinterpolation obtained in equation 13. Putting thistogether with equation 14 and the coordinate expression for the exterior derivative, we obtain anumerical exterior derivative operator d for 1-forms. The case of 0-form and 2-form can be handledsimilarly. We should mention that the case of a 2-form in two dimensions has exterior derivative zerowhich we also impose in our numerical calculations. In this manner, we obtain a numerical operator 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 furtherapplication of other numerical exterior calculus operations.We remark that in practical implimentations of the numerical exterior derivative operator d itis convenient to represent the coordinate conversion in matrix-vector notation as v = GA − ¯v . Thematrix entries [ A − ] jk = a jk correspond to the change from the coordinates of the embedding space tothe local coordinates of the tangent space. With this notation, the derivatives in local coordinates canbe expressed as ∂ s v = ( ∂ s G ) A − ¯v + G ( ∂ s A − ) ¯v + GA − ( ∂ s ¯v ). To avoid differentiating componentsof the inverse matrix A − , we use the identity ∂ s A − = − A − ( ∂ s A ) A − and use a linear algebrasolver to compute the action of A − . This is done at each Lebedev node x (cid:96) with an appropriatechoice made for the coordinate chart that is locally non-degenerate. For additional informationon how we define the coordinate charts and how we perform practical computations from ourrepresentations of the manifold geometry see Appendix B.We approximate next the Hodge star (cid:63) operator on differential forms. We remark that in therelated area of discrete exterior calculus (DEC) efforts are made to preserve geometric structure in Page 8 of 24 he discrete setting, often on triangulated meshes. Interesting issues arise in DEC from the discretegeometry with which one must grapple and extensive studies have been conducted to formulategood approximations for the Hodge star (cid:63) operator [14, 23]. Here we avoid many of these issuessince we treat the operator at the continuum level and have more geometric information availableto us from our spectral representation of both the manifold and the differential forms.We approximate the Hodge star (cid:63) operator on differential forms by a numerical operator (cid:63) which makes use of the representation at the Lebedev nodes. The Hodge star (cid:63) has the featurethat it is a local operation that involves values of the differential form and metric tensor only at anindividual Lebedev node x (cid:96) . We obtain a numerical operator (cid:63) by applying the isomorphisms andmetric tensor using equation 1 and equation 3. The main consideration numerically is to choosewell the coordinate chart so it is locally non-degenerate. The approximation enters through thefidelity of the metric tensor computed from our representation of the manifold geometry. Thegeometry for the radial manifold is determined by the radial function r ( θ, φ ) which is represented asan expansion in a finite number of spherical harmonics up to order L , r ( θ, φ ) = (cid:80) i ˆ r i Y i ( θ, φ ). Thisprovides in a given coordinate chart at location x (cid:96) the metric tensor and curvature tensor alongwith the local coordinate basis vectors ∂ θ , ∂ φ and their derivatives see Appendix B. The numericalHodge star operator (cid:63) maps k -forms defined at the Lebedev nodes to ( n − k )-forms defined at theLebedev nodes. This provides a convenient map between representations for further application ofnumerical exterior calculus operations. In this manner, the numerical exterior derivative d operatorand numerical Hodge star (cid:63) operator can be used through composition to numerically approximatemore complex exterior calculus operations on the manifold. We consider differential equations on the manifold of the form L u = − g. (15)The L denotes a linear differential operator that can be expressed in terms of a composition of theexterior calculus operations d and (cid:63) . For example, in the case of the Laplace-Beltrami equationthis would correspond to the operator L = − δ d = − (cid:63) d (cid:63) d . We discretize the operator L byusing a composition of the numerical exterior calculus operators (cid:63) and d to obtain ˜ L . For theLaplace-Beltrami equation, this corresponds to ˜ L = − ¯ δ ¯d = − (cid:63) d (cid:63) d . We approximate equation 15on the manifold by (cid:104) ˜ L ¯ u, Y i (cid:105) Q = −(cid:104) ¯ g, Y i (cid:105) Q . (16)The ¯ u = (cid:80) j ˆ u j Y j , ¯ g = (cid:80) j ˆ g j Y j denote expansions up to order (cid:98) L/ (cid:99) and (cid:104)· , ·(cid:105) Q denotes the Lebedevinner-product computed up to order L in equation 12. This provides a Galerkin approximationwhere some additional sources of approximation arise from the treatment of the differential operatorby ˜ L . The approximation can be expressed in terms of the solution coefficients ˆu as K ˆ u = − M ˆ g . (17)The ˆu and ˆg denote the collection of coefficients ˆ u j and ˆ g j in the expansion of ¯ u and ¯ g . The K denotes the stiffness matrix with entries K ij = (cid:104) ˜ L Y j , Y i (cid:105) Q and M denotes the mass matrix with Page 9 of 24 ntries M ij = (cid:104) Y j , Y i (cid:105) Q . We again mention that throughout the calculations we use expansions offunctions up to order (cid:98) L/ (cid:99) and Lebedev quadratures of order L . This is important to yield theneeded over-sampling of functions to ensure accurate computation of the inner-product (cid:104)· , ·(cid:105) Q .
4. Convergence Results
We discuss how our numerical methods converge in approximating the fundamental exterior calculusoperations of the exterior derivative d and the Hodge star (cid:63) when applied to 0-forms, 1-forms and2-forms. We then discuss the convergence of our methods for compositions of operators and presentresults for the Laplace-Beltrami equation.We remark that throughout our convergence studies, we describe test functions using for apoint x on the manifold its location within the embedding space. To perform calculations we makeuse of the embedding space coordinates [ x, y, z ] corresponding to x = x ι + y ι + z ι , where ι , ι , ι is the basis for the embedding space. In this manner our test data is not tied to a specific choiceof local coordinates on the manifold. All figures report the relative error (cid:15) rel = (cid:107) ¯w (cid:93) − w (cid:93) (cid:107) / (cid:107) w (cid:93) (cid:107) .The w is the exact result and ¯w is the numerically computed result. For (cid:107) · (cid:107) , we use the L -normof the embedding space. We approximate the Hodge star (cid:63) by the numerical operator (cid:63) using the hyperinterpolation approachdiscussed in Section 3.2. We investigate in practice the accuracy of this approach on a few differentgeometries and differential forms.We first consider a 0-form defined on Manifold B defined in Figure 1. We take the 0-form tobe f = exp( z ) / (3 − y ). We show the accuracy of our numerical operator (cid:63) in approximating theHodge star (cid:63) as the number of Lebedev nodes increases. We find that the main limitation in theaccuracy of the (cid:63) is the resolution of the geometry of the manifold. This is seen in our results whereonce a sufficient number of spherical harmonic modes are reached the relative error rapidly decaysin approximating (cid:63)f . The convergence of (cid:63) as the number of Lebedev nodes is increased and whenthe geometry is varied is shown in Figure 4. Figure 4:
Convergence of the numerical Hodge star operator (cid:63) for -forms. We show for Manifold Bhow the relative error of (cid:63)f in approximating (cid:63)f as the number of Lebedev nodes increases. The -formis f = exp( z ) / (3 − y ) . We investigate how the manifold geometry influences convergence by varying theamplitude r in the range [0 . , . for Manifold B. The amplitude r = 0 . corresponds to a sphere and r = 0 . to the final shape of Manifold B shown in Figure 1. Page 10 of 24 e next consider on Manifold B the 1-form α = (cid:112) | g | exp( z ) dθ + (cid:112) | g | exp( z ) dφ . We again findthat the main limitation in the accuracy of the (cid:63) is the resolution of the geometry of the manifold.In this case we find the error rapidly decreases once a sufficient number spherical harmonic modesare used. The convergence of (cid:63) as the number of Lebedev nodes is increased and when the geometryis varied is shown in Figure 5. Figure 5:
Convergence of the numerical Hodge star operator (cid:63) for -forms. We show for Manifold Bthe relative error of (cid:63)α in approximating (cid:63)α as the number of Lebedev nodes increases. The -form is α = (cid:112) | g | exp( z ) dθ + (cid:112) | g | exp( z ) dφ . We investigate how the manifold geometry influences convergence byvarying the amplitude r in the range [0 . , . for Manifold B. The amplitude r = 0 . corresponds to asphere and r = 0 . to the final shape of Manifold B shown in Figure 1. These results indicate that the numerical operator (cid:63) provides an accurate approximation tothe Hodge star (cid:63) . Given the localized nature of the Hodge star, the main consideration to obtainaccurate results with (cid:63) is to use enough Lebedev nodes to ensure sufficient resolution of the geometryof the manifold.
We investigate the convergence of the numerical operator d in approximating the exterior derivative d when applied to 0-forms and 1-forms. We do not consider 2-forms here since in two dimensionsthe exterior derivative would be zero which we also impose in our numerical calculations [1, 31]. Weconsider the 0-form given by f = exp( z ) and the 1-form given by α = | g | exp( z ) dθ + | g | exp( z ) dφ .We consider for Manifold B the relative error of d in approximating d as the number of Lebedevnodes is increased and the manifold geometry is varied. These results are shown in Figure 6 andFigure 7. Page 11 of 24 igure 6:
Convergence of the numerical exterior derivative operator d for -forms. We show for ManifoldB the relative error of d f in approximating d f as the number of Lebedev nodes increases. The -form is f = exp( z ) . We investigate how the manifold geometry influences convergence by varying the amplitude r inthe range [0 . , . for Manifold B. The amplitude r = 0 . corresponds to a sphere and r = 0 . to the finalshape of Manifold B shown in Figure 1. We see that the numerical operator d converges spectrally in approximating the exteriorderivative d both for the 0-forms and 1-forms. Interestingly, we see that in the case when r = 0 . r (cid:54) = 0. This occurs since r = 0 corresponds to the case when the shape is a sphere where the geometry is relatively simpleand many of the geometric terms to be numerically approximated greatly simplify. The main sourceof error in this case arises primarily from the hyperinterpolation used for computing the derivatives.In the case when r (cid:54) = 0, the isogeometric approach used to compute d approximates the geometryof the manifold using a finite spherical harmonics representation which results in an additionalsource of approximation error. Figure 7:
Convergence of the numerical exterior derivative operator d for -forms. We show for ManifoldB the relative error of d α in approximating d α as the number of Lebedev nodes increases. The -form is α = | g | exp( z ) dθ + | g | exp( z ) dφ . We investigate how the manifold geometry influences convergence by varyingthe amplitude r in the range [0 . , . for Manifold B. The amplitude r = 0 . corresponds to a sphere and r = 0 . to the final shape of Manifold B shown in Figure 1. We investigate how the geometry contributes to convergence by varying r over the range[0 . , . r = 0 . Page 12 of 24 ronounced geometry corresponding to r = 0 .
4. Overall, the results indicate that the numericaloperator d provides for 0-forms and 1-forms an accurate approximation for the exterior derivative d . We consider the composition of exterior calculus operators to represent partial differential equationson manifolds. Taking the approach we discussed in Section 3.3, we consider the Laplace-Beltramioperator L = − δ d = − (cid:63) d (cid:63) d and consider the Poisson problem on the manifold L u = − δ d u = − g. (18)Since there are no boundaries for the manifolds we shall consider, we also impose throughout thatthe zero mode has ˆ u = 0. The Laplace-Beltrami equation 18 provides a test of using compositionsof the numerical exterior calculus operators d and (cid:63) and the associated convergence and accuracy. Page 13 of 24 igure 8:
Solution of the Laplace-Beltrami Equations on Manifold B of Figure 1. For the Laplace-Beltramiequations we show on the manifold surface the source term g (left), solution u (middle), and relative error (cid:15) rel (right) from the vantage point of the coordinate axes associated with the embedding space ι , ι , ι . Shown isthe case with the spherical harmonics resolved with Lebedev nodes.
We investigate how our methods solve the Laplace-Beltrami equation 18 on Manifold B andManifold C of Figure 1. To have a known solution to the Laplace-Beltrami equations with which tocompare, we manifacture a source term for the solution on the manifold given by u = exp( y ) / (3 − z ) ,where the x, y, z refer to the coordinates of the embedding space as discussed in the beginningof Section 4. We compute the source term g = −L u symbolically using the package SymPy andevaluate the expressions numerically when data is needed for the source [22]. In this manner we areable to assess the relative errors of the numerical methods with high-precision. Page 14 of 24 igure 9:
Solution of the Laplace-Beltrami Equations on Manifold B of Figure 1. For the Laplace-Beltramiequations we show on the manifold surface the source term g (left), solution u (middle), and relative error (cid:15) rel (right) from a rotated vantage point. The case shown corresponds to r = 0 . and spherical harmonicsresolved with Lebedev nodes.
We show the source term g , solution u , and relative errors (cid:15) rel in Figure 8 and Figure 9. Weuse the approach discussed for solving PDEs on manifolds in Section 3.3 for the Laplace-Beltramiequations. We show how the composition of the numerical operators d and (cid:63) perform as the numberof Lebedev nodes increases in Figure 10. The results indicate that the numerical methods based oncomposing d and (cid:63) perform well in approximating the true compositions of the exterior calculusoperators (cid:63) and d . Page 15 of 24 igure 10:
Convergence of the numerical methods for the Laplace-Beltrami Equations on Manifold B ofFigure 1. We show for Manifold B the relative error when the Laplace-Beltrami equations are solved using acomposition of d and (cid:63) as discussed in Section 3.3. We investigate how the manifold geometry influencesconvergence by varying the amplitude r in the range [0 . , . for Manifold B. The amplitude r = 0 . corresponds to a sphere and r = 0 . to the final shape of Manifold B shown in Figure 1. We see that when the geometry has r = 0 . r (cid:54) = 0 we see the convergence is comparable in each case. We see theconvergence is slowest for the most pronounced geometry where the representation for the geometryincurs the most approximation error. Overall, we find the compositions of d and (cid:63) perform well andprovide 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 usethe manufactured solution on the surface given by u = exp( y ) / (3 − z ) and g = −L u computedsymbolically. We mention that given the different geometry of Manifold C relative to Manifold Bboth the solution and the source data are different on each of the manifold surfaces. We show on themanifold surface the source term g , solution u , and relative error (cid:15) rel in Figure 11 and Figure 12. Page 16 of 24 igure 11:
Solution of the Laplace-Beltrami Equations on Manifold C of Figure 1. For the Laplace-Beltramiequations we show on the manifold surface the source term g (left), solution u (middle), and relative error (cid:15) rel (right) from the vantage point of the coordinate axes associated with the embedding space ι , ι , ι . Shown isthe case with the spherical harmonics resolved with Lebedev nodes.
Page 17 of 24 igure 12:
Solution of the Laplace-Beltrami Equations on Manifold C of Figure 1. For the Laplace-Beltramiequations we show on the manifold surface the source term g (left), solution u (middle), and relative error (cid:15) rel (right) from a rotated vantage point. The case shown corresponds to r = 0 . and the spherical harmonicsresolved with Lebedev nodes.
For the more geometrically complicated Manifold C, we perform calculations to show how thecomposition of the numerical operators d and (cid:63) perform as the number of Lebedev nodes increasesand r is varied, see Figure 13. Page 18 of 24 igure 13:
Convergence of the numerical methods for the Laplace-Beltrami Equations on Manifold C ofFigure 1. We show for Manifold B the relative error when the Laplace-Beltrami equations are solved using acomposition of d and (cid:63) as discussed in Section 3.3. We investigate how the manifold geometry influencesconvergence by varying the amplitude r in the range [0 . , . for Manifold C. The amplitude r = 0 . corresponds to a sphere and r = 0 . to the final shape of Manifold C shown in Figure 1. We find the numerical methods converge spectrally for Manifold C. As may be expected givenits more complicated geometry, we find the convergence is somewhat slower than for ManifoldB. We find that the numerical methods perform well for Manifold C once a sufficient number ofLebedev nodes are used to resolve the geometry. We again see the trend that the case with smaller r convergence more rapidly. Overall, we find that in practice the composition of the numericaloperator d and (cid:63) performs well in approximating differential operators on the manifold. Thisindicates the approach discussed in Section 3.3 can be applied quite generally in practice for solvingpartial differential equations on manifolds.
5. Conclusion
In summary, we have developed numerical methods for exterior calculus on radial manifolds based onhyperinterpolation and Lebedev quadrature. We have shown our methods are spectrally accurate inapproximating the exterior derivative d , Hodge star (cid:63) , and their compositions. We presented resultsfor the Laplace-Beltrami equations that demonstrate the utility of the approach. The introducednumerical methods can be applied quite generally for approximating exterior calculus operationsand differential equations on radial manifolds. Page 19 of 24 . 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.
References [1] R. Abraham, J.E. Marsden, and T.S. Raiu.
Manifolds, Tensor Analysis, and Applications .Number v. 75. Springer New York, 1988.[2] Congpei An, Xiaojun Chen, Ian H. Sloan, and Robert S. Womersley. Regularized least squaresapproximations on the sphere using spherical designs.
SIAM Journal on Numerical Analysis ,50(3):1513–1534, 2012.[3] Michael Holst Andrew Gillette and Yunrong Zhu. Finite element exterior calculus for evolutionproblems.
Journal of Computational Mathematics , 2017.[4] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Finite element exterior calculus,homological techniques, and applications. 15:1–155–, 2006.[5] Marino Arroyo and Antonio DeSimone. Relaxation dynamics of fluid membranes.
Phys. Rev.E , 79(3):031915–, March 2009.[6] Kendall Atkinson and Weimin Han.
Spherical Harmonics and Approximations on the UnitSphere: An Introduction . Springer, 2010.[7] Fernando de Goes, Mathieu Desbrun, Mark Meyer, and Tony DeRose. Subdivision exteriorcalculus for geometry processing.
ACM Trans. Graph. , 35(4):1–11, 2016.[8] Richard S. Falk Douglas N. Arnold and Ragnar Winther. Finite element exterior calculus: fromhodge theory to numerical stability.
Bull. Amer. Math. Soc. , 2010.[9] J.R. Driscoll and D.M. Healy. Computing fourier transforms and convolutions on the 2-sphere.
Advances in Applied Mathematics , 15(2):202–250, June 1994.[10] James Eells. Geometric aspects of currents and distributions.
Proceedings of the NationalAcademy of Sciences , 41(7):493–496, 1955.[11] Mathieu Desbrun Peter Schroder Fernando de Goes, Keenan Crane. Digital geometry processingwith discrete exterior calculus. In
SIGGRAPH , 2013.[12] D.M. Healy, D.N. Rockmore, P.J. Kostelec, and S. Moore. Ffts for the2-sphere-improvementsandvariations.
Journal of Fourier Analysis and Applications , 9(4):341–385, 2003.[13] Anil N. Hirani.
Discrete Exterior Calculus . PhD thesis, Caltech, 2003.[14] Anil N. Hirani, Kaushik Kalyanaraman, and Evan B. VanderZee. Delaunay hodge star.
Computer-Aided Design , 45(2):540 – 544, 2013. Solid and Physical Modeling 2012.
Page 20 of 24
15] Yuri Bazilevs J. Austin Cottrell, Thomas J. R Hughes.
Isogeometric Analysis: TowardIntegration of CAD and FEA . Wiley, 2009.[16] Eva Kanso, Marino Arroyo, Yiying Tong, Arash Yavari, Jerrold G. Marsden, and MathieuDesbrun. On the geometric character of stress in continuum mechanics.
Zeitschrift fr angewandteMathematik und Physik , 58(5):843–856, 2007.[17] V. I. Lebedev and D. N. Laikov. A quadrature formula for the sphere of the 131st algebraicorder of accuracy.
Dokl. Math. , 59:477481, 1999.[18] V.I. Lebedev. Quadratures on a sphere.
USSR Computational Mathematics and MathematicalPhysics , 16(2):10–24, 1976.[19] David Bau III Lloyd N. Trefethen.
Numerical Linear Algebra . SIAM, 1997.[20] J.E. Marsden and T.J.R. Hughes.
Mathematical Foundations of Elasticity . Dover, 1994.[21] MELVIN LEOK MATHIEU DESBRUN, ANIL N. HIRANI and JERROLD E. MARSDEN.Discrete exterior calculus.
Technical Report , 2003.[22] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondˇrej ˇCert´ık, Sergey B. Kirpichev,Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, ThilinaRathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta,Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, ˇStˇep´anRouˇcka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and AnthonyScopatz. Sympy: symbolic computing in python.
PeerJ Computer Science , 3:e103, January2017.[23] Mamdouh S. Mohamed, Anil N. Hirani, and Ravi Samtaney. Comparison of discrete hodgestar operators for surfaces.
Comput. Aided Des. , 78(C):118–125, September 2016.[24] A. Pressley.
Elementary Differential Geometry . Springer, 2001.[25] Manfred Reimer. Hyperinterpolation on the sphere at the minimal projection order.
Journalof Approximation Theory , 104(2):272–286, 2000.[26] Dzhelil Rufat, Gemma Mason, Patrick Mullen, and Mathieu Desbrun. The chain collocationmethod: A spectrally accurate calculus of forms.
Journal of Computational Physics , 257, PartB:1352–1372, January 2014.[27] Nathanal Schaeffer. Efficient spherical harmonic transforms aimed at pseudospectral numericalsimulations.
Geochem. Geophys. Geosyst. , 14(3):751–758, 2013.[28] Jon Karl Sigurdsson and Paul J. Atzberger. Hydrodynamic coupling of particle inclusionsembedded in curved lipid bilayer membranes.
Soft Matter , 12(32):6685–6707, 2016.[29] Ian H. Sloan and Robert S. Womersley. Constructive polynomial approximation on the sphere.
Journal of Approximation Theory , 103(1):91–118, 2000.
Page 21 of 24
30] I.H. Sloan. Polynomial interpolation and hyperinterpolation over general regions.
Journal ofApproximation Theory , 83(2):238–254, 1995.[31] Michael Spivak.
Calculus On Manifolds: A Modern Approach To Classical Theorems OfAdvanced Calculus . Westview Press, 1971.[32] Ari Stern, Yiying Tong, Mathieu Desbrun, and Jerrold E. Marsden. Geometric computationalelectrodynamics with variational integrators and discrete differential forms. In Dong EuiChang, Darryl D. Holm, George Patrick, and Tudor Ratiu, editors,
Geometry, Mechanics, andDynamics: The Legacy of Jerry Marsden , pages 437–475. Springer New York, New York, NY,2015.[33] G. Strang.
Linear Algebra and Its Applications . Academic Press, Inc., 1980.[34] Ke Wang, Weiwei, Yiying Tong, Mathieu Desbrun, and Peter Schroder. Edge subdivisionschemes and the construction of smooth vector fields. In
ACM SIGGRAPH 2006 Papers , pages1041–1048, Boston, Massachusetts, 2006. ACM.[35] Hassler Whitney.
Geometric Integration Theory . Princeton University Press, 1957.[36] Robert S. Womersley and Ian H. Sloan. How good can polynomial interpolation on the spherebe?
Advances in Computational Mathematics , 14(3):195–226, 2001.
AppendixA. Spherical Harmonics
The spherical harmonics are given by Y mn ( θ, φ ) = (cid:115) (2 n + 1)( n − m )!4 π ( n + m )! P mn (cos( φ )) exp ( imθ ) (19)where m denotes the order and n the degree for n ≥ m ∈ {− n, . . . , n } . The P mn denote the Associated Legendre Polynomials . In our notation, θ denotes the azmuthal angle and φ the polarangle 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 mn = Y − mn . We have found it convenient to represent the sphericalharmonics as Y mn ( θ, φ ) = X mn ( θ, φ ) + iZ mn ( θ, φ ) (20)where X mn and Z mn denote the real and imaginary parts. In our numerical methods we use thissplitting 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) basis elements. In the case of N = 2 we have the basiselements˜ Y = Y , ˜ Y = Z , ˜ Y = Y , ˜ Y = X , ˜ Y = Z , ˜ Y = Z , ˜ Y = Y , ˜ Y = X , ˜ Y = X . (21) Page 22 of 24 e 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 = ˜ Y i / (cid:113) (cid:104) ˜ Y i , ˜ Y i (cid:105) .We compute derivatives of our finite expansions by evaluating analytic formulas for the sphericalharmonics in order to try to minimize approximation error [6]. Approximation errors are incurredwhen sampling the values of expressions involving these derivatives at the Lebedev nodes andperforming quadratures. The derivative in the azmuthal coordinate θ of the spherical harmonics isgiven by ∂ θ Y mn ( θ, φ ) = ∂ θ (cid:115) (2 n + 1)( n − m )!4 π ( n + m )! P mn (cos( φ )) exp ( imθ ) = imY mn ( θ, φ ) . This maps the spherical harmonic of degree n to again a spherical harmonic of degree n . In ournumerics, this derivative can be represented in our finite basis which allows us to avoid projections.This allows for computing the derivative in θ without incurring an approximation error. For thederivative in the polar angle φ we have that ∂ φ Y mn ( θ, φ ) = m cot( φ ) Y mn ( θ, φ ) + (cid:112) ( n − m )( n + m + 1) exp ( − iθ ) Y m +1 n ( θ, φ ) . (22)We remark that the expression can not be represented in terms of a finite expansion of sphericalharmonics. We use this expression for ∂ φ Y mn ( θ, φ ) when we compute values at the Lebedev quadraturenodes in equation 14. This provides a convenient way to compute derivatives of differential formsfollowing the approach discussed in Section 3. We remark that it is the subsequent hyperinterpolationof the resulting expressions where the approximation error is incurred. We adopt the notationalconvention that Y mn = 0 when m ≥ n + 1. For further discussion of spherical harmonics see [6]. B. Differential Geometry of Radial Manifolds
We consider throughout manifolds of radial shape. A radial manifold is defined as a surface whereeach point can be connected by a line segment to the origin without intersecting the surface. Inspherical coordinates, any point x on the radial manifold can be expressed as x ( θ, φ ) = σ ( θ, φ ) = r ( θ, φ ) r ( θ, φ ) (23)where r is the unit vector from the origin to the point on the sphere corresponding to angle θ, φ and r is a positive scalar function.We take an isogeometric approach to representing the manifold M . We sample the scalarfunction r at the Lebedev nodes and represent the geometry using the finite spherical harmonicsexpansion r ( θ, φ ) = (cid:80) i ¯ r i Y i up to the order (cid:98) L/ (cid:99) where ¯ r i = (cid:104) r, Y i (cid:105) Q for a quadrature of order L .We consider two coordinate charts for our calculations. The first is referred to as Chart Aand has coordinate singularities at the north and south pole. The second is referred to as ChartB and has coordinate singularities at the east and west pole. For each chart we use sphericalcoordinates with ( θ, φ ) ∈ [0 , π ) × [0 , π ] but to avoid singularities only use values in the restrictedrange φ ∈ [ φ min , φ max ], where 0 < φ min ≤ π , and π ≤ φ max < π . In practice, one typically takes φ min = 0 . × π and φ max = 0 . × π . For Chart A, the manifold is parameterized in the embeddingspace R as x (ˆ θ, ˆ φ ) = r (ˆ θ, ˆ φ ) r (ˆ θ, ˆ φ ) , r (ˆ θ, ˆ φ ) = (cid:2) sin( ˆ φ ) cos(ˆ θ ) , sin( ˆ φ ) sin(ˆ θ ) , cos( ˆ φ ) (cid:3) (24) Page 23 of 24 nd for Chart B x (¯ θ, ¯ φ ) = r (¯ θ, ¯ φ ) r (¯ θ, ¯ φ ) , ¯r (¯ θ, ¯ φ ) = (cid:2) cos( ¯ φ ) , sin( ¯ φ ) sin(¯ θ ) , − sin( ¯ φ ) cos(¯ θ ) (cid:3) . (25)With these coordinate representations, we can derive explicit expressions for geometric quantitiesassociated with the manifold such as the metric tensor and shape tensor. The derivatives used asthe basis ∂ θ , ∂ φ for the tangent space can be expressed as σ θ ( θ, φ ) = r θ ( θ, φ ) r ( θ, φ ) + r ( θ, φ ) r θ ( θ, φ ) (26) σ φ ( θ, φ ) = r φ ( θ, φ ) r ( θ, φ ) + r ( θ, φ ) r φ ( θ, φ ) . (27)We have expressions for r θ and r φ in the embedding space R using equation 24 or equation 25depending on the chart being used. The first fundamental form I (metric tensor) and secondfundamental form II (shape tensor) are given by I = (cid:20) E FF G (cid:21) = (cid:20) σ θ · σ θ σ θ · σ φ σ φ · σ θ σ φ · σ φ (cid:21) = (cid:20) r θ + r sin( φ ) r θ r φ r θ r φ r φ + r (cid:21) . (28)and II = (cid:20) L MN N (cid:21) = (cid:20) σ θθ · n σ θφ · nσ φθ · n σ φφ · n (cid:21) . (29)The n denotes the outward normal on the surface and is computed using n ( θ, φ ) = σ θ ( θ, φ ) × σ φ ( θ, φ ) (cid:107) σ θ ( θ, φ ) × σ φ ( θ, φ ) (cid:107) . (30)The terms σ θθ , σ θφ , and σ φ,φ are obtained by further differentiation from equation 26 and equation 27.We use the notation for the metric tensor g = I interchangably. In practical calculations wheneverwe 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 notationfor the metric tensor g ij and its inverse g ij which has the formal correspondence g ij = [ I ] i,j , g ij = (cid:2) I − (cid:3) i,j . (31)For the metric factor we also have that (cid:112) | g | = (cid:112) det( I ) = r (cid:113) r θ + ( r φ + r ) sin( φ ) = (cid:107) (cid:126)σ θ ( θ, φ ) × (cid:126)σ φ ( θ, φ ) (cid:107) . (32)To ensure accurate numerical calculations in each of the above expressions the appropriate coordinateseither Chart A or Chart B are used to ensure sufficient distance from coordinate singularities atthe poles. To compute quantities associated with curvature of the manifold we use the Weingartenmap [24] which can be expressed as W = − I − II . (33)To compute the Gaussian curvature K , we use K ( θ, φ ) = det ( W ( θ, φ )) . (34)For further discussion of the differential geometry of manifolds see [1, 24, 31].(34)For further discussion of the differential geometry of manifolds see [1, 24, 31].