A numerical domain decomposition method for solving elliptic equations on manifolds

Elliptic partial differential equations on Riemannian manifolds are of fundamental importance both in analysis and geometry (see e.g., [50, 32]). A simple and important example is

Here $\Delta$ is the Laplace-Beltrami operator, or Laplacian for brevity, defined on a $d$-dimensional Riemannian manifold $M$. Many manifolds are naturally submanifolds of Euclidean spaces, here the “dimension” of a submanifold is referred to as the topological dimension of the manifold, not the dimension of its ambient Euclidean space. For example, the $n$-dimensional unit sphere $S^{n}$ is embedded in $\mathbb{R}^{n+1}$.

When the manifold $M$ is a two-dimensional Riemannian submanifold of $\mathbb{R}^{3}$, i.e. a surface, the numerical methods to solve PDEs, particularly (1), on $M$ have been extensively studied for a long history (see e.g. [42, 41, 6, 19, 20, 22, 44]). Over several decades, among many others, the surface finite element method and its variant have had far-reaching developments (see e.g., [2, 4, 15, 16, 27, 30, 43, 47]), and applications to various PDEs (see e.g. [3, 7, 10, 21, 23, 31]); see also e.g. [14, 22, 8] for surveys and bibliographies, and e.g. [1, 5, 26] for software developments. They have also been widely applied to various areas such as computer graphics (e.g. [25]), surface fitting (e.g. [18]), shape analysis (e.g. [17, 48, 49]), isogeometric topology optimization (e.g. [33, 34]), and medical imaging (e.g. [38]). A basic and common feature of these finite element methods is to construct a global triangulation of $M$, which provides a necessary grid structure for a global finite element space. This triangulation, following [39, Definition 8.3] in the sense of differential topology, is a bijection between a simplicial complex and the manifold $M$, which satisfies certain regularity. In computation, a concrete representation of this triangulation is to approximate $M$ by a polyhedron in $\mathbb{R}^{3}$. Alternatively, $M$ is represented implicitly as a level set of a function $\phi$, say $M=\phi^{-1}(0)$, and an equation on $M$ can be solved by the methods of trace elements or implicit surfaces elements (e.g. [43, 47, 22]).

However, there are several essential difficulties to apply the methods above to general higher dimensional manifolds. First, many important and interesting examples of higher dimensional manifolds are not submanifolds of Euclidean spaces by definition. A notable example is the complex projective spaces $\mathbb{CP}^{n}$, which serve as the foundation for algebraic geometry. Though Whitney Embedding Theorem ([54], see also [29, p. 24-27]) and Nash Embedding Theorem ([40]) do reveal that every smooth manifold can be embedded into a Euclidean space $\mathbb{R}^{k}$ differential-topologically or even geometrically for some large $k$. To the best of our knowledge, there is no literature on a universal algorithm to efficiently construct such an explicit embedding for computation on general higher dimensional manifolds. Second, assuming $M$ is a submanifold of $\mathbb{R}^{k}$, if the codimension of the submanifold $M$ in $\mathbb{R}^{k}$ is greater than $1$, which is often the case, it will be horribly difficult to find an effective polytopal approximation to $M$ in general due to topological and geometrical complexity, meanwhile, $M$ also cannot be represent as a level set of a function.

A global triangulation of a manifold is helpful in numerical computation, since it can provide a global discretization of the target problem. We do acknowledge that J. H. C. Whitehead proved ([53], see also [39, Theorem 10.6]) that every smooth manifold can be globally triangulated in an abstract way. However, to the best of our knowledge, in practice, there has been no algorithm to build such a concrete triangulation or a grid over a high-dimensional manifold in general. To circumvent this difficulty, in [46] which serves as our major inspiration, Qin–Zhang–Zhang proposed a new idea to numerically solve elliptic PDEs on manifolds by avoiding global triangulations completely. Since a $d$-dimensional manifold $M$ has local coordinate charts by definition, $M$ can be decomposed into finitely many subdomains that carry local Cartesian coordinates. Consequently, an elliptic equation on each subdomain can be transformed to one on a domain in $\mathbb{R}^{d}$. Thus an elliptic problem on $M$ can be assembled to either coupled problems on Euclidean domains, or can be solved directly by Domain Decomposition Methods (DDMs). This idea had been numerically verified on the $3$-dimensional unit sphere $S^{3}$ in [46], where $S^{3}$ is decomposed into two subdomains. However, a major drawback of [46] is its lack of flexibility to deal with more general manifolds, whose charts may involve more than two subdomains.

