Quadrilateral Irregular Network for Mesh-Based Interpolation

Abstract

Numerical analysis has been adopted in nearly all modern scientific and engineering fields due to the rapid and ongoing evolution of computational technology, with the number of grid or mesh points in a given data field also increasing. Some values must be extracted from large data fields to evaluate and supplement numerical analysis results and observational data, thereby highlighting the need for a fast and effective interpolation approach. The quadrilateral irregular network (QIN) proposed in this study is a fast and reliable interpolation method that is capable of sufficiently satisfying these demands. A comparative sensitivity analysis is first performed using known test functions to assess the accuracy and computational requirements of QIN relative to conventional interpolation methods. These same interpolation methods are then employed to produce simple numerical model results for a real-world comparison. Unlike conventional interpolation methods, QIN can obtain reliable results with a guaranteed degree of accuracy since there is no need to determine the optimal parameter values. Furthermore, QIN is a computationally efficient method compared with conventional interpolation methods that require the entire data space to be evaluated during interpolation, even if only a subset of the data space requires interpolation.

Introduction

Spatial interpolation is essential in many scientific and engineering fields to ensure the accurate and effective evaluation of physical data in a continuous domain. Numerous spatial interpolation techniques have been introduced, with these approaches possessing various degrees of complexity (Lam, 1983; Myers, 1994; Webster and Oliver, 2007; Li and Heap, 2008); however, selecting the best technique that produces reliable and accurate results is difficult. While various studies have attempted to construct mechanisms for choosing and evaluating the methods that best suit the actual data (Caruso and Quarta, 1998; Li and Heap, 2011), the various factors that affect the performance of spatial interpolations generally indicate that no method is always optimal for specific situations. A method that works well with one dataset may be unsuitable for a different dataset because the interpolated results are strongly dependent on the characteristics of the data. Therefore, comparative studies have recently focused on identifying the optimal method for specific objects of a given study. For example, Sun et al. (2009) compared interpolation methods using data obtained over 22 years on the temporal and spatial variations of groundwater depth in a specific study area. Tabios III and Salas (1985), Keblouti et al. (2012), and Wagner et al. (2012) evaluated various methods to find the most appropriate interpolation method for sparse precipitation data. Piazza et al. (2011) conducted a comparative study on various methods for the estimation of missing data in rainfall time series. Teegavarapu et al. (2012) developed and estimated spatial interpolation weighting methods for the transformation of hourly multisensory precipitation data from one grid to another with different sizes and orientations. Rap et al. (2009) compared interpolation methods for aerosol-cloud parameterization in global climate modeling.

As computational techniques have improved, numerical simulations have grown in popularity in many disciplines and are used as preliminary or predictive tools for specific purposes. Continuous natural objects or characteristics, which are generally obtained from field measurements at scattered points, must be assigned to discrete nodal points for numerical modeling where interpolation is used for data pre-processing. An interpolation may also be required after the simulation to compare the calculated results with observed values. The number of nodal points has increased from dozens or hundreds to thousands or more as more detailed modeling is being undertaken to evaluate datasets. Furthermore, if time is included as a variable, as done in transient modeling, then the data processing requirements can increase nonlinearly. Therefore, a faster and more accurate interpolation method is required to effectively evaluate large datasets. Alternatively, the model domain can be decomposed into several regions, thereby allowing the necessary interpolations to be performed within a targeted sub-region (Mitášová and Mitáš, 1993; Trochu, 1993).

This study proposes a new method using irregular quadrilateral elements to address existing limitations to interpolation. This method is based on the interpolation technique used in finite element methods and can be applied to both rectangular grids and irregular meshes in finite-difference, finite-volume, and finite-element methods. This proposed method, hereafter called the quadrilateral irregular network (QIN), is similar to the triangular irregular network (TIN), which is a well-known conventional method. After introducing the QIN, conventional methods for the linear combination of values, such as TIN and the inverse distance weight method as well as methods based on radial basis functions, such as thin plate spline (TPS), completely regularized spline (CRS), dual kriging (DK), and Hardy (1971)’s multiquadric method (HMQ), are briefly summarized. The QIN and other methods are then estimated using the well-known Franke (1979) test functions, and applied to the results of numerical modeling for a comparison of their effectiveness in interpolating spatial data.

Interpolation Methods

Linear Combination of Values

Linear combination interpolation methods share the same general form as follows:

\(\begin{aligned}u^{*}(x, y)=\sum_{i=1}^{N} \lambda_{i} u\left(x_{i}, y_{i}\right)\end{aligned}\)       (1)

