Unstructured Transfinite Elements Including T-Splines

1. Introduction

Computational methods have undergone several stages of development and continue to evolve. The progression began with global approximation techniques, such as those introduced by Rayleigh [1], Ritz [2], and the Bubnov-Galerkin method (1913-1915), eventually giving way to more localized approaches like the finite element method [3,4]. Among these developments, transfinite interpolation occupies a particularly important role [5,6]. Initially developed to interpolate the geometry of solid objects and surfaces [7], it soon became apparent that the method could also be used to interpolate scalar and even vector fields within three-dimensional domains. Notably, transfinite interpolation allows the construction of finite elements with differing numbers of nodes along opposing edges or surfaces–without requiring transitional elements. This flexibility extends to the merging of dissimilar domains in two or three dimensions.

An effort to integrate ideas from computer-aided geometric design (CAGD) into analysis methods—such as the finite element, boundary element, and collocation methods—was initiated at the National Technical University of Athens (NTUA) in the summer of 1984, within the research group to which the present author belonged; the first publication appeared later in 1989 [8]. A comprehensive monograph, compiled from approximately 60 journal and conference papers produced by this group, along with an extensive literature review until 2018 (here omitted, for the sake of brevity), is presented in reference [9]. In that monograph, the term macroelement is used to describe an entire patch, whether a surface or a volumetric block. In this framework, although B-splines required domain decomposition at individual breakpoints [8], they were initially regarded more as integration cells than as B-spline elements in the modern sense [10]. Concurrently, similar ideas were being explored by other research groups, as documented in [11,12,13,14,15].

