THU-Splines: Highly Localized Refinement on Smooth Unstructured Splines

To start with, a set of $n^{0}$ univariate B-splines $\mathcal{B}^{0}=\{B_{i}^{0}\}_{i=1}^{n^{0}}$ of degree $p$ ($n^{0},\,p\in\mathbb{N}^{+}$), defined on the knot vector $\Xi^{0}=\{\xi_{i}^{0}\}_{i=1}^{n^{0}+p+1}$, is given and initialized as the initial set of THB-spline basis functions, that is,

$$\mathcal{H}^{0}=\mathcal{B}^{0}.$$

(1)

The initial domain serves as the entire domain $\Omega$ of THB-splines, i.e., $\Omega\equiv\Omega^{0}=(\xi_{1}^{0},\xi_{n+p+1}^{0})$.

To facilitate explanation, we introduce a series of nested spaces spanned by the B-splines,

$$\mathcal{B}^{0}\subset\mathcal{B}^{1}\subset\cdots\subset\mathcal{B}^{\ell_{\max}},$$

(2)

where $\ell_{\max}$ is the predefined maximum level, and the set $\mathcal{B}^{\ell}=\{B_{i}^{\ell}\}_{i=1}^{n^{\ell}}$ ($1\leq\ell\leq\ell_{\max}$) is defined on the knot vector $\Xi^{\ell}=\{\xi_{i}^{\ell}\}_{i=1}^{n^{\ell}+p+1}$. Note that the degree $p$ is fixed at all levels. As will become clear shortly, only a small subset of $\mathcal{B}^{\ell}$ is used and the full set $\mathcal{B}^{\ell}$ is never generated in practice. $\Xi^{\ell}$ is obtained by bisecting all the nonzero knot spans of $\Xi^{\ell-1}$, so knot vectors are nested as well,

$$\Xi^{0}\subset\Xi^{1}\subset\cdots\subset\Xi^{\ell_{\max}}.$$

(3)

A hierarchical construction involves a nested sequence of hierarchical domains,

$$\Omega\equiv\Omega^{0}\supseteq\Omega^{1}\supseteq\cdots\supseteq\Omega^{\ell_{\max}}.$$

(4)

This nested sequence implies a series of local refinements and it drives where to perform local refinement. For instance, once a set of THB-splines is initialized, we use $\mathcal{H}^{0}$ to solve the problem of interest. According to a certain a posteriori error estimator, a subdomain of $\Omega^{0}$, which is in fact the domain at Level 1 (i.e., $\Omega^{1}$), is marked to perform local refinement. Once $\Omega^{1}$ is determined, THB-splines $\mathcal{H}^{1}$ can be constructed accordingly.

THB-splines are constructed level by level, so we only need to study a generic two-level construction for selection and truncation: given the THB-splines $\mathcal{H}^{\ell}$ that has been constructed up to Level $\ell$ ($0\leq\ell<\ell_{\max}$) and $\Omega^{\ell+1}$ (i.e., the subdomain of $\Omega^{\ell}$ to be refined), we discuss how to build $\mathcal{H}^{\ell+1}$.

The selection mechanism, first proposed in ref:kraft97 and later generalized in ref:vuong11 , is aimed to ensure linear independence of resulting hierarchical basis functions. It is based on the relation between the support of a candidate B-spline and active subdomains of $\Omega^{\ell}$ and $\Omega^{\ell+1}$. The selection at Level $\ell$ states that a B-spline $B_{i}^{\ell}$ is selected only if its support has a nonzero intersection with $\Omega_{a}^{\ell}:=\Omega^{\ell}\backslash\Omega^{\ell+1}$, the active subdomain at Level $\ell$. All the selected Level-$\ell$ B-splines are collected in $\mathcal{B}_{a}^{\ell}$,

$$\mathcal{B}_{a}^{\ell}=\{B_{i}^{\ell}\in\mathcal{B}^{\ell}:\,\mathrm{supp}B_{i}^{\ell}\cap\Omega_{a}^{\ell}\neq\varnothing\}.$$

(5)

where $\mathrm{supp}B_{i}^{\ell}:=(\xi_{i}^{\ell},\xi_{i+p+1}^{\ell})$ is the support of $B_{i}^{\ell}$. On the other hand, at Level $\ell+1$, a B-spline is selected only if its support is fully contained in $\Omega^{\ell+1}$. All the selected Level-($\ell+1$) B-splines form the set $\mathcal{B}_{a,0}^{\ell+1}$,

$$\mathcal{B}_{a,0}^{\ell+1}=\{B_{i}^{\ell+1}\in\mathcal{B}^{\ell+1}:\,\mathrm{supp}B_{i}^{\ell+1}\subseteq\Omega^{\ell+1}\}.$$

(6)

Selected B-splines are called active, and the complementary set of B-splines (i.e., $\mathcal{B}^{\ell}\backslash\mathcal{B}_{a}^{\ell}$ and $\mathcal{B}^{\ell+1}\backslash\mathcal{B}_{a,0}^{\ell+1}$) are passive. A prefix “H” will be added (i.e., H-active or H-passive) to emphasize the hierarchical refinement when the context is not evident; see Section 4. The subscript “0” in $\mathcal{B}_{a,0}^{\ell+1}$ indicates the initial selection at Level $\ell+1$ with the entire $\Omega^{\ell+1}$ being active. Certain functions in $\mathcal{B}_{a,0}^{\ell+1}$ may latter become passive when Level $\ell+2$ is constructed from Level $\ell+1$ and the active Level-($\ell+1$) subdomain changes from $\Omega^{\ell+1}$ to $\Omega_{a}^{\ell+1}\equiv\Omega^{\ell+1}\backslash\Omega^{\ell+2}$. $\mathcal{I}^{\ell}$, $\mathcal{A}^{\ell}$, $\mathcal{I}^{\ell+1}$ and $\mathcal{A}_{0}^{\ell+1}$ are introduced to denote the index sets of $\mathcal{B}^{\ell}$, $\mathcal{B}_{a}^{\ell}$, $\mathcal{B}^{\ell+1}$ and $\mathcal{B}_{a,0}^{\ell+1}$, respectively.

The truncation mechanism, which enables reduced overlapping and partition of unity, plays the key role in THB-splines. We need the refinability relationship to explain the idea. Refinability states that a Level-$\ell$ B-spline $B_{i}^{\ell}$ can be expressed as a linear combination of certain Level-($\ell+1$) B-splines,

