Discontinuous Galerkin Methods for Isogeometric Analysis for Elliptic Equations on Surfaces

Abstract

We propose a method that combines isogeometric analysis (IGA) with the discontinuous Galerkin (DG) method for solving elliptic equations on 3-dimensional (3D) surfaces consisting of multiple patches. DG ideology is adopted across the patch interfaces to glue the multiple patches, while the traditional IGA, which is very suitable for solving partial differential equations (PDEs) on (3D) surfaces, is employed within each patch. Our method takes advantage of both IGA and the DG method. Firstly, the time-consuming steps in mesh generation process in traditional finite element analysis (FEA) are no longer necessary and refinements, including $h$-refinement and $p$-refinement which both maintain the original geometry, can be easily performed by knot insertion and order-elevation (Farin, in Curves and surfaces for CAGD, 2002). Secondly, our method can easily handle the cases with non-conforming patches and different degrees across the patches. Moreover, due to the geometric flexibility of IGA basis functions, especially the use of multiple patches, we can get more accurate modeling of more complex surfaces. Thus, the geometrical error is significantly reduced and it is, in particular, eliminated for all conic sections. Finally, this method can be easily formulated and implemented. We generally describe the problem and then present our primal formulation. A new ideology, which directly makes use of the approximation property of the NURBS basis functions on the parametric domain rather than that of the IGA functions on the physical domain (the former is easier to get), is adopted when we perform the theoretical analysis including the boundedness and stability of the primal form, and the error analysis under both the DG norm and the ${\mathcal {L}}^{2}$ norm. The result of the error analysis shows that our scheme achieves the optimal convergence rate with respect to both the DG norm and the ${\mathcal {L}}^{2}$ norm. Numerical examples are presented to verify the theoretical result and gauge the good performance of our method.

1 Introduction

In recent years, the attempt to bridge the gap between computer aided design (CAD) and computer aided engineering (CAE) has made important progress through the introduction of Isogeometric Analysis proposed by Hughes et al [2]. What makes the gap noticeable is the time-consuming steps in mesh generation process and difficulties encountered during refinement. For the former, finite element meshes are to be generated for analysis from an analysis suitable geometry (ASG) which is obtained from the CAD model including ambiguities and levels of detail. On the other hand, a series of refinements are often performed to get an accurate solution in the traditional FEA. In the meantime, to make the error decrease as one wish, the geometric approximation to the exact geometry should be better and better along with the refinements, which obviously needs a link between the ASG and the refinement routine. Unfortunately, such a link usually does not exist in practice.

The gap between design and analysis can be closed by IGA in which the numerical solution and the geometry share the same basis functions $-$ usually NURBS (Non-Uniform Rational B-Splines) [1]. The time-consuming steps in mesh generation process in FEA are no longer necessary. Indeed, the mesh itself for describing the geometry is suitable for analysis as well. Additionally, $h$-refinement and $p$-refinement, which both maintain the original geometry, can be easily carried out for the mesh and the corresponding basis functions at any level by knot insertion and order-elevation [1], respectively. Thus, the link between the ASG and the refinement routine is no longer necessary. What is more, with the flexibility of IGA basis functions and the use of multi-patches, we can get more accurate modeling of complex geometries and exactly represent common engineering shapes such as circles, cylinders, spheres, ellipsoids, etc. As a result, the geometrical error is reduced or eliminated at the beginning.

IGA does a perfect job when dealing with problem whose geometry can be described using a single, watertight patch. Both theoretical and practical results have shown that IGA methods are advantageous compared to traditional FEA with respect to convergence and accuracy of the results, see [2–8]. However, it seems difficult or impossible to create a globally smooth parametrization for a complex geometry. Generally an assemblage of multiple patches joining with at least ${\mathrm {C}}^{0}$-continuity is required to represent such complicated geometry. In this case, since each patch has its own NURBS representation, traditional IGA is not feasible, whereas Discontinuous Galerkin-based IGA should play an important role.

Discontinuous Galerkin (DG) methods are a class of finite element methods using completely discontinuous basis functions, which are usually chosen as piecewise polynomials. The first DG method was introduced in 1973 by Reed and Hill [9], in the framework for a time-independent linear hyperbolic equation. Since then, Cockburn, Shu, and their collaborators made a major development of the DG method in a series of papers [10–14] for nonlinear hyperbolic conservation laws. Subsequently, these methods were also applied to purely elliptic problems, see [15–19]. On the other hand, Galerkin methods for elliptic and parabolic equations using discontinuous finite elements were independently proposed in 1970s, and a number of variants were introduced and studied, see, for example, [20–24]. These DG methods were usually called interior penalty (IP) methods, and their development remained independent of the development of the DG methods for hyperbolic equations. In 2002, Arnold et al [25] presented a detailed study of a class of DG methods for second-order elliptic problems which includes all the above-mentioned DG methods and performed a unified error analysis of these methods.

Methods for gluing multiple patches for two-dimensional (2D) plane case and (3D) volume case in IGA were mentioned or studied in many literatures, see [2, 7, 26–32]. However, many applications in practice always need to solve PDEs on complicated (3D) surfaces. Examples of applications are, amongst others, the motion of grain boundaries in alloys, phase transitions, and pattern formation on the surfaces of growing organisms. Specifically, the Laplace–Beltrami operator frequently appears in scientific computations and has a wide range of applications, for example, in shape analysis, segmentation, database retrieval, and so on. Readers can refer to the survey paper [33] for details about different finite element methods for solving PDEs on surfaces. However, the traditional FEA surfers a lot from the time-consuming steps in mesh generation process and the difficulties encountered during refinements. Thus, it is indispensable to perform an efficient and accurate method to solve PDEs on (3D) surfaces.

In this paper, we deal with second-order elliptic problems on (3D) surfaces by combining IGA method with the interior penalty DG (IPDG) method. Since DG method use completely discontinuous basis functions, it can be exactly used to deal with the multi-patch case in IGA, where the discontinuities of the basis functions only appear across the patch interfaces. We remark that the DG ideology is adopted at patch level, i.e., we employ the traditional IGA within each patch, and employ the IPDG method across the patch interfaces to glue the multiple patches. The advantageous features of both IGA and the IPDG method enable us to design a superior finite element formulation. First and most importantly, the use of IGA makes the time-consuming steps in mesh generation process in traditional FEA no longer necessary and refinements, including $h$-refinement and $p$-refinement, can be easily performed by knot insertion and the order-elevation. Secondly, this method can easily and nicely handle cases with non-conforming patches and different degrees across the patch interfaces due to the flexibility of DG method. Thirdly, the geometrical error is eliminated or significantly reduced at the beginning due to the geometric flexibility of IGA basis functions, especially the use of multiple patches. Lastly, this method can be easily formulated and implemented. Recently we are aware of a similar work by Langer et al [34]. Although the idea in [34] is similar to ours, both proofs are different. Furthermore, in our numerical simulations, we can really handle the non-matching mesh.

The remainder of this paper is organized as follows. A brief review of differential operators on surfaces is provided in Sect. 2.1. We generally describe the physical domain and geometrical mappings in Sect. 2.2. Notations for the partitions of the domain and corresponding function spaces are introduced in Sect. 2.3. Additionally, we introduce some notations for the IPDG method in Sect. 2.4. The primal formulation is presented in Sect. 2.5. Theoretical analysis is carried out in Sect. 3. We begin by proving some basic properties of the physical domain and geometrical mappings, and then study the classical properties of boundedness, and stability typically used in finite element analysis after we define appropriate semi-norms and norms for function spaces needed. We obtain our error estimates in Sect. 3.6. Numerical examples are presented in Sect. 4 to verify our theoretical result. We conclude the paper in Sect. 5.

2 DG Method for IGA on Surfaces

In this section, we present our DG formulation after a general description of the physical domain, geometrical mappings, and an introduction of notations for the DG method. A pre-knowledge about differential operators on surfaces is also included here.

We remark that we distinguish the notation of a function on the physical surface from that of a function on the parametric region by putting a tilde over a symbol in Sects. 2 and 3. For example, $f$ is used to denote a function on the parametric region, while $\widetilde{f}$ is used to denote a function on the physical surface.

2.1 Differential Operators on Surfaces

We begin by providing a brief review of the differential operators on surfaces and their computational formulas used later. These can be found in any books about differential geometry, we refer to [35] for more details.

Assume $\mathbf{F }(\xi , \zeta )$ is a parametric representation of a surface ${\mathcal {S}}$, for a function $\Phi (\mathbf{F }) = \Phi (\mathbf{F }(\xi , \zeta ))$ on ${\mathcal {S}}$, let

$$\begin{aligned} \Phi _{\xi } = \frac{\partial \Phi (\mathbf{F }(\xi , \zeta ))}{\partial \xi },\ \ \ \Phi _{\zeta } = \frac{\partial \Phi (\mathbf{F }(\xi , \zeta ))}{\partial \zeta }. \end{aligned}$$

$\Phi _{\xi \xi }, \Phi _{\xi \zeta }$ and $\Phi _{\zeta \zeta }$ can be defined analogously. $\Phi$ is called ${\mathrm {C}}^{1}$-continuous on ${\mathcal {S}}$ if $\Phi _{\xi }, \Phi _{\zeta }$ exist and are both continuous with respect to $\xi , \zeta$. We denote by ${\mathrm {C}}^{1}({\mathcal {S}})$ all the ${\mathrm {C}}^{1}$-continuous functions on ${\mathcal {S}}$. Analogously, we can define ${\mathrm {C}}^{1}({\mathcal {S}}), {\mathrm {C}}^{2}({\mathcal {S}})$ and so on. Let $G$ be the matrix defined as

$$\begin{aligned} G := \left( \begin{array}{c@{\quad }c} g_{11} &{} g_{12} \\ g_{21} &{} g_{22} \end{array} \right) , \end{aligned}$$

(2.1)

where

$$\begin{aligned} \begin{array}{ccc} g_{11} = (\Phi _{\xi }, \Phi _{\xi }), &{} g_{12} = g_{21}= (\Phi _{\xi }, \Phi _{\zeta }), &{} g_{22} = (\Phi _{\zeta }, \Phi _{\zeta }),\\ \end{array} \end{aligned}$$

(2.2)

with $(\cdot , \cdot )$ the vector inner product. Additionally, we let $g := {\mathrm {det}}(G),\ \ [g^{\alpha \beta }] := G^{-1}$.

Tangent gradient operator Assume $\Phi \in {\mathrm {C}}^{1}({\mathcal {S}})$, then the tangent gradient operator $\nabla _{s}$ on $\Phi$ is defined as

$$\begin{aligned} \begin{array}{lll} \displaystyle \nabla _{s}\Phi &{} = &{} \displaystyle \left[ \mathbf{F }_{\xi }, \mathbf{F }_{\zeta }\right] \left[ g^{\alpha \beta }\right] \left[ \Phi _{\xi },\Phi _{\zeta }\right] ^{T} \in {\mathbb {R}}^{3}\\ &{} = &{} \displaystyle g_{\xi }^{\nabla } \Phi _{\xi } + g_{\zeta }^{\nabla } \Phi _{\zeta }, \end{array} \end{aligned}$$

(2.3)

where $[\cdot ]^{T}$ denotes the transposition of a matrix or a vector, and

$$\begin{aligned} g_{\xi }^{\nabla } = \frac{1}{g}\left( g_{22}\mathbf{F }_{\xi } - g_{12}\mathbf{F }_{\zeta }\right) , \ \ \ g_{\zeta }^{\nabla } = \frac{1}{g}\left( g_{11}\mathbf{F }_{\zeta } - g_{12}\mathbf{F }_{\xi }\right) . \end{aligned}$$

From (2.3), we can easily get

$$\begin{aligned} \begin{array}{lll} \displaystyle \nabla _{s}\Psi \cdot \nabla _{s}\Phi &{} = &{} \displaystyle \left[ \Psi _{\xi },\Psi _{\zeta }\right] \left[ g^{\alpha \beta }\right] \left[ \Phi _{\xi }, \Phi _{\zeta }\right] ^{T} \\ &{} = &{} \displaystyle \frac{1}{g}\Big ((g_{22}\Psi _{\xi } - g_{12}\Psi _{\zeta })\Phi _{\xi } + (g_{11}\Psi _{\zeta } - g_{12}\Psi _{\xi })\Phi _{\zeta }\Big ) \ \ \Phi ,\Psi \in {\mathrm {C}}^{1}({\mathcal {S}}). \end{array} \end{aligned}$$