Classical literature on transfinite interpolation (see Refs. [5,6] and references therein) primarily employs linear blending functions, typically Lagrange polynomials. The univariate trial functions defined along the mesh-lines (also referred to as `stations’) are commonly chosen as piecewise linear, piecewise Hermite, or cubic splines. However, in the case of cubic splines, specific implementation details are generally omitted or not fully addressed.

While transfinite interpolation works perfectly well in conjunction with Lagrange polynomials, there is not much experience when it is combined with Bernstein polynomials and B-splines. This is because the former are cardinal polynomials of $[0,1]$-type, in the sense they are interpolatory (take the value of unity) along the stations, whereas the latter are not (take a value less than unity). Nevertheless, both sets hold the partition of unity property. This fact has recently sustained the conjecture that transfinite interpolation can be extended in such a way that blending functions can be Bernstein polynomials or B-splines and at the same time they can also be trial functions which interpolate the primary variable U along the stations of the curvilinear patch [16].

Figure 1 presents a collection of multiple classes of transfinite finite elements, starting from the older Coons element (studied in detail in [8,9]) and ending to a multi-layer element. The well-known case of the T-mesh element [17,18,19] —which is a variation of cases (d,e)– is still missing in this collection, but will be extensively discussed later in Section 6.

Recent studies have shown that a direct substitution of Lagrange polynomials with Bernstein polynomials in the final expressions of the shape functions for all types of transfinite elements is feasible, and sometimes yields equivalent results for both polynomial sets. To date, this direct replacement–accompanied by a pointwise agreement between the two polynomial families–has been verified for the following four classes of transfinite elements [20,21]:

• Coons elements (Figure 1a).
• Tensor product elements (Figure 1b).
• Classical transfinite elements of structured pattern (Figure 1c).
• Deformed tensor-product elements (not shown; similar to Figure 1b, but with boundary nodes arbitrarily positioned relative to the orthogonal projection of the internal nodes).

While the primary variable U is a straightforward physical quantity to interpret, the associated discrete generalized coefficients often cause confusion—particularly among engineers rather than applied mathematicians. In general, however, the coefficients $a_i$ form a set of dual quantities corresponding to the nodal values $U_i$, and the two sets are interconnected through a linear relationship.

First of all, we must report that the expression of tensor-product Bernstein polynomials is not an axiom but a direct sequence of the pre-existed tensor-product of Lagrange polynomials. Clearly, the existing linear relationship between univariate Lagrange and Bernstein polynomials –which share the same monomials– is adequate to offer a mathematical proof for the equivalence between tensor-product Lagrange and tensor-product (non-rational) Bernstein polynomials [9].

Moreover, having established the above explanation for the tensor product of Bernstein polynomials, the interpretation of the tensor-product B-spline follows as a straightforward extension. However, while the rationale behind the tensor-product B-spline is relatively easy to grasp, the same cannot be said for transfinite interpolation, which—even to this day—remains a powerful and versatile tool for constructing novel, large-scale finite elements, as will be demonstrated in the following sections.

The difficulty and delay in advancing the transfinite interpolation method stem from the following issue. Although replacing Lagrange polynomials with Bernstein polynomials is generally straightforward, the challenge lies in substituting the univariate functions $U(\overline{\xi },\eta )$ along each section (station) orthogonal to the $\xi$-axis at point $\overline{\xi }$. This is due to the fact that all the generalized coefficients $a_i$ associated with the two opposite sides (i.e., $AB$ and $DC$) of the quadrilateral patch $ABCD$ (Figure 2b), which are perpendicular to the station under consideration (thus containing its ends), contribute to the determination of function $U(\overline{\xi },\eta )$. Therefore, the classical cardinal blending functions of Lagrange type cannot be directly combined with standard B-splines along the stations, simply because the latter are of global character, that is, it is not possible to interpolate the function $U(\overline{\xi },\eta )$ through a self-contained expression (trial functions) along a $\overline{\xi }$-station. Addressing this incompatibility is the central objective of the present study.

Aiming to address the above issue, a recent preliminary study on transfinite elements has indicated that the classical blending functions –which usually are cardinal functions of [1,0]-type– can be substituted by non-cardinal functions such as Bernstein polynomials and B-splines [16]. The beginning of this inspiration is the 2D analogue of the 1D interpolation using Lagrange polynomials from one side and any other IGA-based functions such as Bernstein polynomials, B-splines, and NURBS. The lowest level of this concept is to think about the quadratic interpolation, which can be implemented either using Lagrange polynomials: $U\left(\xi \right)=L_{1,2}\left(\xi \right)U_1+L_{2,2}\left(\xi \right)U_2+L_{3,2}\left(\xi \right)U_2$ or using Bernstein polynomials: $U\left(\xi \right)=B_{1,2}\left(\xi \right)a_1+B_{2,2}\left(\xi \right)a_2+B_{3,2}\left(\xi \right)a_2$. The issue is as follows: It is quite reasonable that the aforementioned three Lagrange polynomials may serve the role of blending (interpolation) functions in the $\xi$-direction since they fullfil the delta-Kronecker property ($L_i\left({\xi }_j\right)={\delta }_{ij}$) and hold the partition the unity property. Therefore, they can accurately represent a constant, a linear, and a quadratic univariate function. Similarly, if the generalized coefficients ($a_1,a_2,a_3$) are properly chosen each time, the Bernstein polynomials may again accurately represent the same functions (constant, linear, and quadratic).

Clearly, when the blending functions $E_i\left(\xi \right)$ are Bernstein polynomials (or B-splines), the variation in the $\eta$-direction should be performed in terms of a (supposed existing) univariate function $a\left(\eta \right)$, which for $\eta =0$ (i.e., on edge $AB$) coincides with the usual generalized coefficients ($a_1,a_2,a_3$) at points (${\xi }_1,{\xi }_2,{\xi }_3$), respectively. The same holds for a triplet (${\xi }_1^{\prime },{\xi }_2^{\prime },{\xi }_3^{\prime }$) on the opposite edge $DC$ at $\eta =1$.

The critical point is how to describe the abovementioned function $a({\overline{\xi }}_i,\eta )$, which is perpendicularly oriented at point ${\overline{\xi }}_i$. Within the IGA framework, the most easy way is to employ a set of B-spline functions extended in the interval $\eta \in [0,1]$. To this purpose, a safe way is to use clamped splines (for degree $p=3$: $\Xi =[0,0,0,0,\dots ,1,1,1,1]$), whose the end values will be the boundary variable a (and not U). In this way, the previously mentioned difficulty in handling $a\left(\eta \right)$ is removed.

The paper is divided into two main parts. The first part focuses on the construction and performance of simplified elements, both classical and newly introduced, including those composed of layers oriented in a single direction. This serves as a pilot study to validate the approach of treating the discrete coefficients as a bivariate function. The second part addresses the well-known T-spline elements.

The present paper is structured as follows: Section 2 introduces the general formulation of transfinite interpolation, employing either Lagrange or Bernstein polynomials, as well as B-splines, for the representation of blending and trial functions. Section 3 describes the construction of classical transfinite elements. Section 4 presents the development of transfinite elements with nodes arranged in parallel layers. Section 5 focuses on elements that combine tensor-product interior nodes with nonuniform boundary node placement. Section 6 addresses T-spline elements of increasing complexity. Section 7 discusses software-related aspects. Section 8 provides numerical results for all models, Section 9 briefly refers to possible implementation of the collocation method as an alternative, and Section 10 offers a detailed discussion. The paper is supplemented by seven appendices.

2. Transfinite Interpolation Using B-Splines

2.1. Terminology

In classical computer-aided geometric design (CAGD) [22], the term blending function typically refers to univariate interpolation functions used in a specific parametric direction. For instance, consider a Coons patch defined over the quadrilateral $ABCD$, with the origin at vertex A and the $\xi$-axis oriented along edge $AB$ (Figure 2b) [2]. The blending functions $E_1\left(\xi \right)$ and $E_2\left(\xi \right)$ interpolate between the boundary functions $U_{AD}(\eta$) and $U_{BC}(\eta$) in the $\xi$-direction.

This usage of the term differs from its meaning in T-spline literature, where blending functions are often synonymous with basis functions [18]. In this paper, we adopt the classical interpretation: the term blending function will refer specifically to the interpolation functions used in Coons or Gordon patch construction, as adopted in Refs. [2,22,2,22]. Throughout the present paper, the term basis functions will be used uniformly for both B-splines/NURBS and T-splines.

A second key concept in classical transfinite interpolation theory–as well as in finite element methodology, as emphasized in Zienkiewicz’s seminal work–is that of trial functions [2]. In the present study, this term refers to univariate functions defined along boundary edges, or more generally, along horizontal and vertical mesh-lines (stations) within the domain.

A third issue concerns the use of the terms transfinite element and macroelement. In this paper, these terms are used synonymously to refer to the parametric patch defined over the domain $0\le \xi ,\eta \le 1$. Regardless of the choice of blending or trial functions, the resulting basis functions share the same closed-form expression. The main difference lies in the interpolation approach: Lagrange and Bernstein polynomials perform global interpolation over the entire patch, whereas B-splines use piecewise polynomial interpolation between breakpoints, which define the integration cells. Naturally, in the B-spline formulation, not all degrees of freedom (DOFs) are active within every integration cell. However, to maintain a uniform treatment across all three cases (Lagrange, Bernstein, and B-splines), we do not emphasize this distinction.

2.2. State of the Art

The general expression of transfinite interpolation for 2D patches is given in [2]:

$$\begin{array}{cc}\hfill U(\xi ,\eta )& =P_{\xi }\oplus P_{\eta }\hfill \\ \hfill & =P_{\xi }\left\{U\right\}+P_{\eta }\left\{U\right\}-P_{\xi \eta }\left\{U\right\},\hfill \end{array}$$

(1)

where $P_{\xi }\left\{U\right\}$ and $P_{\eta }\left\{U\right\}$ are projectors of the quantity U along the parametric axes $\xi$ and $\eta$, respectively, and $P_{\xi \eta }\left\{U\right\}$ is the corrective projector, defined as the tensor product of nodal values at the intersections of the coordinate stations.

It is well known that, for vertical stations passing through $\xi ={\xi }_1=0,\dots ,{\xi }_{m+1}=1$, and for horizontal stations located at $\eta ={\eta }_1=0,\dots ,{\eta }_{n+1}=1$, we construct univariate blending functions $E_i\left(\xi \right)$ and $E_j\left(\eta \right)$ of degrees m and n, respectively (assuming Lagrange polynomials are used). Based on the univariate functions $U({\xi }_i,\eta )$ and $U(\xi ,{\eta }_j)$ defined along the vertical and horizontal stations, respectively, we obtain:

$$\begin{array}{cc}\hfill P_{\xi }\left\{U\right\}& =\sum _{i=1}^{m+1}{E}_i\left(\xi \right)U({\xi }_i,\eta ),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill P_{\eta }\left\{U\right\}& =\sum _{j=1}^{n+1}{E}_j\left(\eta \right)U(\xi ,{\eta }_j),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill P_{\xi \eta }\left\{U\right\}& =\sum _{j=1}^{n+1}\sum _{i=1}^{m+1}{E}_i\left(\xi \right)E_j\left(\eta \right)U({\xi }_i,{\eta }_j).\hfill \end{array}$$

(4)

On the other hand, along each station, the variable $U(\xi ,\eta )$ is interpolated using trial functions ${\tilde{B}}_i\left(s\right)$, where the variable s denotes either of the parametric directions $\xi$ or $\eta$, as follows:

$$\begin{array}{cc}\hfill \mathrm{Horizontal}\mathrm{station}\text{:}U(\xi ,\overline{\eta })& =\sum _{i=1}^{\tilde{m}+1}{\tilde{B}}_i\left(\xi \right)U({\xi }_i,\overline{\eta }),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \mathrm{Vertical}\mathrm{station}\text{:}U(\overline{\xi },\eta )& =\sum _{j=1}^{\tilde{n}+1}{\tilde{B}}_j\left(\eta \right)U(\overline{\xi },{\eta }_j).\hfill \end{array}$$

(6)

Here, the overbar denotes a prescribed value of a given parameter ($\xi$ or $\eta$) along the corresponding station—vertical for $\overline{\xi }$ or horizontal for $\overline{\eta }$. In general, the number of nodes ($\tilde{m}+1$ or $\tilde{n}+1$) along a station may differ—either smaller or greater—from the number of stations in the same direction. The case of equality (e.g., in a tensor-product configuration) is also included.

Henceforth, for brevity, the notation $\left\{U\right\}$ following each projector is suppressed; for example, $P_{\xi }\left\{U\right\}$ is written simply as $P_{\xi }$, and similarly for the others.

2.3. Extension of the Projectors Using Bernstein Polynomials

The classical interpretation of the well-known projectors ($P_{\xi }$, $P_{\eta }$, and $P_{\xi \eta }$) is founded on their application to the primary variable $U(\xi ,\eta )$, in accordance with the use of Lagrange polynomials as blending functions $E_i\left(\xi \right)$ and $E_j\left(\eta \right)$, as presented in Ref. [2]. It is well established that Coons interpolation represents a direct extension of one-dimensional linear interpolation to two-dimensional domains [9]. Similarly, classical transfinite (Gordon) interpolation generalizes one-dimensional polynomial interpolation to bivariate patches. Since one-dimensional polynomial interpolation can be formulated using either Lagrange or Bernstein polynomials [23], it follows that the classical projectors may likewise be expressed in terms of Bernstein polynomials.

The above claim will be demonstrated—in full detail for the first time—for one projector, namely $P_{\xi }$; the same holds for the remaining ones. Let us consider three points $P_1\left({\xi }_1\right)$, $P_2\left({\xi }_2\right)$, and $P_3\left({\xi }_3\right)$ on a parametric curve $C\left(\xi \right)$ and the associated values $U_1$, $U_2$, $U_3$ of a univariate function $U\left(\xi \right)$ (Figure 2a). These can be interpolated in terms of the nodal values ($U_1$, $U_2$, $U_3$), using the quadratic Lagrange polynomials $E_1\left(\xi \right)=L_{1,2}\left(\xi \right)$, $E_2\left(\xi \right)=L_{2,2}\left(\xi \right)$, $E_3\left(\xi \right)=L_{3,2}\left(\xi \right)$:

$$U\left(\xi \right)=E_1\left(\xi \right)U_1+E_2\left(\xi \right)U_2+E_3\left(\xi \right)U_3.$$

(7)

Alternatively, the same points can be interpolated in terms of the generalized coefficients ($a_1$, $a_2$, $a_3$) using the quadratic Bernstein polynomials ${\tilde{E}}_1\left(\xi \right)=B_{1,2}\left(\xi \right)$, ${\tilde{E}}_2\left(\xi \right)=B_{2,2}\left(\xi \right)$, and ${\tilde{E}}_3\left(\xi \right)=B_{3,2}\left(\xi \right)$:

$$U\left(\xi \right)={\tilde{E}}_1\left(\xi \right)a_1+{\tilde{E}}_2\left(\xi \right)a_2+{\tilde{E}}_3\left(\xi \right)a_3.$$

(8)

The above equations can be extended from one to two dimensions, beginning with the projector $P_{\xi }$, in which the nodal values ($U_1$, $U_2$, $U_3$) are replaced by univariate functions ($U_1\left(\eta \right)$, $U_2\left(\eta \right)$, $U_3\left(\eta \right)$), as illustrated in Figure 2b:

$$P_{\xi }\left\{U\right\}=E_1\left(\xi \right)U_1\left(\eta \right)+E_2\left(\xi \right)U_2\left(\eta \right)+E_3\left(\xi \right)U_3\left(\eta \right).$$

(9)

Similarly, in terms of Bernstein polynomials, the same projector is expressed in the form of a separation of variables, as follows:

$$P_{\xi }\left\{a\right\}={\tilde{E}}_1\left(\xi \right)a_1\left(\eta \right)+{\tilde{E}}_2\left(\xi \right)a_2\left(\eta \right)+{\tilde{E}}_3\left(\xi \right)a_3\left(\eta \right).$$

(10)

Of course, a significant distinction exists between the two formulations presented above. The Lagrange interpolation (Eq. (9)) is a local scheme applied to univariate functions that directly represent the primary variable $U\left(\eta \right)$ at discrete stations located at $\xi =0,{\displaystyle \frac{1}{2}},1$. In contrast, the Bernstein interpolation (Equation (10)) is a global approach, involving univariate functions $a\left(\eta \right)$ constructed from generalized coefficients.

This implies that while the interpolation of the primary variable U is straightforward and intuitive, the corresponding interpolation of the coefficient a is less direct. Nevertheless, the coefficient a embodies a definitive quantity that, in essence, is not fundamentally different from the variable U. For instance, in the case of quadratic interpolation, we observe:

$$a_1=U_1,a_2=-{\displaystyle \frac{1}{2}}{U}_1+2U_2+{\displaystyle \frac{1}{2}}{U}_3,a_3=U_3.$$

(11)

Therefore, if—for instance—the extreme values are equal (i.e., $U_1=U_3$), the coefficient $a_2$ corresponds to $2U_2$. This implies that a duplication of the associated basis function can be employed to ensure that $a_2$ effectively represents the actual variable $U_2$.

Numerical investigations have confirmed that this concept—promoting the use of Bernstein polynomials, and more generally B-splines, as blending functions—yields robust performance, provided that the remaining two projectors ($P_{\xi },P_{\eta }$) are incorporated in a consistent manner [9,16].

2.4. B-Splines Projectors

A comprehensive investigation into the use of B-splines as blending functions was recently presented in Ref. [16], although the treatment therein was primarily speculative. In contrast, the methodology introduced in this paper adopts a more rigorous framework and may be interpreted in analogy with a function of separated variables, namely $U(\xi ,\eta )=X\left(\xi \right)Y\left(\eta \right)$, where the role of $Y\left(\eta \right)$ is assumed by the function $a\left(\eta \right)$, which is intrinsically linked to the generalized coordinates.

In general, the number of mesh lines (stations) along each parametric direction determines the length of the corresponding knot vector used for the blending functions. Likewise, the number of control points defined at each station governs the length of the associated knot vector for the trial functions. Further clarification and illustrative examples will be provided in the subsequent sections.

3. Construction of Classical B-Spline Transfinite Elements

A representative transfinite element $ABCD$, comprising 113 control points, is depicted in Figure 3d. It can be observed that the element is structured with four horizontal and five vertical stations (including the boundary edges), upon which the control points are distributed. The horizontal set of stations is oriented perpendicular to the $\eta$-axis and positioned at $\eta =0,{\displaystyle \frac{1}{3}},{\displaystyle \frac{2}{3}},1$, while the vertical set is oriented perpendicular to the $\xi$-axis and located at $\xi =0,{\displaystyle \frac{1}{4}},{\displaystyle \frac{2}{4}},{\displaystyle \frac{3}{4}},1$.

Based on the aforementioned station configuration, the blending functions in the $\xi$-direction may be chosen as Lagrange or Bernstein polynomials of degree 4, corresponding to five station values. Similarly, in the $\eta$-direction, the blending functions may be Lagrange or Bernstein polynomials of degree 3, associated with four station values.

Alternatively, when employing cubic B-splines, the $\xi$-direction can be represented using the knot vector

$${\Xi }_{\mathrm{blend}}=[0,0,0,0,\frac{1}{2},1,1,1,1],$$

which corresponds to five distinct station values and ensures $C^2$ continuity. In the $\eta$-direction, the blending functions may be defined by the knot vector

$$H_{\mathrm{blend}}=[0,0,0,0,1,1,1,1],$$

which corresponds to four station values and yields Bernstein polynomials of degree 3.

In all cases, the construction of the knot vectors for the blending functions must guarantee that the number of control points aligns with the number of stations, thereby ensuring proper interpolation and continuity across the transfinite element.

Under these conditions, the three projectors defined by Eqs. (2) through () are now reformulated in terms of B-spline-based blending and trial functions, where the trial functions are associated with the generalized coefficients $a_i$.

Upon further elaboration, it is found that the Boolean sum in Equation (1) yields a total of 113 basis functions. These functions can be categorized into two distinct groups, as follows:

(1)

Nodes at section intersections (including corner nodes A, B, C, and D), for example:

$${\psi }_{69}(\xi ,\eta )=E_2\left(\xi \right)N_{9,3}\left(\eta \right)+E_3\left(\eta \right)N_{5,3}\left(\xi \right)-E_2\left(\xi \right)E_3\left(\eta \right).$$

(12)

In this expression, $E_2\left(\xi \right)$ denotes the second blending function in the $\xi$-direction, corresponding to the second vertical station where node 69 is located. Likewise, $N_{9,3}\left(\eta \right)$ refers to the 9th B-spline basis function in the $\eta$-direction, as node 69 is the 9th node from the bottom. The same logic applies to the remaining terms.

(2)

Intermediate nodes on single stations, for example:

$${\psi }_{71}(\xi ,\eta )=E_3\left(\eta \right)N_{7,3}\left(\xi \right).$$

(13)

Here, $E_3\left(\eta \right)$ is the third blending function in the $\eta$-direction, indicating that node 71 lies on the third horizontal station. Similarly, $N_{7,3}\left(\xi \right)$ represents the 7th B-spline basis function in the $\xi$-direction, as node 71 is the 7th node from the left boundary.

By classifying the 113 degrees of freedom into the two aforementioned categories, the complete set of basis functions can be systematically constructed. Specifically, 20 of these functions—corresponding to the black-filled circles in Figure 3d—are composed of three terms, as exemplified by Equation (12), and reflect contributions from all three projectors. The remaining 93 functions—associated with the white-filled circles—are expressed as simple tensor products, as illustrated in Equation (13).

With respect to the blending and trial functions, the formulations presented above are of general applicability. In the following sections, several alternative modeling approaches will be examined and compared.

3.1. Model-1: Lagrange Polynomials

As a representative example, we consider the following configuration (Figure 3d):

(1)

Lagrange polynomials of degree $p_{\xi }=4$ and $p_{\eta }=3$ are employed as blending functions in the $\xi$- and $\eta$-directions, respectively.

(2)

Lagrange polynomials of degree $p_{\xi }^{\prime }=16$ and $p_{\eta }^{\prime }=12$ are used as trial functions in the $\xi$- and $\eta$-directions, respectively.

For this Model-1 setup, the resulting bivariate global shape functions $\varphi (\xi ,\eta )$ are depicted in Figure 4a. Notably, elevated function values are observed along the vertical edges, indicating localized influence. Furthermore, the constructed basis satisfies the partition of unity property:

$$\sum _{i=1}^{113}{\varphi }_i(\xi ,\eta )=1.$$

(14)

3.2. Model-2: Mixed Scheme

In this model, the blending functions are retained from Model-1, namely Lagrange polynomials of degree $p_{\xi }=4$ and $p_{\eta }=3$ in the $\xi$- and $\eta$-directions, respectively. For the trial functions, however, we adopt the classical natural cubic cardinal B-splines, which are directly associated with the nodal values U. This approach mirrors the behavior of Lagrange polynomials in Model-1, thereby eliminating the need for generalized coefficients $a_i$.

This formulation is particularly well-suited to uniform nodal distributions with single knots and facilitates direct comparison with Model-1. It is worth noting that this type of spline interpolation—natural cubic B-splines—was widely used during the mid-1980s (see [24], among others).

In brief, the knot vector for the trial functions in the $\xi$-direction consists of 23 entries:

$$\Xi =[0,0,0,0,\frac{1}{16},\dots ,\frac{15}{16},1,1,1,1],$$

which generates 19 control points. From these, 17 univariate cardinal trial functions are constructed, corresponding to the 17 nodal points along the horizontal stations.

Based on the formulations in Eqs. (12) and (13), the graphical representation of the resulting 113 basis functions is shown in Figure 4b. These basis functions have been verified to satisfy the partition of unity property:

$$\sum _{i=1}^{113}{\psi }_i(\xi ,\eta )=1,$$

(15)

and exhibit unity as their upper bound.

Further details regarding the construction of the cardinal univariate set of natural cubic B-splines can be found in Ref. [9], and relevant implementation code is provided in Appendix A.

3.3. Model-3: Cubic B-Splines

This model closely aligns with the current state-of-the-art practices in Isogeometric Analysis (IGA), where classical B-splines (de Boor formulation) or NURBS are commonly employed. In this framework, B-splines are utilized for constructing both the blending and trial function sets. Classical clamped splines are adopted, reflecting the fact that for a polynomial degree p, the minimum number of control points is $p+1$ in the case of Bernstein polynomials, whereas for pure B-splines, the condition $n_{\mathrm{ctrl}}>p+1$ must be satisfied.

Within this context, the 113-node transfinite element illustrated in Figure 3d should be interpreted as an index space rather than a physical nodal layout. In Models 1 and 2 (see Sects. Section 3.1 and Section 3.2), each of the 113 nodal points was directly associated with a unique bivariate shape function linked to a nodal value $U_i$, leaving no ambiguity in the trial function definition.

Assuming that the 17 points along each horizontal station are treated as single breakpoints, it is well known that the corresponding number of control points becomes 19. Conversely, if only 17 control points per horizontal station are desired, one must define 15 breakpoints, resulting in 14 knot spans in the $\xi$-direction. The associated knot vector is then given by:

$$\Xi =[0,0,0,0,\frac{1}{14},\frac{2}{14},\dots ,\frac{13}{14},1,1,1,1].$$

(16)

Similarly, to obtain 13 control points in the $\eta$-direction (based on 11 breakpoints and 10 knot spans), the corresponding knot vector is:

$$H=[0,0,0,0,\frac{1}{10},\frac{2}{10},\dots ,\frac{9}{10},1,1,1,1].$$

(17)

The graphical representation of the resulting basis functions is shown in Figure 4c.

4. Construction of Multi-Layer Finite Elements with Nodes Arranged in Parallel Layers

We now examine multi-layer finite elements containing internal nodes, where the number of nodes per station (i.e., per parallel layer) may vary. Such configurations arise when nodes—both internal and boundary—are aligned along unidirectional stations, for instance, distributed horizontally, as depicted in Figure 1e.

4.1. Unidirectional Multi-Layer Transfinite Lagrange and Bernstein Elements

The unidirectional analog of Equation (10), corresponding to the configuration shown in Figure 1e where all five stations are parallel to the horizontal $\xi$-axis, is given by the following projector:

$$\begin{array}{cc}\hfill P_{\eta }\left\{a\right\}& ={\tilde{E}}_1\left(\eta \right)a_1\left(\xi \right)+{\tilde{E}}_2\left(\eta \right)a_2\left(\xi \right)+{\tilde{E}}_3\left(\eta \right)a_3\left(\xi \right)\hfill \\ \hfill & +{\tilde{E}}_4\left(\eta \right)a_4\left(\xi \right)+{\tilde{E}}_5\left(\eta \right)a_5\left(\xi \right).\hfill \end{array}$$

(18)

In this formulation, the approximation function is expressed as $U(\xi ,\eta )=P_{\eta }\left\{a\right\}$. Each component function $a_i\left(\xi \right)$, for $i=1,\dots ,5$, can be represented using clamped B-splines. The corresponding knot vectors are constructed such that the number of control points matches the number of nodes shown in Figure 1e.

It is important to note that the interpolatory nature of these B-splines—exhibiting unit values at the endpoints of each horizontal station (i.e., at $\xi =0$ and $\xi =1$)—does not introduce any complications. This is because each function $a_j\left(\xi \right)$ is multiplied by a non-cardinal blending function ${\tilde{E}}_j\left(\eta \right)$, analogous to the behavior observed in one-dimensional problems.

For instance, if $a\left(\xi \right)$ is a constant function, the clamped B-spline basis will represent it exactly. This observation underscores that the value of U is influenced by all coefficients a within the patch, a property particularly relevant for Bernstein polynomials. However, in the present formulation, the focus is on the coefficient functions $a_i\left(\xi \right)$ rather than the potentially problematic direct interpolation of U along each horizontal station.

A previous report, Ref. [16], proposed the conjecture that transfinite interpolation can be effectively employed using a wide variety of blending functions, provided that the partition of unity condition is satisfied and the function set is complete. Within this framework, not only cardinal Lagrange polynomials but also non-cardinal blending functions—such as Bernstein polynomials and B-splines—have been successfully explored.

Furthermore, in the context of trial functions, any well-established univariate functional basis may be considered a viable candidate for use within numerical analysis modules, including Galerkin and collocation methods. This flexibility allows for the integration of diverse approximation schemes while maintaining consistency with the underlying transfinite interpolation principles.

Let us consider a class of transfinite elements with unidirectional nodes, as shown in Figure 5. Let us focus on the 12-node element (Figure 5a). Replacing Equation (9) by the vertical projection $P_{\eta }=E_1\left(\eta \right)U_1\left(\xi \right)+E_2\left(\eta \right)U_2\left(\xi \right)+E_3\left(\eta \right)U_3\left(\xi \right)$, and then expressing the involved univariate functions $U_1\left(\xi \right),U_2\left(\xi \right),U_3\left(\xi \right)$ in terms of Lagrange polynomials $L_{i,p}\left(\xi \right)$ of variable degree p associated with each horizontal station (layer), we obtain the following set of bivariate global shape functions:

$$\begin{array}{cc}{\varphi }_1(\xi ,\eta )& =L_{1,2}\left(\xi \right)L_{1,2}\left(\eta \right),\\ {\varphi }_2(\xi ,\eta )& =L_{2,2}\left(\xi \right)L_{1,2}\left(\eta \right),\\ {\varphi }_3(\xi ,\eta )& =L_{3,2}\left(\xi \right)L_{1,2}\left(\eta \right),\\ {\varphi }_4(\xi ,\eta )& =L_{1,3}\left(\xi \right)L_{2,2}\left(\eta \right),\\ {\varphi }_5(\xi ,\eta )& =L_{2,3}\left(\xi \right)L_{2,2}\left(\eta \right),\\ {\varphi }_6(\xi ,\eta )& =L_{3,3}\left(\xi \right)L_{2,2}\left(\eta \right),\\ {\varphi }_7(\xi ,\eta )& =L_{3,3}\left(\xi \right)L_{2,2}\left(\eta \right),\\ {\varphi }_8(\xi ,\eta )& =L_{1,4}\left(\xi \right)L_{3,2}\left(\eta \right),\\ {\varphi }_9(\xi ,\eta )& =L_{2,4}\left(\xi \right)L_{3,2}\left(\eta \right),\\ {\varphi }_{10}(\xi ,\eta )& =L_{3,4}\left(\xi \right)L_{3,2}\left(\eta \right),\\ {\varphi }_{11}(\xi ,\eta )& =L_{4,4}\left(\xi \right)L_{3,2}\left(\eta \right),\\ {\varphi }_{12}(\xi ,\eta )& =L_{5,4}\left(\xi \right)L_{3,2}\left(\eta \right)\end{array}$$

(19)

Next, replacing Lagrange polynomials with their Bernstein counterparts, or alternatively applying the projector defined by the analog of Equation (18) (including three terms associated with the three layers in Figure 5a) by expressing the involved univariate functions $a_1\left(\xi \right),a_2\left(\xi \right),a_3\left(\xi \right)$ in terms of Bernstein polynomials, we obtain the following set of non-negative basis functions:

$$\begin{array}{cc}{\psi }_1(\xi ,\eta )& =B_{1,2}\left(\xi \right)B_{1,2}\left(\eta \right),\\ {\psi }_2(\xi ,\eta )& =B_{2,2}\left(\xi \right)B_{1,2}\left(\eta \right),\\ {\psi }_3(\xi ,\eta )& =B_{3,2}\left(\xi \right)B_{1,2}\left(\eta \right),\\ {\psi }_4(\xi ,\eta )& =B_{1,3}\left(\xi \right)B_{2,2}\left(\eta \right),\\ {\psi }_5(\xi ,\eta )& =B_{2,3}\left(\xi \right)B_{2,2}\left(\eta \right),\\ {\psi }_6(\xi ,\eta )& =B_{3,3}\left(\xi \right)B_{2,2}\left(\eta \right),\\ {\psi }_7(\xi ,\eta )& =B_{3,3}\left(\xi \right)B_{2,2}\left(\eta \right),\\ {\psi }_8(\xi ,\eta )& =B_{1,4}\left(\xi \right)B_{3,2}\left(\eta \right),\\ {\psi }_9(\xi ,\eta )& =B_{2,4}\left(\xi \right)B_{3,2}\left(\eta \right),\\ {\psi }_{10}(\xi ,\eta )& =B_{3,4}\left(\xi \right)B_{3,2}\left(\eta \right),\\ {\psi }_{11}(\xi ,\eta )& =B_{4,4}\left(\xi \right)B_{3,2}\left(\eta \right),\\ {\psi }_{12}(\xi ,\eta )& =B_{5,4}\left(\xi \right)B_{3,2}\left(\eta \right).\end{array}$$

(20)

The shape functions ${\varphi }_i(\xi ,\eta )$ based on Lagrange polynomials (Equation (19)) are illustrated in Figure 6a, while the basis functions ${\psi }_i(\xi ,\eta )$ corresponding to Bernstein polynomials (Equation (20)) are shown in Figure 6b. This element is not provided for B-spline interpolation, as it contains too few nodes to support such a representation.

4.2. Unidirectional Transfinite B-Spline Elements

Below we focus on the 18-node element shown in Figure 5b. The characteristic of this element is that—in addition to the four edges—it includes two internal stations which are both oriented in the $\xi$-direction (horizontal stations), and thus (as also happened with the 12-node element) the only non-zero projector among the three is $P_{\eta }$, i.e., $P_{\xi }=P_{\xi \eta }=0$. Since the total number of horizontal stations is 4, there is no cubic B-spline to cover this case (unless we select quadratic B-spline), and therefore we resort to cubic Bernstein polynomials (hybrid scheme), as follows:

(1)

The blending functions in vertical ($\eta$)-direction were taken as cubic Bernstein polynomials, $B_{i,3}\left(\eta \right),i=1,2,3,4$, where $B_{1,3}={(1-\eta )}^3,B_{2,3}=3{(1-\eta )}^2\eta ,B_{3,3}=3(1-\eta ){\eta }^2,B_{4,3}={\eta }^3$.

(2)

The base 1-2-3 was modelled by quadratic Bernstein polynomials: $B_{1,2}={(1-\eta )}^2,B_{2,2}=2(1-\eta )\eta ,B_{3,2}={\eta }^2$.

(3)

The second layer, with nodes 4-5-6-7, was modelled by cubic Bernstein polynomials, $B_{i,3}\left(\xi \right),i=1,2,3,4$.

(4)

The third layer, with nodes 8-9-10-11-12, was modelled as a set of five cubic B-splines, $N_{i,3}^{\eta =2/3}\left(\xi \right),i=1,\dots ,5$, with knot vector $[0,0,0,0,{\displaystyle \frac{1}{2}},1,1,1,1]$.

(5)

The fourth (top) layer, with nodes 13-14-15-16-17-18, was modelled as a set of six cubic B-splines, $N_{i,3}^{\eta =1}\left(\xi \right),i=1,\dots ,6$, with knot vector $[0,0,0,0,{\displaystyle \frac{2}{5}},{\displaystyle \frac{3}{5}},1,1,1,1]$.

Writing a projector similar to that of Equation (18), in conjunction with four parallel layers (stations) at $\eta =0,{\displaystyle \frac{1}{3}},{\displaystyle \frac{2}{3}},1$, we obtain the following basis functions:

$$\begin{array}{cc}\hfill {\psi }_1(\xi ,\eta )& =B_{1,2}\left(\xi \right)B_{1,3}\left(\eta \right),\hfill \\ \hfill {\psi }_2(\xi ,\eta )& =B_{2,2}\left(\xi \right)B_{1,3}\left(\eta \right),\hfill \\ \hfill {\psi }_3(\xi ,\eta )& =B_{3,2}\left(\xi \right)B_{1,3}\left(\eta \right),\hfill \\ \hfill {\psi }_4(\xi ,\eta )& =B_{1,3}\left(\xi \right)B_{2,3}\left(\eta \right),\hfill \\ \hfill {\psi }_5(\xi ,\eta )& =B_{2,3}\left(\xi \right)B_{2,3}\left(\eta \right),\hfill \\ \hfill {\psi }_6(\xi ,\eta )& =B_{3,3}\left(\xi \right)B_{2,3}\left(\eta \right),\hfill \\ \hfill {\psi }_7(\xi ,\eta )& =B_{3,3}\left(\xi \right)B_{2,3}\left(\eta \right),\hfill \\ \hfill {\psi }_8(\xi ,\eta )& =N_{1,3}^{\eta =2/3}\left(\xi \right)B_{3,3}\left(\eta \right),\hfill \\ \hfill {\psi }_9(\xi ,\eta )& =N_{2,3}^{\eta =2/3}\left(\xi \right)B_{3,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{10}(\xi ,\eta )& =N_{3,3}^{\eta =2/3}\left(\xi \right)B_{3,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{11}(\xi ,\eta )& =N_{4,3}^{\eta =2/3}\left(\xi \right)B_{3,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{12}(\xi ,\eta )& =N_{5,3}^{\eta =2/3}\left(\xi \right)B_{3,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{13}(\xi ,\eta )& =N_{1,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{14}(\xi ,\eta )& =N_{2,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{15}(\xi ,\eta )& =N_{3,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{16}(\xi ,\eta )& =N_{4,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{17}(\xi ,\eta )& =N_{5,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right),\hfill \\ \hfill {\psi }_{18}(\xi ,\eta )& =N_{6,3}^{\eta =1}\left(\xi \right)B_{4,3}\left(\eta \right).\hfill \end{array}$$

(21)

Obviously, the same form as Equation (21) holds when the blending and trial functions are all of Lagrange type or all of Bernstein type. For all these three cases, the graphs of the basis functions are shown in Figure 7.

Regarding the estimation of the stiffness matrix of the 18-node transfinite element in Figure 5b, the fully-Lagrange and the fully-Bernstein models require an integration scheme of $6\times 4$ Gauss points in the $\xi$- and $\eta$-directions, without local support. This is because the maximum polynomial degrees per single basis function are $p=5$ and $q=3$, respectively, and thus the multiplication in the integrand will contain terms up to $x^{10}{y}^6$. In contrast, for the same patch, the piecewise cubic B-spline interpolation requires $10\times 3$ integration cells (i.e., 30 B-spline elements, as they have been called in standard IGA [10]) with $4\times 4$ Gauss points per element (i.e., a total of 480 integration points, however with local support).

Note: When the number of parallel sections increases to 5 (and thus the number of nodes becomes 25, as shown in Figure 5c), the 25-node transfinite element may be handled using blending functions of B-spline type, based on the knot vector ${\Xi }_{\mathrm{blend}}=[0,0,0,0,\frac{1}{2},1,1,1,1]$. Moreover, regarding the 32-node element made of 6 horizontal sections (partially non-uniform, as shown in Figure 5d), results and discussion are given in Section 8.1.2.

5. Finite Elements Featuring Tensor-Product Interior Nodes and Nonuniform Boundary Node Placement

This is a class of elements (the third in the present paper) in which the boundary nodes are placed arbitrarily, and the interior is populated using a tensor product grid of $m\times n$ internal nodes. A common rule of thumb is to choose the number of internal points in each direction (m or n) based on the average number of intermediate nodes on opposite edges. This class of elements has previously been formulated using Lagrange polynomials as both blending and trial functions [16,20]. In the present paper, the methodology is extended by introducing B-splines as an alternative choice for both blending and trial functions.

5.1. 27-Node Transfinite Element

As an example, we consider a 27-node transfinite element, in which the four edges are uniformly divided into 4, 3, 6, and 5 segments, respectively, as shown in Figure 8. For this configuration, we can construct at least three different formulations (blending and trial functions) of the transfinite element as follows:

• Lagrange polynomials;
• Bernstein polynomials.
• B-splines

Regarding Lagrange and Bernstein polynomials, the associated polynomial degrees per edge are as follows: $p_{AB}=4$, $p_{BC}=3$, $p_{DC}=6$, and $p_{AD}=5$. The internal nodes form a uniform tensor product pattern of scheme $3\times 3$ (nodes 19 to 27 in Figure 8), which means that the degree of the blending functions is $n_{I,\xi }=n_{I,\eta }=4$ (between edges and internal nodes, four node spans per direction are created). Therefore, this element consists of 18 boundary and 9 internal nodes.

The above 27-node B-spline element is defined by a regular tensor-product arrangement of internal nodes in a $3\times 3$ grid. However, each of the four edges features a different discretization: edge $AB$ contains 5 control points, edge $BC$ has 4, edge $DC$ includes 7, and edge $AD$ consists of 6 control points. For clarity, node numbering starts from the boundary and proceeds to the interior. The internal grid lines intersect edge $AD$ at points $(E,G,I)$, edge $BC$ at $(F,H,J)$, and the top edge $DC$ at $(K,11,L)$. Notably, in this example, the bottom edge $AB$ is defined such that its intermediate points coincide exactly with the orthogonal projections of the internal nodes (19 to 27) onto it, for the sake of alignment and illustrative symmetry.

For this element, the transfinite interpolation (Equation (1)) consists of the following projectors (the auxiliary points, which are necessary to define the projectors, are denoted in red ×):

$$\begin{array}{cc}\hfill P_{\xi }& =E_1\left(\xi \right)\left(N_{1,3}^{AD}\left(\eta \right)a_1+N_{2,3}^{AD}\left(\eta \right)a_{18}+N_{3,3}^{AD}\left(\eta \right)a_{17}+N_{4,3}^{AD}\left(\eta \right)a_{16}+N_{5,3}^{AD}\left(\eta \right)a_{15}+N_{6,3}^{AD}\left(\eta \right)a_{14}\right)\hfill \\ \hfill & +E_2\left(\xi \right)\left(N_{1,3}^{L\#2}\left(\eta \right)a_2+N_{2,3}^{L\#2}\left(\eta \right)a_{19}+N_{3,3}^{L\#2}\left(\eta \right)a_{22}+N_{4,3}^{L\#2}\left(\eta \right)a_{25}+N_{5,3}^{L\#2}\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_K\end{array}\right)\hfill \\ \hfill & +E_3\left(\xi \right)\left(N_{1,3}^{L\#3}\left(\eta \right)a_3+N_{2,3}^{L\#3}\left(\eta \right)a_{20}+N_{3,3}^{L\#3}\left(\eta \right)a_{23}+N_{4,3}^{L\#3}\left(\eta \right)a_{26}+N_{5,3}^{L\#3}\left(\eta \right)a_{11}\right)\hfill \\ \hfill & +E_4\left(\xi \right)\left(N_{1,3}^{L\#4}\left(\eta \right)a_4+N_{2,3}^{L\#4}\left(\eta \right)a_{21}+N_{3,3}^{L\#4}\left(\eta \right)a_{24}+N_{4,3}^{L\#4}\left(\eta \right)a_{27}+N_{5,3}^{L\#4}\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_L\end{array}\right)\hfill \\ \hfill & +E_5\left(\xi \right)\left(N_{1,3}^{DC}\left(\eta \right)a_5+N_{2,3}^{DC}\left(\eta \right)a_6+N_{3,3}^{DC}\left(\eta \right)a_7+N_{4,3}^{DC}\left(\eta \right)a_8\right.,\hfill \end{array}$$

(22)

also

$$\begin{array}{cc}\hfill P_{\eta }& =E_1\left(\eta \right)\left(N_{1,3}^{AB}\left(\xi \right)a_1+N_{2,3}^{AB}\left(\xi \right)a_2+N_{3,3}^{AB}\left(\xi \right)a_3+N_{4,3}^{AB}\left(\xi \right)a_4+N_{5,3}^{AB}\left(\xi \right)a_5\right)\hfill \\ \hfill & +E_2\left(\eta \right)\left(N_{1,3}^{L\#2}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_E\end{array}+N_{2,3}^{L\#2}\left(\xi \right)a_{19}+N_{3,3}^{L\#2}\left(\xi \right)a_{20}+N_{4,3}^{L\#2}\left(\xi \right)a_{21}+N_{5,3}^{L\#2}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_F\end{array}\right)\hfill \\ \hfill & +E_3\left(\eta \right)\left(N_{1,3}^{L\#3}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_G\end{array}+N_{2,3}^{L\#3}\left(\xi \right)a_{22}+N_{3,3}^{L\#3}\left(\xi \right)a_{23}+N_{4,3}^{L\#3}\left(\xi \right)a_{24}+N_{5,3}^{L\#3}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_H\end{array}\right)\hfill \\ \hfill & +E_4\left(\eta \right)\left(N_{1,3}^{L\#4}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_I\end{array}+N_{2,3}^{L\#4}\left(\xi \right)a_{25}+N_{3,3}^{L\#4}\left(\xi \right)a_{26}+N_{4,3}^{L\#4}\left(\xi \right)a_{27}+N_{5,3}^{L\#4}\left(\xi \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_J\end{array}\right)\hfill \\ \hfill & +E_5\left(\eta \right)\left(N_{1,3}^{DC}\left(\xi \right)a_{14}+N_{2,3}^{DC}\left(\xi \right)a_{13}+N_{3,3}^{DC}\left(\xi \right)a_{12}+N_{4,3}^{DC}\left(\xi \right)a_{11}+N_{5,3}^{DC}\left(\xi \right)a_{10}\right.\hfill \\ \hfill & \left.+N_{6,3}^{DC}\left(\xi \right)a_9+N_{7,3}^{DC}\left(\xi \right)a_8\right),\hfill \end{array}$$

(23)

and

$$\begin{array}{cc}\hfill P_{\xi \eta }& =E_1\left(\xi \right)E_1\left(\eta \right)a_1+E_2\left(\xi \right)E_1\left(\eta \right)a_2+E_3\left(\xi \right)E_1\left(\eta \right)a_3+E_4\left(\xi \right)E_1\left(\eta \right)a_4+E_5\left(\xi \right)E_1\left(\eta \right)a_4\hfill \\ \hfill & +E_1\left(\xi \right)E_2\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_E\end{array}+E_2\left(\xi \right)E_2\left(\eta \right)a_{19}+E_3\left(\xi \right)E_2\left(\eta \right)a_{20}+E_4\left(\xi \right)E_2\left(\eta \right)a_{21}+E_5\left(\xi \right)E_2\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_F\end{array}\hfill \\ \hfill & +E_1\left(\xi \right)E_3\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_G\end{array}+E_2\left(\xi \right)E_3\left(\eta \right)a_{22}+E_3\left(\xi \right)E_3\left(\eta \right)a_{23}+E_4\left(\xi \right)E_3\left(\eta \right)a_{24}+E_5\left(\xi \right)E_3\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_H\end{array}\hfill \\ \hfill & +E_1\left(\xi \right)E_4\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_I\end{array}+E_2\left(\xi \right)E_4\left(\eta \right)a_{25}+E_3\left(\xi \right)E_4\left(\eta \right)a_{26}+E_4\left(\xi \right)E_4\left(\eta \right)a_{27}+E_5\left(\xi \right)E_4\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_J\end{array}\hfill \\ \hfill & +E_1\left(\xi \right)E_5\left(\eta \right)a_{14}+E_2\left(\xi \right)E_5\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_K\end{array}+E_3\left(\xi \right)E_5\left(\eta \right)a_{11}+E_4\left(\xi \right)E_5\left(\eta \right)\begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_L\end{array}+E_5\left(\xi \right)E_5\left(\eta \right)a_8.\hfill \end{array}$$

(24)

Furthermore, if we consider that for either $s=\xi ,\eta$ we have:

$$N_{i,3}^{L\#j}\left(s\right)=E_i\left(s\right),i=1,\dots ,5\mathrm{for}\mathrm{all}\mathrm{internal}\mathrm{layers},$$

(25)

the substitution of Equation (22) to Equation (25) into the Boolean sum given by Equation (1) cancels all the auxiliary red-colored terms, and thus the bivariate function $U(\xi ,\eta )$ is approximated by:

$$\begin{array}{cc}\hfill U(\xi ,\eta )& =\left(N_{1,3}^{AD}\left(\eta \right)E_1\left(\xi \right)+N_{1,3}^{AB}\left(\xi \right)E_1\left(\eta \right)-E_1\left(\xi \right)E_1\left(\eta \right)\right)a_1\hfill \\ \hfill & +N_{2,3}^{AB}\left(\xi \right)E_1\left(\eta \right)a_2+N_{3,3}^{AB}\left(\xi \right)E_1\left(\eta \right)a_3+N_{4,3}^{AB}\left(\xi \right)E_1\left(\eta \right)a_4\hfill \\ \hfill & +\left(N_{1,3}^{BC}\left(\eta \right)E_5\left(\xi \right)+N_{5,3}^{AB}\left(\xi \right)E_1\left(\eta \right)-E_5\left(\xi \right)E_1\left(\eta \right)\right)a_5\hfill \\ \hfill & +N_{2,3}\left(\eta \right)E_5\left(\xi \right)a_6+N_{3,3}\left(\eta \right)E_5\left(\xi \right)a_7\hfill \\ \hfill & +\left(N_{4,3}^{BC}\left(\eta \right)E_5\left(\xi \right)+N_{7,3}^{DC}\left(\xi \right)E_5\left(\eta \right)-E_5\left(\xi \right)E_5\left(\eta \right)\right)a_8\hfill \\ \hfill & +N_{6,3}\left(\xi \right)E_5\left(\eta \right)a_9+N_{5,3}\left(\xi \right)E_5\left(\eta \right)a_{10}+N_{4,3}\left(\xi \right)E_5\left(\eta \right)a_{11}\hfill \\ \hfill & +N_{3,3}\left(\xi \right)E_5\left(\eta \right)a_{12}+N_{2,3}\left(\xi \right)E_5\left(\eta \right)a_{13}\hfill \\ \hfill & +\left(N_{6,3}\left(\eta \right)E_1\left(\xi \right)+N_{1,3}\left(\xi \right)E_5\left(\eta \right)-E_1\left(\xi \right)E_5\left(\eta \right)\right)a_{14}\hfill \\ \hfill & +N_{5,3}\left(\eta \right)E_1\left(\xi \right)a_{15}+N_{4,3}\left(\eta \right)E_1\left(\xi \right)a_{16}+N_{3,3}\left(\eta \right)E_1\left(\xi \right)a_{17}+N_{2,3}\left(\eta \right)E_1\left(\xi \right)a_{18}\hfill \\ \hfill & +E_2\left(\xi \right)E_2\left(\eta \right)a_{19}+E_3\left(\xi \right)E_2\left(\eta \right)a_{20}+E_4\left(\xi \right)E_2\left(\eta \right)a_{21}\hfill \\ \hfill & +E_2\left(\xi \right)E_3\left(\eta \right)a_{22}+E_3\left(\xi \right)E_3\left(\eta \right)a_{23}+E_4\left(\xi \right)E_3\left(\eta \right)a_{24}\hfill \\ \hfill & +E_2\left(\xi \right)E_4\left(\eta \right)a_{25}+E_3\left(\xi \right)E_4\left(\eta \right)a_{26}+E_4\left(\xi \right)E_4\left(\eta \right)a_{27}.\hfill \end{array}$$

(26)

The above model offers flexibility in the choice of the univariate blending and trial functions. Therefore, we can proceed as follows:

(1)

The five blending functions per direction, horizontal or vertical, can be ensured by the knot vector $\Xi =[0,0,0,0,{\displaystyle \frac{1}{2}},1,1,1,1]$ which guarantees five control points.

(2)

The trial functions along the edge $AB$ are chosen to match the blending functions in the $\xi$-direction, since all associated control points are orthogonal projections of the internal points onto this edge.

(3)

The four control points along the edge $BC$ are ensured by the knot vector $H=[0,0,0,0,1,1,1,1]$, which effectively corresponds to a set of four Bernstein polynomials of degree 3.

(4)

The seven control points along the edge $DC$ are ensured by the knot vector $\Xi =[0,0,0,0,{\displaystyle \frac{1}{4}},{\displaystyle \frac{2}{4}},{\displaystyle \frac{3}{4}},1,1,1,1]$.

(5)

The six control points along the edge $AD$ are determined by the knot vector $\Xi =[0,0,0,0,\frac{1}{3},\frac{2}{3},1,1,1,1]$.

Under these assumptions, the 27 basis functions ${\psi }_j(\xi ,\eta )$ are defined as follows:

Regardless of the element type, the formulas for the shape or basis functions remain identical. According to the Boolean sum formulation, the corner nodes $(A,B)$ on the bottom edge—discretized identically to the interior—are not influenced by the blending functions, whereas the corner nodes $(C,D)$ on the top edge are affected.

For instance, using cubic B-splines for the blending functions:

$$E_1=N_{1,3},E_2=N_{2,3},E_3=N_{3,3},E_4=N_{4,3},E_5=N_{5,3},$$

and also employing cubic B-splines for the trial functions along the edges, the resulting set of bivariate basis functions is:

$$\begin{array}{cc}\hfill & {\psi }_1(\xi ,\eta )=N_{1,3}\left(\xi \right)E_{1,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_2(\xi ,\eta )=N_{2,3}\left(\xi \right)E_{1,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_3(\xi ,\eta )=N_{3,3}\left(\xi \right)E_{1,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_4(\xi ,\eta )=N_{4,3}\left(\xi \right)E_{1,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_5(\xi ,\eta )=N_{5,3}\left(\xi \right)E_{1,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_6(\xi ,\eta )=E_{5,3}\left(\xi \right)N_{2,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_7(\xi ,\eta )=E_{5,3}\left(\xi \right)N_{3,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_8(\xi ,\eta )=B_{4,3}\left(\eta \right)E_{5,3}\left(\xi \right)-E_{5,3}\left(\xi \right)E_{5,3}\left(\eta \right)+E_{5,3}\left(\eta \right)B_{7,3}\left(\xi \right),\hfill \\ \hfill & {\psi }_9(\xi ,\eta )=N_{6,3}\left(\xi \right)E_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{10}(\xi ,\eta )=N_{5,3}\left(\xi \right)E_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{11}(\xi ,\eta )=N_{4,3}\left(\xi \right)E_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{12}(\xi ,\eta )=N_{3,3}\left(\xi \right)E_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{13}(\xi ,\eta )=N_{2,3}\left(\xi \right)E_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{14}(\xi ,\eta )=N_{1,3}\left(\xi \right)E_{5,3}\left(\eta \right)-E_{1,3}\left(\xi \right)E_{5,3}\left(\eta \right)+E_{1,3}\left(\xi \right)N_{6,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{15}(\xi ,\eta )=E_{1,3}\left(\xi \right)N_{5,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{16}(\xi ,\eta )=E_{1,3}\left(\xi \right)N_{4,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{17}(\xi ,\eta )=E_{1,3}\left(\xi \right)N_{3,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{18}(\xi ,\eta )=E_{1,3}\left(\xi \right)N_{2,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{19}(\xi ,\eta )=E_{2,3}\left(\xi \right)E_{2,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{20}(\xi ,\eta )=E_{3,3}\left(\xi \right)E_{2,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{21}(\xi ,\eta )=E_{4,3}\left(\xi \right)E_{2,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{22}(\xi ,\eta )=E_{2,3}\left(\xi \right)E_{3,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{23}(\xi ,\eta )=E_{3,3}\left(\xi \right)E_{3,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{24}(\xi ,\eta )=E_{4,3}\left(\xi \right)E_{3,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{25}(\xi ,\eta )=E_{2,3}\left(\xi \right)E_{4,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{26}(\xi ,\eta )=E_{3,3}\left(\xi \right)E_{4,3}\left(\eta \right),\hfill \\ \hfill & {\psi }_{27}(\xi ,\eta )=E_{4,3}\left(\xi \right)E_{4,3}\left(\eta \right).\hfill \end{array}$$

(27)

The graph of the 27 basis functions ${\psi }_j(\xi ,\eta )$, which are involved in Equation (26) and described by Equation (27), is shown in Figure 9. It has also been verified that the functional set $\left\{{\psi }_j(\xi ,\eta )\right\}$ satisfies the partition of unity property, i.e.,

$$\sum _{i=1}^{27}{\psi }_i(\xi ,\eta )=1.$$

(28)

It is noted that the basis functions can be categorized into three groups as follows:

(1)

Internal nodes: Defined by a simple tensor product of blending functions.

(2)

Intermediate boundary nodes: Constructed as local tensor products of trial functions along the edge and blending functions in the perpendicular direction.

(3)

Corner nodes: Formulated using a Boolean sum of three terms. Two terms correspond to projectors perpendicular to the edges connected to the corner node, while the third term is a tensor product of the associated blending functions, serving as a correction.

Although the symbol E in Equation (27) is identical to N, it was used for the sake of readability, allowing one to distinguish between blending functions (E) and trial functions (N). A graph of the set of basis functions is illustrated in Figure 9.

6. Handling T-Spline Elements

Although the previous sections referred to T-spline-like elements, this section systematically addresses the classical T-spline element. This element is defined by a T-index, which includes m knots $\{{\xi }_1,\dots ,{\xi }_m\}$ in the $\xi$-direction and n knots $\{{\eta }_1,\dots ,{\eta }_n\}$ in the $\eta$-direction of the parametric domain $ABCD$, where $0\le \xi ,\eta \le 1$.

In the common case of a cubic T-spline, the anchors of the element form a tensor-product grid of $m\times n$ potential anchor positions. However, some of these anchors may be missing, depending on the local refinement and topology of the T-mesh.

6.1. The Two Approaches

In general, two equivalent approaches can be developed, as follows [16]:

(1)

Transfinite interpolation

(2)

Background tensor-product

Similar techniques have previously been applied in conjunction with Lagrange polynomials [21] and Bernstein polynomials [20]. In those cases, it was possible to globally eliminate the degrees of freedom (DOFs) associated with the missing nodes in the $m\times n$ tensor-product grid.

Since the central approach of this paper is to first control the stations (sections or isolines) and then blend them, it is immaterial whether the focus is on a boundary edge or an internal station (inter-boundary). Naturally, when operating along a boundary edge (e.g., edge $AB$), the computations are strictly one-dimensional, as each edge is independent of the remaining data associated with the patch $ABCD$.

In contrast, knot insertion in the interior of a cubic T-spline typically affects a surrounding region spanning two steps in each direction. However, in our numerical scheme—regarded as a first approximation—this influence is intentionally neglected.

6.2. Transfinite Interpolation Versus Tensor-Product

Transfinite interpolation applies to cases where each internal node lies on a single line—either horizontal or vertical—so that the associated degree of freedom (DOF) is involved exclusively in either the $P_{\eta }$ or the $P_{\xi }$ projector, respectively. However, if the internal point lies at the intersection of a horizontal and a vertical line, an alternative approach must be adopted—for example, by taking the average of the two lines.

As an example, let us consider a T-mesh with a maximum of 9 control points per section (Figure 10). Suppose that two knots are missing on the bottom edge $AB$, and thus the function $U(\xi ,0)$ is approximated as a B-spline associated with 7 univariate functions:

$$U(\xi ,0)=\sum _{i=1}^7{N}_{i,3}\left(\xi \right)a_i,$$

as illustrated by the black-colored nodes 1 to 7 in Figure 10b (red-colored nodes 77 and 78 at $\xi =\frac{3}{6},\frac{5}{6}$ are missing). Since this set of B-splines is involved in the projector $P_{\eta }$, it follows that the associated bivariate basis functions are given by:

$${\psi }_i(\xi ,\eta )=N_{i,3}\left(\xi \right)E_1\left(\eta \right),i=1,\dots ,7.$$

The seven B-spline functions along the edge $AB$ are also plotted as solid lines in Figure 11.

The same result along the edge $AB$ can also be obtained via an alternative approach. We virtually assume that all nine control points—fully occupying the nine uniformly spaced positions—contribute to the interpolation of the primary variable U along the edge $AB$. Furthermore, we consider this complete configuration as resulting from a previously incomplete one, following knot insertion at the previously missing locations: $\xi =\frac{3}{6},\frac{5}{6}$.

The relationship between the assumed initial (incomplete) and the final complete configuration is detailed in Appendices Appendix B and Appendix C. More specifically, the complete set (${\mathbf{a}}_Q$ of size $9\times 1$) is related to the incomplete set (${\mathbf{a}}_P$ of size $7\times 1$) through the transformation matrix ${\mathbf{T}}_{\mathbf{AB}}$, as follows:

$$\begin{array}{c}\hfill \underset{\begin{array}{c}\left({\mathbf{a}}_Q\right)\\ (\mathrm{Tensor}-\mathrm{product})\end{array}}{\underbrace{\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4^{tp}\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{77}\end{array}\\ a_5^{tp}\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{78}\end{array}\\ a_6^{tp}\\ a_7\end{array}\right]}}=\underset{\left({\mathbf{T}}_{AB}^t\right)}{\underbrace{\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1/4& 3/4& 0& 0& 0\\ 0& 0& 0& 3/5& 2/5& 0& 0\\ 0& 0& 0& 3/20& 53/80& 3/16& 0\\ 0& 0& 0& 0& 1/4& 3/4& 0\\ 0& 0& 0& 0& 0& 1/2& 1/2\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]}}\underset{\begin{array}{c}\left({\mathbf{a}}_P\right)\\ \left(\mathrm{Primary}\right)\end{array}}{\underbrace{\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\\ a_7\end{array}\right]}}.\end{array}$$

(29)

It is noted that the matrix ${\mathbf{T}}_{AB}^t$ was obtained by calling the subroutine cited in Appendix C twice. The first call corresponds to the insertion of the knot at $\xi =\frac{3}{6}$, resulting in the transformation matrix ${\mathbf{T}}_1$, while the second call corresponds to the insertion of the knot at $\xi =\frac{5}{6}$, yielding the transformation matrix ${\mathbf{T}}_2$. The cumulative effect of inserting both knots leads to the overall transformation matrix ${\mathbf{T}}_{AB}={\mathbf{T}}_1{\mathbf{T}}_2$, and consequently, its transpose is given by ${\mathbf{T}}_{AB}^t={\mathbf{T}}_2^t{\mathbf{T}}_1^t$.

Substituting Equation (29) into the full B-spline sum expansion, we obtain:

$$\begin{array}{cc}\hfill U\left(\xi \right)& =\left[\begin{array}{cccc}{N}_{1,3}\left(\xi \right)& N_{2,3}\left(\xi \right)& \dots & N_{9,3}\left(\xi \right)\end{array}\right]{\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{77}\end{array}\\ a_6\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{78}\end{array}\\ a_8\\ a_9\end{array}\right]}_{\mathrm{Tensor}-\mathrm{product}}\hfill \\ \hfill & =\left[\begin{array}{cccc}{N}_{1,3}\left(\xi \right)& N_{2,3}\left(\xi \right)& \dots & N_{9,3}\left(\xi \right)\end{array}\right]\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1/4& 3/4& 0& 0& 0\\ 0& 0& 0& 3/5& 2/5& 0& 0\\ 0& 0& 0& 3/20& 53/80& 3/16& 0\\ 0& 0& 0& 0& 1/4& 3/4& 0\\ 0& 0& 0& 0& 0& 1/2& 1/2\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]{\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\\ a_7\end{array}\right]}_{\mathrm{Primary}}\hfill \end{array}$$

(30)

Obviously, Equation (30) dictates that the set of B-splines based on the seven control points (incomplete system: ${\overline{N}}_{i,3}\left(\xi \right),i=1,\dots ,7$) is a linear combination of the basis functions of the complete system ($N_{i,3}\left(\xi \right),i=1,\dots ,9$), which is defined by the product of a row vector and a transformation matrix:

$$\begin{array}{cc}\hfill {\overline{N}}_{1,3}\left(\xi \right)& =N_{1,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{2,3}\left(\xi \right)& =N_{2,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{3,3}\left(\xi \right)& =N_{3,3}\left(\xi \right)+\frac{1}{4}{N}_{4,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{4,3}\left(\xi \right)& =\frac{3}{4}{N}_{4,3}\left(\xi \right)+\frac{3}{5}{N}_{5,3}\left(\xi \right)+\frac{3}{20}{N}_{6,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{5,3}\left(\xi \right)& =\frac{2}{5}{N}_{5,3}\left(\xi \right)+\frac{53}{80}{N}_{6,3}\left(\xi \right)+\frac{1}{4}{N}_{7,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{6,3}\left(\xi \right)& =\frac{3}{16}{N}_{6,3}\left(\xi \right)+\frac{3}{4}{N}_{7,3}\left(\xi \right)+\frac{1}{2}{N}_{8,3}\left(\xi \right),\hfill \\ \hfill {\overline{N}}_{7,3}\left(\xi \right)& =\frac{1}{2}{N}_{8,3}\left(\xi \right)+N_{9,3}\left(\xi \right).\hfill \end{array}$$

(31)

Based on the full set and using Equation (31), the B-spline functions are also plotted as small circles in Figure 11. One may observe the coincidence between the B-splines directly generated by the knot vector $\Xi =[0,0,0,0,1/6,2/6,4/6,1,1,1,1]$ and those obtained via Equation (31).

6.3. Application of Knot Insertion Functions on the Entire Boundary of T-Spline

The above procedure is straightforward and can be extended from edge $AB$ to the remaining boundary edges: $BC$, $CD$, and $DA$. From a numerical standpoint, substituting the incomplete basis with the complete one does not require matrix inversion, but only a simple matrix multiplication: $\mathbf{T}{\mathbf{N}}_Q={\mathbf{N}}_P$, where ${\mathbf{N}}_P$ denotes the incomplete vector of basis functions, $\mathbf{T}$ is the transformation matrix, and ${\mathbf{N}}_Q$ is the complete vector associated with the tensor-product basis. Since the background tensor-product basis satisfies the partition of unity property, the same property holds after the substitution described by Equation (A19) in Appendix B.

Let us now complete the discussion of Figure 10b, which illustrates a T-spline patch composed of 76 control points. This configuration is derived from a full tensor-product grid of 81 points, with five control points missing from the boundary. In this extreme case, the patch interior is fully populated, while the boundary is incomplete—lacking five control points (indicated by red crosses in Figure 10b and highlighted in red below) required to complete the full tensor structure. By considering the 76 actual (or primary) control points—shown in black—and the five missing (secondary) ones, we conceptually reconstruct a complete tensor-product grid in the background. The expansion of the primary variable over this virtual tensor-product grid (in which the updated variables—after elimination using Equation (29) and subsequent steps—are denoted by the superscript `tp’) is then expressed as:

$\begin{array}{cc}\hfill U(\xi ,\eta )=& \left(E_{1x}{E}_{1y}\right)a_1+\left(E_{2x}{E}_{1y}\right)a_2+\left(E_{3x}{E}_{1y}\right)a_3+\left(E_{4x}{E}_{1y}\right)a_4^{tp}\hfill \\ \hfill +& \left(E_{5x}{E}_{1y}\right)\begin{array}{c}\hfill 1,0,0a_{77}\end{array}+\left(E_{6x}{E}_{1y}\right)a_5^{tp}+\left(E_{7x}{E}_{1y}\right)\begin{array}{c}\hfill 1,0,0a_{78}\end{array}+\left(E_{8x}{E}_{1y}\right)a_6^{tp}+\left(E_{9x}{E}_{1y}\right)a_7\hfill \\ \hfill +& \left(E_{1x}{E}_{2y}\right)a_8^{tp}+\left(E_{2x}{E}_{2y}\right)a_9+\left(E_{3x}{E}_{2y}\right)a_{10}+\left(E_{4x}{E}_{2y}\right)a_{11}\hfill \\ \hfill +& \left(E_{5x}{E}_{2y}\right)a_{12}+\left(E_{6x}{E}_{2y}\right)a_{13}+\left(E_{7x}{E}_{2y}\right)a_{14}+\left(E_{8x}{E}_{2y}\right)a_{15}+\left(E_{9x}{E}_{2y}\right)a_{16}\hfill \\ \hfill +& \left(E_{1x}{E}_{3y}\right)\begin{array}{c}\hfill 1,0,0a_{80}^{tp}\end{array}+\left(E_{2x}{E}_{3y}\right)a_{17}+\left(E_{3x}{E}_{3y}\right)a_{18}+\left(E_{4x}{E}_{3y}\right)a_{19}\hfill \\ \hfill +& \left(E_{5x}{E}_{3y}\right)a_{20}+\left(E_{6x}{E}_{3y}\right)a_{21}+\left(E_{7x}{E}_{3y}\right)a_{22}+\left(E_{8x}{E}_{3y}\right)a_{23}+\left(E_{9x}{E}_{3y}\right)a_{24}\hfill \\ \hfill +& \left(E_{1x}{E}_{4y}\right)a_{25}^{tp}+\left(E_{2x}{E}_{4y}\right)a_{26}+\left(E_{3x}{E}_{4y}\right)a_{27}+\left(E_{4x}{E}_{4y}\right)a_{28}\hfill \\ \hfill +& \left(E_{5x}{E}_{4y}\right)a_{29}+\left(E_{6x}{E}_{4y}\right)a_{30}+\left(E_{7x}{E}_{4y}\right)a_{31}+\left(E_{8x}{E}_{4y}\right)a_{32}+\left(E_{9x}{E}_{4y}\right)a_{33}\hfill \\ \hfill +& \left(E_{1x}{E}_{5y}\right)a_{34}+\left(E_{2x}{E}_{5y}\right)a_{35}+\left(E_{3x}{E}_{5y}\right)a_{36}+\left(E_{4x}{E}_{5y}\right)a_{37}\hfill \\ \hfill +& \left(E_{5x}{E}_{5y}\right)a_{38}+\left(E_{6x}{E}_{5y}\right)a_{39}+\left(E_{7x}{E}_{5y}\right)a_{40}+\left(E_{8x}{E}_{5y}\right)a_{41}+\left(E_{9x}{E}_{5y}\right)a_{42}^{tp}\hfill \\ \hfill +& \left(E_{1x}{E}_{6y}\right)a_{43}^{tp}+\left(E_{2x}{E}_{6y}\right)a_{44}+\left(E_{3x}{E}_{6y}\right)a_{45}+\left(E_{4x}{E}_{6y}\right)a_{46}\hfill \\ \hfill +& \left(E_{5x}{E}_{6y}\right)a_{47}+\left(E_{6x}{E}_{6y}\right)a_{48}+\left(E_{7x}{E}_{6y}\right)a_{49}+\left(E_{8x}{E}_{6y}\right)a_{50}+\left(E_{9x}{E}_{6y}\right)\begin{array}{c}\hfill 1,0,0a_{79}\end{array}\hfill \\ \hfill +& \left(E_{1x}{E}_{7y}\right)\begin{array}{c}\hfill 1,0,0a_{81}^{tp}\end{array}+\left(E_{2x}{E}_{7y}\right)a_{51}+\left(E_{3x}{E}_{7y}\right)a_{52}+\left(E_{4x}{E}_{7y}\right)a_{53}\hfill \\ \hfill +& \left(E_{5x}{E}_{7y}\right)a_{54}+\left(E_{6x}{E}_{7y}\right)a_{55}+\left(E_{7x}{E}_{7y}\right)a_{56}+\left(E_{8x}{E}_{7y}\right)a_{57}+\left(E_{9x}{E}_{7y}\right)a_{58}^{tp}\hfill \\ \hfill +& \left(E_{1x}{E}_{8y}\right)a_{59}^{tp}+\left(E_{2x}{E}_{8y}\right)a_{60}+\left(E_{3x}{E}_{8y}\right)a_{61}+\left(E_{4x}{E}_{8y}\right)a_{62}\hfill \\ \hfill +& \left(E_{5x}{E}_{8y}\right)a_{63}+\left(E_{6x}{E}_{8y}\right)a_{64}+\left(E_{7x}{E}_{8y}\right)a_{65}+\left(E_{8x}{E}_{8y}\right)a_{66}+\left(E_{9x}{E}_{8y}\right)a_{67}\hfill \\ \hfill +& \left(E_{1x}{E}_{9y}\right)a_{68}+\left(E_{2x}{E}_{9y}\right)a_{69}+\left(E_{3x}{E}_{9y}\right)a_{70}+\left(E_{4x}{E}_{9y}\right)a_{71}\hfill \\ \hfill +& \left(E_{5x}{E}_{9y}\right)a_{72}+\left(E_{6x}{E}_{9y}\right)a_{73}+\left(E_{7x}{E}_{9y}\right)a_{74}+\left(E_{8x}{E}_{9y}\right)a_{75}+\left(E_{9x}{E}_{9y}\right)a_{76}\hfill \end{array}$