$$B_{i}^{\ell}=\sum_{j\in\mathcal{C}_{i}^{\ell+1}}h_{ij}^{\ell+1}B_{j}^{\ell+1},\quad\text{with}\quad\mathcal{C}_{i}^{\ell+1}:=\{j\in\mathcal{I}^{\ell+1}:\,\mathrm{supp}B_{j}^{\ell+1}\subset\mathrm{supp}B_{i}^{\ell}\},$$

(7)

where $h_{ij}^{\ell+1}$ are coefficients obtained from the knot insertion algorithm ref:boehm80 , and $B_{j}^{\ell+1}$ ($j\in\mathcal{C}_{i}^{\ell+1}$) are called the children (or H-children, children in the hierarchical refinement) of $B_{i}^{\ell}$. In fact, refinability is the reason for Eq. (2) to hold.

The truncation mechanism applies to $B_{i}^{\ell}$ when some of its children are active, and it discards the active children from Eq. (7). As a result, the truncated function of $B_{i}^{\ell}$ is defined as

$$\mathrm{trun}_{H}B_{i}^{\ell}=\sum_{j\in\mathcal{C}_{i}^{\ell+1}\backslash\mathcal{A}_{0}^{\ell+1}}h_{ij}^{\ell+1}B_{j}^{\ell+1},$$

(8)

which is also referred to as the inter-level truncation in this paper. Equivalently speaking, the coefficient $h_{ij}^{\ell+1}$ in Eq. (7) is set to be $0$ by truncation if $B_{j}^{\ell+1}$ is active. Note that when all the children are passive (i.e., $\mathcal{C}_{i}^{\ell+1}\cap\mathcal{A}_{0}^{\ell+1}=\varnothing$), we have $\mathcal{C}_{i}^{\ell+1}\backslash\mathcal{A}_{0}^{\ell+1}\equiv\mathcal{C}_{i}^{\ell+1}$ and thus $\mathrm{trun}_{H}B_{i}^{\ell}\equiv B_{i}^{\ell}$. In other words, Eq. (8) can represent a generic definition including both truncated and non-truncated B-splines.

In summary, after selection of $\mathcal{B}_{a}^{\ell}$ and $\mathcal{B}_{a,0}^{\ell+1}$, truncation is applied to certain B-splines in $\mathcal{B}_{a}^{\ell}$, leading to a set of truncated (active) B-splines at Level $\ell$,

$$\mathcal{T}_{a}^{\ell}=\{\mathrm{trun}_{H}B_{i}^{\ell}:\,B_{i}^{\ell}\in\mathcal{B}_{a}^{\ell}\}.$$

(9)

The THB-spline basis functions up to Level $\ell+1$ are now obtained as

$$\mathcal{H}^{\ell+1}=\mathcal{H}^{\ell-1}\cup\mathcal{T}_{a}^{\ell}\cup\mathcal{B}_{a,0}^{\ell+1},$$

(10)

where note that $\mathcal{H}^{\ell-1}:=\varnothing$ when $\ell=0$.

In Eq. (10), the local refinement is realized by enriching basis functions and it solely depends on $\mathcal{B}_{a,0}^{\ell+1}$. $\mathcal{H}^{\ell-1}$ and $\mathcal{T}_{a}^{\ell}$ are merely adjustments of existing basis functions to accommodate $\mathcal{B}_{a,0}^{\ell+1}$. Therefore, local refinement is effective only when $\mathcal{B}_{a}^{\ell+1}\neq\varnothing$; otherwise we end up with $\mathcal{H}^{\ell+1}\equiv\mathcal{H}^{\ell}$. A non-empty $\mathcal{B}_{a,0}^{\ell+1}$ implies that $\Omega^{\ell+1}$ must be able to cover the support of at least one Level-($\ell+1$) B-spline; see Eq. (6). This is indeed the only requirement for local refinement of THB-splines. It has a definite termination condition and it indicates a highly localized refinement with no propagation beyond the region of interest.

We first introduce several terminologies of unstructured meshes with the reference to Fig. 1. An element is called an irregular element if any of its vertices²²2We use “vertex”, “edge” and “face” to emphasize mesh connectivity (or topology), whereas using “point” to carry the position or geometry information. “Element” and “face” are used interchangeably. is an EV; otherwise it is regular. An element is a boundary element if any of its vertices lies on the boundary; otherwise it is an interior element. The one-ring neighborhood of a vertex includes all the elements sharing this vertex; recursively, the $n$-ring ($n\geq 2$) neighborhood is the ($n-1$)-ring neighborhood together with those sharing vertices with the ($n-1$)-ring neighborhood. The neighborhood of an element can be defined similarly. We adopt the following assumptions regarding the input mesh to simplify discussion.

• •

  An irregular element can only have one EV;

• •

  A boundary vertex can only be shared by either one or two elements; and

• •

  An irregular element must not lie in the three-ring neighborhood of any boundary vertex.

Each vertex corresponds to a control point, and we call such points V-points (i.e., vertex points). Each edge is assigned with a non-negative real number called the knot span or knot interval, which is used to define element (parametric) domains as well as spline functions. For simplicity, we adopt a semi-uniform configuration of knot intervals where all edges have a unit knot interval except for those perpendicular to the boundary, which have a zero knot interval. As a result, edges around EVs all have the same knot intervals, which indeed is a necessary assumption for the D-patch construction ref:reif97 .

Note that given the third assumption regarding the input mesh, Bézier extraction is never needed for elements near the boundary. Without loss of generality, we next explain how to extract a Bézier representation for an irregular element $\Omega_{e}$ ref:scott13 . Fig. 2(a) shows as a generic irregular element $\Omega_{e}$, where $P_{e,i}$ ($i=1,\ldots,4$) are its 4 vertices, and $P_{e,1}$ is an EV. In the bicubic case, Bézier extraction is to find the $4\times 4$ Bézier control points $Q_{e,j}$ ($j=1,\ldots,16$) corresponding to $\Omega_{e}$. They can be divided into face, edge and vertex Bézier points, and they are called FB-points, EB-points, and VB-points, respectively; see Fig. 2(b). The 4 FB-points are computed as convex combinations of the 4 element corners, that is,

