Modified PHT-splines

TL;DR

This paper introduces modified PHT-splines using split-in-half mesh refinement to enhance flexibility and efficiency in surface fitting and isogeometric analysis, especially for anisotropic features.

Contribution

The paper proposes modified hierarchical T-meshes and PHT-splines with split-in-half refinement to improve local refinement and flexibility over traditional methods.

Findings

Enhanced surface fitting accuracy with modified PHT-splines.

Improved isogeometric analysis performance on anisotropic features.

Numerical examples demonstrate the advantages of the new splines.

Abstract

The local refinement of PHT-splines (polynomial splines over hierarchical T-meshes) is achieved by a simple cross insertion, which may introduce superfluous control points or coefficients. By allowing split-in-half in mesh refinement, modified hierarchical T-meshes are defined. Using this approach, polynomial splines defined over the modified hierarchical T-meshes (modified PHT-splines) are introduced to increase the flexibility of PHT-splines. Numerical examples demonstrate the advantages of our new splines when applied to surface fitting and isogeometric analysis problems with anisotropic features.

Tables

Table 1. Table 1: Comparison of results

	original method			new method
Level 0	#cp	max error	average error	#cp	max error	average error
Level 1	4356	0.0174641	0.00161962	520	0.00486232	0.00160813
Level 2	4356	0.0174641	0.00162178	2244	0.0174388	0.00219253
Level 3	4356	0.0222689	0.00222854	3420	0.0209686	0.00248527
Level 4	6952	0.0546684	0.01776762	5552	0.0336198	0.0193418

Equations

$$\mathcal{S}(m,n,\alpha ,\beta ,\mathbb{T}):=\{s(x,y)\in C^{\alpha ,\beta }(\Omega )\text{arrowvert}s(x,y){|}_{\varphi }\in {\mathbb{P}}_{mn}\text{mbox}forany\varphi \in \mathcal{F}\}$$

$$\dim \mathcal{S}(3,3,1,1,\mathbb{T})=4(V^b+V^{+}),$$

$$\Xi =[s_0,s_0,s_1,s_1,s_2,s_2,\dots ,s_{m-1},s_{m-1},s_m,s_m],$$

$$[s_{i-1},s_{i-1},s_i,s_i,s_{i+1}]\times [t_{i-1},t_{i-1},t_i,t_i,t_{i+1}],$$

$$[s_{i-1},s_i,s_i,s_{i+1},s_{i+1}]\times [t_{i-1},t_{i-1},t_i,t_i,t_{i+1}],$$

$$[s_{i-1},s_{i-1},s_i,s_i,s_{i+1}]\times [t_{i-1},t_i,t_i,t_{i+1},t_{i+1}],$$

$$[s_{i-1},s_i,s_i,s_{i+1},s_{i+1}]\times [t_{i-1},t_i,t_i,t_{i+1},t_{i+1}]$$

$$\bar{b}_{i}^{k+1}(s,t)=\mathfrak{M}^{k+1}\big{(}b_{i}^{k}(s,t)\big{)}$$

$${\mathcal{S}}_A^{k+1}=\{{\mathfrak{M}}^{k+1}(b_i^k(s,t)):b_i^k(s,t)\in {\mathcal{S}}^k\},$$

$${\mathcal{S}}_B^{k+1}=\{b_i^{k+1}(s,t):\text{mbox}newbasisfunctionsatlevelk+1\}.$$

$$f(s,t)=\sum _{j=0}^{4n-1}{c}_j{b}_j(s,t).$$

$$f(s,t)\equiv 0\Rightarrow c_j=0,j=0,\dots ,4n-1.$$

$$\mathfrak{L}(f(s,t))=\big{(}f(s,t),f_{s}(s,t),f_{t}(s,t),f_{st}(s,t)\big{)}.$$

$$\mathfrak{L}(f({\mathbf{v}}_i))=\sum _{j=0}^3{c}_{4i+j}\mathfrak{L}(b_{4i+j}({\mathbf{v}}_i))=(c_{4i},c_{4i+1},c_{4i+2},c_{4i+3})\cdot \mathbf{B}({\mathbf{v}}_i)=0,$$

$$\mathbf{B}(s,t)=$$

$$\mathbf{B}(\mathbf{v}_{i})=\left(\begin{array}[]{cccc}(1-\lambda)(1-\mu)&\lambda(1-\mu)&\lambda\mu&\mu(1-\lambda)\\ \lambda(1-\mu)&\alpha(1-\mu)&-\beta\mu&\alpha\beta\\ (1-\lambda)\mu&\alpha\mu&\beta(1-\lambda)&\alpha\beta\\ \lambda\mu&\alpha\mu&\beta\lambda&\alpha\beta\\ \end{array}\right),$$

$$f_1(s,t)=\sum _{j=0}^{4V_{k+1}-1}{b}_j^{k+1}(s,t),f_2(s,t)=\sum _{j=0}^{d_k-1}{\bar{b}}_j^{k+1}(s,t),$$

$$\mathbf{S}(s,t)=\sum _{j=0}^d{\mathbf{C}}_j{b}_j(s,t),(s,t)\in [0,1]\times [0,1],$$

$$\min \sum _{i=0}^{M-1}\Vert \mathbf{S}(s_i,t_i)-{\mathbf{P}}_i{\Vert }^2=\underset{{\mathbf{C}}_j}{\min }\sum _{i=0}^{M-1}\Vert \sum _{j=0}^d{\mathbf{C}}_j{b}_j(s_i,t_i)-{\mathbf{P}}_i{\Vert }^2,$$

$$\mathfrak{L}(\mathbf{S}({\mathbf{v}}_i))=\sum _{j=0}^d{\mathbf{C}}_j\mathfrak{L}{b}_j({\mathbf{v}}_i)=\sum _{k=0}^3{\mathbf{C}}_{4i+k}\mathfrak{L}{b}_{4i+k}({\mathbf{v}}_i)=\mathbf{C}\cdot \mathbf{B}({\mathbf{v}}_i),$$

$${\kappa_{s}}_{j}=\frac{\|\frac{\partial\mathbf{S}_{k}}{\partial s}\times\frac{\partial^{2}\mathbf{S}_{k}}{\partial s^{2}}\|}{\|\frac{\partial\mathbf{S}_{k}}{\partial s}\|^{3}}\bigg{|}_{(\overline{s}_{j},\overline{t}_{j})},\ \ \ {\kappa_{t}}_{j}=\frac{\|\frac{\partial\mathbf{S}_{k}}{\partial t}\times\frac{\partial^{2}\mathbf{S}_{k}}{\partial t^{2}}\|}{\|\frac{\partial\mathbf{S}_{k}}{\partial t}\|^{3}}\bigg{|}_{(\overline{s}_{j},\overline{t}_{j})},\ j=1,\ldots,l.$$