(29)

The main concept is that the basis functions associated with the 76 primary DOFs of this T-spline element will be produced by eliminating the red-colored secondary DOFs. This idea has been successfully applied in conjunction with Lagrange polynomials, where the values of the incomplete polynomials were calculated at the missing points and then substituted into virtual polynomials of higher degree in the tensor product [21], in a global manner for the entire station under consideration.

However, in the context of B-splines, if we wish to preserve the same polynomial degree, a different technique is required to substitute the red-colored variables.

Since the tensor-product configuration is fully controlled and the associated basis functions are standard B-splines, it can be considered as a state produced by virtual knot insertion into a coarser mesh—namely, the given T-mesh (defined by its index space). In this example, we must eliminate two DOFs along edge $AB$, one DOF along edge $BC$, and two DOFs along edge $AD$.

In total, the model includes 76 primary nodes and 5 auxiliary (secondary) nodes—located on the boundary—which merely complete (or close) the tensor-product structure in the background.

More specifically, the treatment of edge $AB$ has already been discussed in Equation (29).

Similarly, for the edge $BC$ we have:

$$\begin{array}{c}\hfill {\left[\begin{array}{c}{a}_7\\ a_{16}\\ a_{24}\\ a_{33}\\ a_{42}\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{79}\end{array}\\ a_{58}\\ a_{67}\\ a_{76}\end{array}\right]}_{\mathrm{Tensor}}-\mathrm{product}=\underset{\left({\mathbf{T}}_{BC}^T\right)}{\underbrace{\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 1/4& 3/4& 0& 0& 0\\ 0& 0& 0& 0& 1/2& 1/2& 0& 0\\ 0& 0& 0& 0& 0& 2/3& 1/3& 0\\ 0& 0& 0& 0& 0& 2/3& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]}}{\left[\begin{array}{c}{a}_7\\ a_{16}\\ a_{24}\\ a_{33}\\ a_{42}\\ a_{58}\\ a_{67}\\ a_{76}\end{array}\right]}_{\mathrm{Primary}}.\end{array}$$

(32)

Since the top edge $DC$ is complete, we keep it intact and proceed to the final edge $AD$, for which we have:

$$\begin{array}{c}\hfill {\left[\begin{array}{c}{a}_1\\ a_8\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{80}\end{array}\\ a_{25}\\ a_{34}\\ a_{43}\\ \begin{array}{c}\hfill named-contentcolortypergbcontent-type1,0,0a_{81}\end{array}\\ a_{59}\\ a_{68}\end{array}\right]}_{\mathrm{Tensor}}-\mathrm{product}=\underset{\left({\mathbf{T}}_{AD}^T\right)}{\underbrace{\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 1/2& 1/2& 0& 0& 0& 0& 0\\ 0& 2/3& 1/3& 0& 0& 0& 0\\ 0& 0& 3/4& 1/4& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1/4& 3/4& 0& 0\\ 0& 0& 0& 0& 1/3& 2/3& 0\\ 0& 0& 0& 0& 0& 1/2& 1/2\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]}}{\left[\begin{array}{c}{a}_1\\ a_8\\ a_{25}\\ a_{34}\\ a_{43}\\ a_{59}\\ a_{68}\end{array}\right]}_{\mathrm{Primary}}.\end{array}$$

(33)

Substituting the three sets of equations—Equation (29), Equation (32), and Equation (33)—each corresponding to an active edge, into the tensor product defined by Equation (Section 6.3), and performing factorization, we obtain:

$$U(\xi ,\eta )=\sum _{i=1}^{76}{\psi }_i(\xi ,\eta )a_i,$$

(34)

where the basis functions are given by:

$$resizeboxwidth298.8987pt\begin{array}{cccc}\hfill {\psi }_1& =E_{1x}{E}_{1y}+\frac{1}{2}{E}_{1x}{E}_{2y}\hfill & \hfill {\psi }_2& =E_{2x}{E}_{1y}\hfill \\ \hfill {\psi }_3& =E_{3x}{E}_{1y}+\frac{1}{4}{E}_{4x}{E}_{1y}\hfill & \hfill {\psi }_4& =\frac{3}{4}{E}_{4x}{E}_{1y}+\frac{3}{5}{E}_{5x}{E}_{1y}+\frac{3}{20}{E}_{6x}{E}_{1y}\hfill \\ \hfill {\psi }_5& =\frac{2}{5}{E}_{5x}{E}_{1y}+\frac{53}{80}{E}_{6x}{E}_{1y}+\frac{1}{4}{E}_{7x}{E}_{1y}\hfill & \hfill {\psi }_6& =\frac{3}{16}{E}_{6x}{E}_{1y}+\frac{3}{4}{E}_{7x}{E}_{1y}+\frac{1}{2}{E}_{8x}{E}_{1y}\hfill \\ \hfill {\psi }_7& =\frac{1}{2}{E}_{8x}{E}_{1y}+E_{9x}{E}_{1y}\hfill & \hfill {\psi }_8& =\frac{1}{2}{E}_{1x}{E}_{2y}+\frac{2}{3}{E}_{1x}{E}_{3y}\hfill \\ \hfill {\psi }_9& =E_{2x}{E}_{2y}\hfill & \hfill {\psi }_{10}& =E_{3x}{E}_{2y}\hfill \\ \hfill {\psi }_{11}& =E_{4x}{E}_{2y}\hfill & \hfill {\psi }_{12}& =E_{5x}{E}_{2y}\hfill \\ \hfill {\psi }_{13}& =E_{6x}{E}_{2y}\hfill & \hfill {\psi }_{14}& =E_{7x}{E}_{2y}\hfill \\ \hfill {\psi }_{15}& =E_{8x}{E}_{2y}\hfill & \hfill {\psi }_{16}& =E_{9x}{E}_{2y}\hfill \\ \hfill {\psi }_{17}& =E_{2x}{E}_{3y}\hfill & \hfill {\psi }_{18}& =E_{3x}{E}_{3y}\hfill \\ \hfill {\psi }_{19}& =E_{4x}{E}_{3y}\hfill & \hfill {\psi }_{20}& =E_{5x}{E}_{3y}\hfill \\ \hfill {\psi }_{21}& =E_{6x}{E}_{3y}\hfill & \hfill {\psi }_{22}& =E_{7x}{E}_{3y}\hfill \\ \hfill {\psi }_{23}& =E_{8x}{E}_{3y}\hfill & \hfill {\psi }_{24}& =E_{9x}{E}_{3y}\hfill \\ \hfill {\psi }_{25}& =\frac{1}{3}{E}_{1x}{E}_{3y}+\frac{3}{4}{E}_{1x}{E}_{4y}\hfill & \hfill {\psi }_{26}& =E_{2x}{E}_{4y}\hfill \\ \hfill {\psi }_{27}& =E_{3x}{E}_{4y}\hfill & \hfill {\psi }_{28}& =E_{4x}{E}_{4y}\hfill \\ \hfill {\psi }_{29}& =E_{5x}{E}_{4y}\hfill & \hfill {\psi }_{30}& =E_{6x}{E}_{4y}\hfill \\ \hfill {\psi }_{31}& =E_{7x}{E}_{4y}\hfill & \hfill {\psi }_{32}& =E_{8x}{E}_{4y}\hfill \\ \hfill {\psi }_{33}& =E_{9x}{E}_{4y}+\frac{1}{4}{E}_{9x}{E}_{5y}\hfill & \hfill {\psi }_{34}& =\frac{1}{4}{E}_{1x}{E}_{4y}+E_{1x}{E}_{5y}+\frac{1}{4}{E}_{1x}{E}_{6y}\hfill \\ \hfill {\psi }_{35}& =E_{2x}{E}_{5y}\hfill & \hfill {\psi }_{36}& =E_{3x}{E}_{5y}\hfill \\ \hfill {\psi }_{37}& =E_{4x}{E}_{5y}\hfill & \hfill {\psi }_{38}& =E_{5x}{E}_{5y}\hfill \\ \hfill {\psi }_{39}& =E_{6x}{E}_{5y}\hfill & \hfill {\psi }_{40}& =E_{7x}{E}_{5y}\hfill \\ \hfill {\psi }_{41}& =E_{8x}{E}_{5y}\hfill & \hfill {\psi }_{42}& =\frac{3}{4}{E}_{9x}{E}_{5y}+\frac{1}{2}{E}_{9x}{E}_{6y}\hfill \\ \hfill {\psi }_{43}& =\frac{3}{4}{E}_{1x}{E}_{6y}+\frac{1}{3}{E}_{1x}{E}_{7y}\hfill & \hfill {\psi }_{44}& =E_{2x}{E}_{6y}\hfill \\ \hfill {\psi }_{45}& =E_{3x}{E}_{6y}\hfill & \hfill {\psi }_{46}& =E_{4x}{E}_{6y}\hfill \\ \hfill {\psi }_{47}& =E_{5x}{E}_{6y}\hfill & \hfill {\psi }_{48}& =E_{6x}{E}_{6y}\hfill \\ \hfill {\psi }_{49}& =E_{7x}{E}_{6y}\hfill & \hfill {\psi }_{50}& =E_{8x}{E}_{6y}\hfill \\ \hfill {\psi }_{51}& =E_{2x}{E}_{7y}\hfill & \hfill {\psi }_{52}& =E_{3x}{E}_{7y}\hfill \\ \hfill {\psi }_{53}& =E_{4x}{E}_{7y}\hfill & \hfill {\psi }_{54}& =E_{5x}{E}_{7y}\hfill \\ \hfill {\psi }_{55}& =E_{6x}{E}_{7y}\hfill & \hfill {\psi }_{56}& =E_{7x}{E}_{7y}\hfill \\ \hfill {\psi }_{57}& =E_{8x}{E}_{7y}\hfill & \hfill {\psi }_{58}& =\frac{1}{2}{E}_{9x}{E}_{6y}+\frac{2}{3}{E}_{9x}{E}_{7y}\hfill \\ \hfill {\psi }_{59}& =\frac{2}{3}{E}_{1x}{E}_{7y}+\frac{1}{2}{E}_{1x}{E}_{8y}\hfill & \hfill {\psi }_{60}& =E_{2x}{E}_{8y}\hfill \\ \hfill {\psi }_{61}& =E_{3x}{E}_{8y}\hfill & \hfill {\psi }_{62}& =E_{4x}{E}_{8y}\hfill \\ \hfill {\psi }_{63}& =E_{5x}{E}_{8y}\hfill & \hfill {\psi }_{64}& =E_{6x}{E}_{8y}\hfill \\ \hfill {\psi }_{65}& =E_{7x}{E}_{8y}\hfill & \hfill {\psi }_{66}& =E_{8x}{E}_{8y}\hfill \\ \hfill {\psi }_{67}& =\frac{1}{3}{E}_{9x}{E}_{7y}+E_{9x}{E}_{8y}\hfill & \hfill {\psi }_{68}& =\frac{1}{2}{E}_{1x}{E}_{8y}+E_{1x}{E}_{9y}\hfill \\ \hfill {\psi }_{69}& =E_{2x}{E}_{9y}\hfill & \hfill {\psi }_{70}& =E_{3x}{E}_{9y}\hfill \\ \hfill {\psi }_{71}& =E_{4x}{E}_{9y}\hfill & \hfill {\psi }_{72}& =E_{5x}{E}_{9y}\hfill \\ \hfill {\psi }_{73}& =E_{6x}{E}_{9y}\hfill & \hfill {\psi }_{74}& =E_{7x}{E}_{9y}\hfill \\ \hfill {\psi }_{75}& =E_{8x}{E}_{9y}\hfill & \hfill {\psi }_{76}& =E_{9x}{E}_{9y}\hfill \end{array}$$

(35)

In Equation (35), we have assumed that $E_{1x}$ through $E_{9x}$ represent the B-splines $N_{i,3}\left(\xi \right)$ for $i=1,\dots ,9$, while $E_{1y}$ through $E_{9y}$ correspond to the B-splines $N_{i,3}\left(\eta \right)$ for $i=1,\dots ,9$. The symbol E is also consistent with the interpretation of blending functions. However, the tensor formulation—whether based on Lagrange, Bernstein, or B-spline functions—is fundamentally governed by one of the three projectors: $P_{\xi }$, $P_{\eta }$, or $P_{\xi \eta }$, as thoroughly discussed in Ref. [9].

The graph of these basis functions is shown in Figure 12. One may observe that:

(1)

None of the degrees of freedom (DOFs) associated with the top edge are affected, as this edge is fully occupied.

(2)

Among the 76 DOFs, the only affected ones are those associated with nodes $(1,3,4,5,6,7,8,25,33,34,42,43,58,59,67,68)$. All of these nodes lie on the affected boundary, comprising the three edges: $AB$, $BC$, and $AD$.

(3)