where u^* is the interpolated value of a variable at a target point (x, y), u is the known value at the i-th sampled point (x_i, y_i), λ_i is the weight assigned to the i-th sampled point, and N is the number of sampled points used for the interpolation.

Quadrilateral Irregular Network (QIN)

Finite element processes generally transform each quadrilateral element in a global coordinate system to a square isoparametric element that is defined in a ξ-η local coordinate system, as shown in Fig. 1. These local coordinates, ξ and η, are used to determine the sampling points of the Gaussian quadrature, which is a way to numerically approximate the value of an integration. A value (u) within a quadrilateral element is estimated from the nodal values of the isoparametric element as follows:

Fig. 1. Transformation of quadrilateral and triangular elements to isoparametric elements.

\(\begin{aligned}u=\sum_{i=1}^{N_{e}} N_{i}(\xi, \eta) u_{i}\end{aligned}\)       (2)

where N_e (= 4 for linear quadrilateral elements) is the number of nodes in an element, u_i is the known value at the i-th nodal point, and N_i(ξ, η) is the i-th shape function for a local element, which is defined in bilinear form (Becker et al., 1981; Heinrich and Pepper, 1990) as follows:

\(\begin{aligned}\begin{array}{l}N_{1}(\xi, \eta)=\frac{1}{4}(1-\xi)(1-\eta) \\ N_{2}(\xi, \eta)=\frac{1}{4}(1+\xi)(1-\eta) \\ N_{3}(\xi, \eta)=\frac{1}{4}(1+\xi)(1+\eta) \\ N_{4}(\xi, \eta)=\frac{1}{4}(1-\xi)(1+\eta)\end{array}\end{aligned}\)       (3)

Eqs. (1) and (2) are of the same form, except for the number of points used for the interpolation. The weights of the quadrilateral elements are determined by the ratio of the areas diagonally opposite the corresponding nodes to the total area, for example, N₁ = A₁/4 (see Fig. 1). Using Eq. (2), points in the local coordinate system can be mapped to corresponding points in the actual global coordinate system as follows:

\(\begin{aligned}x=\sum_{i=1}^{N_{e}} N_{i}(\xi, \eta) x_{i}\end{aligned}\)       (4)

\(\begin{aligned}y=\sum_{i=1}^{N_{e}} N_{i}(\xi, \eta) y_{i}\end{aligned}\)       (5)

Therefore, the values of x, y, and variables in the global coordinate system can be easily estimated at a given point (ξ, η) within an isoparametric element through Eqs. (2)~(5). However, for a point (x, y) within an irregular-shaped quadrilateral element in the global coordinate system, estimation of variables at a given point is difficult because the local coordinate values (ξ, η) corresponding to the point are not known. Therefore, inverse mapping is necessary to estimate the variables.

Substitution of Eq. (3) into Eqs. (4) and (5) yields:

x = a₁ + a₂ξ + a₃η + a₄ξη       (6)

y = b₁ + b₂ξ + b₃η + b₄ξη,       (7)

respectively, where the coefficients are a₁ = (x₁ + x₂ + x₃ + x₄)/4, a₂ = (-x₁ + x₂ + x₃ - x₄)/4, a₃ = (-x₁ - x₂ + x₃ + x₄)/4, a₄ = (x₁ - x₂ + x₃ - x₄)/4, b₁ = (y₁ + y₂ + y₃ + y₄)/4, b₂ = (-y₁ + y₂ + y₃ - y₄)/4, b₃ = (-y₁ - y₂ + y₃ + y₄)/4, and b₄ = (y₁ - y₂ + y₃ - y₄)/4. Coefficients ai and b_i are known values from four nodal coordinates of a quadrilateral element. Therefore, if the coordinate values (x, y) of an arbitrary point are given, then Eqs. (6) and (7) become non-linear equations of ξ and η. Although the numerical solution for ξ and η can be obtained through the iterative Newton-Raphson method, there is a possibility that a convergent solution cannot be obtained due to numerical error. Therefore, an analytic solution for ξ and η should be considered instead of an iterative numerical solution.

The quadratic equation for ξ such as Eq. (9) can be obtained by rearranging Eq. (6) for η and substituting the resultant Eq. (8) into Eq. (7).

\(\begin{aligned}\eta=\frac{x-a_{1}-a_{2} \xi}{a_{3}+a_{4} \xi}\end{aligned}\)       (8)

Aξ² + Bξ + C = 0       (9)

where the coefficients A = a₂b₄ - a₄b₂, B = a₂b₃ - a₃b₂ + (a₁ - x)b₄ - a₄(b₁ - y), and C = (a₁ - x)b₃ - a₃(b₁ - y). The solution for ξ is as follows:

\(\begin{aligned}\xi=\frac{1}{2 A}\left(-B \pm \sqrt{B^{2}-4 A C}\right)\end{aligned}\).       (10)

If a quadrilateral element is defined in a counterclockwise direction, then the positive square root should be applied to Eq. (10). Conversely, if a quadrilateral element is defined in a clockwise direction, the negative square root should be chosen as the solution for ξ. In addition, if coefficient A in Eq. (9) equals zero, then the quadratic equation is reduced to a linear equation, and the solution for ξ becomes:

ξ = -C/B       (11)

Coefficient A = 0 is equal to u₁v₂ - u₂v₁ = 0 when coefficient A is expressed by the components of vector u and v in Fig. 2a for a detailed geometrical configuration. This means that u × v = 0, that is, the two vectors are parallel. Once ξ is obtained using Eq. (10) or Eq. (11), η can be estimated from Eq. (8). However, η cannot be calculated when the denominator in Eq. (8) equals zero. In this case, another approach for estimating η should be considered.

Fig. 2. Vector representations for (a) the quadrilateral element, and (b) the center of a circumcircle.

The first condition that results in a zero denominator is when a₃ = a₄ = 0. By applying this condition to Eq. (6), a simplified equation such as x = a₁ + a₂ ξ can be obtained. Therefore, substituting ξ = (x - a₁)/a₂ into Eq. (7) and rearranging the equation yield the following:

\(\begin{aligned}\eta=-\frac{\left(a_{1}-x\right) b_{2}-a_{2}\left(b_{1}-y\right)}{-a_{2} b_{3}+\left(a_{1}-x\right) b_{4}}\end{aligned}\).       (12)

To explain the geometrical meaning of this condition, the coefficients a₃ and a₄ can be expressed as components of vectors p and q in Fig. 2a:

\(\begin{aligned}\begin{array}{l}a_{3}=\frac{1}{4}\left[\left(x_{4}-x_{1}\right)+\left(x_{3}-x_{2}\right)\right]=\frac{1}{4}\left(p_{1}+q_{1}\right) \\ a_{4}=\frac{1}{4}\left[-\left(x_{4}-x_{1}\right)+\left(x_{3}-x_{2}\right)\right]=\frac{1}{4}\left(-p_{1}+q_{1}\right) .\end{array}\end{aligned}\)       (13)

The condition a₃ = a₄ = 0 is only satisfied when p₁ and q₁ both equal zero, which means that p and q are parallel to both each other and the vertical axis.

The second condition is a₃ = ξ = 0 and a₄ ≠ 0. Substituting this condition into Eq. (6) shows that this situation occurs when x = a₁. Substituting ξ = 0 into Eq. (7) and rearranging the equation yield as follows:

\(\begin{aligned}\eta=-\frac{b_{1}-y}{b_{3}}\end{aligned}\).       (14)

This is the same as Eq. (12) when a₁ - x = 0. Geometrically, only a₃ may be zero when components p₁ and q₁ of vectors p and q in Fig. 2a are equal in magnitude and in opposite directions, as can be inferred from Eq. (13).

The third condition is when a₃ ≠ 0 and a₄ ≠ 0, but a₃ + a₄ ξ = 0. Applying ξ = -a₃/a₄ to Eq. (6) and Eq. (7), respectively, yield x = a₁ - (a₂a₃)/a₄ and the relation for η as follows:

\(\begin{aligned}\eta=-\frac{a_{3} b_{2}-a_{4}\left(b_{1}-y\right)}{a_{3} b_{4}-a_{4} b_{3}}\end{aligned}\).       (15)

It may be inefficient and time-consuming to compute ξ and η in all elements to find the element in which an arbitrary point is located. This is especially true when the number of elements is very large. Therefore, either rectangles or circles enclosing an element are utilized in this search. As shown in Fig. 3a, two opposite corners of a rectangle enclosing an element are compared with the coordinate (x, y) of an arbitrary point. It is then temporarily determined that an arbitrary point may be located within the element when the following inequality is satisfied:

Fig. 3. Methods for the discrimination of a temporarily internal point a quadrilateral element with (a) surrounding box, and (b) circumcircles.

(x _min, y _min) ≤ (x, y) ≤(x _max, y _max).       (16)

Alternatively, the center coordinates of the circumcircle are necessary when using circles that enclose an element. As shown in Fig. 2b, vectors u and v₁ and vectors w and v₂ are perpendicular to each other, respectively, based on the definition of the circumcircle. That is, u ∙ v₁ = 0 and w ∙ v₂ = 0, and the following system of linear equations is derived.

