Discontinuous Galerkin Isogeometric Analysis for Elliptic Problems with Discontinuous Coefficients on Surfaces

TL;DR

This paper extends discontinuous Galerkin isogeometric analysis to handle elliptic surface problems with discontinuous coefficients and non-matching meshes, providing theoretical error estimates and numerical validation.

Contribution

It generalizes a priori error estimates to non-matching meshes and discontinuous coefficients on surfaces, advancing the applicability of dGIGA in practical scenarios.

Findings

Error estimates are validated through numerical experiments.

The method effectively handles discontinuities across patch interfaces.

The approach accommodates non-matching meshes in surface problems.

Abstract

This paper is concerned with using discontinuous Galerkin isogeometric analysis (dGIGA) as a numerical treatment of Diffusion problems on orientable surfaces $\Omega \subset \mathbb{R}^3$. The computational domain or surface considered consist of several non-overlapping sub-domains or patches which are coupled via an interior penalty scheme. In Langer and Moore U. Langer and S. E. Moore,2014, we presented a priori error estimate for conforming computational domains with matching meshes across patch interface and a constant diffusion coefficient. However, in this article, we generalize the \textit{a priori} error estimate to non-matching meshes and discontinuous diffusion coefficients across patch interfaces commonly occurring in industry. We construct B-Spline or NURBS approximation spaces which are discontinuous across patch interfaces. We present \textit{a priori} error estimate for the symmetric discontinuous Galerkin scheme and numerical experiments to confirm the theory.

Equations

$$-{\text{div}}_{\Omega }(\alpha {\nabla }_{\Omega }u)+u=f\text{in}\Omega ,u=0\text{on}{\Gamma }_D,\mathbf{n}\cdot (\alpha {\nabla }_{\Omega }u)=g_N\text{on}{\Gamma }_N,$$

$$a(u,v)$$

$$a(u,v)={\int }_{\Omega }\alpha {\nabla }_{\Omega }u\cdot {\nabla }_{\Omega }v+uvdx\text{and}\ell (v)={\int }_{\Omega }fvdx+{\int }_{{\Gamma }_N}{g}_{N}vds,$$

$$\Vert u{\Vert }_{H^2(\Omega )}\le \Vert f{\Vert }_{L_2(\Omega )},$$

$$\mathrm{\Phi }:\hat{\Omega }\to \Omega \subset {\mathbb{R}}^3,\xi \to x=\mathrm{\Phi }(\xi ),$$

$$\hat{J}$$

$$\hat{F}(\xi )$$

$$\varphi (x)=\hat{\varphi }(\xi )\circ {\Phi }^{-1}(x),x\in \Omega ,$$

$${\nabla }_{\Omega }\varphi (x):=\hat{J}(\xi ){\hat{F}}^{-1}(\xi )\nabla \hat{\varphi }(\xi )\circ {\Phi }^{-1}(x).$$

$${\nabla }_{\Omega }\cdot \varphi (x):=\frac{1}{\hat{g}\xi )}\nabla \cdot \left[\hat{g}(\xi ){\hat{F}}^{-1}(\xi ){\hat{J}}^T(\xi )\hat{\varphi }(\xi )\right]\circ {\Phi }^{-1}(x).$$

$${\Delta }_{\Omega }\varphi (x)=\frac{1}{\hat{g}(\xi )}\nabla \cdot \left[\hat{g}(\xi ){\hat{F}}^{-1}(\xi )\nabla \hat{\varphi }(\xi )\right]\circ {\Phi }^{-1}(x).$$

$$\hat{\mathbf{n}}(\xi ):=\frac{{\hat{t}}_1(\xi )\times {\hat{t}}_2(\xi )}{\Vert {\hat{t}}_1(\xi )\times {\hat{t}}_2(\xi )\Vert },$$

$$J(x)=\hat{J}(\xi )\circ {\mathrm{\Phi }}^{-1}(x),F(x)=\hat{F}(\xi )\circ {\mathrm{\Phi }}^{-1}(x)\text{and}g(x)=\hat{g}(\xi )\circ {\mathrm{\Phi }}^{-1}(x).$$

$${\hat{B}}_{i,0}(\xi )$$

$${\hat{B}}_{i,p}(\xi )$$

$${\hat{B}}_{\mathbf{i},\mathbf{p}}(\xi ):=\prod _{k=1}^2{\hat{B}}_{i_k,p_k}({\xi }_k),$$

$$\mathbf{\Phi }(\xi ):=\sum _{\mathbf{i}\in \overline{\mathcal{I}}}{C}_{\mathbf{i}}{\hat{B}}_{\mathbf{i},\mathbf{p}}(\xi ),$$

$${\mathbb{V}}_h=\text{span}\{B_{\mathbf{i},\mathbf{p}}:\mathbf{i}\in \overline{\mathcal{I}}\}.$$

$$a(u_h,v_h)=\ell (v_h),\forall v_h\in {\mathbb{V}}_h,$$

$$H^{\mathbf{s}}(\Omega ,{\mathcal{T}}_h):=\{v\in L_2(\Omega ):v{|}_{{\Omega }_i}\in H^{s_i}({\Omega }_i),i=1,\dots ,N\},$$