$$f(u,v)$$

$$+B_{1,7}(u)(2\sin (120u)\sin (2\pi u))$$

$$+B_{7,7}(u)(2-2(1+0.4\sin (60v))|\cos (2\pi v)|)),$$

$$-\Delta u=f\text{mbox}in\Omega ,$$

$$u=0\text{mbox}on{\Gamma }_D,$$

$$\frac{\partial u}{\partial \mathbf{n}}=g_N\text{mbox}on{\Gamma }_N$$

$$\mathbf{G}:$$

$$(s,t)\mapsto (x,y)=\sum _{i=1}^d{\mathbf{P}}_i{b}_i(s,t),$$

$$H^1(\Omega )=\{v\in L^2(\Omega ):|\nabla v|\in L^2(\Omega )\}$$

$$a(u,v)=⟨f,v⟩,\forall v\in V$$

$$a(u,v)={\int }_{\Omega }\nabla u\cdot \nabla vd\Omega ,$$

$$⟨f,v⟩={\int }_{\Omega }fvd\Omega +{\int }_{{\Gamma }_N}{g}_{N}vd\Gamma$$

$$a(u^h,v^h)=⟨f,v^h⟩,\forall v^h\in V^h,$$

$$V^h=\text{mbox}span\{{\psi }_i(x,y)=b_i\circ {\mathbf{G}}^{-1},{\psi }_i(x,y){|}_{{\Gamma }_D}=0,i=1,2,\dots ,n\}.$$

Full text

1 Introduction

The multivariate splines used in geometric design (GM) and isogeometric analysis (IGA) are often based on the tensor-product of B-splines or NURBS. However, the standard tensor-product mechanism prevents a localized editing of the mesh. As such, a large number of superfluous control points or coefficients may be contained under the tensor-product representations when we deal with problems having anisotropic features. To overcome this weakness and provide more flexibility, locally refinable splines are introduced, which generalize the B-spline model in different ways. Hence the development of locally refinable splines has recently become an active research topic in geometric modeling and isogeometric analysis. There are a number of different locally refinable splines, and among them, a) T-splines, b) HB-splines and THB-splines, c) LR-splines, d) PHT-splines are very popular. We shall review these splines briefly below.

1. By allowing T-junctions in the control meshes, T-splines [31, 30] were introduced to remove the redundant control points in geometric modeling. The potential of T-splines in isogeometric analysis was reported in [1, 8]. Additionally, AST-splines (analysis-suitable T-splines) were introduced to fix the linear dependency problem of T-splines [24, 29], and are important for their applications in isogeometric analysis. Recently, AS++ T-splines (analysis-suitable++ T-splines), which include analysis-suitable (AS) T-splines as a special case, were also introduced [23].

2. As a classical extension of tensor-product representations, hierarchical B-splines (HB-splines) allow an effective local control of refinement [13]. The basis functions of HB-splines are locally supported, linearly independent and non-negative. The application of HB-splines in geometric modeling and isogeometric analysis was reported in [11, 17, 28, 37]. Unfortunately, the standard hierarchical construction does not preserve the partition of unity property. By adequately truncating basis functions in the construction of HB-splines, truncated hierarchical B-splines (THB-splines), which form a convex partition of unity, were introduced by Giannelli et al. in [14, 16]. In addition, based on their good stability [16] and approximation abilities [33], THB-splines are well suited for applications in geometric modeling, computer aided design, and isogeometric analysis [15, 16, 20, 35]. In addition, based on the techniques of truncation, the normalized quasi-hierarchical basis was proposed for hierarchical Powell-Sabin splines in [34]. Recently, truncated T-splines are also introduced in [40].

3. Based on the idea of splitting basis functions, locally refined (LR) splines were proposed by Dokken et al. in [7]. LR splines possess the nested spaces property while they have difficulties with linear independence. The properties of LR splines and their hierarchical construction in 2D were further discussed in [2, 3]. LR splines were also applied in isogeometric analysis [19]. Recently, Johannessen et al. analyzed the corresponding stiffness and mass matrices for HB-splines, THB-splines and LR-splines, in terms of sparsity patterns and conditioning numbers [18]. The comparison shows that the basis functions of these splines in general do not span the same space, and that conditioning numbers are comparable.

4. Polynomial splines over hierarchical T-meshes (PHT-splines), i.e., bi-cubic splines over hierarchical T-meshes with reduced $C^{1}$ regularity, use a modification mechanism to obtain nested spline spaces [5]. PHT-splines modify their basis functions in coarse levels directly by resetting the Bézier ordinates corresponding to the new basis vertices. The modification mechanism ensures the properties like partition of unity, nonnegativity, and linear independence. Therefore, PHT-splines can be applied successfully in surface model reconstruction [5, 22, 38], adaptive finite element method for elliptic equations [36], and isogeometric analysis for solving elastic problems [26, 27, 39]. The local refinement for hierarchical T-meshes of PHT-splines is achieved by cross insertion [5], i.e., splitting each candidate cell into four subcells by inserting a cross. This simple refinement rule may still introduce redundant control points or coefficients when applied to the models or problems with anisotropic features, i.e., the problems or models are directionally dependent.

The present paper is devoted to increasing flexibility to PHT-splines, which can handle models or problems with anisotropic features efficiently. Instead of performing refinement exclusively by cross insertion, modified hierarchical T-meshes are obtained by allowing for the splitting of cells into halves (horizontal or vertical). Similar to the construction of classical PHT-splines, polynomial splines over modified hierarchical T-meshes (modified PHT-splines) are also presented. Numerical experiments show that our new splines have advantages when applied to problems with anisotropic features. It is worth noting that using anisotropic refinement for surface construction has been reported recently in [9, 10].

The remainder of the paper is organized as follows. Section 2 reviews the definition of T-meshes and hierarchical T-meshes. The definition of modified hierarchical T-meshes and the dimension formula for polynomial spline spaces over modified hierarchical T-meshes are also provided. In Section 3, we discuss the construction of basis functions for modified PHT-splines. Section 4 presents three examples for fitting open meshes with modified PHT-splines. In Section 5, the adaptive isogeometric analysis based on modified PHT-splines is provided. Section 6 concludes the paper with a summary and some future work.

2 Polynomial splines over T-meshes

In this section, we briefly review some notations of T-meshes and hierarchical T-meshes. The definition of modified hierarchical T-meshes is provided and the dimension formula for polynomial splines space over modified T-meshes is also given.

2.1 T-meshes