In this paper, we shall develop the idea in [46] further to solve problems on general manifolds. Similar to [46], global triangulations or grids shall be completely avoided. In fact, we shall solve such problems by an overlapping domain decomposition method (DDM). The development of DDM has a long history. It was first invented by H. A. Schwarz [51]. The version of DDM we mainly follow was originally proposed by P. L. Lions in [36, I. 4] for solving continuous problems in Euclidean spaces. It has been well-developed and is later named as Multiplicative Schwarz Method (c.f. [52]). The later development usually takes this DDM as a preconditioner for a globally discretized problem (see e.g. [9, 24], and pure algebraic versions in e.g. [11, 35]). However, a globally discretized problem should be based on a global grid on a manifold $M$, which is not accessible in our case. Therefore, we shall more closely follow Lions’ original approach rather than the later development. This original approach is a very simple iteration scheme such that a local problem in a subdomain is solved in each step. We found that this DDM can be well-adapted to solve problems on manifolds. As in [46], a problem in each subdomain can be converted to one in a Euclidean domain. It is simpler to obtain a grid and hence a discretization over each Euclidean domain. To minimize the difficulty of coding, these Euclidean domains can be even chosen as (high-dimensional) rectangles. The transition of information among these Euclidean domains is by virtue of the transition maps of coordinates.

Overall, our approach is a discrete imitation of Lions’ method with additional transition maps of coordinates. It is necessary to point out that the transition maps in the proposed method are not required to preserve nodes or grids. More precisely, in computations, suppose there are two Euclidean domains $D_{i}$ and $D_{j}$. Both $D_{i}$ and $D_{j}$ carry a grid, respectively. If $\phi_{ji}:D_{i}\rightarrow D_{j}$ is such a transition map between $D_{i}$ and $D_{j}$, the image of the grid on $D_{i}$ under $\phi_{ji}$ may not match the grid on $D_{j}$. Furthermore, for a node $\xi\in D_{i}$, $\phi_{ji}(\xi)$ may not be a node anymore on $D_{j}$ (see Fig. 1 for an example). Nevertheless, this incompatibility shows the high flexibility of our approach. In fact, if all transition maps preserved grids, we would obtain a global grid on $M$ which is practically almost impossible. The transition maps in [46, Section 4] do not preserve grids but do preserve boundary nodes on each subdomain, i.e., the $\phi_{ji}$ above maps nodes on $\partial D_{i}$ to nodes of $D_{j}$. As a result, the methods therein have some flexibility to solve PDEs on $S^{3}$. In order to handle more general high-dimensional manifolds, this paper improves the flexibility further. Particularly, interpolation techniques are employed now to handle nonmatching grids.

The proposed approach is numerically verified on closed manifolds of dimension $4$. They are $4$-dimensional unit sphere $S^{4}$, the complex projective space $\mathbb{CP}^{2}$, and the product manifold $S^{2}\times S^{2}$. The numerical results show that our method solves equation (1) on these manifolds in a natural manner.

The outline of this paper is as follows. In Section 2, we shall extend P. L. Lions’ method in [36] to a DDM for continuous problems on manifolds and prove its convergence. In Section 3, we shall propose our discrete imitation of Lions’ method. Some specialty of product manifolds will be explained in Section 4. Finally, some numerical results will be presented in Section 5.

In this section, we shall first formulate a second order elliptic model problem on general manifolds. Then we introduce a domain decomposition method (DDM) generalized from [36, I. 4] to solve the model problem (see Algorithm 1 below). We shall formulate and prove Theorem 2.1 below on the convergence of this iterative procedure. This procedure also motivates us to propose Algorithm 2 in next section, which gives numerical approximations to the solution to the model problem.