(2.4)

For vector-valued function ${\varvec{\Phi }} = [\Phi _{1},\cdots , \Phi _{m} ]^{T}$, we define that

$$\begin{aligned} \nabla _{s}{\varvec{\Phi }} = \left[ \nabla _{s}\Phi _{1},\cdots ,\nabla _{s}\Phi _{m} \right] \in {\mathbb {R}}^{3\times m} \ \ \ {\varvec{\Phi }} \in \left[ {\mathrm {C}}^{1}({\mathcal {S}})\right] ^{m}. \end{aligned}$$

(2.5)

Tangent divergence operator Assume ${\varvec{v}}$ is a ${\mathrm {C}}^{1}$ vector field on ${\mathcal {S}}$, then the tangent divergence operator ${\mathrm {div}}_{s}$ on ${\varvec{v}}$ is defined as

$$\begin{aligned} {\mathrm {div}}_{s}({\varvec{v}}) = \frac{1}{\sqrt{g}}\bigg [\frac{\partial }{\partial \xi }, \frac{\partial }{\partial \zeta }\bigg ] \bigg [\sqrt{g}[g^{\alpha \beta }][\mathbf{F }_{\xi }, \mathbf{F }_{\zeta }]^{T}{\varvec{v}}\bigg ]. \end{aligned}$$

Laplace–Beltrami (L–B) operator Let $\Phi \in {\mathrm {C}}^{2}({\mathcal {S}})$, then $\nabla _{s}\Phi$ is a smooth vector field on ${\mathcal {S}}$. The L-B operator $\Delta _{s}$ on $\Phi$ is defined as

$$\begin{aligned} \Delta _{s}\Phi = {\mathrm {div}}_{s}(\nabla _{s}\Phi ). \end{aligned}$$

From the definition of $\nabla _{s}$ and ${\mathrm {div}}_{s}$, one can easily derive that

$$\begin{aligned} \Delta _{s}\Phi = \frac{1}{\sqrt{g}}\bigg [\frac{\partial }{\partial \xi }, \frac{\partial }{\partial \zeta }\bigg ] \bigg [\sqrt{g}[g^{\alpha \beta }][\Phi _{\xi }, \Phi _{\zeta }]^{T}\bigg ]. \end{aligned}$$

Lemma 2.1

(Green’s second identity for L-B operator [35]) Assume ${\mathcal {S}}$ is an orientable surface, and $\Omega$ a subsection of ${\mathcal {S}}$ with piecewise smooth boundary $\partial \Omega$. Denote by ${\varvec{n}}$ the outward unit normal for $\Omega$ along $\partial \Omega$, then for all ${\mathrm {C}}^{1}$ vector field ${\varvec{v}}$ on ${\mathcal {S}}$ and $\Phi \in {\mathrm {C}}^{1}({\mathcal {S}})$, we have

$$\begin{aligned} \int _{\Omega }\Big ( {\varvec{v}}\cdot \nabla _{s}\Phi + \Phi \ {\mathrm {div}}_{s}({\varvec{v}})\Big ) \, {\mathrm {d}}A = \int _{\partial \Omega }\Phi \ {\varvec{v}}\cdot {\varvec{n}} \, {\mathrm {d}} \lambda . \end{aligned}$$

2.2 Physical Domain and Geometrical Mappings

The physical domain ${\mathfrak {D}}$, composed of ${\mathfrak {M}}$ non-overlapping NURBS patches $\{{\mathfrak {D}}_{k}\}_{k = 1}^{{\mathfrak {M}}}$, is at least first-order geometrically continuous (${\mathrm {G}}^{1}$-continuous) [1] across the patch interfaces, i.e., the outward normals of different patches along patch interfaces are parallel.

Fig. 1

The physical domain, composed of several colored NURBS patches $\{\mathfrak {D}_{k}\}_ {k=1}^{\mathfrak {M}}$, is $\mathrm {G}^{1}$- continuous across the patch interfaces. $\varvec{n}_{i}$ and $\varvec{n}_{j}$ denote the unit outward normals of adjacent patches along the common interface, $\varvec{n}_{i} = - \varvec{n}_{j}$

For each $k \in \{1, \ldots , {\mathfrak {M}}\}$, the geometrical mapping $\mathbf{F }_{k}: D_{k} \mapsto {\mathfrak {D}}_{k}$ takes the following form:

$$\begin{aligned} \mathbf{F }_{k}(\xi , \zeta ) = \sum _{i=1}^{m_{k}}\sum _{j=1}^{n_{k}} \mathbf{P}_{i,j}^{k} R_{i,j}^{p_{k},q_{k}}(\xi , \zeta ) \quad (\xi , \zeta ) \in D_{k}, \end{aligned}$$

(2.6)

where $D_{k} = [a_{k}, b_{k}]^2$, referred to as a parametric patch, is an square which is always taken as $[0,1]^{2}$ in literatures on computer aided geometric design (CAGD). $\{\mathbf{P}_{i,j}^{k} \in {\mathbb {R}}^{3}\}_{i=1,j=1}^{m_{k},n_{k}}$ are called control points, and $\{R_{i,j}^{p_{k},q_{k}}(\xi , \zeta )\}_{i=1,j=1}^{m_{k},n_{k}}$ represent bivariate NURBS basis functions of bi-degree $(p_{k}, q_{k})$. $m_{k}, n_{k}$ denote the numbers of bivariate NURBS basis functions on $D_{k}$ in $\xi$-direction and $\zeta$-direction, respectively. We assume each $\mathbf{F }_{k}$ is smooth enough with a smooth inverse. Figure 1 gives an illustration of the physical domain.

Without loss of generality, we make a fundamental assumption that

$$\begin{aligned} \frac{\partial {\varvec{F}}_{k}}{\partial \xi } \times \frac{\partial {\varvec{F}}_{k}}{\partial \zeta } \ne {\varvec{0}} \quad \forall (\xi , \zeta )\in D_{k},\ k =1, \ldots , {\mathfrak {M}}, \end{aligned}$$

(2.7)

that is, for each $k\in \{1,\ldots ,{\mathfrak {M}} \}$ and $\forall (\xi , \zeta )\in D_{k}$, vector $\frac{\partial {\varvec{F}}_{k}}{\partial \xi }$ and $\frac{\partial {\varvec{F}}_{k}}{\partial \zeta }$ are linearly independent.

2.3 Notations of Partitions and Function Spaces for the Parametric Domain and Physical Domain

Notations in this subsection will make a more rigorous description of the whole paper. The objection of this subsection is to introduce the partitions and function spaces for the whole parametric domain and the physical domain. However, it will inevitably involve the partitions and function spaces over one parametric patch, which will make this subsection very tedious. Thus, we will omit the details when introducing the partitions and function spaces over one parametric patch which have already been defined in detail in [3]. In addition, we will make use of the classical Lebesgue space $L^{2}(\cdot )$ endowed with the norm $||\cdot ||_{L^{2}(\cdot )}$ throughout the analysis. The classical Hilbert spaces $H^{l}(\cdot )$, endowed with the semi-norms $\{|\cdot |_{H^{j}(\cdot )}\}_{j=0}^{l}$ and norm $||\cdot ||_{H^{l}(\cdot )}$ with $l$ a non-negative integer, is also used.

We denote the parametric domain by $D$. Without loss of generality, we simply assume in the following analysis that $D$ is a convex polygonal domain directly composed of all the parametric patches $\{D_{k}\}_{k = 1}^{{\mathfrak {M}}}$. In our method, each parametric patch $D_{k}$ will be further partitioned into non-overlapping elements $\{d_{k,i}\}$, see Fig. 2 for the description of $D$.

Fig. 2

Description of $D$. Each parametric patch, e.g., $D_{k}$ (green), is partitioned into non-overlapping elements. (Color figure online)

For each $k\in \{1,\ldots ,{\mathfrak {M}}\}$, we denote by ${\mathcal {T}}_{k}$ the partition of the parametric patch $D_{k}$ into 2-dimensional elements $\{d\}$ such that $\overline{D_{k}} = \cup _{d\in {\mathcal {T}}_{k}} \overline{d}$ and $d_{1}\cap d_{2} = \emptyset , \forall d_{1}, d_{2}\in {\mathcal {T}}_{k}$. The “bent” Sobolev space [3] of order $l \in {\mathbb {N}}$ over one parametric patch $D_{k}$ associated with ${\mathcal {T}}_{k}$ is denoted by ${\mathbb {H}}^{l,k}({\mathcal {T}}_{k})$, which is a well-defined Hilbert space over $D_{k}$, endowed with its semi-norms and norm [3]. The restriction, denoted by ${\mathbb {H}}^{l,k}(\breve{d})$, of ${\mathbb {H}}^{l,k}({\mathcal {T}}_{k})$ to a given support extension [3] $\breve{d}$ of $d$ is also defined in [3], and it is useful in the local approximation property (3.29) in Sect. 3.5.

Now we can give the partition, denoted by ${\mathcal {T}}$, of the whole parametric domain $D$ for the multi-patch case as follows:

$$\begin{aligned} {\mathcal {T}} := \bigcup _{k = 1}^{{\mathfrak {M}}} \{d\in {\mathcal {T}}_{k}\}. \end{aligned}$$

The set of all the edges of ${\mathcal {T}}$ is given by $\Gamma _{{\mathcal {T}}} := \big (\cup _{d_{1},d_{2}\in {\mathcal {T}}}(\partial d_{1} \cap \partial d_{2})\big ) \bigcup \big ( \cup _{d \in {\mathcal {T}}}(\partial d \cap \partial D) \big )$. We further introduce a subset of interest, denoted by $\Gamma _{D}$, of $\Gamma _{{\mathcal {T}}}$ as follows:

$$\begin{aligned} \Gamma _{D} := \left\{ e\in \Gamma _{{\mathcal {T}}} : e \subset \cup _{k=1}^{{\mathfrak {M}}}(\partial D_{k}) \right\} . \end{aligned}$$

Let $h_{d}, h_{e}$ denote the diameter of the element $d\in {\mathcal {T}}$ and the length of the edge $e\in \Gamma _{D}$, respectively. The global mesh size of the parametric domain is defined as

$$\begin{aligned} h := \max \left\{ h_{d}: d \in {\mathcal {T}} \right\} . \end{aligned}$$

The partition ${\mathcal {T}}$ is assumed to be shape regular, that is, the ratio between the smallest edge of $d\in {\mathcal {T}}$ and its diameter $h_{d}$ is bounded, uniformly with respect to $d$ and $h$.

With the notations above, we can now define our “broken” spaces on parametric domain $D$ associated with ${\mathcal {T}}$ as follows:

$$\begin{aligned} {\mathbb {H}}^{l}({\mathcal {T}}) := \big \{ v\in L^{2}(D) : v|_{D_{k}} \in {\mathbb {H}}^{l,k}({\mathcal {T}}_{k}), k = 1,\ldots ,{\mathfrak {M}} \big \}, \end{aligned}$$

and simply define its semi-norms $\{|v|_{{\mathbb {H}}^{j}({\mathcal {T}})}\}_{j=0}^{l}$ and norm $||v||_{{\mathbb {H}}^{l}({\mathcal {T}})}$ from those of ${\mathbb {H}}^{l,k}({\mathcal {T}}_{k})$ by a summation over the index $k\in \{1,\ldots ,{\mathfrak {M}}\}$. In particular, we let ${\mathbb {L}}^{2}({\mathcal {T}}) := {\mathbb {H}}^{0}({\mathcal {T}})$.