$$\begin{bmatrix}Q_{e,6}\\ Q_{e,7}\\ Q_{e,11}\\ Q_{e,10}\\ \end{bmatrix}=\begin{bmatrix}a&b&c&b\\ b&a&b&c\\ c&b&a&b\\ b&c&b&a\\ \end{bmatrix}\begin{bmatrix}P_{e,1}\\ P_{e,2}\\ P_{e,3}\\ P_{e,4}\\ \end{bmatrix},\quad\text{with}\quad a=\frac{4}{9},\,b=\frac{2}{9},\,c=\frac{1}{9}.$$

(11)

As shown in Fig. 2(c, d), EB- and VB-points are simply an average of their neighboring FB-points, for examples,

	$\displaystyle Q_{e,9}$	$\displaystyle=Q_{e+1,3}=\frac{1}{2}Q_{e,10}+\frac{1}{2}Q_{e+1,7},$		(12)
	$\displaystyle Q_{e,1}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}Q_{i,6},$

where $N\in\mathbb{N}^{+}$ is the number of faces sharing $P_{e,1}$, also called the valence of $P_{e,1}$.

Now let $\bm{Q}_{e}:=[Q_{e,1},\ldots,Q_{e,16}]^{T}$ denote the vector of Bézier points of $\Omega_{e}$, and $\bm{P}_{v}$ and $\bm{P}_{f}$ contain the V-points and FB-points of the one-ring neighborhood of $\Omega_{e}$, respectively. From Eqs. (11, 12), we observe that $\bm{Q}_{e}$ can be obtained through either

$$\bm{Q}_{e}=\bm{A}_{e}\bm{P}_{f},$$

(13)

or

$$\bm{Q}_{e}=\bm{E}_{e}\bm{P}_{v},$$

(14)

where the matrix $\bm{A}_{e}$ contains the averaging coefficients in Eq. (12), and $\bm{E}_{e}$ is the extraction matrix with coefficients from both Eqs. (11, 12). When $N=4$, Eqs. (13, 14) indicate Bézier extraction for a regular bicubic B-spline element.

The Bézier representation of $\Omega_{e}$ is written as

$$\bm{S}(u,v)|_{\Omega_{e}}=\bm{Q}_{e}^{T}\bm{b}(u,v),\quad(u,v)\in[0,1]^{2},$$

(15)

where $\bm{b}(u,v):=[b_{1}(u,v),\ldots,b_{16}(u,v)]^{T}$ is the vector of bicubic Bernstein polynomials. $b_{j}(u,v)$ is expressed as,

$$b_{j}(u,v)=N_{(j-1)\%4}(u)N_{(j-1)/4}(v),\quad(u,v)\in[0,1]^{2},$$

(16)

where “$\%$” and “/” denote remainder and quotient operations, respectively; and

$$N_{k}(t)=\binom{3}{k}t^{k}(1-t)^{3-k},\quad k=0,\ldots,3.$$

(17)

Alternatively, the surface patch $\bm{S}(u,v)|_{\Omega_{e}}$ can also be expressed in terms of $\bm{P}_{v}$ and the associated spline functions $\bm{B}_{v}(u,v)$,

$$\bm{S}(u,v)|_{\Omega_{e}}=\bm{P}_{v}^{T}\bm{B}_{v}(u,v).$$

(18)

On the other hand, Eqs. (15, 18) must be equivalent, leading to

$$\bm{P}_{v}^{T}\bm{B}_{v}(u,v)\equiv\bm{Q}_{e}^{T}\bm{b}(u,v)=\bm{P}_{v}^{T}\bm{E}_{e}^{T}\bm{b}(u,v),$$

(19)

where we have used Eq. (14) to arrive at the last term. According to the arbitrariness of $\bm{P}_{v}$, we conclude that

$$\bm{B}_{v}(u,v)=\bm{E}_{e}^{T}\bm{b}(u,v).$$

(20)

Similarly, the surface patch can also be written with FB-points $\bm{P}_{f}$,

$$\bm{S}(u,v)|_{\Omega_{e}}=\bm{P}_{f}^{T}\bm{B}_{f}(u,v),\quad\text{with}\quad\bm{B}_{f}(u,v)=\bm{A}_{e}^{T}\bm{b}(u,v).$$

(21)

The central idea of Bézier extraction is to express each spline function as a linear combination of Bernstein polynomials, which enables general spline functions to be easily integrated with existing finite element and isogeometric codes. While in regular elements, $\bm{B}_{v}$ and $\bm{B}_{f}$ are actually bicubic $C^{2}$ and $C^{1}$ B-splines, respectively, they are only $C^{0}$-continuous in irregular elements. The D-patch construction essentially adjusts $\bm{E}_{e}^{T}$ (or $\bm{A}_{e}^{T}$) to yield smooth spline functions, which we discuss in the following.

In a local D-patch construction, the D-patch method is only adopted in regions around EVs, where we need to add necessary face-based control points and define their associated splines. Face-based control points are referred to as F-points in the following. Splines associated with V-points and F-points are referred to as V-splines and F-splines, respectively.

Four F-points are then added to each enriched element. In fact, F-points are equivalent to FB-points, and they are obtained according to Eq. (11) in the input mesh. As a result, the control mesh of a local D-patch representation consists of both V-points and F-points; see Fig. 3. However, while all the F-points of enriched elements are active in the sense that they are used in the local D-patch representation, not all V-points are active. A V-point is passive if all its one-ring neighboring elements are enriched elements; otherwise it is active. A prefix “D” will be added to emphasize the context of local D-patch representation when necessary to distinguish those in hierarchical refinement. Clearly, all EVs are passive according to this definition. Essentially, passive V-points are replaced by their neighboring F-points. Note that such a control mesh corresponds to the analysis space discussed in ref:toshniwal17 .

In the D-patch method, a smooth representation is achieved by adjusting all $\bar{Q}_{e,k}^{11}$, $e\in\{1,\ldots,N\}$, $k\in\{1,2,\ldots,16\}\backslash\{11,12,15,16\}$, which are collected into a vector $\bar{\bm{Q}}_{s}^{11}$. This is done through a linear projection matrix $\bm{\Pi}$,

$$\tilde{\bm{Q}}_{s}^{11}=\bm{\Pi}\bar{\bm{Q}}_{s}^{11}.$$

(23)