Let $M$ be a $d$-dimensional compact smooth manifold without or with boundary $\partial M$. Equipping $M$ with a Riemannian metric $g$, the Laplacian $\Delta$ can be defined on $M$. In general, neither $g$ nor $\Delta$ can be expressed by coordinates globally because $M$ does not necessarily have a global coordinate chart. In a local chart with coordinates $(x_{1},\dots,x_{d})$, the Riemannian metric tensor $g$ is expressed as

where the matrix $(g_{\alpha\beta})_{d\times d}$ is symmetric and positive definite. The Laplacian $\Delta$ can be then expressed in this chart as

where $G=\det\left((g_{\alpha\beta})_{d\times d}\right)$ is the determinant of the matrix $(g_{\alpha\beta})_{d\times d}$ and $(g^{\alpha\beta})_{d\times d}$ is the inverse of $(g_{\alpha\beta})_{d\times d}$. It is well-known that $\Delta$ is an elliptic differential operator of second order.

We consider the following model problem on $M$

where $b\geq 0$ is a constant and $f\in L^{2}(M)$. A weak solution to (3) is a solution to the following problem: Find a $u\in H^{1}_{0}(M)$ such that, $\forall v\in H^{1}_{0}(M)$,

Here $\nabla u$ and $\nabla v$ are the gradients of $u$ and $v$ with respect to $g$, $\langle\nabla u,\nabla v\rangle$ is the inner product of $\nabla u$ and $\nabla v$, and $\mathrm{dvol}$ is the volume form. In terms of local coordinates,

and

For brevity, we shall omit the symbol $\mathrm{dvol}$ in integrals on $M$ throughout this paper.

Note that it suffices to solve (4) on each component of $M$. Therefore, without loss of generality, $M$ is assumed to be connected. In addition, if $\partial M=\emptyset$, the condition above $u|_{\partial M}=0$ is vacuously satisfied in regards to $H^{1}_{0}(M)=H^{1}(M)$. To guarantee (4) is well-posed, we further assume $b>0$ if $\partial M=\emptyset$. (Actually, if $b=0$, one may imposed additional conditions such as $\int_{M}u\mathrm{dvol}=0$ to guarantee the well-posedness. On the other hand, the numerical algorithm would be more complicated in this situation. We shall study the case later.)

Now we describe a domain decomposition iterative procedure to solve (4). This method was originally proposed by P. L. Lions in [36, I. 4], in which the classical Schwarz Alternating Method is extended to an iterative procedure with many subdomains. The manifold nature of this method is intrinsically adapted to solve PDEs on manifolds. More precisely, suppose $M$ is decomposed into $m$ subdomains, i.e.

Here $M_{i}$ is a closed subdomain (submanifold with codimension $0$) of $M$ with Lipschitz boundary, and $\mathrm{Int}M_{i}$ is the interior of $M_{i}$ in the sense of point-set topology of $M$. Clearly, an element in $H^{1}_{0}(M_{i})$ can be naturally considered as an element in $H^{1}_{0}(M)$ by zero extension. Thus, $H^{1}_{0}(M_{i})\subseteq H^{1}_{0}(M)$. The iterative procedure to solve (4) is the following Algorithm 1.

To obtain the $u^{n,i}$ in (7), one first solve an elliptic problem on $M_{i}$ with Dirichlet boundary condition $u^{n,i}|_{\partial M_{i}}=u^{n,i-1}|_{\partial M_{i}}$, then extend the solution to a function $u^{n,i}$ on $M$ by defining $u^{n,i}|_{M\setminus M_{i}}=u^{n,i-1}|_{M\setminus M_{i}}$. Thus the $u^{n,i}$ is well-defined and hence Algorithm 1 is well-posed.

We have the following theorem on the geometrical convergence of Algorithm 1.

Theorem 2.1 was originally proved by P. L. Lions ([36, Theorem. I.2]) in the case that $M$ is a Euclidean domain. We shall adapt his proof to the case of manifolds.

An elementary proof of the following lemma can be found in [36, p. 17]. It is also an immediate corollary of [56, (1.2)].

