Adaptive sampling in higher dimensions for point-wise experimental measurement techniques

When performing point-wise measurements within a pre-defined domain, the experimentalist is faced with the problem of defining the spatial locations where to collect data based on an a priori unknown underlying signal. While structured sampling grids are most common, these are rarely optimal from a time-efficiency perspective. In this work an adaptive process is presented for point-wise measurement techniques to guide the spatial distribution of sampling locations in two-dimensions. Thin plate splines of varying degree of smoothing are utilised to iteratively obtain surrogate models of the spatial distribution of the quantities of interest. Adaptive sampling criteria reflecting typical human decisions are combined into a unique objective function, allowing the derivation of suitable measurement positions. The total number of sampling locations thus no longer requires to be pre-defined, as the sampling process is terminated automatically once the addition of measurements is deemed not to provide additional information. The overall implementation is elaborated and the methodology is assessed on the basis of numerical simulations and an experimental test case. Results support the new method’s advocated superiority compared to traditional full-factorial sampling in terms of accuracy and reliability.


Introduction
Whereas full-field image-based measurement techniques such as particle image velocimetry (PIV) are capable of providing velocity data within planar or even volumetric sections of the scrutinised flow (Kähler et al 2016), highest spatial resolutions (excluding advanced PIV alternatives e.g. micro-PIV (Meinhart et al 1999)) continue to be attainable with more traditional point-wise methodologies e.g. hot-wire anemometry (HWA) or laser Doppler anemometry (LDA) (Lavoie et al 2007). Contrary to PIV where data rates are defined by the laser's pulse repetition and camera frame rate (and therefore associated costs), metrologies such as HWA and LDA also offer cheaper alternatives with high frequency responses. This makes such techniques extremely suitable for spectral analyses in turbulent flow research. Although the limited measurement volume of the point-wise techniques enables the detection of smaller flow scales, because measurement volumes are typically hundreds of orders of magnitude smaller compared to the flow region of interest, some form of traversing is required when performing a flow survey. Typical HWA wire diameters for example are in the order of 5 μm compared to flow-related length scales in the order of centimeters (Gao et al 2015). Unless flow-related time scales are much larger than the Measurement Science and Technology Adaptive sampling in higher dimensions for point-wise experimental measurement techniques overall measurement time, full-field data analyses with pointwise measurement techniques are therefore often limited to temporal averages, i.e. time-average flow fields.
Point-wise probes are commonly mounted on traverse mechanisms and are therefore able to be, in principle, positioned continuously in space. Ignoring limitations imposed by the finite probe dimensions, the level of flow detail captured is thus directly defined by the spacing of spatial sampling locations considered through the well-known Nyquist sampling criterion (Petersen and Middleton 1962). This criterion dictates the smallest measurable length scale to be at least twice the minimum sample spacing. The experimentalist is thus forced to compromise between data content and time efficiency: the more samples of the flow, the longer the measurement time, the higher the running costs of the experimental facility yet the more flow detail can be obtained. In turn, when faced with the task of characterising an unknown flow phenom enon, every experimentalist will be ultimately faced with the same conundrum: where to measure and how many samples to take.
Of the various sampling strategies, space-filling is the most common and includes, among others, the distribution of equidistant samples across the measurement domain (fullfactorial sampling: FF), Latin Hypercube (LH) sampling, Poisson disc sampling, etc. The number of samples adopted in such methods will always vary as a power of the ratio between domain extent L and user-defined minimum sample spacing h. For full factorial and Poisson disk sampling in d-dimensional space, the number of samples will for example be in the order of O((L/h) d ), whereas, because of the permutations, LH typically results in O(L/h) samples. As the user refines the sampling grid, the number of measurement points clearly becomes rapidly impractical. Assuming the data collection spans ten seconds at each measurement location (common in LDA), characterisation of a three-dimensional velocity field with only ten samples along each coordinate axis will require a total run-time of approximately three hours for 10 3 samples while yielding a still under-resolved flow field. It is for this reason that more intelligent sampling approaches are needed, capable of focusing, and thereby minimising, the data accrual to flow regions in which a higher number of measurements are beneficial in the reconstruction and subsequently understanding, of the underlying a priori unknown flow field.
All space filling methods blindly adopt the identical sampling strategy irrelevant of the scrutinised process (LH does involve the probability density function of each variable in theory, yet for spatial sampling each coordinate has a uniform probability, eliminating any probability dependency). Instead, adaptive processes whereby the sampling strategy is adjusted in situ on the basis of collected observations, have the potential to yield results more efficiently if the stopping criteria are quality-defined, or with higher accuracy, when limiting the permitted number of accrued data points. It is therefore not surprising that adaptive sampling has been widely applied in e.g. surface reconstruction (Chen and Peng 2017), survey studies (Yu et al 2012), computer vision (Bolin and Meyer 1998), computer graphics (Chen et al 2011), computer aided design (Li 1995), optimisation (Mackman et al 2013) and in the prediction of suitable measurement locations for PIV (Theunissen et al 2007) or one-dimensional velocity surveys (Theunissen et al 2015).
The typical adaptive process outline is depicted in figure 1 and starts with collecting data at initial sampling locations, defining the domain outline in which an evaluation grid is specified. This grid defines the locations at which the surrogate model is evaluated and must subsequently be sufficiently small to capture the model's behaviour. The surrogate model itself is constructed from the available data and needs to be updated each time new data is accrued. Various routes exist to construct surrogate or meta-models such as for example neural networks (Holeňa et al 2010), support vector machines (Ciccazzo et al 2014), splines (Grimstad et al 2015), polynomial regression (Gano et al 2006), Kriging or Gaussian process regression (Denimal et al 2016) and radial basis functions (RBF) (Durantin et al 2017).
Adaptivity is finally based on some form of criterion, the so-called objective, enabling potential samples to be ranked according to their quantified suitability. In numerical grid refinement for example, this criterion can be based on a second derivative error norm of the surrogate model (Li 2010), local field curvature (Mackman and Allen 2010) or a combination of various functionals. When dealing with sampling optimisation, such cost functions are generally based on estimates of the spatial distributions of the scrutinized parameters (Forrester et al 2008).
In the work of Theunissen et al (2015), which forms the basis of the current study, a procedure was presented whereby sequential sample locations were selected adaptively thus eliminating user-dependency and/or a priori knowledge of the investigated flow. Weighted criteria involving local measurement error estimates, space filling heuristics, local surrogate improvement and curvature were combined into a unique objective function such that its peak location indicated the most suitable sampling location, replicating the exper imentalist's decision processes. A C 2 -continuous Wendland function, extended with third-order polynomial interpolation, was adopted because of its ability to provide second order derivative estimates with relative computational ease. The autonomous methodology was shown to be more efficient and able to capture flow features of interest more accurately, although restricted to one-dimensional, single component velocity profile extraction. In the current work such concepts are extended to higher dimensions with a focus on implementation.
The aim of the current work is to devise an autonomous routine capable of guiding the spatial sampling process in a typical experiment in which point-wise measurement techniques are utilized to survey a prescribed domain in terms of scalar quantities such as e.g. pressure, velocity, temperature, etc. At each sampling location the statistical average of multiple readings serves as the quantity of interest in the characterisation of the underlying phenomenon. Given the exponential increase of the number of sampling positions with dimensionality of the problem (and consequently the experiment run time), irrespective of adaptation, considered applications are limited to two-dimensions. The positioning of the data extraction locations within the 2D domain is performed through an adaptive process similar to Theunissen et al (2015) but extended into higher dimensions. To this extent the current study starts with a review of the objective function followed by a discussion on the choice of RBF kernel for interpolation. Additional computational procedures are presented to ensure viability and reduce computational effort in each of the different procedural steps. Preceding the conclusions, a numerical assessment as well as an experimental application with varying degree of data uncertainty corroborate the advocated advantages.
It is important to note that the importance of efficient collection of meaningful data spreads beyond the boundaries of experimental fluid mechanics and has a broader application referred to as the design of experiments. As such, the proposed methodology and associated implementation is expected to benefit a wider range of research.