As a result, for every $e\in\{1,\ldots,N\}$, $\tilde{Q}_{e,1}^{11}=\tilde{Q}_{e,2}^{11}=\tilde{Q}_{e,5}^{11}=\tilde{Q}_{e,6}^{11}\equiv\bar{Q}_{EV}$, and the five points $\tilde{Q}_{e,3}^{11}$, $\tilde{Q}_{e,7}^{11}$, $\tilde{Q}_{e,9}^{11}$, $\tilde{Q}_{e,10}^{11}$ and $\bar{Q}_{EV}$ are coplanar; see indices in Fig. 4. The reminder of $\bar{\bm{Q}}_{e}^{ij}$ stays the same, and the resulting Bézier surface patches are smoothly joined together around the EV. The smoothing matrix $\bm{\Pi}$ corresponds to the idempotent version of the analysis space ref:casquero20 ; ref:toshniwal17 . Interested readers are referred there for further details. Note that although $\bm{\Pi}$ contains negative coefficients and thus leads to slightly negative spline functions, it guarantees refinability, which is the fundamental property required in a hierarchical spline construction.

Now we discuss the F-splines defined on an irregular element $\Omega_{e}$, which are associated with the F-points of the one-ring neighborhood of $\Omega_{e}$. They are piecewise-defined on the four subelements ($\Omega_{e}^{ij}$, $i,j=1,2$) as

$$\tilde{\bm{B}}_{f}(u,v)|_{\Omega_{e}^{ij}}=\bm{A}_{e}^{T}\bm{D}_{ij}^{T}\bm{\Pi}_{e,ij}^{T}\bm{b}(u,v).$$

(24)

where the entries of $\bm{\Pi}_{e,ij}$ come from $\bm{\Pi}$ with an proper arrangement for $\Omega_{e}^{ij}$, and in particular, $\bm{\Pi}_{e,22}$ is an identity matrix. Comparing Eqs. (21, 24), we observe that by introducing matrices $\bm{D}_{ij}$ and $\bm{\Pi}_{e,ij}$, $C^{1}$ continuity is built into $\tilde{\bm{B}}_{f}$. Moreover, we have $\tilde{\bm{B}}_{f}\equiv\bm{B}_{f}$ outside of the irregular element $\Omega_{e}$ thanks to the $2\times 2$ split, or equivalently speaking, the presence of $\bm{D}_{ij}$. Mostly importantly, it has also been proved that $\tilde{\bm{B}}_{f}$ are refinable in the sense of Eq. (7); see ref:tnguyen16 . In the reminder of the paper, we omit “$\sim$” over $\tilde{\bm{B}}_{f}$ to keep notations simple.

Next, we introduce how to smoothly join enriched and non-enriched elements in a local D-patch representation. The following terms are needed to facilitate the discussion. An interior edge is called an interface edge if it is shared by an enriched element and a non-enriched element. Vertices of interface edges are interface vertices. An element is then called a transition element if it contains an interface vertex. Note that transition elements can be regular or irregular; see Fig. 1. Now let $\mathcal{V}$ and $\mathcal{F}$ denote the index sets of V-points and F-points, respectively, and $\mathcal{V}_{a}$ be the active subset of $\mathcal{V}$.

We now discuss the splines associated with the active V-points ($\mathcal{V}_{a}$), which include interface points ($\mathcal{V}_{t}$) and the active regular points ($\mathcal{V}_{r}$) that do not lie on the interface, that is, $\mathcal{V}_{a}=\mathcal{V}_{t}\cup\mathcal{V}_{r}$. Splines associated with $\mathcal{V}_{r}$ are simply bicubic $C^{2}$ B-splines, whereas those associated with $\mathcal{V}_{t}$ need truncation to enable a smooth transition between enriched and non-enriched elements. It is related but not the same as that introduced in the context of THB-splines. It was first proposed in blended B-splines ref:wei18 to recover the optimal convergence property when using unstructured hexahedral meshes. To distinguish, we refer to the truncation in a local D-patch representation as the same-level truncation. Recall that it is called the inter-level truncation in the context of THB-splines.

For the V-spline $B_{v,i}$ associated with an interface point $P_{i}$ ($i\in\mathcal{V}_{t}$), we can express it in terms of its one-ring neighboring F-splines $B_{f\!,j}$,

$$B_{v,i}=\sum_{j\in\mathcal{F}_{i}}d_{ij}B_{f\!,j},$$

(25)

where the coefficients $d_{ij}$ are obtained from Eq. (11) and they are illustrated in Fig. 5(a), $\mathcal{F}_{i}\subseteq\mathcal{F}$ is the index set of the F-points on the one-ring neighborhood of $P_{i}$, and $B_{f\!,j}$ are the children (or D-children, children in the local D-patch representation) of $B_{v,i}$. $B_{f\!,j}$ is defined by Eq. (24) if its F-point lies in an irregular element; otherwise it is defined by Eq. (21). Note that not all $B_{f\!,j}$ are associated with the F-points of enriched elements. As $P_{i}$ is an interface point, there must exist at least one non-enriched element in its one-ring neighborhood. F-points on non-enriched elements are not part of the control mesh in a local D-patch representation, and thus they are passive. In contrast, F-points of enriched elements are active.

The same-level truncation for $B_{v,i}$ states that $d_{ij}$ is set to be 0 if $B_{f\!,j}$ is active. In other words, active children of $B_{v,i}$ are discarded,

$$\mathrm{trun}_{D}B_{v,i}=\sum_{j\in\mathcal{F}_{i}\backslash\mathcal{F}_{a}}d_{ij}B_{f\!,j},$$

(26)

where $\mathcal{F}_{a}$ is the index set of active F-points. As a result of truncation, $\mathrm{trun}_{D}B_{v,i}$ is always a linear combination of regular $C^{1}$ B-splines. Fig. 5(b) shows an example of $\mathrm{trun}_{D}B_{v,i}$, where the red squares are active F-points as the gray element is enriched whereas the blue squares are passive and only serve as an auxiliary purpose to define $\mathrm{trun}_{D}B_{v,i}$. Therefore, the splines associated with interface points are truncated B-splines defined according to Eq. (26). From an element-wise perspective, truncation has no influence outside of transition elements in the sense that $\mathrm{trun}_{D}B_{v,i}\equiv B_{v,i}$.