Define the energy bilinear form on $H^{1}_{0}(M)$ as

In this section, we propose a numerical DDM iterative procedure (Algorithm 2 below) to obtain approximations to the solution of (4). The procedure is as follow. First, $M$ is decomposed into overlapping subdomains $M_{i}$ ($1\leq i\leq m$), and each $M_{i}$ is in a coordinate chart. Second, a DDM iterative procedure, which serves as a discrete counterpart of Algorithm 1 is applied. Due to $M_{i}$ being a coordinate chart, an elliptic problem on $M_{i}$ can be naturally converted to one on a domain in a Euclidean space. Then, this problem on $M_{i}$ can be solved approximately using conventional finite element methods. The transition of information among subdomains is by interpolation.

For the purpose of presentation, only manifolds without boundary are considered in the numerical examples. A forthcoming work will study in detail the numerical implementation on manifolds with boundaries. In that more general case, we shall have to apply some special technique to deal with the boundary.

Suppose $M$ and $M^{\prime}$ are compact manifolds with dimensions $d$ and $d^{\prime}$ respectively. The Cartesian product $M\times M^{\prime}$ is a compact manifold with dimension $d+d^{\prime}$. A decomposition of $M$ and another one of $M^{\prime}$ canonically result in a decomposition of $M\times M^{\prime}$. Actually, suppose $M=\bigcup_{i=1}^{m}\mathrm{Int}M_{i}$ and $M^{\prime}=\bigcup_{i^{\prime}=1}^{m^{\prime}}\mathrm{Int}M^{\prime}_{i^{\prime}}$, where $M_{i}$ (resp. $M^{\prime}_{i^{\prime}}$) are subdomains of $M$ (resp. $M^{\prime}$). Then

where $M\times M^{\prime}$ is decomposed into $mm^{\prime}$ subdomains $M_{i}\times M^{\prime}_{i^{\prime}}$ for $1\leq i\leq m$ and $1\leq i^{\prime}\leq m^{\prime}$.

This canonical decomposition of $M\times M^{\prime}$ reflects another advantage of the spaces of rectangular finite elements. More precisely, suppose there are diffeomorphisms $\phi_{i}:D_{i}\rightarrow M_{i}$ and $\phi^{\prime}_{i^{\prime}}:D^{\prime}_{i^{\prime}}\rightarrow M^{\prime}_{i^{\prime}}$, where each $D_{i}\subset\mathbb{R}^{d}$ (resp. $D^{\prime}_{i^{\prime}}\subset\mathbb{R}^{d^{\prime}}$) is a $d$-rectangle (resp. $d^{\prime}$-rectangle). Then we have the diffeomorphisms

where each $D_{i}\times D^{\prime}_{i^{\prime}}\subset\mathbb{R}^{d+d^{\prime}}$ is a $(d+d^{\prime})$-rectangle. The transition of coordinates between $D_{i}\times D^{\prime}_{i^{\prime}}$ and $D_{j}\times D^{\prime}_{j^{\prime}}$ is

If rectangular grids are created over both $D_{i}$ and $D^{\prime}_{i^{\prime}}$, a rectangular grid over $D_{i}\times D^{\prime}_{i^{\prime}}$ follows automatically.

In summary, the procedures of the decomposition and discretization of factor manifolds are helpful for those of product manifolds.

We perform several numerical tests of the proposed method on manifolds $S^{4}$, $\mathbb{CP}^{2}$ and $S^{2}\times S^{2}$. They are $4$-dimensional compact manifolds without boundary. While the proposed method applies to problems in all dimensions, the sizes of the linear systems derived from subdomains will increase exponentially with respect to the dimension. On one hand, we would have trouble in the storage of data. On the other hand, we would struggle to find solutions to these linear systems with desired accuracy. This difficulty is so called “the curse of dimensionality”. It is actually a typical phenomenon of Euclidean spaces rather than of general manifolds. For the sake of presentation and due to the constraint of computing resources, the numerical examples in this paper consider manifolds with dimension no more than $4$. The numerical challenge in higher dimensions will be tackled in a forthcoming future work.