Transfinite Patches for Isogeometric Analysis

1. Introduction

Classical textbooks in computer-aided geometric design (CAGD) [1, 2], suggest that one of the first CAGD formulations in Computational Geometry is the Coons’ interpolation formula [3]. Despite its high historical importance, we could say that the latter is merely the extension of the linear interpolation in two-dimensional quadrilateral domains (patches). Even in the first Coons’ report [3], it is clearly written that this interpolation refers to both the geometry of a patch as well as to the estimation of a function U(x,y) when its values are prescribed along the four edges of a quadrilateral patch ABCD.

The disadvantage of Coons’ interpolation formula is its incapability of considering information related to inner points. For example, if the values of the function U(x,y) vanish on the entire boundary, then Coons’ interpolation formula will provide zero values everywhere in the domain ABCD. This shortcoming was resolved by William Gordon at General Motors Research Laboratories (GMRL), who proposed that information may be given of course along the boundary (as also happened with Coons’ interpolation) but also along internal lines parallel to the edges of the quadrilateral patch (sometimes called inter-boundaries or stations). In more detail, Gordon preserved the meaning of blending functions, which had been previously proposed by Coons, but extended it to include the inter-boundaries as well [4-7].

It is worth mentioning that the same team at GMRL applied Coons’ interpolation for the implementation of the Ritz-Galerkin method to the numerical solution of a boundary value problem (BVP) regarding Laplace equation within a complex domain in the form of a five-sided patch [8]. This fact can be rephrased as follows. The utilization of the same global interpolation (within the entire patch) for both the geometry and the numerical approximation of the solution were proposed for the first time by the GMRL group, for practical engineering purposes. This is simply the isoparametric consideration, well known in finite element analysis, in which the global trial function was taken equal to that which describes the geometry of the geometric patch.

The few papers published by the abovementioned GMRL group were written in a very dense style, so their decryption was not fully performed widely. The impression a reader takes from their reading is that the blending functions must be Lagrange polynomials which vanish at the internal stations, whereas as the trial functions along them may be either Lagrange polynomials or Hermite splines. The latest papers published by the GMRL group on this issue are [9,10].

For over thirty years, the above idea of using the same CAGD-based interpolation for both the geometry and the physical function U(x,y), which is involved into the governing PDE D(u) – f = 0, did not attract much attention [11-13], until the advent of the so-called isogeometric analysis (IGA), which changed the establishment [14,15]. In [15](p.8), it is recognized that the same interpolation was used in Coons and transfinite elements as well, but IGA dominated because it closely follows the most recent rules of CAGD-theory in which the shape remains unaltered after knot insertion, degree elevation, or any combination of them. The history of implementing transfinite interpolation, as an alternative to FEM, includes a volume of 60 papers written from 1989 till 2018, of which the resumé has been included in a monograph published in 2019 by Provatidis [16].

Although there is no doubt that IGA is an extraordinary technique itself, this paper shows that it can be considered as a particular case of transfinite interpolation. Clearly, transfinite interpolation is a superset of tensor-product Lagrange polynomials, tensor-product Bernstein polynomials, tensor-product B-splines, and tensor-product NURBS. The novel key point to this consideration is to revise the meaning of the so-called blending functions including non-cardinal sets as well, and thus to obtain hybrid schemes, with the same or different functional sets within the transfinite element which describes an entire quadrilateral patch. Preliminary but rather empirical thoughts on this issue have been recently presented in Ref. [17]. In the present paper a new conjecture is proposed.

The paper is structured as follows. In section 2 we begin with the classical linear blending functions which are involved in Coons’ interpolation formula and discuss the meaning of the three projectors involved in the Boolean sum. In section 3 we continue with higher-order blending functions which are involved in Gordon’s transfinite interpolation formula, explain the meaning of the stations, and prove that the tensor product is a special case. In section 4 we present the most useful univariate functions, which may play the role of either a blending or a trial function. In section 5 we show that the non-rational tensor-product Bernstein polynomials is an equivalent tensor-product Lagrange polynomials, and we provide a formula for the determination of the α-coefficients therein. In section 6 we show that the tensor-product B-splines are produced using B-splines for both the blending and trial functions. Based on this finding, we further conjecture the utilization for blending functions, of any linear independent and complete functional set, including non-cardinal (de Boor) B-splines and NURBS. In section 7 we present variations of unstructured transfinite elements. In section 8 we discuss mesh refinement using a single element, whereas in section 9 we deal with some simple T-spline elements. In general, except for finding the closed-form analytical expressions, we also interpolate within numerous transfinite macro-elements under given boundary conditions of Dirichlet type and determine the L2-error norm.

2. Linear Blending Functions in Coons’ interpolation

2.1 Linear interpolation in the physical space

Given, U₁ and U₂, which denote the values of a function U(x) the end-points of an one-dimensional interval [0,a], at any position x∈[0,a] the value of the function may be linearly interpolated by:

U(x) = (1-x/a)U₁ + (x/a)U₂.

(1)

Figure 1. The quadrilateral ABCD.

Figure 1. The quadrilateral ABCD.

Equation (1) can be extended to a rectangle [0,a] × [0,b], considering the tensor product, i.e., through the expression:

U(x) = (1-x/a)(1-y/b)U₁ + (x/a)(1-y/b)U₂ + (x/a)(y/b)U₃ + (1-x/a)(y/b)U₄.

(2)

2.2. Linear interpolation in the parametric space

If the position x∈[0,a] is normalized to its length ’a’ by introducing the parameter:

ξ = x/a,

(3)

so that ξ∈[0,1], Equation (1) is re-written as follows:

U(x) = (1-ξ)U₁ + ξU₂.

(4)

Moreover, regarding the above rectangle [0,a] × [0,b], the position (x,y) within it can be normalized, so as the x-coordinate is replaced according to Equation (3), whereas the y-coordinate is also normalized as:

η = y/b.

(5)

Then, Equation (2) is re-written as:

U(x) = (1-ξ)(1-η)U₁ + ξ(1-η)U₂ + ξηU₃ + (1-ξ)ηU₄.

(6)

2.3. Projectors within a quadrilateral

Definition of parameters:

The position of a point P(x,y) within a planar quadrilateral ABCD is usually defined by its two Cartesian coordinates, x and y. Alternatively, the same point P(ξ,η) can be uniquely determined using two parameters, ξ and η.

When the quadrilateral is rectangular and the axis origin is at corner A, the definition of the parameters (ξ,η) is direct, according to Equation (3) and Equation (5).

In general, given the parameter ξ, we define the unique points E and F on the edges AD and BC (see, Figure 1), respectively, so that:

AE/AD = BF/BC = ξ.

(7)

The points E and F, which divide the opposite edges AD and BC in the same ratio ξ according to Equation (7), are called ‘corresponding points’ in the ξ-direction.

Similarly, for the parameter η we define the unique points G and H on the edges AB and DC, respectively, so that:

AG/AB = DH/DC = η.

(8)

The points G and H, which divide the opposite edges AB and DC in the same ratio η according to Equation (8), are called ‘corresponding points’ in the η -direction.

Having determined the two couples of points, (E,F) and (G,H), the intersection of the lines EF and GH determines the unique point P(ξ,η). Inversely, for any point P internal to the quadrilateral ABCD, there is a unique line EF passing through it which intersects the opposite edges AD and BC in equal ratios, which is the parameter η. The same holds true for the intersection of the opposite edges AB and DC at the same ratio, which is the parameter ξ.

To understand the meaning of the projector P_ξ and P_η, we introduce the following Theorem-1.

Theorem-1: The point P(ξ,η) within a quadrilateral ABCD divides the segment of its corresponding points (E,F) and (G,H) into the same ratios as those on the opposite edges:

EP/EF = ξ and GP/GH = η.

(9)

Proof: Let us consider point P(x,y) within the quadrilateral ABCD. Obviously, there is a unique line passing through this point which intersects the opposite edges AD and BC at the points E and F, respectively, having the same ratio, i.e., η = AE/AD = BF/BC, according to Equation (7).

Then, the position vector of the point I_ξ which divides the segment EF into the ratio EI_ξ/EF = ξ, will be given as (see, Figure 1)

$$\begin{array}{ll}{\overrightarrow{r}}_{I_{\xi }}& =(1-\xi ){\overrightarrow{r}}_E+\xi {\overrightarrow{r}}_F\\ & =(1-\xi )\left[(1-\eta ){\overrightarrow{r}}_A+\eta {\overrightarrow{r}}_D\right]+\xi \left[(1-\eta ){\overrightarrow{r}}_B+\eta {\overrightarrow{r}}_C\right]\\ & =(1-\xi )(1-\eta ){\overrightarrow{r}}_A+\xi (1-\eta ){\overrightarrow{r}}_B+\xi \eta {\overrightarrow{r}}_C+(1-\xi )\eta {\overrightarrow{r}}_D.\end{array}$$

(10)

Similarly, the position vector of the point I_η which divides the segment GH into the ratio E I_η /EF = η, will be given as (see, Figure 1)

$$\begin{array}{ll}{\overrightarrow{r}}_{I_{\eta }}& =(1-\eta ){\overrightarrow{r}}_G+\eta {\overrightarrow{r}}_H\\ & =(1-\eta )\left[(1-\xi ){\overrightarrow{r}}_A+\xi {\overrightarrow{r}}_B\right]+\eta \left[(1-\xi ){\overrightarrow{r}}_D+\xi {\overrightarrow{r}}_C\right]\\ & =(1-\xi )(1-\eta ){\overrightarrow{r}}_A+\xi (1-\eta ){\overrightarrow{r}}_B+\xi \eta {\overrightarrow{r}}_C+(1-\xi )\eta {\overrightarrow{r}}_D.\end{array}$$

(11)

Since the right-hand sides of Equation (10) and (11) are identical, the vectors ${\overrightarrow{r}}_{I_{\xi }}$and ${\overrightarrow{r}}_{I_{\eta }}$ must be identical, i.e., ${\overrightarrow{r}}_{I_{\xi }}\equiv {\overrightarrow{r}}_{I_{\eta }}$, which means that the points I_ξ and I_η coincide (I_ξ = I_η). Therefore, since (I_ξ ∈ EF ∧ I_η ∈ GH ∧ I_ξ = I_η), we derive that either of the points (I_ξ, I_η) is the unique intersection of the lines EF and GH, which completes the proof. ■

2.4. Interpolation by a Boolean sum

The above ‘Theorem-1’ has probably inspired the so-called “Coons’ interpolation formula”, which refers to an arbitrary curvilinear patch ABCD. In general, the parameterization of this patch depends on the choice of the four corner points (A,B,C,D) on the true patch in the physical space (x,y). Given these corners, the boundary part AB defines the body-fitted (curvilinear) ξ-axis whereas the part AD defines the η-axis. Given the parameters (ξ,η), we define the points (E,F,G,H) on the edges (AD,BC,AB,DC), respectively, and then determine the abovementioned points I_ξ and I_η. In general, these two points are different (I_ξ ≠ I_η), and thus it becomes necessary to find a way to harmonize them. Note that the obvious remedy to take the mean average (I_ξ + I_η)/2 does not work well, because it cannot accurately represent the prescribed boundary conditions. The answer to this problem was given by Steven Coons in 1964, and then in his final report at MIT [3]. For the smooth transition to the problem and for the sake of completeness, we express the essence of Coons’ interpolation in the form of a modified theorem (Theorem-2).

Theorem-2: Given a parametric curvilinear patch ABCD with prescribed boundaries [U(ξ,0), U(1,η), U(ξ,1), U(0,η)], each couple of parameters (ξ,η) determines the position of a point P(ξ,η) which fulfils the boundary conditions, when we apply the following Boolean sum:

$$U(\xi ,\eta )=P_{\xi }\left\{U\right\}+P_{\eta }\left\{U\right\}-P_{\xi \eta }\left\{U\right\}.$$

(12)

In terms of the definitions cited in section 2.3, the three projectors involved in Equation (12) are defined as follows [3]:

$$\begin{array}{l}{P}_{\xi }\left\{U\right\}=(1-\xi )U_E+\xi U_F,\\ P_{\eta }\left\{U\right\}=(1-\eta )U_G+\eta U_H,\\ P_{\xi \eta }\left\{U\right\}=(1-\xi )(1-\eta )U_A+\xi (1-\eta )U_B+\xi \eta U_C+(1-\xi )\eta U_D.\end{array}$$

(13)

The blending functions involved in Equation (13) are defined as:

$$\begin{array}{l}{E}_0(\xi )=(1-\xi ),\hspace{0.33em}{E}_1(\xi )=\xi ,\\ E_0(\eta )=(1-\eta ),\hspace{0.33em}{E}_1(\eta )=\eta .\end{array}$$

(14)

Obviously, the blending function $E_0(\xi )$ is associated to data along the edge AD at ξ=0, whereas $E_1(\xi )$ is associated to data along the edge BC at ξ=1. Similarly, the blending function $E_0(\eta )$ is associated with the edge AB at η=0, whereas $E_1(\eta )$ is associated with the edge DC at η=1. In other words, the subscript of each blending function is the parameter of the edge to which it is referred.

Although the above formula was labelled as ‘Theorem-2’, in fact it is only a heuristic conjecture proposed by Steven Coons, which (again) identically fulfils the boundary conditions. The latter issue can be proved as follows.

Proof: Let us assume that point P tends to the boundary and particularly the edge AB. Then, points E and F tend to the corners A and B, respectively, point G coincides with P, and thus (omitting the notation {U}), since η = 0, Equation (13) shows that the three projectors will be:

$$\begin{array}{}{\left.P_{\xi }\right|}_{AB}=(1-\xi )U_A+\xi U_B,\\ {\left.P_{\eta }\right|}_{AB}=(1-0)U_P+0U_H=U_P,\\ {\left.P_{\xi \eta }\right|}_{AB}=(1-\xi )(1-0)U_A+\xi (1-0)U_B+\xi 0U_C+(1-\xi )0U_D=(1-\xi )U_A+\xi U_B\end{array}$$

(15)

Substituting Equation (15) into Equation (12), we receive:

$$U(\xi ,0)=U_P,\hspace{1em}\mathrm{Q}.\mathrm{E}.\mathrm{D}.$$

(16)

Therefore, Coons’ interpolation formula [Equation (12)], fulfils the prescribed boundary condition along the edge AB, and similarly all the other edges. ■

3. Higher-order Blending Functions in Gordon’s (transfinite) interpolation

3.1. From Coons’ patch to Gordon’s patch theory

The disadvantage of Coons’ interpolation formula is that it cannot deal with internal data at all. For example, if the boundary values are homogeneous (i.e., U = 0), all the three projectors in Equation (13) vanish since they interpolate zero values, and therefore Equation (12) leads to zero values of the interpolated function U(ξ,η), everywhere within the patch. This shortcoming was overcome by William Gordon, who extended the meaning of the three projectors involved in Equation (12) [4,5].

In more detail, the shape of the patch ABCD, i.e. of the function U(ξ,η), is controlled through (m+1) ‘stations’ perpendicular to the ξ-axis and (n+1) stations perpendicular to the η-axis, where the three projectors involved are given as

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

(17)

The functions $E_i(\xi ),\hspace{0.17em}\hspace{0.17em}i=1,\dots ,m+1$ are again called “blending functions” toward the ξ-direction and generally are of higher degree than the liner ones used in Coons’ interpolation (see Equation (14)). They interpolate the functions $U({\xi }_i,\eta )$ at $m+1$ stations, and usually fulfil the Kronecker’s delta condition $E_i({\xi }_j)={\delta }_{ij}$.

According to Gordon’s theory, considering the data along all the stations and then calculating the three projectors according to Equation (17), the interpolation is given by the Boolean sum:

$$U(\xi ,\eta )=P_{\xi }\otimes P_{\eta }=P_{\xi }\left\{U\right\}+P_{\eta }\left\{U\right\}-P_{\xi \eta }\left\{U\right\}.$$

(12)

repeated

To express the bivariate function $U(\xi ,\eta )$ in terms of discrete values, $U_i$, it is necessary to interpolate the $(m+1)$ univariate functions $U({\xi }_i,\eta )$ at each station ${(S_i)}_{\xi ={\xi }_i=const.}$ using the so-called “trial functions”. Similarly, each of the $(n+1)$ univariate functions $U(\xi ,{\eta }_j)$ at the corresponding station ${(T_j)}_{\eta ={\eta }_j=const.}$ is interpolated using a separate set of trial functions.

In summary, the transfinite interpolation requires a set of blending functions associated with stations perpendicular to the parameter axes ξ and η, as well as sets of trial functions which interpolate the univariate data along each station.

3.2. Transfinite interpolation

The term “transfinite interpolation” is strictly related to Gordon’s interpolation formula, which -in addition to the boundary- includes internal points as well. Transfinite interpolation usually refers to structured patches in which the (m+1) stations are similar one another, and the same holds for the other set of (n+1) stations. In this configuration, the nodes are distributed along the stations. This, in turn, means that between two parallel neighboring stations, there may be intermediate nodes, as shown in Figure 2. One may observe therein that each of the 21-node, 27-node and 33-node elements have three stations in each direction, at the locations ξ = 0,1/2,1 and η = 0,1/2,1. Moreover, the 113-node element sown in Figure 2(d) includes 5 vertical (at ξ = 0, 1/4, 2/4, 3/4, 1) and 4 horizontal (at η = 0, 1/3, 2/3, 1) stations.

The advantage of such a setup is that we save computer resources since we avoid the ideal tensor product, which is full of nodes.

Another important issue in the implementation of the transfinite interpolation is the choice of the univariate interpolation along the horizontal and vertical stations. This is the reason that a special part of the present paper (section 4) is devoted to it.

3.3. Tensor-product under the umbrella of transfinite interpolation

In general, the data are given at discrete positions such as the vertical stations (S1 to S5) and the horizontal ones (T1 to T5) shown in Figure 3. This makes it necessary to interpolate the univariate functions $U({\xi }_i,\eta )$ and $U(\xi ,{\eta }_j)$, which appear in Equation (17), along these stations.

Before we go on with Lagrange and Bernstein polynomials, for the sake of completeness we provide some historical facts about interpolation.

3.3.1. Historical note on univariate interpolation

Lagrange interpolation was published in the work of Joseph Louis Lagrange “Théorie des fonctions analytiques” (Theory of Analytical Functions) in 1797, which is available in a later edition [18]. This topic was later discussed by a few authors such as Jenkins and Greville, as reported in 1946 by Schoenberg [19], while a more recent review paper is [20]. Today, Lagrange interpolation is cited in every textbook of interpolation as [21] and numerical analysis (e.g., [22]).

Polynomials in Bernstein form were first used by Sergei Bernstein, in a constructive proof for the Weierstrass approximation theorem [23]. The use of Bernstein polynomials in the treatment of differential equations was continued by others (see Davis [21] and papers therein). While other forms of polynomials such as Chebyshev and Legendre attracted the researchers in computational mathematics due to simplicity in their differentiation compared to Lagrange polynomials (e.g., Lanczos [24]), Bernstein polynomials were rewarmed under the seminal contribution of Pierre Bézier, a French engineer at the automobile industry Renault, to describe curves and surfaces [25]. As a result, today the term Bernstein-Bézier polynomials is widely used.

3.3.2. Tensor product using Lagrange polynomials

The tensor product of Lagrange polynomials emerged naturally as a result of combining two key ideas: Lagrange interpolation and the tensor product framework (which is widely used in mathematics, particularly in multivariate settings). This method became popular in the development of numerical analysis for interpolating multivariate functions by extending the Lagrange polynomial from one dimension to higher dimensions. It gained prominence in the 20th century with the rise of finite element methods and numerical techniques in computational mathematics.

3.3.3. Generalization of tensor product

In the continuous process of creating controllable curves and surfaces, the transfinite interpolation goes beyond the tensor product of Lagrange polynomials, alternatively using Bernstein polynomials, B-splines and NURBS for trial functions (see section 4). Moreover, transfinite interpolation includes the tensor product of any other trial functions (not only Lagrange polynomials) as a special case, as will be described by the following Theorem-3.

Theorem-3: In tensor product, like Fig. 3(a), all the three projectors are equal one another, i.e.:

$$P_{\xi }=P_{\eta }=P_{\xi \eta },$$

(18a)

and therefore the Boolean sum becomes

$$P_{\xi }\otimes P_{\eta }=P_{\xi }+P_{\eta }-P_{\xi \eta }\equiv P_{\xi \eta }.$$

(18b)

Proof: Given the (m+1) stations perpendicular to ξ-axis and the (n+1) ones perpendicular to the η-axis, and assuming the same set of interpolation functions (also called “trial” functions), ${\tilde{B}}_i(\xi ),i=1,\dots ,m+1$ and ${\tilde{B}}_j(\eta ),j=1,\dots ,n+1$, we have:

$$U(\xi ,\overline{\eta })={\displaystyle \sum _{i=1}^{m+1}{\tilde{B}}_i(\xi )U_i(\xi ,\overline{\eta })},\hspace{1em}U(\overline{\xi },\eta )={\displaystyle \sum _{j=1}^{n+1}{\tilde{B}}_j(\eta )U_j(\overline{\xi },\eta )},$$

(19)

where the bar denotes constant values.

The characteristic of a tensor product is that the trial functions are the same with the blending functions toward the same direction, i.e.

$${\tilde{B}}_i(\xi )\equiv E_i(\xi ),i=1,\dots ,m+1 \mathrm{and} {\tilde{B}}_j(\eta )\equiv E_j(\eta ),j=1,\dots ,n+1$$

(20)

Substituting Equation (20) into Equation (17), Equation (18a) is derived, and this completes the proof. ■

4. Trial functions

Trial functions are the univariate functions ${\tilde{B}}_i(\xi )$ and${\tilde{B}}_j(\eta )$, which are involved in Equation (19), aiming at interpolating or approximating the value of the function $U(\xi ,\eta )$ along the stations of the patch. In principle, they can be any known functional set, such as Lagrange and Bernstein polynomials, B-splines, NURBS, Fourier series, etc. In this paper we restrict ourselves to the first four types only, which are discussed in detail below.

4.1. Lagrange polynomials

It is well known that when a univariate function $u(x)$ is represented by a polynomial of degree $m$,

$$u(x)=a_1+a_{2}x+\cdots +a_{m+1}{x}^m,$$

(21)

it can be also expressed in terms of Lagrange polynomials, $L_i(x)$:

$$u(x)={\displaystyle \sum _{i=1}^{m+1}{L}_{i,m}(x)}{U}_i,$$

(22)

where

$$\begin{array}{ll}{L}_{i,m}(x)& ={\displaystyle \frac{(x-x_1)(x-x_2)\dots (x-x_i)(x-x_{i+2})\dots (x-x_{m+1})}{(x_i-x_1)(x_i-x_2)\dots (x_i-x_i)(x_i-x_{i+2})\dots (x_i-x_{m+1})}}\\ & ={\displaystyle \prod _{\begin{array}{l}j=1\\ j\ne i\end{array}}^{m+1}\frac{(x-x_j)}{(x_i-x_j)}}\end{array}$$

(23)

is the Lagrange polynomial of degree m, and $U_i,i=1,\dots ,m+1$ is any set of nodal values of the polynomial at the positions $x_1,\dots ,x_{m+1}$.