$$\Vert v{\Vert }_{H^{\mathbf{s}}(\Omega ,{\mathcal{T}}_h)}:={\left(\sum _{i=1}^N\Vert v{\Vert }_{H^{s_i}({\Omega }_i)}^2\right)}^{1/2}\text{and}|v{|}_{H^{\mathbf{s}}(\Omega ,{\mathcal{T}}_h)}:={\left(\sum _{i=1}^N|v{|}_{H^{s_i}({\Omega }_i)}^2\right)}^{1/2},$$

$$[[v]]:=v_i-v_j,and\{v\}:=\frac{1}{2}\left(v_i+v_j\right),F_{ij}\in {\mathcal{F}}_I,$$

$$a_h(u,v)={\ell }_h(v),\forall v\in V,$$

$$a_{h}(u,v)=\sum_{i=1}^{N}\big{(}a_{i}(u,v)+s_{i}(u,v)+p_{i}(u,v)\big{)},$$

$$a_i(u,v)$$

$$s_i(u,v)$$

$$p_i(u,v)$$

$${\ell }_h(v)$$

$$V_h:=\{v\in {\mathbb{V}}_{h,i},i=1,\dots ,N\},$$

$$a_h(u_h,v_h)=\ell (v_h),\forall v_h\in V_h.$$