A T-mesh is a rectangular grid that allows T-junctions [6, 31]. It is assumed that except for the corner case, the endpoints of each grid line in the T-mesh must be on two other grid lines, and each cell or facet in the grid must be a rectangle. Figure 1(a) shows an example of a T-mesh, and Figure 1(b) shows an example of a non-T-mesh. A grid point in a T-mesh is also called a vertex of the T-mesh. If a vertex is on the boundary of the domain, then it is called a boundary vertex. Otherwise, it is called an interior vertex. For example, $\mathbf{b}_{i}$, $i=1,\ldots,12$ in Figure 1(a) are boundary vertices, while $\mathbf{v}_{i}$, $i=1,\ldots,9$ are interior vertices. Interior vertices have two types: one is crossing vertex, e.g., $\mathbf{v}_{1}$, $\mathbf{v}_{3}$, $\mathbf{v}_{5}$, $\mathbf{v}_{7}$, and $\mathbf{v}_{9}$ in Figure 1(a), and the other is T-vertex, e.g., $\mathbf{v}_{2}$, $\mathbf{v}_{4}$, $\mathbf{v}_{6}$, and $\mathbf{v}_{8}$ in Figure 1(a). The line segment connecting two adjacent vertices on a grid line is called an edge of the T-mesh. If an edge is on the boundary of the T-mesh, it is called a boundary edge. Otherwise, it is called an interior edge. A cell is called an interior cell if all four edges are interior. Otherwise, it is called a boundary cell.

In terms of neighbor relations of cells in a T-mesh, a pair of cells are called adjacent if they share one common edge. Moreover, a pair of adjacent cells are called horizontally (vertically) aligned-adjacent if they share only one vertical (horizontal) edge, and the knot intervals of these two cells are identical in the vertical (horizontal) direction. A pair of cells are aligned-adjacent if they are horizontally or vertically aligned-adjacent.

2.2 Hierarchical T-meshes and modified hierarchical T-meshes

As a special type of T-mesh, hierarchical T-mesh $\mathbb{T}$ is defined in a hierarchical manner [5]. To describe the level of hierarchical T-mesh, we fix the following notations in the rest of the paper. $\mathbb{T}_{k}$ denotes the T-mesh at level $k$. Let $\mathcal{F}_{k}$ be the set consisting of all cells of level $k$. Note that the level of cells and the level of mesh are defined separately; not every cell in $\mathbb{T}_{k}$ is of level $k$. In this paper, we focus on a special type of hierarchical T-mesh, i.e., the cells to be subdivided are chosen from $\mathcal{F}_{k}$. Let $\Theta_{k}\subset\mathcal{F}_{k}$ be the set consisting of all cells in $\mathcal{F}_{k}$ that are to be subdivided at level $k$. Then, a hierarchical T-mesh $\mathbb{T}$ is defined as follows:

Start with a tensor product mesh $\mathbb{T}_{0}$. Let $\mathcal{F}_{0}$ be the set of all cells in $\mathbb{T}_{0}$. 2. 2.

Recursive steps: Subdivide each cell in $\Theta_{k}$ by inserting a cross to get a new T-mesh $\mathbb{T}_{k+1},k=0,1,\ldots,N-1$. Let $\mathcal{F}_{k+1}$ be the set consisting of all new cells, which emerge at level $k+1$ after subdivision. 3. 3.

$\mathbb{T}=\mathbb{T}_{N}$.

Remark 1

Based on the above (ii), if a cell $\theta\notin\Theta_{k}$, i.e., $\theta$ is not subdivided at level $k$, the cell $\theta$ will be excluded from subdivision henceforth.

The level of a vertex is defined as follows. Given a vertex $\mathbf{v}$ in $\mathbb{T}$, if $\mathbf{v}\in\mathbb{T}_{k}$ but $\mathbf{v}\not\in\mathbb{T}_{l},l=0,\ldots,k-1$, then it is called a vertex of level $k$. Specifically, every vertex in $\mathbb{T}_{0}$ is a vertex of level 0.

In Figure 2, an example is provided to illustrate the dynamic refinement process of a hierarchical T-mesh.

To improve the flexibility of hierarchical T-meshes, we modify hierarchical T-meshes as follows.

Definition 1

A modified hierarchical T-mesh is generated as follows. In the recursive step (ii) of hierarchical T-meshes, besides cross insertion, cells in $\Theta_{k}$ can be subdivided by inserting a single (horizontal or vertical) edge in half.

Figure 3 is an example of the refinement process of a modified hierarchical T-mesh.

It should be noted that a strategy of type selection for cells will be provided in Section 3.1.

2.3 Spline spaces over modified hierarchical T-meshes

Given a T-mesh $\mathbb{T}$, $\mathcal{F}$ represents the set of all the cells in $\mathbb{T}$ and $\Omega$ represents the region occupied by $\mathcal{F}$. The polynomial spline space over $\mathbb{T}$ is defined as

[TABLE]

where $\mathbb{P}_{mn}$ is the space of all the polynomials with bi-degree $(m,n)$, and $C^{\alpha,\beta}(\Omega)$ is the space consisting of all the bivariate functions that are continuous in $\Omega$ with order $\alpha$ in the $x$-direction and with order $\beta$ in the $y$-direction.

An explicit formula of the spline space $\mathcal{S}(m,n,\alpha,\beta,\mathbb{T})$ in the case of $m\geq 2\alpha+1$, $n\geq 2\beta+1$ is provided in [6]. For the spline space $\mathcal{S}(3,3,1,1,\mathbb{T})$ with $\mathbb{T}$ being a modified hierarchical T-mesh, the dimension formula can be simplified into

[TABLE]

where $V^{b}$ and $V^{+}$ denote the number of boundary vertices and interior crossing vertices in $\mathbb{T}$ respectively. The dimension formula (1) implies that every boundary vertex or interior crossing vertex corresponds to four basis functions. Therefore, following the method in [5], we call a boundary vertex or an interior crossing vertex a basis vertex.

For brevity, the polynomial spline in $\mathcal{S}(3,3,1,1,\mathbb{T})$ defined over a modified hierarchical T-mesh is called modified PHT-spline.

3 Basis functions of polynomial splines over modified hierarchical T-meshes

Inspired by [5], basis functions of the modified PHT-splines can be constructed in a level-by-level approach.

For the initial level $\mathbb{T}_{0}$, the standard bicubic $C^{1}$ continuous tensor-product B-splines are used as basis functions. Let the knot vector of a $C^{1}$ continuous cubic spline be

[TABLE]

where $s_{i}<s_{i+1},1\leq i\leq m-2$, $s_{0}=s_{1}$, and $s_{m-1}=s_{m}$. For every interior knot $s_{i}$, there are two basis functions with two knot vectors $[s_{i-1},s_{i-1},s_{i},s_{i},s_{i+1}]$, $[s_{i-1},s_{i},s_{i},s_{i+1},s_{i+1}]$. Hence these two basis functions have support $[s_{i-1},s_{i+1}]$. Extending this fact to the tensor-product case of the surface, every vertex of $\mathbb{T}_{0}$ is associated with four basis functions with knot vectors