Changing the positions of the nodal points, a new set of Lagrange polynomials is derived. The most usual case is the uniform location of the $x_i$’s.

4.2. Bernstein polynomials

On the other hand, the same function u(x) can be also expressed in terms of Bernstein-Bézier basis polynomials, $B_{i,m}(x)$:

$$u(x)={\displaystyle \sum _{i=1}^{m+1}{B}_{i,m}(x)}{\alpha }_i,$$

(24)

where

$$\begin{array}{ll}{B}_{i,m}(x)& =\left(\begin{array}{c}m\\ i-1\end{array}\right){(1-x)}^{m-i+1}{x}^{i-1}\\ & ={\displaystyle \frac{m!}{(i-1)!(m-i+1)!}}{(1-x)}^{m-i+1}{x}^{i-1},\hspace{0.33em}i=1,\dots ,m+1\end{array}$$

(25)

and ${\alpha }_i$ are generalized coordinates, i.e., coefficients without physical meaning.

4.3. Relationship between Lagrange and Bernstein polynomials

Since in Eq. (22) and Eq. (24), the function u(x) shares the same monomials $(x^i,i=0,\dots ,x^m)$, the identical equation of these two polynomials demands the equality of all the coefficients in the monomials (powers) of the same degree $(x^i,i=0,\dots ,x^m)$. This is turn gives a system of $(m+1)$ equations, which relates the $(m+1)$ Lagrange polynomials with the $(m+1)$ Bernstein basis polynomials.

Another way is to consider $(m+1)$ test points in the domain $\left[x_1,x_{m+1}\right]$ and apply to them the identity between Eq. (22) and Eq. (24):

$$\underset{(m+1)\times 1}{\underbrace{U_{test}}}=\underset{(m+1)\times (m+1)}{\underbrace{\left[L_{test}\right]}}\underset{(m+1)\times 1}{\underbrace{U_{nodal}}}=\underset{(m+1)\times (m+1)}{\underbrace{\left[B_{test}\right]}}\underset{(m+1)\times 1}{\underbrace{A}}$$

(26)

If it happens to consider the test points exactly at the nodal points, $x_1,\dots ,x_{m+1}$, the matrix $\left[L_{test}\right]$ equals to the unity matrix ($\left[L_{test}\right]=I$) while the other becomes $\left[B_{test}\right]=\left[B_{nodal}\right]$, and hence

$$A={\left[B_{nodal}\right]}^{-1}{U}_{nodal}$$

(27)

Substituting Eq. (27) in the last equality of Eq. (28) below, which equates the two polynomial expressions,

$$\begin{array}{ll}u(x)& =\left[\begin{array}{ccc}{L}_1(x)& \dots & L_{m+1}(x)\end{array}\right]\cdot \left\{U_{nodal}\right\}\\ & =\left[\begin{array}{ccc}{B}_1(x)& \dots & B_{m+1}(x)\end{array}\right]\cdot \left\{A\right\}\end{array}$$

(28)

and then deleting the common vector, $\left\{U_{nodal}\right\}$, we eventually derive:

$$L=T\cdot B$$

(29)

where

$$T={\left[B_{nodal}^{-1}\right]}^T,$$

(30a)

is the transformation matrix,

$$L={\left[\begin{array}{ccc}{L}_1(x)& \dots & L_{m+1}(x)\end{array}\right]}^T,$$

(30b)

also

$$B={\left[\begin{array}{ccc}{B}_1(x)& \dots & B_{m+1}(x)\end{array}\right]}^T,$$

(30c)

and

$$B_{nodal}=\left[\begin{array}{cccc}{B}_1(x_1)& B_2(x_1)& \cdots & B_{m+1}(x_1)\\ B_1(x_2)& B_2(x_2)& \cdots & B_{m+1}(x_2)\\ ⋮& ⋮& \ddots & ⋮\\ B_1(x_{m+1})& B_2(x_{m+1})& \cdots & B_{m+1}(x_{m+1})\end{array}\right].$$

(30d)

4.3.1. Uniform distribution of nodal points

Using either of the two abovementioned methodologies to determine T, for the uniform polynomials of low degree, we can derive the following transformation matrices, which fulfil Equation (29):

$$\begin{array}{}\mathrm{Linear}, p=1:\hspace{0.33em}{T}_1=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],\\ \mathrm{Quadratic}: p=2:\hspace{0.33em}{T}_2=\left[\begin{array}{ccc}1& -1/2& 0\\ 0& 2& 0\\ 0& -1/2& 1\end{array}\right],\\ \mathrm{Cubic}: p=3:\hspace{0.33em}{T}_3=\left[\begin{array}{cccc}1& -5/6& 1/3& 0\\ 0& 3& -3/2& 0\\ 0& -3/2& 3& 0\\ 0& 1/3& -5/6& 1\end{array}\right],\\ \mathrm{Quartic}: p=4:\hspace{0.33em}{T}_4=\left[\begin{array}{ccccc}1& -13/12& 13/18& -1/4& 0\\ 0& 4& -32/9& 4/3& 0\\ 0& -3& 20/3& -3& 0\\ 0& 4/3& -32/9& 4& 0\\ 0& -1/4& 13/18& -13/12& 1\end{array}\right],\end{array}$$

(30e)

and so on.

4.3.2 Non-uniform distribution of nodal points

For each set of internal nodal points [x₂, …, x_m] within the domain [x₁, x_m+1], a different set of Lagrange polynomials, {L₁(x), …, L_m+1(x)}, is produced. Always, the basis of Lagrange polynomials is associated to the nodal values U_i at the nodes [x₁, x_m+1]. Sometimes the position of the points is critical for the accurate representation of non-polynomial functions u(x), e.g. those including high gradients, singularities etc. For each set of the nodal points, {x₁, x_m+1}, Equation (30a) determines a corresponding transformation matrix T_p.

On the other hand, in the representation of the same functions u(x) using Bernstein polynomials, we do not use internal points, nor do we change the form of the involved basis functions. To achieve the same accuracy as node shifting in Lagrange interpolation, the only matter which is influenced in Bernstein interpolation is the value of the generalized coefficients ${\alpha }_i$ involved in Equation (24). Looking the same issue from the scope of computer-aided geometric design (CAGD), the function u(x) is replaced by the form of a curve C(x), and the values of the coefficients ${\alpha }_i$ are substituted by the so-called control points $P_i$.

In this context, the uniform position of nodal points used in Lagrange representation is equivalent to the uniform location of the control points, which also lie in the positions where the Bernstein polynomials take their maximum values. Moreover, shifting the nodal points at non-uniform positions and at the same time preserving the Bernstein polynomials, if we demand that the approximation curves of degree $p$ (Lagrange and Bernstein) pass through these (non-uniform) points, it is found that the representations are identical (see, Appendix A).

4.4 B-splines

In one-dimensional problems, the variable $u(x)$ can be approximated by B-splines, merely by substituting the Bernstein-Bézier polynomials in Equation (24) by the piecewise polynomial B-spline (basis) functions $N_{i,p}$ of degree $p$ (order $k=p+1$), as follows:

$$u(x)={\displaystyle \sum _{i=1}^{m+1}{N}_{i,p}(x)}{\alpha }_i,$$

(39)

where $m+1$ is the number of control points in the $x$-direction.

The functions $N_{i,p}$ are very versatile, dependent on the knot vector $\Xi$ and the polynomial degree $p$.

For example, considering the nodes $\left\{x_1,x_2,\dots ,x_{m+1}\right\}$ and $p=1$, if the knot vector becomes $\Xi =\left[\underset{double}{\underbrace{x_1,x_1}},x_2,x_3,\dots ,x_{m-1},x_m,\underset{double}{\underbrace{x_{m+1},x_{m+1}}}\right]$, we eventually derive $(m+1)$ shape functions, exactly as those piecewise-linear in the conventional linear finite elements.

Another extreme situation is when the knot vector includes only its ends with multiplicity $\lambda =p+1$, and then the set of B-splines coincides with the set of non-rational Bernstein-polynomials, which −as previously discussed− is a linear transformation (and thus equivalent) to the set of Lagrange polynomials.

Remark: Obviously, B-splines includes the class of Bernstein polynomials, as a subclass. Moreover, since Bernstein polynomials constitute a linear transformation of Lagrange polynomials, this in turn means that the space of B-splines includes the subspace of Lagrange polynomials.

4.5. Non-uniform Rational B-splines (NURBS)

In one-dimensional problems concerning the representation of a curve $C(\xi )$, it is possible to use weights aiming at strengthening or weakening the involved B-splines, and thus to accurately represent it.

In general, the 1D analogue of Eq. (39) is

$$u(x)={\displaystyle \sum _{i=1}^{m+1}{R}_{i,p}(x)}{\alpha }_i,$$

(41)

with

$$R_{i,p}(x)={\displaystyle \frac{w_i{B}_i(x)}{{\displaystyle \sum _{j=1}^{m+1}{w}_j{B}_j(x)}}},$$

(42)

where $m+1$ is the number of control points.

5. Tensor-products of trial functions

Having introduced all the important univariate (trial) functions in section 4, in this section we proceed with the tensor product of them.

Figure 4. The quadrilateral ABCD. Arrangement of (n +1) × (m +1) nodal points in a 2D patch.

Figure 4. The quadrilateral ABCD. Arrangement of (n +1) × (m +1) nodal points in a 2D patch.

5.1. Lagrange polynomials

Let a 2D patch, on the plane $Oxy$, interpolated through a net of $(m+1)$ vertical and $(n+1)$ horizontal lines (also called stations). The vertical lines intersect the $x$-axis at the nodal points $\left\{x_1,x_2,\dots ,x_{m+1}\right\}$, whereas the horizontal lines intersect the $y$-axis at the nodal points $\left\{y_1,y_2,\dots ,y_{n+1}\right\}$, as shown in Figure 4. In general, the $i$-th vertical line (which passes through the point $x=x_i$) will intersect the $j$-th horizontal line (which passes through the point $y=y_j$) at the point $P_{i,j}(x_{i,j},y_{i,j})=(x_i,y_j)$, which is associated to a value $U_{i,j}$. In other words, the first index “i” is associated to the $x$-direction, while the second index “j” is associated to the $y$-direction.

When the variable is bivariate, $u(x,y)$, Equation (22) is replaced by the tensor product:

$$u(x,y)={\displaystyle \sum _{j=1}^{n+1}{\displaystyle \sum _{i=1}^{m+1}{L}_i(x)}{L}_j(y)U_{i,j}},$$

(31)

where $U_{i,j}$ are the nodal values associated to the Lagrange basis functions $L_i(x)$ and $L_j(y)$. In our formulation, the totality of nodal points is stored into a matrix $U$ of size $(m+1)\times (n+1)$, arranged in $(m+1)$ matrix rows and $(n+1)$ matrix columns, as follows:

$$U=\left[\begin{array}{ccc}{U}_{1,1}& \cdots & U_{1,(n+1)}\\ ⋮& \ddots & ⋮\\ U_{(m+1),1}& \cdots & U_{(m+1),(n+1)}\end{array}\right]$$

(32)

In other words, the $j$-th column of matrix U in Eq. (12) represents the $j$-th row in the spatial arrangement of Figure 4, measured from bottom. Similarly, the $i$-th row of matrix U corresponds to the $i$-th column of Figure 4, measured from left.

Given the abovementioned net of nodal points stored in $U$, the tensor product in this patch is written in the form:

$$u(x,y)=L_x^T\cdot U\cdot L_y,$$

(33)

with

$$L_x={\left[\begin{array}{ccc}{L}_1(x)& \cdots & L_{m+1}(x)\end{array}\right]}^T$$

(34a)

and

$$L_y={\left[\begin{array}{ccc}{L}_1(y)& \cdots & L_{n+1}(y)\end{array}\right]}^T$$