Next, we introduce the partition of the physical domain ${\mathfrak {D}}$, followed by an introduction of our “broken” spaces on the physical domain ${\mathfrak {D}}$. For each $k \in \{1, \ldots , {\mathfrak {M}}\}$, each element $d\in {\mathcal {T}}_{k}$ is mapped into an element ${\mathfrak {d}} = \mathbf{F }_{k}(d) :=\{ \mathbf{F }_{k}(\xi , \zeta ) : (\xi , \zeta )\in d \}.$ Thus, we could naturally introduce the partition, denoted by ${\mathfrak {T}}_{k}$, of the physical patch ${\mathfrak {D}}_{k}$ as ${\mathfrak {T}}_{k}:= \{ {\mathfrak {d}} = \mathbf{F }_{k}(d) : d\in {\mathcal {T}}_{k} \}.$ Analogously, let ${\mathfrak {T}}$ denote the partition of the whole physical domain ${\mathfrak {D}}$ such that

$$\begin{aligned} {\mathfrak {T}} := \cup _{k = 1}^{{\mathfrak {M}}} \{{\mathfrak {d}} \in {\mathfrak {T}}_{k}\}. \end{aligned}$$

The set of all the edges of ${\mathfrak {T}}$ is given by $\Gamma _{{\mathfrak {T}}} := \big (\cup _{{\mathfrak {d}}_{1},{\mathfrak {d}}_{2}\in {\mathfrak {T}}}(\partial {\mathfrak {d}}_{1} \cap \partial {\mathfrak {d}}_{2})\big ) \bigcup \big ( \cup _{{\mathfrak {d}} \in {\mathfrak {T}}}(\partial {\mathfrak {d}} \cap \partial {\mathfrak {D}}) \big )$. We define the subset of interest of $\Gamma _{{\mathfrak {T}}}$ by

$$\begin{aligned} \Gamma _{{\mathfrak {D}}} := \{ \widetilde{e}\in \Gamma _{{\mathfrak {T}}} : \widetilde{e} \subset \cup _{k=1}^{{\mathfrak {M}}}(\partial {\mathfrak {D}}_{k}) \}, \end{aligned}$$

and again let $h_{\widetilde{e}}$ denote the length of the curved edge $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$.

Let ${\mathcal {H}}^{l,k}({\mathfrak {T}}_{k})$ be the push-forward of the space ${\mathbb {H}}^{l,k}({\mathcal {T}}_{k})$ through the geometric mapping $\mathbf{F }_{k}$, thus it is given by ${\mathcal {H}}^{l,k}({\mathfrak {T}}_{k}) := \big \{ \tilde{v} : \tilde{v} = v \circ \mathbf{F }_{k}^{-1}, v\in {\mathbb {H}}^{l,k}({\mathcal {T}}_{k}) \big \}, k = 1,\ldots ,{\mathfrak {M}}.$ We now could define our “broken” spaces on the physical domain ${\mathfrak {D}}$ associated with ${\mathfrak {T}}$ as follows:

$$\begin{aligned} {\mathcal {H}}^{l}({\mathfrak {T}}) := \big \{ \tilde{v}\in L^{2}({\mathfrak {D}}) : \tilde{v}|_{{\mathfrak {D}}_{k}} \in {\mathcal {H}}^{l,k}({\mathfrak {T}}_{k}), \quad k = 1,\ldots ,{\mathfrak {M}} \big \}, \end{aligned}$$

and we let ${\mathcal {L}}^{2}({\mathfrak {T}}) := {\mathcal {H}}^{0}({\mathfrak {T}})$.

Further, we naturally define the ${\mathcal {L}}^{2}$-inner product for piecewise functions $\widetilde{w_{1}}, \widetilde{w_{2}}$ and piecewise vector-valued functions ${\widetilde{\varvec{w_{1}}}}, {\widetilde{\varvec{w_{2}}}}$ in ${\mathcal {L}}^{2}({\mathfrak {T}})$ as follows:

$$\begin{aligned} (\widetilde{w_{1}}, \widetilde{w_{2}})_{{\mathfrak {T}}} := \sum _{{\mathfrak {d}}\in {\mathfrak {T}}} \int _{{\mathfrak {d}}} \widetilde{w_{1}} \ \widetilde{w_{2}} \, {\mathrm {d}}A,\ \ \ ({\widetilde{\varvec{w_{1}}}}, {\widetilde{\varvec{w_{2}}}})_{{\mathfrak {T}}} := \sum _{{\mathfrak {d}}\in {\mathfrak {T}}} \int _{{\mathfrak {d}}} {\widetilde{\varvec{w_{1}}}} \cdot {\widetilde{\varvec{w_{2}}}} \, {\mathrm {d}}A, \end{aligned}$$

from which we could define the ${\mathcal {L}}^{2}$ norm as $||\cdot ||_{{\mathcal {L}}^{2}({\mathfrak {T}})} : = \sqrt{(\cdot ,\cdot )_{{\mathfrak {T}}}}$. In addition, for any set $\gamma _{{\mathfrak {D}}} \subset \Gamma _{{\mathfrak {D}}}$, the ${\mathcal {L}}^{2}$-inner product for piecewise vector-valued functions ${\widetilde{\varvec{w_{1}}}}, {\widetilde{\varvec{w_{2}}}}$ over $\gamma _{{\mathfrak {D}}}$ is defined by

$$\begin{aligned} \langle {\widetilde{\varvec{w_{1}}}}, {\widetilde{\varvec{w_{2}}}} \rangle _{\gamma _{{\mathfrak {D}}}} := \sum _{\tilde{e}\in \gamma _{{\mathfrak {D}}}} \int _{\tilde{e}} {\widetilde{\varvec{w_{1}}}} \cdot {\widetilde{\varvec{w_{2}}}} \, {\mathrm {d}}\lambda . \end{aligned}$$

2.4 Notations for DG Method

In this subsection, we introduce some trace operators that will appear in our formulation in next section.

For $\widetilde{q} \in {\mathcal {H}}^{1}({\mathfrak {T}})$, we define the average $\{\widetilde{q}\}$ and the jump $[\![\widetilde{q}]\!]$ of $\tilde{q}$ on $\Gamma _{{\mathfrak {D}}}$ as follows. Let $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$ be an interior edge shared by elements ${\mathfrak {d}}_{\delta (1)} \in {\mathfrak {T}}$ and ${\mathfrak {d}}_{\delta (2)} \in {\mathfrak {T}}$. Let $\varvec{\widetilde{n}_{i}}$ denote the unit normal vectors of $\mathfrak {D}$ on $\widetilde{e}$ pointing exterior to $\mathfrak {d}_{\delta (i)}$, $i = 1,2$. With $\widetilde{q}_{i} :=\widetilde{q}|_{\partial {\mathfrak {d}}_{\delta (i)}}$, we set

$$\begin{aligned} \{\widetilde{q}\} = \frac{1}{2}(\widetilde{q}_{1} + \widetilde{q}_{2}), \ \ [\![ \widetilde{q} ]\!] = \widetilde{q}_{1} {\widetilde{\varvec{n}}_{1}} + \widetilde{q}_{2} {\widetilde{\varvec{n}}_{2}}\ \ \ {\mathrm {on}} \ \widetilde{e}. \end{aligned}$$

For ${\widetilde{\varvec{\varphi }}} \in ({\mathcal {H}}^{1}({\mathfrak {T}}))^{3}$, we define ${\widetilde{\varvec{\varphi }}_{1}}$ and ${\widetilde{\varvec{\varphi }}_{2}}$ analogously and set

$$\begin{aligned} \{ {\widetilde{\varvec{\varphi }}} \} = \frac{1}{2}({\widetilde{\varvec{\varphi }}_{1}} + {\widetilde{\varvec{\varphi }}_{2}}), \ \ [\![ {\widetilde{\varvec{\varphi }}} ]\!] = {\widetilde{\varvec{\varphi }}_{1}} \cdot {\widetilde{\varvec{n}}_{1}} + {\widetilde{\varvec{\varphi }}_{2}} \cdot {\widetilde{\varvec{n}}_{2}} \ \ \ {\mathrm {on}} \ \widetilde{e}. \end{aligned}$$

Notice that the jump $\,[\![\widetilde{q}]\!]$ of the scalar function $\widetilde{q}$ is a vector-valued function, and the jump$\,[\![ {\widetilde{\varvec{\varphi }}} ]\!]$ of the vector-valued function ${\widetilde{\varvec{\varphi }}}$ is a scalar quantity. The advantage of these definitions is that they do not depend on assigning an ordering to the elements ${\mathfrak {d}}_{\delta (i)}$.

For a boundary edge $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$, each $\widetilde{q} \in {\mathcal {H}}^{1}({\mathfrak {T}})$ and ${\widetilde{\varvec{\varphi }}} \in ({\mathcal {H}}^{1}({\mathfrak {T}}))^{3}$ has a uniquely defined restriction on $\widetilde{e}$:

$$\begin{aligned}{}[\![\widetilde{q}]\!] = \widetilde{q}\ {\widetilde{\varvec{n}}}, \ \ \{ {\widetilde{\varvec{\varphi }}} \} = {\widetilde{\varvec{\varphi }}} \ \ \ {\mathrm {on}} \ \ \widetilde{e}, \end{aligned}$$

where ${\widetilde{\varvec{n}}}$ is the outward unit normal of ${\mathfrak {D}}$ along $\widetilde{e}$.

2.5 Formulation of the Problem

In this paper, we display the most distinctive features of our method using as simple a setting as possible. Thus, when formulating the problem, we just consider the model problem

$$\begin{aligned} \left\{ \begin{array}{llll} -\Delta _{s}\widetilde{u} &{} = &{} \widetilde{f} &{}\quad \text {in}\, {\mathfrak {D}}, \\ \widetilde{u} &{} = &{} 0 &{}\quad \text {on}\,\partial {\mathfrak {D}}, \end{array} \right. \end{aligned}$$

(2.8)

with $\widetilde{f} : {\mathfrak {D}} \mapsto {\mathbb {R}}$ a given function.

Assume $\widetilde{f}\in {\mathcal {L}}^{2}({\mathfrak {T}})$, we begin by presenting, without derivation, the weak formulation of (2.8) in the “broken” function spaces ${\mathcal {H}}^{2}({\mathfrak {T}})$ as follows: Find $\widetilde{u}\in {\mathcal {H}}^{2}({\mathfrak {T}})$ such that

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{u}, \widetilde{v}) = ( \widetilde{f}, \widetilde{v})_{{\mathfrak {T}}} \quad \forall \widetilde{v}\in {\mathcal {H}}^{2}({\mathfrak {T}}), \end{aligned}$$

(2.9)

with

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{u}, \widetilde{v}) = ( \nabla _{s} \widetilde{u}, \nabla _{s} \widetilde{v} )_{{\mathfrak {T}}} - \langle \{\nabla _{s} \widetilde{u}\}, [\![\widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}} - \langle [\![\widetilde{u}]\!], \{\nabla _{s} \widetilde{v}\} \rangle _{\Gamma _{{\mathfrak {D}}}} + \langle \widetilde{\mu } [\![\widetilde{u}]\!], [\![\widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}}, \end{aligned}$$

(2.10)

where the penalty weighting function $\widetilde{\mu } :\Gamma _{{\mathfrak {D}}} \mapsto {\mathbb {R}}$ is given by $\eta _{\widetilde{e}}h_{\widetilde{e}}^{-1}$ on each $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$ with $\eta _{\widetilde{e}}$ a positive number.

Remark

$\eta _{\widetilde{e}}$ should be taken large enough to guarantee the stability of ${\mathfrak {B}}_{{\mathfrak {T}}}$, which will be discussed in Sect. 3.4 and indicated when presenting numerical examples. We call ${\mathfrak {B}}_{{\mathfrak {T}}}$ in (2.10) the primal form. Additionally, we omit the derivation of the weak formulation, and point out that one can easily achieve it using Lemma 2.1 and the notations introduced in Sect. 2.4.