[TABLE]

respectively.

3.1 Mesh refinement at level $k$

Assume $\Theta_{k}$ is given at level $k$ first. For each cell in $\Theta_{k}$, its anisotropic information is also provided as labels, which indicates the subdivision type of the cell that tends to be selected. Namely, for each cell $\theta\in\Theta_{k}$, it has been labeled by ‘H’, ‘V’, or ‘C’ based on anisotropic estimation. Here, ‘H’, ‘V’, and ‘C’ represent horizontal subdivision, vertical subdivision, and cross insertion respectively, which are three different subdivision types that $\theta$ tends to choose. For brevity, the label of the cell $\theta$ is denoted by label$(\theta)$ in the rest of the paper. For more details about the estimation of anisotropic information, see Section 4.3 and Section 5.3. If we refine these cells directly as their labels indicate, then the following two issues may occur, which increase the complexity in the construction of the basis functions.

Some T-vertices in $\mathbb{T}_{l}$ $(l<k)$ change into crossing vertices in $\mathbb{T}_{k}$ (see Figure 4 for an example). 2. 2.

No basis vertex appears for some cells belonging to $\Theta_{k}$ after subdivision (see Figure 5 for an example).

With the first issue, the basis construction is slightly complicated, and we need to revisit the levels at which these T-vertices first appear, which is required in the modification of basis functions at a later stage. With the second issue, invalid refinement will occur after subdivision, which contradicts the most common situation in geometric modeling and adaptive isogeometric analysis applications.

Hence, in the current paper, we explore the construction of polynomial spline in $\mathcal{S}(3,3,1,1,\mathbb{T})$, which is defined over modified hierarchical T-meshes in the absence of the above two issues. To resolve the above two issues, we provide a refinement strategy to guide subdivision in the following Algorithm 1. The refinement strategy is based on anisotropic features and neighbor relations of cells.

Algorithm 1. Anisotropic Refinement Strategy

Input: $\mathbb{T}_{k}$, $\Theta_{k}$, $\{\mbox{label}(\theta)\mid\theta\in\Theta_{k}\}$

Output: $\mathbb{T}_{k+1}$

Classify $\Theta_{k}$ into several connected groups $G_{1}^{k},\cdots,G_{M}^{k}$ (See Remark 2 for more details). 2. 2.

For each connected group $G_{i}^{k}$, check label$(\theta)$ for all $\theta\in G_{i}^{k}$.

1. (a)

When labels are identical, check whether new basis vertices will appear for each $\theta\in G_{i}^{k}$ or not.

If no, for each cell $\theta_{0}\in G_{i}^{k}$ without new basis vertices, find all cells that are in $G_{i}^{k}$ and connected to $\theta_{0}$ by a sequence of vertically (or horizontally, if there exists) aligned-adjacent cells in $G_{i}^{k}$. Relabel them with ‘C’. 2. (b)

When labels are nonidentical, check neighbor cells for each $\theta\in G_{i}^{k}$.

1. i.

If $\theta$ only has horizontally aligned-adjacent cells $\overline{\theta}\in G_{i}^{k}$, $\mbox{label}(\theta)=$ ‘H’; 2. ii.

If $\theta$ only has vertically aligned-adjacent cells $\overline{\theta}\in G_{i}^{k}$, $\mbox{label}(\theta)=$ ‘V’; 3. iii.

if $\theta$ has both vertically and horizontally aligned-adjacent cells or has not any aligned-adjacent cell $\overline{\theta}\in G_{i}^{k}$, $\mbox{label}(\theta)=$ ‘C’. 3. 3.

Subdivide cells as new labels indicates. For each cell $\theta\in G_{i}^{k}$,

1. (a)

if $\theta$ has the label ‘H’, split $\theta$ in half horizontally; 2. (b)

if $\theta$ has the label ‘V’, split $\theta$ in half vertically; 3. (c)

if $\theta$ has the label ‘C’, subdivide $\theta$ by inserting a cross. 4. 4.

Output $\mathbb{T}_{k+1}$.

Remark 2

In Step 1) of Algorithm 1, the classification is achieved by ”flood-fill” through aligned-adjacent relations. Namely, two cells of level $k$ will be classified into one connected group, 1) if they are adjacent, they must be aligned-adjacent; 2) if they are not adjacent, there exists a sequence of aligned-adjacent cells in $\Theta_{k}$ connecting them. Namely, for two cells $\theta_{s},\theta_{e}\in\Theta_{k}$, if they are not adjacent, there exist $\hat{\theta}_{j}^{k}\in\Theta_{k},j=1,\ldots,L$ such that $\hat{\theta}_{1}^{k}=\theta_{s}$, $\hat{\theta}_{L}^{k}=\theta_{e}$, $\hat{\theta}_{j}^{k}$ and $\hat{\theta}_{j+1}^{k}$ are aligned-adjacent, $j=1,\ldots,L-1$.

Remark 3

In Step 2.1) of Algorithm 1, cells in $G_{i}^{k}$ are relabeled based on anisotropic information if their original labels are identical. Otherwise, cells in $G_{i}^{k}$ are relabeled in Step 2.2) based on neighbor relations.

An example in Figure 6 illustrates our subdivision strategy. All cells in an initial tensor-product mesh $\mathbb{T}_{0}$ have labels ‘H’(see Figure 6(a)). After executing Algorithm 1 once, we get $\mathbb{T}_{1}$ in Figure 6(b). Cells to be subdivided in $\mathbb{T}_{1}$ are classified into three connected groups (see Figure 6(c)). If we subdivide the cells directly as their labels indicate, no basis vertex will appear in the upper right cell in the yellow group, and T-vertices may change into basis vertices after subdividing cells several times (blue parts). By Step 2) and Step 3) of our strategy, we get $\mathbb{T}_{2}$ in Figure 6(d).

The connected group plays an important role in Algorithm 1. Here we shall provide the properties of the connected groups.

Lemma 2

Suppose $G$ is a connected group. Subdivide all cells in $G$ as in Algorithm 1, and let $G^{\prime}$ be the set consisting of all new cells. Then:

There is no T-vertex on the interior edges of the group $G$. 2. 2.

$G^{\prime}$* forms a new connected group.* 3. 3.

There is no T-vertex on the interior edges of the group $G^{\prime}$. 4. 4.

Let $G_{1}$ be a subset of $G^{\prime}$. Assume $G_{1}$ is classified into several connected groups (based on Remark 2). Then, there is no common edge between any two connected groups.

Proof** 1**

1. There is no T-vertex on the interior edges of $G$ before subdivision. Otherwise, three cells around the interior T-vertex, which are in one connected group, are adjacent to each other but not aligned-adjacent, which contradicts the assumption that they are from one connected group.