(34b)

5.2 Bernstein polynomials

Based on the matrix form of tensor-product Lagrange basis polynomials (Equation (33)), it is easy to show that a similar tensor-product Bernstein basis polynomials is derived as well.

By substituting Equation (29) into Equation (33), we receive:

$$\begin{array}{ll}u(x,y)& ={\left(T_x{B}_x\right)}^T\cdot U\cdot \left(T_y{B}_y\right)\\ & \hspace{0.17em}=B_x^T\cdot \left(T_x^{T}U{T}_y\right)\cdot B_y\end{array}$$

(35)

By introducing the matrix:

$$A=T_x^{T}U{T}_y,$$

(36)

Equation (15) is written as:

$$u(x,y)=B_x^T\cdot A\cdot B_y$$

(37)

Since Equation (37) has the same form with Equation (33) which is a tensor product of Lagrange polynomials, it is obvious that the bivariate function can be written as a tensor product of Bernstein polynomials

$$u(x,y)={\displaystyle \sum _{j=1}^{n+1}{\displaystyle \sum _{i=1}^{m+1}{B}_{i,m}(x)}{B}_{j,n}(y){\alpha }_{ij}},$$

(38)

associated to arbitrary coefficients ${\alpha }_{i,j}$ which form a matrix A given by Equation (16).

In other words, the non-rational tensor-product of Bernstein-Bézier basis polynomials is only a linear transformation of the tensor-product of Lagrange basis polynomials.

5.3 . B-splines

Since the univariate B-spline is a generalization of the Lagrange and Bernstein polynomials, simply by inserting a piecewise character in the interpolation, it is reasonable that bivariate interpolation may be carried by the analogy of Equation (38), that is by the formula:

$$u(x,y)={\displaystyle \sum _{j=1}^{n+1}{\displaystyle \sum _{i=1}^{m+1}{N}_{i,p}(x)}{N}_{j,q}(y){\alpha }_{ij}},$$

(40)

where p and q are the polynomial degrees in the x- and y-direction, respectively.

5.4. NURBS approximation

By the analogy of Equation (40), the interpolation in a curvilinear NURBS patch is given by the formula:

$$u(x,y)={\displaystyle \sum _{j=1}^{n+1}{\displaystyle \sum _{i=1}^{m+1}{R}_{ij,pq}(x,y)}{\alpha }_{ij}},$$

(43)

where $R_{ij,pq}(x,y)$ is the bivariate NURBS basis function:

$$R_{ij,pq}(x,y)={\displaystyle \frac{w_{ij}{N}_{i,p}(x)N_{j,q}(y)}{{\displaystyle \sum _{j=1}^{n+1}{\displaystyle \sum _{i=1}^{m+1}{w}_{ij}{N}_{i,p}(x)N_{j,q}(y)}}}},$$

(44)

The latter is a true tensor-product only when the weight $w_{ij}=w_i{w}_j$, otherwise the denominator of Equation (44) includes the sum of all the $(m+1)\times (n+1)$ terms.

6. The proposed Conjecture

Considering the material which was systematically presented and discussed in the previous sections, this paper makes the innovative step beyond the state-of-the-art and tests the hypothesis (or conjecture) that:

CONJECTURE: In the transfinite interpolation, the blending functions may be any complete and linear independent set provided the Partition of Unity Property is preserved.

It is worthy to mention that the original transfinite interpolation formula (by Coons and later by Gordon) is heuristic, and thus cannot be mathematically sustained (in absolute sense). According to Ref. [7](p. 114), ”. . . it is guided by experience, geometric intuition and analysis, and is accomplished with the aid of computer graphics, i.e. actual visual inspection”.

Since the original source is heuristic, the novel idea is obviously heuristic as well.

The proposed conjecture was progressively established in author’s mind, step-by-step, as follows. In the beginning, it was found that it is sufficient to ‘mechanically’ replace all the Lagrange polynomials involved in the transfinite interpolation of a curvilinear patch (i.e. the blending and the trial functions) by Bernstein-Bézier polynomials of the same degree. In this way, it was made possible to determine the weights associated with the control points so that the patch may accurately represent a circular or a spherical part [17]. The amazing issue related to the previously mentioned mechanical replacement, is that the produced shape functions are fully symmetrical with respect to the geometry of the patch (see, [17]), and this was precisely the cause of the inspiration for the present conjecture.

The argumentation in favour of the above conjecture continues below.

6.1. Bernstein tensor-product elements

In this context, first let us consider the tensor product Bernstein elements (e.g. Fig. 3a). Obviously, here the assumption that blending and trial functions are both Bernstein polynomials denoted by the symbol $B$ $(E_i=B_i,{\tilde{B}}_j=B_j)$, leads to three equal projectors in the transfinite interpolation, i.e. $P_{\xi }=P_{\eta }=P_{\xi \eta }$, and therefore the Boolean sum becomes $P_{\xi }+P_{\eta }-P_{\xi \eta }\equiv P_{\xi \eta }$ (compare with Theorem-1). Since the assumption of Bernstein polynomial as blending functions lead to the anticipated tensor-product, the proposed Conjecture is verified. As an example, the case of uniform control points in the tensor-product Bernstein-Bézier patch is shown in Fig. 3a. We recall that, according to section 2.2, the uniform tensor-product of Bernstein polynomials of fourth degree (shown in Fig. 3a) is equivalent to the unform tensor-product of Lagrange polynomials of the same degree (shown in Fig. 2a).

6.2. B-spline tensor-product elements

Let us now consider a B-spline patch of degree p. Again, the assumption of both blending and trail functions as univariate B-splines, leads to a Boolean sum equal to the anticipated tensor-product of B-splines. Therefore, again, our hypothesis is still reliable.

6.3. NURBS tensor-product elements

Considering now a NURBS patch with weights $w_{ij}=w_i{w}_j$, and repeating the same procedure, the Boolean sum leads to a tensor-product NURBS patch.

6.4. Overall evidence for tensor-product elements

Therefore, summarizing the above findings and the review of this section, we report:

• Regarding the tensor-product, for any of the four kinds of trial functions (i.e. Lagrange, Bernstein, B-splines, NURBS), this is the outcome of transfinite interpolation in which the blending functions are identical with the trial functions.
• Piecewise-linear blending functions are admissible, since when combined with piecewise-linear trial functions lead to an assembly of conventional bilinear finite elements (see, Ref. [26](p.956)).
• Pure transfinite elements, with nodal points distributed along stations (i.e. different than tensor-product arrangement), gave the same result when the Lagrange polynomials were mechanically replaced by Bernstein polynomials. This is equivalent to considering Bernstein polynomial for both blending and interpolating purposes.
• All the above successful cases are characterized by completeness, linear independence and all of them have the Partition of Unity Property.

6.5. Applications

6.5.1. The 21-node transfinite element

Regarding the 21-node transfinite element (Figure 2a), since there are three stations per direction (ends and middle), the blending functions may be Lagrange polynomials of degree 2, whereas the trial functions may be Lagrange polynomials of degree 4 (because there are five nodes along each edge and any internal station), as shown in Figure 5a.

Implementing the conjecture of the present paper, if we assume that the blending functions are Bernstein polynomials of degree 2 whereas the trial functions are Bernstein polynomials of degree 4, the final closed-form expression of the global shape functions is merely the replacement of the Lagrange by the Bernstein polynomials, as shown in Figure 5b. The latter has been previously performed by intuition in Ref. [17], but the present paper generalizes this idea and promotes the implementation of other trial functions apart from Bernstein polynomials as well (e.g., B-splines, NURBS).

In addition to that, we can preserve the blending function as a Bernstein polynomial of degree 2, whereas the trial function must be consistent to the five nodal points along each station. The latter can be accomplished by five control points which are consistent with either of the following choices for the trial functions:

• A Bernstein polynomial of degree p = 4 and the obvious knot vector U = [0,0,0,0,0, 1,1,1,1,1]. This choice was already mentioned above and is shown in Figure 5b.
• A cubic B-spline (p = 3) with knot vector U = [0,0,0,0, 1/2, 1,1,1,1], shown in Figure 5c.
• A quadratic B-spline (p = 2) with knot vector U = [0,0,0, 1/3, 2/3, 1,1,1], shown in Figure 5d.

To give an idea about the performance of these elements, at the 21 nodes of Figure 2a we collocate the function $f(x,y)=e^x\cos (\pi y)$, and then interpolate in the interior of the element using the 21 global shape functions. The results are shown in Table 1, where one may observe that Lagrange and Bernstein polynomials of the same degree led to the same result. The replacement of the fully polynomial trial functions with B-splines decreased accuracy but not that much.

6.5.2. The 113-node transfinite element

The concern with the 113-node element (Figure 2d), is whether the involved Lagrange and Bernstein polynomials of higher degree (p_x = 16, p_y = 12) lead to numerical oscillations. Actually, one may observe the high values in Figure 6(a) when Lagrange polynomals are used.

This shortcoming is overcome by the proposed novel procedure, according to which we can consider that the blending functions be B-splines of degree (say) $p=3$, in conjunction with the proper knot vector which produces so many control points as the number of nodal points per station.

In contrast, Figure 6(b) shows that when both blending and trial functions are chosen as Bernstein polynomials, the positive high values of the basis functions fully disappear, but small negative values still remain.

One would expect that the non-negative Bernstein polynomials would also lead to non-negative basis functions, which of course is true when only tensor products are involved. To explain this issue, below we isolate the first two out of the 113 basis functions:

$$\begin{array}{l}{\varphi }_1(\xi ,\eta )=E_0(\xi )L_0^{(12)}(\eta )+L_0^{(16)}(\xi )E_0(\eta )-E_0(\xi )E_0(\eta ),\\ {\varphi }_2(\xi ,\eta )=L_{\frac{1}{16}}^{(16)}(\xi )E_0(\eta ).\end{array}$$

(45a)

The first equality in the Lagrange-polynomial based Equation (45a) is associated to the node ‘1’ and corresponds to a typical intersection between two stations. The second equality is associated to node ‘2’ and refers to a typical intermediate node of a station. Based on the number of stations per direction, the blending functions $E_0(\xi )$ and $E_0(\eta )$ are Lagrange polynomials of degree 4 and 3, respectively. Moreover, $L_i^{(16)}(\xi )$ are Lagrange polynomials of degree 16 for interpolating along the horizontal stations, whereas $L_j^{(12)}(\eta )$ are Lagrange polynomials of degree 12 for interpolating along the vertical stations. Since Lagrange polynomials take positive and negative values, the shape of Figure 6(a) is fully justified.

Furthermore, even if the Lagrange polynomials are replaced by the non-negative Bernstein polynomials as follows: $E_0(\xi )\to B_{0,4}(\xi )$, $E_0(\eta )\to B_{0,3}(\eta )$, etc, $L_0^{(16)}(\xi )\to B_{0,16}(\xi )$, etc, the basis functions ${\psi }_1(\xi ,\eta ),{\psi }_2(\xi ,\eta )$ and the similars may obviously obtain negative values:

$$\begin{array}{l}{\psi }_1(\xi ,\eta )=B_{0,4}(\xi )B_{0,12}(\eta )+B_{0,16}(\xi )B_{0,3}(\eta )-B_{0,4}(\xi )B_{0,3}(\eta ),\\ {\psi }_2(\xi ,\eta )=B_{1,16}(\xi )B_{0,3}(\eta ).\end{array}$$

(45b)

To evaluate the quality of the approximation using transfinite interpolation, in both cases (Lagrange-Lagrange and Bernstein-Bernstein), after collocating the function $f(x,y)=e^x\cos (\pi y)$ at the 113 nodal points, Table 2 shows that the L2-error norm equals to only (1.2326 × 10^-5) %. Since the maximum polynomial degree involved in the L2-norm is 2 × 16 = 32, a 17 × 17 Gaussian scheme was imposed for the entire patch.