On the other hand, given an active non-interface point $P_{k}$ ($k\in\mathcal{V}_{r}$), its associated spline $B_{v,k}$ is a $C^{2}$ B-spline and can be also expressed in the form of Eq. (26), where, however, truncation is trivial as $\mathcal{F}_{k}\backslash\mathcal{F}_{a}=\mathcal{F}_{k}$, because we have $\mathcal{F}_{k}\cap\mathcal{F}_{a}=\varnothing$ by definition. In other words, $\mathrm{trun}_{D}B_{v,k}\equiv B_{v,k}$ holds for $k\in\mathcal{V}_{r}$.

To this end, all spline functions of a local D-patch construction have been discussed, which form a partition of unity and are linearly independent ref:casquero20 . The local D-patch representation is written as

$$S=\sum_{i\in\mathcal{F}_{a}}Q_{i}\,B_{f\!,i}+\sum_{i\in\mathcal{V}_{t}}P_{i}\,\mathrm{trun}_{D}B_{v,i}+\sum_{i\in\mathcal{V}_{r}}P_{i}\,B_{v,i}=\sum_{i\in\mathcal{F}_{a}}Q_{i}\,B_{f\!,i}+\sum_{i\in\mathcal{V}_{a}}P_{i}\,\mathrm{trun}_{D}B_{v,i},$$

(27)

where we recall that $\mathcal{V}_{a}=\mathcal{V}_{r}\cup\mathcal{V}_{t}$, and $Q_{i}$ and $P_{i}$ are (active) F-points and V-points, respectively. The corresponding matrix form is given as,

$$\bm{S}=\begin{bmatrix}\bm{Q}^{T}&\bm{P}^{T}\end{bmatrix}\begin{bmatrix}\bm{B}_{f}\\ \bm{B}_{v}\end{bmatrix},$$

(28)

where $\bm{B}_{v}$ and $\bm{B}_{f}$ are vectors of V-splines and F-splines, respectively, and $\bm{B}_{v}$ includes both regular $C^{2}$ B-splines ($\mathcal{V}_{r}$) and truncated B-splines ($\mathcal{V}_{t}$).

Finally, we discuss how to perform global $h$-refinement for a local D-patch representation. In the absence of EVs, $h$-refinement of a B-spline patch is straightforward through the knot insertion algorithm, where the underlying geometric mapping stays the same during refinement. We call such a refinement the consistent refinement. However, consistent refinement is usually not available when EVs are involved. Currently, only two methods³³3We only focus on methods based on quadrilaterals. have this property: the D-patch method and subdivision. In fact, subdivision implies a large family of methods that feature an infinite piecewise definition of splines in irregular regions. Interested readers are referred to ref:xli19 ; ref:wei20 for some of the latest developments that aim to recover optimal convergence rates.

Consistent global refinement of a local D-patch representation consists of a topological step and a geometric step. The topological step deals with mesh connectivity as well as inheritance of enriched elements, whereas the geometric step computes control point coordinates of the refined mesh. Elementwise definitions of basis functions are discussed as well after the topological step.

The computation of point coordinates has two cases: B-spline refinement for non-enriched elements, and Bézier refinement for enriched elements. First, for a surface patch restricted to a non-enriched element (no matter if the element is transition or not), it is simply a B-spline patch when we only focus on its V-points. All the V-points of refined patches are obtained through the knot insertion algorithm applied to such a B-spline patch. Moreover, when the non-enriched element is a transition element, recall that some of its child elements are enriched elements, whose F-points are added according to Eq. (11).

Second, for a surface patch restricted to an enriched element, the corresponding local D-patch representation is first converted to its Bézier representation. Refinement is then applied to this Bézier patch using de Casteljau’s algorithm. In each of the resulting child Bézier patches, the 4 interior Bézier control points are added to the corresponding element (which must be an enriched element by definition) as the F-points.

To this end, both basis functions and control points are defined for a refined mesh. Its local D-patch representation has the same form as Eq. (28).

THU-splines are initialized with the local D-patch construction on an input unstructured quad mesh $\mathcal{M}^{0}$. For the ease of discussion, we introduce a series of consecutively refined meshes, $\{\mathcal{M}^{\ell}\}_{\ell=0}^{\ell_{\max}}$, where $\ell_{\max}\in\mathbb{N}^{+}$ is the predefined maximum number of refinements, and $\mathcal{M}^{\ell}$ is a global refinement of $\mathcal{M}^{\ell-1}$ ($1\leq\ell\leq\ell_{\max}$) according to Section 3.4. Note that global refinement is never actually performed and these globally refined meshes are introduced only to simplify explanation. Correspondingly, a local D-patch basis $\mathcal{B}^{\ell}$, which only consists of D-active splines, is constructed on each $\mathcal{M}^{\ell}$. The initial THU-spline basis $\mathcal{H}^{0}$ is given by $\mathcal{B}^{0}$. Moreover, with a closely related proof (Proposition 4.4, ref:toshniwal17 ), we conjecture that refinability holds for $\{\mathcal{B}^{\ell}\}$, for which we will show numerical evidence in Section 5.

Let $\mathcal{F}^{\ell}$ and $\mathcal{V}_{a}^{\ell}$ denote the index sets of F-points and active V-points of $\mathcal{M}^{\ell}$, respectively. Note that $\mathcal{F}^{\ell}=\mathcal{F}_{a}^{\ell}\cup\mathcal{F}_{p}^{\ell}$, where passive F-points ($\mathcal{F}_{p}^{\ell}$) are introduced as an auxiliary means to define truncated B-splines associated with interface points, and splines associated with $\mathcal{F}_{p}^{\ell}$ do not belong to $\mathcal{B}^{\ell}$.

The same as in Section 2, we focus on constructing $\mathcal{H}^{\ell+1}$ from $\mathcal{H}^{\ell}$. Following Eq. (4), we assume that a series of nested hierarchical domains $\{\Omega^{\ell}\}_{\ell=0}^{\ell_{\max}}$ is given to guide local refinement. However, we need to define the “$\supset$” operation in Eq. (4) for unstructured meshes. Let $\Omega^{0}:=\{\Omega_{e}^{0}\}_{e=1}^{n^{0}}$ be a collection of all $n^{0}$ quad elements in $\mathcal{M}^{0}$. $\Omega^{\ell}$ ($1\leq\ell\leq\ell_{\max}$) is defined similarly but does not necessarily contain all elements of $\mathcal{M}^{\ell}$. $\Omega^{\ell}\supset\Omega^{\ell+1}$ means that every element in $\Omega^{\ell+1}$ is a child element of a certain element in $\Omega^{\ell}$. In other words, $\Omega^{\ell+1}$ is obtained by locally refining certain elements in $\Omega^{\ell}$; after refinement, the parent elements of those in $\Omega^{\ell+1}$ are no longer used and thus they are marked as passive. The active subdomain of $\Omega_{a}^{\ell}$ consists of all active elements at Level $\ell$.