Next, we introduce the DG finite element space associated with ${\mathfrak {T}}$. In the standard IGA, basis functions over one physical patch ${\mathfrak {D}}_{k}$ are defined from $\{R_{i,j}^{p_{k},q_{k}}(\xi , \zeta )\}_{i=1,j=1}^{m_{k},n_{k}}$ by considering a composition with the inverse of the geometrical mapping $\mathbf{F }_{k}$. In our multi-patch case, we could therefore define the IGA basis functions over ${\mathfrak {D}}$, denoted by $\{ \widetilde{R}_{i,j}^{k}(x,y,z) \}_{k = 1,i=1,j=1}^{{\mathfrak {M}},m_{k},n_{k}}$, as follows: For each $k\in \{1,\ldots ,{\mathfrak {M}}\},\ i\in \{1,\ldots , m_{k}\},\ j\in \{1,\ldots ,n_{k}\}$,

$$\begin{aligned} \widetilde{R}_{i,j}^{k}(x,y,z) := \left\{ \begin{array}{ll} R_{i,j}^{p_{k},q_{k}}\circ \mathbf{F }_{k}^{-1}(x,y,z) &{} \quad \forall \, (x,y,z)\in {\mathfrak {D}}_{k},\\ 0 &{} \quad {\mathrm {otherwise}}. \end{array} \right. \end{aligned}$$

Since the employed NURBS basis functions are continuous within each parametric patch, discontinuities of the IGA basis functions over ${\mathfrak {D}}$ could only appear across the patch interfaces. The DG finite element space is defined as follows:

$$\begin{aligned} \widetilde{ W_{h} } = \Big \{ \widetilde{ w }(x,y,z) \in {\mathcal {L}}^{2}({\mathfrak {T}}) : \widetilde{ w }\big |_{{\mathfrak {D}}_{k}} \in {\mathrm {span}}\big \{\, \cup _{i=1}^{m_{k}} \cup _{j=1}^{n_{k}} \widetilde{R}_{i,j}^{k}(x,y,z)\, \big \}, \ k = 1,\ldots ,{\mathfrak {M}} \Big \}. \end{aligned}$$

(2.11)

With the notations above, we present our DG scheme for problem (2.8) as follows: Find $\widetilde{u_{h}} \in \widetilde{ W_{h} }$ such that

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{u_{h}}, \widetilde{v}) = ( \widetilde{f}, \widetilde{v})_{{\mathfrak {T}}}\quad \forall \widetilde{v} \in \widetilde{ W_{h} }, \end{aligned}$$

(2.12)

which is called the primal formulation. Let $\widetilde{u}$ solve the boundary value problem (2.8), then we have

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{u}, \widetilde{v}) = ( \widetilde{f}, \widetilde{v})_{{\mathfrak {T}}}\quad \forall \widetilde{v} \in \widetilde{ W_{h} }. \end{aligned}$$

(2.13)

We remark that (2.12) is consistent, which is defined to mean that (2.13) holds, at least for all $\widetilde{v} \in \widetilde{ W_{h} }$. In view of (2.12) it means that the Galerkin orthogonality holds

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{u} - \widetilde{u_{h}}, \widetilde{v}) = 0 \quad \forall \,\widetilde{v} \in \widetilde{ W_{h} }. \end{aligned}$$

(2.14)

Although we have only considered the model problem (2.8) with homogeneous Dirichlet boundary conditions, extension to more general boundary conditions can be easily carried out.

3 Theoretical Analysis

The objective of this section is to perform theoretical analysis for our method. We discuss separately the boundedness and stability of the primal form ${\mathfrak {B}}_{{\mathfrak {T}}}$ with respect to an appropriate norm which will be defined later. We will then be ready to present an error analysis incorporating with the approximation properties in the parameter space.

Let $p$ denote the lowest degree of the bivariate NURBS basis functions for geometrical mappings, thus it is given by

$$\begin{aligned} p := \min _{k=1}^{{\mathfrak {M}}} \{p^{k}, q^{k}\}. \end{aligned}$$

(3.1)

Here, we assume that $p\ge 3$ in this section. $c$ and $C$ denote generic positive constants which are independent of $h$.

3.1 Properties of the Geometrical Mappings

Denoting by $J_{k}$ the Jacobian matrix of ${\varvec{F}}_{k}$, and by $G_{k}$ the matrix whose elements are the coefficients of the first fundamental form of ${\mathfrak {D}}_{k}$, $k =1, \ldots , {\mathfrak {M}}$, we have

$$\begin{aligned} J_{k} = \Big [\frac{\partial {\varvec{F}}_{k}}{\partial \xi }, \frac{\partial {\varvec{F}}_{k}}{\partial \zeta }\Big ]_{3\times 2}, \ \ G_{k} = J_{k}^{T} J_{k} \ \ \ k =1, \ldots , {\mathfrak {M}}. \end{aligned}$$

(3.2)

With the assumption (2.7), we summarize some properties which will be used below in the following results.

Lemma 3.1

Assume (2.7) holds, we claim that there exist positive constants $c$ and $C$ such that

$$\begin{aligned} \begin{array}{rcl} c \le \Big \Vert \frac{\partial {\varvec{F}}_{k}}{\partial \xi }\Big \Vert \le C,\ c \le \Big \Vert \frac{\partial {\varvec{F}}_{k}}{\partial \zeta }\Big \Vert \le C &{} &{}\quad {\mathrm {on}}\ d\in {\mathcal {T}}_{k},\ k =1, \ldots , {\mathfrak {M}},\\ c \le ||J_{k}||_{L^{\infty }(d)} \le C &{} &{}\quad \forall d\in {\mathcal {T}}_{k},\ k =1, \ldots , {\mathfrak {M}},\\ ||(J_{k})_{\xi }||_{L^{\infty }(d)} \le C, \ ||(J_{k})_{\zeta }||_{L^{\infty }(d)} \le C &{} &{}\quad \forall d\in {\mathcal {T}}_{k},\ k =1, \ldots , {\mathfrak {M}},\\ c \le det(G_{k}) \le C &{} &{}\quad {\mathrm {on}}\ d\in {\mathcal {T}}_{k},\ k =1, \ldots , {\mathfrak {M}},\\ ||G_{k}^{-1}||_{L^{\infty }(d)} \le C &{} &{}\quad \forall d\in {\mathcal {T}}_{k},\ k =1, \ldots , {\mathfrak {M}},\\ \end{array} \end{aligned}$$

(3.3)

where $||\cdot ||$ denotes the 2-norm of a vector and $c\mathrm{,}\, C$ depend only on the dimensions of the related matrices. (Each row of the inequalities in (3.3) is referred to as the $i^{th}$ row, $i = 1,\ldots , 5$.)

Proof

First, for bivariate NURBS basis functions of bi-degree $(p,q)$ with $p,q\ge 3$, the geometrical mapping ${\varvec{F}}_{k}$ is then ${\mathrm {C}}^{2}$-continuous on $D_{k}$, $k = 1,\ldots , {\mathfrak {M}}$. Thus, the $3^{th}$ row, as well as the upper bound estimates in the $1^{th}$ row and the $2^{th}$ row, obviously hold. We claim that the lower bound estimates in the $1^{th}$ row also hold, otherwise, there will exist a $k_{0}\in \{1,\ldots ,{\mathfrak {M}}\}$ and $(\xi _{0},\zeta _{0})\in d \in {\mathcal {T}}_{k_{0}}$ such that $\Big \Vert \frac{\partial {\varvec{F}}_{k_{0}}}{\partial \xi }(\xi _{0},\zeta _{0})\Big \Vert = 0$ or $\Big \Vert \frac{\partial {\varvec{F}}_{k_{0}}}{\partial \zeta }(\xi _{0},\zeta _{0})\Big \Vert = 0$, i.e., $\frac{\partial {\varvec{F}}_{k_{0}}}{\partial \xi }(\xi _{0},\zeta _{0}) = {\varvec{0}}$ or $\frac{\partial {\varvec{F}}_{k_{0}}}{\partial \zeta }(\xi _{0},\zeta _{0}) = {\varvec{0}}$, which are in conflict with (2.7). The lower bound estimations in the $2^{th}$ row follows directly from the lower bound estimations in the $1^{th}$ row, incorporating with the equivalency of vector norms and matrix norms. $c, C$ are positive constants depending on the dimensions of the related vectors and matrices.

Next, we consider about properties of matrices $\{G_{k}\}_{k = 1}^{{\mathfrak {M}}}$. For each $k \in \{1, \ldots , {\mathfrak {M}}\}$, since vector $\frac{\partial {\varvec{F}}_{k}}{\partial \xi }$ and $\frac{\partial {\varvec{F}}_{k}}{\partial \zeta }$ are linearly independent, $G_{k}$ is then symmetric positive definite according to (3.2). Thus, there exists an orthogonal matrix $T_{k}$ such that

$$\begin{aligned} G_{k} = T_{k}^{T} \left( \begin{array}{cc} \lambda _{k,1} &{} 0 \\ 0 &{} \lambda _{k,2} \end{array} \right) T_{k}, \end{aligned}$$

(3.4)

where $\lambda _{k,1}, \lambda _{k,2}$ are the eigenvalues of $G_{k}$ which satisfy

$$\begin{aligned} c \le \lambda _{k,2} \le \lambda _{k,1}\quad \forall \, (\xi , \zeta )\in d \in {\mathcal {T}}_{k}, \end{aligned}$$

(3.5)

with $c>0$ a positive constant independent of $k$. With the equivalency and consistency of matrix norms, we have

$$\begin{aligned} \lambda _{k,1} = ||G_{k}||_{2} \le C ||G_{k}||_{L^{\infty }(d)} \le C ||J_{k}||_{L^{\infty }(d)}^{2} \quad \forall \, (\xi , \zeta )\in d \in {\mathcal {T}}_{k}, \end{aligned}$$

where we use (3.2) for the second inequality and $C$ is a positive constant depending only on the dimensions of the related matrices. From the $2^{th}$ row, we know that

$$\begin{aligned} \lambda _{k,1} \le C \ \ \ \ \forall \, (\xi , \zeta )\in d \in {\mathcal {T}}_{k}, \end{aligned}$$

(3.6)

where $C$ is a positive constant depending only on the dimensions of the related matrices. Thus, it follows that the $4^{th}$ row holds from (3.5) and (3.6) since the fact that

$$\begin{aligned} \lambda _{k,2} ^{2} \le det(G_{k}) \le \lambda _{k,1} ^{2}\ \ \ \ \forall \, (\xi , \zeta )\in d \in {\mathcal {T}}_{k}. \end{aligned}$$

Actually, the $5^{th}$ row also holds from the fact that

$$\begin{aligned} ||G_{k}^{-1}||_{2} = \frac{1}{\lambda _{k,2}} \le \frac{1}{C_{0}}\ \ \ \ \forall \, (\xi , \zeta )\in d \in {\mathcal {T}}_{k}, \end{aligned}$$

and a combination with the consistency of matrix norms. $\square$

3.2 Mesh-dependent Semi-norms and Norms

In this subsection, we define appropriate mesh-dependent semi-norms and norms, associated with the DG method, for spaces of functions on the parametric domain $D$ and the physical domain ${\mathfrak {D}}$, respectively.

We define trace operators, average and jump, for functions in ${\mathbb {H}}^{1}({\mathcal {T}})$ in the same way as we did in Sect. 2.4 and adopt the same notations, just noting that the related valued-functions and outward unit normal vectors here are 2-dimensional. Just like (2.11), we define the space of functions on $D$ as

$$\begin{aligned} W_{h} = \Big \{ w(\xi , \zeta ) \in {\mathbb {L}}^{2}({\mathcal {T}}) : w\big |_{ D_{k} } \in {\mathrm {span}}\big \{\ \cup _{i=1}^{m_{k}} \cup _{j=1}^{n_{k}} R_{i,j}^{k}(\xi , \zeta )\ \big \}, \ k = 1,\ldots ,{\mathfrak {M}} \Big \}. \end{aligned}$$

Further on, we let

$$\begin{aligned} W(h) = W_{h} + {\mathbb {H}}^{2}({\mathcal {T}})\cap {\mathbb {H}}^{1}_{0}({\mathcal {T}}), \end{aligned}$$

(3.7)

where the subscript 0 of ${\mathbb {H}}^{1}_{0}({\mathcal {T}})$ means that functions in ${\mathbb {H}}^{1}_{0}({\mathcal {T}})$ vanish on $\partial D$. Now, we define the following semi-norms and the norm for $w\in W(h)$:

$$\begin{aligned} |w|^{2}_{1,{\mathcal {T}}}= & {} \sum _{d\in {\mathcal {T}}} |w|^{2}_{H^{1}(d)}, \ \ |w|^{2}_{*} = \sum _{e\in \Gamma _{D} } \int _{e} \frac{1}{h_{e}}\, [\![ w ]\!] \cdot [\![ w ]\!] \, {\mathrm {d}}\gamma ,\nonumber \\ |\!|\!|w|\!|\!|^{2}= & {} |w|^{2}_{1, {\mathcal {T}}} + \sum _{d\in {\mathcal {T}}} h_{d}^{2} |w|^{2}_{H^{2}(d)} + |w|^{2}_{*}. \end{aligned}$$

(3.8)

The norm (3.8) is the natural one for obtaining boundedness of the bilinear form ${\mathfrak {B}}_{{\mathfrak {T}}}$. On the other hand, the DG norm

$$\begin{aligned} ||w||_{DG, {\mathcal {T}}} = ( |w|^{2}_{1, {\mathcal {T}}} + |w|^{2}_{*} )^{1/2}, \end{aligned}$$

(3.9)

is the natural one for analyzing the stability. Restricted to $w\in W_{h}$, these two norms are equivalent, i.e., there exist positive constants $c$ and $C$, such that

$$\begin{aligned} c||w||_{DG, {\mathcal {T}}} \le |\!|\!|w|\!|\!| \le C||w||_{DG, {\mathcal {T}}}, \end{aligned}$$

(3.10)

as is evident from a local inverse inequality, see [25].

Next, we consider function space on ${\mathfrak {D}}$. We let

$$\begin{aligned} \widetilde{W(h)} = \widetilde{W_{h}} + {\mathcal {H}}^{2}({\mathfrak {T}})\cap {\mathcal {H}}^{1}_{0}({\mathfrak {T}}), \end{aligned}$$

and define the following norm for $\widetilde{w} \in \widetilde{W(h)}$:

$$\begin{aligned} ||\widetilde{w}||^{2}_{DG,{\mathfrak {T}}} = \sum _{{\mathfrak {d}}\in {\mathfrak {T}}} |\widetilde{w}|_{1,{\mathfrak {d}}}^{2} + |\widetilde{w}|^{2}_{\bigstar }, \end{aligned}$$

(3.11)

where for each ${\mathfrak {d}}\in {\mathfrak {T}}$,

$$\begin{aligned} \begin{array}{cc} \displaystyle |\widetilde{w}|_{1,{\mathfrak {d}}}^{2} = \displaystyle \int _{{\mathfrak {d}}} \nabla _{s}\widetilde{w} \cdot \nabla _{s}\widetilde{w}\, {\mathrm {d}}A,&\displaystyle |\widetilde{w}|^{2}_{\bigstar } = \displaystyle \sum _{\widetilde{e} \in \Gamma _{{\mathfrak {D}}} } \int _{\widetilde{e}} \frac{1}{h_{\widetilde{e}}}\, [\![ \widetilde{w}]\!] \cdot [\![ \widetilde{w}]\!] \, {\mathrm {d}} \lambda . \end{array} \end{aligned}$$

Now, we show in a lemma the relationship between $||\cdot ||_{{\mathcal {L}}^{2}({\mathfrak {T}})}$ and $||\cdot ||_{{\mathbb {L}}^{2}({\mathcal {T}})}$, $||\cdot ||_{DG,{\mathfrak {T}}}$ and $||\cdot ||_{DG,{\mathcal {T}}}$. Since each patch of the physical domain has a parameter representation, a function in $\widetilde{W(h)}$ can also be regarded as a function in $W(h)$, i.e., a function on the parameter domain $D$.

Lemma 3.2

For all $\widetilde{w} \in \widetilde{W(h)}$, there exist positive constants $c$ and $C$ such that

$$\begin{aligned} c ||w||_{{\mathbb {L}}^{2}({\mathcal {T}})}\le & {} ||\widetilde{w}||_{{\mathcal {L}}^{2}({\mathfrak {T}})} \le C ||w||_{{\mathbb {L}}^{2}({\mathcal {T}})}, \end{aligned}$$

(3.12)

$$\begin{aligned} c||w||_{DG,{\mathcal {T}}}\le & {} ||\widetilde{w}||_{DG,{\mathfrak {T}}} \le C ||w||_{DG,{\mathcal {T}}}, \end{aligned}$$

(3.13)

where $w \in W(h)$ denotes the pull-back of $\widetilde{w}$ via the geometrical mappings.

Proof

For (3.12), a straightforward computation gives that for each ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$,

$$\begin{aligned} \int _{{\mathfrak {d}}} \widetilde{w} \widetilde{v}\ {\mathrm {d}}A = \iint _{ d } wv \ \sqrt{det(G_{k})}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta , \end{aligned}$$

(3.14)

where $d = \mathbf{F }_{k}^{-1}({\mathfrak {d}})$. Using the $4^{th}$ row in (3.3), (3.12) holds from (3.14) by a summation over ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$ and the index $k\in \{1,\ldots , {\mathfrak {M}}\}$.

For (3.13), we consider each term that appears in $||\cdot ||_{DG,{\mathfrak {T}}}$. On one hand, we claim that

$$\begin{aligned} c|w|_{1,{\mathcal {T}}}^{2} \le \sum _{{\mathfrak {d}}\in {\mathfrak {T}}} |\widetilde{w}|_{1,{\mathfrak {d}}}^{2} \le C |w|_{1,{\mathcal {T}}}^{2}, \end{aligned}$$

which will be proved in Sects. 3.3 and 3.4, see (3.16) and (3.28), respectively. On the other hand, we claim that

$$\begin{aligned} c|w|_{*}^{2} \le |\widetilde{w}|^{2}_{\bigstar } \le C |w|_{*}^{2}, \end{aligned}$$

which will also be proved in Sects. 3.3 and 3.4, see (3.20) and (3.26), respectively. Thus, (3.13) holds from the above two inequalities.

The proof is finished. $\square$

3.3 Boundedness

We now show as a lemma that the primal form ${\mathfrak {B}}_{{\mathfrak {T}}}$ in (2.10) is bounded with respect to the norm $|\!|\!|\cdot |\!|\!|$.

Lemma 3.3

[Boundedness]

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{w}, \widetilde{v}) \le C |\!|\!|w|\!|\!|\, |\!|\!|v|\!|\!| \ \ \ \ \ \forall \widetilde{w}, \widetilde{v} \in \widetilde{W(h)}, \end{aligned}$$

where $C$ is a positive constant independent of $h$ and $w, v\in W(h)$ denote the pull-back of $\widetilde{w}, \widetilde{v}$ via the geometrical mappings, respectively.

Proof

We show this by bounding each of the integral terms that appears in ${\mathfrak {B}}_{\mathfrak {T}}$. For the first term, a straightforward computation and scaling give that for each ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$ with $d = \mathbf{F }_{k}^{-1}({\mathfrak {d}})$,

$$\begin{aligned} \displaystyle \int _{{\mathfrak {d}}} \nabla _{s}\widetilde{w} \cdot \nabla _{s}\widetilde{v}\, {\mathrm {d}}A= & {} \displaystyle \iint _{ d} (\nabla w)^{T} G_{k}^{-1} (\nabla v) \ \sqrt{det(G_{k})}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta \nonumber \\\le & {} \displaystyle \iint _{ d} ||\nabla w|| \ ||G_{k}^{-1}||_{F} \ ||\nabla v|| \ \sqrt{det(G_{k})}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta \nonumber \\\le & {} \displaystyle C \iint _{ d} ||\nabla w|| \ ||G_{k}^{-1}||_{L^{\infty }(D_{k})} \ ||\nabla v|| \ \sqrt{det(G_{k})}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta ,\quad \quad \quad \end{aligned}$$

(3.15)

where the first inequality comes from the consistency of the Frobenius norm $||\cdot ||_{F}$ of a matrix, while the second one makes use of the equivalency of matrix norms. $C$ is a positive constant depending only on the dimension of the related matrix. After using the $4^{th}$ row, the $5^{th}$ row in (3.3) and Hölder’s inequality, we get

$$\begin{aligned} ( \nabla _{s}\widetilde{w} , \nabla _{s}\widetilde{v})_{{\mathfrak {T}} } \le C |w|_{1,{\mathcal {T}}}|v|_{1,{\mathcal {T}}}, \end{aligned}$$

(3.16)

from a simple summation over ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$ and the index $k\in \{1,\ldots , {\mathfrak {M}}\}$.

Next, we consider about the penalty term. Before doing that, let’s introduce some notations that will be used here. Let $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$ be an interior edge shared by elements ${\mathfrak {d}}_{\delta (1)}\in {\mathfrak {T}}$ and ${\mathfrak {d}}_{\delta (2)}\in {\mathfrak {T}}$. Let $\varvec{\widetilde{n}_{i}}$ denote the unit normal vectors on $\widetilde{e}$ pointing exterior to $\mathfrak {d}_{\delta (i)}$, $i = 1,2$. We know that there is a corresponding interior edge $e \in \Gamma _{ D }$ shared by $d_{\delta (1)}\in {\mathcal {T}}$ and $d_{\delta (2)}\in {\mathcal {T}}$. Similarly, we define the unit normal vectors ${\varvec{n_{1}}}$ and ${\varvec{n_{2}}}$ on $e$ pointing exterior to $d_{\delta (1)}$ and $d_{\delta (2)}$, respectively. With $\widetilde{w}_{i} = \widetilde{w}|_{\partial {\mathfrak {d}}_{\delta (i)}}$, $\widetilde{v}_{i} = \widetilde{v}|_{\partial {\mathfrak {d}}_{\delta (i)}}$ , $w_{i} = w|_{\partial d_{\delta (i)}}$, and $v_{i} = v|_{\partial d_{\delta (i)}}, i = 1,2$, we have

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{3} : &{} = &{} \displaystyle \int _{\widetilde{e}}\ \widetilde{\mu }\ [\![\widetilde{w}]\!]\cdot [\![\widetilde{v}]\!] \, {\mathrm {d}} \lambda \\ &{} = &{} \displaystyle \int _{\widetilde{e}}\ \widetilde{\mu }\ \left( \widetilde{w}_{1} {\widetilde{\varvec{n}}_{1}} + \widetilde{w}_{2} {\widetilde{\varvec{n}}_{2}}\right) \cdot \left( \widetilde{v}_{1} {\widetilde{\varvec{n}}_{1}} + \widetilde{v}_{2} {\widetilde{\varvec{n}}_{2}}\right) \, {\mathrm {d}} \lambda . \end{array} \end{aligned}$$

(3.17)

Noting that we have ${\widetilde{\varvec{n}}_{2}} = - {\widetilde{\varvec{n}}_{1}}$ from the first-order geometric continuity of the physical domain ${\mathfrak {D}}$, we get

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{3} &{} = &{} \displaystyle \int _{\widetilde{e}}\ \widetilde{\mu }\ (\widetilde{w}_{1} - \widetilde{w}_{2} ) \left( \widetilde{v}_{1}- \widetilde{v}_{2} \right) \, {\mathrm {d}} \lambda \\ &{} = &{} \displaystyle \int _{ e }\ \widetilde{\mu }\ ( w_{1} - w_{2} ) (v_{1}- v_{2} ) \ L \, {\mathrm {d}} \gamma , \end{array} \end{aligned}$$

(3.18)

where $L$, the length scale, is $||\frac{\partial {\varvec{F}}_{i}}{\partial \xi }||$ or $||\frac{\partial {\varvec{F}}_{i}}{\partial \zeta }||,\ i \in \{\delta (1),\delta (2)\}$, on $e$. $\widetilde{\mu }$ should be written as $\widetilde{\mu }(\xi , \zeta )$ in the second equality in (3.18). Using the fact ${\varvec{n_{2}}} = -{\varvec{n_{1}}}$, we have

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{3} &{} = &{} \displaystyle \int _{ e }\ \widetilde{\mu }\ ( w_{1}{\varvec{n_{1}}} + w_{2}{\varvec{n_{2}}} ) \cdot (v_{1}{\varvec{n_{1}}} + v_{2}{\varvec{n_{2}}} )\ L \, {\mathrm {d}} \gamma \\ &{} = &{} \displaystyle \int _{ e }\ \widetilde{\mu }\ [\![ w ]\!] \cdot [\![ v ]\!] \ L \, {\mathrm {d}} \gamma \\ &{} = &{} \displaystyle \int _{ e }\ \eta _{\widetilde{e}} \frac{1}{h_{\widetilde{e}}}\ [\![ w ]\!] \cdot [\![ v ]\!] \ L \, {\mathrm {d}} \gamma . \end{array} \end{aligned}$$