In the case of the combination of blending functions using either Lgrange or Bernstein polynomials in conjunction with cubic B-splines based on 14 and 10 uniform knot spans in x- and y-directions (thus producing 17 and 13 control points), respectively, the basis functions are shown in Figure 6(c). In the latter case, the domain was divided into 14 × 10 uniform cells (elements), and on each element a 4 × 4 Gaussian scheme was imlemented. The L2-error norm was found equal to 0.0013%.

7. Unstructured transfinite elements

Now we leave the structured stencils shown in Figure 2 (developed at General Motors Research Laboratories), and deal with patches in which the nodes (or the control points) do not lie along isolines only, as was the case discussed in section 6. There are many gradations of this group, at least as follows:

• The internal nodes form a tensor product, but are not consistent with the boundary nodes (neither the position nor the quantity);
• The internal nodes have been produced after a number of subdivisions (see Section 8);
• The internal nodes are arranged like the well-known T-mesh (see Section 9).

All the above three categories may be tackled using univariate blending and trial functions which may be (i) Lagrange polynomials, (ii) Bernstein polynomials, (iii) B-splines, and (iv) NURBS (when curvilinear).

The treatment of all these cases may be faced with a well-programmed computer program, but for the sake of clariness, below we present typical cases separately.

7.1. Tensor-product-like 25-node element

To this purpose, let us consider the 25-node element shown in Figure 3b, which is a distrorted tensor-product of which some boundary nodes have been shifted from their uniform position. Therein, one may observe 25 nodes (16 on the boundary in form of black circles, and 9 inner nodes denoted by red-coloured circles). The characteristic of this element is that some of the normal projections of the inner nodes on the boundary are shifted with respect to the boundary positions. To depict this deviation, these projections (called artificial points) are denoted by smaller auxiliary numbers in green (nodes 26,…,35), and green-coloured crosses (×).

The first question to answer below is, whether the Lagrange-polynomial and the Bernstein-polynomial based numerical solution is the same, and how it is performed and justified.

7.1.1. Lagrange-polynomial based elements

The number of layers of internal nodes per direction defines the degree of the blending functions, which in principle are Lagrange polynomials. Applying the transfinite (Gordon) interpolation to the assemblage of the internal nodes and the artificial points, after all these artificial nodes are interpolated in terms of the actual nodes along the corresponding edge, we eventually derive two types of global shape functions as follows:

• Those along the boundary, depending on both the trial functions and the associated blending function in the normal direction. By construction, they vanish at all the boundary nodes (except for one) and vanish at all the internal nodes.
• Those global shape functions, associated with the inner nodes. They are the well-known classical tensor product Lagrangian polynomials (i.e., the blending functions).

For the distorted setup of nodes shown in Figure 3b, the closed-form expressions of the Lagrange-polynomial based global shape functions have been derived in Appendix B and are illustrated in Figure 7(a).

7.1.2. Bernstein-polynomial based elements

While the Lagrange-based global basis functions depend on whether the associated nodal point is at the corner (and then is an actual Boolean sum of blending and trial functions) or is intermediate on the boundary, Bernstein polynomials operate in a slightly different way as is discussed below.

For the distorted setup of twenty-five nodes shown in Figure 3b, the stations to which the Bernstein polynomials refer are uniformly distributed, and thus blending and trial sets coincide one another. Consequently, in Equation (B.5) (of Appendix B) the Lagrange polynomials must be replaced by Bernstein polynomials of the same degree, and the result is illustrated in Figure 7(b). Obviously, this means that for all the internal nodes and the intermediate boundary nodes, the global basis functions will be tensor-product Bernstein polynomials. Furthermore, if the same replacement is performed for a corner node (say the bottom left No. 1), the result will be again a perfect tensor product:

$$\begin{array}{ll}{\left[N_1(\xi ,\eta )\right]}_{\mathrm{Bernstein}}& =B_{0,4}(\xi )B_{0,4}(\eta )+{\bcancel{B}}_{0,4}(\eta )B_{0,4}(\xi )-{\bcancel{B}}_{0,4}(\xi )B_{0,4}(\eta )\\ & =B_{0,4}(\xi )B_{0,4}(\eta ).\end{array}$$

(46)

In other words, the analogue of a distorted tensor-product Lagrangian transfinite element is an ideal tensor-product Bernstein polynomials (with symmetrically placed control points), and this is possible because any deviation from symmetry is absorbed by the generalized coefficients α.

Interestingly, setting the values of the function $u(x,y)=e^x\cos (\pi y)$ at the same 25 non-uniform places, in both formulations (Lagrange, Bernstein) the L2-error norm within the square domain was found equal to L2 = 0.6361%.

7.2. Tensor-product-like 27-node element

Another instructive example of a generalized 27-node transfinite element, a little more complex than that described above, is shown in Figure 8a. One may observe that the opposite edges are divided into a uniform but different number of subdivisions to one another.

The three stations of internal nodes parallel to the horizontal axis depict blending functions $E_i(\xi )$ of fourth degree in $\xi$ (with $i=1,\dots ,5$). Similarly, the three stations of internal nodes parallel to the vertical axis (thus five stations, in total, including the edges AD and BC) depict blending functions $E_j(\eta )$ again of fourth degree in $\eta$.

Regarding the trial functions, the edge AB requires a polynomial of fourth degree (i.e., $L_i^{(4)AB}(\xi )$ with $i=1,\dots ,5$), the edge BC of third degree (i.e., $L_j^{(3)BC}(\eta )$ with $j=1,\dots ,4$), the edge DC of sixth degree, whereas the edge AD needs a polynomial of fifth degree.

As usual, we distinguish three categories of nodes, i.e. those at corners, those at intermediate boundary nodes, and finally those in the interior. The totality of the basis functions is shown in Equation (47), where $L$ represent the Lagrange polynomials performing as trial functions.

When working with Lagrange polynomials for blending and trial functions, not only the degree of the polynomial but also the location of the involved nodes plays a certain role. We recall that in the previous example of 25-node element (section 7.1), although the inner nodes were uniformly distributed, the nonuniform location of the nodes 2, 3, and 4 was introduced in the set of the five Lagrange polynomials along the bottom edge AB. The same holds true in the 27-node element under consideration:

$$\begin{array}{l}Boundary\hspace{0.17em}nodes:\\ N_1(\xi ,\eta )=E_1(\xi )L_1^{AD}(\eta )+E_1(\eta ){\bcancel{L}}_1^{AB}(\xi )-{\bcancel{E}}_1(\xi )E_1(\eta ), \\ N_2(\xi ,\eta )=E_1(\eta )L_2^{AB}(\xi ),\\ N_3(\xi ,\eta )=E_1(\eta )L_3^{AB}(\xi ),\\ N_4(\xi ,\eta )=E_1(\eta )L_4^{AB}(\xi ),\\ N_5(\xi ,\eta )=E_5(\xi )L_1^{BC}(\eta )+E_1(\eta ){\bcancel{L}}_5^{AB}(\xi )-{\bcancel{E}}_5(\xi )E_1(\eta ), \\ N_6(\xi ,\eta )=E_5(\xi )L_2^{BC}(\eta ),\\ N_7(\xi ,\eta )=E_5(\xi )L_3^{BC}(\eta ),\\ N_8(\xi ,\eta )=E_5(\xi )L_4^{BC}(\eta )+E_5(\eta )L_7^{DC}(\xi )-E_5(\xi )E_5(\eta ),\\ N_9(\xi ,\eta )=E_5(\eta )L_6^{DC}(\xi ),\\ N_{10}(\xi ,\eta )=E_5(\eta )L_5^{DC}(\xi ),\\ N_{11}(\xi ,\eta )=E_5(\eta )L_4^{DC}(\xi ),\\ N_{12}(\xi ,\eta )=E_5(\eta )L_3^{DC}(\xi ),\\ N_{13}(\xi ,\eta )=E_5(\eta )L_2^{DC}(\xi ),\\ N_{14}(\xi ,\eta )=E_1(\xi )L_6^{AD}(\eta )+E_5(\eta )L_1^{DC}(\xi )-E_1(\xi )E_5(\eta ), \\ N_{15}(\xi ,\eta )=E_1(\xi )L_5^{AD}(\eta ),\\ N_{16}(\xi ,\eta )=E_1(\xi )L_4^{AD}(\eta ),\\ N_{17}(\xi ,\eta )=E_1(\xi )L_3^{AD}(\eta ),\\ N_{18}(\xi ,\eta )=E_1(\xi )L_2^{AD}(\eta ),\\ Inner\hspace{0.17em}nodes:\\ N_{19}(\xi ,\eta )=E_2(\xi )E_2(\eta ),\\ N_{20}(\xi ,\eta )=E_3(\xi )E_2(\eta ),\\ N_{21}(\xi ,\eta )=E_4(\xi )E_2(\eta ),\\ N_{22}(\xi ,\eta )=E_2(\xi )E_3(\eta ),\\ N_{23}(\xi ,\eta )=E_3(\xi )E_3(\eta ),\\ N_{24}(\xi ,\eta )=E_4(\xi )E_3(\eta ),\\ N_{25}(\xi ,\eta )=E_2(\xi )E_4(\eta ),\\ N_{26}(\xi ,\eta )=E_3(\xi )E_4(\eta ),\\ N_{27}(\xi ,\eta )=E_4(\xi )E_4(\eta ).\end{array}$$

(47)

In contrast, when working with Bernstein polynomials, only the polynomial degrees are of interest, since the nodal points do not contribute to the formulation of the basis functions.

Focusing on Figure 8a, the blending functions are polynomials of degree four, either uniform Lagrange or Bernstein, accordingly.

Moreover, in both formulations (i.e., Lagrange and Bernstein), the corners nodes 1 and 5 eventually appear only in the $P_{\xi }$ projection, whereas the rest corner nodes 8 and 14 are influenced by all the three projections.

The result shown in Figure 8b has been obtained using Lagrange polynomials as blending functions of fourth degree, and for the edges it is as follows: The edges AB (bottom), BC, CD and DA are interpolated by polynomials or degree 4, 3, 6, and 5, respectively.

Regarding the Bernstein formulation, the bivariate basis functions are given by replacing the Lagrange with Bernstein polynomials of the same degree in Equation (47) and are illustrated in Figure 8(c).

It is also noted that, in both formulations, the determinant of the Jacobian takes everywhere its accurate value, and the partition of unity holds everywhere.

Now, setting the values of the function $u(x,y)=e^x\cos (\pi y)$ at the 27 uniform places on the boundary and the interior (nodes 1 to 27), the 27 coefficients α of the distrorted tensor-product Bernstein-Bézier element were calculated. Based on these two sets of 27 values separately (U or α values), a small difference was found for the L2-error norm within the square domain (Lagrange formulation according to Equation (47): L2 = 0.6408%, Bernstein formulation: L2 = 0.6154%).

8. Mesh Refinement

There are many engineering problems in which the exact solution varies rapidly at a certain point of the domain. Due to this high gradient, it becomes necessary to refine the mesh at this point of the domain, but then the problem of hanging nodes appears.

A simple method to overcome this problem is to subdivide the initial patch ABCD, using progressively one-fourth of the previous subdomain closest to the ‘singularity’ at the corner D, as shown in Figure 9. One may observe the hanging nodes 5 and 10 after the first refinement, as well as 15 and 18 after the second refinement. Below, we will see that in the transfinite formulation they do not constitute a difficulty to deal with.

Therefore, starting from the classical quadratic element shown in Figure 9a, we subdivide the upper left fourth, as shown in Figure 9b. In this case, we have four vertical stations at $\xi =0,\frac{1}{4},\frac{1}{2},1$ in the horizontal and $\eta =0,\frac{1}{2},\frac{3}{4},1$ in the vertical direction.