It is left to define the support of every basis function to apply the selection mechanism. Given a generic basis function $B_{i}^{\ell}$ (whose control point is denoted as $P_{i}^{\ell}$), its support, $\mathrm{supp}B_{i}^{\ell}$, is defined as a collection of certain elements in $\mathcal{M}^{\ell}$. There are four cases depending on the type of $P_{i}^{\ell}$: (1) a V-point that is not on the interface, (2) an F-point on a regular element, (3) an F-point on an irregular element, and (4) an interface point; see Fig. 6.

In Case (1), $B_{i}^{\ell}$ is a bicubic $C^{2}$ B-spline, so $\mathrm{supp}B_{i}^{\ell}$ is the two-ring neighborhood of $P_{i}^{\ell}$. In Case (2), $B_{i}^{\ell}$ is a bicubic $C^{1}$ B-spline with doubly repeated knots, and $\mathrm{supp}B_{i}^{\ell}$ is the one-ring neighborhood of $P_{i}^{\ell}$. When $P_{i}^{\ell}$ is an F-point, its one-ring neighborhood means the one-ring neighborhood of its closest V-point. As a special case, the F-point $Q_{e,11}$ in Fig. 2(b), although being on an irregular element, also falls into Case (2), whereas the other three ($Q_{e,6}$, $Q_{e,7}$ and $Q_{e,10}$) belong to Case (3).

Case (3) must involve a neighboring EV, and $B_{i}^{\ell}$ is defined as a D-patch spline. $\mathrm{supp}B_{i}^{\ell}$ consists of two parts. First, the same as Case (2), $B_{i}^{\ell}$ has support on the one-ring neighborhood of $P_{i}^{\ell}$. Moreover, $B_{i}^{\ell}$ is defined using the smoothing matrix such that it also has support on the one-ring neighborhood of the EV. The union of the two one-ring neighborhoods yields $\mathrm{supp}B_{i}^{\ell}$.

In Case (4), $B_{i}^{\ell}$ (more precisely, $\mathrm{trun}_{D}B_{i}^{\ell}$) is a truncated B-spline that is expressed as a linear combination of certain D-passive $C^{1}$ B-splines. Therefore, its support is the union of all the one-ring neighborhoods of these D-passive F-points.

With the hierarchical domains and the support of splines well defined, the selection mechanism follows the same as those in Eqs. (5, 6) to collect all the H-active splines as $\mathcal{B}_{a}^{\ell}$ and $\mathcal{B}_{a,0}^{\ell+1}$; see an example in Fig. 7. Note that during selection, only D-active splines at each level can be selected to be H-active or H-passive, whereas D-passive ones only serve as an auxiliary purpose to help define truncated B-splines and they never enter the selection procedure. In other words, D-passive functions are neither H-active nor H-passive.

Finally, recall that $\Omega^{\ell+1}$ needs to be large enough such that at least one Level-($\ell+1$) spline function can be added to $\mathcal{H}^{\ell+1}$. In the THU-spline construction, this condition can be easily guaranteed if the parent elements of those in $\Omega^{\ell+1}$ are a collection of one-ring neighborhoods of certain Level-$\ell$ points. In this manner, $\Omega^{\ell+1}$ at least contains the two-ring neighborhood of a Level-($\ell+1$) point. As the largest support in the above four cases is a two-ring neighborhood, a Level-($\ell+1$) spline can always be added.

The inter-level truncation follows the same as in Eq. (8). Certain H-active Level-$\ell$ functions are truncated. As a result, $\mathcal{B}_{a}^{\ell}$ is changed to $\mathcal{T}_{a}^{\ell}$, and $\mathcal{H}^{\ell+1}$ is constructed according to Eq. (10). Although the expressions stay the same as those in THB-splines, evaluation of THU-splines becomes much more involving because there are various types of splines rather than a single type as in THB-splines.

Without loss of generality, we study the THU-splines defined on a transition element $\Omega_{e}^{\ell+1}$, where we consider two levels (Levels $\ell$ and $\ell+1$) to simplify explanation. Following ref:bornermann13 , we evaluate all levels of H-active functions ($\bm{H}_{a}^{\ell}$ and $\bm{H}_{a}^{\ell+1}$) with the functions (both H-active and H-passive) at the highest level, and we write it in the matrix form,

$$\bm{H}_{a}:=\begin{bmatrix}\bm{H}_{a}^{\ell}\\ \bm{H}_{a}^{\ell+1}\\ \end{bmatrix}=\begin{bmatrix}\bm{0}&\bm{C}\\ \bm{I}&\bm{0}\\ \end{bmatrix}\begin{bmatrix}\bm{H}_{a}^{\ell+1}\\ \bm{H}_{p}^{\ell+1}\\ \end{bmatrix},$$

(30)

where the coefficients of $\bm{C}$ come from Eq. (8), and spline functions in $\bm{H}_{a}^{\ell}$ and $\bm{H}_{*}^{\ell+1}$ ($*\in\{a,p\}$) are from $\mathcal{B}_{a}^{\ell}$ and $\mathcal{B}^{\ell+1}$, respectively. Note that the inter-level truncation has been incorporated for $\bm{H}_{a}^{\ell}$ through the zero submatrix corresponding to $\bm{H}_{a}^{\ell+1}$.

We further split $\bm{H}_{a}^{\ell}$ and $\bm{H}_{*}^{\ell+1}$ as V-splines ($\bm{H}_{v,a}^{\ell}$, $\bm{H}_{v,*}^{\ell+1}$) and F-splines ($\bm{H}_{f\!,a}^{\ell}$, $\bm{H}_{f\!,*}^{\ell+1}$). As a result, Eq. (30) becomes