2. Consider any two adjacent cells, denoted by $\theta_{0}$ and $\theta_{1}$, in $G$ first. By Remark 2, $\theta_{0}$ and $\theta_{1}$ are also aligned-adjacent. Let the common edge of $\theta_{0}$ and $\theta_{1}$ be $\mathbf{e}$. Assume that $\theta_{0}$ and $\theta_{1}$ are horizontally (vertically) adjacent, then a T-vertex will appear on the common edge $\mathbf{e}$ if

one is subdivided by inserting a cross, and the other is split in half vertically (horizontally), or 2. 2.

one is split in half vertically, and the other is split in half horizontally.

However, by Step 2.1) and Step 2.2) in Algorithm 1, neither of the above cases will occur after subdivision. Hence, any new cells acquired by subdividing two adjacent cells will be aligned-adjacent if they are adjacent. If new cells are not adjacent, there exists a sequence of new aligned-adjacent cells connecting them. Consequently, all these new cells will be classified into the same connected group.

Repeat the above process for any two adjacent cells in $G$. We get that all new cells, which are generated by subdividing cells in $G$ as Algorithm 1, are in the same connected group.

3. This result is a direct consequence of property 1 and property 2 in Lemma 2.

4. Note that there is no common cell between any two different connected groups. Otherwise, these two connected groups will merge into a bigger connected group through the common cells. Hence, suppose there exists a common edge $\mathbf{e}$ between two different connected groups. Then, the common edge $\mathbf{e}$ is shared by two cells, which are from different connected groups. By property 2 in Lemma 2, $G^{\prime}$ is a connected group, and it follows that these two cells are aligned-adjacent. Therefore, these two connected groups should be classified into the same connected group, which contradicts our assumption that they are classified into different connected groups. \qed

Theorem 3

Start from a tensor-product mesh $\mathbb{T}_{0}$, let $\mathbb{T}_{N}$ be the mesh output by performing Algorithm 1 $N$ times ($N>0$). Then, $\mathbb{T}_{N}$ is a modified hierarchical T-mesh. Hence, for each level,

no T-vertex changes into a crossing vertex; 2. 2.

new basis vertices will appear for each cell that is to be subdivided.

Proof** 2**

In Algorithm 1, cells are allowed to be subdivided by splitting in half or inserting crosses. It is straightforward to check that $\mathbb{T}_{N}$ is a modified hierarchical T-mesh.

1. For any connect group $G_{i}^{k}$ at level $k$, by property 1 of Lemma 2, there is no T-vertex on the interior edges of the group $G_{i}^{k}$ before subdivision. Therefore, the change from a T-vertex to a crossing vertex on the interior edges of the group is impossible.

On the other hand, since we start with a tensor product mesh $\mathbb{T}_{0}$, it follows by Property 2 of Lemma 2 that there is no common edge for any two different connected groups at level 0. Note that by Remark 1, if a cell is not subdivided at level $k$, it will be excluded from subdivision henceforth. Hence, by property 4 of Lemma 2, for each time we perform Algorithm 1, any two different connected groups at the same level have no common edge. Consequently, no T-vertex on the boundary edges of connected groups changes into a crossing vertex at each level.

Therefore, no T-vertex will change into a crossing vertex if we subdivide cells as Algorithm 1.

2. At each level, when the labels of cells in a connected group $G_{i}^{k}$ are identical, Algorithm 1 will subdivide these cells directly as their labels indicate when new basis vertices will appear for each cell. The exceptions are

for those cells with label ‘H’, there are no horizontally neighbor cells in $G_{i}^{k}$; 2. 2.

for those cells with label ‘V’, there are no vertically neighbor cells in $G_{i}^{k}$.

With the case (a) (case (b)), Algorithm 1 will relabel these cells and their neighbor cells in the same row (column) with ‘C’ as Step 2.1), which ensures the new basis vertices for these cells.

When the labels of cells in a connected group are not identical, Algorithm 1 will subdivide these cells based on neighboring relations. It is straightforward to verify that for each cell in $G_{i}^{k}$, at least one basis vertex will appear on the edge of this cell. Thus, the appearance of new basis vertices for each cell, which is to be subdivided, is guaranteed in Algorithm 1.\qed

3.2 Basis construction

The concept in constructing basis functions relies on the modifying mechanism that operates on the Bézier form of basis functions directly.

Definition 4

Given a basis function $b_{i}^{k}(s,t)$ at level $k$, for each cell $\theta$ of $\Theta_{k}$ with \theta\subset\mbox{supp}\big{(}b_{i}^{k}(s,t)\big{)}. The modification of $b_{i}^{k}(s,t)$ from level $k$ to level $k+1$

[TABLE]

is defined as follows.

Represent $b_{i}^{k}(s,t)$ over $\theta$ by specifying 16 Bézier ordinates. 2. 2.

Generate new Bézier ordinates by applying de Casteljau’s algorithm according to the different subdivision type (see Figure 7 for an illustration). 3. 3.

Set all of the Bézier ordinates that are associated with the new basis vertices to zero.

Note that the function $b_{i}^{k}(s,t)$ remains unchanged if there is no cell $\theta$ in $\Theta_{k}$ such that \theta\subset\mbox{supp}\big{(}b_{i}^{k}(s,t)\big{)}. For each cell that is to be subdivided at level $k+1$, the 16 Bézier ordinates are divided into four or two parts according to its subdivision type. Hence, each part is associated with a cell corner vertex, which is illustrated in Figure 8. By setting all of the Bézier ordinates that are associated with the new basis vertices to zero (see Figure 9 for an illustration), the basis function $b_{i}^{k}(s,t)$ at level $k$ is modified into a basis function $\bar{b}_{i}^{k}(s,t)$ at level $k+1$. In addition, the Bézier ordinates around the basis vertex (located in neighbor cells) are reset simultaneously. Hence, the conditions of $C^{1}$ continuity still hold and $\bar{b}_{i}^{k+1}(s,t)\in\mathcal{S}(3,3,1,1,\mathbb{T}_{k+1})$.

Next, we shall discuss adding new basis functions at level $k+1$.

By the definition of modified hierarchical T-mesh, a new interior basis vertex appears at level $k+1$ if

a cell at level $k$ is subdivided by cross insertion, which introduces a crossing vertex at the center of the cell; 2. 2.

the common edge of two aligned-adjacent cells is split after subdivision, which introduces a crossing vertex on the common edge.