$$\displaystyle\|v\|_{h}^{2}:=\sum_{i=1}^{N}\bigg{(}\alpha_{i}\|\nabla_{\Omega}v_{i}\|_{L_{2}(\Omega_{i})}^{2}+\|v\|_{L_{2}(\Omega_{i})}^{2}$$

$$\Vert \nabla v{\Vert }_{L_2({\Omega }_i)}\le C_{inv,1,u}{h}_i^{-1}\Vert v{\Vert }_{L_2({\Omega }_i)},$$

Full text

1 Introduction

In this paper, we consider the second-order elliptic boundary value problem on a open, smooth, connected and oriented two dimensional surface $\Omega\subset\mathbb{R}^{3}$ as follows: find $u:\overline{\Omega}\rightarrow\mathbb{R}$ such that

[TABLE]

where the diffusion coefficient $\alpha$ is uniformly bounded i.e. $\alpha_{min}\leq\alpha\leq\alpha_{max}$ with positive constants $\alpha_{max}$ and $\alpha_{min}$, $f$ and $g_{N}$ are given sufficiently smooth data. The physical or computational domain $\Omega\subset\mathbb{R}^{3}$ is compact, connected and positively oriented surface with boundary $\partial\Omega.$ The boundary of the computational domain consists of the Dirichlet part $\Gamma_{D}$ with positive boundary measure and a Neumann part $\Gamma_{N}$ such that $\partial\Omega:=\Gamma_{D}\bigcup\Gamma_{N}.$ The operators $\text{div}_{\Omega}$ and $\nabla_{\Omega}$ are the surface divergence and surface gradient respectively, and will be defined in Section 2.

Partial Differential Equations (PDEs) on surfaces arise in many fields of application like material science, fluid mechanics, electromagnetics, biology and image processing, see e.g.[7] for several interesting discussions on applications. For several years, numerical methods dedicated to the solutions of PDEs on manifolds including conforming and non-conforming finite element methods (FEM) have been well studied and applied to compute the solution of elliptic and parabolic evolution problems on fixed and evolving computational domains, see, e.g., [7, 5]. We note that there are however some drawbacks to the standard surface FEM. The standard surface FEM has two main sources of error: the error due to the approximation of the infinite dimensional spaces with finite dimensional spaces in the variational problem and the geometric error resulting from the approximation of the surface. These drawbacks are due to the discrete variational formulation of the PDE that is constructed on a triangulated surface which contains the finite elements space as discussed by Dzuik and Elliott in [7].

As an alternative aproach to the surface FEM, we resort to Isogeometric Analysis (IGA). The numerical scheme is based on B-splines and Non-Uniform B-splines (NURBS) and was proposed to approximate solutions of PDEs, see e.g. [10]. The method uses the same class of basis functions for representing both the geometry of the computational domain and also approximating the solution of problems modeled by PDEs. By using the exact representation of the geometry, the geometrical errors introduced by approximation of computational domains in the surface FEM are eliminated. This is especially of importance in the discretization of PDEs on surfaces. However, we note that IGA can also have geometry-related failures such as holes, singularities, etc. see e.g.[18]. Such failures or features are beyond the scope of this article. IGA uses B-splines or Non-Uniform Rational B-Splines (NURBS) basis functions which are standard in Computer Aided Design (CAD). The NURBS basis functions have several advantages making them suitable for analysis, see [10]. The mathematical analysis of the approximation properties, the stability and discretization error estimates of NURBS spaces and analysis of several refinement strategies, i.e., $h$-$p$-$k$ refinements can be found in [2]. In many practical applications, the computational domains cannot be represented by a single B-spline or NURBS domain but by several patches or sub-domains. In this sense, single patch IGA and multi-patch IGA have been addressed in [3].

Alternatively, multi-patches can also be coupled via interior penalty Galerkin methods. In our earlier articles, see e.g., [13, 12, 15], we analyzed the multi-patch discontinuous Galerkin IGA (dGIGA) for diffusion and biharmonic problems and presented several convincing numerical results for conforming domains with matching meshes. However, in this paper, we will generalize the analysis to include non-matching meshes with jumping diffusion coefficients across patch boundaries and present a priori error estimates for diffusion problems. Our analysis follows the monograph [6] and requires three main ingredients; discrete stability, consistency and boundedness of the discrete bilinear form. Then using the approximation estimates, see e.g., [2], we finally derive a priori error estimate. The linear system obtained from the discretization of the problem is solved by means of a preconditioned conjugate gradient (PCG) with a scaled Dirichlet preconditioner as presented in primal isogeometric tearing and interconnecting (dG-IETI-DP) see e.g., [9].

The rest of the paper is organized as follows; Section 2 gives a brief introduction to function spaces, weak formulation, NURBS surfaces and geometrical mappings and isogeometric analysis. We present the dGIGA scheme in Section 3. In Section 4, we present the multi-patch dGIGA and the analysis of the dGIGA scheme. The a priori error estimate is presented in Section 5. We present numerical results for an open surface and a closed surface with non-matching meshes respect to the jumping diffusion coefficient in Section 6. Finally, we conclude and give an outlook.

2 Preliminaries

In this section, we introduce briefly introduce Sobolev spaces, NURBS surfaces and isogeometric analysis method, see e.g. [1, 10] for detailed study. Firstly, let the computational domain $\Omega$ be a compact smooth and oriented surface with boundary $\partial\Omega$. We introduce the Sobolev space $H^{s}(\Omega):=\{v\in L_{2}(\Omega)\,:\,D^{\kappa}v\in L_{2}(\Omega),\,\,\text{for}\,\,0\leq|\kappa|\leq s\}$, where $L_{2}(\Omega)$ denote the space of square integrable functions and $\kappa=(\kappa_{1},\ldots,\kappa_{d})$ be a multi-index with non-negative integers $\kappa_{1},\ldots,\alpha_{d}$, and $|\kappa|=\kappa_{1}+\ldots+\kappa_{d},$ $D^{\kappa}:=\partial^{|\kappa|}/\partial x^{\kappa}.$ We associate the Sobolev space $H^{s}(\Omega)$ with the norm $\|v\|_{H^{s}(\Omega)}=\left(\sum_{0\leq|\kappa|\leq s}\|D^{\kappa}v\|^{2}_{L_{2}(\Omega)}\right)^{1/2}.$

The variational formulation of the surface diffusion problem (1.1) reads: find $u\in V_{0}$ such that

[TABLE]

where the bilinear and linear forms are given by

[TABLE]

with $V_{0}:=\{v\in H^{1}(\Omega):v=0\quad\text{on}\quad\Gamma_{D}\}.$ The existence and uniqueness of the solution of such a variational problem (2.1) follows the standard arguments of Lax-Milgram lemma if $u\in H^{2}(\Omega)$ satisfies

[TABLE]

see e.g. [17] for further details.

2.1 NURBS Geometrical Mapping and Surfaces

Let $\xi=(\xi_{1},\xi_{2})\in\mathbb{R}^{2}$ be a vector-valued independent variable in the parameter domain $\widehat{\Omega}$. By means of a smooth and invertible geometrical mapping $\mathrm{\Phi}$, the computational domain $\Omega\subset\mathbb{R}^{3}$ is defined as

[TABLE]

where the parameter domain $\widehat{\Omega}\subset\mathbb{R}^{2}$ as illustrated in Fig. 1.

We introduce briefly some important mathematical tools necessary for the analysis of surface PDEs. The following objects are obtained by means of the geometrical mapping (2.4) in the parameter domain. The Jacobian $\widehat{J},$ first fundamental form $\widehat{F}$ and the determinant $\widehat{g}$ of the geometrical mapping are respectively given by

[TABLE]

Next, we present some differential operators by using notations in the parameter domain. We consider a smooth function $\phi$ defined on the manifold $\Omega,$ by using the invertible geometrical mapping (2.4) to obtain

[TABLE]

where $\widehat{\phi}(\xi)=\phi(\mathrm{\Phi}(\xi)).$ Using the gradient operator in the parameter space $\nabla\widehat{\phi}$, the tangential gradient of the manifold is given by

[TABLE]

The divergence operator for the vector-valued function can be written as

[TABLE]

The Laplace-Beltrami operator on the manifold $\Omega$ is defined for a twice continuously differentiable function $\phi:\Omega\rightarrow\mathbb{R}$ as

[TABLE]

The unit normal vector on the manifold $\Omega\subset\mathbb{R}^{3}$ is obtained by the geometrical mapping of

[TABLE]

where $\widehat{t}_{l}(\xi):=\partial\Phi(\xi)/\partial\xi_{l}$ is the tangent vector to a curve in $\mathbb{R}^{3}$ with $l=1,2.$ The manifold $\Omega$ has a tangent plane at $\xi$ if the tangent vectors are linearly independent.

Finally, by means of the geometrical mapping (2.4), we can write the Jacobian, first fundamental form and the determinant on the computational domain $\Omega$ as follows

[TABLE]

2.2 NURBS and Isogeometric Analysis

We begin by introducing the univariate B-splines since they are usually the industry standard.Given positive integers $p$ and $n,$ we define a vector $\mathrm{\Xi}:=\left\{0=\xi_{1},\ldots,\xi_{n+p+1}=1\right\}$ with a non-decreasing sequence of real numbers in the unit interval or parameter domain $\widehat{\Omega}=[0,1]$ called a knot vector. Given $\mathrm{\Xi},p\geq 1,$ and $n$ the number of basis functions, the univariate B-spline functions are defined by the following recursion formula

[TABLE]

where a division by zero is defined to be zero. We note that a basis function of degree $p$ is $(p-m)$ times continuously differentiable across a knot value with the multiplicity $m$. If all internal knots have the multiplicity $m=1$, then B-splines of degree $p$ are globally $(p-1)-$continuously differentiable.

The bivariate B-spline basis functions are tensor products of the univariate B-spline basis functions (2.2). Let $\mathrm{\Xi}_{k}=\left\{\xi_{1,k},\ldots,\xi_{n_{k}+p_{k}+1,k}\right\}$ be the knot vectors for every direction $k=1,2$. Let $\mathbf{i}:=(i_{1},i_{2}),\mathbf{p}:=(p_{1},p_{2})$ and the set $\overline{\mathcal{I}}=\{\mathbf{i}=(i_{1},i_{2}):i_{k}=1,2,\ldots,n_{k};\;k=1,2\}$ be multi-indicies. Then the tensor product B-spline basis functions are defined by

[TABLE]

where $\xi=(\xi_{1},\xi_{2})\in\widehat{\Omega}=(0,1)^{2}.$ The univariate and bivariate B-spline basis functions are defined in the parametric domain by means of the corresponding B-spline basis functions $\{\widehat{B}_{\mathbf{i},\mathbf{p}}\}_{\mathbf{i}\in\overline{\mathcal{I}}}.$

The distinct values $\xi_{l},l=1,\ldots,n$ of the knot vectors $\mathrm{\Xi}$ provides a partition of $(0,1)^{2}$ creating a mesh $\widehat{\mathcal{K}}_{h}$ in the parameter domain where $\widehat{K}$ is a mesh element. The computational domain is described by means of a geometrical mapping $\mathbf{\Phi}$ such that $\Omega=\mathbf{\Phi}(\widehat{\Omega})$ and

[TABLE]

where $C_{\mathbf{i}}$ are the control points.

We define the basis functions in the computational domain by means of the geometrical mapping as $B_{\mathbf{i},\mathbf{p}}:=\widehat{B}_{\mathbf{i},\mathbf{p}}\circ\mathbf{\Phi}^{-1}$ and the discrete function space by

[TABLE]

Finally, the NURBS isogeometric analysis scheme reads as folows; Find $u_{h}\in\mathbb{V}_{h}$ such that

[TABLE]

with $\mathbb{V}_{h}\subset V_{0}.$

However, for many practical applications, the physical domain $\Omega$ consists of non-overlapping domains $\Omega_{i},i=1,\dots,N$ called subdomains denoted by $\mathcal{T}_{h}:=\{\Omega_{i}\}_{i=1}^{N}$ such that $\overline{\Omega}=\bigcup_{i=1}^{N}\overline{\Omega}_{i}$ and $\Omega_{i}\cap\Omega_{j}=\emptyset$ for $i\neq j.$ Each patch is the image of an associated geometrical mapping $\mathbf{\Phi}_{i}$ such that $\mathbf{\Phi}_{i}(\widehat{\Omega})=\Omega_{i},i=1,\ldots,N,$ see Fig. 2.

We denote by $F_{ij}=\partial\Omega_{i}\cap\partial\Omega_{j},i\neq j,$ the interior facets of two patches see Fig. 3. We assume the $F_{ij}\subset\partial\Omega_{i}$ for the interior facets. Let $F_{i}=\partial\Omega_{i}\cap\partial\Omega$ denote an edge of $\partial\Omega_{i}.$ $\mathcal{F}:=\mathcal{F}_{I}\cup\mathcal{F}_{D}.$

We assume that for each patch $\Omega_{i},i=1,\ldots,N,$ the underlying mesh $\mathcal{K}_{h,i}$ is quasi-uniform i.e. $h_{K}\leq h_{i}\leq C_{u}h_{K},$ for all $K\in\mathcal{K}_{h,i},\quad i=1,\ldots,N,$ where $C_{u}\geq 1$ and $h_{i}=\max\{h_{K},K\in\mathcal{K}_{h,i}\}$ is the mesh size of $\Omega_{i}$ and $h_{K}$ is the diameter of of the mesh element $K$ see e.g. [14].

3 Discontinuous Galerkin IGA Scheme formulation

We recall some function spaces required for the derivation of interior penalty Galerkin schemes. We assign to each patch $\Omega_{i}$ a a real number $s_{i}$ and collect them in the vector $\mathbf{s}=\{s_{1},\ldots,s_{N}\}.$ Let us now define the broken Sobolev space

[TABLE]

and the corresponding broken Sobolev norm and semi-norm

[TABLE]

respectively.

We denote the restrictions of the function $v$ on patches $\Omega_{i}$ and $\Omega_{j},$ by $v_{i}$ and $v_{j}$ respectively. For the interior facets $F_{ij}\subset\partial\Omega_{i},$ let $\mathbf{n}_{i}$ be the outward unit normal vector with respect to $\Omega_{i},$ which coincides with the outward unit normal $\mathbf{n}$ on $\partial\Omega,$ see Fig. 3. We define the jump and average across the interior facets $F_{ij}$ of a smooth function $v\in H^{1}(\Omega,\mathcal{T}_{h})$ by

[TABLE]

whereas the jump and average functions on the facets $F_{i}$ are given by $\llbracket{v}\rrbracket:=v_{i},$ and $\{v\}:=v_{i}.$

Now, we present the dGIGA variational scheme as follows: find $u\in V=H^{1+\mathbf{s}}(\Omega,\mathcal{T}_{h})$ with $\mathbf{s}\in(1/2,1],$ such that

[TABLE]

where the dG bilinear and linear forms considered throughout this paper are defined by the relationships where the bilinear form is given as

[TABLE]

with

[TABLE]

where $\delta$ is a non-zero positive real number. We have used a harmonic mean for the edges on the interface i.e. $h_{ij}=2h_{i}h_{j}/(h_{i}+h_{j})$ with $h_{ij}\leq 2h_{i}$ and $h_{ij}\leq 2h_{j}$ and similarly for the diffusion coefficient i.e. $\alpha_{ij}=2\alpha_{i}\alpha_{j}/(\alpha_{i}+\alpha_{j})$ with $\alpha_{ij}\leq 2\alpha_{i}$ and $\alpha_{ij}\leq 2\alpha_{j}.$ The linear form is given by

[TABLE]

where $\mathcal{F}_{N}$ is the collection of all edges on the Neumann Boundary parts.

Remark 1

The choice of the penalty parameter $\delta$ depends on B-spline or NURBS degree $p$ and the dimension of the computational domain $\Omega\subset\mathbb{R}^{d},$ for example in FEM $\delta=2(p+1)(p+d)/d,$ see e.g. [16].

4 Analysis of the dGIGA Scheme

For each subdomain $\Omega_{i},i=1,\ldots,N,$ we will consider the discrete space $\mathbb{V}_{h,i},i=1,\ldots,N$ where $\mathbb{V}_{h}$ is given by (2.16). We define the discrete space corresponding to the domain $\Omega$ as

[TABLE]

which allows discontinuities across the patch interface. The discrete dGIGA scheme then reads as: find $u_{h}\in V_{h}$ such that

[TABLE]

The existence and uniqueness of the bilinear form $a_{h}(\cdot,\cdot)$ follow the popular Lax-Milgram theorem by showing the coercivity and boundedness. Next, we show that the bilinear form $a_{h}(\cdot,\cdot)$ is $V_{h}-$coercive with respect to the dG-norm

[TABLE]

Remark 2

The discrete norm $\|\cdot\|_{h}$ from (4.2) is a norm on $V_{h}.$ Indeed, if $\|v\|_{h}=0$ for some function $v\in V_{h}$, then $\nabla_{\Omega}v=0$ in each subdomain $\Omega_{i}$. This means that the function $v$ is a constant on each patch $\Omega_{i}$, $i=1,\ldots,N$. Furthermore, $\|v\|_{h}=0$ yields that $\llbracket{v}\rrbracket=0$ across the internal facets $F_{ij}\subset\partial\Omega_{i}$ are zero, i.e., $v$ is constant in $\overline{\Omega}$. Finally, $v_{i}=0$ on the boundary $\partial\Omega$ implies that this constant must be zero. Thus, $v=0$ in $\overline{\Omega}$. The other norm axioms are obviously fulfilled.

To analyze the multi-patch interior penalty Galerkin scheme, the following discrete inverse and trace inequalities are required.

Lemma 1

Let $v\in V_{h},$ then the following inverse inequalities hold;

[TABLE]

and

[TABLE]

where $C_{inv,1,u}$ and $C_{inv,0,u}$ are positive constants, which are independent of $h_{i}$ and $\Omega_{i}$.

We conclude with the continuous trace inequality,

Lemma 2

Let $\Omega_{i}=\Phi_{i}(\widehat{\Omega})$ for $i=1,\dots,N$. Then the patch-wise scaled trace inequality

[TABLE]

holds for all $v\in H^{1/2+\epsilon}(\Omega_{i}),\epsilon\in(0,1/2],$ where $h_{i}$ denotes the maximum mesh size in the physical domain, and $C_{t,u}$ is a positive constant that only depends on the shape regularity of the mapping $\Phi_{i}.$

The proofs of Lemma 1 and Lemma 2 follows the standard procedure see e.g., [8] and [14]. from the Finite Element.

Lemma 3

For an arbitrary positive $\varepsilon$ and for $F_{ij}\subset\partial\Omega_{i}$ the estimates

[TABLE]

holds for all $v_{h,i},v_{h,j}\in V_{h}$, a positive constant $C_{inv,0,u}$ and $\alpha_{i}>0.$

Proof

Following the Cauchy Schwarz inequality, since $F_{ij}\subset\partial\Omega_{i},$ by using the trace inequality (4.4), we have

[TABLE]

where we use $h_{ij}\leq 2h_{i}$ and $\alpha_{ij}\leq 2\alpha_{i},$ together with the inequality $ab\leq\varepsilon a^{2}/2+b^{2}/(2\varepsilon),\forall a,b\in\mathbb{R}$ with $\varepsilon>0$ to complete the proof. ∎

Using the above result, we proceed to show the coercivity of the bilinear form $a_{h}(\cdot,\cdot).$

Lemma 4 (Coercivity)

Let $a_{h}(\cdot,\cdot):V_{h}\times V_{h}\rightarrow\mathbb{R}$ be the bilinear form (3.5). There exists $\delta_{0}>0$ and $\mu_{c}>0$ such that $\delta>\delta_{0}$ and the discrete bilinear form $a_{h}(\cdot,\cdot)$ is $V_{h}-$coercive with respect to the norm $\|\cdot\|_{h}$, i.e.

[TABLE]

where $\mu_{c}$ is independent of $\alpha_{i},h_{i}$ and $N.$

Proof

By using Cauchy-Schwarz’s inequality, we proceed as follows

[TABLE]

Using Cauchy Schwarz’s inequality and Lemma 3, we have

[TABLE]

For example, for $\mu_{c}=1/2,$ we choose $\varepsilon=1$ and $\delta\geq C^{2}_{inv,0,u}.$ ∎

Next, we prove the uniform boundedness for the bilinear form $a_{h}(\cdot,\cdot)$ on $V_{h,*}\times V_{h},$ where $V_{h,*}=V_{0}\cap H^{1+\mathbf{s}}(\Omega,\mathcal{T}_{h})+V_{h}$ with $\mathbf{s}>1/2$ is equipped with the norm

[TABLE]

To prove a priori estimates, we need to show the uniform boundedness of the bilinear form. We need the following two auxiliary lemmata to proof for the boundedness of the bilinear form.

Lemma 5

For a positive parameter $\delta$ and for $F_{ij}\subset\partial\Omega_{i},i=1,\ldots,N$ and diffusion coefficients $\alpha_{i}$ and $\alpha_{ij},$ the estimates

[TABLE]

hold for all $u\in V_{h,*}$ and for all $v_{h}\in V_{h}.$

Proof

For $F_{ij}\subset\partial\Omega_{i},$ using Cauchy-Schwarz’s inequality, we obtain

[TABLE]

We conclude the proof since $h_{ij}\leq 2h_{i}$ and $\alpha_{ij}\leq 2\alpha_{i}$ For the second inequality, we apply the Cauchy Schwarz inequality to obtain

[TABLE]

Since $h_{ij}\leq 2h_{i}$ and $\alpha_{ij}\leq 2\alpha_{i},$ by applying the inequality (4.4) for $v_{h,i}\in V_{h},$ we complete the proof. ∎

Next, we proceed with the boundedness of the bilinear form $a_{h}(\cdot,\cdot)$ as follows ;

Lemma 6 (Boundedness)

The discrete bilinear form $a_{h}(\cdot,\cdot)$ is uniformly bounded on $V_{h,*}\times V_{h},$ i.e. there exists a mesh-independent positive constant $\mu_{b}$ such that

[TABLE]

Proof

The first term of the bilinear form (3.5) is estimated using Cauchy-Schwarz’s inequality as follows

[TABLE]

Using Cauchy Schwarz’s inequality together with Lemma 5, the second term yields

[TABLE]

Also, the last term of the bilinear form is estimated by applying Cauchy-Schwarz’s inequality to obtain

[TABLE]

Combining all the terms (4.12) – (4.14), we conclude the proof with the positive constant $\mu_{b}=2\sqrt{\max\{1,(1+C_{inv,0,1}^{2}/\delta)\}}.$ ∎

We note that the discrete norms $\|\cdot\|_{h}$ and $\|\cdot\|_{h,*}$ are uniformly equivalent on the dicrete space $V_{h}.$ In the next lemma, we present this equivalence of the discrete norms since the convergence analysis is considered in the discrete norm $\|\cdot\|_{h}.$

Lemma 7

The norms $\|\cdot\|_{h}$ and $\|\cdot\|_{h,*}$ are uniformly equivalent on the discrete space $V_{h}$ such that

[TABLE]

where $C_{e}$ is mesh independent.

Proof

The proof of the upper bound follows immediately. However, the proof of the lower bound follows by using the definition of the norm (4.8) together with the trace inequality (4.4) with $C_{e}=\left(1+C_{inv,0,1}^{2}/\delta\right)^{-1}.$ ∎

A consequence of Lemma 7 yields the boundedness of the bilinear form $a_{h}(\cdot,\cdot)$ with a positive constant $\tilde{\mu}_{b}$ as follows

[TABLE]

where $\tilde{\mu}_{b}=\mu_{b}\left(1+C_{inv,0,1}^{2}/\delta\right)^{-1/2}.$

5 Error Analysis of dGIGA Discretization

Finally, we present the approximation estimates required to obtain a priori error estimates. For patch $\Omega_{i},i=1,\ldots,N,$ let $\Pi_{h,i}:L_{2}(\Omega_{i})\rightarrow\mathbb{V}_{h,i}$ denote a quasi-interpolant that yields optimal approximation results. Of course, such an interpolant is known to exist and has been well studied and presented in [2, 3] as follows

Lemma 8

Let $l$ and $s$ be integers with $0\leq l\leq s\leq p+1$ and $K\in\mathcal{K}_{h,i}$. Then there exist an interpolant $\Pi_{h,i}v\in\mathbb{V}_{h,i}$ for all $v\in H^{s}(\Omega_{i})$ and a constant $C_{s}>0$ such that the following inequality holds

[TABLE]

where $h_{i}$ is the mesh size in the physical domain, and $p$ denotes the underlying polynomial degree of the B-spline or NURBS.

For patch $\Omega_{i},i=1,\ldots,N,$ the local estimate (5.1) yields a global estimate if the multiplicity of the inner knots is not larger than $p+1-l$ and $\Pi_{h,i}v\in\mathbb{V}_{h,i}\cap H^{l}(\Omega_{i}).$

Proposition 1

Let $v\in H^{s}(\Omega_{i})$ be a function defined in the physical domain $\Omega_{i}.$ Given an integer $l$ such that $0\leq l\leq p+1,l\leq s,$ and $p\leq s+1,$ where $s$ is the smoothness of the considered B-Spline basis. Then there exists a projection operator $\Pi_{h,i}:L_{2}(\Omega_{i})\rightarrow\mathbb{V}_{h,i}$ such that the approximation error estimates

[TABLE]

where $\beta=\min\{p+1,l\},$ $h_{i}$ denotes the maximum mesh-size parameter in the physical domain and the generic constant $C_{s}$ only depends on $l,s$ and $p$, the shape regularity of the physical domain $\Omega_{i}$ described by the mapping $\Phi$ and, in particular, $\nabla_{\Omega}\mathrm{\Phi}.$

Proof

See [4, Proposition 3.2].

For the error analysis, we assume that the patches have the same regularity such that $\mathbf{s}=\{s_{1},s_{2},\ldots,s_{N}\}=s$ and $H^{1+s}(\Omega_{i}),i=1,\ldots,N.$

Lemma 9

Let $v\in V_{0}\cap H^{1+s}(\Omega_{i})$ with $s>1/2$ and $p\geq 1.$ By assuming quasi-uniform meshes, then there exists a projection $\Pi_{h}v\in V_{h}$ and generic positive constants $C_{0}$ and $C_{1}$ such that the following error estimates hold

[TABLE]

where $r=\min\{s,p\},$ the constant $C_{0}$ and $C_{1}$ $\Omega_{i},$ are independent of $h_{i}$ and $h_{j}$ and $F_{ij}\subset\partial\Omega_{i}$ are the interior facets.

Proof

By using (3.3) with Lemma 2 and Proposition 1, we estimate the first term as follows

[TABLE]

Similarly, the second term

[TABLE]

Now, we complete the proof by summing with respect to the interior facets $F_{ij}\subset\partial\Omega_{i}$ and $i=1,2,\ldots,N,$ to obtain

[TABLE]

where $C_{0}=\delta C_{t,u}^{2}C_{s}.$ The proof of (5.4) follows by using Lemma 2 and the approximation estimate of Proposition 1 as follows

[TABLE]

where $C_{1}=2C_{s}C_{t,u}^{2}.$ ∎

To derive the a priori error estimate, we show that the interpolant yields the optimal approximation estimate in the discrete norms.

Lemma 10

Let $v\in V_{0}\cap H^{1+s}(\Omega_{i})$ with $s>1/2$ and $p\geq 1.$ Then there exists a projection $\Pi_{h}v\in V_{h}$ and generic positive constants $C_{2}$ and $C_{3}$ such that

[TABLE]

where $F_{ij}\subset\partial\Omega_{i}$ are the interior facets, of $\Omega_{i},$ $r=\min\{s,p\}$ and $C_{2}$ and $C_{3}$ only depend on $s$ and $p.$

Proof

Following from the definition of the discrete norms (4.2) and (4.8) together with Lemma 9, we complete the proof. ∎

Finally, we prove the main result in this section, namely a priori error estimate for surfaces. We will present the results for the discrete norm $\|\cdot\|_{h}$ and the $\|\cdot\|_{L_{2}(\Omega)}-$norm.

Theorem 5.1

Let $u\in V_{0}\cap H^{1+s}(\Omega_{i})$ with $s>1/2$ be the exact solution of the model (2.1) and $u_{h}\in V_{h}$ with $p\geq 1$ be the discrete solution of the dGIGA scheme (4.1). For the penalty parameter $\delta$ chosen as in Lemma 4, then the discretization error estimate

[TABLE]

holds true, where $r=\min\{s,p\}$ and $p$ denotes the underlying NURBS degree of the patch $\Omega_{i},$ and $C$ is a positive constant independent of the $h_{i}$ and $h_{j}.$

Proof

By using the coercivity result Lemma 4, Galerkin orthogonality (5) and the boundedness of the discrete bilinear form, Lemma 6, we obtain

[TABLE]

Thus, we have

[TABLE]

Using Lemma 9, we get

[TABLE]

where $C=(C_{3}+\left(\mu_{b}/\mu_{c}\right)^{2}C_{4}).$ Lemma 6. ∎

Remark 3

If, we assume matching meshes i.e. $h_{i}=h_{j},$ then the a priori error estimate (5.11) yields

[TABLE]

which has been studied and presented in [13].

6 Numerical Results

In this section, we present numerical results for the dGIGA scheme and a priori error estimate of Theorem 5.1. All the numerical experiments have been performed in G**+**Smo see [11]. We solve the linear system arising from the dGIGA formulation by means a preconditioned conjugate gradient (PCG) algorithm with a scaled Dirichlet preconditioner where we choose vertex evaluation and edge averages as primal variables in the so-called dual-primal isogeometric tearing and interconnecting (dG-IETI-DP) solver. The solver is known to be robust with respect to diffusion coefficient see e.g. [9]. A reduction of the initial residual factor of $10^{-6}$ is used as a stopping criterion together with a zero initial guess. In the examples, we present non-matching grid of ratio $h_{i}/h_{j}=2^{q},$ where $q$ is the mesh refinement. The ratio $h_{i}/h_{j}$ denotes the relative number of refinement on the neighboring patches and $h_{i},h_{j}$ are the maximum mesh sizes of patches $\Omega_{i}$ and $\Omega_{j}.$ The penalty parameter is chosen to be $\delta=2(p+2)(p+1),$ where $p$ is the NURBS degree. The convergence rate is computed using the formula $rate=\log_{2}\left(e_{i+1}/e_{i}\right),$ where $e_{i+1}=\|u-u_{h,i+1}\|_{h}$ and $e_{i}=\|u-u_{h,i}\|_{h}$ to study the discrete solution of the model problem. We consider as computational domains a quarter cylinder and a torus for the open and closed surfaces respectively, see Fig. 4.

6.1 Open Surface

We consider a diffusion problem with homogeneous Dirichlet boundary condition on an open surface $\Omega$ that is given by a quarter cylinder in the first quadrant i.e. $x\geq 0$ and $y\geq 0$ with unitary radius and height $L=4.$ The computational domain $\Omega$ is decomposed into 4 patches, with each of the patches having height of $L=1$ and depicted by different color as seen on the left-hand side of Fig. 4 (left). The knot vectors representing the geometry of each patch are given by $\Xi_{1}=\{0,0,0,1,1,1\}$ and $\Xi_{2}=\{0,0,1,1\}$ in the $\xi_{1}-$direction and $\xi_{2}-$direction respectively. Let $f(\phi,z)=\varrho\left(\frac{\sigma^{2}\pi^{2}}{L^{2}}g_{\phi,1}(\phi)-g_{\phi,2}(\phi)\right)g_{z}(z),$ where $\phi:=\arctan\left(\frac{x}{y}\right),$ $g_{\phi,1}(\phi):=(1-\cos(\phi))(1-\sin(\phi)),$ $g_{\phi,2}(\phi):=(\cos(\phi)+\sin(\phi)-4\sin(\phi)\cos(\phi)),$ and $g_{z}(z):=\sin\left(\sigma\pi\frac{z}{L}\right)$ for $\sigma\in\mathbb{N}_{0}$ and $\varrho>0.$ The exact solution of the problem is $u(\phi,z)=\varrho g_{\phi,1}(\phi)g_{z}(z).$ In our numerical experiments, we set $\sigma=3,\varrho=1/\left(3/2-\sqrt{2}\right).$ Fig. 5. We present the convergence behavior of the dGIGA scheme with respect to the discrete norm $\|\cdot\|_{h}$ in Fig.. 5 by successive mesh refinement of ratio $h_{i}/h_{j}=2^{q},$ where $q=1,2,3$ and $q=4$ are the refinement level using NURBS of degree $p=2$ and $p=4.$ We observe the optimal convergence rate as theoretically predicted in Theorem 5.1 for smooth functions.

6.2 Closed Surface

We consider the closed surface

[TABLE]

that is nothing but a torus decomposed into 4 patches, see Fig. 4 (right). The knot vectors describing the NURBS used for the geometrical representation of the patches $\Xi_{1}=\{0,0,0,0.25,0.25,0.50,0.50,0.75,0.75,1,1,1\}$ and $\Xi_{2}=\{0,0,0,1,1,1\}$ . Let us consider the surface Poisson equation with the right-hand side

[TABLE]

where $\phi=\arctan(y/x)$, $\theta=\arctan(z/(\sqrt{x^{2}+y^{2}}-R))$, $R=2$ and $r=1$. The exact solution is given by $u=\sin(3\phi)\cos(3\theta+\phi)$. The functions $u$ and $f$ are chosen such that the zero mean compatibility condition holds. We present the convergence behavior of the dGIGA scheme with respect to the discrete norm $\|\cdot\|_{h}$ by successive mesh refinement of ratio $h_{i}/h_{j}=2^{q},$ where $q=1,2,3$ and $q=4$ are the number of mesh refinements using NURBS degrees $p=2$ and $p=4$ see Fig. 6. We observe the optimal convergence rate as theoretically predicted in Theorem 5.1 for smooth functions.

Conclusion

In this article, we considered the discontinuous Galerkin isogeometric analysis (dGIGA) for the surface diffusion problem with jumping coefficient and geometrically non-matching meshes. We analyzed the well-posedness and presented a priori error estimates. Finally, we presented numerical results confirming the theory presented. In solving the linear system arising from the dGIGA scheme, we applied the dual-primal discontinuous Galerkin isogeometric tearing and interconnecting method (dG-IETI-DP). This involved a Preconditioned Conjugate Gradient (PCG) algorithm with the scaled Dirichlet preconditioner which is known to be robust with respect to jumping diffusion coefficient. An extension of the results to non-orientable surfaces as well evolving surfaces will be considered in our next article.

Acknowledgement

The author acknowledges the Horizon 2020 Programme (2014-2020) under grant agreement number 678727.