$$\begin{bmatrix}\bm{H}_{v,a}^{\ell}\\ \bm{H}_{f\!,a}^{\ell}\\ \bm{H}_{v,a}^{\ell+1}\\ \bm{H}_{f\!,a}^{\ell+1}\\ \end{bmatrix}=\begin{bmatrix}\bm{0}&\bm{0}&\bm{C}_{11}&\bm{C}_{12}\\ \bm{0}&\bm{0}&\bm{C}_{21}&\bm{C}_{22}\\ \bm{I}&\bm{0}&\bm{0}&\bm{0}\\ \bm{0}&\bm{I}&\bm{0}&\bm{0}\\ \end{bmatrix}\begin{bmatrix}\bm{H}_{v,a}^{\ell+1}\\ \bm{H}_{f\!,a}^{\ell+1}\\ \bm{H}_{v,p}^{\ell+1}\\ \bm{H}_{f\!,p}^{\ell+1}\\ \end{bmatrix},$$

(31)

where $\bm{C}_{21}=\bm{0}$ because an F-spline can only have F-splines as its children. Moreover, $\bm{H}_{v,*}^{\ell+1}$ can be evaluated according to Eq. (29), that is, $\bm{H}_{v,*}^{\ell+1}=\bm{T}_{*}^{\ell+1}\bm{B}_{f\!,p}^{\ell+1}$, where $\bm{B}_{f\!,p}^{\ell+1}$ is a vector of D-passive F-splines that do not belong to $\mathcal{B}^{\ell+1}$. Substituting this into Eq. (31), we have

$$\bm{H}_{a}=\begin{bmatrix}\bm{H}_{v,a}^{\ell}\\ \bm{H}_{f\!,a}^{\ell}\\ \bm{H}_{v,a}^{\ell+1}\\ \bm{H}_{f\!,a}^{\ell+1}\\ \end{bmatrix}=\begin{bmatrix}\bm{0}&\bm{0}&\bm{C}_{11}&\bm{C}_{12}\\ \bm{0}&\bm{0}&\bm{0}&\bm{C}_{22}\\ \bm{I}&\bm{0}&\bm{0}&\bm{0}\\ \bm{0}&\bm{I}&\bm{0}&\bm{0}\\ \end{bmatrix}\begin{bmatrix}\bm{0}&\bm{0}&\bm{T}_{a}^{\ell+1}\\ \bm{I}&\bm{0}&\bm{0}\\ \bm{0}&\bm{0}&\bm{T}_{p}^{\ell+1}\\ \bm{0}&\bm{I}&\bm{0}\\ \end{bmatrix}\begin{bmatrix}\bm{H}_{f\!,a}^{\ell+1}\\ \bm{H}_{f\!,p}^{\ell+1}\\ \bm{B}_{f\!,p}^{\ell+1}\\ \end{bmatrix}=:\bm{M}_{f}\bm{B}_{f}^{\ell+1},$$

(32)

where $\bm{B}_{f}^{\ell+1}$ contains the full set of F-splines (including both D-active ones $\bm{B}_{f\!,a}^{\ell+1}:=[\bm{H}_{f\!,a}^{\ell+1,T},\bm{H}_{f\!,p}^{\ell+1,T}]^{T}$, and D-passive ones $\bm{B}_{f\!,p}^{\ell+1}$) defined on the transition element at Level $\ell+1$. Eq. (32) implies that THU-splines $\bm{H}_{a}$ defined on a transition element can be evaluated in this unified manner as linear combinations of F-splines. The key is to obtain elementwise $\bm{M}_{f}$. With the two kinds of truncation, it can be shown that each column sum of $\bm{M}_{f}$ equals to one. Therefore, THU-splines $\bm{H}_{a}$ form a partition of unity. While we postpone the related proofs to a forthcoming paper, we will show numerical evidence in Section 5.

When $\Omega_{e}^{\ell+1}$ is a non-transition enriched element (i.e., type-E), only F-splines $\bm{H}_{f\!,*}^{\ell+1}$ are involved. In other words, we have $\bm{H}_{v,*}^{\ell+1}=\bm{0}$ in Eq. (31). In contrast, when $\Omega_{e}^{\ell+1}$ is a type-R element, only $C^{2}$ B-splines $\bm{H}_{v,*}^{\ell+1}$ are involved and they not truncated. Moreover, $\bm{H}_{f\!,*}^{\ell+1}=\bm{0}$ as no F-splines are involved. In fact, this latter case coincides with standard THB-splines.

We discuss how THU-splines can achieve a consistent refinement in this section. We approach it in a constructive manner to show how the two kinds of truncation work together to retain this property.

We assume that only two levels (i.e., $\ell$ and $\ell+1$) are involved to provide insights of construction while keeping notations simple. Globally refined meshes ($\mathcal{M}^{\ell}$ and $\mathcal{M}^{\ell+1}$) and and their bases ($\mathcal{B}^{\ell}$ and $\mathcal{B}^{\ell+1}$) are in place to help discussion. Also recall the following notations for index sets at Level $L\in\{\ell,\ell+1\}$:

• •

  $\mathcal{F}_{a}^{L}$ and $\mathcal{F}_{p}^{L}$: D-active and D-passive F-points/splines, respectively;

• •

  $\mathcal{V}_{a}^{L}$: D-active V-points/splines; and

• •

  $\mathcal{A}^{L}$: H-active points/splines (which can be F- or V-points/splines).

In what follows, we start with Level $\ell+1$ alone and then add back Level $\ell$. Without loss of generality, we study a transition element at Level $\ell+1$, $\Omega_{e}^{\ell+1}$, which can be enriched or non-enriched. Its surface patch can be written as follows with a full set of Level-($\ell+1$) F-splines,