For the case (a), the neighboring four cells around the new interior vertex are just four subcells of the original cell. For the case (b), the neighboring four cells around the new interior vertex are also aligned-adjacent. Therefore, for each new basis vertex, if it is an interior vertex, its neighboring four cells form a $2\times 2$ tensor-product mesh. On the other hand, new boundary basis vertices are introduced by subdividing of the boundary cells at level $k$. Therefore, for each new boundary basis vertex, its neighboring two cells are two aligned-adjacent subcells of the original boundary cell at level $k$. Hence, these two subcells also form a tensor-product mesh. Thus, for each new basis vertex $\mathbf{v}_{i}^{k+1}$, four new basis functions associated with the basis vertex $\mathbf{v}_{i}^{k+1}$ can be constructed.

We are now ready to construct the basis functions of modified PHT-splines.

Definition 5

Given a modified hierarchical T-mesh $\mathbb{T}$ in the 2D plane. A set of basis functions $b_{i}(s,t)$ for the spline space $\mathcal{S}(3,3,1,1,\mathbb{T})$ is recursively constructed as follows:

Initialization: $\mathcal{S}^{0}=\{b_{i}^{0}(s,t):$ tensor-product bicubic B-spline basis function over $\mathbb{T}_{0}$}. 2. 2.

Recursive case: $S^{k+1}=\mathcal{S}_{A}^{k+1}\cup\mathcal{S}_{B}^{k+1}$ for $k=0,\ldots,N-1$, where

[TABLE] 3. 3.

$\mathcal{S}=\mathcal{S}^{N}$**

Next we shall discuss the properties of the basis functions constructed according to the above definition.

Theorem 6

Let $\mathbb{T}$ be a modified hierarchical T-mesh. The modified PHT-spline basis $\mathcal{S}$ is constructed according to Definition 5. Then:

The functions in $\mathcal{S}$ are linearly independent; 2. 2.

$\mathcal{S}$* spans the modified PHT-spline space $\mathcal{S}(3,3,1,1,\mathbb{T})$.* 3. 3.

$\mathcal{S}$* forms a partition of unity.*

Proof** 3**

Assume that there are $n$ basis vertices in $\mathbb{T}$, denoted by $\mathbf{v}_{i},i=0,1,\ldots,n-1$. Let $b_{4i+j}(s,t),j=0,1,2,3$ be the four basis functions that associated with $\mathbf{v}_{i}$. By the construction of $\mathcal{S}$, it follows that there are $4n$ basis functions in $\mathcal{S}$.

a) Let $f(s,t)$ be a linear combination of all functions in $\mathcal{S}$,

[TABLE]

The linearly independence of $\mathcal{S}$ can be proven if

[TABLE]

Define a linear operator $\mathfrak{L}(\cdot)$ as

[TABLE]

Substituting $\mathbf{v}_{i}$ into (2) yields

[TABLE]

where $c_{4i+j}$ is the coefficient for $b_{4i+j}(s,t)$ and

[TABLE]

Suppose that the basis vertex $\mathbf{v}_{i}=(s_{i},t_{i})$ is of level $k,k=0,1,\ldots,N$. Then its four neighbor cells at level $k$ form a rectangle with four vertices $(s_{i}-3\Delta s_{i},t_{i}-3\Delta t_{i})$, $(s_{i}+3\Delta s_{i+1},t_{i}-3\Delta t_{i})$, $(s_{i}-3\Delta s_{i},t_{i}+3\Delta t_{i+1})$, $(s_{i}+3\Delta s_{i+1},t_{i}+3\Delta t_{i+1})$ (from left to right and bottom to top). Hence $\mathbf{B}(\mathbf{v}_{i})$ can be expressed by

[TABLE]

where $\alpha=\frac{1}{\Delta s_{i}+\Delta s_{i+1}},\beta=\frac{1}{\Delta t_{i}+\Delta t_{i+1}}$, $\lambda=\alpha\Delta s_{i}$, $\mu=\beta\Delta t_{i}$. It is straightforward to check that $\mathbf{B}(\mathbf{v}_{i})$ is invertible. Hence $(c_{4i},c_{4i+1},c_{4i+2},c_{4i+3})=\mathbf{0}$.

Iterating through each basis vertex in $\mathbb{T}$, and finally we have $c_{i}=0,i=0,\ldots,4n-1$. Hence the functions in $\mathcal{S}$ are linearly independent.

b) Since $b_{i}(s,t)\in\mathcal{S}(3,3,1,1,\mathbb{T})$, $i=0,\ldots,4n-1$ and dim $\mathcal{S}(3,3,1,1,\mathbb{T})=4n$, it follows directly by a) that $\mathcal{S}$ spans the modified PHT-spline space $\mathcal{S}(3,3,1,1,\mathbb{T})$.

c) From level $k$ to level $k+1$, we can define

[TABLE]

where $b_{j}^{k+1}(s,t),j=0,\ldots 4V_{k+1}-1$ are new basis functions associated with new basis vertices at level $k+1$, and $\bar{b}_{j}^{k+1}(s,t),j=0,\ldots d_{k}-1$ are basis functions after modification according to Definition 4 above. Check the Bézier ordinates of $f_{1}(s,t)$ and $f_{1}(s,t)$ in every cell. For $f_{1}(s,t)$, the Bézier ordinates associated with new basis vertices are one, while all the others are zero; for $f_{2}(s,t)$, the Bézier ordinates associated with new basis vertices are zero, while all of the others are one. Thus $f_{1}(s,t)+f_{2}(s,t)\equiv 1.$ \qed

In addition, it is straightforward to check that the basis functions in $\mathcal{S}$ have other properties, such as nonnegativity, and local support.

Remark 4

*By setting certain Bézier ordinates to zero (as done in Definition 4), we modify the basis functions, while THB-splines use similar techniques, which truncates the basis functions for splines with a more general setting [14]. *

Definition 7

Let $\mathbb{T}$ be a modified hierarchical T-mesh, and $b_{j}(s,t)$, $j=0,1,\ldots,d$ be the basis functions constructed as Definition 5. A modified PHT-spline surface over $\mathbb{T}$ is defined by

[TABLE]

where $\mathbf{C}_{j},j=0,1,\ldots,d$, are control points.

Similar to PHT-spline surfaces in [5], modified PHT-spline surfaces have beneficial properties such as convex hull, affine invariant, and local support.

4 Fitting open meshes

Surface fitting is a fundamental task in computer aided geometric design and computer graphics that has been discussed in many papers (see [4] for a review of the literature). To show the potential of modified PHT-splines, we present a scheme to fit open mesh models based on modified PHT-spline surfaces in this section.

Given an open mesh $\mathcal{M}$ with vertices $\mathbf{P}_{j}$, $j=0,1,\ldots,M-1$. Let the corresponding parameter values of $\mathbf{P}_{j}$ be $(s_{j},t_{j})$, $j=0,1,\ldots,M-1$, which are calculated from some parameterization of the mesh (uniform weights are taken in barycentric mapping for the current paper). The parameter domain is $[0,1]\times[0,1]$. For more details about mesh parameterization, we refer readers to a specific review on this topic in [12].