As was shown in the previous sections, we are flexible in choosing the blending and trial functions. Below we present the global shape functions when Lagrange polynomials have been implemented everywhere. Then, the blending functions are Lagrange polynomials of third degree but are based on non-uniform nodes (not the same in the two directions, and thus the subscripts denote the position of each associated station):

x-direction: $E_0(\xi ),E_{\frac{1}{4}}(\xi ),E_{\frac{2}{4}}(\xi ),E_1(\xi )$

y-direction: $E_0(\eta ),E_{\frac{2}{4}}(\eta ),E_{\frac{3}{4}}(\eta ),E_1(\eta )$

Regarding the trial functions, we distinguish three sets of Lagrange polynomials, as follows:

- Those along the boundary stations AB and BC, which are uniform quadratic polynomials $(p=2)$: $L_0^{(2)}(s)=(2s-1)(s-1)$, $L_{\frac{1}{2}}^{(2)}(s)=4s(1-s)$, and $L_1^{(2)}(s)=s(2s-1)$, where $s$ denotes either of $\xi$ or $\eta$.

- Those along the horizontal stations 4-5-6-7 and 11-12-13-14, which are non-uniform Lagrange polynomials, identical with the blending functions toward the $\xi$-axis: $L_0^{(3)}(\xi )=E_0(\xi ),$$L_{\frac{1}{4}}^{(3)}(\xi )=E_{\frac{1}{4}}(\xi ),$ $L_{\frac{2}{4}}^{(3)}(\xi )=E_{\frac{2}{4}}(\xi ),$ and $L_1^{(3)}(\xi )=E_1(\xi )$.

- Those along the vertical stations 1-4-8-11 and 2-6-10-13, which are identical with the blending functions toward the $\eta$-axis: $L_0^{(3)}(\eta )=E_0(\eta ),$ $L_{\frac{2}{4}}^{(3)}(\eta )=E_{\frac{2}{4}}(\eta ),$ $L_{\frac{3}{4}}^{(3)}(\eta )=E_{\frac{3}{4}}(\eta ),$ and $L_1^{(3)}(\eta )=E_1(\eta )$.

Based on the above three functional sets of trial functions, the application of the Boolean sum in conjunction with artificial nodes which are the intersections of the line 9-5 with 1-2 and the line 9-10 with 7-14, leads to the following global shape functions for the entire patch (for details, see Appendix C):

$$\begin{array}{l}{N}_1(\xi ,\eta )=E_0(\xi )L_0^{(3)}(\eta )+E_0(\eta )L_0^{(2)}(\xi )-E_0(\xi )E_0(\eta ),\\ N_2(\xi ,\eta )=E_{\frac{2}{4}}(\xi )L_0^{(3)}(\eta )+E_0(\eta )L_{\frac{1}{2}}^{(2)}(\xi )-E_{\frac{2}{4}}(\xi )E_0(\eta ),\\ N_3(\xi ,\eta )=E_1(\xi )L_0^{(2)}(\eta )+E_0(\eta )L_1^{(2)}(\xi )-E_1(\xi )E_0(\eta ),\\ N_4(\xi ,\eta )=E_0(\xi )L_{\frac{2}{4}}^{(3)}(\eta )+E_{\frac{2}{4}}(\eta )L_0^{(3)}(\xi )-E_0(\xi )E_{\frac{2}{4}}(\eta ),\\ N_5(\xi ,\eta )=E_{\frac{1}{4}}(\xi )L_{\frac{2}{4}}^{(3)}(\eta )+E_{\frac{2}{4}}(\eta )L_{\frac{1}{4}}^{(3)}(\xi )-E_{\frac{1}{4}}(\xi )E_{\frac{2}{4}}(\eta ),\\ N_6(\xi ,\eta )=E_{\frac{2}{4}}(\xi )L_{\frac{2}{4}}^{(3)}(\eta )+E_{\frac{2}{4}}(\eta )L_{\frac{2}{4}}^{(3)}(\xi )-E_{\frac{2}{4}}(\xi )E_{\frac{2}{4}}(\eta ),\\ N_7(\xi ,\eta )=E_1(\xi )L_{\frac{1}{2}}^{(2)}(\eta )+E_{\frac{2}{4}}(\eta )L_1^{(3)}(\xi )-E_1(\xi )E_{\frac{2}{4}}(\eta ),\\ N_8(\xi ,\eta )=E_0(\xi )L_{\frac{3}{4}}^{(3)}(\eta )+E_{\frac{3}{4}}(\eta )L_0^{(3)}(\xi )-E_0(\xi )E_{\frac{3}{4}}(\eta ),\\ N_9(\xi ,\eta )=E_{\frac{1}{4}}(\xi )L_{\frac{3}{4}}^{(3)}(\eta )+E_{\frac{3}{4}}(\eta )L_{\frac{1}{4}}^{(3)}(\xi )-E_{\frac{1}{4}}(\xi )E_{\frac{3}{4}}(\eta ),\\ N_{10}(\xi ,\eta )=E_{\frac{2}{4}}(\xi )L_{\frac{3}{4}}^{(3)}(\eta )+E_{\frac{3}{4}}(\eta )L_{\frac{2}{4}}^{(3)}(\xi )-E_{\frac{2}{4}}(\xi )E_{\frac{3}{4}}(\eta ),\\ N_{11}(\xi ,\eta )=E_0(\xi )L_1^{(3)}(\eta )+E_1(\eta )L_0^{(3)}(\xi )-E_0(\xi )E_1(\eta ),\\ N_{12}(\xi ,\eta )=E_{\frac{1}{4}}(\xi )L_1^{(3)}(\eta )+E_1(\eta )L_{\frac{1}{4}}^{(3)}(\xi )-E_{\frac{1}{4}}(\xi )E_1(\eta ),\\ N_{13}(\xi ,\eta )=E_{\frac{2}{4}}(\xi )L_1^{(3)}(\eta )+E_1(\eta )L_{\frac{2}{4}}^{(3)}(\xi )-E_{\frac{2}{4}}(\xi )E_1(\eta ),\\ N_{14}(\xi ,\eta )=E_1(\xi )L_1^{(2)}(\eta )+E_1(\eta )L_1^{(3)}(\xi )-E_1(\xi )E_1(\eta ).\end{array}$$

(48)

The shape functions of Equation (48) are illustrated in Figure 10a, whereas after replacing the Lagrange with Bernstein polynomials, the updated basis functions are shown in Figure 10b.

Considering the almost singular function $u(x,y)=1/\ln (r+0.1)$, in which the position $r$ is measured from the top left corner D (node 11) of Figure 9, for both formulations (Lagrange and Bernstein polynomials), the 9-node element gave L2-error norm equal 10.7304%, whereas the 14-node decreased to 7.6496%.

Moreover, following the same procedure as above, the 19-node element which is shown in Figure 9(c), can be treated through the global basis functions:

$$\begin{array}{l}{N}_1(\xi ,\eta )=E_0(\xi )L_0^{(4)}(\eta )+E_0(\eta )L_0^{(2)}(\xi )-E_0(\xi )E_0(\eta ),\\ N_2(\xi ,\eta )=E_{\frac{4}{8}}(\xi )L_0^{(3)}(\eta )+E_0(\eta )L_{\frac{1}{2}}^{(2)}(\xi )-E_{\frac{4}{8}}(\xi )E_0(\eta ),\\ N_3(\xi ,\eta )=E_1(\xi )L_0^{(2)}(\eta )+E_0(\eta )L_1^{(2)}(\xi )-E_1(\xi )E_0(\eta ), \\ N_4(\xi ,\eta )=E_0(\xi )L_{\frac{4}{8}}^{(4)}(\eta )+E_{\frac{4}{8}}(\eta )L_0^{(3)}(\xi )-E_0(\xi )E_{\frac{4}{8}}(\eta ),\\ N_5(\xi ,\eta )=E_{\frac{2}{8}}(\xi )L_{\frac{4}{8}}^{(4)}(\eta )+E_{\frac{4}{8}}(\eta )L_{\frac{1}{4}}^{(3)}(\xi )-E_{\frac{2}{8}}(\xi )E_{\frac{4}{8}}(\eta ),\\ N_6(\xi ,\eta )=E_{\frac{4}{8}}(\xi )L_{\frac{4}{8}}^{(3)}(\eta )+E_{\frac{4}{8}}(\eta )L_{\frac{2}{4}}^{(3)}(\xi )-E_{\frac{4}{8}}(\xi )E_{\frac{4}{8}}(\eta ),\\ N_7(\xi ,\eta )=E_1(\xi )L_{\frac{1}{2}}^{(2)}(\eta )+E_{\frac{4}{8}}(\eta )L_1^{(3)}(\xi )-E_1(\xi )E_{\frac{4}{8}}(\eta ),\\ N_8(\xi ,\eta )=E_0(\xi )L_{\frac{6}{8}}^{(4)}(\eta )+E_{\frac{6}{8}}(\eta )L_0^{(4)}(\xi )-E_0(\xi )E_{\frac{6}{8}}(\eta ),\\ N_9(\xi ,\eta )=E_{\frac{2}{8}}(\xi )L_{\frac{6}{8}}^{(4)}(\eta )+E_{\frac{6}{8}}(\eta )L_{\frac{2}{8}}^{(4)}(\xi )-E_{\frac{2}{8}}(\xi )E_{\frac{6}{8}}(\eta ),\\ N_{10}(\xi ,\eta )=E_{\frac{4}{8}}(\xi )L_{\frac{6}{8}}^{(3)}(\eta )+E_{\frac{6}{8}}(\eta )L_{\frac{4}{8}}^{(4)}(\xi )-E_{\frac{4}{8}}(\xi )E_{\frac{6}{8}}(\eta ),\\ N_{11}(\xi ,\eta )=E_0(\xi )L_1^{(4)}(\eta )+E_1(\eta )L_0^{(4)}(\xi )-E_0(\xi )E_1(\eta ),\\ N_{12}(\xi ,\eta )=E_{\frac{2}{8}}(\xi )L_1^{(4)}(\eta )+E_1(\eta )L_{\frac{2}{8}}^{(4)}(\xi )-E_{\frac{2}{8}}(\xi )E_1(\eta ),\\ N_{13}(\xi ,\eta )=E_{\frac{4}{8}}(\xi )L_1^{(3)}(\eta )+E_1(\eta )L_{\frac{4}{8}}^{(4)}(\xi )-E_{\frac{4}{8}}(\xi )E_1(\eta ),\\ N_{14}(\xi ,\eta )=E_1(\xi )L_1^{(2)}(\eta )+E_1(\eta )L_1^{(4)}(\xi )-E_1(\xi )E_1(\eta ),\\ N_{15}(\xi ,\eta )=E_{\frac{1}{8}}(\xi )L_{\frac{6}{8}}^{(4)}(\eta )+E_{\frac{6}{8}}(\eta )L_{\frac{1}{8}}^{(4)}(\xi )-E_{\frac{1}{8}}(\xi )E_{\frac{6}{8}}(\eta ),\\ N_{16}(\xi ,\eta )=E_0(\xi )L_{\frac{7}{8}}^{(4)}(\eta )+E_{\frac{7}{8}}(\eta )L_0^{(4)}(\xi )-E_0(\xi )E_{\frac{7}{8}}(\eta ),\\ N_{17}(\xi ,\eta )=E_{\frac{1}{8}}(\xi )L_{\frac{7}{8}}^{(4)}(\eta )+E_{\frac{7}{8}}(\eta )L_{\frac{1}{8}}^{(4)}(\xi )-E_{\frac{1}{8}}(\xi )E_{\frac{7}{8}}(\eta ),\\ N_{18}(\xi ,\eta )=E_{\frac{2}{8}}(\xi )L_{\frac{7}{8}}^{(4)}(\eta )+E_{\frac{7}{8}}(\eta )L_{\frac{2}{8}}^{(4)}(\xi )-E_{\frac{2}{8}}(\xi )E_{\frac{7}{8}}(\eta ),\\ N_{19}(\xi ,\eta )=E_{\frac{1}{8}}(\xi )L_1^{(4)}(\eta )+E_1(\eta )L_{\frac{1}{8}}^{(4)}(\xi )-E_{\frac{1}{8}}(\xi )E_1(\eta ).\end{array}$$