(3.19)

For the case in which $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$ is a boundary edge, one can deal with it similarly. Then after using the $1^{th}$ and $2^{th}$ rows in (3.3), we can easily get

$$\begin{aligned} \langle \widetilde{\mu }\, [\![\widetilde{w}]\!], [\![\widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}} \le C ( \sup \limits _{\widetilde{e}} \eta _{\widetilde{e}}) |w|_{*}|v|_{*}, \end{aligned}$$

(3.20)

with $C$ a positive constant, from a summation over $\widetilde{e} \in \Gamma _{{\mathfrak {D}}}$. We have used the fact that $h_{\widetilde{e}} \sim {\mathcal {O}}(h_{e})$ to get (3.20).

At last, we bound the terms $\langle \{\nabla _{s}\widetilde{w}\} , [\![ \widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}}$ and $\langle [\![\widetilde{w}]\!] , \{ \nabla _{s}\widetilde{v} \} \rangle _{\Gamma _{{\mathfrak {D}}}}$. On one hand, for each $d\in {\mathcal {T}}$, if $w\in H^{2}(d)$ and $e$ is an edge of $d$, we can get from [21, eq. (2.5)] that,

$$\begin{aligned} \int _{e} ||\nabla w||^{2} \, {\mathrm {d}} \gamma \le C(h_{e}^{-1}|w|^{2}_{H^{1}(d)} + h_{e} |w|^{2}_{H^{2}(d)}), \end{aligned}$$

(3.21)

where $C$ depends only on the minimum angle of $d$. It follows that, for every $q\in L^{2}(e)$,

$$\begin{aligned} \int _{e} ||\nabla w||\ q \, {\mathrm {d}} \gamma \le C\left( |w|^{2}_{H^{1}(d)} + h_{e}^{2} |w|^{2}_{H^{2}(d)}\right) ^{1/2}\left( h_{e}^{-1}\int _{e} q^{2} d\gamma \right) ^{1/2}. \end{aligned}$$

On the other hand, a straightforward computation gives

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{2} : &{} = &{} \displaystyle \int _{\widetilde{e} } \nabla _{s}\widetilde{w} \cdot [\![ \widetilde{v}]\!] \, {\mathrm {d}} \lambda \\ &{} = &{} \displaystyle \int _{ e } \ (\nabla w )^{T} G_{k}^{-1} J_{k}^{T}\ {\widetilde{\varvec{n}}_{1}}\ (v_{1}-v_{2})\ L \, {\mathrm {d}} \gamma , \end{array} \end{aligned}$$

(3.22)

where $L$ again denotes the length scale as above, and ${\widetilde{\varvec{n}}_{1}}$ should be written as ${\widetilde{\varvec{n}}_{1}}(\xi ,\zeta )$. A simple scaling gives

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{2} &{} \le &{} \displaystyle \int _{ e } \ ||\nabla w ||\ || G_{k}^{-1}||_{F} ||J_{k}^{T}||_{F}\ |v_{1}-v_{2}|\ L \, {\mathrm {d}} \gamma \\ &{} \le &{} \displaystyle \int _{ e } \ ||\nabla w ||\ || G_{k}^{-1}||_{L^{\infty }(d)} ||J_{k}^{T}||_{L^{\infty }(d)}\ |v_{1}-v_{2}|\ L \, {\mathrm {d}} \gamma . \end{array} \end{aligned}$$

(3.23)

After we use the 1th, 2th and 5th rows in (3.3) and the fact that $|v_{1}-v_{2}|$ can be written as $||[\![ v ]\!]||$, it follows that,

$$\begin{aligned} \begin{array}{lll} \displaystyle I_{2} &{} \le &{} \displaystyle C \int _{ e } \ ||\nabla w || \ ||[\![ v ]\!] || \, {\mathrm {d}} \gamma \\ &{} \le &{} \displaystyle C\,\left( |w|^{2}_{H^{1}(d)} + h_{e}^{2} |w|^{2}_{H^{2}(d)}\right) ^{1/2}\left( h_{e}^{-1}\int _{e} [\![ v ]\!]^{2} \, {\mathrm {d}}\gamma \right) ^{1/2}. \end{array} \end{aligned}$$

(3.24)

This implies that

$$\begin{aligned} \begin{array}{lll} \displaystyle \langle \{\nabla _{s}\widetilde{w}\} , [\![ \widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}} &{} = &{} \sum \limits _{ \widetilde{e} \in \Gamma _{{\mathfrak {D}}} } \int _{\widetilde{e} } \{\nabla _{s}\widetilde{w}\} \cdot [\![\widetilde{v}]\!]\, {\mathrm {d}} \lambda \\ &{} \le &{} \displaystyle C \Big ( \sum _{d\in {\mathcal {T}}} (|w|^{2}_{H^{1}(d)} + h_{d}^{2} |w|^{2}_{H^{2}(d)}) \Big )^{1/2} \Big ( \sum _{e \in \Gamma _{D}} h_{e}^{-1}\int _{e} [\![ v ]\!]^{2}\, {\mathrm {d}}\gamma \Big )^{1/2}\\ &{} \le &{} \displaystyle C\, |\!|\!|w|\!|\!| \ |v|_{*}. \end{array} \end{aligned}$$

(3.25)

Similarly, $\langle [\![\widetilde{w}]\!] , \{ \nabla _{s}\widetilde{v} \} \rangle _{\Gamma _{{\mathfrak {D}}}} \le C |w|_{*}|\!|\!|v|\!|\!|$.

We have thus established the bound and the proof is finished. $\square$

3.4 Stability

Now we show that our method satisfies the stability condition in the following lemma.

Lemma 3.4

(Stability)

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{v}, \widetilde{v}) \ge C |\!|\!|v|\!|\!|^{2} \ \ \ \ \ \forall \, \widetilde{v} \in \widetilde{W_{h}}, \end{aligned}$$

where $C$ is a positive constant independent of $h$ and $v\in W_{h}$ denotes the pull-back of $\widetilde{v}$ via the geometrical mappings.

Proof

For the sake of convenience for writing, we let

$$\begin{aligned} b_{1}(\widetilde{v}, \widetilde{v}) := \left( \nabla _{s}\widetilde{v} , \nabla _{s}\widetilde{v}\right) _{{\mathfrak {T}} } , \quad b_{3}(\widetilde{v}, \widetilde{v}) := \langle \widetilde{\mu }\, [\![\widetilde{v}]\!], [\![\widetilde{v}]\!] \rangle _{\Gamma _{{\mathfrak {D}}}}, \end{aligned}$$

and $b_{2}(\widetilde{v}, \widetilde{v})$ gathers up all the remaining terms of the primal form ${\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{v}, \widetilde{v})$. Thus we can write ${\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{v}, \widetilde{v})$ as

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{v}, \widetilde{v}) = b_{1}(\widetilde{v}, \widetilde{v}) + b_{2}(\widetilde{v}, \widetilde{v}) + b_{3}(\widetilde{v}, \widetilde{v}). \end{aligned}$$

From the bounds we obtained in the last section, we know that $|b_{2}(\widetilde{v}, \widetilde{v})| \le C_{2} |\!|\!|v|\!|\!|\ |v|_{*}$ for some constant $C_{2}$ and all $\widetilde{v} \in \widetilde{W(h)}$.

When present, the penalty term $b_{3}(\widetilde{v}, \widetilde{v})$ contributes to the stability of our method. In view of (3.19) and noting $h_{\widetilde{e}} \sim {\mathcal {O}}(h_{e})$, we have

$$\begin{aligned} b_{3}(\widetilde{v}, \widetilde{v}) \ge C_{3} \eta _{0} |v|_{*}^{2} \ \ \ \ \forall \, \widetilde{v} \in \widetilde{W_{h}}, \end{aligned}$$

(3.26)

with $\eta _{0} := \inf \limits _{\widetilde{e}} \eta _{\widetilde{e}}$, from the 1th and 2th rows in (3.3).

For $b_{1}(\widetilde{v}, \widetilde{v})$, we get from (3.4) that for each ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$ with $d = \mathbf{F }_{k}^{-1}({\mathfrak {d}})$,

$$\begin{aligned} \begin{array}{lll} \displaystyle \int _{{\mathfrak {d}}} \nabla _{s}\widetilde{v} \cdot \nabla _{s}\widetilde{v}\, {\mathrm {d}}A &{} = &{} \iint _{ d} (\nabla v)^{T} G_{k}^{-1} (\nabla v) \ \sqrt{det(G_{k})}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta \\ &{} = &{} \displaystyle \iint _{d} (\nabla v)^{T} T_{k}^{T} \left( \begin{array}{cc} \lambda _{k,1}^{-1} &{} 0 \\ 0 &{} \lambda _{k,2}^{-1} \end{array} \right) T_{k}(\nabla v) \ \sqrt{\lambda _{k,1}\lambda _{k,2}}\, {\mathrm {d}}\xi {\mathrm {d}}\zeta \\ &{} \ge &{} \displaystyle (\lambda _{k,1}^{-1})^{2} \lambda _{k,2} \iint _{d} (\nabla v)^{T} T_{k}^{T} T_{k} (\nabla v) \, {\mathrm {d}}\xi {\mathrm {d}}\zeta \\ &{} \ge &{} \displaystyle C |v|^{2}_{H^{1}(d)}, \end{array} \end{aligned}$$

(3.27)

where the last inequality uses the boundedness of $\lambda _{k,1}, \lambda _{k,2}$ obtained in the $4^{th}$ and $5^{th}$ row in (3.3) and that $T_{k}$ is an orthogonal matrix. Then a simple summation over ${\mathfrak {d}}\in {\mathfrak {T}}_{k}$ and the index $k\in \{1,\ldots , {\mathfrak {M}}\}$ gives

$$\begin{aligned} b_{1}(\widetilde{v}, \widetilde{v}) \ge C_{1} |v|^{2}_{1,{\mathcal {T}}}\ \ \ \ \forall \, \widetilde{v} \in \widetilde{W_{h}}. \end{aligned}$$

(3.28)

Thus, we have

$$\begin{aligned} {\mathfrak {B}}_{h}(\widetilde{v}, \widetilde{v}) \ge C_{1} |v|^{2}_{1,{\mathcal {T}}} + C_{3} \eta _{0} |v|_{*}^{2} - C_{2} |\!|\!|v|\!|\!|\ |v|_{*} \ \ \ \ \forall \, \widetilde{v} \in \widetilde{W_{h}}. \end{aligned}$$

Now, we may use the arithmetic-geometric mean inequality $(2ab\le a^{2}\epsilon + b^{2}/\epsilon )$ on the last term, and then the equivalence (3.10) of the norms (3.8) and (3.9), to show that the stability holds for large enough $\eta _{0}$. $\square$

Obviously, the stability result in Lemma 3.4 implies the uniqueness since $|\!|\!|\cdot |\!|\!|$ is a norm. In addition, (2.12) leads to a linear equation system in finite dimension where the number of equations equals to the number of unknowns, existence is equivalent to uniqueness. Thus, the well-posedness of (2.12) is guaranteed.

3.5 Approximation in the Parameter Space

Let $\widetilde{u} \in \widetilde{W(h)}$ be the exact solution for problem (2.8), and $u \in W(h)$ denote the pull-back of $\widetilde{u}$ via the geometrical mappings. The last ingredient in the error analysis is a bound on the approximation error $|\!|\!|u - u_{I}|\!|\!|$, where $u_{I} \in W_{h}$ is a suitable approximation of $u$.