4.1 Outline of the fitting process

The surface fitting scheme repeats the sequential steps 2, 3 and 4 until the fitting error in each cell is less than some tolerance $\varepsilon$.

Initialize a tensor product mesh $\mathbb{T}_{0}$. Set level $k=0$. 2. 2.

Keep unchanged the control points associated with the old basis functions. Compute the control points for the new basis functions on the $k$th level mesh $\mathbb{T}_{k}$ to obtain a modified PHT-spline surface $\mathbf{S}_{k}(s,t)$ (in the beginning, every basis function is new). See the following subsection 4.2 for details. 3. 3.

Find the cells of level $k$ whose fitting errors are greater than $\varepsilon$, and denote them as $\Theta_{k}$. The fitting error over the cell $\theta$ is defined to be $\max_{(s_{j},t_{j})\in\theta}\|\mathbf{P}_{j}-\mathbf{S}_{k}(s_{j},t_{j})\|$. 4. 4.

Label each cell in $\Theta_{k}$ with marks based on discrete curvature information (see the following Subsection 4.3 for details). Subdivide these cells as Algorithm 1. Set $k=k+1$.

4.2 Compute control points

To get a modified PHT-spline surface that fit the given open mesh $\mathcal{M}$, one needs to determine the control points in (4). As a standard way to evaluate the control points, we may solve the following least square optimization problem

[TABLE]

which needs to solve $\mathbf{C}_{j}$ globally.

Similar to [5, 35], an efficient method based on local geometric information can be provided as follows.

Based on the proof of Theorem 6, for each basis vertex $\mathbf{v}_{i}$ and its associated four basis functions $b_{4i+k}(s,t),k=0,1,2,3$,

[TABLE]

where $\mathbf{C}=(\mathbf{C}_{4i},\mathbf{C}_{4i+1},\mathbf{C}_{4i+2},\mathbf{C}_{4i+3})$ is a $3\times 4$ matrix and $\mathbf{B}(\mathbf{v}_{i})$ is defined as (3). On the other hand, by fitting the set of points around $\mathbf{v}_{i}$ with a quadratic surface, $\mathfrak{L}(\mathbf{S}(\mathbf{v}_{i}))$ can be approximated. Hence by (5), the control points $\mathbf{C}$ can be estimated by $\mathcal{L}(\mathbf{S}(\mathbf{v}_{i}))\cdot\mathbf{B}^{-1}.$

4.3 Label the cells

To capture the anisotropic information of the fitting surface, discrete curvature information of $\mathbf{S}_{k}(s,t)$ over those cells is evaluated and used as follows.

For any cell $\theta\in\Theta_{k}$, choose $l$ parametric points $(\overline{s}_{j},\overline{t}_{j})\in\theta$, $j=1,\ldots,l$. Curvatures along the $s$ and $t$ directions at these points $\mathbf{S}_{k}(\overline{s}_{j},\overline{t}_{j})$ are

[TABLE]

Set $K_{s}=\frac{1}{l}\sum_{j=1}^{l}{\kappa_{s}}_{j}$ and $K_{t}=\frac{1}{l}\sum_{j=1}^{l}{\kappa_{t}}_{j}$.

Assume that $K_{t}\neq 0$, then the ratio $\rho=K_{s}/K_{t}$ is chosen to characterize the anisotropic feature. When $\rho>\delta$ or $\rho<1/\delta$, the fitting surface is considered to change sharply along one direction but stay flat along the other direction, which is assumed to possess the anisotropic feature over the cell $\theta$. Hence when $\rho>\delta$, label the cell $\theta$ with ‘V’. Similarly, label the cell $\theta$ with ‘H’ when $\rho<1/\delta$. For the other cases, label the cell $\theta$ with ‘C’. Note that the choice of the threshold $\delta$ may influence the label of anisotropic information on each cell. Therefore, $\delta$ can be assigned depending on the problem being addressed.

4.4 Numerical Examples

Figure 10 illustrates four examples, in which fitting open mesh models with modified PHT-splines are shown. The first two examples are open meshes generated by sampling points on a conic surface and a paraboloid, respectively. The fittings of these two open meshes work as prototypes to verify the efficiency of our modified PHT-splines when dealing with open meshes with anisotropic features. The third model is a practical model. The fourth model is referred to as Example 6 in [9]. We sampled $101\times 101$ data points on a uniform grid from the function

[TABLE]

with $B_{i,7}(u)$ being the $i$-th Bernstein polynomial of degree $7$.

For all four examples, the initial T-meshes are set to $2\times 2$ uniform tensor-product meshes for the first two examples and $5\times 5$ tensor-product meshes for the latter two examples on the square domain $[0,1]\times[0,1]$. The tolerance of the fitting error is $\epsilon=0.1\%$, which refers to the size of the corresponding model’s bounding box. The threshold $\delta=0.5$ is chosen to characterize anisotropic features for all of these four examples. The statistical data for comparison are listed in Table 1.

5 IGA based on modified PHT-splines

In this section, we attempt to use our modified PHT-splines in an isogeometric method to solve elliptic partial differential equations.

5.1 Model problem

Suppose the model problem is an elliptic partial differential equation defined as

[TABLE]

where $\Omega\in\mathbb{R}^{2}$ is a connected, bounded domain with a Lipschitz-continuous boundary $\Gamma=\Gamma_{D}\bigcup\Gamma_{N}$, $\Gamma_{D}\bigcap\Gamma_{N}=\emptyset$, $\mathbf{n}$ is the outward unit normal to $\Gamma_{N}$, $\Gamma_{D}$ is assumed to be closed relative to $\Gamma$ and has a positive length, while $f$ and $h$ are square-integrable on $\Omega$ and $\Gamma_{N}$, respectively.

5.2 Discretization

The framework of isogeometric analysis based on modified PHT-splines is as follows. For more details about isogeometric analysis based on splines, see [25, 39] and references therein.

Suppose the parametrization $\mathbf{G}$ of the physical domain $\Omega$ is defined by

[TABLE]

where $\Omega_{0}=[0,1]\times[0,1]$, $\mathbf{P}_{i}\in\mathbb{R}^{2}$, $b_{i}(s,t)$ is a modified PHT-spline function, $d$ is the number of basis functions.

Let $V(\Omega)=H^{1}(\Omega)$ be the underlying Hilbert space for both the space of test functions and the solution space. The function space $H^{1}(\Omega)$ is defined by

[TABLE]

The weak form solution of problem (5.1) is to find $u\in V=\{v\in H^{1}(\Omega):v|_{\Gamma_{D}}=0\}$ such that

[TABLE]

where $a(,)$ is a bilinear form and $\langle,\rangle$ is a linear functional defined by