Although node 2 lies on edge $AB$, it is not affected by the nearest inserted knot at point `77’, since it is located three units away.

(4)

None of the interior tensor-product DOFs are influenced by the coarse mesh along the boundary.

(5)

The shape and its parametric representation remain unchanged.

(6)

The partition of unity is satisfied a priori for all $(\xi ,\eta )\in \Omega$:

$$\sum _{i=1}^{76}{\varphi }_i(\xi ,\eta )=1,$$

(36)

and therefore, no normalization is required.

(7)

The $L_2$ error for the tensor-product element (81-node) is $5.6581\times {10}^{-4}\%$.

(8)

The $L_2$ error for the 76-node element was found to be $0.0011\%$, which is approximately double that of the tensor-product element.

6.4. Treatment of Internal Nodes

Regarding the internal nodes, any missing knot can be treated in an analogous manner to the boundary case. The key difference is that when a gap is filled by inserting a knot, the influence may extend in two parametric directions. The simplest scenario occurs when the missing knot affects only a single station—e.g., a horizontal line. This concept is directly applicable when using Lagrange or Bernstein polynomials.

Let us consider a simple case inspired by Ref. [25]. Specifically, three internal vertices are missing along three horizontal lines at $\xi =\frac{1}{2}$, corresponding to the parametric positions $\eta =\frac{1}{6},\frac{2}{6},\frac{3}{6}$, as illustrated in Figure 10c. For each of these missing vertices, we define the following transformation matrix:

$$\mathbf{T}=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 1/4& 0& 0& 0& 0& 0\\ 0& 0& 0& 3/4& 1/2& 0& 0& 0& 0\\ 0& 0& 0& 0& 1/2& 3/4& 0& 0& 0\\ 0& 0& 0& 0& 0& 1/4& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right],$$

(37)

so that the tensor-product degrees of freedom (DOFs), denoted by ${\mathbf{a}}_Q$, are related to the actual ones, ${\mathbf{a}}_P$, through the transformation

$${\mathbf{a}}_Q={\mathbf{T}}^t{\mathbf{a}}_P.$$

In more detail, we impose the following sets of constraints for each missing internal vertex:

$$\mathrm{Isoline}\eta =\frac{1}{6}:a_{19}^{\mathrm{tp}}=\frac{1}{4}{a}_{18}+\frac{3}{4}{a}_{19},a_{77}^{\mathrm{tp}}=\frac{1}{2}{a}_{19}+\frac{1}{2}{a}_{20},a_{20}^{\mathrm{tp}}=\frac{3}{4}{a}_{20}+\frac{1}{4}{a}_{21}.$$

(38)

$$\mathrm{Isoline}\eta =\frac{2}{6}:a_{27}^{\mathrm{tp}}=\frac{1}{4}{a}_{26}+\frac{3}{4}{a}_{27},a_{78}^{\mathrm{tp}}=\frac{1}{2}{a}_{27}+\frac{1}{2}{a}_{28},a_{28}^{\mathrm{tp}}=\frac{3}{4}{a}_{28}+\frac{1}{4}{a}_{29}.$$

(39)

$$\mathrm{Isoline}\eta =\frac{3}{6}:a_{35}^{\mathrm{tp}}=\frac{1}{4}{a}_{34}+\frac{3}{4}{a}_{35},a_{79}^{\mathrm{tp}}=\frac{1}{2}{a}_{35}+\frac{1}{2}{a}_{36},a_{36}^{\mathrm{tp}}=\frac{3}{4}{a}_{36}+\frac{1}{4}{a}_{37}.$$

(40)

In Eqs. (38)–(40), the superscript “tp” (standing for tensor product) in $a_k^{\mathrm{tp}}$ indicates that the k-th degree of freedom (DOF), shown in Figure 10c and subject to the constraint, represents the final value after the elimination process. This final value is subsequently employed in the expression of the initial tensor-product form, whereas the variables on the right-hand side correspond to the values prior to elimination.

After accounting for the five missing boundary nodes (Eqs. (29), (32), and (33)), along with the three internal nodes (Eqs. (38) to (40)), the resulting graph of the basis function set is illustrated in Figure 13.

6.5. Intersected Isolines

While handling missing vertices along the isolines of a T-mesh is relatively straightforward—as demonstrated in Section 6.4—complications arise when a missing vertex is located at the intersection of two isolines. In such cases, the numerical results may vary depending on the choice of isoline along which the elimination is performed.

Extensive numerical experimentation suggests a more robust strategy: perform the elimination independently along both isolines and then take the arithmetic mean of the two results as the final expression for elimination.

A simple test case is illustrated in Figure 14, where the central node (located at $\xi =\eta =\frac{1}{2}$) has been omitted from the tensor-product grid of 81 nodes.

The imposed constraints are:

Regarding the horizontal $\xi$-direction:

$$\mathrm{Isoline}\eta =\frac{1}{2}:a_{40}^{\mathrm{tp}}=\frac{1}{4}{a}_{39}+\frac{3}{4}{a}_{40},a_{81,x}^{\mathrm{tp}}=\frac{1}{2}{a}_{40}+\frac{1}{2}{a}_{41},a_{41}^{\mathrm{tp}}=\frac{3}{4}{a}_{41}+\frac{1}{4}{a}_{42}.$$

(41)

Regarding the vertical $\eta$-direction:

$$\mathrm{Isoline}\xi =\frac{1}{2}:a_{32}^{\mathrm{tp}}=\frac{1}{4}{a}_{23}+\frac{3}{4}{a}_{32},a_{81,y}^{\mathrm{tp}}=\frac{1}{2}{a}_{32}+\frac{1}{2}{a}_{49},a_{49}^{\mathrm{tp}}=\frac{3}{4}{a}_{49}+\frac{1}{4}{a}_{58}.$$

(42)

After substituting Equation (41) and Equation (42) into the tensor-product expression, we find that only the eight blue-colored vertices in Figure 14 are affected. Using the MATLAB code provided in Appendix D, the reader may verify that the corresponding (modified) basis functions are given by:

$$\begin{array}{cccccc}\hfill {\psi }_1& =E_{1x}{E}_{1y}\hfill & \hfill {\psi }_2& =E_{2x}{E}_{1y}\hfill & \hfill {\psi }_3& =E_{3x}{E}_{1y}\hfill \\ \hfill {\psi }_4& =E_{4x}{E}_{1y}\hfill & \hfill {\psi }_5& =E_{5x}{E}_{1y}\hfill & \hfill {\psi }_6& =E_{6x}{E}_{1y}\hfill \\ \hfill {\psi }_7& =E_{7x}{E}_{1y}\hfill & \hfill {\psi }_8& =E_{8x}{E}_{1y}\hfill & \hfill {\psi }_9& =E_{9x}{E}_{1y}\hfill \\ \hfill {\psi }_{10}& =E_{1x}{E}_{2y}\hfill & \hfill {\psi }_{11}& =E_{2x}{E}_{2y}\hfill & \hfill {\psi }_{12}& =E_{3x}{E}_{2y}\hfill \\ \hfill {\psi }_{13}& =E_{4x}{E}_{2y}\hfill & \hfill {\psi }_{14}& =E_{5x}{E}_{2y}\hfill & \hfill {\psi }_{15}& =E_{6x}{E}_{2y}\hfill \\ \hfill {\psi }_{16}& =E_{7x}{E}_{2y}\hfill & \hfill {\psi }_{17}& =E_{8x}{E}_{2y}\hfill & \hfill {\psi }_{18}& =E_{9x}{E}_{2y}\hfill \\ \hfill {\psi }_{19}& =E_{1x}{E}_{3y}\hfill & \hfill {\psi }_{20}& =E_{2x}{E}_{3y}\hfill & \hfill {\psi }_{21}& =E_{3x}{E}_{3y}\hfill \\ \hfill {\psi }_{22}& =E_{4x}{E}_{3y}\hfill & \hfill {\psi }_{23}& =E_{5x}{E}_{3y}+\frac{1}{4}{E}_{5x}{E}_{4y}\hfill & \hfill {\psi }_{24}& =E_{6x}{E}_{3y}\hfill \\ \hfill {\psi }_{25}& =E_{7x}{E}_{3y}\hfill & \hfill {\psi }_{26}& =E_{8x}{E}_{3y}\hfill & \hfill {\psi }_{27}& =E_{9x}{E}_{3y}\hfill \\ \hfill {\psi }_{28}& =E_{1x}{E}_{4y}\hfill & \hfill {\psi }_{29}& =E_{2x}{E}_{4y}\hfill & \hfill {\psi }_{30}& =E_{3x}{E}_{4y}\hfill \\ \hfill {\psi }_{31}& =E_{4x}{E}_{4y}\hfill & \hfill {\psi }_{32}& =\frac{3}{4}{E}_{5x}{E}_{4y}+\frac{1}{4}{E}_{5x}{E}_{5y}\hfill & \hfill {\psi }_{33}& =E_{6x}{E}_{4y}\hfill \\ \hfill {\psi }_{34}& =E_{7x}{E}_{4y}\hfill & \hfill {\psi }_{35}& =E_{8x}{E}_{4y}\hfill & \hfill {\psi }_{36}& =E_{9x}{E}_{4y}\hfill \\ \hfill {\psi }_{37}& =E_{1x}{E}_{5y}\hfill & \hfill {\psi }_{38}& =E_{2x}{E}_{5y}\hfill & \hfill {\psi }_{39}& =E_{3x}{E}_{5y}+\frac{1}{4}{E}_{4x}{E}_{5y}\hfill \\ \hfill {\psi }_{40}& =\frac{3}{4}{E}_{4x}{E}_{5y}+\frac{1}{4}{E}_{5x}{E}_{5y}\hfill & \hfill {\psi }_{41}& =\frac{1}{4}{E}_{5x}{E}_{5y}+\frac{3}{4}{E}_{6x}{E}_{5y}\hfill & \hfill {\psi }_{42}& =\frac{1}{4}{E}_{6x}{E}_{5y}+E_{7x}{E}_{5y}\hfill \\ \hfill {\psi }_{43}& =E_{8x}{E}_{5y}\hfill & \hfill {\psi }_{44}& =E_{9x}{E}_{5y}\hfill & \hfill {\psi }_{45}& =E_{1x}{E}_{6y}\hfill \\ \hfill {\psi }_{46}& =E_{2x}{E}_{6y}\hfill & \hfill {\psi }_{47}& =E_{3x}{E}_{6y}\hfill & \hfill {\psi }_{48}& =E_{4x}{E}_{6y}\hfill \\ \hfill {\psi }_{49}& =\frac{1}{4}{E}_{5x}{E}_{5y}+\frac{3}{4}{E}_{5x}{E}_{6y}\hfill & \hfill {\psi }_{50}& =E_{6x}{E}_{6y}\hfill & \hfill {\psi }_{51}& =E_{7x}{E}_{6y}\hfill \\ \hfill {\psi }_{52}& =E_{8x}{E}_{6y}\hfill & \hfill {\psi }_{53}& =E_{9x}{E}_{6y}\hfill & \hfill {\psi }_{54}& =E_{1x}{E}_{7y}\hfill \\ \hfill {\psi }_{55}& =E_{2x}{E}_{7y}\hfill & \hfill {\psi }_{56}& =E_{3x}{E}_{7y}\hfill & \hfill {\psi }_{57}& =E_{4x}{E}_{7y}\hfill \\ \hfill {\psi }_{58}& =\frac{1}{4}{E}_{5x}{E}_{6y}+E_{5x}{E}_{7y}\hfill & \hfill {\psi }_{59}& =E_{6x}{E}_{7y}\hfill & \hfill {\psi }_{60}& =E_{7x}{E}_{7y}\hfill \\ \hfill {\psi }_{61}& =E_{8x}{E}_{7y}\hfill & \hfill {\psi }_{62}& =E_{9x}{E}_{7y}\hfill & \hfill {\psi }_{63}& =E_{1x}{E}_{8y}\hfill \\ \hfill {\psi }_{64}& =E_{2x}{E}_{8y}\hfill & \hfill {\psi }_{65}& =E_{3x}{E}_{8y}\hfill & \hfill {\psi }_{66}& =E_{4x}{E}_{8y}\hfill \\ \hfill {\psi }_{67}& =E_{5x}{E}_{8y}\hfill & \hfill {\psi }_{68}& =E_{6x}{E}_{8y}\hfill & \hfill {\psi }_{69}& =E_{7x}{E}_{8y}\hfill \\ \hfill {\psi }_{70}& =E_{8x}{E}_{8y}\hfill & \hfill {\psi }_{71}& =E_{9x}{E}_{8y}\hfill & \hfill {\psi }_{72}& =E_{1x}{E}_{9y}\hfill \\ \hfill {\psi }_{73}& =E_{2x}{E}_{9y}\hfill & \hfill {\psi }_{74}& =E_{3x}{E}_{9y}\hfill & \hfill {\psi }_{75}& =E_{4x}{E}_{9y}\hfill \\ \hfill {\psi }_{76}& =E_{5x}{E}_{9y}\hfill & \hfill {\psi }_{77}& =E_{6x}{E}_{9y}\hfill & \hfill {\psi }_{78}& =E_{7x}{E}_{9y}\hfill \\ \hfill {\psi }_{79}& =E_{8x}{E}_{9y}\hfill & \hfill {\psi }_{80}& =E_{9x}{E}_{9y}\hfill & \hfill .\end{array}$$

(43)

One may observe that only the blue-colored vertices are influenced—a well-known characteristic in standard T-spline practice. The complete set of basis functions is illustrated in Figure 15.

6.6. Continuity Issues

All bivariate basis functions are ultimately constructed as tensor products of blending functions and trial functions—both of which are a priori continuous. Consequently, the numerical solution, being an algebraic sum of such tensor products, inherits this continuity.

For B-splines, where the blending functions are piecewise polynomials of degree $p_{\mathrm{blend}}$ in the $\xi$-direction and the trial functions are of degree $p_{\mathrm{trial}}$ in the $\eta$-direction, the continuity at internal breakpoints is $C^{p_{\mathrm{blend}}-1}$ and $C^{p_{\mathrm{trial}}-1}$, respectively.

In contrast, tensor-product Lagrange (or Bernstein) polynomials of degrees $(p,q)$ yield bivariate shape (or basis) functions that are infinitely differentiable within the open patch—effectively $C^p\times C^q$ continuous throughout the interior, which is treated as a single macroelement. However, globally, the assembled surface exhibits only $C^0\times C^0$ continuity across element (patch) interfaces.

7. Software Issues

All computations were performed within the standard MATLAB® environment (version R2024b), which includes the Spline Toolbox. When applicable, the evaluation of B-splines, $N_{i,3}(\xi ,\eta )$, and their derivatives was carried out using the built-in function spcol, originally developed in Fortran77 by Carl de Boor and later ported to MATLAB [26].

In addition to built-in tools, a suite of self-contained in-house codes was developed to support the workflow, comprising: (i) pre-processing modules for mesh generation and data preparation, (ii) analysis routines for solving the governing equations, and (iii) post-processing utilities for visualization and error assessment.

Given that the primary objective was proof-of-concept—specifically, to evaluate the accuracy of transfinite elements under various unstructured configurations—the computational strategy was tailored accordingly.

7.1. Lagrange and Bernstein Polynomials

A MATLAB® function named lagrange, which computes univariate Lagrange polynomials (both uniform and non-uniform) along with their first derivatives, is available in [23]. Additionally, a tensor-product implementation of uniform Lagrange polynomials, named shape_full_lagrange, is provided in [23]. The output of this function includes the bivariate shape functions and their partial derivatives (stored in the array shp), as well as the determinant of the Jacobian matrix (variable xsj). With minor modifications, this function can be extended to handle non-uniform polynomials.

An auxiliary in-house function, Bernstein(tau, p), for computing Bernstein–Bézier polynomials is listed in Appendix E. Gaussian integration was performed using the subroutine lgwt (see, [27]).

For each pair of parameters $(\xi ,\eta )$, a subroutine shapeX was invoked. This subroutine utilizes either closed-form expressions from the main text or formulations from Appendix D, depending on context. Its general form is:

[shp, xsj] = shapeX(xi, eta, XL, NEL, ...)

where:

- shp contains the basis functions and their partial derivatives:

shp(1,i) = $\partial N/\partial x$,

shp(2,i) = $\partial N/\partial y$,

shp(3,i) = N, for $i=1,\dots ,\mathtt{NEL}$.

- xsj is the determinant of the Jacobian matrix.

- xi, eta are the parametric coordinates $(\xi ,\eta )$.

- XL contains the Cartesian coordinates of the element nodes.

- NEL is the number of nodes in the element.

For the simplest case—a 12-node transfinite element described by Equation (19)—the complete subroutine is provided in Appendix F.

This subroutine structure aligns with standard finite element implementation practices, as discussed in Ref. [2]. The full computer code is included in Appendix G.

7.2. B-Spline Interpolation

The MATLAB® function spcol, part of the Spline Toolbox, generates a collocation matrix for B-splines by evaluating spline basis functions and their derivatives at specified sites.

In the context of B-spline interpolation, the function shapeX is enhanced to accept additional input arguments: the polynomial degrees px, py, and the knot vectors knotx, knoty. This extended formulation enables direct comparison of different blending and trial functions within a unified in-house finite element framework, without requiring special treatment for their local support characteristics.

Numerical integration was performed over cells defined by adjacent breakpoints—i.e., cubic B-spline elements in the sense of isogeometric analysis (IGA)—using a $4\times 4$ or $6\times 6$ Gaussian quadrature scheme for rectangular or curvilinear elements, respectively. Notably, both Lagrange and Bézier elements do not require domain subdivision and can be integrated using a global quadrature over the entire patch.

The determination of integration cells (elements) can be efficiently performed by applying the MATLAB function unique to the knot vectors. Although optimizing the code for computational efficiency using a connectivity vector IEN is straightforward, this enhancement has been deferred to future work.

7.3. T-Spline

In this context, a dedicated routine was developed to compute the transformation matrix $\mathbf{T}$ associated with knot insertion (see, for example, Equation (37), and for further details, refer to Appendix B and Appendix C). These formulas are ultimately incorporated manually into the implementation provided in Appendix D.

8. Numerical Solution of Boundary-Value Problems

To illustrate the proposed method, four representative examples are presented below. The first involves a rectangular domain, while the second considers a semicircular ring (i.e., a half-annulus), both addressing the solution of the Laplace equation.

For these two cases, the $L_2$ error norm (expressed as a percentage) is computed over the entire domain $\Omega$ using the following formula:

$$L_2={\displaystyle \frac{{\int }_{\Omega }{\left(U_{\mathrm{calculated}}-U_{\mathrm{exact}}\right)}^{2}d\Omega }{{\int }_{\Omega }{\left(U_{\mathrm{exact}}\right)}^{2}d\Omega }}\times 100(\%).$$

(44)

The third example concerns an eigenvalue problem, for which details on the error are provided in the dedicated subSection 8.3. Finally, the fourth example is a patch test in plane elasticity (Section 8.4).

8.1. Example 1

A square plate $ABCD$ is subjected to steady-state heat conduction, with boundary conditions illustrated inFigure 16. The exact solution is given by:

$$U(x,y)=U_m{\displaystyle \frac{sinh\left(\frac{\pi y}{2a}\right)}{sinh\left(\frac{\pi b}{2a}\right)}}cos\left({\displaystyle \frac{\pi x}{2a}}\right).$$

(45)

This problem was solved using the following nine transfinite element types:

• Classical B-spline transfinite element (21-27-33-113 nodes; see Figure 3).
• 12-node element (see Figure 5a).
• 18-node element (see Figure 5b).
• 25-node element (see Figure 5c).
• 27-node element (see Figure 8).
• 81-node element (tensor product; see Figure 10a).
• 80-node element (central node No. 41 is missing; see Figure 14).
• 76-node element (five boundary nodes are missing; see Figure 10b).
• 73-node element (five boundary nodes plus three internal nodes are missing; see Figure 10c).

For the sake of comparison, in addition to B-spline interpolation—which serves as the reference in the present paper—Lagrange and Bernstein polynomials have also been employed in selected cases.

The numerical results are summarized below:

8.1.1. Classical Transfinite Elements (21-27-33-113 Nodes)

These elements belong to the first class studied in the present paper. The overall results of the three models mentioned in Section 3 (each with 21, 27, 33, and 113 DOFs, shown in Figure 3) are presented in Table 1, which clearly illustrates the convergence trend.

For the sake of brevity, below we discuss only the case of the 113-DOF element (Figure 3d):

• The implementation of Model-1 (Section 3.1), based entirely on Lagrange polynomials for both blending and trial functions (of degree $p_{\xi }=16,p_{\eta }=12$), yields an excellent result: $L_2=1.9140\times {10}^{-5}\%$, using a $17\times 13$ Gaussian quadrature scheme. The same result was obtained using Bernstein polynomials.
• For Model-2 (Section 3.2), using $16\times 12$ integration cells (also referred to as natural B-spline elements), and applying a $4\times 4$ Gaussian quadrature per cell, the error was found to be $L_2=0.0324\%$.
• Finally, Model-3 (Section 3.3) requires $14\times 10$ integration cells (also referred to as Cox–de Boor spline elements), with $4\times 4$ Gauss points per cell, and yields an error ($L_2=6.2379\times {10}^{-5}\%$) of the same order of magnitude when using either Lagrange or Bernstein polynomials.

8.1.2. Unidirectional 12-18-25-32-Node Elements

These elements belong to the second class studied in the present paper (Figure 5). The results for all three element types are presented in Table 2. As the number of nodes increases (from 12 to 25), a monotonic convergence is observed across all types of polynomial and B-spline approximations.

The 32-node element (shown in Figure 5d) is derived from the 25-node element by subdividing the upper strip into two equal parts, thereby producing non-uniform blending functions. When Lagrange or Bernstein polynomials are employed, the six horizontal station positions are uniquely determined at $\eta =0,\frac{1}{4},\frac{2}{4},\frac{3}{4},\frac{7}{8},1$. In contrast, the use of B-splines as blending functions is not unique, since six control points can be generated from multiple combinations of breakpoints.

Nevertheless, despite the non-uniform placement of horizontal stations near the top (nodes 19 to 25), adopting the uniform knot vector $H=[0,0,0,0,\frac{1}{3},\frac{2}{3},1,1,1,1]$ yields an error of $L_2=0.0173\%$ (last column in Table 2), indicating convergence for the B-spline formulation as well.

Moreover, when the station at $\eta =\frac{7}{8}$ is defined with nodes 19 to 25 placed non-uniformly—e.g., following an arithmetic progression measured from the left at $\xi =\{0,\frac{1}{12},\frac{1}{5},\frac{7}{20},\frac{8}{15},\frac{3}{4},1\}$—the Lagrange and Bernstein formulations were affected only in the seventh and sixth decimal places, respectively. In contrast, using the knot vector $\Xi =[0,0,0,0,\frac{1}{10},\frac{3}{10},\frac{6}{10},1,1,1,1]$, the B-spline formulation showed a change in the third decimal place, yielding $L_2=0.0182\%$ instead of the previous $L_2=0.0173\%$.

8.1.3. Arbitrary-Noded 27-Node Transfinite Element

This element belongs to the third class examined in the present study (Section 5.1, Figure 8). The $L_2$ error norms of the numerical solution obtained from the three models applied to the 27-node element are presented in Table 3.

To demonstrate the effect of mesh non-uniformity, boundary nodes 15–18 and 9–13 were clustered toward the concealed vertex D (node 14) according to an arithmetic progression. The resulting difference in the error norm $L_2$ (%) was observed only beyond the fifth decimal place.

Similarly, regardless of the location of the boundary nodes, the tensor product of the internal nodes may be either uniform or non-uniform, for example, arranged as shifted Legendre roots or as Gauß–Lobatto–Legendre (GLL) points, as discussed later in Section 8.2.1.

8.1.4. T-Spline Elements

I. Proposed formulation

This element belongs to the fourth class examined in this paper (Section 6). We begin with a tensor-product T-spline element comprising 81 degrees of freedom (DOFs) arranged in a $9\times 9$ configuration. To investigate the impact of vertex reduction on solution accuracy, the number of control points is gradually decreased. The models considered include: one with a missing central node (80 DOFs, Figure 14), one with five missing boundary nodes (76 DOFs, Figure 10b), and another with five boundary and three internal nodes removed (73 DOFs, Figure 10c). In all cases, numerical integration was performed within the 36 cells (in $6\times 6$ setup) which are determined by the seven unique knots (breakpoints) per direction (at $\xi ,\eta =[0,1/6,2/6,3/6,4/6,5/6,1]$). In each case, the control points were selected such that the parameterization coincides with the physical coordinates, i.e., $x(\xi ,\eta )=\xi$ and $y(\xi ,\eta )=\eta$, and thus the determinant of the Jacobian was a constant ($det\mathbf{J}=1$). The corresponding results of the proposed model are presented in Table 4, which demonstrate a progressive increase in error as the number of DOFs decreases. In contrast, the condition number oscillates in the interval $[38.1,61.2]$.

II. Conventional T-spline

For completeness, the same T-index was employed for the conventional T-spline [17,18]. For the case of 73 nodes (Figure 10c), two matrices, Uvec and Vvec, each of size $73\times 5$, were constructed to represent the vertex topology. Using this local vertex information, the de Boor tensor-product B-spline basis function

$$T_k(\xi ,\eta )=N_{i,3}\left(\xi \right)N_{j,3}\left(\eta \right)$$

(46)

was evaluated over a uniform grid of points. The spatial distribution of the sum ${\sum }_{k=1}^{73}{T}_k(\xi ,\eta )$ is illustrated in Figure 17, revealing that the partition of unity (PU) property is not satisfied everywhere in the conventional T-spline formulation. In contrast, the proposed method inherently preserves this property.

To impose the PU property in the conventional T-spline, at each pair of parameters $(\xi ,\eta )$, Equation (46) is normalized as follows:

$${\tilde{T}}_i(\xi ,\eta )={\displaystyle \frac{w_i{T}_i(\xi ,\eta )}{{\sum }_{k=1}^{73}{w}_k{T}_k(\xi ,\eta )}},i=1,\dots ,73.$$

(47)

Obviously, although the denominator in Equation (47) theoretically includes 73 terms, in reality only a few of them are different than zero (local support property).

Remark: The comparison between the basis functions of this conventional formulation and the abovementioned proposed T-spline is as follows:

(1)

Before normalization: 59 out of the 73 basis functions (Equation (46)) are identical between the two formulations. The remaining 14 functions—specifically those numbered 18, 19, 20, 21, 26, 27, 28, 29, 34, 35, 36, 37, 46, and 53—exhibit differences within the interval $[-0.12,0.15]$.

(2)

After normalization: Considering the weights in Equation (47) equal to 1 (i.e., $w_i=1,i=1,\dots ,73$), 32 of the 73 basis functions are identical in both formulations. These include the basis functions associated with vertices numbered 1–17, 24, 25, 32, 33, 40, 41, 42, 48, 49, 56, 57, 58, 65, 66, and 67.

For a fair comparison between conventional T-splines (Equation (47)) and the proposed method, we employed the same control points, namely those of the latter (as explained in the first paragraph of this subsection). The convergence behavior of the conventional T-spline is summarized in Table 5. Similar to the proposed T-spline formulation, the results exhibit monotonic convergence up to a certain level, followed by a progressive increase in error as the number of DOFs decreases. In contrast, the condition number varies in the interval $[39.8,61.2]$.

Comparing Table 4 with Table 5, one may observe that—for the same control points—the proposed method is competitive to the conventional T-spline in terms of the error norm $L_2$, while also exhibiting a smaller condition number.

8.1.5. Coons-Patch Elements

This represents a type of element in addition to the four classes mentioned above. Although Coons-patch macroelements have been extensively discussed in numerous papers reviewed in the monograph [9], this model is included here for direct comparison with the models presented in Example 1. For the boundary-node configurations shown in Figure 5a–d (with 10, 13, 16, and 18 nodes, respectively), all illustrated in Figure 18, the computed results are presented in the left portion of Table 6.

For completeness, we examined both the standard linear blending functions, $E_1\left(\xi \right)=1-\xi$, $E_2\left(\xi \right)=\xi$ (first set), and the cosine-like blending functions, $E_1\left(\xi \right)=cos\left({\displaystyle \frac{\pi \xi }{2}}\right)$, $E_2\left(\xi \right)=1-cos\left({\displaystyle \frac{\pi \xi }{2}}\right)$ (second set). The results indicate that, with linear blending functions, the numerical solution converges to an erroneous value of approximately $L_2=2.3\%$, which is substantially larger than the competing values reported in Table 2. In contrast, the second set of blending functions—although heuristic in nature and inspired by the boundary conditions imposed on the top edge—leads rapidly to the exact solution.

Overall, the results of Example 1 suggest the following:

(1)

In a fully two-dimensional problem without any symmetry, the conventional Coons interpolation—implemented with linear blending functions—is insufficient for accurately representing the exact solution. Therefore, the inclusion of internal nodes is generally necessary to enhance accuracy.

(2)

These internal nodes may be arranged along horizontal and/or vertical stations and can follow a transfinite interpolation formula that is applied globally across the entire patch.

(3)

One practical approach is to place internal nodes along horizontal layers (stations), i.e., parallel to the $\xi$-axis. The station locations may correspond to either uniform or non-uniform $\eta$-values, and the associated blending functions will be constructed accordingly—based on uniform or non-uniform nodal points (breakpoints).

(4)

Regardless of whether the blending functions are uniform or non-uniform, the trial functions along each station may also be chosen to be uniform or non-uniform.

(5)

Once the closed-form expressions for the global bivariate shape functions have been derived using Lagrange polynomials, they can be readily extended to Bernstein polynomials and B-splines. While the local support property of B-splines affects the resulting bivariate basis functions, splines should be viewed as a specific choice of trial functions for univariate interpolation along each station. The same flexibility applies to the blending functions, which need not be restricted to cardinal types; in addition to Lagrange polynomials, Bernstein polynomials and B-splines are also valid options.

8.2. Example 2

Steady-State Conduction in a Cylindrical Wall (Half-Annulus): This example analyzes a long, hollow cylinder subjected to steady-state heat conduction, with a uniform inner surface temperature of $T_i={1000}^{\circ }$C and a uniform outer surface temperature of $T_o=0^{\circ }$C. The cylinder, defined by inner and outer radii $R_i=1$ and $R_o=32$, is assumed to be insulated to prevent axial heat flow (see Figure 19).

According to [28], the steady-state temperature distribution in the radial direction is given by:

$$T\left(r\right)=T_i+{\displaystyle \frac{(T_o-T_i)}{log\left(\frac{R_o}{R_i}\right)}}log\left({\displaystyle \frac{r}{R_i}}\right)$$

(48)

The radii are chosen to produce a sufficiently steep temperature gradient, making this example suitable for a convergence test. This test case was previously presented using single Coons macroelements based on cardinal natural cubic B-splines [2], formulated via truncated power series (i.e., not the procedure shown in Appendix A, but mathematically equivalent).

In the present study, we introduce several additional models and clarify all relevant implementation details. Unlike Example 1, where the Coons-patch element failed to converge accurately, here the Coons formulation does converge to the exact solution. Therefore, we adopt a different strategy: we begin with the Coons model in several variations (Section 8.2.1) and progressively introduce internal nodes (Section 8.2.2).

8.2.1. Coons-Patch Element

In this subsection, we investigate the convergence behavior of several models constructed using a single Coons-patch element, assuming the standard linear blending functions: $E_1\left(s\right)=1-s$, $E_2\left(s\right)=s$, where s denotes either of the parametric coordinates $\xi$ or $\eta$.

As reviewed in Ref. [9], previous studies have employed various cardinal trial functions, including:

(i) piecewise-linear functions,

(ii) cardinal natural cubic B-splines,

(iii) Lagrange polynomials.

In the present work, we extend this framework by also implementing non-uniform Lagrange polynomials and Cox–de Boor B-splines [29,30,31].

Although each edge of the Coons patch may consist of an arbitrary number of nodal or control points (each associated with degrees of freedom), for brevity of the presentation we assume equal number of points on opposite edges, as follows:

– The number of spans along the parallel edges $(AB,CD)$ is denoted $n_{\xi }$.

– The number of spans along the edges $(BC,AD)$ is denoted $n_{\eta }$.

For global Lagrange and Bernstein polynomials, the polynomial degrees in each parametric direction are $p=n_{\xi }$ and $q=n_{\eta }$, respectively. Interpolation using piecewise-linear or Lagrange polynomials requires $n_{\xi }+1$ and $n_{\eta }+1$ nodal points per side in the $\xi$ and $\eta$ directions, respectively. These nodal points correspond to the breakpoints used in the Cox–de Boor formulation. For cubic Cox–de Boor B-splines, the number of control points per edge exceeds the number of breakpoints by two (e.g., $n_{c,\xi }=n_{\xi }+3$) [26].

The parameterization of the patch is illustrated in Figure 19:

- Point A is the origin of the axis

- Arc $AB$ defines the $\xi$-axis

- Segment $AD$ defines the $\eta$-axis

- Edges $AB$ and $CD$ are subjected to Dirichlet boundary conditions

- Edges $BC$ and $AD$ are subjected to Neumann boundary conditions ($\partial T/\partial n=0$), due to symmetry.

In the numerical setup, each circular arc is uniformly subdivided into six breakpoint spans ($n_{\xi }=6$). The number of uniform breakpoint spans along the radial direction (e.g., edge $AD$) is varied incrementally: $n_{\eta }=1,2,3,\dots$

Model C1: Piecewise-Linear. The easiest case for implementing a Coons-patch macroelement is the adoption of piecewise-linear trial functions. This implies minimal cost in the estimation of these functions and also ensures local support. The vector of unknowns includes $2(n_{\xi }+n_{\eta })$ nodal values, $T_i$, all of them located on the boundary of the patch. Although this model is generally not equivalent to a set of $n_{\xi }\times n_{\eta }$ bilinear finite elements, the quality of the numerical solution is sometimes similar. In the particular case of Example 2, in which the solution does not depend on the polar angle $\theta$ (axisymmetric problem with no angular dependence), the numerical solution of the Coons-patch element with piecewise-linear trial functions coincides with that of the FEM solution using bilinear 4-node elements—provided the same mesh density is used. The result is shown in Figure 20 (blue line).

Model C2: Cardinal Natural Cubic B-Spline. This model briefly reproduces part of Ref. [2], in which cardinal natural cubic B-splines were previously used. Before the imposition of boundary conditions (BCs), the number of nodes is $n_{\mathrm{tot}}=2(n_{\xi }+n_{\eta })$. Since the $2(n_{\xi }+1)$ nodal points on edges $AB$ and $DC$ are restricted under Dirichlet-type BCs, the number of free degrees of freedom becomes $n_{\mathrm{free}}\_\mathrm{dofs}=2(n_{\eta }-1)$. The vector of unknowns includes only nodal values, $T_i$, which was the primary reason for adopting this model in the mid-1980s. Numerical integration is performed within the $n_{\xi }{n}_{\eta }$ cells created by nodal breakpoints in both directions. Given the piecewise polynomial degrees $p=q=3$, the number of Gauss points per integration cell is $6\times 6$ for curvilinear patches (and $4\times 4$ for rectangular ones). The results, shown in Figure 20 (red line), are identical whether the old-fashioned methodology of Ref. [2] (based on truncated powers) or the procedure of Appendix A is used—namely, the Curry-Schoenberg formulation [29], computed via the Cox-de Boor algorithm [30,31], with vanishing second derivatives at the ends of each edge.

Model C3: Truncated Power Series. This model extends the previous Model C2 on the same single Coons-patch element, but now additionally incorporates rotational degrees of freedom (DOFs) associated with the second derivatives at the ends of each edge—as illustrated, for example, in Figure A1c. In this case, although an edge such as $AB$ is discretized into $n_{\xi }$ breakpoint spans (yielding $n_{\xi }+1$ breakpoints), the number of associated control points along that edge becomes $n_{\xi }+3$. Consequently, after imposing Dirichlet boundary conditions, the final number of free DOFs becomes $n_{\mathrm{free}}\_{\mathrm{dofs}}^{\prime }=2(n_{\eta }+1)$. In this quantity, eight additional ’rotational’ DOFs are considered, though four of them—those along edges $AB$ and $DC$ under Dirichlet BCs—are eliminated. The vector of unknowns includes both nodal values, $T_i$, and the four curvatures $T{\left(\eta \right)}^{\prime \prime }$ at the ends of edges $AD$ and $BC$. The layout and number of integration cells remain the same as in Model C2, defined by the breakpoints. The corresponding results are presented in Figure 20 (green line).

Model C4: Cox–de Boor B-spline. This model continues to use the same Coons-patch element. Let ${\xi }_1=0,{\xi }_2,\dots ,{\xi }_{n_{\xi }+1}=1$ be the parameters associated with the breakpoints along edge $AB$. For a piecewise-cubic polynomial, we construct the well-known knot vector$\Xi =[0,0,0,0,{\xi }_1,\dots ,{\xi }_{n_{\xi }},1,1,1,1]$, which defines $n_{\xi }+3=n+2$ basis functions: $N_{1,3},\dots ,N_{n+2,3}$, where $n=n_{\xi }+1$ is the number of breaks along the edge, following the notation of Appendix A. The number of DOFs, both before and after applying boundary conditions, is the same as in Model C3, but here they consist of generalized coefficients ${\alpha }_i$ rather than nodal values. The integration cells are also the same, now with the additional property of local support. Despite this difference, the results are identical. In other words, although Model C3 lacks compact support (and deals with DOFs: $T_i$ and $T_{end}^{\prime \prime }$) while Model C4 possesses it (and deals with DOFs: $a_i$), the two functional sets are equivalent, as confirmed by spline theory and also sustained by numerous test cases. It should be noted that most literature discussing Coons patches in conjunction with B-splines and NURBS refers to the present Model C4, as it is based on the modern spline formulation introduced by Curry and Schoenberg (1966) [29] and implemented via the recursive Cox–de Boor algorithm (1972) [30,31]. Note that, since the linear blending functions used in Coons interpolation are identical across all formulations (Lagrange, Bernstein, and B-splines), the B-spline version of the Coons element follows directly.

Model C5: Uniform Lagrange polynomials. This model continues with Coons elements, using uniform global Lagrange polynomials as trial functions along each entire edge. Each bivariate function $N_i(\xi ,\eta )$ is continuous, influences the entire quadrilateral patch, and is associated with the nodal value $T_i$. Let p and q denote the polynomial degrees in the $\xi$- and $\eta$-directions, respectively. It can be readily verified that the integrand of each entry in the stiffness matrix, $k_{ij}={\int }_{\Omega }\nabla N_i\nabla N_j\det Jd\Omega$, is of degree $4p-1$ per direction—specifically, $2p$ from $\nabla N_i\nabla N_j$ and $2p-1$ from $\det J$. Consequently, for a rectangular patch $ABCD$, the stiffness matrix and error norm are estimated using a global scheme of $(p+1)\times (q+1)$ Gauss points, whereas curvilinear patches such as that in Figure 19 require $2p\times 2q$ Gauss points spanning the entire patch. Due to this, the terms “patch” and “Coons element” are used interchangeably. Closely related terminology such as “macroelement” or “Coons-patch element” has also appeared in previous literature [9]. From the convergence diagram in Figure 20 (orange color), it becomes evident that increasing the polynomial degree p in the $\eta$-direction steadily improves numerical accuracy up to $p_{\mathrm{cr}}=21$, beyond which the solution becomes unstable.

Model C6: Non-uniform Lagrange polynomials. The sixth model studies Coons elements, using non-uniform global Lagrange polynomials as trial functions along each entire edge. Following the practice of spectral methods by Karniadakis and Sherwin [32], one of the best choices is the set of Gauß–Lobatto–Legendre (GLL) points, which are defined as follows:

$${\xi }_i^{\mathrm{GLL}}=\left\{\begin{array}{cc}-1\hfill & \mathrm{for}i=1\hfill \\ {\widehat{\xi }}_i\hfill & \mathrm{for}i=2,3,\dots ,p\hfill \\ +1\hfill & \mathrm{for}i=p+1,\hfill \end{array}\right.$$

(49)

where ${\widehat{\xi }}_i$ denotes the roots of the Lobatto polynomial $L{o}_{p-1}\left(\xi \right)$ of order $p-1$, which is defined as the first derivative of the Legendre polynomial $L_p\left(\xi \right)$ of order p:

$$L{o}_{p-1}\left(\xi \right)={\displaystyle \frac{d{L}_p\left(\xi \right)}{d\xi }}.$$

(50)

Practically, for any given degree p, the $(p+1)$ GLL points ${\xi }_i^{\mathrm{GLL}}\in [-1,1]$, defined in Equation (49), can be numerically determined by setting i = p - 1 in the following MATLAB command:

   roots = vpasolve((legendreP(i,x) - legendreP(i+2,x)) == 0);

Based on the above nodal points (images of the GLL points), we can easily construct non-uniform Lagrange polynomials and use them as trial functions for each edge of the quadrilateral patch $ABCD$.

In our case, we kept the uniform subdivision of the circular arcs at $n_{\xi }=6$, and adopted GLL-based nodal points along the edges $BC$ and $AD$. The polynomial degree varied from $p_{\min }=1$ to $p_{\max }=99$. It was found that, up to degree $p_{\mathrm{cr}}=21$, the results are practically the same in both formulations—i.e., the uniform and non-uniform (GLL-based). By further increasing the polynomial degree up to $p_{\max }=99$, the non-uniform formulation yielded a monotonically decreasing error norm, converging to the value $L_2=0.0102053\%$ (Figure 19). This remaining error is attributed to the discretization of the circular arcs, which are based on seven nodal points per arc (i.e., six uniform nodal spans: $n_{\xi }=6$).

Remark: Let us divide the interval $\xi \in [0,1]$ into n segments (i.e., $n+1$ nodal points) in the $\xi$-direction. Regarding global interpolation, one possibility is to introduce $(n+1)$ Lagrange polynomials $L_i\left(\xi \right)$ of degree $p=n$, which leads to integrands in the stiffness and mass matrix entries

$$k_{ij}=\int L_i^{\prime }\left(\xi \right)L_j^{\prime }\left(\xi \right)d\xi ,m_{ij}=\int L_i\left(\xi \right)L_j\left(\xi \right)d\xi$$

of degree $p_{\mathrm{integrand}}=2n$; therefore, the required number of Gauss points for the whole interval is $n_g^{\mathrm{Lagrange}}=n+1$.

On the other hand, cubic B-spline interpolation over the aforementioned n elements leads to matrix element integrands of degree six, and thus requires four Gauss points per element; this results in a total of $n_g^{\mathrm{B}}-\mathrm{spline}=4n$ Gauss points.

Of course, the stiffness and mass matrices in the Lagrange formulation will be of size $(n+1)\times (n+1)$, whereas in the B-spline formulation they will be of size $(n+3)\times (n+3)$ (because $n+1$ breakpoints give $n+3$ control points). In any case, the B-spline formulation requires more Gauss points than the Lagrange formulation.

Model C7: Finite Element Method (FEM). Based on the pair $(n_{\xi },n_{\eta })$, the FEM model using $\left(n_{\xi }{n}_{\eta }\right)$ bilinear 4-node elements consists of $(n_{\xi }+1)\times (n_{\eta }+1)$ nodal points. In the general case in which the solution does not follow the Coons interpolation (like Example-1), the model of piecewise-linear Coons-patch element—with $2(n_{\xi }+n_{\eta })$ nodes—differs from the FEM solution. In other cases such as in Example-2, the piecewise-linear based Coons element is equivalent to the FEM model, and thus both models have the same numerical error $L_2$.

Overall, the results shown in Figure 20 suggest:

(1)

All models converge to different values. This fact is mainly due to the incapability of accurately representing the circular arc (for $n_{\xi }=6$). Below we start from the less accurate model and end with the most accurate.

(2)

Piecewise-linear Coons interpolation coincides with the FEM solution. This is because Example-2 is an axisymmetric problem with no angular dependence.

(3)

The cardinal natural cubic B-spline is characterized by a smaller error than the above piecewise-linear model and FEM.

(4)

Uniform Lagrange polynomials lead to a smaller error but, after the value $n_{\eta }=21$, diverge.

(5)

De Boor cubic B-splines closely follow the accuracy of uniform Lagrange polynomials up to $n_{\eta }=5$, then become less accurate until $n_{\eta }\approx 25$, and eventually converge to the value $L_2=0.0160\%$.

(6)

The non-uniform Lagrange polynomials based on GLL points rapidly converge to the accurate solution. A small error ($L_2=0.0102\%$) remains, due to the incapability of accurately representing the circular arc using $n_{\xi }=6$ spans.

In all the above models, it is anticipated that the error norm $L_2$ will be decreased when the circular arcs are represented through more nodal points or NURBS.

8.2.2. Transfinite Elements

A. Classical Transfinite Elements

Each of the four classical transfinite elements shown in Figure 3 is mapped onto the entire half-annulus depicted in Figure 19. Given the parametric coordinates $\xi$ and $\eta$ of the nodal points, the corresponding polar coordinates r and $\theta$ in the physical domain are obtained as:

$$\begin{array}{cc}\hfill r& =R_i+\eta \left(R_o-R_i\right),\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill \theta & =(1-\xi )\pi .\hfill \end{array}$$

(52)

The numerical results are reported in Table 7, where convergence behavior can be observed. As in Example 1, the same result was obtained when using Bernstein polynomials.

Moreover, it is instructive to compare the above results with those obtained for the previously mentioned Model C5 (uniform Lagrange or Bernstein polynomials for a single Coons element) under similar conditions, i.e., with $n_{\eta }=4,5,6,12$ subdivisions along the radial direction.

From Figure 20 (orange curve), the corresponding errors are 8.5851%, 4.3559%, 2.3578%, and 0.1127%, respectively.

It is evident that, in Example 2, the aforementioned accuracy of the Coons model is very close to that of the classical transfinite elements (Table 7), primarily due to the independence of the solution from the polar angle $\theta$.

To some extent, the small differences can be attributed to the different discretization of the semi-circular edges in the transfinite patch. It is likely that, if NURBS were employed, these differences would be further reduced.

The aforementioned similarity between classical transfinite elements and the Coons element is observed for all models discussed in Section 8.2.1. For instance, both the 113-node transfinite element (with 56 DOFs on the boundary) and its counterpart Coons element (model C4, based on Cox–de Boor cubic B-splines with a total of 56 DOFs) yield the same error norm, $L_2=0.2777\%$. This value is approximately 2.5 times larger than the error obtained when using Lagrange polynomials as trial functions (see $L_2=0.1127\%$ in Table 7).

B. T-spline Elements

The half annulus shown in Figure 19 was discretized using a uniform tensor-product of $7\times 7$ breakpoints, resulting in a grid of $9\times 9$ control points (81 DOFs), as illustrated in Figure 21a. More specifically, a tensor-product of Cox–de Boor cubic B-splines ($p=3$) was constructed using a uniform knot vector

$$\Xi =[0,0,0,0,1/6,2/6,3/6,4/6,5/6,1,1,1,1].$$

The parametric space $(\xi ,\eta )\in [0,1]\times [0,1]$ was uniformly divided into a grid of $9\times 9$ test points $({\xi }_i,{\eta }_i)$, for $i=1,\dots ,81$, and a corresponding uniform mesh of test points ${\mathbf{x}}_{\mathrm{test}}$ was generated in the physical domain.

The system of 81 equations for $(k,l=1,\dots ,9)$:

$${\mathbf{x}}_{\mathrm{test}}({\xi }_k,{\eta }_l)=\sum _{i=1}^9\sum _{j=1}^9{N}_{i,p}\left({\xi }_k\right)N_{j,p}\left({\eta }_l\right){\mathbf{x}}_c^{ij},k=1,\dots ,9;l=1,\dots ,9,$$

(53)

determines the positions of the 81 control points ${\mathbf{x}}_c$, illustrated by red crosses in Figure 21a. In this way, the B-spline curve (with unit weights), corresponding to each circular arc, passes through nine uniformly distributed points along the arc. Consequently, a small deviation from the exact representation of the arc occurs between these points. Extending this formulation to NURBS is a straightforward task, but it is left to the interested reader.

It was found that the above B-spline tensor product yields the same numerical result, i.e., $L_2=1.0494\%$, as that obtained using the corresponding Coons-patch element with the same number of breakpoint spans (see Section 8.2.1, Model C4). As previously mentioned, this outcome arises because the problem is axisymmetric and exhibits no angular dependence.

Starting from the abovementined B-spline tensor product with 81 control points, we proceed as follows. First, we remove the central node, resulting in a T-mesh with 80 nodes (see Figure 14 in the parameter space, and Figure 21b in the physical space). Second, we retain the central node but eliminate five boundary nodes, yielding a model with 76 nodes (see Figure 10b in the reference square, and Figure 21c in the physical domain). Third, we further remove three internal nodes, ultimately producing a T-mesh with 73 nodes (see Figure 10c in parametric space and Figure 21d for the physical annular domain). For all these cases, the results are presented in Table 8, where it can be observed that as the number of removed nodes increases, the error also increases.

For completeness, Table 9 also presents the numerical solution based on the conventional normalized T-spline formulation [18] (Equation (47)). It can be observed that as the number of removed nodes increases, the error also increases.

8.3. Example 3

Eigenvalues of acoustic cavity: This example deals with the natural frequencies of a rectangular acoustic cavity of size $a\times b=2.5\times 1.1$, under Neumann boundary conditions ($\partial p/\partial n=0$). The governing wave equation is

$${\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^{2}p}{\partial t^2}}-{\nabla }^{2}p=0,$$

(54)

where p denotes the acoustic pressure. According to Courant and Hilbert [33], the exact eigenvalues are given by:

$${\lambda }_{mn}={\omega }_{mn}^2={\pi }^2{c}^2\left({\displaystyle \frac{m^2}{a^2}}+{\displaystyle \frac{n^2}{b^2}}\right),m,n=0,1,\dots ,\infty$$

(55)

The application of the Galerkin procedure to Equation (54) for the eigenvalue problem, leads to the following matrix formulation:

$$\left[\mathbf{M}\right]\left\{\ddot{\mathbf{p}}\right\}+\left[\mathbf{K}\right]\left\{\mathbf{p}\right\}=0,$$

(56)

where the entries of mass matrix ($\left[\mathbf{M}\right]$) and stiffness matrix ($\left[\mathbf{K}\right]$) are given by:

$$m_{ij}={\displaystyle \frac{1}{c^2}}{\int }_{\Omega }{N}_i{N}_{j}d\Omega ,k_{ij}={\int }_{\Omega }\nabla N_i·\nabla N_{j}d\Omega$$

(57)

The computational mesh of 73 nodes was that of Figure 10c, adapted from unit square to a rectangle of dimensions $2.5\times 1.1$. The first eigenvalue (Mode-1) was found equal to zero (modes $m=n=0$; ${\lambda }_{00}=0$) as should be, whereas for the next modes the error (in percent) was calculated by the formula:

$$E_{mn}={\displaystyle \frac{{\lambda }_{mn}^{\mathrm{calculated}}-{\lambda }_{mn}^{\mathrm{exact}}}{{\lambda }_{mn}^{\mathrm{exact}}}}\times 100(\%)$$

(58)

Using the proposed formulation of the T-spline (according to Section 6) and also that of the conventional T-spline (Equation (47))–both of them associated with the 73-node element that is illustrated in Figure 10c–the errors of calculated eigenvalues are shown in Table 10. One may observe that the proposed T-spline outperforms.

For the sake of completeness, in addition we present the results obtained using the B-spline tensor-product of 81 control points (Figure 10a) for degree $p=3$, and also the assembly of 64 uniform bilinear elements (81 nodes). The results are shown in Table 11, where one may observe that the tensor product B-spline model of 81 nodes is slightly better than both T-spline models of 73 nodes (with the exception of Mode-3 where it outperforms). In addition, all B-spline and T-spline models substantially outperformed compared with the conventional bilinear FEM model of 81 nodes.

8.4. Example 4

Patch test in plane elasticity: Among several patch tests [34], we present here the case of a unit square under plane stress, with Dirichlet boundary conditions $u(\xi ,\eta )=0$, $v(\xi ,\eta )=\eta$, where $0\le \xi ,\eta \le 1$. For the 12-node element (Figure 5a), the implementation of Lagrange polynomials (Equation (19)) or Bernstein polynomials (Equation (20)) yields exactly the theoretical values at the internal nodes 5 and 6 (i.e., horizontal displacements $u=0$, and vertical displacements $v=0.5$). Moreover, for elastic modulus $E=1$, the stress throughout the element coincides exactly with the theoretical values, namely

$${\sigma }_{xx}={\displaystyle \frac{\nu E}{1-{\nu }^2}}=0.3297,{\sigma }_{yy}={\displaystyle \frac{E}{1-{\nu }^2}}=1.0989,{\tau }_{xy}=0.$$

Similarly, we tested the case where node 2 (Figure 5a) was fully supported, while nodes 1 and 3 were constrained to roll horizontally. The top edge 8–9–10–11–12 was uniformly loaded in tension, and the vertical edges were free-free. In this setting, both the Lagrange (Equation (19)) and Bernstein (Equation (20)) polynomials yielded the theoretical values, i.e., ${\sigma }_{xx}=0$, ${\sigma }_{yy}=1$, and ${\tau }_{xy}=0$.

Regarding the implementation of B-splines on the 12-node element (Figure 5a), the same accurate result as above was obtained, despite its hybrid nature. One illustrative example is the following. The blending function along the vertical ($\eta$) direction was taken as a Bernstein polynomial of degree $p_y=2$, since a cubic B-spline is not applicable. The trial functions were chosen as follows:

(1)

On the bottom layer (consisting of 3 nodes: 1–2–3), a non-cardinal Bernstein polynomial of degree $p_x=2$ was used.

(2)

On the middle layer (consisting of 4 nodes: 4–5–6–7), a Bernstein polynomial of degree $p_x^{\prime }=3$ was employed. Alternatively, one could adopt a quadratic B-spline ($p_x^{\prime \prime }=2$) with knot vector $\Xi =[0,0,0,\frac{1}{2},1,1,1]$.

(3)

On the top layer (consisting of 5 nodes: 8–9–10–11–12), a cubic B-spline ($p_x=3$) was used with knot vector ${\Xi }^{\prime }=[0,0,0,0,\frac{1}{2},1,1,1,1]$.

Further results will be reported in a future paper, dedicated to elasticity problems.

9. Extension to Collocation Method

Retaining the bivariate global shape (or basis) functions of the transfinite elements discussed in the previous sections, the Galerkin method can be readily replaced by the Collocation Method. This alternative has been explored in a series of studies since 2007 (see [9] for an overview, including IGA-based papers since 2010), with representative applications to potential and elasticity problems presented in [35] and [36], respectively.

In the context of transfinite elements using Lagrange polynomials for both blending and univariate trial functions, collocation points can be selected as either Gauss or Chebyshev points, defined globally over the entire patch. This formulation offers the flexibility to choose a number of collocation points equal to the number of unknowns. Notably, it has been demonstrated that the computed eigenvalues remain unchanged whether a uniform nodal mesh is used in conjunction with Gauss/Chebyshev collocation points, or a non-uniform mesh is constructed by placing nodes at the mapped roots of Legendre or Chebyshev polynomials and applying nodal collocation directly.

In contrast, B-spline-based collocation requires a standard elementwise distribution of collocation points within each B-spline element across the patch. In general, this results in a number of collocation points exceeding the number of unknowns, thereby necessitating a least-squares solution procedure. Alternatively, following the pioneering work of de Boor [26], Greville or Demko abscissae can be globally selected as collocation points across the entire patch.

An additional challenge inherent to the collocation method is the imposition of mixed boundary conditions, where segments of the boundary are governed by Dirichlet conditions while others are subject to Neumann conditions. This complicates the consistent enforcement of boundary constraints. The optimal handling of this issue—along with the assessment of the performance of alternative basis functions such as trigonometric or Hermite B-splines [37,38,39]—remains an open area for future investigation.

10. Discussion

This paper demonstrates that transfinite elements are not limited to Lagrange polynomials but can also be formulated using B-splines, which may serve as both blending and trial functions. In general, all three formulations—Lagrange polynomials, Bernstein polynomials, and B-splines—yield basis functions with identical closed-form expressions. Consequently, the term macroelement has been adopted to describe the unified framework represented by a collection of sixty papers documented in Ref. [9]. Within this framework, the entire patch can be treated as a single macroelement, regardless of the specific type of blending or trial functions employed. This paper supports the conjecture that the classical blending functions used in Coons–Gordon patches can be replaced by Bernstein polynomials or B-splines, thereby enabling the seamless transfer of established techniques into the isogeometric analysis (IGA) domain.

The use of B-splines is inherently characterized by their local support, meaning that within each integration cell—defined between adjacent breakpoints—only a limited number of basis functions are non-zero. However, in this paper, we did not focus on bandwidth optimization; rather, our primary objective was to establish a proof of concept demonstrating that B-splines are fully applicable within the framework of transfinite elements. This applicability was verified across at least four distinct classes, ranging from classical transfinite elements, such as those of Gordon and Hall [2], to more advanced constructs like T-splines. For completeness, the Coons–patch macroelement was studied as well.

The outcome of the present study is the culmination of a long-term effort and accumulated experience that dates back to mid-1984. At that time, the author was co-supervising a PhD thesis—which was ultimately never defended—and subsequently a Master’s thesis in 1988, both conducted within the CAD/CAE group of the Mechanical Engineering Department at the National Technical University of Athens. Further historical context is provided in Ref. [9]. The group’s first publication appeared in 1989 [2]; however, due to various circumstances, the second publication—part of the collection mentioned in the first paragraph of this section—was delayed until 1999 (see Ref. [9]).

Motivated by concerns related to the Runge phenomenon, the group initially developed a Coons macroelement based on cardinal natural cubic B-splines defined along each edge of the Coons patch (explained in Appendix A). Considerable time and deliberation were required before adopting Lagrange polynomials, which ultimately proved to offer excellent performance, even when constructed using uniformly spaced nodal points. Based on our experience, we propose that potential and elasticity problems can be effectively addressed using Gordon’s transfinite interpolation, with polynomial degrees of up to 8 and 10, respectively. Today, we also advocate for the use of spectral methods, which mitigate the excessive end-point oscillations often encountered in high-degree polynomial approximations (see Ref. [40] for further details, and also refer to Model C6 around Equation (49)).

The deeper motivation behind this research stems from our conviction that, since isogeometric analysis (IGA) operates over entire patches, it is fundamentally aligned with the general principles of transfinite interpolation. This perspective suggests that IGA is not an entirely distinct methodology, but rather a natural evolution or extension of transfinite concepts. Preliminary reflections on this connection were outlined in our earlier monograph (see Ref. [9], among others), while a more recent review paper is due to Provatidis et al. [41].

An early attempt to adapt transfinite interpolation for the accurate representation of conic sections and quadrics was presented in Ref. [42], where weights were determined via nonlinear optimization. In that study, it was—almost intuitively—found that a direct substitution of Lagrange polynomials with their Bernstein counterparts produced nearly identical and highly accurate numerical results. This observation further supported the idea that transfinite formulations are inherently compatible with the foundational principles of isogeometric analysis (IGA).

Since then, significant progress has been made:

First, it was proven that the substitution of Lagrange polynomials with Bernstein polynomials of the same degree is fully compatible in both Coons elements and classical transfinite elements, such as the 113-node element shown in Figure 3. Specifically, the existence of a transformation matrix between the two sets of basis functions was revealed [20].

Second, it was shown that elements with structured interior nodes but irregular boundary node distributions do not admit the aforementioned transformation matrix. Nevertheless, replacing the Lagrange polynomials with Bernstein polynomials of the same degree yields slightly improved results. This observation was explained by extending the formulation of the transfinite interpolation formula with respect to the choice of blending functions and trial functions [16].

Third, a mechanism was identified that enables the application of the transfinite interpolation formula to T-meshes, but only with Lagrange [21] and Bernstein [20] polynomials. Notably, it was revealed that this can be achieved by enforcing constraints both through the transfinite interpolation formula and with the aid of a tensor product defined on the background mesh.

The present paper completes the aforementioned efforts by demonstrating the use of B-spline interpolation both as blending functions and as trial functions, treating the entire patch as a single transfinite element. This approach is illustrated across four classes of elements: classical transfinite elements, one-directional elements with nodes arranged in parallel layers, elements with structured interiors but irregular boundaries, and, finally, T-splines.

We now refer to these four classes:

First class: The classical 21-, 27-, 33-, and 113-node transfinite elements were successfully treated using B-spline interpolation, together with integration cells automatically determined from the mesh density along the horizontal and vertical stations containing the nodal points.

Second class: The unidirectional transfinite element is also handled automatically by considering a single subdivision in the vertical direction, while in the horizontal direction—along which the layers of nodes are aligned—the union of knot spans is readily determined using the MATLAB function unique.

Third class: The transfinite element with a structured interior arranged in a tensor-product pattern but with arbitrarily distributed boundary nodes extends the previous approach by considering stations in both directions. The advantage of using the MATLAB function unique to identify knot spans for numerical integration of the stiffness matrix remains applicable.

Fourth class: Regarding the class of transfinite elements known as T-splines—possibly first considered in Ref. [2] using Hermite splines—prior studies have addressed the enforcement of constraints at the missing nodes within the T-mesh. Within this framework, Lagrange polynomials can be systematically treated on each segment containing a missing node, with an analogous procedure applicable to Bernstein polynomials. Notably, in both cases, the constraint formulation—linking the state of an incomplete station to that of a complete station—remains identical [16,20].

The present study extends this approach by employing concepts from computer-aided geometric design (CAGD), specifically one-dimensional knot insertion theory, to impose analogous constraints. For boundaries featuring missing nodes, the transformation matrix is readily derived (cf. Appendix B), given the independence of each edge from the remainder of the patch; consequently, knot insertion on the curve representing the edge is well-defined and effective.

In this work, each station is treated as an individual curve with associated data, which leads to the assumption that the same type of transformation matrix applies uniformly to all horizontal and vertical stations (cf. Appendix B). Furthermore, at vertices corresponding to the intersection of two stations, the mean average of the contributing stations was employed. These assumptions, previously validated for both Lagrange and Bernstein polynomials [20,21], form the foundation of the present methodology.

The advantage of the proposed transfinite element—compared to the conventional T-spline method [18,19] (nowadays of great interest, e.g., [43,44])—is that it a priori satisfies the partition of unity property, thus eliminating the need for normalization of the resulting basis functions. In fact, as noted at the end of Section 8.1, the proposed methodology produces basis functions that are closer to the raw (unnormalized) set than to the normalized one. Nevertheless, the error norm obtained using the proposed method for the solution of Laplace equation is of the same order of magnitude with the conventional normalized T-spline, but the numerical solution is characterized by a smaller condition number. In the solution of the eigenvalue problem, the proposed T-spline outperforms. A potential disadvantage, however, is that the proposed method does not a priori guarantee the initial parameterization of object’s shape under analysis.

In all computations reported here, when the patch is a rectangle of size $L_x\times L_y$, the control points were chosen so that at the vertices of the parametric unit square, the physical coordinates satisfy $x(\xi ,\eta )=L_x\xi$ and $y(\xi ,\eta )=L_y\eta$, which leads to a constant Jacobian everywhere. A similar technique—using the mapping between uniform test points on both the reference unit square and the true domain—was applied for the half-annulus of Example 2 (Section 8.2). Therefore, the implemented computer codes perform equally well for curvilinear patches as well.

The treatment of circular (or possibly elliptic) parts is better accomplished using NURBS with predefined control points and associated weights, rather than the inverse engineering approach applied here for estimating control points that accurately represent these curves only at distinct uniform positions. However, this is a straightforward task that may be carried out by the interested reader to observe firsthand the improvement resulting from this modification.

Although the present paper provides self-contained MATLAB codes, it is not difficult to incorporate the proposed formulations—Bernstein, B-spline, and the NURBS version—into popular academic software such as IGAFEM [45] and GeoPDEs [46].

The proposed procedure leads to the construction of basis functions that can be directly applied to the estimation of mass and stiffness matrices for the numerical solution of potential and elasticity problems in the continuum. Both static, dynamic, and coupled problems (e.g., thermoelasticity, fluid–structure interaction) can be solved by following standard FEM/IGA algorithms. A common feature of these problems is that only $C^0$ continuity is required between the proposed elements in an assembly. In contrast, for shells and plates requiring $C^1$ continuity, the formulation presented in this paper must be extended. A brief summary of past research, along with ongoing work, is presented in the following three paragraphs.

The original report by Coons [2] introduced two interpolation formulas for the function $w(\xi ,\eta )$ over a curvilinear patch. The first formula is based on prescribed values of w along the four boundaries of the patch; this case was thoroughly discussed in the present paper. The second formula involves prescribed values of w and its partial derivatives $\partial w/\partial n$ on the four edges of the patch. Based on this second interpolation formula, Provatidis and Angelidis [47] constructed a rectangular Coons-patch macroelement (with nodes located only on the boundary) suitable for thin plate-bending analysis. Following Coons [2], the four blending functions per direction were chosen as cubic Hermite polynomials (two cardinal and two non-cardinal, as in beam bending), while the trial functions were Hermite polynomials of degree $(2n-1)$, with n denoting the number of nodes along an edge. This element was tested in six examples with simply supported or clamped edges, subjected to concentrated or distributed loads. Despite using only boundary information, in most cases the deflection error was reported to be less than 0.2%. Furthermore, it was shown that the lowest-order version of this transfinite element—using cubic Hermite polynomials as both blending and trial functions with two nodes per edge $(n=2)$—reduces to the well-known Bogner–Fox–Schmit (BFS) element [48].

In a later study ([9]), it was demonstrated that by selecting the same degree of Hermite polynomials for both blending and trial functions, transfinite interpolation degenerates to conventional tensor-product Hermitian interpolation. In this case, the results for both static and eigenvalue problems were found to be exceptionally accurate. By analogy to the present paper, where the treatment of internal nodes deviating from the tensor-product structure has been addressed, Hermite-based transfinite plate-bending elements with slightly unstructured interiors (such as the one shown in Figure 1e) can be readily constructed in a similar manner, although no published report currently exists. Despite this flexibility, Hermite interpolation cannot naturally accommodate non-rectangular plate elements, which motivates the use of a different approach, as discussed in the next paragraph.

To overcome the limitations of Hermite polynomials and to accommodate smooth curvilinear patches such as circular plates, the use of B-splines and NURBS has been proposed [9]. The central idea is to express the deflection as a tensor product of B-splines (or NURBS), namely

$$w(\xi ,\eta )=\sum _{i=0}^n\sum _{j=0}^m{N}_{i,n}\left(\xi \right)N_{j,m}\left(\eta \right)a_{ij},$$

where $a_{ij}$ denote generalized coefficients. This formulation has been shown to provide excellent accuracy for both static and dynamic problems. Nevertheless, within this framework, the treatment of non-tensor-product B-spline/NURBS elements using non-cardinal blending functions (either B-splines or NURBS) remains an active area of research. From an industrial perspective, the important aspects of the proposed formulation that favor its practical adoption are: (i) the partition of unity property without the need for normalization, (ii) the capability to handle unstructured grids, and (iii) the straightforward joining of patches. Moreover, after the analytical determination of the basis functions, the proposed method was found to be ten to twenty times faster than the conventional T-spline; however, this improvement may be related to the choice of B-spline subroutines (e.g., MATLAB function spcol) or to specific programming implementations. These features have recently attracted considerable interest from the research community (e.g., [49,50,51]).

This study is also connected to recent advances in multi-material and isogeometric-inspired topology optimization, where some remarkable contributions include [52,53,54,55,56].

In summary, the present study advances the theory and implementation of transfinite interpolation by integrating B-spline and NURBS-based formulations within a unified framework that accommodates a wide range of element classes, including classical transfinite elements, structured and irregular patches, and T-splines. By leveraging concepts from CAGD and knot insertion theory, the methodology ensures compatibility with isogeometric principles while preserving desirable properties such as partition of unity. Although the initial parameterization is not inherently guaranteed, practical strategies have been demonstrated to achieve consistent mappings and accurate geometric representations. Future work may focus on extending this framework to three-dimensional patches, refining control point estimation for complex geometries, and integrating the proposed formulations into broader isogeometric analysis pipelines. The flexibility and modularity of the approach make it a promising candidate for further exploration in both academic and industrial settings.

11. Conclusions

This work demonstrates that classical transfinite finite elements can be reformulated using Bernstein polynomials and B-splines, enabling a unified framework for both structured and partially unstructured grids, including unidirectional layouts and T-meshes. By replacing Lagrange polynomials with Bernstein or B-spline bases and expressing transfinite interpolation projectors in B-spline form, the method allows independent control of horizontal and vertical stations and facilitates the imposition of linear constraints along patch boundaries and internal mesh-lines using established knot insertion techniques. Two complementary strategies—transfinite interpolation and tensor-product background grids—lead to automatic construction of blending functions that satisfy the Partition of Unity Property without normalization. The approach further simplifies the system by eliminating secondary degrees of freedom in favor of primary ones, resulting in a flexible and computationally efficient formulation. Convergence studies in two examples of potential problems, eigenvalue acoustic analysis in a third example, and a patch test in plane elasticity have shown that the proposed approach performs robustly for both uniform and non-uniform grids. The proposed formulation can be extended to plate-bending and shell analysis as well as to three-dimensional problems.