Here, we choose $u_{I}$ such that the restriction of $u_{I}$ to each parametric patch $D_{k}$ is the projection of $u$ defined in [3], eq. (3.27)]. For $u\in {\mathbb {H}}^{p+1}({\mathcal {T}})$, the space of bivariate NURBS functions with lowest degree $p$ gives the following local approximation property (see [3], Lemma 3.4]:

$$\begin{aligned} |u-u_{I}|_{H^{\nu }(d)} \le C h_{d}^{p+1-\nu } ||u||_{{\mathbb {H}}^{p+1}(\breve{d})} \ \ \ \ \forall \, d\in {\mathcal {T}}, \ \ \nu = 0,1,2, \end{aligned}$$

(3.29)

where $C$ is independent of $h$ and $\breve{d}$ is the support extension of $d$ as said in Sect. 2.3.

From the definition of the triple norm in (3.8), we have

$$\begin{aligned} |\!|\!|u-u_{I}|\!|\!|^{2} \le C \left( \ \sum _{d\in {\mathcal {T}}} |u\!-\!u_{I}|_{H^{1}(d)}^{2} \!+\! \sum _{d\in {\mathcal {T}}} h_{d}^{2} |u\!-\!u_{I}|_{H^{2}(d)}^{2} \!+\! \sum _{d\in {\mathcal {T}}} h_{d}^{-2} ||u\!-\!u_{I}||_{L^{2}(d)}^{2} \ \right) , \end{aligned}$$

(3.30)

where $C$ is independent of $h$ and we have used the fact (see, e.g., [21, eq. (2.4)]) that

$$\begin{aligned} ||v||^{2}_{L^{2}(e)} \le C\left( h_{e}^{-1}||v||^{2}_{L^{2}(d)} + h_{e}|v|^{2}_{H^{1}(d)} \right) \quad \forall v\in H^{1}(d), \end{aligned}$$

when dealing with the penalty term in (3.8). From (3.29) and (3.30), we can obtain that

$$\begin{aligned} |\!|\!|u-u_{I}|\!|\!|\le C h^{p} ||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.31)

with $C$ a positive constant independent of $h$. Similarly, we can repeat the above process and just omit handling terms containing second derivatives to get that

$$\begin{aligned} ||u-u_{I}||_{DG,{\mathcal {T}}} \le C h^{p} ||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.32)

with $C$ a positive constant independent of $h$. Note that we have used the fact that $h_{e},\, h_{d}\sim {\mathcal {O}}(h)$ when getting (3.31) and (3.32).

3.6 Error Estimates

We now present error estimates for our method using the properties of boundedness, stability and approximation discussed in Sect. 3.5. We will show our method converges with optimal order with respect to the DG norm $||\cdot ||_{DG, {\mathfrak {T}}}$ and ${\mathcal {L}}^{2}$ norm $||\cdot ||_{{\mathcal {L}}^{2}({\mathfrak {T}})}$.

Theorem 1

(Error estimates in $||\cdot ||_{DG, {\mathfrak {T}}}$ and $||\cdot ||_{{\mathcal {L}}^{2}({\mathfrak {T}})}$) Assume $\widetilde{u_{h}} \in \widetilde{W_{h}}$ is the numerical solution of problem (2.8) with a sufficient smooth exact solution $\widetilde{u}\in \widetilde{W(h)}$, $W_{h}$ is the corresponding bivariate NURBS function space with lowest degree $p$ defined in (3.1). Then we have

$$\begin{aligned} ||\widetilde{u } - \widetilde{u_{h}}||_{DG, {\mathfrak {T}}} \le C h^{p} ||u||_{H^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.33)

further, if $D$ is convex we obtain

$$\begin{aligned} ||\widetilde{u } - \widetilde{u_{h}}||_{{\mathcal {L}}^{2}({\mathfrak {T}})} \le C h^{p+1} ||u||_{H^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.34)

where $C$ is a positive constant independent of $h$ and $u$ denotes the pull-back of $\widetilde{u}$ via the geometrical mappings.

Proof

We denote by $u_{h}$ the pull-back of $\widetilde{u_{h}}$ via the geometrical mappings. Let $u_{I} \in W_{h}$, which satisfies (3.29), be an suitable approximation of $u$ as described in Sect. 3.5.

Using the stability in Lemma 3.4, Galerkin orthogonality (2.14), boundedness in Lemma 3.3 and the approximation property (3.31), we have

$$\begin{aligned} \begin{array}{ll} &{} c|\!|\!|u_{I}-u_{h}|\!|\!|^{2} \\ &{}\quad \le {\mathfrak {B}}_{{\mathfrak {T}}}\left( \widetilde{u_{I}}- \widetilde{u_{h}}, \widetilde{u_{I}}- \widetilde{u_{h}}\right) \\ &{}\quad = {\mathfrak {B}}_{{\mathfrak {T}}}\left( \widetilde{u_{I}}- \widetilde{u}, \widetilde{u_{I}}- \widetilde{u_{h}}\right) \\ &{}\quad \le C |\!|\!|u - u_{I} |\!|\!|\ |\!|\!|u_{I}-u_{h}|\!|\!| \\ &{}\quad \le C h^{p}||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})} |\!|\!|u_{I}-u_{h}|\!|\!|, \\ \end{array} \end{aligned}$$

where $\widetilde{u_{I}}$ denotes the push-forward of $u_{I}$ via the geometrical mappings. Thus, we have

$$\begin{aligned} |\!|\!| u_{I}-u_{h}|\!|\!| \le C h^{p}||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

and then

$$\begin{aligned} || u_{I}-u_{h}||_{DG, {\mathcal {T}}} \le C h^{p}||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

using (3.10). Hence, the triangle inequality gives

$$\begin{aligned} |\!|\!| u - u_{h} |\!|\!|\le & {} C h^{p}||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.35)

$$\begin{aligned} || u - u_{h} ||_{DG, {\mathcal {T}}}\le & {} C h^{p}||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

(3.36)

from (3.31) and (3.32), respectively. Thus, (3.33) follows from (3.36) incorporating with (3.13) in Lemma 3.2.

Next, we show that we can obtain the optimal error estimate in the ${\mathcal {L}}^{2}$ norm by using the standard duality argument. As usual, we define the auxiliary function $\widetilde{\psi }$ as the solution of the problem

$$\begin{aligned} - \Delta _{s} \chi = \widetilde{u } - \widetilde{u_{h}} \ \ \ in \ {\mathfrak {D}}, \ \ \ \ \chi = 0 \ \ \ on \ \partial {\mathfrak {D}}, \end{aligned}$$

(3.37)

and again let $\psi$ denote the pull-back of $\widetilde{\psi }$ via the geometrical mappings. In view of (2.13) and the fact that our primal form ${\mathfrak {B}}_{{\mathfrak {T}}}(\cdot , \cdot )$ is symmetric, we have

$$\begin{aligned} {\mathfrak {B}}_{{\mathfrak {T}}}(\widetilde{v}, \widetilde{\psi }) = (\widetilde{u} - \widetilde{u_{h}} , \widetilde{v})_{{\mathfrak {T}}} \quad \forall \, \widetilde{v}\in \widetilde{W_{h}}. \end{aligned}$$

(3.38)

We take $\psi _{I}$ to be a piecewise linear interpolation of $\psi$, and let $\widetilde{\psi _{I}}$ be the push-forward of $\psi _{I}$ via the geometrical mappings. Now, taking $\widetilde{v} = \widetilde{u } - \widetilde{u_{h}}$ in (3.38) and using Galerkin orthogonality (2.14), boundedness in Lemma 3.3 and the approximation property (3.31), we obtain

$$\begin{aligned} \begin{array}{ll} &{} c||u-u_{h}||^{2}_{{\mathbb {L}}^{2}({\mathcal {T}})} \\ &{}\quad \le ||\widetilde{u } - \widetilde{u_{h}}||^{2}_{{\mathcal {L}}^{2}({\mathfrak {T}})} \\ &{}\quad = {\mathfrak {B}}_{{\mathfrak {T}}}\left( \widetilde{u } - \widetilde{u_{h}}, \widetilde{\psi }\right) \\ &{}\quad = {\mathfrak {B}}_{{\mathfrak {T}}}\left( \widetilde{u } - \widetilde{u_{h}}, \widetilde{\psi } - \widetilde{\psi _{I}}\right) \\ &{}\quad \le C |\!|\!|u-u_{h}|\!|\!| \ |\!|\!|\psi - \psi _{I}|\!|\!| \\ &{}\quad \le C h ||\psi ||_{{\mathbb {H}}^{2}({\mathcal {T}})}|\!|\!|u-u_{h}|\!|\!|, \\ \end{array} \end{aligned}$$

(3.39)

where for the first inequality, we have used (3.12) in Lemma 3.2.

Further on, if $D$ is convex, the elliptic regularity gives $||\psi ||_{{\mathbb {H}}^{2}({\mathcal {T}})} \le C_{r} ||u-u_{h}||_{{\mathbb {L}}^{2}({\mathcal {T}})}$ with $C_{r}$ depending only on $D$. Thus, (3.39) gives

$$\begin{aligned} ||u-u_{h}||_{{\mathbb {L}}^{2}({\mathcal {T}})} \le C h |\!|\!|u-u_{h}|\!|\!|, \end{aligned}$$

(3.40)

with $C$ a positive constant independent of $h$. Using (3.12), we get from (3.35) the optimal estimation in the ${\mathcal {L}}^{2}$ norm:

$$\begin{aligned} ||\widetilde{u } - \widetilde{u_{h}}||_{{\mathcal {L}}^{2}({\mathfrak {T}})} \le C h^{p+1} ||u||_{{\mathbb {H}}^{p+1}({\mathcal {T}})}, \end{aligned}$$

with $C$ a positive constant independent of $h$. $\square$

4 Numerical Examples

In this section, we present some numerical results of our method. We first consider problem (2.8) on a disk and hemispherical surface with conforming patches. Next, cases with non-conforming patches and different degrees across patch interfaces will be considered on a cylindrical surface. Then we compute problem (2.8) on a graph to show that our method is not only capable for parametric surfaces but also capable for graphs of functions over flat regions. In addition, our method can be easily extended to linear elasticity problem. We present an example which is considered as a benchmark in IGA literatures at last.

For the parameter $\eta _{\widetilde{e}}$ mentioned in Sect. 2.5 and discussed in Sect. 3.4, it is easy to choose in the implementation since our scheme is not sensitive to it at all. We just need to take it large enough to guarantee the stability of the bilinear form ${\mathfrak {B}}_{{\mathfrak {T}}}$. As observed from numerical experiments, $\eta _{\widetilde{e}}$ can take the range from 1 to $10^6$, and we take $\eta _{\widetilde{e}} = 10$ in our following numerical tests.

In the following examples, each parameter patch $D_{k}\in \{D_{k}\}_{k = 1}^{{\mathfrak {M}}}$ is further uniformly partitioned by a $N_{k}\times N_{k}$ mesh. The bivariate NURBS functions over $D_{k}$ are taken to be of bi-degree $(p_{k}, p_{k})$. We let $N = (N_{1},\ldots , N_{k}, \ldots , N_{{\mathfrak {M}}})$ and ${\mathcal {P}} = (p_{1},\ldots , p_{k}, \ldots , p_{{\mathfrak {M}}})$, which will be specified in each example.

4.1 Poisson Equation

In this example, we consider the Poisson equation (2.8) on different domain ${\mathfrak {D}}$ with the right hand side $f$ given by the exact solution, which will be specified in each case for test.

Fig. 3

Unit disk. Partition of the domain ${\mathfrak {D}}$; Azimuth $= 0^{\circ }$, Elevation $= 90^{\circ }$

4.1.1 Unit Disk: Conforming Patches and Same Degree Across Patches

In this case, domain ${\mathfrak {D}} = \{(x,y,z) : x^{2}+y^{2} \le 1, z = 0 \}$ is divided into five patches, $\{{\mathfrak {D}}_{i}\}_{i = 1}^{5}$, as depicted in Fig. 3. The exact solution is taken as $\sin (x^{2} + y^{2})\pi$ for test.