[TABLE]

Based on Galerkin’s principle, finite dimensional function space $V_{h}$ is set up to solve the following problem:

Find $u^{h}\in V^{h}$ such that

[TABLE]

where $V^{h}\subset V$ is

[TABLE]

The approximation solution $u^{h}$ can be written as

[TABLE]

where $c_{i},i=1,\ldots,n$ are coefficients need to be determined. Define the stiffness matrix $\mathbf{A}$ by

[TABLE]

and the load vector $\mathbf{F}$ by

[TABLE]

Thus problem (7) is equivalent to the following linear system

[TABLE]

where $\mathbf{c}=(c_{1},c_{2},\ldots,c_{n})$ is the coefficient vector.

5.3 Solving in the adaptive process with modified PHT-splines

Consider the flexibility of the T-meshes and local refinement algorithm for modified PHT-splines, the adaptive procedure consists of the following successive loops

[TABLE]

The essential part of the loops is the estimate and mark step. Follow the well-developed way in IGA [21, 39], the posteriori error is considered here. It should be noted that error estimators for anisotropic refinement have been discussed in the literature on finite element methods [32]. The posteriori error on a cell is

[TABLE]

where $h_{\theta}$ is the diameter of cell $\theta$ and $\Delta=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$. The posteriori error on the modified hierarchical T-mesh $\mathbb{T}$ is the sum of the posteriori errors over all the cells, that is,

[TABLE]

The strategy for the estimate and mark, and refine steps in the loops is as follows.

Find the cells of the current level that need to be subdivided, which are indicated by the the posteriori error estimator (8). Here we use a prescribed threshold to decides the cells. 2. 2.

For each cell that to be subdivided, choose $l$ parametric points $(\overline{s}_{j},\overline{t}_{j})$, $j=1,\ldots,l$ in the cell and evaluate

[TABLE]

Set $\hat{K}_{s}=\sum_{j=1}^{l}{\hat{k}_{sj}}/l,\hat{K}_{t}=\sum_{j=1}^{l}{\hat{k}_{tj}}/l$. If $\hat{K}_{t}\neq 0$, the ratio $\hat{K}_{s}/\hat{K}_{t}$ is adopted to describe the anisotropic features of this cell. If $\hat{K}_{s}/\hat{K}_{t}>\delta_{0}$ or $\hat{K}_{s}/\hat{K}_{t}<\delta_{1}$, where $\delta_{0}>\delta_{1}$ be two prescribed values, the solution $u^{h}(G(s,t))$ is assumed to change sharply along one direction but stay flat along the other direction. In addition, we can treat the case $\hat{K}_{t}=0$ similarly. 3. 3.

Use Algorithm 1 to refine the mesh.

One of the important steps of constructing modified PHT-splines is to label cells in T-meshes, which indicates the subdivision type the cells tend to choose. Similar to the method we used for fitting open mesh in subsection 4.3, for each cell that needs to be subdivided, the discrete second-order partial derivatives along two parametric directions at parametric points are computed and the ratio $\rho$ is prescribed to characterize the anisotropic feature of the solution $u^{h}$.

The boundary conditions in IGA for modified PHT-splines are imposed as follows. 1)In the case of a homogenous boundary condition on $\Gamma_{D}$, for each boundary vertex that is mapped into a point in $\Gamma_{D}$, find its associated four basis functions. Set the coefficients of those functions that not vanish on $\Gamma_{D}$ to zero. 2)In the case of a nonhomogeneous boundary condition on $\Gamma_{D}$, for each boundary vertex that is mapped into a point in $\Gamma_{D}$, find its associated four basis functions and solve the coefficients of these functions by minimizing errors (the least square method may be an option). 3)In the case of a Neumann boundary condition on $\Gamma_{D}$, for each boundary vertex that is mapped into a point in $\Gamma_{D}$, find its adjacent vertices. Construct the equations for the coefficients based on the information of directional derivatives.

5.4 Numerical experiments

Here, an IGA problem based on modified PHT-splines is illustrated. For convenience, DOF is used as the abbreviation for degree of freedom. A cell $\theta$ is marked for refinement, if $\eta_{\theta}>0.0001$. In addition, $\delta=0.5$ is chosen as the threshold to characterize the anisotropic feature.

The example was also considered by Dörfel et al [8] and Kleiss et al. [21]. We solve the Laplace equation (5.1) on the $L$-shaped domain $\Omega=(-1,1)^{2}\backslash[0,1]^{2}$, with boundaries $\Gamma_{D}=\{0\}\times\{0,1\}\cup[0,1]\times\{0\}$, and $\Gamma_{N}=\partial\Omega\backslash\Gamma_{D}$, that is illustrated in Figure 13. Double control points are used to model the corners at $(0,0)$ and $(-1,-1)$, the complete geometry data is referred to Appendix A.2 in [21].

The following function $u$ solves (5.1) and is used as our exact solution:

[TABLE]

Hence, $g_{D}$ and $g_{N}$ in (5.1) are determined by the exact solution $u$.

We start from a $3\times 3$ tensor-product mesh. The refined meshes on the physical domains solved by PHT-splines and modified PHT-splines are shown in Figure 13 and Figure 13 respectively.

A comparison of the solution of the PHT-splines and modified PHT-splines is provided in Figure 14. The comparison shows that we can achieve the better accuracy with the same convergence rate.

6 Conclusion and future work

In this paper, we extend PHT-splines to modified PHT-splines, a special type of polynomial spline defined over modified hierarchical T-meshes. Simultaneously, the basis functions of modified PHT-splines that we constructed inherit the beneficial properties of PHT-splines, i.e., nonnegativity, partition of unity, local support, and linear independence. A refinement strategy based on neighborhood relations and labels of cells is also presented, in which the labels can be estimated from the geometric information in practical problems.

Modified PHT-splines have been applied to fitting open meshes and isogeometric analysis for the elliptic partial differential equation. Numerical results show that modified PHT-splines have advantages when applied to problems with anisotropic features.

It should be pointed out that only cells of lasted level are marked for refinement in our scheme. This assumption may limit our splines to be used in further applications. Especially, for the problems of isogeometric analysis, there is no guarantee that only cells of lasted level need to be marked in some practical problem, although our refinement scheme has good performance in many numerical examples. In the future, we will focus on improving the algorithm when it is applied to isogeometric analysis. In addition, we will further explore the applications of modified PHT-splines in isogeometric analysis. To exactly represent common geometric objects such as circles, cylinders, spheres, and ellipsoids, rational forms of modified PHT-splines may be required. Furthermore, the generalization of modified PHT-splines in three-dimensional space is worthy of consideration.

Acknowledgements

We would like to thank the anonymous referees for providing us with constructive comments and suggestions.