RBF interpolation
In its simplest form the surrogate model in d-dimensions consists of an interpolation of available data f (x t ) with x = (x 1 , . . . , x d ) onto the evaluation grid. In view of required adaptivity to local curvature, this interpolant must be at least second order continuous (C 2 ). Because of their simplicity and ease of implementation, polynomial regression and RBF interpolation are the most common methods (Fasshauer 2007). Moreover, RBF functions have been shown to substantially reduce effects of measurement noise when interpolating (Casa and Krueger 2013). Given N (i) data points (t = 1, . . . , N (i) ), the value of the surrogate model s (i) (x s ) in the ith iteration at a random location x s in d-dimensional space is formulated through an augmented RBF model with polynomial summands as where ψ(r) accounts for the contribution of the tth basis function in function of the radial distance r = x t − x s between the new location and the tth data point. It is evident that with every iteration the number of samples N (i) can change and the metamodel s (i) (x) needs to be re-evaluated. Without loss of clarity, the superscript (i) will be omitted for brevity hereafter. The polynomial part q(x s ) is added to guarantee polynomial accuracy up to degree m and is given by where a possible choice of polynomial basis functions is given by Scaling coefficients α t , β t ∈ R are obtained in a straightforward manner by solving equation (1) for the known data points, i.e. the interpolation conditions s(x t ) = f (x t ) with the M additional constraints N t=1 α t p j (x t ) = 0, j = 1, . . . , M. This is equivalent to solving the matrix system The additional term λI is referred to as ridge-regression and introduces a degree of smoothing, which can be controlled through λ.
Having calculated the coefficients for each data point, these are then substituted into equation (1) to solve for the interpolation at any random location. Especially the Wendland kernels have received great attention as choices for ψ(r) because these can be easily extended to higher dimensions and any degree of smoothness (Wendland 2004). More importantly, these kernels introduce sparsity in the matrix systems to solve (see equation (2)) because of their compactness (Schaback and Wendland 1999), thus reducing computational effort. However, RBF functions require a well-defined support radius since enlarging the RBF kernel generally leads to enhanced interpolation accuracy at the expense of denser and ill-conditioned matrix systems (Boyd and Gildersleeve 2011). Too large values of the support cause the RBF interpolant to exhibit spurious overshoots reminiscent of the Gibbs phenomenon for polynomial interpolants. To this extent (Schaback and Wendland 1999) suggested the support radius to be chosen such that each RBF support contains approximately the same number of samples, which in turn depends on the sample density (i.e. the fill distance defined as the largest of the nearest neighbor sample distances). Especially given the expected non-uniformity in sample density, a unique choice of support radius may prove troublesome and much continuing research is devoted to the proper choice of a potentially adaptive scale parameter (Rippa 1999, Scheuerer 2011, Uddin 2017. Thin Plate Splines (TPS) on the other hand do not require shape parameters (Fasshauer 2007) and are more intuitive as they provide an interpolation minimising the overall 'bending' energy. They are inhibitive though as they are intrinsically not compact (Li et al 2007) and thus involve solving large matrix systems (Wood 2003). However, experimental applications targeted in the current work will concern a number of sampling locations ordinarily in the order of 10 3 -10 4 , rendering the matrix system (depending on the choice of evaluation grid density) of only low to moderate complexity. Considerable work has also been done in optimising the ridge regression parameter λ (see equation (2)) to accommodate splines for noisy data (Wahba 1990). This will prove to be advantageous as it allows the sampling process to consider recursively refined spatial scales of the underlying function through iterative reduction of λ (see section 2.3). Moreover, the TPS kernel has a well-defined integral, contrary to the Wendland functions, which will be shown to reduce the computational complexity of spatial filtering based on moving averages (see section 3.3). For this reason, in the following TPS radial basis functions of second (n = 2) and third order (n = 3) (see table 1) will be adopted.
In terms of polynomial order m, the inclusion of high order polynomials has been reported to worsen the conditioning problem as well as increasing the cost of the matrix solution (Costin and Allen 2013). It is therefore advisable to retain the polynomial order to a minimum. Following (Duchon, 1977) the order of the polynomial will be at most m = n − 1 (m = 1 for TPS2 and m = 2 for TPS4).

Preconditioning
The definition of TPS basis functions in two dimensions (d = 2) and associated derivatives are tabulated in table 1. The most common TPS is that of order two (TPS2), despite its second derivative exhibiting an unavoidable singularity at the origin. In view of the requirement to calculate the surrogate model's curvature (see section 3.1), the order three variant (TPS4) is more suitable. On uniform grids, the TPS4 has a condition number which increases as O(N 7 ) compared to O(N 4 ) for TPS2 (Boyd and Gildersleeve 2011). This translates in ill-conditioned, and subsequently, unstable solutions for the coefficient vectors α and β in equation (2). To overcome this problematic, interpolation matrices are preconditioned and TPS2 kernels are adopted to construct the surrogate models whereas the TPS4 kernels are applied in the calculation of curvature heuristics. The latter will be elaborated in section 3.1.
The preconditioned equivalent of matrix system Ax = y is given as CAx = Cy such that CA is better conditioned. For two dimensional problems the preconditioning strategy proposed by Dyn et al (1986) is specifically tailored to TPS2 and has been extended here for TPS4. Contrary to equation (2), the system considered in Dyn did not involve ridge regression. However, the addition of λ along the diagonal elements of Ψ does not invalidate the presumptions that the iterated Laplacian of the radial basis function (including the TPS kernels) tends towards zero as the radial distance increases and approaches infinity in the limit of zero radial distance. Preconditioning also intends to increase the diagonal elements relative to the off-diagonal ones in Ψ. Ridge regression enforces such differences and as such, the described preconditioning remains valid.
The preconditioning matrix C in Dyn is obtained through the The summation can be represented on an analogue (Delaunay) triangulated domain containing N ∆ triangles on the basis of the N datapoints; Functions u and v can be approximated on the jth triangle, which has an area T j , as nth order polynomials ũ j = n t=0 a t,j x t n t=0 b t,j y t = p+q n ξ pq,j x p y q = X pq ξ j and similarly ṽ j = X pq η j . Vectors ξ j and η j contain the N ξ = 1 2 (n + 2)(n + 1) polynomial coefficients for ũ and ṽ respectively, which can be solved for in a least-squares manner ξ j = X T pq X pq −1 X T pq u j = S j u j and alike η j = S j u j . Vectors u j and v j are N v,j × 1 in size and contain the u and v-values respectively of the N v vertices neighbouring the jth triangle, i.e. N v,j . These matrix systems can be rewritten such that all N vertices are included; S j u j = Bu and S j v j = Bv. Matrix B is then a N ξ × N sparse matrix with copied entries of the N ξ × N v,j matrix S j and vectors u and v contain the u and v values of all the N vertices. Using the polynomial representations, each of the derivative terms in equation (3) can now be expressed in terms of ξ j and η j and consequently with reference to discrete values u and v. Matrix C is constructed by summing up the different terms pertaining the contributing polynomial coefficients for all triangles. The process is clarified for the case of n = 3 in appendix A.1. It was shown in Dyn et al (1986) that this type of preconditioner reduces the condition number by several orders of magnitude. Solutions for the TPS coefficients α are retrieved through an iterative conjugate gradient scheme, which is presented in appendix A.2 for completeness.

Scale-adaptive TPS-based modelling
Measurement data will unavoidably contain some level of uncertainty due to noise originating from temporal fluctuations in the data, limited spatial resolution of the measurement Table 1. 2D TPS basis functions of order n, ψ(r) = r 2n−2 log(r), with associated derivatives.
system, electronic noise, etc. In view of adjusting the sampling process to signal curvature, such regions of jitter will unnecessarily attract samples and continue to do so, as with each new sample the noise component will alter the surrogate model and amplify higher order spatial derivatives. A degree of filtering is therefore required. Data can be low-pass filtered prior to the interpolation process, although this requires an appropriate choice of smoothing kernel and extent (Carr et al 2003). Both are dependent on the data landscape a priori unknown to the user. Alternatively, ridge regression presented in equation (2) renders the RBF interpolation an approximant. With increasing λ, the surrogate model equation (1) becomes a smoother representation of the underlying data, effectively modulating smaller scale oscillations. To this extent the choice of λ must be in relation to the noise level. However, over-smoothing must be avoided as oscillations (in space) may genuinely originate from smaller scale phenomena such as e.g. turbulence. For this reason the adaptive sampling routine is nested in a primary process in which the degree of smoothing, λ (k) , is progressively reduced over three iterations (figure 2).
In the first iteration, values of λ (1) are obtained on the basis of ensemble average data scatter σ t = 1 N N t=1 σ t (Forrester et al 2008) and are updated with each new adaptively placed data sample. Multiple readings are taken at every sampling location enabling the estimation of local signal scatter as the signal standard deviation σ t (see also section 3.2). Once the adaptive sampling process has converged, the second iteration is initiated in which a smoothed TPS2 spline is estimated recursively using a generalised cross-validation (GCV) scheme as suggested and implemented by Bates et al (1987). The procedure to obtain the coefficients α θ and β θ of the smoothed spline (with polynomial precision) is incorporated in appendix B for completeness. While the smoothed spline is not obtained through ridge regression, it can be expressed in an equivalent form, which will prove useful when evaluating relative improvements in the underlying surrogate model (section 3.3) ; representing the vector containing the values of the smoothed spline evaluated at the sampling locations by f θ , the approximation utilising the coefficients of the smoothed spline is given by Ψα θ + Pβ θ = f θ . This can be written analogously in terms of ridge regression (equation (2)), (Ψ + λI) α + Pβ = f, such that Ψα + Pβ = f θ . Equating the two expressions implies α = α θ and β = β θ and the corre sponding ridge parameter is obtained as To ensure a gradual decrease in smoothing and allow for sampling of smaller scales in case of low data uncertainty, λ (k) -values are iteratively scaled down with the following intermediate conditional parameter value (5) An example of the influence of λ is illustrated in figure 3. In the first iteration λ (1) is based on the ensemble averaged data scatter σ t and the surrogate model reflects the overall signal behaviour. As the value of λ (k) is gradually reduced  2 cos(20πx) is sampled with 50 randomly positioned points (N = 50), each superimposed with a local uncertainty (error bar width) σ t = |5γ t cos(100πx t ) sin(160πx t )| where γ t are random numbers between 0 and 1 with uniform probability. Measurement noise n t is defined as n t = 3τ t with τ t ∈ N (0, 1). The inset shows a close-up of the graph portion bounded by the rectangle. towards λ GCV (in the current example λ GCV = 3 × 10 −4 ), the underlying surrogate model tends to reflect smaller scale oscillations, though a degree of smoothing is retained.

Adaptive sampling objective function
The objective function J (i) (x s ) defined in Theunissen et al (2015) consisted of various individual objective factors converting the typical sampling criteria used by experimentalists in analytical formats. With new information added every iteration the objective function is iteratively updated (also here the subscripts (i) will be omitted for brevity) and maxima in the function then indicate new suitable sampling locations.

Curvature dependent sampling
The first criterion J C ensures an increased number of samples in regions of higher curvature, similar to mesh-refinement heuristics in computational fluid dynamics e.g. Mavriplis (1995). Curvature is estimated on the basis of the Laplacian Δ of the RBF-based surrogate model s, yielding the criterion J L ; The Laplacian of the TPS kernels is tabulated in table 1. For TPS2 kernels the Laplacian of the polynomial term vanishes. Whereas the surrogate model is constructed using the 2nd order TPS, its singularity in the second-order derivative is avoided by basing the curvature on the 3rd order TPS (∆q(x s ) = 2(β 5 + β 6 )). The objective function aims to indicate regions of higher curvature, not to be an accurate representation. As such, the use of multiple kernels is permitted.
To compensate for the detrimental effect of stochastic noise, J L is combined with a measure for local signal change where s MA (x s ) constitutes the moving average filtered response evaluated with a filter of approximately 30% domain width. While a moving average operation is computationally relatively inexpensive in one dimension when selecting neighbouring data points through a kd-tree, the process becomes too time consuming when extended into two dimensions. Advantage can be taken though of the existence of explicit analytical expressions for the integral of the TPS functions over a 2a × 2b rectangular domain; which for the case of the 2nd order TPS (n = 2) gives In one dimension the smoothed kernels reduce to The evaluation of the spatially smoothed surrogate hence reverts to a re-evaluation of equation (1), with a novel kernel and polynomial but readily available pre-determined coefficients α and β (recall that the TPS2 kernel is used to generate the surrogate model while the TPS4 kernel is only used in the calculation of curvature). Following a unity-based normalisation, referenced with subscript n, (·) n = [(·) − (·) min ]/[(·) max − (·) min ], both objectives are combined into the curvature objective J C ;

Uncertainty-based sampling
The second criterion ensures regions of higher uncertainty are adequately sampled by introducing the local level of data scatter as a component in the objective function. Regions in which accrued data shows larger scatter will undoubtedly demand additional samples not only from the perspective of data validation but also from a flow physics perspective. Additional samples will thus verify whether the increased scatter is due to an isolated outlier event or to a coherent region of higher turbulence intensity. While each variable is assigned one unique value at each sampling location, these values typically correspond to the average of N i independent measurements taken at each sampling station. Accordingly, the relative measurement error at the tth sampling location for 95% confidence level in variable (·) is given by where (·) t refers to the mean statistic (time-average) of N i values and σ (·)t the standard deviation in the data. The number of readings N i is typically set to 10 000 to retrieve reliable mean statistics or can be locally adjusted according the local data scatter to ensure a constant relative error across the domain. Subsequently, the objective component pertaining error can be expressed as J E (x t ) = ( (·)t ) n . The discrete J E values at the data points can also be interpolated using thin plate splines as per section 2.1 yielding a continuous function J E (x s ) for each of the considered variables.

Improvement-driven sampling
Sampling locations whose addition does not alter the surrogate model are evidently of low importance. Conversely, regions where the addition of new data yields a significant change in the previous underlying approximant merit further exploration. This thought is captured in the improvement objective J (i) To retain historic model change, the improvement objective combines the four most recent Ω's in a weighted average, with all weights w i−m currently equal to 1/4; Changes in surrogate model are however not necessarily linked to the addition of novel data, but can be due to alterations in the ridge-regression parameter, uncertainty in the model due to uncertainty in the acquired data and/or limitations in interpolation accuracy. This has been accounted for through three parameters quantifying each of these sources of variation in surrogate model; E λ (x s ), E σ (x s ) and E TPS .
1. With each ith iteration, the ridge regression parameter λ (k) is updated as new data is added. As such, the surrogate model can show variations which do not necessarily relate to an enhanced representation of the underlying signal, but are merely due to a different degree in smoothing, i.e. changes in ridge regression parameter. To quantify these changes, the matrix system in equation (2) is first rewritten as The first N rows and columns in I * constitute the identity matrix and the trailing M rows and columns are zero. With L = [Ψ P] evaluated at a location x s , the surrogate model value in the ith iteration is given as and the change in surrogate model due to this variation in In the second step of the derivation the Woodbury matrix identity was applied. It should be noted that E λ is not simply equivalent to s (i−1) (x s ) − s (i) (x s ) due to the addition of new sampling locations with every iteration. The difference between s (i−1) and s (i) thus encompasses a change in ridge regression as well as new data and would require re-evaluation of s (i−1) on the ith evaluation grid as this can change between iterations (section 4.2). Instead, in the final expression of equation (14), K −1 has already been calculated in the ith iteration to evaluate the TPS interpolation coefficients α (i) for the surrogate model s (i) (x s ). With each of the terms readily available, the change in surrogate model E λ (x s ) caused purely by a change in ridge regression parameter is easily obtained. 2. When performing point-wise measurements, multiple readings are taken at each station and local data scatter can be suitably modelled through a normal distribution N (0, σ). The local average of the readings is subsequently used to fit the surrogate model. When imposing a regression parameter to smoothen too rapid changes in underlying signal merely due to noise, the surrogate model becomes sensitive to the fitted data values as these might vary as a result of the underlying scatter. This subsequently introduces an uncertainty in the thin plate spline approximant, which has been dealt with in Wahba (1990). In this case the variation in surrogate model is caused by data uncertainty. Let σ t symbolise the data scatter in the tth sampling point and σ t the ensemble average of all σ t=1...N values. The suggested Bayesian confidence interval at 95% for the tth sample is given as where H tt is the tth entry on the diagonal of the influence matrix H.
i.e. the largest nearest-neighbour distance based on the avail- able data points within the domain, and a scaled native space norm F; The authors would like to stress that they do not claim E TPS to be a proven or mathematically underpinned error bound. The obtained error bound only provides a conservative guesstimate above which changes in surrogate model can no longer be ascribed to interpolation errors.
Values of Q (i) (x s ) exceeding the three thresholds; E λ (x s ), E σ (x s ) and E TPS , are finally considered sufficiently high to be changes in the surrogate model, genuinely caused by the addition of newly accrued data and are taken into account in the calculation of Ω (i) .

Iteration dependent exploratory sampling
New measurement locations will be prescribed by the objective function. However, each region within the domain should ideally be sampled at least once to ensure a representative objective function. To this extent oversampling of local extrema in objective function should be avoided, which in this work is achieved through a separation function J h (x s ). This function attains zero values at each existing measurement location and tends to unity away from the sample location. The objective function surrounding the points is thus given less importance, promoting exploration of the domain. The influence of each point should be radially symmetric and finite, making the compactly supported Wendland functions ideal kernels (Wendland 2004). Function J h (x s ) is therefore obtained through RBF interpolation as per section 2.1 with ψ(r) = (1 − ηr) 6 + 35(ηr) 2 + 18ηr + 3 , f = 1 and omitting the polynomial term. The cutoff function is defined as (ζ) + = ζ for ζ 0 and (ζ) + = 0 else. The radial extent of the kernel is controlled through η. Since increasingly smaller scales are envisaged within the iterative process (section 2.3), samples must be allowed to be positioned in closer proximity to each other as the iterations continue. This can be achieved by adjusting the Wendland kernel's support in function of the changing fill distance h s and iteration k in ridge regression (figure 4(a)); Starting with a radial extent equal to the fill distance, the kernel reduces asymptotically to h (i) s /10 ( figure 4(a)). The authors would like to draw the reader's attention to the definition of η as the scaling factor preceding h (i) s is altered with every ridge-regression iteration k while the fill-distance itself is updated with every sample addition (iteration (i)).
In line with the concept behind domain exploration, it would not be efficient to sample the domain in sequential regions. Instead, areas of interest must receive a higher number of samples, while avoiding global under-sampling. This is accounted for by introducing a secondary weighting to the objective function, particular to the tth sample point and dependent on the iteration number i; where i t designates the iteration number the tth sample point was introduced. The weighting starts from 0 and gradually reaches a value of 1 after N κ iterations ( figure 4(a)). Parameter N κ essentially dictates how fast a region can be re-sampled and has been set to N κ = 15 in the current study. The w t value at each of the sampling locations is subsequently interpolated again using the TPS2 kernel producing a continuous  (17) and (18) respectively setting N κ = 15 (equation (18)). (b) Illustration of the influence of the iterative weighting factor on the Franke function ( figure 8(a)). Insets depict the obtained sampling distribution corresponding to the curves in the graphs aside, superimposed on the contour levels of the underlying Laplacian of the signal.
weighting J w (x s ) for the objective function. J w will avoid recently sampled regions from being re-sampled, enforces domain exploration, yet gradually allows further investigation of already explored domain sections.
The effect of the new weighting term is demonstrated in figure 4(b) for a translated Franke function (see section 7), considering the first (k = 1) and last (k = 3) data smoothening iteration (see section 2.3). The objective function was based solely on curvature and seperation as per equations (11) and (17) respectively. Depicted is the coverage with increasing number of iterations, whereby five samples ( N S = 5) are added each iteration. A circular area of π/4N was centered on each sample. Considering N congruent, non-overlapping circles of diameter D in a square L × L domain such that √ ND = L (in the current approach the domain will be normalised, i.e. L ≡ 1-section 4.1), the ratio between the total area of the circles and domain is given by N πD 2 4 1 L 2 = π 4 . The coverage was subsequently calculated as the ratio between the unified domain area covered by all sampling points and the ideal maximum (π/4). In the first scale-adaptive smoothing iteration (k = 1) homogeneity in the sampling distribution (coverage tending to unity) is driven by J h (x) and J w is of lesser importance. Indeed, with increasing number of samples the omission of J w (or equivalently J w ≡ 1) yields a distribution quasi identical to that obtained with active J w . When η (i) is decreased in the third smoothing step (k = 3), the effect of J w becomes more prominent. Sampling locations become clustered in regions of high curvature when J w ≡ 1. A spatially varying J w on the other hand returns a distribution more populated in those regions, while simultaneously exploring their vicinity, which is reflected in the higher coverage attained.

Combined objective
At each sampling station multiple variables can be measured. In the case of two-dimensional LDA for example, the velocity magnitudes in two directions are obtained. Each of the N var variables will therefore have an associated J C (x s ), J E (x s ) and J I (x s ). For each variable these elementary functions are combined as The empirical offset constants κ were set as κ C = 0.5, κ I = 0.1 and κ E = 0.5. The definitive encompassing objective function is finally given by

Domain normalisation
While the above procedures are generally valid, a domain normalisation is introduced for simplicity, i.e. data locations are rescaled to the unit domain such that x s ∈ [0, 1] d . This is achieved by initiating the adaptive process by covering the d-dimensional domain with equispaced samples along each dimension. The minimum and maximum along each dimension thus define the extent of the observed flow region. The choice of initial number of samples is left user dependent. However, it should be noted that the total number of samples increases as a power of the dimension. Selecting for example 20 samples in each direction of a two-dimensional domain quickly yields a total of 400 initial samples, while considering a typical total of 1000 measurement locations. Moreover, the numerical assessments presented hereafter will show that the performances of adaptive sampling are relatively insensitive to the choice of initial sampling density. To that extent a modest 7-by-7 sampling grid is deemed a good compromise between effort and information content.

Adaptive evaluation grid
The evaluation grid defines the underlying nodes at which the surrogate models, and consequently the objective function defined in equation (20), is computed. While coarse grids hamper the detection of smaller structures in the objective function, a too fine grid spacing will needlessly augment the computation time. The authors circumvented this problem by partitioning the global domain into regular blocks and adapting the local evaluation grid spacing in each block in response to the local sample spacing. This concept is illustrated in figure 5. The Nyquist sampling criterion implies that an oscillation in a signal can only be reconstructed when sampled at least twice. The smallest nearest neighbour distance between samples is obtained as and this will be representative of the smallest measurable spatial fluctuation. Proper characterisation of any fluctuation at this scale will necessitate at least two evaluation grid points and for this reason the spacing is further reduced to fit five grid points. The local grid spacing of the equispaced evaluation grid is thus defined as δ = 1 5 (h t ) MA where (h t ) MA is the smoothed distribution of h t obtained by means of the moving averaging filtering described in section 3.1. Representing, in the ith iteration, the number of nodes in the jth block, which has a square size l × l, as N

Objective function guided sample allocation
While peaks in the objective function J O (x s ) (equation (20)) indicate the most suitable sampling locations, such a process would rapidly favour local extrema and hence clustering of the newly to-be-added measurement stations within an iteration. To avoid over-sampling, an alternative process has been implemented in this work, basing the sample locations no longer on simply the highest J O value, but on local peak geometry. Having defined the objective function, a threshold T JO is calculated based on the average of all detected local maxima ( figure 6(a)). Potential regions for further investigation are identified where J O (x s ) exceeds the threshold and sorted on the basis of their parameter Γ. This heuristic is evaluated for each region V j and quantifies the local extent of the objective function (in case of two dimensional domains this equals the volume underneath the objective function) and height (maximum value); Each region is then attributed n j = max 1, Γj Γ samples with Γ = 1 NS j Γ j and N S equals the number of new samples to be distributed each iteration. Of all the identified regions, only the N S regions with highest Γ are retained and each is subsequently sub-divided into 2 √ n j × 2 √ n j patches. For each region, patches with the highest J O average are then assigned one of the n j samples in their corresponding local centers of gravity (based on J O (x s )). This process tends to deliver an evenly-spread distribution of points across each region, while ensuring samples are spatially assigned in accordance with the objective function as illustrated in figure 6(b). The minimum sample spacing is thus also limited, allowing for a well-conditioned interpolation matrix.
The number of new sampling points N S to allocate each iteration is highly dependent on the underlying function. Generally, adding too many new samples each iteration will automatically enforce space exploration and negate any benefit of adaptive sampling. On the other hand, incrementing the number of samples by only one will yield a lengthy and computationally intense sampling process since the surrogate model needs re-evaluation with each newly accrued data point. The evolution in error and iteration number with N S is presented in section 7.1 and the authors found N S = 5 to provide a reasonable compromise while for one-dimensional adaptivity N S = 1 is more conducive (Theunissen et al 2015). Once the sampling locations have been identified, the order in which they are visited is optimised analogous to the Traveling Salesman Problem to minimise traverse time.

Surrogate model convergence
While the number of repeated measurement at a certain sampling location for a Gaussian process can be related to confidence level, local fluctuation level in the data and a user-imposed relative error (equation (12)), the required number of sampling locations to retrieve a reliable estimate of the underlying spatial data distribution is less evident and a priori unknown. Too few samples will yield a poor representation of the underlying data landscape, while oversampling is an inefficient use of computational and experimental effort. For this reason the adaptive sampling process should be considered converged when the addition of new data points yields a minimal change in surrogate model over a sufficient number of iterations. The experimentalist is therefore left with the problem of deciding for how many iterations the improvement must remain below a user-defined threshold. In this work a modified version of the adaptive stopping procedure of Theunissen et al (2015) has been employed.
While convergence in the surrogate model would imply a near-zero discrepancy, small model changes due to interpolation inaccuracies or modifications caused by uncertainty in the data should actually be discarded. Convergence in the surrogate model is thus essentially equivalent to demanding the previously introduced improvement to be below a certain level, which has been discussed already in section 3.3. The number of consecutive converged iterations N CCI in which the maximum improvement in surrogate model of each of the measured variables is negligible, i.e. Q (i) (x s ) max (E λ (x s ), E σ (x s ), E TPS ) ∀x s , has been defined as Parameter µ 1 quantifies the relative importance of the measurement uncertainty with respect to the overall sur- t=1 t as per equation (12). The contrast in surrogate model is essentially , adopting the previous definition of J ∆ (equation (7)). An overall large data scatter with respect to the underlying signal will result in high values of µ 1 and subsequently, an increase in N CCI . Such a more stringent convergence criterion will thus lead to an increase in number of samples placed and concomitantly an improved confidence in the obtained surrogate model's representativeness of the underlying signal. High values of µ 2 imply the surrogate landscape to contain few, yet very distinct (high amplitude) peaks. In the presence of such isolated features, a higher number of samples is required to encourage further domain exploration, which is analogous to a higher N CCI . Conversely, a relatively flat signal would correspond to low values of µ 2 and requires only few samples to be appropriately modelled. Note that N CCI is adjusted with every addition of new samples. In the ideal situation of no scatter in the data (µ 1 = 0) and uniform signal (hence J (i) ∆ = 0 and µ 2 = 0) the minimum number of consecutively converged iterations will equal N CCI = 11 (equation (22)), which was empirically found to be a suitable number.

Computational overhead
It is clear from the foregoing that adaptive sampling will be accompanied by an additional overhead compared to traditional sampling. In view of the numerical performance assessment in the following paragraph, figure 7 depicts the contribution of the different stages in adaptive sampling to the overall time (T total ) for the case where the underlying signal is defined by Franke's function. All processes pertaining adaptive sampling were implemented in the commercial software MATLAB (2014) on a 64-bit standard PC (Intel Core i7-4790, 3.6GHz CPU) with 16GB RAM. Regarding data accrual, at each sampling location N i = 10 4 samples are assumed to be collected at a typical data rate f s of 1 kHz to retrieve meaningful statistical moments.
Neither the allocation of samples (section 4.3), nor the establishment of an adaptive evaluation grid (section 4.2), nor convergence checks (section 5) signify a high computational load as seen in figure 7(a). The total time for adaptive sampling is dominated by data acquisition since in each iteration 50 seconds (= NS×Ni fs ) are allocated to data collection. With increasing number of samples the construction of the surrogate model and objective function start to take up a significant portion of the total processing time. The establishment of the surrogate model involves the evaluation of the interpolation coefficients α and β by means of an iterative conjugate gradient method (appendix A.2), the RBF matrix Φ, as well as a preconditioning and generalised cross-validation (section 2). The objective function involves the interpolation of the model onto an adaptively refined evaluation grid (section 4.2). As the number of samples grows, so do the involved matrices and hence the computational effort in their manipulation. Exactly these two subroutines constitute the additional overhead compared to the full factorial approach as illustrated by the red lines in figure 7. At approximately 600 samples, adaptive sampling takes approximately 20% longer than FF. However, as will be shown hereafter, with the benefit of obtaining consistent data quality and lower error levels compared to the structured sampling approach.
In figure 7(b) the initial sampling grid consists of 15 nodes in each direction compared to the 7 × 7 grid in figure 7(a). The evaluation grid thus contains a larger number of blocks (section 4.2) and, consequently, each block will accommodate fewer evaluation grid points. Indeed, for a given number of total samples, comparison between figures 7(a) and (b) indicates a reduction in additional computational effort. Moreover, the establishment of the surrogate model, which is independent of the evaluation grid, also shows to take slightly less effort, implying that optimised internal memory management may bring the computational cost down further. This is however beyond the scope of the current work.

Extension to higher dimensions
In describing the adaptive routines proposed, no limitations were imposed in terms of dimensionality of the problem and as such, the concept underlying the proposed adaptive process remains valid in higher dimensions. Some implementation aspects can be expected to require more attention. The preconditioning described in section 2.2 for example was tailored for surface fitting and hence no longer suitable when sampling in a volumetric domain. In that case, alternative RBF functions, such as the Wendland functions mentioned afore, may prove to be more appropriate. Especially given the rapid increase in data points, the sparsity promoted by these RBF kernels will be beneficial in reducing the computational overhead and negating the need for preconditioning. However, these kernels will require a proper definition of the inherent support radius. Also the optimal ridge-regression parameter λ (k) will necessitate a review of its estimation when changing RBF kernel (Schaback 2007) as this further influences the improvementdriven sampling. Given the vast application of RBF interpolation, the majority of these problematics has already been solved, making the adaptation of the presented strategy for volumetric sampling less cumbersome.

Numerical assessment
To test the performances of the proposed methodology, the adaptive sampling strategy has been tested on two analytical test functions; Franke's bivariate function (Franke 1979) and the numerically simulated flow field underlying a 2D lid driven cavity at a cavity-height-based Reynolds number of 100.
Franke's function was translated by 0.6 and 0.3 in x and y direction respectively and coordinates were rescaled by a factor five; f (5(x − 0.6), 5(y − 0.3)). The surface consists of localised peaks and a trough superimposed on a global exponential ( figure 8(a)). The lid driven cavity is chosen as it involves adaptivity to multiple velocity components and presents a complex flow field with a central vortex and two additional regions of recirculating flow in the cavity corners ( figure 11(a)).
With no uncertainty in the data, the µ 1 component in equation (22) equaled zero for all cases while µ 2 = 1 as the minimum in J D was 0. Accordingly, all adaptive sampling processes, each consisting of three iterations in smoothing parameter λ (k) , were considered converged after N CCI = 11 consecutive iterations in which the improvement in surrogate model was negligible.
Two initial structured grids were considered to construct the very first estimate of the surrogate model and initiate the adaptive process; a 7 × 7 and 15 × 15 grid. Once the adaptive approach converged, the uniform sample spacing in the FF was defined as h FF = 1/ √ N at which the function values were evaluated. Note that the number of samples with FF will always satisfy N FF = h −2 FF N. To juxtapose the 'goodness' of the adaptive sampling (AS) distribution with equispaced sampling, a Laplacian-based quality heuristic G was introduced defined as G = 1 The authors would like to re-iterate the fact that the total number of samples was not pre-defined, but automatically determined by the algorithm based on the stopping criteria discussed in section 5. Figure 8 gives an overview of the algorithm's performance with respect to the traditional full factorial sampling (FF) approach in the analysis of the shifted Franke function. As the number of samples imposed per iteration N S is increased, fewer iterations are needed to reach convergence ( figure 8(b)). Augmenting the initial number of samples from 49 to 225 shows the process to converge faster because of the better initial estimate of the underlying signal. However, the total number of samples continues to increase in a near linear fashion (figure 8(c)) and is nearly equal for both initial grids. The standard deviation in error (figure 8(d)) and maximum error (figure 8(e)) initially decrease, reaching asymptotic values equal for both FF and adaptive sampling. Up to N S ≈ 30 AS reduces both errors by a factor 3 compared to FF, independent of the initial sampling density used. As such, the adaptive approach shows to be independent of the initial sampling grid used, provided of course the initial surrogate model is at least to some degree representative of the underlying function. Also the sampling quality is consistently higher for adaptive sampling as the samples are located in appropriate regions ( figure 8(f)). This can be clearly seen in figure 9. With N S = 71 the regions of higher curvature are clearly visible through the heightened sample concentrations. However, these regions become increasingly oversampled when N S is enlarged as the interpolation accuracy improves due to the reduced sample spacing. Driven by the objective J I (x s ) (section 3.3) additional samples are placed to further enhance the surrogate model. For N S = 5 the denser initial (structured) grid provides an overall, better initial estimate of the underlying function ( figure 9(d)), yet, the function's extremities are sampled to a lesser amount compared to the sparser initial grid, yielding a comparable maximum error (figure 8(e)). For both initial grids tested, the addition of 5 new samples ( N S = 5) with each iteration shows to result in an appropriate sampling. In view of the time required to perform the experiment, a total number of samples N ∼ O(10 3 ) is desirable, in which case the typical number of samples to be added each iteration should be less than 20 according to figure 8(c). Comparing the error distributions (figure 10) it can be seen that little is gained going from N S = 5 to N S = 11 despite a near doubling of the total number of samples. Both AS and FF have the largest errors in the regions of largest spatial variance, but already for N S = 5 AS reaches error levels comparable to those achievable with N S = 11. It is as such proposed that N S = 5 is most suitable. For the lid driven cavity figures 11(b) and (c) depict the evolution in maximum and random error with respect to the total number of samples adding with each iteration five new samples ( N S = 5). It should be noted that the sampling is based on surrogate models of both the horizontal (U) and vertical (V) velocity component. All heuristics, calculated for each of the components, are therefore arithmetically averaged and displayed in figure 11. Initially the error plots for both AS and FF approaches are very similar. However, structured sampling yields a higher variation in maximum error and increasing the number of samples does not guarantee a reduction in maximum error ( figure 11(b)). The adaptive process on the other hand yields a typical reduction in (δ ) max and random error by a factor two compared to FF. In other words, for FF to achieve the same error level as AS, at least double the number of samples are required. Even then, based on figure 11(b) this still does not ensure FF to attain the same error level.
The smoothing iterations can be clearly identified in figures 11(b) and (c); for an initial grid of 7 × 7 samples the value λ (1) , based on data scatter, is used for the first 38 iterations ( N = 7 × 7 + 38 × 5 = 239). Although the maximum error for λ (2) remains largely constant, this tendency does not reflect the local improvements in surrogate model. As can be seen in figure 11(c), the random error continues to decrease. In approximately the 144th iteration ( N = 7 × 7 + 144 × 5 = 769) the final smoothing is applied.
In this stage smaller flow scales are considered, which, especially for j N (0) e,j = 7 × 7 yields a sudden reduction in maximum error.
Increasing the initial sampling density allows AS to reach lower maximum errors during the converging process. However, as can be seen from figures 11(b) and (c), once the process converges, final error levels are identical. Higher sampling densities with AS are observed in regions of strong gradients, notably the corner vortices, the top of the cavity and the central inflection point ( X ≈ 0.6, Y ≈ 0.75). This clearly leads to a drastic reduction in local error (figures 11(e) versus (f)). With FF exactly these regions consistently present the highest errors.

Experimental assessment
The final assessment consisted of the application of the adaptive methodology to an experimental test. The case of an open cavity flow with a length-to-depth (L/H) ratio of 2 was considered at a Reynolds number Re H ≈ 37 × 10 3 . Experiments were performed in the open-return low-speed wind tunnel  of the University of Bristol at a flow speed of approximately 10.33 m s −1 and a turbulence intensity of approximately 2%. 2D2C Particle Image Velocimetry measurements were performed on a cavity 110 mm long and 55 mm deep using a Dantec Dynamics system consisting of a 200 mJ Litron laser and 4MP FlowSense EO camera. A sketch of the setup and exemplary PIV image illustrating the location of the cavity walls is provided in figure 12. The 2072 × 2072 px 2 sensor corresponded to a field of view of 115 × 115 mm 2 and images were analysed using initial interrogation windows of 65 × 65 px 2 , recursively reduced to 17 × 17 px 2 . A mutual window overlap of 75% yielded a final vector spacing of approximately 0.2 mm. In total 1000 instantaneous velocity fields were considered ( N i = 10 3 ) in the calculation of the mean and fluctuating velocity components, (u, v) and (u 2 , v 2 ) respectively, where The latter stresses then enabled the determination of the local standard deviations, e.g. σ u = Ni Ni−1 u 2 , to establish relative measurement estimates as per equations (12) and (15). The extracted PIV data subsequently served as input into the adaptive sampling process adopting linear data interpolation from the structured PIV data to the intermediate locations.
The adaptive process was initiated by taking measurements at 49 (7 × 7) equispaced samples. Each iteration five new samples ( N S = 5) were added. Typical LDA traverse increments are in the order of 0.1 mm. In terms of the unit domain, this implies that sample locations must be located on multiples of 0.1 115 . This was taken into account when attributing new data positions. Useful data can only be extracted within the flow region. To this extent cavity walls are to be excluded when attributing new sampling locations. This was accounted for by setting the objective function (equation (20)) to zero within the corresponding regions as evidenced when overlaying the FOV with the predicted sampling locations ( figure 12(c)).
The number of consecutive iterations to establish convergence in each of the smoothing phases through ridge regression was automatically set to N CCI = 16, N CCI = 15 and N CCI = 15 respectively. A higher concentration of samples is allocated within the cavity compared to the outer flow region. A zone of recirculation can be noted in the upstream corner of the cavity as visualised by the streamlines and this region, as well as the upstream and downstream cavity stagnation points are adaptively allocated more samples. Contrary to the numerical cases, which are void of any fluctuation in signal, the distribution of new samples is now influenced by the local data uncertainty, i.e. J E . This is manifested when samples are superimposed onto the product of local Reynolds stresses u 2 × v 2 ( figure 12(d)). Whereas the adaptive approach yields a consistent reduction in maximum error with increasing number of samples ( figure  12(e)), the standard approach shows stronger fluctuations. Not only does AS attain error levels inferior to those of FF (figures 12(e) and (f)), but, as per the numerical cavity flow, adding more data does not ensure either the maximum, or the rms error of FF to reduce. On the contrary, adding samples can actually increase the maximum error as samples are not optimally placed in flow regions of spatially rapid changes such as the cavity corners. This leads to poor representations of the flow as manifested in the error distribution plots (figures 12(g)-(j)). In the bulk of the flow where spatial gradients in time-averaged velocity are low, differences between full factorial and adaptive sampling are marginal. However, in the close vicinity of the cavity walls, adaptive sampling clearly yields a representation in better agreement with the PIV results. Both this test case and foregoing numerical assessments corroborate AS to be more reliable and efficient in the measurement of an a priori unknown flow field; 1. Error levels in the reconstructed signal with data sampled through AS are lower than those sampled by means of FF. 2. Without knowing the underlying signal, it is impossible to preset the required number of samples to collect. In the FF approach N is thus rendered a critical user-defined parameter. The adaptive routine on the other hand automatically stops when sufficient convergence is reached and does not require the experimentalist to make this decision. 3. The observed instability in error with FF, from an experimental perspective, is a highly unwanted anomaly. When using the traditional equispaced sampling, the assessments have indicated that the experimentalist cannot be certain that the obtained result would be more accurate if more data points were added. He/she is therefore required to gradually increase the number of samples and make decisions based on observed changes in the reconstructed flow field. Such a process however intrinsically replicates the concept behind adaptive sampling. 4. Adaptive sampling shows a more consistent decrease in error (both maximum and rms) with total number of samples. Aborting the adaptive process before reaching convergence will subsequently still yield an overall result more reliable than FF.

Conclusions
With the aim of improving experimental efficiency when surveying e.g. flow regions by means of point-wise metrologies, an adaptive strategy has been proposed to properly sample 2D areas of interest in an automated manner. The methodology incorporates sampling criteria based on local signal curvature, improvement, uncertainty and domain exploration. In addition, initial smoothing is applied and progressively reduced as to mitigate the effect of local data jitter due to measurement noise. The different processes are elaborated within this paper. Thin Plate Splines finally allow the reconstruction of the surrogate model onto a self-adapting evaluation grid. Automated stopping criteria have been implemented to imply convergence is reached and the numerical and computational implications are discussed within the document. The performance of the adaptive sampling strategy is assessed and juxtaposed with results obtained by traditional, equispaced sampling (full factorial sampling) based on Franke's function and the lid driven cavity flow. In addition, both approaches are applied to an open cavity experiment. All test cases attest regions of interest to be consistently attributed a higher number of samples through adaptive sampling. The optimum number of samples to add with each iteration was found to be 5. This yields a drastic reduction (typically factor 2-3) in error compared to traditional full factorial sampling, independent of the initial number of samples. To reach the same level of error as the adaptive approach, traditional factorial sampling thus requires at least twice the number of samples. More importantly, whereas adaptive sampling yields lower error levels with increasing number of samples, full factorial sampling does not guarantee such a tendency. However, without any a priori knowledge of the flow field, the experimentalist is unable to determine an appropriate number of samples when adopting equispaced sampling. Traditional sampling thus demands the experimentalist to frequently reevaluate the reconstructed signal and augment the number of samples until presumed convergence is reached. Such a process intrinsically replicates the concept behind adaptive sampling and, combined with the fact that the automatic stopping criteria implemented have shown to be effective in the assessments provided, corroborates the adaptive approach to be more reliable and efficient in the measurement of an unknown signal/flow field. ξ j = Bu and η j = Bv. The preconditioning matrix C then follows from approximating the Laplacian in terms of the derivatives of the polynomials; (36ξ 30 η 30 + 12ξ 12 η 12 + 12ξ 21 η 21 + 36ξ 03 η 03 ) j T j ≈ N∆ j=0 (36B 1, * uB 1, * v + 12B 4, * uB 4, * v + 12B 3, * uB 3, * v + 36B 2, * uB 2, * v) j T j (A. 2) The notation B i,* refers to the ith row of matrix B including all columns.
(A.3) Superscript (q) indicates the iteration number. The process is stopped when convergence is reached in α (q) or the number of iterations exceeds a maximum.