(49)

9. T-mesh like transfinite elements

The capability of transfinite elements to deal with mesh refinement is encouraging, and therefore it is important to check whether it can also treat unstructured domains in the form of T-mesh patches.

EXAMPLE-1: Let us consider the simple T-mesh shown in Figure 11. Below we test and compare several cases, such as polynomials (Lagrange and Bernstein) as well as cubic B-splines.

9.1. Lagrange and Bernstein polynomials

Regarding the application of the conventional transfinite interpolation, one may recognize in Figure 11(a) five vertical and another five horizontal uniform stations. Therefore, each of them will be polynomial of the fourth degree ($p_x=p_y=4$).

To build the tree projectors, again we need to complete the tensor product by adding the two points (G and F) shown in Figure 11(b). Each of the three edges (BC, DC, and AD) contain five nodes and thus are represented by quartic polynomials, whereas the bottom edge AB includes only four nodes (1, 2 , 3, and 4). Clearly, the function U(ξ,0) along AB will be represented using a non-uniform Lagrange polynomial of degree three, and this affects the projector P_η. In contrast, the auxiliary point G (between 2 and 3) is involved in the interpolation of the function U(ξ = 1/2, η) which participates as a third additive term in the projector P_ξ.

The complete procedure will be given below, but the symmetry of the mesh in Figure 11 with respect to the y-axis dictates, in advance, that the transfinite interpolation degenerates to only the vertical projector P_η.

Actually, in more detail we have:

$$\begin{array}{ll}{P}_{\xi }& =E_0(\xi )\left[L_0^{(4)}(\eta )U_1+L_{\frac{1}{4}}^{(4)}(\eta )U_5+L_{\frac{2}{4}}^{(4)}(\eta )U_9+L_{\frac{3}{4}}^{(4)}(\eta )U_{14}+L_1^{(4)}(\eta )U_{19}\right]\\ & +E_{\frac{1}{4}}(\xi )\left[L_0^{(4)}(\eta )U_2+L_{\frac{1}{4}}^{(4)}(\eta )U_6+L_{\frac{2}{4}}^{(4)}(\eta )U_{10}+L_{\frac{3}{4}}^{(4)}(\eta )U_{15}+L_1^{(4)}(\eta )U_{20}\right]\\ & +E_{\frac{2}{4}}(\xi )\left[L_0^{(4)}(\eta )U_G+L_{\frac{1}{4}}^{(4)}(\eta )U_H+L_{\frac{2}{4}}^{(4)}(\eta )U_{11}+L_{\frac{3}{4}}^{(4)}(\eta )U_{16}+L_1^{(4)}(\eta )U_{21}\right]\\ & +E_{\frac{3}{4}}(\xi )\left[L_0^{(4)}(\eta )U_3+L_{\frac{1}{4}}^{(4)}(\eta )U_7+L_{\frac{2}{4}}^{(4)}(\eta )U_{12}+L_{\frac{3}{4}}^{(4)}(\eta )U_{17}+L_1^{(4)}(\eta )U_{22}\right]\\ & +E_1(\xi )\left[L_0^{(4)}(\eta )U_4+L_{\frac{1}{4}}^{(4)}(\eta )U_8+L_{\frac{2}{4}}^{(4)}(\eta )U_{13}+L_{\frac{3}{4}}^{(4)}(\eta )U_{18}+L_1^{(4)}(\eta )U_{23}\right],\end{array}$$

(50)

also

$$\begin{array}{ll}{P}_{\eta }& =E_0(\eta ){\left[L_0^{(3)AB}(\xi )U_1+L_{\frac{1}{4}}^{(3)AB}(\xi )U_2+L_{\frac{2}{4}}^{(3)AB}(\xi )U_3+L_1^{(3)AB}(\xi )U_4\right]}_{BC}\\ & +E_{\frac{1}{4}}(\eta )\left[L_0^{(3)\eta =1/4}(\xi )U_5+L_{\frac{1}{4}}^{(3)\eta =1/4}(\xi )U_6+L_{\frac{2}{4}}^{(3)\eta =1/4}(\eta )U_7+L_1^{(3)\eta =1/4}(\eta )U_8\right]\\ & +E_{\frac{2}{4}}(\eta )\left[L_0^{(4)}(\eta )U_9+L_{\frac{1}{4}}^{(4)}(\eta )U_{10}+L_{\frac{2}{4}}^{(4)}(\eta )U_{11}+L_{\frac{3}{4}}^{(4)}(\eta )U_{12}+L_1^{(4)}(\eta )U_{13}\right]\\ & +E_{\frac{3}{4}}(\eta )\left[L_0^{(4)}(\eta )U_{14}+L_{\frac{1}{4}}^{(4)}(\eta )U_{15}+L_{\frac{2}{4}}^{(4)}(\eta )U_{16}+L_{\frac{3}{4}}^{(4)}(\eta )U_{17}+L_1^{(4)}(\eta )U_{18}\right]\\ & +E_1(\eta )\left[L_0^{(4)}(\eta )U_{19}+L_{\frac{1}{4}}^{(4)}(\eta )U_{20}+L_{\frac{2}{4}}^{(4)}(\eta )U_{21}+L_{\frac{3}{4}}^{(4)}(\eta )U_{22}+L_1^{(4)}(\eta )U_{23}\right],\end{array}$$

(51)

whereas the tensor product is:

$$\begin{array}{ll}{P}_{\xi \eta }& =E_0(\xi )E_0(\eta )U_1+E_{\frac{1}{4}}(\xi )E_0(\eta )U_2+E_{\frac{2}{4}}(\xi )E_0(\eta )U_G+E_{\frac{3}{4}}(\xi )E_0(\eta )U_3+E_1(\xi )E_0(\eta )U_4\\ & +E_0(\xi )E_{\frac{1}{4}}(\eta )U_5+E_{\frac{1}{4}}(\xi )E_{\frac{1}{4}}(\eta )U_6+E_{\frac{2}{4}}(\xi )E_{\frac{1}{4}}(\eta )U_H+E_{\frac{3}{4}}(\xi )E_{\frac{1}{4}}(\eta )U_7+E_1(\xi )E_{\frac{1}{4}}(\eta )U_8\\ & +E_0(\xi )E_{\frac{2}{4}}(\eta )U_9+E_{\frac{1}{4}}(\xi )E_{\frac{2}{4}}(\eta )U_{10}+E_{\frac{2}{4}}(\xi )E_{\frac{2}{4}}(\eta )U_{11}+E_{\frac{3}{4}}(\xi )E_{\frac{2}{4}}(\eta )U_{12}+E_1(\xi )E_{\frac{2}{4}}(\eta )U_{13}\\ & +E_0(\xi )E_{\frac{3}{4}}(\eta )U_{14}+E_{\frac{1}{4}}(\xi )E_{\frac{3}{4}}(\eta )U_{15}+E_{\frac{2}{4}}(\xi )E_{\frac{3}{4}}(\eta )U_{16}+E_{\frac{3}{4}}(\xi )E_{\frac{3}{4}}(\eta )U_{17}+E_1(\xi )E_{\frac{3}{4}}(\eta )U_{18}\\ & +E_0(\xi )E_1(\eta )U_{19}+E_{\frac{1}{4}}(\xi )E_1(\eta )U_{20}+E_{\frac{2}{4}}(\xi )E_1(\eta )U_{21}+E_{\frac{3}{4}}(\xi )E_1(\eta )U_{22}+E_1(\xi )E_1(\eta )U_{23}.\end{array}$$

(52)

Since the trial functions $L_i^{(4)}(\eta )$ are identical with the blending functions, i.e. $L_i^{(4)}(\eta )=E_i(\eta )$, one may observe that in this example we have:

$$P_{\xi }=P_{\xi \eta },$$

(53)

whence

$$\begin{array}{ll}U(\xi ,\eta )& ={\bcancel{P}}_{\xi }+P_{\eta }-{\bcancel{P}}_{\xi \eta }\\ & =P_{\eta }.\end{array}$$

(54)

The cancelation of the difference $(P_{\xi }-P_{\xi \eta })$ also means the cancelation of the involved auxiliary variables ($U_G$ and $U_H$). Eventually, the factors of the 23 nodal values $U_i,i=1,\dots ,23$ in Equation (49) denote the desired shape functions of the macroelement associated with the patch ABCD, which are illustrated in Figure 12(a) using either Lagrange polynomials. When the Lagrange polynomials are replaced by Bernstein ones, we receive the set of 23 basis functions shown in Figure 12(b).

The quality of the above obtained 23 shape functions is evaluated as follows. Setting the nodal values equal to the function $f(x,y)=e^x\cos (\pi y)$, then the L₂-norm was found equal to 0.4968% (in both Lagrange and Bernstein polynomials).

9.2. B-spline blending and trial functions

As an innovative step beyond state-of-the-art T-splines, we consider that the blending functions are B-splines and the same holds for the trial functions as well.

9.2.1 Example 9-1

In this context, the existence of five stations per direction (see Figure 11 and its analogue Figure 13) is consistent with a univariate B-spline formulation with five control points which span the entire domain [0,1]. This can be accomplished as follows:

• Quadratic B-spline (p = 2) in conjunction with knot vector U = [0,0,0, 1/3,2/3, 1,1,1];
• Cubic B-spline (p = 3) in conjunction with knot vector U = [0,0,0,0, 1/2, 1,1,1,1].

The results for the L2-error norm are shown in Table 3.

Clearly, cases 7 and 8 of Table 3 were obtained using the transfinite interpolation formula given by Equation (51), in which the Lagrange polynomials were replaced by B-splines. A question here arises whether the so produced basis functions are related to those that can be obtained through the conventional rule of T-splines. To answer this question, we form the arrays (each of size 23 × 5) which for each anchor provide the local knot vectors for the T-mesh shown in Figure 13 ($U_{\xi }$ toward ξ-direction, $V_{\eta }$ toward η-direction):

$$U_{\xi }=\left[\begin{array}{ccccc}0& 0& 0& 0& 1\\ 0& 0& 0& 1& 1\\ 0& 0& 1& 1& 1\\ 0& 1& 1& 1& 1\\ 0& 0& 0& 0& 1\\ 0& 0& 0& 1& 1\\ 0& 0& 1& 1& 1\\ 0& 1& 1& 1& 1\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0.5& 1& 1& 1& 1\end{array}\right],\hspace{1em}{V}_{\eta }=\left[\begin{array}{ccccc}0& 0& 0& 0& 0.5\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0& 0.5\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0& 0.5& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0& 0.5& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0& 0.5& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0.5& 1& 1& 1& 1\\ 0.5& 1& 1& 1& 1\end{array}\right]$$

(55)

One may observe in Figure 14 that the basis functions which were produced through transfinite interpolation using cubic B-splines for both blending and trial functions, are identical with those obtained using the classical T-spline formulation based on the local knot vectors given by Equation (55). Numerically, the differences were less than 10^-16.

EXAMPLE-2: Missing boundary knots

