Discontinuous Galerkin Isogeometric Analysis of Convection Problem on Surface

: The objective of this work is to study ﬁnite element methods for approximating the solution of convection equations on surfaces embedded in R 3 . We propose the discontinuous Galerkin (DG) isogeometric analysis (IgA) formulation to solve convection problems on implicitly deﬁned surfaces. Three numerical experiments shows that the numerical scheme converges with the optimal convergence order.


Introduction
Surface partial differential equations (SPDEs) arise in natural sciences and applied areas, such as minimal surface equation, Willmore flow, transport of surfactants along interfaces in multiphase fluids [1], and lipid interactions in cell membranes [2]. In this paper, we consider the DG-IgA methods for the following model problem: where ∇ Γ · = div Γ denotes the surface divergence, Γ is a smooth two-dimensional surface without boundary embedded in R 3 . β = (b 1 , b 2 , b 3 ) is the advective velocity with b i ∈ W 1 ∞ (Γ), γ ∈ L ∞ (Γ) is reaction coefficient. We call this model problem the convection equation on surface Γ. The problem (1) has a unique solution u ∈ H(div, Γ) satisfying Γ udσ = 0, where dσ is surface measure. Moreover, for source term f ∈ L 2 (Γ) we have Γ f dσ = 0.
To ensure the existence and uniqueness of a solution u ∈ H 1 (Γ) to the model problem (1.1), we adopt the following hypothesis: There exists a positive constant γ 0 such that Finite element methods of SPDEs have been extensively developed. G. Dziuk [3] first used the finite element method for solving the surface elliptic problem. A. Bonito and J. E. Pasciak [4] designed and analyzed multigrid algorithms for the Laplace-Beltrami operator on a smooth and closed surface. DG methods for solving the first-order hyperbolic problems were defined in [5][6][7][8]. Since DG methods have well-known stability and local conservation in numerical application of partial differential equation problems [9,10] and capture solution blow-ups, it is natural to extend the DG approximation to PDEs on surfaces or manifolds. Ref. [11] extended the discontinuous Galerkin (DG) methods to the first-order hyperbolic and advection-dominated problems on surfaces.
IgA has been introduced in [12] as a new tool for solving numerically PDEs with the complicated geometric domains, in particular, surfaces. DG-IgA methods have been applied for approximating solutions of elliptic problems on surfaces [13,14]. Ref. [15] developed DG-IgA numerical schemes for solving problems on segmentations with gaps. Ref. [16] studied NURBS-based isogeometric analysis for the computation of flows about rotating components. In [17], a new stabilized symmetric Nitsche method was proposed for enforcement of Dirichlet boundary conditions for elliptic problems in Cut-IgA. The remainder of the paper is organized as follows: First, we introduce some preliminaries about IgA and discrete NURBS finite space V(h). Then we derive the DG-IgA scheme of convection problem (1.1). Finally, we present three numerical experiments to illustrate the discrete formulation.

Differential Operators on Surfaces
Let us consider a two-dimensional surface Γ ⊂ R 3 defined in the physical space R 3 . Assume that the surface is characterized by a geometrical mapping from a parameter space Γ ⊂ R 2 . Let X ∈ C 2 ( Γ, R 3 ) be a local parameterization of θ ∈ Γ, where the vector-valued independent variable θ = (θ 1 , θ 2 ) is called parametric coordinate. Indeed, X defines the geometrical mapping as follows: To define surface gradient, we introduce the Jacobi matrix J(θ) = (J ij ) 3×2 as So, the metric tensor G(θ) = (g ij ) 2×2 is represented by The inverse of G(θ) is denoted by G −1 (θ) = (g ij ) 2×2 . The tangential or surface gradient is given by So, we can deduce the tangential or surface divergence. For υ is a C 1 vector field on Γ, we have Remark that the vector field υ may not be tangential to surface Γ. However, the 2-dimensional vector G −1 J T υ is the projection of υ in parametric coordinate θ, which is a tangent vector field of Γ. So, the definition of surface divergence is well-posed.
Let L 2 (Γ) denote the usual L 2 −space on the surface Γ with norm Furthermore, we use standard notations W m p (Γ) and H m (Γ) for Sobolev spaces on Γ with norm and semi-norm To present the weak form of problem (1), we introduce the Stokes theorem on surface Γ. Let υ be a C 1 vector field on surface Γ, q ∈ C 1 (Γ), then the following identity holds:

Isogeometric Analysis
To apply the IgA methodology of the problem (1), the domain Γ is partitioned into some closed sub-patches such that We denote the set of sub-patches as . Without loss of generality, we simply assume a parametric domain D of unit length, i.e., D = [0, 1] 2 . For each patch Γ i , we associate the knot vectors V i on D, which defines a partition T i . Any patch Γ i can be represented by a parametrization map as follows The set of all the edges of partition T h is denoted by E as follows, We denote the faces of all patches as F defined by Let h K and h e denote the size of element K ∈ T h and the length of the face e ∈ E respectively.
We assume that the shape of the elements is regular and quasi-uniform, i.e., h e ∼ h and h K ∼ h Γ i .
Next, we define the broken space on the physical domain Γ associated with T h by using the introduced push-forward function Now we can define the broken Sobolev space H k (T h ): To apply DG method to problem (1), where n i denotes the unit normal vector of Γ i on F ij pointing exterior to Γ j .