\(\begin{aligned}\left[\begin{array}{ll}\Delta x_{21} \Delta y_{21} \\ \Delta x_{31} \Delta y_{31}\end{array}\right]\left[\begin{array}{l}x_{c} \\ y_{c}\end{array}\right]=\left[\begin{array}{l}\bar{x}_{21} \Delta x_{21}+\bar{y}_{21} \Delta y_{21} \\ \bar{x}_{31} \Delta x_{31}+\bar{y}_{31} \Delta y_{31}\end{array}\right]\end{aligned}\)       (17)

where Δv_i1 = v_i - v₁, \(\begin{aligned}\bar{v}_{i 1}=\left(v_{\mathrm{i}}+v_{1}\right) / 2\end{aligned}\), and x_c and y_c are the center coordinates of the circumcircle. Two circumcircle radii can be obtained through Eq. (17) when a quadrilateral is divided into two triangles, as shown in Fig. 3b. If the distance between one of the center points of the circumcircles and an arbitrary point is less than or equal to the radii of the circumcircles, it is temporarily determined that an arbitrary point may be located within the element, and the exact values of ξ and η can be estimated analytically. If both ξ and η are in the range of [-1, 1], then an arbitrary point is definitely within the quadrilateral element, and the variables can be estimated by using Eq. (2).

Triangular Irregular Network (TIN)

It is simpler to find local coordinates within triangular elements than within quadrilateral elements. The local coordinates ξ and η within a triangular element can be estimated directly from global coordinates as follows (Becker et al., 1981; Heinrich and Pepper, 1990):

\(\begin{aligned}\xi=\frac{\frac{1}{2}\left|\begin{array}{ccc}x & y & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1\end{array}\right|}{\frac{1}{2}\left|\begin{array}{lll}x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1\end{array}\right|}=\frac{\text { Area of } \Delta P P_{2} P_{3}}{\text { Area of } \Delta P_{1} P_{2} P_{3}}\end{aligned}\)       (18)

\(\begin{aligned}\eta=\frac{\frac{1}{2}\left|\begin{array}{ccc}x_{1} & y_{1} & 1 \\ x & y & 1 \\ x_{3} & y_{3} & 1\end{array}\right|}{\frac{1}{2}\left|\begin{array}{lll}x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1\end{array}\right|}=\frac{\text { Area of } \Delta P_{1} P P_{3}}{\text { Area of } \Delta P_{1} P_{2} P_{3}}\end{aligned}\).       (19)

As with quadrilateral elements, the areal ratio is also used for triangular elements. An arbitrary point is within a triangular element when ξ and η satisfy the following inequality:

0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1, and ξ + η ≤ 1.       (20)

If so, some values at an arbitrary point are interpolated from the nodal values of an isoparametric triangular element by using Eq. (2). The shape functions for triangular elements are as follows:

N₁ (ξ, η) = ξ

N₂ (ξ, η) = η

N₃ (ξ, η) = 1 - ξ - η       (21)

which are the same as the quadrilateral elements in that the ratio of the areas diagonally opposite from the corresponding nodes to the total area is used as a weight (see Fig. 1). As in the quadrilateral, it can be preliminarily evaluated whether or not an arbitrary point is located within the element by using either a rectangle or circle that encloses a triangular element.

Inverse Distance Weight Method

In distance weighting interpolation methods, the weight of the i-th sampled point is defined as a function of the distance from the target point to the sampled points, given in general as follows:

\(\begin{aligned}\lambda_{i}=\frac{f\left(d_{i}\right)}{\sum_{j=1}^{N} f\left(d_{j}\right)}\end{aligned}\)       (22)

where function f(d_i) is defined as d_i^p , and d_i is the distance from the target point to the i-th sampled point. When the power parameter p of -2 or -1 is selected, this is called the inverse square distance (ISD) or reciprocal distance weight (RDW) method, respectively. Therefore, the weights diminish as the distance increases, and the values sampled nearby have more influence on the estimation. When the target point coincides with one of the sampled points, the value of the target point should take the value of the sampled point because the weight becomes infinite. For wide areas or for large data sets, the searching radius or threshold distance (d_max) may be set (Caruso and Quarta, 1998), and the function of distance is defined as follows:

\(\begin{aligned}f\left(d_{i}\right)=\left\{\begin{array}{ll}d_{i}^{p} & \text { for } d_{i} \leq d_{\max } \\ 0 & \text { for } d_{i}>d_{\max }\end{array}\right.\end{aligned}\)       (23)

This situation is called localized ISD or localized RDW (Wagner et al., 2012).

Radial Basis Functiom

Interpolation functions based on the radial basis function have the following general form (Talmi and Gilat, 1977; Mitášová and Mitáš, 1993):

\(\begin{aligned}u^{*}(x, y)=T(x, y)+\sum_{j=1}^{N} \beta_{j} R\left(d_{j}\right)\end{aligned}\)       (24)

\(\begin{aligned}T(x, y)=\sum_{k=1}^{M} \alpha_{k} f_{k}(x, y)\end{aligned}\)       (25)

where R(∙) is a radial basis function, T(∙) is called the trend function, f_k(∙) is a set of linearly independent monomial functions, and α_k and β_j are coefficients. Once the radial basis function is determined, the coefficients α_k and β_j can be obtained by solving the following system of linear equations.

u*(x_i, y_i) = u(x_i, y_i) for i = 1, 2, ⋯, N

\(\begin{aligned}\sum_{j=1}^{N} \beta_{j} f_{k}\left(x_{j}, y_{j}\right)=0 \text { for } k=1,2, \cdots, M\end{aligned}\)       (26)

In a geostatistical framework, Eq. (25) is called a drift, the second term on the right hand side of Eq. (24) is called a fluctuation, and the function R(∙) within the second term on the right hand side of Eq. (24) is replaced with the function K(∙), which is called the generalized covariance. The final M equations in Eq. (26) are called the no-bias condition, which shows that the trend function or drift represents the mean of sampled values. Therefore, the fluctuation term is adjusted so that the interpolation coincides with the data points (Trochu, 1993).

Thin Plate Spline (TPS)

Duchon (1976) proposed the following trend and radial basis functions for the TPS interpolation method.

T(x, y) = α₁ + α₂x + α₃y       (27)

R(d_j) = d²_j ln d_j       (28)

Therefore, coefficients α_k (1 ≤ k ≤ 3) and β_j (1 ≤ j ≤ N) in Eq. (24) are characterized by the following linear system.

u*(x_i, y_i) = u(x_i, y_i) for i = 1, 2, ⋯, N

\(\begin{aligned}\sum_{j=1}^{N} \beta_{j}=\sum_{j=1}^{N} \beta_{j} x_{j}=\sum_{j=1}^{N} \beta_{j} y_{j}=0\end{aligned}\)       (29)

Completely Regularized Spline (CRS)

Mitášová and Mitáš (1993) derived a new interpolation method called the completely regularized spline as an alternative method to regularized TPS with tension which was developed to compensate for the shortcomings of TPS by Mitáš and Mitášová (1988). The trend and radial basis function of CRS are as follows:

T(x, y) = α₁       (30)

\(\begin{aligned}R\left(d_{j}\right)=-\left\{\ln \left[\left(\frac{\phi d_{j}}{2}\right)^{2}\right]+E_{1}\left[\left(\frac{\phi d_{j}}{2}\right)^{2}\right]+C_{E}\right\}\end{aligned}\)       (31)

where E₁(∙) is the exponential integral function found in Abramowitz and Stegun (1972), C_E is the Euler constant of 0.5772156649, and φ is a generalized tension parameter that should be determined empirically. Consequently, coefficients α_k (k = 1) and β_j (1 ≤ j ≤ N) in Eq. (24) can be obtained from the N + 1 linear equation.

Dual Kriging (DK)

Trochu (1993) proposed the following generalized covariance for the 2-D case derived from TPS in order to incorporate the distance of influence in his dual kriging model.

\(\begin{aligned}R\left(d_{j}\right)=\left\{\begin{array}{cc}1+\frac{1}{c_{1}}\left(\frac{d_{j} c_{0}}{d_{\max }}\right)^{2} \ln \left(\frac{d_{j} c_{0}}{d_{\max }}\right) & \text { for } 0 \leq d_{j} \leq d_{\max } \\ 0 & \text { for } d_{j}>d_{\max }\end{array}\right.\end{aligned}\)       (32)

where c₁ is a constant of 0.18393972, and c₀ is a constant of 0.60653066, which is the value when the function d² lnd attains its minimum; -c₁ = c²₀lnc₀. The function R(∙) in Eq. (32) continuously decreases from 1 when d_j = 0 to 0 when d_j > d_max. As in TPS, the values of coefficients α_k and β_j in Eq. (24) can be obtained by solving N + 3 linear equations.

Hardy’s Multiquadric (HMQ) Method

Hardy (1971) proposed the following relatively simple function for the interpolation of irregular surfaces.

R(d_j) = (d²_j + Δ²)^1/2       (33)

Eq. (33) represents the multiquadric surface, which is a single cone or hyperboloid dependent upon whether Δ² equals zero or is a positive constant, respectively. In HMQ, the trend function of Eq. (25) is zero. Therefore, β_j (1 ≤ j ≤ N) in Eq. (24) is obtained by solving linear equations of order N. For reference, the method of applying a power of -1/2 instead of 1/2 to the right side of Eq. (33) is called the reciprocal multiquadric (RMQ) method, and successful applications of HMQ in various disciplines can be found in Hardy (1990).

Comparison of Interpolarion Methods

Test Functions and Data Sets

The comparison of interpolation methods described in the preceding chapter is based on the six test functions used by Franke (1979) in his critical comparison of interpolants for scattered data, reported as:

\(\begin{aligned} f_{1}(x, y)= & 0.75 \exp \left[-\frac{(9 x-2)^{2}+(9 y-2)^{2}}{4}\right]+0.75 \exp \left[-\frac{(9 x+1)^{2}}{49}-\frac{9 y+1}{10}\right] \\ & +0.5 \exp \left[-\frac{(9 x-7)^{2}+(9 y-3)^{2}}{4}\right]-0.2 \exp \left[-(9 x-4)^{2}-(9 y-7)^{2}\right] \\ f_{2}(x, y)= & \frac{\tanh (9 y-9 x)+1}{9} \\ f_{3}(x, y)= & \frac{1.25+\cos (5.4 y)}{6\left[1+(3 x-1)^{2}\right]} \\ f_{4}(x, y)= & \frac{1}{3} \exp \left\{-\frac{81}{16}\left[\left(x-\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}\right]\right\} \\ f_{5}(x, y)= & \frac{1}{3} \exp \left\{-\frac{81}{4}\left[\left(x-\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}\right]\right\} \\ f_{6}(x, y)= & \frac{1}{9} \sqrt{64-81\left[\left(x-\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}\right]}-\frac{1}{2}\end{aligned}\)       (34)

Four of the six test functions (test functions 2 and 4-6) were originally provided by McLain (1974), then slightly modified and scaled down from [1, 10] to [0, 1] by Franke (1979). A total of three grid networks were constructed as shown in Fig. 4. The first case is an 11 × 11 square grid, which has an origin of (0, 0) with both a horizontal (Δx) and vertical (Δy) interval of 0.1, whereby 121 total points were used to generate the test function values in Eq. (34). The known values in the sparse square grid were then interpolated into points on a fine grid that was rotated 60° counterclockwise from the horizontal axis and consisted of 20 rows and 31 columns, with its origin at (0.35, 0.02), Δx = 0.03 and Δy = 0.02 (Fig. 4a). The second case is a 13 × 13 square grid that was rotated 15° counterclockwise from the horizontal, with its origin at (0.1, -0.25), and Δx = Δy = 0.1. The 169 generated values were then interpolated to points on a 19 × 19 fine square grid, with its origin at (0.05, 0.0), Δx = Δy = 0.05, and no rotation (Fig. 4b). The third case is a mesh grid consisting of 120 randomly generated points. The triangular elements were constructed from the distribution of these points over the entire area, with some triangles then combined into quadrilateral elements for the QIN application. Therefore, the third case was divided into a fully triangular element mesh (Fig. 4c) and a mixed mesh with triangles and quadrilaterals (Fig. 4d). The nodal values were interpolated to points on an 11 × 11 square grid, with its origin at (0.0, 0.0), Δx = Δy = 0.1, and no rotation. Therefore, both pure TIN and a combination of QIN and TIN (hereafter QIN&TIN) were applied to the third case, whereas only QIN was used for the first and second cases.

Fig. 4. Target points (dots) and neworks (solid lines) for (a) case 1 and (b) case 2, and mesh networks (solid line) for cade 3 composed of (c) triangular elements and (d) mixed elements.

Results of Test and Discussion

The sensitivity analysis results, presented as a function of a given parameter for each interpolation method and test case, are shown in Figs. 5~7. The lower horizontal axis is the searching radius (d_max) for DK, RDW, and ISD and the parameter value (Δ²) for HMQ, and the upper horizontal axis is a generalized tension parameter value (φ) for CRS. The vertical axis represents the corresponding RMSEs on a logarithmic scale. The resultant QIN (cases 1 and 2), and QIN&TIN and TIN (case 3) RMSEs, which appear as constant values, are provided for comparison with the other methods.

Fig. 5. RMS errors versus parmeter values in the first case for (a) test function 2, (c) test function 3, (d) test function 4, (e) test function 5, and (f) test function 6.

Fig. 6. RMS errors versus parameter values in the second case for (a) test function 1, (b) test function 2, (c) test function 3, (d) test function 4, (e) test function 5, and (f) test function 6.

Fig. 7. RMS errors versus parameter values in the third case for (a) test function 1, (b) test function 2, (c) test function 3, (d) test function 4, (e) test function 5, and (f) test function 6.

At first, the methods using a searching radius are discussed. With RDW or ISD, the most accurate results were obtained with d_max of about 0.1 for the first and second cases, and about 0.2 for the third. As d_max increases from these values, the accuracy of the interpolants decreases. Eventually, the effect of the searching radius disappears from around 0.6. It is difficult to find the optimal d_max when RDW or ISD is applied to complex natural phenomena. Therefore, in many practical applications of RDW and ISD, weighting coefficients are generally calculated using all data points without considering the searching radius. Although RDW and ISD are well-known and popular methods due to their practical simplicity, they seem to be the least sophisticated among all of the methods tested in this study when the searching radius is not considered. Trochu (1993) introduced the searching radius for geostatistical kriging in order to account for the distance of influence and, as a result, reduced the computational effort due to the derivation of banded structures in the kriging matrix. He mentioned that d_max should be selected carefully because it controls the smoothness of the interpolation, the loss in accuracy may be unacceptable when d_max is small, and there is no difference between the results with and without the distance of influence when d_max greater than a certain threshold value is selected. However, he did not indicate a specific value or range as the optimal d_max for DK. As shown in Figs. 5~7, the accuracy of DK is improved as d_max increases, reaching an asymptotic constant equal to the RMSE of TPS. Therefore, to avoid the ambiguity of d_max as a parameter and the effort of trying to find the proper value for d_max, it would be better to use TPS instead of DK in practice.

The influence of Δ² on the HMQ results is now discussed. Carlson and Foley (1991) stated that larger Δ² values are necessary in HMQ when using the data from spherical or quadratic-shaped functions such as test functions 4 and 6, whereas very small Δ² value are required for rapidly varying or highly oscillatory sampled values. As shown in Figs. 5~7, HMQ yields better results than the other methods for test function 4, even with relatively large Δ² values, although this is not necessarily true for test function 6. Several researchers have attempted to find an optimal Δ² value. The formula of Hardy (1971) yields Δ² values of 0.011, 0.011, and 0.015, and the Franke (1979) formula yields Δ² values of 0.026, 0.027, and 0.034 for cases 1, 2, and 3, respectively. However, selecting these values as the optimal Δ² values is not appropriate when considering the results in Figs. 5~7. Carlson and Foley (1991) argued that the optimal Δ² value should be estimated with the values of the sampled data as well as the number and location of data points, and proposed their own method for optimal Δ² in which the data points and values were scaled to the unit cube and interpolated by applying their formula for an effective Δ² value, and then the resulting interpolants were scaled back to the original domain. Kansa (1990) and Kansa and Carlson (1992) suggested another method in which various Δ² values for all data points were assessed rather than a constant. However, some ambiguity still remains with their method since the user should provide a minimum and maximum constraint and Δ² is not dependent on the sampled values but only on the data points. Nonetheless, it is apparent that when a value in the range of about 0.02 to 0.3 is selected for Δ² , better results than the others can be obtained for all functions, with the exception of test function 4, based on the results of Figs. 5~7.

Mitášová and Mitáš (1993) reported that a generalized tension parameter φ in CRS adjusts the ratio between the weights of the first- and higher-order derivatives in the smooth seminorm and thus controls the behavior of the resulting surface from a membrane to a thin plate. They also remarked that the value of φ should be determined empirically and can be found within a few trials. However, as shown in Figs. 5~7, the CRS accuracy varies the most depending on the selected φ value, thereby highlighting the potential difficulty in selecting the appropriate φ value. These test results indicate that a φ value in the range from 10 to 20 is recommended to obtain more stable and accurate CRS results.

As stated above, the results of interpolations may be greatly influenced by the selected parameter value. In addition, the suggested optimal Δ² and φ ranges in HMQ and CRS, respectively, are limited to the test functions that correspond to the true value in this study, and can also be changed according to the number and distribution of sampled points. Through cross-validation methods, it is possible to find the optimal parameter value that minimizes the error of interpolation. However, this approach may be time-consuming and computationally intensive. Furthermore, the accuracy of the interpolated results for a natural phenomenon, especially time-series data which are complex and difficult to predict, cannot be guaranteed. Therefore, it would be preferable to use interpolation methods with no parameters and a certain degree of guaranteed accuracy, such as TPS and QIN.

Application of Interpolation Methods

Fig. 8 shows an example of the numerical results of flow characteristics in a U-shaped channel. The width of the channel is 1.7 m, and the radius of the center line in a curved section is 4.25 m. The slope of the channel bed is 1.76471 × 10^-3. Meshes are roughly composed of quadrilateral elements for which Δx and Δy in straight reaches are 0.34 m and 0.17 m, respectively, and Δr and Δθ in a curved reach are 0.17 m and 4.5°, respectively (here, r and θ represent the cylindrical coordinate system). The total number of nodes and elements is 1551 and 1400, respectively. When numerical modeling was performed with specific boundary conditions for the upstream and downstream ends, a slip condition in the lateral side walls, and a fixed bed condition, longitudinal mean flow velocity and water surface elevation were obtained as shown in Fig. 8. Assuming that the observational data at 61 and 15 points were regularly obtained along the A-A' and B-B' routes in Fig. 8a, respectively, the observational points did not coincide with the nodal points. In this case, it is necessary to apply interpolation methods for comparison of the observed values with the numerical results.

Fig. 8. Distribution of (a) longitudinal mean flow velocity (m/s) and (b) water surface elevation (m) using a quadrilateral mesh for a U-shaped channel.

The searching radius of the RDW and ISD was not considered, and a generalized tension parameter φ of 15 in CRS and Δ² of 0.1 in HMQ was chosen for this application. Fig. 9 shows the interpolated results for channel bed, water surface, and mean velocity. Since the results for CRS, TPS, and HMQ are almost identical, they are represented by CRS dotted lines in Fig. 9a. The mean velocity with the QIN along A-A' shows slight differences from those obtained with CRS, TPS, and HMQ, while the results for channel bed and water surface are almost the same for QIN, CRS, TPS, and HMQ, as shown in Fig. 9a. The CRS, TPS, and HMQ show smoothed curves based on radial basis functions, while QIN shows a slight oscillation in the mean velocity. The RDW and ISD show more oscillatory results than the QIN, as well as underestimated values for the mean velocity compared with the others. On the other hand, the water surface and channel bed results along B-B' using HMQ in Fig. 9b show some oscillations unlike those of QIN, CRS, and TPS, while the interpolated mean velocity is the same for QIN, HMQ, CRS, and TPS. Since the results from the QIN are almost identical to those from CRS and TPS, the results from CRS and TPS are overlapped with the representative solid line for the QIN in Fig. 9b.

Fig. 9. Interpolated results of bed elevation, water surface elevation, and longitudinal mean flow velocity along (a) section A-A' and (b) section B-B' shown in Fig. 8a.

Comparing the processing speed and results, it is clear that the ISD and RDW are superior in terms of calculation time, but produce inferior results when the searching radius is not considered. The CRS and HMQ with properly selected parameters and TPS are fairly reliable, but they are computationally intensive because the entire matrix must be resolved for each variable. The QIN is as fast as the ISD and RDW, and has reliable accuracy to some extent. Therefore, the QIN is proposed as an alternative to traditional interpolation methods for large data fields that are composed of quadrilateral elements.

Conclusion

Analytical solutions for the inverse mapping of local coordinates in an isoparametric square element corresponding to the global coordinates of arbitrarily given points were proposed in this study. In addition, a method for preliminarily determining whether an arbitrary point is within an element with a rectangle or circumcircles enclosing an element was proposed for fast calculation before the estimation of exact local coordinates. Once it was determined that an arbitrary point was located within the quadrilateral element, the values at a target point were interpolated using shape functions of the finite element method when the evaluated local coordinates were between -1 and +1. The details of this proposed algorithm, the quadrilateral irregular network (QIN), were presented, and interpolation schemes of conventional linear combination methods (TIN, RDW, and ISD) and methods based on radial basis functions (TPS, CRS, DK, and HMQ) were briefly summarized. All of these interpolation methods were applied to six test functions and three mesh types to evaluate their accuracy. The results showed that parameters such as searching radius in RDW, ISD, and DK, as well as Δ² in HMQ, and tension in CRS have a significant influence on the interpolation results, and that inverse distance methods without consideration of a searching radius produce inferior results compared with the other methods. It is difficult to find the optimal parameter value in practical applications because it greatly depends on sampled values as well as the number and location of data points. Therefore, we cautiously recommend the use of methods such as TPS or QIN in order to avoid ambiguity and minimize the amount of effort required for determining non-physical parameter values, although not all of the known interpolation methods were compared in this study. The accuracy of QIN is guaranteed compared with the other methods and there is no need to determine the optimal parameter values. If an initial mesh is maintained without changes, coordinates ξ and η within an isoparametric element for a target point can be saved and reused in subsequent interpolations. Consequently, the QIN is computationally more efficient than conventional methods that require the interpolation values to be reevaluated using all of the nodal values, even when only one of the nodal values changes.