In Figure 15(a,b) we present a trivial case of 14 nodes where the nodal points G and F are missing to complete the tensor product. The four stations per direction suggest blending functions of degree three. The third vertical station measured from left requires the artificial points G and interpolation of third degree. Similarly, the third horizontal station measured from bottom requires the artificial point F and polynomial of third degree. Using Lagrange polynomials, the three projectors and their Boolean sum are given below.

$$\begin{array}{ll}{P}_{\xi }& =E_0(\xi )\left[E_0(\eta )U_1+E_{\frac{1}{3}}(\eta )U_4+E_{\frac{2}{3}}(\eta )U_8+E_1(\eta )U_{11}\right]\\ & +E_{\frac{1}{3}}(\xi )\left[E_0(\eta )U_2+E_{\frac{1}{3}}(\eta )U_5+E_{\frac{2}{3}}(\eta )U_9+E_1(\eta )U_{12}\right]\\ & +E_{\frac{2}{3}}(\xi )\left[E_0(\eta )U_G+E_{\frac{1}{3}}(\eta )U_6+E_{\frac{2}{3}}(\eta )U_{10}+E_1(\eta )U_{13}\right]\\ & +E_1(\xi ){\left[L_0^{(2)BC}(\eta )U_3+L_{\frac{1}{3}}^{(2)BC}(\eta )U_7+L_1^{(2)BC}(\eta )U_{14}\right]}_{BC},\end{array}$$

(56)

also

$$\begin{array}{ll}{P}_{\eta }& =E_0(\eta ){\left[L_0^{(2)AB}(\xi )U_1+L_{\frac{1}{3}}^{(2)AB}(\xi )U_2+L_1^{(2)AB}(\xi )U_3\right]}_{BC}\\ & +E_{\frac{1}{3}}(\eta )\left[E_0(\eta )U_4+E_{\frac{1}{3}}(\eta )U_5+E_{\frac{2}{3}}(\eta )U_6+E_1(\eta )U_7\right]\\ & +E_{\frac{2}{3}}(\eta )\left[E_0(\eta )U_8+E_{\frac{1}{3}}(\eta )U_9+E_{\frac{2}{3}}(\eta )U_{10}+E_1(\eta )U_F\right]\\ & +E_1(\eta )\left[E_0(\eta )U_{11}+E_{\frac{1}{3}}(\eta )U_{12}+E_{\frac{2}{3}}(\eta )U_{13}+E_1(\eta )U_{14}\right],\end{array}$$

(57)

whereas the tensor product is:

$$\begin{array}{ll}{P}_{\xi \eta }& =E_0(\xi )E_0(\eta )U_1+E_{\frac{1}{3}}(\xi )E_0(\eta )U_2+E_{\frac{2}{3}}(\xi )E_0(\eta )U_G+E_1(\xi )E_0(\eta )U_3\\ & +E_0(\xi )E_{\frac{1}{3}}(\eta )U_4+E_{\frac{1}{3}}(\xi )E_{\frac{1}{3}}(\eta )U_5+E_{\frac{2}{3}}(\xi )E_{\frac{1}{3}}(\eta )U_6+E_1(\xi )E_{\frac{1}{3}}(\eta )U_7\\ & +E_0(\xi )E_{\frac{2}{3}}(\eta )U_8+E_{\frac{1}{3}}(\xi )E_{\frac{2}{3}}(\eta )U_9+E_{\frac{2}{3}}(\xi )E_{\frac{2}{3}}(\eta )U_{10}+E_1(\xi )E_{\frac{2}{3}}(\eta )U_F\\ & +E_0(\xi )E_1(\eta )U_{11}+E_{\frac{1}{3}}(\xi )E_1(\eta )U_{12}+E_{\frac{2}{3}}(\xi )E_1(\eta )U_{13}+E_1(\xi )E_1(\eta )U_{14}.\end{array}$$

(58)

One may observe that:

• The interpolation along the right vertical edge BC is fully preserved using non-uniform polynomials of degree 2, and this justifies the subscript BC (boundary conditions) in the last square bracket of the projector $P_{\xi }$ (see Equation (56)).
• The interpolation along the bottom edge AB is fully preserved using non-uniform polynomials of degree 2, and this justifies the subscript BC (boundary conditions) in the first square bracket of the projector $P_{\eta }$ (see Equation (57)).
• The artificial values ($U_G,U_F$) appear in the projectors $P_{\xi }$ and $P_{\eta }$, respectively, whereas both of them appear in the tensor product $P_{\xi \eta }$. As a result, they are eventually eliminated.
• After performing the Boolean sum $P_{\xi }\otimes P_{\eta }=P_{\xi }+P_{\eta }-P_{\xi \eta }$, the transfinite interpolation offers the following 14 global shape functions in the domain $[0,1]\times [0,1]$:

  $$\begin{array}{ll}U(\xi ,\eta )=& \left[L_0^{(2)AB}(\xi )E_0(\eta )\right]\cdot U_1\\ & +\left[L_{\frac{1}{3}}^{(2)AB}(\xi )E_0(\eta )\right]\cdot U_2\\ & +\left[E_1(\xi )L_0^{(2)AB}(\eta )+L_1^{(2)AB}(\xi )E_0(\eta )-E_1(\xi )E_0(\eta )\right]\cdot U_3\\ & +\left[E_0(\xi )E_{\frac{1}{3}}(\eta )\right]\cdot U_4\\ & +\left[E_{\frac{1}{3}}(\xi )E_{\frac{1}{3}}(\eta )\right]\cdot U_5\\ & +\left[E_{\frac{2}{3}}(\xi )E_{\frac{1}{3}}(\eta )\right]\cdot U_6\\ & +\left[E_1(\xi )L_{\frac{1}{3}}^{(2)BC}(\eta )\right]\cdot U_7\\ & +\left[E_0(\xi )E_{\frac{2}{3}}(\eta )\right]\cdot U_8\\ & +\left[E_{\frac{1}{3}}(\xi )E_{\frac{2}{3}}(\eta )\right]\cdot U_9\\ & +\left[E_{\frac{2}{3}}(\xi )E_{\frac{2}{3}}(\eta )\right]\cdot U_{10}\\ & +\left[E_0(\xi )E_1(\eta )\right]\cdot U_{11}\\ & +\left[E_{\frac{1}{3}}(\xi )E_1(\eta )\right]U_{12}\\ & +\left[E_{\frac{2}{3}}(\xi )E_1(\eta )\right]\cdot U_{13}\\ & +\left[E_1(\xi )L_1^{(2)BC}(\eta )\right]\cdot U_{14}\\ & ={\displaystyle \sum _{j=1}^{14}{\varphi }_j(\xi ,\eta )U_j}.\end{array}$$

  (59)

The fourteen shape functions of Equation (59) are illustrated in Figure 16(a), and thus they obtain positive and negative values. Although at the nodes they take the value of unity, in between them they may exceed the unity.

Furthermore, if blending and trial functions are Bernstein polynomials of the same degree as the previously mentioned Lagrange polynomials, the form of Equation (59) is again derived, and the updated results are illustrated in Figure 16(b). One may recognize the Bernstein polynomials on the boundaries and may observe only positive values.

It is noted that in both cases, i.e. using either Lagrange of Bernstein polynomials, the partition of unity is fulfilled:

$${\displaystyle \sum _{j=1}^{14}{\varphi }_j(\xi },\eta )=1,\hspace{1em}\forall (\xi ,\eta )\in {\left[0,1\right]}^2.$$

(60)

However, to obtain the above Partition of Unity Property (from which the transfinite interpolation originated), we must consider polynomials of second degree as well, which are denoted in Equation (59) with red color.

Interestingly, when the T-spline local vectors were calculated only in terms of degree p = 3, the T-spline formulation failed to ensure unity basis functions at the corner B and C, simply because there were not four discrete nodes therein. In contrast, the transfinite interpolation could treat this issue in an effective way, according to Equation (59).

10. Discussion

Currently, the tensor product is probably the most usual formula for interpolation and approximation in well-shaped patches, although in several scientific fields such as geo-informatics there are more complex interpolations which allow for global or local control [27-29].

In computational mechanics, one disadvantage of the tensor product is that it requires the same number and relative location of nodes on the opposite edges of a large element or patch, with the internal nodes which should be along discrete isolines passing through the boundary nodes. This in turn leads to connectivity problems when two dissimilar patches must be joined together. As a result, the buildup of “transition elements”, which can ensure the inter-element compatibility, is imperative. Such elements can be automatically generated using the transfinite interpolation [30].

In addition to the inter-element compatibility, the existence of high gradients requires local mesh refinement leading to “hanging nodes”, but this issue is easily treated through the transfinite interpolation, using a larger number of stations close to the singularity. The hanging nodes are treated by implicitly extending the lines to which they belong until they reach the boundary. However, the produced nodes are ‘artificial’ (i.e. auxiliary), in the sense that they are eventually eliminated, and thus not involved in the final expressions of the basis functions.

Aiming at treating the above two needs, while most of the transfinite elements can be developed using Lagrange polynomials, when the number of the involved nodes in a station becomes high, the interpolation may be amenable to numerical oscillations, and therefore it may be useful to replace them with B-splines. In other words, if the number of stations is relatively small, the blending functions may be selected as Lagrange polynomials whereas the trial functions along the stations may be chosen to be B-splines or NURBS.

In addition to the Lagrange-polynomial based blending functions, it has been previously found that they can be chosen as any other kind of cardinal functions such as piecewise-linear [26], cosine-like [16](p.125), cardinal natural cubic B-spline, and others. For example, the choice of piecewise-linear blending functions in conjunction with piecewise-linear trial functions, lead to an assembly of conventional bilinear finite elements [26], and so on.

Recently, it was found that Lagrange polynomials may be replaced by Bernstein polynomials of the same degree, and this idea was combined with weights which make the transfinite formulation capable of accurately representing conics and spherical cups [17].

Beyond the above state-of-the-art, the present paper proposes a still further effective replacement of the higher order blending functions, not only by Bernstein polynomials as previously, but also by non-cardinal B-splines. From the theoretical point of view, in the light of Theorem-3 which states that every tensor product equals either of the three projectors in transfinite interpolation, the choice of Bernstein polynomials for both the blending and the trial functions, makes the transfinite interpolation degenerate to tensor-product Bernstein polynomials. Similarly, the choice of B-splines for both the blending and the trial functions, makes the transfinite interpolation degenerate to tensor-product B-splines. As a result of this finding, it was conjectured that B-splines can be generally applied as blending functions, despite the fact they are not cardinal (i.e., not of [0-1]-type).

From the practical point of view, the utilization of B-splines as blending functions assists in the buildup of many mixed transfinite elements, which may combine even different polynomial degrees between blending and trial functions, and in some cases (where some internal lines do not intersect the boundary) can treat patches like T-splines.

On-going research on the general T-spline reveals difficulties when there are many intersecting points, because not always all the artificial points are eliminated. This topic will be presented in full detail in a future paper.

11. Conclusions

Within a geometric curvilinear patch, the transfinite interpolation formula is applicable in conjunction with several sets of blending functions, which do not necessarily are of cardinal (i.e., [0-1])-type. In this sense, any tensor-product of well-known functions such as Lagrange and Bézier-Bernstein polynomials, B-splines, NURBS and others, may be considered under the umbrella (prism) of the transfinite interpolation. For example, the NURBS-based isogeometric analysis may be considered as a special case where both the blending and trial functions are chosen as weighted B-splines, while some aspects of T-splines can be derived from the general theory, which allows for the simultaneous utilization of different polynomial degrees. Although the present paper was limited to two dimensions under Dirichlet boundary conditions, it can be extended to any kind of mixed conditions and three-dimensional problems as well. The closed form of the global basis functions (or the easy generic production of them) makes the methodology applicable to the whole spectrum of computational mechanics as well as the numerical solution of boundary-value and eigenvalue problems governed by any type of partial differential equations.