A Novel Finite Element Discretization of Domains with Spheroidal Geometry Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. Necibe Tuncer Certificate of Approval: Paul Schmidt Professor Mathematics and Statistics Amnon J. Meir, Chair Professor Mathematics and Statistics Georg Hetzer Professor Mathematics and Statistics Richard Zalik Professor Mathematics and Statistics George T. Flowers Interim Dean Graduate School A Novel Finite Element Discretization of Domains with Spheroidal Geometry Necibe Tuncer A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 10, 2007 A Novel Finite Element Discretization of Domains with Spheroidal Geometry Necibe Tuncer Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Dissertation Abstract A Novel Finite Element Discretization of Domains with Spheroidal Geometry Necibe Tuncer Doctor of Philosophy, May 10, 2007 (M.A.,Dokuz Eylul University, 2001) (B.S., Dokuz Eylul University, 1999) 66 Typed Pages Directed by Amnon J. Meir We describeand analyze a newfiniteelement discretizations for domainswith spheroidal geometry. In particular, we describe how the method can be used to approximate solutions as well as eigenvalues and eigenfunctions of partial differential equations posed on the sphere, ellipsoidal shell, and cylindrical shell. These novel, so-called, ?radially projected finite ele- ments? are particularly attractive for numerical simulations since the resulting finite element discretization is ?logically rectangular? and may be easily implemented or incorporated into existing finite element codes. iv Style manual or journal used Journal of Approximation Theory (together with the style known as ?aums?). Bibliograpy follows van Leunen?s A Handbook for Scholars. Computer software used The document preparation package TEX (specifically LATEX) together with the departmental style-file aums.sty. v Table of Contents List of Figures vii 1 Introduction 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Our Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 The Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Discretization of the Domain and Basis Elements . . . . . . . . . . . 7 1.2.2 Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 Smooth Manifolds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 Eigenvalue Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2 The Radial Projection 24 2.1 Properties of the Radial Projection . . . . . . . . . . . . . . . . . . . . . . . 26 2.2 Some Examples of Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.1 Mesh on the Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.2 Meshes on Ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.3 Mesh on the Cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.2.4 Mesh on the Disc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3 Analysis 37 3.1 Finite Element Construction . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4 Numerical Experimetns 50 4.1 Finite elements on the sphere . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.4 Example 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5 Conclusion 57 Bibliography 58 vi List of Figures 1.1 A uniform triangulation of the square . . . . . . . . . . . . . . . . . . . . . 9 1.2 Linear basis element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 A transition map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Mesh on the box and the sphere. . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Faces of the box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.3 F+1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4 F+1 and P (F+1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Geodesic distance between a and b. . . . . . . . . . . . . . . . . . . . . . . . 30 2.6 Geodesic and Euclidean distance of a and b . . . . . . . . . . . . . . . . . . 32 2.7 Mesh on the ellipsoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.8 Mesh on the cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.9 Mesh on the disc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.1 Spherical triangle, K and planar triangle Kh. . . . . . . . . . . . . . . . . . 42 4.1 Basis functions on the sphere . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Numerical approximation to ?(a) = cosa1. . . . . . . . . . . . . . . . . . . . 53 vii Chapter 1 Introduction 1.1 Overview Mathematical models of many engineering and scientific problems usually consist of partial differential equations and some given conditions. The solution of such problems can be approximated by numerical solution techniques. Scientific computing gained significant importance in studies of physical problems. The Finite Element Method (FEM) is a nu- merical approach to approximate the solution of a partial differential equation. Compared to the Finite Difference Method, the FEM is a fairly new method. FEM has become a favorite method immediately, since it is very easy to use even for complex domains, for general boundary value problems, and also has a strong theoretical basis which yields easy derivation of error estimates. The FEM is very nicely structured so that a general computer program can be used to approximate solutions of various problems easily. Despite all the ad- vantages of the FEM, it is mainly convenient for domains that have a polyhedral geometry. Thus, it has mainly been applied to the polygonal and polyhedral domains. Many studies have been performed for FEM disretization of curved domains or domains with spherical geometry. However, the results of these studies have focused on approximating the domain by a polygonal domain, rather than applying FEM to these exact domains. Recently, there has been growing interest in finite element discretization of partial differential equations de- fined on the sphere, since these equations have many applications in areas such as climate modeling or weather forecast and various engineering problems. 1 1.1.1 Background Mathematical weather models or climate models consist of partial differential equations posed on a domain which is assumed to be a sphere for the obvious reason. Lately, the study of numerically approximating the solution of partial differential equations defined on the sphere has become a favorite subject. However, this is a challenging problem no matter which numerical method is being applied to approximate the solutions of the equations because of the geometry of the sphere. Finite difference methods, spectral methods, finite element methods, finite volume methods are all numerical techniques that have been applied to approximate the solutions of the partial differential equations on the sphere. Grid gener- ation is an important part of approximating solutions of partial differential equations, since the accuracy of numerical solutions depend on the quality of the grid. Various researchers have used different methods for mesh generations on the sphere, such as, constructing the finite element basis directly on spherical triangles by using barycentric coordinate systems [10], approximating, not just the sphere, but any surface with a polyhedral surface, [9] also [6] and [7]. Both of these approaches, described in [10],[9], yield the optimal convergence rate (which is 2 if linear finite elements are used to approximate the solution, the error is measured in the L2 norm, and the solution is smooth enough). Although all these ap- proaches result in a nice mesh, they suffer from some weaknesses, such as not discretizing the sphere exactly or requiring multiple computations for each refinement step. Our work aims to overcome these difficulties by developing a new method that discretizes the sphere exactly and that can also be easily implemented and incorporated into existing finite element codes. 2 1.1.2 Our Contribution Our work is motivated by a desire for exact finite element discretization of the sphere which is yielding a conforming finite element method. We develop and analyze a method which can also be easily implemented and incorporated into an existing code. In our ap- proach, we disceretize the sphere exactly, rather than approximate it, which has not been done before. The method is also very easy to implement if one knows how to implement finite element methods on planar domains. It requires very little computational effort since all computations are done on a logically rectangular domain. Our primary objective is to develop a finite element discretization of a sphere or a ball. Our method can easily be gen- eralized to other domains, such as cylinders, ellipsoids or tori, and can be used for problems posed on domains with spherical, cylindrical, or toroidal holes. In this work we present our approach for the discretization of the sphere, and the same method can be also used for cylinders, ellipsoids and tori. Basically, one can apply this discretization to any domain that has a spherical, cylindrical, ellipsoidal, or toroidal geometry, which we call spheroidal geometry. We develop and analyze a novel finite element discretization of domains with spheroidal geometry .We present some computational examples related to some of these domains. We firstgive an introduction to finite element method in section 1.2, introduce smooth manifolds and Laplace Beltrami operator in section 1.3. In chapter 2 we introduce our new method, in chapter 3 we analyze the method, derive error estimates, and in chapter 4 we present numerical results from computer experiments. 3 1.2 The Finite Element Method In this section, we introduce the finite element method and illustrate it with the classical Laplace?s equation on a bounded, connected, open domain ? in Rn. Consider the following boundary value problem, ??u = f in ? (1.1) u = 0 on ?? (1.2) which is also known as Poisson?s equation. In (1.1), ? is the Laplacian, ?? is the boundary of the domain ?, and u is the unknown function. The right hand side f is a real valued function defined on the open subset ? of Rn, we say that f ? Cr(?), if all the kth order (classical) partial derivatives exist and are continuous for k ? r, where r,k are positive integers. We call any function u ? C2(?) that satisfies the differential equation and the boundary condition a classical solution. Consider the Lebesgue measure of Rn, let ? be a Lebesgue measurable domain, and let f be a Lebesgue-integrable function, then we denote the Lebesgue integral by integraldisplay ? fd? where d? denotes the Lebesgue measure. We say that f is in Lp(?) , for 1 ? p < ? ifparenleftbiggintegraldisplay ? |f|pd? parenrightbigg1/p 0 such that a(v,v) ??bardblvbardbl2H1(?) ?v ?H10(?). We say that the the bilinear form is continuous or bounded if there exists a constant C > 0 such that a(u,v) ?CbardblubardblH1 0(?) bardblvbardblH1 0(?) . Clearly, a(?,?) induces a norm, which is equivalent to the norm of the Hilbert space H10(?). This norm bardblvbardbla = radicalbig a(v,v) is called the energy norm. The question we need to answer is the relationship between (1.1) and its variational problem (1.3). Let f ? L2(?) and u ? H2(?) then u satisfies (1.3) if and only if u is a solution of (1.1). Within the theory of functional analysis, one can prove the following theorem, which is known as Riesz Representation theorem in the literature [12]. 6 Theorem 1.1 Let H be a Hilbert space and let L be a continuous linear functional L : H ? R, s.t. L(u) = (u,v) ?u?H then L is an isomorphism. Theorem 1.2 There exists a unique solution to (1.3). Proof: The bilinear form a(?,?) is a bounded, elliptic such that (H10(?),a(?,?))is a Hilbert space. Then the theorem follows as a consequence of Riesz representation theorem. square To approximate the solution of (1.3), we construct the discrete variational problem in a finite dimensional subspace of H10(?). Let ? be the finite dimensional subspace of H10(?), then the discrete analog of problem (1.3) is : Find uh ?? such that integraldisplay ? ?uh ??vhd? = integraldisplay ? f ?vhd? ?vh ??. (1.4) It is worth noting that by (1.4), we construct a discrete scheme for approximating the solution of (1.1). Also, at this point we need to remark that (1.4) has a unique solution, since ? together with the bilinear form, ah(uh,vh) :=integraltext??uh ??vhd? is a Hilbert space. 1.2.1 Discretization of the Domain and Basis Elements Finite element method can be used to solve the discrete problem (1.4). One important step is to define, ?, the finite dimensional subspace of H10(?) . We define ? by first parti- tioning the domain. For simplicity, hereinafter let ? be a convex polygonal domain in R2. Let K1,K2,...,Kn be triangles or quadrilaterals such that they form a partition of ?, i.e ? = nuniondisplay i=1 Ki and ?Ki intersectiondisplay ? Kj = ?. 7 For any intersection of two triangles or quadrilaterals Ki and Kj, if the intersection consists of one point then it is a common vertex, if the intersection consists of more than one point then it is a common edge. Let T j = {Kj1,Kj2,...,Kjnj} ,j ? N denote the family of triangulation where j denotes the refinement step. We say that the mesh size of the triangulation T at the jth refinement step is hj where hj = max Kji?T j hKj i , and hKj i = diam(Kji), here diam(Kji) = sup x,y?Kji bardblx?ybardbl. The quantity hj is a measure of how refined the mesh is. The smaller hj is, the finer the mesh. Let ?Kj i denote the radius of inscribed circle in Kji. We say that a family of triangulation T j, j ? N is shape regular if there exists a constant ? > 0, independent of j such that for each Kji ? T j, we have hKj i ?Kj i ??. Figure 1.1 shows a uniform triangulation of the square. Usually we choose the finite dimensional subspace to be the space of piecewise polynomials on the domain ?. Piece- wise linear or quadratic functions are widely used in defining these subspaces. Figure 1.2 shows the linear basis element. For instance, define ? to be the space of piecewise linear polynomials, that is: ? = {? :? is piecewise linear continuous polynomials and ? = 0 on ??}. 8 ?1 ?0.5 0 0.5 1?1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 0.8 1 Figure 1.1: A uniform triangulation of the square Let {?i}ni be a basis for ?, then for any uh ? ? we have uh = nsummationdisplay i ui?i. Thus, (1.4) reduces to a linear system of equations of the following form Au = f where A is a n?n matrix whose entries are Aij = ah(?j,?i), and u and f are n?1 vectors of the form u(i) = ui, f(i) = (f,?i).Since the bilinear form ah(?,?) is symmetric and elliptic, the matrixA is symmetric and positive definite, which assures the existence and uniqueness of a solution of (1.4). The matrix A is also a sparse matrix, if the basis functions are chosen to have compact support. In the following section we talk about the error estimates in the finite element method. 9 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 10 0.5 1 Figure 1.2: Linear basis element 1.2.2 Error Estimates Let uh be the solution of (1.4) and u be the exact solution of (1.3), we first observe the orthogonality between u?uh and uh in terms of the energy inner product. Simply, subtracting, (1.4) from (1.3), we get, a(u?uh,vh) = 0, ?vh ??. (1.5) The following lemma, which is known as Cea?s lemma, plays a significant role in estab- lishing error estimates in finite element method. Lemma 1.1 (Cea?s Lemma) Let uh be the solution of (1.4) and u be the exact solution of (1.3), then bardblu?uhbardblH1(?) ? C? infv h?? bardblu?vhbardblH1(?). 10 Proof: Thebilinear formis elliptic,a(u?uh,u?uh) ??bardblu?uhbardbl2H1(?) anda(u?uh,u?uh) = a(u?uh,u?vh), since a(u?uh,uh?vh) = 0 from (1.5) and the bilinear form is continuous, that is, a(u?uh,u?vh) ? Cbardblu?uhbardblH1(?)bardblu?vhbardblH1(?) ?bardblu?uhbardbl2H1(?) ? Cbardblu?uhbardblH1(?)bardblu?vhbardblH1(?) bardblu?uhbardblH1(?) ? C?bardblu?vhbardblH1(?). square Cea?s Lemma states that the solution of (1.4) is a best approximation to the solution of (1.3) by an element of ?. Since, ? is a convex polygonal domain in R2, by Sobolev Embed- ding theorem, H2(?) is continuously imbedded in C(??). Thus, for every u?H2(?), there exists a uniquely defined interpolation operator, Iu, associated with the nodal points. In the theory of finite elements one can prove the following theorem which states the convergence rates of the finite element method [11]. Theorem 1.3 Let u ? H2(?) and T be a shape regular triangulation of ?, and I be an interpolation operator I : H2(?) ? Pt?1, where Pt?1 denotes the space of piecewise polynomial functions of degree ?t?1, then there exists a constant c such that bardblu?IubardblHm(Ki) ?cht?m|u|Ht(Ki) for u?H2(?), 0 ?m?t, t = 2. 11 1.3 Smooth Manifolds We are interested in finite element discretization of spheres, ellipsoids, cylinders and tori or geometries that we generally call spheroidal geometries. Spheres are smooth manifolds in R3, or genereally in Rn. In this section, we give a brief introduction to manifolds and define the Laplace-Beltrami operator on smooth manifolds. Before introducing smooth manifolds, let?s first define topological manifolds. Definition 1.1 A set X together with the collection of its subsets, J, form a topology if ? X,? ?J, ? union of any sets that are in J is also in J ? the intersection of any two sets that are in J is also in J. We call these sets that are in J the open sets. A special case of topological spaces is a Hausdorff space. It is a topological space with the property that for every two distinct points x,y in X, there exists an open neighborhood U of x and an open neighborhood V of y, which are distinct as well. The definition of k-dimensional topological manifold, M negationslash= ?, M ? Rn, is given in [1] as a Haussdorff, second countable topological space that is locally homeomorphic (one-to-one, onto, continuous and has a continuous inverse) to Rk. In other words, topological manifolds are the spaces that locally look like Euclidean spaces. Example The unit sphere, S2, in R3 is a 2-dimensional topological manifold. The unit sphere, S2, is the set of unit vectors in R3, that is S2 = {x = (x1,x2,x3), x?R3 : |x| = 1}. 12 The unit sphere, S2, is by definition a topological subspace of R3. Thus, it is Haussdorff and second countable. We need to show that it is locally Euclidean. To show that it is locally homeomorphic to R2, we need to show that for every point p in S2, there exists a homeomorphism, ?: U ? S2 ?V ? R2 between the open neighborhood U of p and a open subset V of R2. Let U+i denote the subset of S2 such that ith coordinate is positive, for i = 1,2,3. Similarly, let U?i denote the subset of S2 such that ith coordinate is negative, for i = 1,2,3, U+i = {(x1,x2,x3) ? S2 : xi > 0}i = 1,2,3, U?i = {(x1,x2,x3) ? S2 : xi < 0}i = 1,2,3. Let?U+ i : U+i ? R2 denote the homeomorphisms fromU+i to?U+ i (U+i ) ? R2, fori = 1,2,3, which are defined as ?U+ 1 (x1,x2,x3) = (x2,x3), ??1U+ 1 (x2,x3) = ( radicalBig 1?(x22 +x23),x2,x3), ?U+ 2 (x1,x2,x3) = (x1,x3), ??1U+ 2 (x1,x3) = (x1, radicalBig 1?(x21 +x23),x3), ?U+ 3 (x1,x2,x3) = (x1,x2), ??1U+ 3 (x1,x2) = (x1,x2, radicalBig 1?(x21 +x22)), ?U? 1 (x1,x2,x3) = (x2,x3), ??1U? 1 (x2,x3) = (? radicalBig 1?(x22 +x23),x2,x3), ?U? 2 (x1,x2,x3) = (x1,x3), ??1U? 2 (x1,x3) = (x1,? radicalBig 1?(x21 +x23),x3), ?U? 3 (x1,x2,x3) = (x1,x2), ??1U? 3 (x1,x2) = (x1,x2,? radicalBig 1?(x21 +x22)). Thus, the unit sphere in R3 is 2-dimensional topological manifold. 13 Let U be an open subset of M and let ? :U ??(U) ? Rk be a homeomorphism then we call the pair (U,?) a chart. If p is a point in U we say that the pair (U,?) is a chart at p, and we call U the coordinate neighborhood of p, and we call the Cartesian coordinates of p, ?(p) = (?1(p),?2(p),...,?k(p)) = (x1,x2,...,xk) ? Rk the coordinates of p. Definition 1.2 A function f :M ? R is continuous at p?M, if for some chart (U,?) at p, f ???1 : ?(U) ? R is continuous at ?(p). The function f : M ? R is continuous on the open set U ?M, if it is continuous at all points p in U. If M is a topological manifold, the notion of a continuous function f : M ? R makes sense, but the notion of a differentiable function does not, since differentiablity property of a function is not invariant under homeomorphism. We need to add more features to a topological manifold and get a smooth manifold. We say that f is smooth, if f ? C?(U), that is if f ? Cr(U) for every r. The function f is said to be a diffeomorphism, if f and its inverse f?1 are smooth. Definition 1.3 Let (U,?) and (V,?) be two charts, such that UintersectiontextV negationslash= ?, the charts (U,?) and (V,?) are called C?-compatible if the transition map ????1 : ?(UintersectiontextV) ??(UintersectiontextV) is a diffeomorphism. If we index the charts, (Ui,?i), i ? I where i ranges over some index set I, such that the sets Ui cover M, then we call such a collection of charts an atlas and denote it by A. An atlas A is called a differential structure if (U,?) is a chart in M such that it is compatible with every chart in M, then (U,?) is in A. Definition 1.4 The pair (M,A) is a smooth manifold, where M is a topological manifold and A is a differential structure [2]. 14 M U V ? ? ?(U) ? Rn ?(V) ? R n ????1 Figure 1.3: A transition map We usually omit the differential structure and say M is a smooth manifold. A smooth manifold is a topological manifold whose transitions maps are all smooth. For more detailed information on smooth manifolds see [2]. Example The unit sphere, S2, in R3 is a 2-dimensional smooth manifold. Since, it is clear from the previous example that {U+i ,U?i ,i = 1,2,3} is an open cover for S2. For simplicity, let?s consider, U+1 intersectiontextU?1 = {(x1,x2,x3) ? S2 : x1 > 0,x2 > 0}, then the transi- tion maps, ?U+ 1 ???1U+ 2 = (radicalbig1?(x21 +x22),x3) and ?U+ 2 ???1U+ 1 = (radicalbig1?(x22 +x23),x3) are compatible. Any other pair can be shown to be compatible, in a similar way. Furthermore, the unit sphere is a compact smooth manifold. Definition 1.5 A function f : M ? R is smooth (C?(M)) if and only if f is continuous and for every p in M, and there exists a chart (U,?) at p such that f ???1 : ?(U) ? R is smooth (C?((?(U)). 15 This definition is well defined, since it is independent of the choice of the the chart at p. We can generalize this definition easily to smooth maps between smooth manifolds. Let M negationslash= ?, M ? Rn, be k-dimensional smooth manifold and N negationslash= ?, N ? Rn be l-dimensional smooth manifold, the mapping F : M ? N is a smooth mapping, if and only if F is continuous and for every p ? M, there exists a chart (U,?) in M and a chart (V,?) in N such that F(U) ?V, and the mapping ??F ???1 : ?(U) ??(V) is smooth. Tangent Spaces Hereinafter M denotes k-dimensional smooth submanifold of Rn. We first, define the partial derivatives of a differentiable function f : M ? R with respect to the coordinate system (U,?). At this point, it is worth recalling the limit definition of partial derivatives of the function f : Rk ? R. Let x = (x1,...,xk) ? Rk, then we denote by Dif(x1,....xk), the number lim h?0 f(x1,...,xi +h,...,xk)?f(x1,....xk) h . For a function f : M ? R and local charts (U,?), we define Dif(p) = Diparenleftbigf ???1)(?(p)parenrightbig. (1.6) Thus, Dif(p) is the number lim h?0 fparenleftbig??11 (x),...,??1i (x) +h,...,??1k (x)parenrightbig?f(??1(x)) h . 16 If we define the curve ?i : (??,?) ? M by ?i(h) = (??11 (x),...,??1i (x) +h,...,??1k (x)), then (1.6) is just (f ??i)?(0) = lim h?0 f(?i(h))?f(?i(0)) h . Define the collection of all such C1 curves in M by Cp = {? ? C1 ((?,?),M),?> 0,?(0) = p}. We say that x? Rn is tangent to M at the point p, if and only if there exists a ? ? Cp such that x = ??(0). Definition 1.6 We define the tangent space at p as TpM = {(??(0)) :? ? Cp}, and tangent bundle TM as TM := uniondisplay p?M TpM. An element of TM is an ordered pair (p,x) with p ? M x ? TpM. The projection map ? : TM ? M maps each vector in TpM to the point at which it is tangent as ?(p,x) = p. Remark 1.1 For any p?M, TpM is a k-dimensional vector space, and if (U,?) is chart containing p, then (D1??1,D2??1,...,Dk??1) forms a basis for TpM. Remark 1.2 Tangent bundle TM is a 2k-dimensional manifold and the ? : TM ? M is a smooth map. Definition 1.7 The tangent map TF :TM ?TN for a C1-map F :M ?N is defined by TF(?(0),??(0)) = (?(0),(F ??)?(0)). 17 The differential of F, is a linear mapping DF :TpM ?TF(p)N which is defined as: F(??(0)) = (F ??)?(0). Define gp(?,?) :TpM ?TpM ? R as gp(x,y) = x?y, where x?TpM y ?TpM and x?y is the inner product TpM inherited from Rn. Recall that (D?1??1,D2??1,...Dk??1) is a basis for TpM, and define the metric tensor G = gij as gij = Di??1 ?Dj??1 and g = det(G) is called the Gramm determinant and is strictly positive. Let f : M ? R be a differentiable function, then the differential of f is a linear mapping TpM to R, thus it is a functional on TpM. So, by the Riesz Representation (Theorem 1.1), there exists a unique vector in TpM, call it gradient of f, and denote it by ?f, such that Df(p)v = ?f ?v for every v ?TpM. Remark 1.3 Let f :M ? R and (U,?) a chart containg p, then ?f(p) is defined as ?f(p) = ksummationdisplay i=1 ksummationdisplay j=1 gijDj(f ???1)(?(p)). (1.7) Definition 1.8 A vector field is a continuous map X : M ?TM, such that ??X = IdM, that is X(p) = (p,x(p)) ?p?M such that x(p) is a tangent vector to p. 18 Let X be C1-vector field on M with local representation X ?? = ksummationdisplay i=1 Xi?i?, define a scalar function ??X :M ? R as (??X)?? = 1?g ksummationdisplay i=1 ?iparenleftbigXiparenrightbig. For any C2 function f :M ? R define Laplace-Beltrami operator as ?f = ???f. Let ? be a local chart then the local representation is (?f)?? = (???f)?? = 1radicalbig|g| ksummationdisplay j=1 ?j( radicalbig |g| ksummationdisplay i=1 gij?i(f ??)) (1.8) where gij are the components of g?1. The Laplace Beltrami Operator on the Sphere Let ?(?1,?2) = (sin?1 cos?2,sin?1 sin?2,cos?1) 0 ? ?1 ? ?, 0 ? ?2 < 2? be a local chart of the sphere, then g = ? ?? ??1????1? ??1????2? ??2????1? ??2????2? ? ??= ? ?? 1 0 0 sin2?1 ? ?? and also, g?1 = ? ?? 1 0 0 1sin2 ? 1 ? ?? and |g| = sin2? 1. 19 We obtain the local representation of the Laplace-Beltrami ?Sf of any C2 function f on the sphere as (?Sf)?? = 1radicalbig|g| 2summationdisplay j=1 ?j( radicalbig |g| 2summationdisplay i=1 gij?i(f ??)) = 1sin? 1 [??1(sin?1??1)+??2( 1sin? 1 ??2)](f ??) (1.9) = 1sin? 1 ??1(sin?1??1) + 1sin2? 1 ??2(??2)f ??. 1.4 Eigenvalue Problem Let ?S be the Laplace operator on the sphere, then the eigenvalue problem is to find a scalar ? such that ??S? = ??. If the local representation of Laplace-Beltrami operator is obtained using polar coordinate charts (1.9), then the eigenvalue problem is given by 1 sin?1??1(sin?1??1?) + 1 sin2?1??2(??2?) = ???. (1.10) Writing ? in separable form ?(?1,?2) = ?(?1)?(?2), we obtain the following two eigenvalue problems d2? d?22 +??= 0 ?(0) =?(2?), d? d?2(0) = d? d?2(2?) (1.11) and sin?1 dd? 1 (sin?1 d?d? 1 ) +(?sin2?1 ??)? = 0, lim ?1?0 ? 0, c2 > 0 such that |J| ?c1 and |J?1| ? 1c2.Thus bardbl?bardbl2L2(K) = integraldisplay K |?|2 = integraldisplay Kh |??P|2|J| ? c1 integraldisplay Kh |v|2 = c1bardblvbardbl2L2(Kh). Similarly, bardblvbardbl2L2(Kh) ? 1c 2 bardbl?bardbl2L2(K). For any x ? Kh, there exists constants such that d1 ? 1bardblxbardbl ? d2, then by (3.14) we get d1|?aE?(a)| ? |?xE?(x)| ?d2|?aE?(a)|. Thus, integraldisplay K |?S?|2 = integraldisplay K |?E?|2 = integraldisplay Kh |?E??P)|2|J| ? c1d 1 integraldisplay Kh |?v|2. So, there exists some constants c3 > 0, and c4 > 0 such that c3bardblvbardblH1(Kh) ? bardbl?bardblH1(K) ?c4bardblvbardblH1(Kh). 47 Taking the gradient of both sides of (3.14), we get ?2xE?(x) = ?x parenleftbigg 1 bardblxbardbl parenrightbigg ?aE?(a) + 1bardblxbardbl2A?a (?aE?(a)). So for some constant c5 we get the following estimate, |v|H2(Kh) ?c5bardbl?bardblH2(K). square Let I denote the linear interpolation operator for continuous functions defined on S, then I?(a) = msummationdisplay i=1 ?(?i)?i(a) is uniquely determined by the interpolation conditions. De- note the restriction of I to the spherical triangle K by IT = I|T , then IT? = 3summationdisplay i=1 ?(?i)?i. Let Ih denote the linear interpolation operator for continuous functions defined on ?, then Ihu(x) = msummationdisplay i=1 u(?i)?i(x), is also uniquely determined by the interpolation conditions. Simi- larly, let IKh = Ih|Kh denote the restriction of Ih onto the planar triangle Kh. The linear interpolation operator Ih on ? and interpolation operator I on S are related to each other, in particular I = Ih ?P?1. Theorem 3.2 Let ??H2(S), then we have the following estimates bardbl???hbardblL2(K) ?h2bardbl?bardblH2(K), and bardbl???hbardblH1(K) ?hbardbl?bardblH2(K) 48 Proof: bardbl???hbardblH1(K) ? bardbl??IK?bardblH1(K), ? cbardblu?IKhubardblH1(Kh), ? ch|u|H2(Kh), ? chbardbl?bardblH2(K). Similarly, bardbl???hbardblL2(K) ? chbardbl???hbardblH1(K), ? ch2bardbl?bardblH2(K). square 49 Chapter 4 Numerical Experimetns 4.1 Finite elements on the sphere Consider the unit sphere, S of R3, centered at the origin and consider the box B in R3, with side length 2 which is also centered at the origin. Let {K1,...,Kn} denote a triangulation of the box and recall that the box is radially projected onto the sphere. Let {T1,...,Tn} denote the corresponding finite element discretization of the sphere. For simplicity, we use linear finite elements for the box, so?h denotes the finite element space of the box such that it consists of piecewise linear continuous polynomials: ?h = {v : B ? R s.t. v is a piece wise continuous linear polynomial }. Let {?i}mi=1 denote nodes of the finite element discretization of the box and let {?i}mi=1 denote the basis functions of the finite element space ?h with ?i(?j) = ?ij where ?ij is the Kronecker-delta. The support of the basis function ?i consists of all the triangles which share the node ?i. Since, ?i are linear basis functions, for any point x = (x1,x2,x3) on the box, we have ?i(x1,x2,x3) = kx1 +mx2 +nx3 +d for some constants k,m,n and d. For example, if x = (x1,x2,x3) is on F1, then ?i(1,x2,x3) = mx2 +nx3 +c Let ? denote the finite element space on the unit sphere. We denote the basis function of?by {?i}mi=1. Letabe any point on the sphere then there exists a pointxon the box with P(x) = a, we define the basis functions on the sphere as ?i = ?i?P?1. For example, let the inverse projection of a = (a1,a2,a3) be on the face F1, then ?i(a1,a2,a3) = ?i(1, a2a1, a3a1) = ma2a1 +na3a1 +c. So, the basis functions on the sphere are rational functions. Figure (4.1) shows a basis function on the sphere. 50 0 1 2 ?1 ?0.5 0 0.5 1 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 0.8 Figure 4.1: Basis functions on the sphere 4.2 Example 1 To illustrate our method and to verify the theoretical results, we present numerical experiments for the following model problem ??S?+? = f, on S. We choose a classical solution of the above problem to be ?(a1,a2,a3) = cosa1, then we calculate the right hand side, f by f = ??S ??S +?, ?S = ???(???n)n ??S ??S = ???S ? 3summationdisplay i=1 (?S,i ?n)ni (4.1) and get that f(a) = a41cos(a1) ? (1 + a21)a21cos(a1) + 2cos(a1) + 2a31sin(a1) ? (2 + 2a? 12)a1sin(a1). We denote by ?h the approximate solution, and by h the largest diameter of 51 the spherical triangles at the jth refinement step. The experimentally observed convergence rate, p is p = ln parenleftBig Ej Ej+1 parenrightBig ln parenleftBig hj hj+1 parenrightBig, where Ej denotes the error in the L2-norm at jth refinement step. Similarly, q is the experimentally observed convergence rate when the error is measured in the H1-norm. We present the results in the following table (4.2). j n h bardbl???hbardblL2(S) p bardbl???hbardblH1(S) q 1 48 0.9553 0.0439 - 0.3466 - 2 192 0.6155 0.0300 0.8661 0.2763 0.5157 3 768 0.3398 0.0087 2.0837 0.1478 1.0531 4 3072 0.1750 0.0023 2.0049 0.0757 1.0083 5 12288 0.0882 0.0005848 1.9986 0.0382 0.9982 6 49152 0.0441 0.000147 1.9987 0.0191 1.0033 7 196,608 0.0221 0.000036799 1.9981 0.0096 0.9925 Table 4.1: Numerical results. In this table j denotes the refinement step, h denotes the largest diameter of the spherical triangles measured using geodesic distance formula (2.5), n deontes the number of spherical triangles, p denotes the experimentally observed convergence rate in the L2-norm, and q denotes the experimentally observed convergence rate in the H1-norm . Figure (4.2) is a plot of the approximate solution ?h. 52 ?1 ?0.5 0 0.5 1 ?1 ?0.5 0 0.5 1?1 ?0.5 0 0.5 1 y x z 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Figure 4.2: Numerical approximation to ?(a) = cosa1. 4.3 Example 2 Next consider the following second order elliptic equation defined on the ellipsoid, whose equation is given by a21 + a224 + a234 = 1 ??E?+? = f, on E. We chose the exact solution to be ?(a) = cos(a1), then the right hand side f is calculated by using (4.1), where n = (2a21, 2a224 , 2a234 ). As in the previous example, ?h denotes the approximate solution, and h denotes the largest diameter of the finite element at the jth refinement step. We calculate the experimentally observed convergence rate and the results are presented in table (4.3). 53 j n h bardbl???hbardblL2(S) p bardbl???hbardblH1(S) q 1 48 0.9553 0.0439 - 0.3466 - 2 192 0.6155 0.0300 0.8661 0.2763 0.5157 3 768 0.3398 0.0087 2.0837 0.1478 1.0531 4 3072 0.1750 0.0023 2.0049 0.0757 1.0083 5 12288 0.0882 0.0005848 1.9986 0.0382 0.9982 6 49152 0.0441 0.000147 1.9987 0.0191 1.0033 7 196,608 0.0221 0.000036799 1.9981 0.0096 0.9925 Table 4.2: Numerical results 4.4 Example 3 Consider the following eigenvalue problem, ??S? = ??. (4.2) Eigenvalues of the problem (4.2) are ?k = k(k + 1) k = 0,1,... and the eigenvalue ?k has a multiplicity of 2k+ 1. Galerkin approximation of the variational problem is find ?h ? R and 0 negationslash= ?h ?? s.t. ah(?h,?h) = ?h(?h,?h) ??h ??. Here ah(?,?) is as in (1.15) and ? is as in (3.5). Exact eigenvalues 2,6,12,20 and their approximates with the finite element method is presented in table (4.4). Here n denotes the number of elements in each refinement step. 54 n 2 6 12 20 48 2.2325 6.5237 14.8632 29.3408 2.2325 6.5237 18.3974 30.1297 2.2325 9.015 18.3974 30.1297 9.015 18.3974 35.6114 9.015 22.2659 35.6114 22.2659 35.6114 22.2659 35.6114 51.311 51.311 51.311 192 2.0568 6.3224 13.437 24.192 2.0568 6.3224 13.656 24.192 2.0568 6.5956 13.656 24.385 6.5956 13.656 25.108 6.5956 14.362 25.108 14.362 25.108 14.362 26.483 26.483 26.483 768 2.0146 6.0912 12.416 21.248 2.0146 6.0912 12.444 21.248 2.0146 6.1491 12.444 21.341 6.1491 12.444 21.341 6.1491 12.581 21.341 12.581 21.371 12.581 21.502 21.502 21.502 3072 2.0037 6.0239 12.107 20.324 2.0037 6.0239 12.115 20.324 2.0037 6.0375 12.115 20.341 6.0375 12.115 20.341 6.0375 12.146 20.341 12.146 20.358 12.146 20.376 20.376 20.376 Table 4.3: Exact eigenvalues and their approximates 55 2 6 12 20 emin/emax pmin/pmax emin/emax pmin/ pmax emin/emax pmin/pmax emin/emax pmin/pmax 0.2325/ 0.2325 0.5237 / 3.015 2.8632/10.266 9.3408/31.311 0.0568/0.0568 3.2061/3.2061 0.3224/0.5956 1.1036/3.6893 1.437/2.362 1.5682/3.3425 4.192 /6.483 1.8226/3.5824 0.0146/ 0.0146 2.2868/2.2868 0.0912/0.1491 2.1255 / 2.3313 0.416 /0.581 2.0866/ 2.3608 1.248/1.502 2.0395/2.4616 0.0037/ 0.0037 2.0686/2.0686 0.0239/0.0375 2.0181/2.0801 0.107/0.146 2.0463/2.0814 0.324/0.376 2.0323/2.0871 Table 4.4: Experimentally calculated convergence rates. 56 Chapter 5 Conclusion We developed and analyzed a new method for exactly discretizing spheroidal domains and constructing finite element spaces on such domains, thus yielding conforming finite element discretizations. This method may be used to approximate solutions of partial differential equations, as well as eigenvalues and eigenfunctions of differential operators de- fined on such domains (spherical, ellipsoidal, cylindrical, and toroidal shells). This method which was described for the Laplace-Beltrami equation can be easily generalized to approx- imate solutions of other partial differential equations and eigenvalue problems defined on the spheroidal domains. The methods is easy to implement and can be easily incorporated into existing finite element codes. Moreover, the method can be generalized to other convex domains by means of generalized projections, provided these are available analytically. 57 Bibliography [1] Denis Barden and Charles Thomas, ?An Introduction to Differential Manifolds?, Im- perial College Press, 2005. [2] John M. Lee, ?Introduction to Smooth Manifolds?, Springer, 2002. [3] J. Wloka, ?Partial differential equations?, Cambridge University Press, 1987. [4] P.E. Ciarlet, ?The Finite Element Method for Elliptic Problems,? SIAM, 2002. [5] Gilbert Strang, George Fix, ?An Analysis of the Finite Element Method ,? Prentice- Hall, 1973. [6] Ronchi, C.; Paolucci, P. S., The ?cubed sphere?: a new method for the solution of partial differential equations in spherical geometry, J. Comput. Phys., 124(1996), 93- 114. [7] Qiang Du, Lili Ju, ?Finite volume methods on spheres and spherical centroidal voronoi meshes?, SIAM Journal on Numerical Analysis, 43, (2005),1673-1692. [8] D. Gilbarg, N.S. Trudinger ? Elliptic partial differential equations of second or- der?,Springer, 1977. [9] Gerhard Dziuk, Finite elements for the Beltrami operator on arbitrary surfaces.Partial differential equations and calculus of variations, 142-155, Lecture Notes in Math., 1357, Springer, Berlin, 1988. [10] John R. Baumgardner; Paul O. Frederickson, Icosahedral Discretization of the Two- Sphere, SIAM Journal on Numerical Analysis, 22 (1985) 1107-1115. [11] Dietrich Braess, ?Finite elements. Theory, fast solvers, and applications in solid me- chanics?, Cambridge, 2001. [12] Susanne Brenner, L. Ridgway Scott ? The mathematical theory of finite element meth- ods?, Springer, Verlag, 1994. [13] Andrew V. Knyazev, John E. Osborn New a priori FEM error estimates for eigenvalues, SIAM Journal on Numerical Analysis, 43(2006),2647-2667. [14] I. Babuska, J.E. Osborn Estimated for the errors in eigenvalue and eigenvector ap- proximation by Galerikin methods, with particular attention to the case of multiple eigenvalues. SIAM Journal on Numerical Analysis, 24(1987)1249-1276. 58 [15] I. Babuska, J.E. Osborn Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems. Mathematics of Computation,52(1989)275-297. [16] Thierry Aubin Nonlinear Analysis on manifolds. Monge-Amp`ere Equations, Springer- Verlag, 1982. 59