We consider conforming patches in this case, i.e., meshes of different patches match across patch interfaces $(N_{i} = N_{j}, i \ne j)$. As for ${\mathcal {P}}$, here we take $p_{i} = p_{j}, i \ne j$.

Numerical results are presented with different $N$ and ${\mathcal {P}}$. Convergence results in the ${\mathcal {L}}^{2}$ norm are shown in Table 1, from which we can see that the ${\mathcal {L}}^{2}$ convergence rate is approximately $p+1$, where $p$ is the lowest degree in each case. The data verify the theoretical result achieved in last section.

Table 1 Unit disk. $L^{2}$ error and convergence order

4.1.2 Hemispherical Surface: Conforming Patches and Same Degree Across Patches, Singularities in Geometrical Mappings

In this case, we use four patches to represent the domain surface ${\mathfrak {D}} = \{(x,y,z) : x^{2}+y^{2} +z^{2} = 1, z \ge 0 \}$, see Fig. 4 for the partition. We take the exact solution as $u(x,y,z) = (\sin x + \sin y)z + \sin z$ for test.

Fig. 4

Hemispherical surface. Partition of the domain ${\mathfrak {D}}$; Azimuth $= 75^{\circ }$, Elevation $= 50^{\circ }$

Conforming patches are still employed in this case. The configurations of $N$ and ${\mathcal {P}}$ are the same with those described in last example in Sect. 4.1.1.

We show the convergence results in the ${\mathcal {L}}^{2}$ norm in Table 2. One point we wish to make clear for this instance is that each geometrical mapping has a singular point (the pole of the hemispherical surface). As a result, we can still achieve the optimal ${\mathcal {L}}^{2}$ convergence rate $p+1$ with $p$ the lowest degree in each case when the degrees of NURBS are low, for example, $p = 3$ and $p = 4$. However, convergence rates are not that good when $p$ gets higher due to the singular point of each geometrical mapping.

Table 2 Hemispherical surface. $L^{2}$ error and convergence order

4.1.3 Cylindrical Surface: Non-conforming Patches and Different Degrees Across Patches

The physical domain for this numerical example is a cylindrical surface, whose radius is 1 and height is 2, i.e.,

$$\begin{aligned} {\mathfrak {D}} = \left\{ (x, y, z) : x^{2} + y^{2} = 1, 0\le z \le 2 \right\} , \end{aligned}$$

which is uniformly partitioned into four patches, see Fig. 5. We test this example using the exact solution $u(x,y,z) = 10\cos ^{2}{x}\, \sin ^{2}{y}\, \sin {\frac{\pi z}{2}}$.

Fig. 5

Cylindrical surface. Partition of the domain ${\mathfrak {D}}$

In this case, we consider non-conforming patches , i.e., meshes of different patches do not match across patch interfaces ($N_{i} \ne N_{j}$ for adjacent patches), see Fig. 6 for non-conforming meshes of two patches. For a positive integer $n$, $N$ is taken to be $(n, 2n, n, 2n)$ in this example. As for ${\mathcal {P}}$, degrees differ from patch to patch, i.e., $p_{i} \ne p_{j}$ for adjacent patches. ${\mathcal {P}} = (2,3,2,3)$, ${\mathcal {P}} = (3,4,3,4)$ and ${\mathcal {P}} = (2,3,4,5)$ are considered here.

Fig. 6

Non-conforming meshes of two patches

Convergence results in the ${\mathcal {L}}^{2}$ norm are depicted in Table 3. As can be seen, again we obtain that the ${\mathcal {L}}^{2}$ convergence rate is approximately $p+1$ with $p$ the lowest degree in each case, which is the best according to our theoretical result. Thus, cases with non-conforming patches and different degrees across patch interfaces can be handled very well.

Table 3 Cylindrical surface. $L^{2}$ error and convergence order

4.2 Graph: Conforming Patches and Same Degree Across the Patches

Our method is not only capable for parametric surfaces but also capable for graphs over some plat regions. In this subsection, we consider the model problem (2.8) on graphs over a plat region. Conforming patches and the same degree across the patches are employed in the implementation.

The domain surface ${\mathfrak {D}}$ of this example is the graph of class ${\mathrm {C}}^{1,\alpha }$ given by

$$\begin{aligned} z(x,y) = (0.75-x^{2}-y^{2})_{+}^{1+\alpha }, \end{aligned}$$

(4.1)

over the flat domain $\Omega = [0,1]^{2}$ and $f$ is given by the exact solution which is taken for test as $u(x,y,z) = \sin {(\pi x)}\sin {(\pi y)}e^{1-z}.$

We wish to emphasize here that the essential feature for IGA is the isoparametric concept rather than the NURBS representation of the geometrical mapping. Actually, the geometrical mappings adopted in this example take the following form:

$$\begin{aligned} \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{c} x(\xi , \zeta ) \\ y(\xi , \zeta ) \\ z(x(\xi , \zeta ), y(\xi , \zeta )) \end{array} \right) , \end{aligned}$$

(4.2)

where $\left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} x(\xi , \zeta ) \\ y(\xi , \zeta ) \end{array} \right)$ is a NURBS representation of a mapping from an unit square to the corresponding plane domain, while $z(x,y)$ is directly given by (4.1).

Two cases, $\alpha = 1$ and $\alpha = 2/5$, will be considered here. Figure 7 shows how ${\mathfrak {D}}$ is divided into four patches according to its characteristics for the two cases.

Fig. 7

Partition of the domain ${\mathfrak {D}}$; Azimuth $= 64^{\circ }$, Elevation $= 14^{\circ }$. Left $\alpha = 1$; Right $\alpha = 2/5$

Table 4 Hemispherical surface, $\alpha = 1.\,L^{2}$ error and convergence order

We present the convergence results in the ${\mathcal {L}}^{2}$ norm for case $\alpha = 1$ and case $\alpha = 2/5$ in Tables 4 and 5, respectively. When $\alpha = 1$, it is a case where the geometrical mappings are smooth enough to guarantee the optimal ${\mathcal {L}}^{2}$ convergence rate, approximately $p+1$ with $p$ the lowest degree in each case, as can be observed from Table 4. This conforms to our theory presented in Sect. 3.

Since we adopt (4.2) as the geometrical mappings and $z(x,y)$ in it is directly given by (4.1), the geometrical mappings for some surface patches are no longer smooth enough to satisfy all the inequalities in Lemma 3.1 when $\alpha = 2/5$. Thus, convergence results in the ${\mathcal {L}}^{2}$ norm like previous examples can not be achieved in this case. We, in turn, compute the error $||\nabla _{s}(\widetilde{u}- \widetilde{u_{h}})||_{{\mathcal {L}}^{2}({\mathfrak {T}})}$. Here, we just consider the case where ${\mathcal {P}} = (2,2,2,2)$, and point out that errors do not decrease apparently when using higher degree NURBS functions. Values of $||\nabla _{s}(\widetilde{u}- \widetilde{u_{h}})||_{L^{2}({\mathfrak {T}})}$ with different $N$ and the corresponding convergence rates are presented in Table 5. We can obtain from Table 5 that the convergence rate is approximately 0.4, which is as good a result as shown for the same example in [36], in which an adaptive finite element method (AFEM) is presented for several applications governed by geometric PDEs.

Table 5 $\alpha = 2/5$. $||\nabla _{s}(\widetilde{u}- \widetilde{u_{h}})||_{L^{2}({\mathfrak {T}})}$ and convergence order

4.3 Linear Elasticity: Infinite Plate with Circular Hole Under Constant In-plane Tension in the $x$-direction

Similar methods can be easily applied to numerically approximate the solutions of more complex problems. Here we show an example for the system of linear elasticity. It is a benchmark problem in IGA literatures. The infinite plate is modeled by a finite quarter plate, which is partitioned into two patches in our implementation, see Fig. 8 (left).

Fig. 8

Linear elasticity. Partition of the domain ${\mathfrak {D}}$ (left) and the exact solution (right)

Conforming meshes are used in the implementation. The exact solution, given in [37], is applied as a Neumann boundary condition. We put it here for reference:

$$\begin{aligned} \sigma _{rr}(r,\theta )= & {} \frac{T_{x}}{2} \left( 1-\frac{R^{2}}{r^{2}}\right) + \frac{T_{x}}{2} \left( 1-4\frac{R^{2}}{r^{2}} +3\frac{R^{4}}{r^{4}} \right) \cos {(2\theta )},\nonumber \\ \sigma _{\theta \theta }(r,\theta )= & {} \frac{T_{x}}{2} \left( 1+\frac{R^{2}}{r^{2}}\right) - \frac{T_{x}}{2} \left( 1 +3\frac{R^{4}}{r^{4}} \right) \cos {(2\theta )},\nonumber \\ \sigma _{r\theta }(r,\theta )= & {} -\frac{T_{x}}{2} \left( 1+2\frac{R^{2}}{r^{2}} -3\frac{R^{4}}{r^{4}} \right) \sin {(2\theta )}, \end{aligned}$$

(4.3)

where $T_{x}$ is the magnitude of the stress applied for the infinity plate case.

The image of the exact solution is presented in Fig. 8 (right). The setup is illustrated in Fig. 9 [2]. $R$ is the radius of the hole, $L$ is the length of the finite quarter plate, $E$ is Young’s modulus, and $\nu$ is Poisson’s ration.

Fig. 9

Setup of an infinite elastic plate with a circular hole

Contour plots of the numerical solutions with ${\mathcal {P}} = (3,3)$ are shown in Fig. 10. The applied stress is $T_{x} = 10$ and the contours show that the stress concentration of $\sigma _{xx} = 30$ at the edge of the hole (i.e., at $r = R$, $\theta = 3\pi /2$) is obtained as the mesh is refined.

Table 6 gives the convergence results in the ${\mathcal {L}}^{2}$ norm of stress. As can be seen, the $L^{2}$ convergence rate of stress is approximately $p$ with $p$ the lowest degree in each case. This is the best according to our theory, since $\sigma _{xx}$ is a linear combination of the first-order derivatives of the displacement vector components, with which we are actually dealing in the implementation.

Fig. 10

Linear elasticity. Contour plots of the Numerical solutions with ${\mathcal {P}} = (3,3)$

Table 6 Elasticity. $L^{2}$ error and convergence order

5 Conclusions

In this paper, we proposed a method that combines IGA with the DG method together for solving second-order elliptic problem on (3D) surfaces. DG ideology is adopted at patch level, i.e., we employ the traditional IGA within each patch, and employ the DG method across the patch interfaces to glue the multiple patches. Our method enjoys the following advantages. First and most significantly, the time-consuming steps in mesh generation process in traditional FEA is no longer necessary and refinements, including $h$-refinement and $p$-refinement, can be easily performed at any level for the mesh and the corresponding basis functions. Moreover, the flexibility of DG method makes this method suitable for handling cases with non-conforming patches and different degrees across patch interfaces. In addition, the geometric flexibility of IGA basis functions, especially the use of multi-patch, enables a more accurate representation of a much larger class of complicated objects than standard finite element technology. In particular, all conic sections can be represented exactly, thus eliminating the geometrical errors at the beginning. Finally, this method can be easily formulated and implemented.

A new ideology is adopted when we perform the error analysis, in which we directly make use of the approximation property of the NURBS basis functions on the parametric domain rather than that of the IGA functions on the physical domain. Obviously, the former is easier to get. Theoretical analysis demonstrates that our method achieves the optimal convergence rate with respect to both the DG norm $||\cdot ||_{DG, {\mathfrak {T}}}$ and the ${\mathcal {L}}^{2}$ norm $||\cdot ||_{{\mathcal {L}}^{2}({\mathfrak {T}})}$. Numerical examples are presented to verify the theoretical result and gauge the good performance of our method.

For the future research, one possible modification is to make our method self-adaptive. There are two possible way to reach self-adaptivity. One is to make use of the flexibility of the DG method in adaptivity, we just need to cut patches into smaller ones whenever necessary by adaptivity. The other one is to adopt more powerful IGA basis function to achieve local refinement within each patch. Additionally, the thought proposed in this paper could obviously be used for numerically approximating the solutions of more complex problems with appropriate DG methods.