$$S_{e}^{\ell+1}=\sum_{i\in\mathcal{F}_{a}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{i\in\mathcal{F}_{p}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1},$$

(33)

which is essentially the same as Eq. (21). Next, we express $S_{e}^{\ell+1}$ in terms of D-active points/splines only. A D-passive F-point is obtained through its neighboring D-active V-points at the same level, so for $i\in\mathcal{F}_{p}^{\ell+1}$, we have

$$Q_{i}^{\ell+1}=\sum_{j\in\mathcal{V}_{a}^{\ell+1}}d_{ji}^{\ell+1}P_{j}^{\ell+1},\quad\sum_{j\in\mathcal{V}_{a}^{\ell+1}}d_{ji}^{\ell+1}=1,$$

(34)

which is in fact the same as Eq. (13), with $d_{ji}^{\ell+1}\in\{4/9,2/9,1/9,0\}$. In fact, there are only four of $d_{ji}^{\ell+1}$ that are nonzero. Substituting Eq. (34) to Eq. (33), we have

	$\displaystyle S_{e}^{\ell+1}$	$\displaystyle=\sum_{i\in\mathcal{F}_{a}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{i\in\mathcal{F}_{p}^{\ell+1}}B_{f\!,i}^{\ell+1}\sum_{j\in\mathcal{V}_{a}^{\ell+1}}d_{ji}^{\ell+1}P_{j}^{\ell+1}$		(35)
		$\displaystyle=\sum_{i\in\mathcal{F}_{a}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{j\in\mathcal{V}_{a}^{\ell+1}}P_{j}^{\ell+1}\sum_{i\in\mathcal{F}_{p}^{\ell+1}}d_{ji}^{\ell+1}B_{f\!,i}^{\ell+1}$
		$\displaystyle=\sum_{i\in\mathcal{F}_{a}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{j\in\mathcal{V}_{a}^{\ell+1}}P_{j}^{\ell+1}\mathrm{trun}_{D}B_{v,j}^{\ell+1},$

where note that we have used the definition of the same level in the last equation; see Eq. (26). Eq. (35) is the local D-patch representation and it is indeed equivalent to Eq. (27).

Now we add to Eq. (35) the refinement relation between Levels $\ell+1$ and $\ell$. A Level-($\ell+1$) D-active V-point $P_{j}^{\ell+1}$ ($j\in\mathcal{V}_{a}^{\ell+1}$) is computed with the Level-$\ell$ D-active V-points,

$$P_{j}^{\ell+1}=\sum_{j\in\mathcal{V}_{a}^{\ell}}h_{ji}^{\ell+1}P_{j}^{\ell},\quad\sum_{j\in\mathcal{V}_{a}^{\ell}}h_{ji}^{\ell+1}=1,$$

(36)

where $h_{ji}^{\ell+1}$ are the same coefficients as Eq. (7) and only a small number of them are nonzero. On the other hand, a Level-($\ell+1$) D-active F-point is computed from Level-$\ell$ F-points (both D-active and D-passive), so for $i\in\mathcal{F}_{a}^{\ell+1}$, we have

$$Q_{i}^{\ell+1}=\sum_{j\in\mathcal{F}_{a}^{\ell}}h_{ji}^{\ell+1}Q_{j}^{\ell}+\sum_{j\in\mathcal{F}_{p}^{\ell}}h_{ji}^{\ell+1}Q_{j}^{\ell},\quad\sum_{j\in\mathcal{F}_{a}^{\ell}\cup\mathcal{F}_{p}^{\ell}}h_{ji}^{\ell+1}=1,$$

(37)

which is merely applying the B-spline knot insertion to doubly repeated knot vectors. When $j\in\mathcal{F}_{p}^{\ell}$, $Q_{j}^{\ell}$ is computed according to Eq. (34) (replacing $\ell+1$ with $\ell$).

We further divide $\mathcal{F}_{a}^{L}$ and $\mathcal{V}_{a}^{L}$ ($L\in\{\ell,\ell+1\}$) into H-active and H-passive subsets, leading to $\mathcal{A}_{f}^{L}:=\mathcal{F}_{a}^{L}\cap\mathcal{A}^{L}$, $\mathcal{A}_{v}^{L}:=\mathcal{V}_{a}^{L}\cap\mathcal{A}^{L}$, $\mathcal{N}_{f}^{L}:=\mathcal{F}_{a}^{L}\backslash\mathcal{A}^{L}$, and $\mathcal{N}_{v}^{L}:=\mathcal{V}_{a}^{L}\backslash\mathcal{A}^{L}$. Next, starting from Eq. (35), we have

	$\displaystyle S_{e}^{\ell+1}$	$\displaystyle=\sum_{i\in\mathcal{F}_{a}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{j\in\mathcal{V}_{a}^{\ell+1}}P_{j}^{\ell+1}\mathrm{trun}_{D}B_{v,j}^{\ell+1}$		(38)
		$\displaystyle=\sum_{i\in\mathcal{A}_{f}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{i\in\mathcal{N}_{f}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}$
		$\displaystyle+\sum_{j\in\mathcal{A}_{v}^{\ell+1}}P_{j}^{\ell+1}\mathrm{trun}_{D}B_{v,j}^{\ell+1}+\sum_{j\in\mathcal{N}_{v}^{\ell+1}}P_{j}^{\ell+1}\mathrm{trun}_{D}B_{v,j}^{\ell+1}.$

We keep the terms corresponding to H-active subsets unchanged but further expand H-passive ones using Eqs. (36, 37),

	$\displaystyle S_{e,\mathcal{N}}^{\ell+1}$	$\displaystyle:=\sum_{i\in\mathcal{N}_{f}^{\ell+1}}B_{f\!,i}^{\ell+1}Q_{i}^{\ell+1}+\sum_{j\in\mathcal{N}_{v}^{\ell+1}}P_{j}^{\ell+1}\mathrm{trun}_{D}B_{v,j}^{\ell+1}$		(39)
		$\displaystyle=\sum_{i\in\mathcal{N}_{f}^{\ell+1}}B_{f\!,i}^{\ell+1}\left(\sum_{j\in\mathcal{F}_{a}^{\ell}}h_{ji}^{\ell+1}Q_{j}^{\ell}+\sum_{j\in\mathcal{F}_{p}^{\ell}}h_{ji}^{\ell+1}Q_{j}^{\ell}\right)+\sum_{i\in\mathcal{N}_{v}^{\ell+1}}\left(\sum_{j\in\mathcal{V}_{a}^{\ell}}h_{ji}^{\ell+1}P_{j}^{\ell}\right)\mathrm{trun}_{D}B_{v,i}^{\ell+1}$
		$\displaystyle=\sum_{i\in\mathcal{N}_{f}^{\ell+1}}B_{f\!,i}^{\ell+1}\left(\sum_{j\in\mathcal{F}_{a}^{\ell}}h_{ji}^{\ell+1}Q_{j}^{\ell}+\sum_{j\in\mathcal{F}_{p}^{\ell}}h_{ji}^{\ell+1}\sum_{k\in\mathcal{V}_{a}^{\ell}}d_{kj}^{\ell+1}P_{k}^{\ell}\right)$
		$\displaystyle+\sum_{i\in\mathcal{N}_{v}^{\ell+1}}\left(\sum_{j\in\mathcal{V}_{a}^{\ell}}h_{ji}^{\ell+1}P_{j}^{\ell}\right)\mathrm{trun}_{D}B_{v,i}^{\ell+1}.$