DG-IgA Discretization
Next, we introduce the finite element space associated with the partition T h . In general, the basis function B The DG-IgA approximation is formulated as follows: where where η e is the stability parameter, defined as |β · n|, otherwise, The bilinear map a h (u h , v h ) can be modified. Indeed, by using the continuity of β and the function of V h in each patch, i.e., Similarly, we can define the outflow faces F + . So, the bilinear map a h (u h , v h ) can be modified as follows We can derive the discrete coercivity, stability and consistency of above numerical scheme with a similar technique in [19]. Due to the discontinuity of β, it needs a few skills to prove the discrete coercivity. So there exists a unique discrete solution u h ∈ V h satisfying DG-IgA scheme (5). Here we briefly show the a priori error result without detailed proof. We introduce the dual or adjoint weak form: for all v ∈ V(h). With the above adjoint weak form, inverse inequality and the approximation of the interpolant operator, we can prove the following result. (1) and (5), respectively. Assume that u| K ∈ H k+1 (K), ω ∈ H k+1 (K), ∀K ∈ T h . Then we have u − u h L 2 (Γ) ≤ Ch k+1 u H k+1 (Γ) .

Numerical Experiments
In this section, we present some numerical experiments of convection problems on surface. Numerical examples are presented for a sphere and a quarter of a cylinder.
We divide the unit sphere Γ into 6 patches. For each patch, the knot vector is taken as [0, 0, 0, 1, 1, 1] × [0, 0, 0, 1, 1, 1] to represent the geometry of each patch. We generate the mesh by refining the parameter element of each patch, whose mesh size is denoted by h. We show the patches and the uniform meshes of the sphere for h = 1/4 in Figure 1(Left). The numerical L 2 errors and convergence results are given in Table 1 for k = 1, 2 and 3, respectively. In Figure 2, we present the convergence histories of L 2 errors. These results show that the IgA-DG method yields O(h k+1 ) convergent solution. We present the numerical solution of convection problem (1) on the sphere in Figure 1(Right).

Numerical Experiment 2
Here we continue to consider the model problem on the surface of torus subject to the compatibility condition Γ f dσ = 0, where β = [1, 1, 1] T and γ = 1. The torus is the surface with r = 1 and R = 2. We take coordinates (φ, θ) as and select the source function f such that the exact solution is u = sin(θ) sin(3φ). We divide the torus into 8 patches. For each patch, we take the knot vector as [0, 0, 0, 1/2, 1/2, 1, 1, 1] × [0, 0, 0, 1, 1, 1] to give the geometrical representation. We plot the patches and the uniform meshes of the torus for h = 1/4 in Figure 3(Left). The numerical L 2 errors and convergence rates of this problem for k = 1, 2 and 3 are shown in Table 2. Table 3 indicate that the rates are also O(h k+1 ) for L 2 norm. Figure 4 shows the convergence history of errors. Finally, we plot the numerical solution of convection problem (1) on torus in Figure 3(Right).

Numerical Experiment 3
Next, we solve on the surface of a square of the cylinder the model problem where γ = 1. The domain Γ is the surface of a quarter of the cylinder, shown in Figure 5. In contrast to the case of sphere, this problem needs the boundary condition Γ − determined by β. Then, writing (r = 1, θ, z) to denote the system of cylindrical coordinate, we impose an appropriate boundary condition g for u and source function f so that the exact solution is u = sin(zθ).
We consider the model problem on the surface of a square of the cylinder with continuous and discontinuous coefficient β.
(1) We take continuous coefficient as β = [−1, 1, 1] T . We divide the cylinder into 9 patches. For each patch, we take the knot vector as [0, 0, 0, 1, 1, 1] × [0, 0, 1, 1] to give the geometrical representation. Similarly, we plot the patches and the uniform meshes of the cylinder for h = 1/4 in Figure 5(Left).  We present numerical L 2 errors and convergence rates of this problem for k = 1, 2 and 3 in Table 3. Table 3 indicate that the rates are also O(h k+1 ) for L 2 norm. Figure 6 shows the convergence history of errors. Finally, we plot the numerical solution of convection problem (1) on cylinder in Figure 5(Right). (2) We consider the discontinuous advective velocity β as where the index (i) is the patch number. Observe that the source function is discontinuous across the patches according to the choice of β.
According to the discontinuity of β, we divide the cylinder into 4 patches. For each patch, we take the knot vector as [0, 0, 0, 1, 1, 1] × [0, 0, 1, 1] to give the geometrical representation. Similarly, we plot the patches and the uniform meshes of the cylinder for h = 1/4 in Figure 7(Left), where patch 1 and patch 4 are painted blue, and patch 2, 3 are painted yellow. The numerical L 2 errors and convergence rates of this problem for k = 1, 2 and 3 are shown in Table 4. Table 4 indicate that the rates are also O(h k+1 ) for L 2 norm. Figure 8 shows the convergence history of errors. Finally, we plot the numerical solution of convection problem (1) on cylinder in Figure 7(Right).

Conclusions
In this paper, we present the new penalty discontinuous Galerkin (DG) isogeometric analysis(IgA) methods to solve convection problems with continuous or discontinuous coefficient on implicitly defined surfaces. For further purpose, it is worthy studying the stability and error analysis of this method and more practical problems on surfaces.

Conflicts of Interest:
The authors declare no conflict of interest.