Adaptive Isogeometric Digital Height Correlation: Application to Stretchable Electronics

A novel adaptive isogeometric digital height correlation (DHC) technique has been developed in which the set of shape functions, needed for discretization of the ill‐posed DHC problem, is autonomously optimized for each specific set of profilometric height images, without a priori knowledge of the kinematics of the experiment. To this end, an adaptive refinement scheme is implemented, which refines the shape functions in a hierarchical manner. This technique ensures local refinement, only in the areas where needed, which is beneficial for the noise robustness of the DHC problem. The main advantage of the method is that it can be applied in experiments where the deformation mechanisms are unknown in advance, thereby complicating the choice of suitable shape functions. The method is applied to a virtual experiment in order to provide a proof of concept. A second virtual experiment is executed with stretchable electronics interconnects, which entail localized buckles upon deformation with complex kinematics. In both cases, accurate results were obtained, demonstrating the beneficial aspects of the proposed method. Moreover, the technique performance on profilometric images of a real experiment with stretchable interconnects was demonstrated.


Introduction
Digital image correlation (DIC) is nowadays an almost indispensable technique in experimental mechanics [1,2]. Subsequent images of a test specimen taken during an experiment are correlated to determine displacement fields. Initially, DIC was developed to analyse in-plane, two-dimensional displacements. Meanwhile, several extensions of the method have been developed. The method has been extended to 3D by correlating images from two cameras in stereo to obtain the displacement of a surface in three dimensions: stereo-DIC [1,3,4]. Furthermore, the method has been advanced to a fully 3D method: digital volume correlation (DVC), where also the internal kinematics of the sample is tracked, rather than solely the surface deformation [5,6]. This is achieved by using for example X-ray tomography scans instead of planar photographic images.
Another recent development concerns the correlation of profilometric images in order to identify both in-plane and out-of-plane deformation fields: digital height correlation (DHC). In Figure 1A-C, examples from literature are shown where this technique is applied. In all three cases, accurate results were obtained on the microscale. Han et al. included mode-I crack displacement fields in the algorithm to describe the kinematics of a propagating crack in a glass specimen [7]. Bergers et al. included the function describing the shape of a single-clamped beam in the algorithm, which was required to calculate the curvature of a bending microbeam with a resolution of $ 1.5 Á 10 À 6 μm À 1 [8]. Neggers et al. used globally defined, continuously differentiable polynomial functions to accurately capture the local strain and curvature fields of a bulging membrane under pressure [9], which enables the measurement of the local plane strain and biaxial elastic moduli within $2% accuracy. The results in these cases were convincing; however, in all cases, the discretization of the DHC problem (necessary because of the intrinsic ill-posed nature of a DHC formulation) was adapted to the specific mechanics of the considered experiment. However, in most mechanics problems, it is not possible to assess the kinematics of the unknown displacement field a priori, for instance, in the case of Figure 1D, where a copper stretchable electronic interconnect delaminates from the rubber substrate and buckles in specific local areas, which is an active field of research [10,11]. Therefore, a more generic DHC framework is called for, which preferably autonomously adjusts to the kinematics, without using prior knowledge.
For the case of in-plane, two-dimensional DIC, a novel adaptive isogeometric global digital image correlation (iso-GDIC) scheme was recently developed [12]. Isogeometric shape functions for the discretization, both of the sample domain and the unknown displacement field, were used, i.e. non-uniform rational B-splines (NURBS). It was shown that this type of shape function is versatile and able to capture a wide range of kinematics. This type of shape functions is used increasingly in a DIC setting, both in 2D [12][13][14] and 3D (stereo-vision) [15][16][17]. NURBS have proven to be less sensitive to noise than finite elements (FE) [18]. Moreover, NURBS originate from CAD-modelling and are able to describe many shapes exactly. In the case of stretchable electronics, this is of particular interest. In the three examples shown in Figure 1A-C, a rectangular mesh for the discretization of the DIC problem sufficed for the purpose of the correlation procedure. Yet, to be able to describe both the complex shape of the interconnects and the buckle pattern, which occurs mainly at the edges of the interconnect, to accommodate the sample's edges, a more advanced mesh must be used. Furthermore, the continuity of NURBS functions across the element edge is adjustable by inter alia choosing the polynomial order of the shape functions. This is beneficial in the case where one, for example, is interested in the calculation of the curvature field, which requires at least C 1 continuity. In [12], the DIC algorithm was combined with an adaptive refinement procedure, in order to autonomously optimize the shape functions. The number of degrees of freedom (DoFs) thereby remains limited, which is beneficial for the noise robustness. Refinement is executed locally, such that a finer discretization is only used in areas where this is necessary for accurately capturing the kinematics of the unknown displacement field, thereby compromising noise robustness [19], while retaining a coarser mesh in areas where the kinematics allow for this, preserving noise robustness. Furthermore, problem-specific choices on the discretization are not required in advance. In this paper, a generic nearly autonomous DHC framework is developed, which requires adaptation of the 2D adaptive iso-GDIC formulation towards quasi-3D: DHC. An advanced meshing framework is thereby required, which is generic for a myriad of shape function types, polynomial orders and mesh generating interfaces. In comparison to most global DIC formulations, the proposed method is less user dependent, since the most important choice, for the set of shape functions, is automated. The potential of the novel method is demonstrated on both virtual and real experiments with interface delamination of stretchable interconnects. In the "Methodology" section, the methodology used is explained: first NURBS shape functions are introduced, the discretization of the specimen shape is shown, and also, the refinement procedure of the discretization is clarified. Furthermore, the DHC algorithm is defined. In "Demonstration: Virtual Experiment" Section, the novel adaptive isogeometric global digital height correlation (iso-DHC) technique is applied to two virtual experiments in order to provide a proof of concept. Also, noise is included in the analysis. The method is applied to a real experiment with stretchable electronics in "Experiment: Application to Stretchable Electronics" Section. Finally, conclusions are drawn in "Conclusion" Section.

Methodology
In this section, the methodology for digital height correlation is detailed. First, the parametrization procedure for the sample geometry is clarified in "Isogeometric shape functions and parametrization" Section. The required shape functions are thereby introduced. These are also used for the DHC algorithm, which is addressed in "Digital height correlation algorithm" Section. Finally, the adaptive refinement procedure is described in "Adaptive refinement" Section.

Isogeometric shape functions and parametrization
The shape functions that are used to parametrize and regularize both the sample geometry and the unknown displacement field are NURBS: non-uniform rational B-splines. This type of shape function originates from CAD (computer aided design) modelling but is used increasingly in the computational analysis of mechanical problems: isogeometric analysis (IGA) [20]. Both B-splines and their generalization, NURBS, have been used in digital image correlation methods [12,14].
In this work, the CAD representation of the sample is used directly in the DIC analysis. The procedure to generate the initial discretization is explained by using the commercial CAD program Autodesk AutoCAD [21]. The geometry parametrization could also have been retrieved from other CAD programs. An image of the undeformed sample, the reference image, is loaded in AutoCAD to act as a reference for the creation of the mesh. Thereafter, a NURBS surface is inserted with the desired order and number of degrees of freedom in both directions. Associated with this surface are control vertices, or control points, that can be translated to make the surface fit the sample geometry, utilizing the inserted background image. The geometry parametrization is described bȳ where N i ðξÞ are the two-dimensional NURBS shape functions defined on parametric coordinatesξ,p i are the control points, andxðξÞ represents the mapping of the created surface, or mesh, to the physical domain. Since the DHC algorithm concerns images that are defined on a regular grid of pixels, it is necessary to evaluate the shape functions at the pixel locations, instead of the local element grid resulting from equation 1. This can either be done by interpolation, or, faster, by a nearest neighbours search algorithm, which couples the shape function value of the nearest local coordinate to the pixel coordinate. A choice for the second option is made because the loss of accuracy was found to be negligible for a fine local element grid, while still keeping a considerable gain in speed. This mapping procedure is illustrated in Figure 2. The fundamental building block of a NURBS surface is the univariate B-spline [20], which is a piecewise polynomial function of order p that is defined over a knot vector Ξ = {ξ 1 , ξ 2 , …, ξ k }, where each knot determines an element boundary in the domain. For NURBS, these knots are not necessarily uniformly distributed. Also each knot can occur more than once; the continuity of the shape functions across an element border (knot location) is controlled by the multiplicity m of the knot: C p À m . The number of shape functions (and hence number of degrees of freedom, DoFs) is determined by the number of knots k and the polynomial order ( Figure 3).
In the isogeometric GDIC approach, it is necessary to reconstruct the CAD shape functions and geometry in the DIC code. Here, this information is obtained by extracting the required data from the AutoCAD file (drawing exchange format): (i) the chosen polynomial order of the NURBS surface; (ii) the unique knot values in both directions and their corresponding multiplicities; and (iii) the control points and, possibly, their weights. Using the knot information, the Bézier extraction procedure [22] is applied to compute the element extraction operatorsC ē . With these operators, the spline basis functions on an element can be constructed as a linear combination of a canonical set of shape functions, in this case, the Bernstein polynomialsB: This process is illustrated in Figure 4. The reader is referred to [22] for details on the extraction procedure. The resulting shape function N i ðxÞ is composed of the contributions from all elements: Some of the resulting shape functions for second-order NURBS are plotted in Figure 3.
It is emphasized that this extraction is not restricted to B-splines or NURBS but can also be used to construct, e.g. T-splines [23]. From the perspective of the DIC algorithm, this extraction process provides a unified interface for the implementation of a variety of spline technologies from different CAD interfaces.
The shape functions are not only important for the parametrization of the sample geometry but also important for the discretization of the DHC problem, as will be shown in the next section.

Digital height correlation algorithm
The shape functions are used not only for the parametrization of the specimen geometry but also for regularizing the displacement field,ŪðxÞ, in DIC, i.e. in the correlation of the images of a deforming test specimen. The first image, f, is After this mapping, the shape functions are known at a local element grid of coordinatesx e ðξÞ, as indicated in black in the zoom of an element in the right image. For the DHC algorithm, it is required that the shape functions are evaluated at the pixel locations, indicated by the blue grid in the same image. Therefore, a nearest neighbour search algorithm is employed, which for each pixel centre inside an element (blue dots) finds the nearest element grid point (green dots) and assigns the value of the shape function at this point to the pixel. Note that for illustration purpose the pixel grid and element grid are depicted coarse. In reality, the element grid is significantly finer than the pixel grid, such that the loss of accuracy of the nearest neighbour mapping method is negligible In regular, two-dimensional DIC, the images are characterized by the grey-scale intensities measured at the pixel locations, and the corresponding brightness at the material points is assumed to remain constant upon deformation of the underlying material, i.e. brightness conservation holds: f ðxÞ À g∘ΦðxÞ ¼ rðxÞ≈0; (4) ΦðxÞ ¼x þŪ xy ðxÞ; where rðxÞ is the residual image, which is zero in the absence of noise, andΦðxÞ is a vector function which maps the reference coordinatex to the deformed coordinate. Note that throughout this article the same notation is followed as in Ref. [24], i.e. the coordinatex refers to the (Lagrangian) reference coordinates, while the deformed coordinates are consistently expressed using the mapping functionΦðxÞ. The residual is minimized to achieve optimal correlation, thereby obtaining the two-dimensional, in-plane displacement field U xy ðxÞ (see, e.g. [12,25]). For those cases where the out-of-plane deformation field WðxÞ is also desired, the DIC algorithm can be extended to digital height correlation (DHC) [9]. In that case, the images are not defined by the grey-scale intensities, but each pixel contains a quantitative measurement of the height of the surface, obtained with, e.g. optical profilometry, atomic force microscopy, or scanning tunnelling microscopy. The conservation relation therefore transforms to surface height conservation, i.e.: f ðxÞ À ðg∘ΦðxÞ þ WðxÞÞ ¼ rðxÞ≈0; ΦðxÞ ¼x þŪ xy ðxÞ; whereŪ xy ðxÞ is now the in-plane component of the total, three-dimensional displacement fieldŪðxÞ ¼ UðxÞē x þ VðxÞē y þ WðxÞē z , which is a function of the two-dimensional position vectorx ¼ xē x þ yē y . Identifying the displacement field that will satisfy equation 7 is an ill-posed problem, which deteriorates through the inevitably present additional noise field. Therefore, DIC methods approximate the true displacement field with a field represented by a finite and limited set of unknowns, and where is a column of degrees of freedom (DoFs), i.e. = [a 1 , a 2 , …, a 3n ] T , with n DoFs for each of the three components of the displacement fieldŪðxÞ. Applying more pixels per DoF allows for attenuation of acquisition noise (e.g. [26]), provided that the discretized displacement field can adequately describe the true displacement field.
As is commonly done in DIC, the displacement field is approximated as a linear summation of DIC basis functions,φ i ðxÞ: Note, however, that these basis functions are threedimensional vector-valued fields. For this purpose, the After multiplication with extraction operatorsC e , the shape functionsN e on an element e (with parametric coordinate ξ) are obtained (right). Note that the extraction process also involves a mapping from parent coordinate b ξ to local, parametric coordinate ξ same NURBS shape functions, N j ðxÞ , that are used for parametrization of the sample geometry (see "Isogeometric shape functions and parametrization" Section), are implemented. Note that the NURBS functions are two-dimensional, scalar-valued functions, and each NURBS function is used three times to describe the three components of the displacement field with independent DoFs: Similar to regular 2D-DIC, a cost function Ψ(a) is defined as the L 2 (Ω) norm of the residual and a minimization problem is formulated to solve for the DoFs and thus the optimal approximate solution to the displacement field: The conventional derivation of the DIC solution scheme to determine the DoFs a in Eq. 13 involves, first, linearization of the conservation equation, followed by an iterative optimization algorithm (usually Gauss-Newton), resulting in a two-step linearization with a number of implicit assumptions. This two-step linearized system of equations is then iteratively solved to retrieve the optimal unknowns . In Ref. [24], however, it was demonstrated that the non-linear conservation equation can be minimized in a consistent mathematical setting, yielding a one-step linearization, thereby highlighting the implicit assumptions made. Here, the same one-step linearization is followed, resulting in a system of equations that is iteratively solved for the unknowns where is the iterative update of the DoFs. As argued in Ref. [24], the tangent operator contains three terms, of which the second term is zero because the adopted basis is here linearly independent, while the third term is neglected as it contains the second gradient of the image making it highly sensitive to measurement noise. In that case, only one tangent operator term remains: whereby the right-hand member of Eq. 14 is given by Here,Ḡ is the true image gradient, i.e. the gradient of the image g evaluated at the deformed coordinates: The term in the z-direction is added for DHC, to correctly deal with the out-of-plane displacements. As detailed in Ref. [24], using the deformation gradient tensorF , the true gradient can be related to the gradient of the backtransformed image, Therefore, to simplify the true gradient to the one typically found in literature, first, small deformations are assumed, i.e.F T ≡Ī . Second, grad e g ð Þ is replaced with grad f ð Þ, which has been justified in the literature because e g is updated at each iteration and converges towards f (see, e.g. [27]), For the present case, this simplification to grad f ð Þ was found to have a negligible effect on the accuracy, in agreement with the guide lines given in Ref. [24]. Therefore, grad f ð Þ was implemented to reduce computational costs; however, extension to the true gradient is trivial. In order to determine in Eq. 19, the surface height values in the deformed image gðxÞ need to be determined at the deformed planar positions , which requires interpolation. Interpolation is a source of error and to minimize its impact a cubic spline interpolation scheme is here implemented, as suggested by Schreier et al. [28].

Adaptive refinement
In order to be able to accurately describe the kinematics of the displacement field, the regularized displacement field should be sufficiently rich, i.e. should have enough degrees of freedom. However, if too many DoFs are used, the solution becomes highly sensitive to noise [13,29]. It is therefore important that the number of degrees of freedom is balanced. In certain cases, where the kinematics of the problem is known in advance, it is possible to make a good estimation on the number and distribution of DoFs (i.e. the configuration of the mesh) to obtain accurate results, e.g. in the examples shown in Figures 1A-C. However, in some cases, including the experiments with stretchable interconnects considered in this work, shown in Figure 1D, it is more difficult to assess the kinematics in advance. Furthermore, displacement fields might be rather complex and exhibit localized features, which calls for a more detailed regularization in specific areas.
With an adaptive refinement algorithm, the mesh can be optimized autonomously. This can be achieved by either p-refinement, where the polynomial order of the shape functions is elevated in designated elements of the mesh [30], or by h-refinement, where the elements themselves are refined [12], or isogeometric k-refinement [20]. Both h-refinement and p-refinement can be done in an adaptive fashion, where the algorithm autonomously determines in which area refinement is required based on the error or residual rðxÞ in that area. In that way, the solution is not dependent on user experience. Furthermore, both methods can be implemented efficiently, since only the refined shape functions have to be added in case of p-refinement, or substituted in case of h-refinement, in the columnφðxÞ that contains all shape functions. Hence, there is no need to rebuild the entire set of shape functions. Which type of refinement is used is a matter of taste. Nonetheless, in the case of the three-dimensional deformation of structures, including buckles and delamination, the curvature of the material may be of interest, which requires second-order derivatives of the displacement field. Therefore at least C 1 , continuity across element borders is required. Since for NURBS shape functions continuity across element borders is C p À m , with p the polynomial order and m the multiplicity of knots at the element border, C 1 continuity can be achieved by using second-order B-splines with multiplicity 1. With prefinement-based finite elements only C 0 continuity can be achieved, which is the reason why in this work it was decided to employ isogeometric analysis. In combination with h-refinement, this allows for the refinement of the computational grid, while preserving the necessary continuity properties.
The approach adopted in this paper is a hierarchical refinement [31,32] scheme, identical to the technique used in [12]. In this method, multiple bases of shape functions are defined, which represent subsequent levels of uniform refinement. If refinement of a certain shape function is desired, this shape function is replaced by shape functions from the underlying basis that lie in the same support area. The result of this concept is that refinement occurs in a local fashion, in contrast to knot insertion, where the tensor product structure induces refinement of an entire row and an entire column of elements. The idea of hierarchical refinement is depicted in Figure 5.

Refinement indicator
The selection of shape functions for refinement is based on the residual rðxÞ, since the residual is also used as an error estimator in the DHC procedure itself. The residual gives full-field information, which makes it possible to distinguish between areas where correlation of the deformed image to the reference image is successful, and areas where it is not possible to approximate the displacement field accurately. For each shape function, an averaged value of the residual in the area of its support is calculated and weighted: First, the residual is weighted with the shape function N j itself, in order to couple the residual in a certain area to the shape function with the largest influence in that area, i.e. a larger value. Furthermore, this scaling assists in preventing larger shape functions always being favoured for refinement at the expense of basis functions with a smaller support area. Additionally scaling with the mean , which is a measure for contrast, is applied. This is because the residual is influenced not only by non-exact correlation but also by changes in the contrast in the pattern of the sample. Imagine two neighbouring pixels with a different value (either grey scale intensity or height). Now, correlation is slightly inaccurate, and the value of one pixel is assigned to the other pixel in the back-transformed image e g. Since the residual is defined as the difference between the original value of this pixel and the value in the backtransformed image, the residual will be larger if the difference in value between the two neighbouring pixels is larger, i.e. if the contrast in that area is larger. To compensate for this the refinement indicator C j is scaled with the relative mean intensity gradient δ f,global /δ f,j , where δ f,j represents the contrast in the area of support of shape function j and δ f,global in the entire region of interest. Finally, scaling with the range of pixel values f, here the range of height values, is implemented. This makes it possible to base the refinement criterium on the level of acquisition noise of the images. A shape function is selected for refinement if the refinement indicator exceeds a certain threshold: C j > T (see Figure 6). This threshold is set to C is the average of the refinement indicator C j of all shape functions, and σ is the standard deviation. This threshold ensures that only shape functions are refined for which this is necessary, i.e. for which the refinement indicator is significantly large with respect to the other shape functions. If the differences in refinement indicator between the shape functions are too small, e.g. when the displacement field is homogeneous and no refinement is required, this threshold assures that no shape functions are selected for refinement. The value T nl is an absolute threshold which corresponds to the noise level of the images, which has to be determined for each set of experimental images separately, such that refinement is not triggered by artefacts caused by image noise. This choice for the threshold is in correspondence to Ref. [12].
Since the correlation becomes highly sensitive to noise if the number of degrees of freedom becomes too large compared with the number of pixels, a minimum is set for the number of pixels within an element. This threshold is based on the correlation length, ℓ c , which is a measure for the in-plane length scale of the pattern, i.e. the average size of the pattern features. At least several pattern features should be present inside an element, otherwise noise is dominant for correlation. Therefore, the minimum element size is set to 10ℓ c × 10ℓ c , corresponding to Ref. [12]. The correlation length is determined automatically for each experiment. If an element becomes smaller than this threshold, the shape functions corresponding to this element are excluded from further refinement.
The mathematical formulation allows to use a different set of shape functions for each direction and thus only refine the shape functions for one direction, e.g. only for displacement in the height direction. Especially in this particular example of buckling of a stretchable electronics interconnect a different set of shape functions in the out-of-plane direction and refinement of only this set would make sense, considering the more complex nature of the out-of-plane deformation with respect to the in-plane deformation. However, this buckling case is a specific example and implementation of a scheme that relies on the known kinematics of the particular problem (by either choosing a different set of shape functions for one direction or refining only in this direction or applying both) would imply a loss of generality. Furthermore, the residual field is a result of the correlation of the displacement fields in all directions, and hence, it is not possible to distinguish between the accuracy of the correlation in the different directions. Therefore, the refinement is carried out in the shape functions in all three directions, and thus, the same set of shape functions is used for all x, y and z directions.

Demonstration: Virtual Experiment
The developed adaptive isogeometric digital height correlation algorithm is applied in a virtual setting in order to demonstrate the method. First, a proof of concept is given with a virtual displacement field that Figure 5: A graphical representation of the hierarchical refinement process: in the left figure, the initial mesh is depicted, plotted on top of the undeformed image of the stretchable interconnect. One of the shape functions is shown, and the maxima of all shape functions are indicated by blue dots. For this example, only the depicted shape function is selected for refinement, and the resulting refined mesh is shown in the middle figure. Note that all elements that form the support of the initial shape function, i.e. the top two-by-two elements in the left figure, are refined. Again, one shape function is shown, and the maxima of all shape functions are indicated by blue dots. The refinement process is repeated for the shape function of the middle figure, and the result is shown in the right figure represents a localized buckle pattern. Here, a rectangular mesh is used. It is demonstrated that the developed adaptive iso-DHC method is easily used for different orders of the NURBS shape functions. Also the influence of noise is investigated. In the subsequent example, a more complex geometry is used: a stretchable electronics interconnect, which requires the advanced meshing framework introduced in "Isogeometric shape functions and parametrization" Section.

Localized buckles
In this experiment, a virtual height profile is analytically deformed. The height profile is in this case generated analytically and contains both coarse and fine in-plane features, making it suitable for DHC analysis [33]. The applied out-of-plane displacement field represents a localized buckling pattern, with two sinusoidal peaks of different size, of which one is pointing upwards and the other pointing downwards (see Figure 9A). The in-plane displacement is zero in both x and y direction in the entire domain and therefore not discussed in the results. In the next example, a virtual experiment will be shown where also the in-plane displacement is taken into account. The out-of-plane displacement is applied in four increments in which the amplitude of the sinusoidal peaks increases from 1 to 4μm. The reference image f, the final image g 4 and intermediate images g 1 and g 2 are shown in Figure 7. The developed DHC method is applied, using second-order (p = 2) NURBS shape functions. The nearly autonomous refinement algorithm adaptively refines the mesh in the areas where refinement is required. The resulting meshes are also shown in Figure 7.
Since refinement is based on the residual field, it is interesting to analyse these fields. In Figure 8, the residuals are shown for each refinement step, corresponding to the meshes shown in Figure 7. In the first figure, it is observed that the residual is high in the area where the peaks occur, this means that the original set of shape functions, with the mesh of Figure 7A, is not able to capture the kinematics of these out-of-plane displacement peaks. The shape functions which span the region where the residual is high are refined, and the residual decreases (see Figure 8B). After the last refinement step, the residual has decreased to almost zero ( Figure 8D) in the entire region, indicating that the new set of shape functions is successful in describing the displacement field.
The resulting calculated displacement field is shown next to the analytical displacement field in Figure 9B. In case of a virtual experiment, it is possible to determine the exact error in the calculated displacement field, which is the difference between the exact and calculated displacement fields: ε w ¼ w ref ðxÞ À wðxÞ . The error field is displayed in Figure 11B. In the error field, small 'wiggles' appear that are characteristic for polynomial shape functions. However, looking at the values of the error field, compared to the amplitude of the sinusoidal peaks, the error is reasonably small. The method is therefore able to calculate sufficiently accurate results due to the autonomous refinements.
The adaptive iso-DHC method can readily be used with other polynomials orders of the NURBS shape functions. In the meshing procedure, the polynomial order is an input setting in AutoCAD, which can be set to the desired order. In the DHC algorithm, the Bézier functions for different orders Figure 6: The refinement indicator C j plotted for each shape function j (left). The position (top) of each shape function is shown in the right figure. The red line in the left image represents the threshold T ¼ max C þ σ; T nl À Á . All shape functions above this threshold are selected for refinement, as indicated by the red circle. In this case, this concerns only one shape function, namely shape function nr. 6, which corresponds to the refinement step between Figure 5A and 5B need to be implemented in order to use them. In this example, we repeat correlation of the above virtual experiment with firstand third-order shape functions. The same initial mesh ( Figure 7A) is used. Note that the number of DoFs is not the same for the three cases, since a set of higher order shape functions consists of more functions than a lower order basis. The refined meshes for the firstand third-order NURBS are shown in Figure 10. As can be seen, the refinement for the first-order shape functions remains more local than the second-order, while the thirdorder shape functions refine in a less local fashion. This is because NURBS shape functions overlap multiple elements, depending on their order. A first-order NURBS (not at the edge of the domain) covers two-by-two elements, while a third-order NURBS occupies four-by-four elements. As explained in "Adaptive refinement" Section, the entire support of the selected shape functions is refined, resulting in less local refinement for higher order shape functions. The resulting decrease in residual corresponding to the refined meshes of Figure 10 is similar to that of the second-order shape functions, shown in Figure 8, and therefore not shown.
The error fields resulting from the calculation of the displacement field with the adaptive iso-DHC method with firstand third-order NURBS are plotted in Figure 11A and C. Especially for the third-order shape functions, the characteristic 'wiggles' are again recovered; however, they now spread out over a larger region, originating from the larger support of the higher order functions. The level of the error is similar to the error of the second-order shape functions. The error field for the first-order shape functions exhibits more local features due to the more local nature of the lower order shape functions; however, its features have significantly higher amplitude, indicating that, despite refinement, these shape functions are not able to describe the displacement field as accurately as the higher order shape functions.
In this work, we demonstrate the method on stretchable interconnects, which have a slender geometry, thereby limiting the number of elements in the width direction. Therefore, the refinement process should be local. Firstorder shape functions were found to refine locally, but are not optimally suited for capturing the kinematics of localized buckling, while also providing only C 0 continuity on the element boundaries. Third-order shape functions were found to be less local. Therefore, second-order shape functions form an adequate compromise, with the preferred element boundary continuity of C 1 . They will be used for the remainder of the paper.

Noise analysis
The virtual experiment is repeated with shape functions of the second order for a case where noise is present. From real experimental data, the noise level is determined by subtracting multiple images taken subsequently, with no deformation, and calculating the RMS value of the residual.
The noise level is assessed at about 1.5%. A safety factor of 2 is administered, and random noise of 3% of the height level range is artificially applied to the images f and g. The same algorithm is applied, starting from the initial mesh that is shown in Figure 7A, and the resulting mesh refinement and residual field after refinement are shown in Figure 12.  From this image and the corresponding root mean square (RMS) value, it is observed that the residual decreases to about the level of the noise (which has an RMS value of 0.0964 μm), indicating that an optimal correlation has been obtained. The refined mesh is essentially the same as the case where no noise is present. The calculated displacement field and corresponding error field, including the RMS value, are similar to Figures 9B and 11B respectively and hence are not displayed here. The results illustrate the noise robustness of the proposed method.

Stretchable interconnect
Although the first virtual experiment provides a proof of concept for the adaptive iso-DHC algorithm, a more complex sample geometry and mesh are considered next. The method is applied to a virtual experiment on a stretchable interconnect (SI). To this end, a real height profile ( Figure 13A) from a stretchable electronics structure, measured using a Sensofar PLμ2300 confocal optical profilometer, is analytically deformed. The applied displacement field again represents localized buckles, as depicted in Figure 1D. The buckles are represented by two sinusoidal peaks that are cut off at the edge of the SI geometry. In this experiment, also in-plane deformation is considered, namely uniaxial stretching in x-direction and rigid body translation in y-direction. Like for the previous case, the final displacement field is applied in four increments, for which the resulting images in steps 2 and 4 are shown in Figure 13B and C.
The initial mesh is built with AutoDesk AutoCAD, as described in "Isogeometric shape functions and parametrization" Section. In this case, 6 × 2 elements are the minimum to accurately describe the sample's contour (see Figure 13A). The DHC algorithm is solved, and the mesh refines in the regions around the peaks (see Figure 13B and C).
From the residual images, shown in Figure 14, it is observed that with the initial mesh, it is not possible to accurately capture the kinematics of the buckles; therefore, the residual is high in the area surrounding the buckles. After refining the mesh, the residual decreases. Note that the residual does not decrease to the low level achieved in the virtual experiment of the previous section. This is due to the limited amount of pixels, only 768 × 558 used here  Notice that the RMS value of the error for the 1st order NURBS (0.04496 μm) is much larger than for the 2nd (0.01670 μm) and 3rd-order (0.01653 μm) NURBS compared to the image with a more common size of 1200 × 1000 pixels of the previous experiment. To more accurately describe the displacement field with localized buckles, the mesh should presumably be refined one more step. However, this is not allowed, since the number of pixels in an element would become too small, and the problem becomes more sensitive to noise. Therefore, this result is the best that can be obtained with this image.
A quantitative measure for the accuracy of the method is shown in Figure 15. The applied displacement fields uðxÞ , vðxÞ and wðxÞ, or reference fields, are shown in Figure 15A, and the same fields calculated with the DHC algorithm are depicted in Figure 15B. Since this is a virtual experiment, we are able to calculate the exact error field, which is the difference between the reference and calculated field. This error field is shown for all three directions in Figure 15C. It is observed that the displacement fields for both the inplane directions x and y and the out-of-plane direction z are captured well. They are calculated rather accurately, although the error is larger than in the previous test case. However, this is partly due to the lower amount of pixels in the image, which can be optimized by using a profilometer with a high-resolution camera. The novel method, with a complex mesh that was constructed using the proposed meshing procedure, is adequately able to autonomously determine the degrees of freedom that optimally describe the localized displacement field with a representative buckling profile.

Experiment: Application to Stretchable Electronics
A real experiment concerning a stretchable electronics interconnect is executed in order to demonstrate the adaptive isogeometric DHC method's performance in a real situation. If the interconnect is stretched, it buckles in certain regions. The objective of this experiment is to calculate the displacement field describing these buckles as well as the in-plane displacements, with an autonomously optimized set of shape functions.

Specimen and test set-up
The specimen used for this experiment is a stretchable electronics interconnect consisting of a 10-μm-thick polyimide substrate with a 1-μm-thick aluminium meander.   The residual images are shown after correlation using the corresponding meshes from Figure 13. It can be seen that the buckles cannot be described accurately with the initial mesh, resulting in a high residual in this area. After refinement in this region, the residual decreases significantly. The RMS values are reported below the figures This S-shaped aluminium interconnect structure has width 20 μm, inner radius 20 μm and a 40-μm rectilinear segment between the curvilinear sections (see Figure 16). For DHC, a certain contrast in height values on the sample, or pattern, is required. For this purpose, an Indium-Tin (In-Sn) layer is deposited using a planar magnetron sputtering system. In situ heating of the sample to 80°C, close to the melting temperature of In-Sn, in combination with a high deposition rate is used to prevent deposition of a homogeneous In-Sn layer, but instead achieve distinct height features. Since the temperature during this pattern deposition procedure is significantly lower than the processing temperature of the sample, it is not degraded using this technique.
The experimental set-up consists of a Kammarath&Weiss uniaxial tensile/compression module placed underneath a Sensofar Pl μ2300 confocal surface profilometer equipped with a 150X objective. After deposition of the indium-tin layer, the aluminium/polyimide interconnect is glued onto disposable grippers that are clamped in the tensile stage and stretched. After each elongation increment, a topographic image is acquired.

Results
The dimensions of the specimen are very small and therefore application of a pattern with sufficiently distinct features and accurate, reproducible measurement of this height profile with a profilometer is a known challenge, as discussed for instance in [8]. For instance, when comparing the profilometric images of different increments, it is observed that the pattern features change between the images (see Figure 17). This might be due to a relatively high signal-to-noise ratio caused by steep edges of the tiny pattern features, thereby working at the limits of the profilometer. Measurement of the features from a slightly different position and angle, due to in-plane deformation and especially out-of-plane rotation of the underlying sample, further decreases the measurement reproducibility.
These discrepancies between the incremental images make correlation difficult, because the residual will not reduce to (almost) zero for the correct displacement field. The detrimental effects are somewhat reduced by applying some Gaussian blurring (kernel size 10, σ = 2) to the images before correlation, which is a known technique to reduce bias error [34]. Still it was found that shape functions with a small support are sensitive to these measurement artefacts, especially the shape functions in the corners of the domain. To complicate the test further, it was observed that the correlation length of the pattern is small. This, in combination with the change in pattern features, makes any DIC algorithm sensitive to a good initial guess. It is a known feature of DIC algorithms that there is a possibility of correlation to a local minimum instead of the global minimum [2]. Starting from a coarse mesh with limited degrees of freedom reduces this risk, while using the correlation result of the coarse mesh as initial guess for correlation for a refined mesh. In the virtual cases, this was sufficient, and an initial guess of zero displacement everywhere was acceptable, but in this experiment a good initial guess is inevitable. This good initial guess for all images was obtained by a correction for rigid body motion of the specimen centre and running the algorithm first over all images with the refinement option turned off, i.e. with the large-area NURBS shape functions shown in Figure 18A.
For the correlation, we zoom in on one of the rectilinear parts, as this gives the most interesting displacement field, since it buckles upon stretching. The images before and after deformation are shown in Figure 18, with the initial and refined mesh (after two refinement steps) plotted on top. The buckles that occur are about 3.5 μm high and the mesh refines in the area of the buckles. However, the sensitivity of the corner shape functions to the measurement artefacts is clear in the refined mesh. The mesh is refined in the corner elements, while there is no kinematic reason for it. Refinement leads to more freedom in this area, causing an even higher sensitivity to experimental uncertainties and further refinement in the same area. Also note that in the second image, it might appear that the refined mesh does not correctly conform the sample anymore, i.e. the mesh appears smaller than Figure 16: Scanning electron microscopic image of the aluminium stretchable interconnect on a polyamide substrate, taken after deformation. The interconnect has delaminated from the substrate, which exposes regions of the substrate that are not covered with an Indium-Tin layer. The In-Sn layer is characterized by the granular texture on top of the entire sample. In the rectilinear parts, the lighter regions indicate the location of buckles the interconnect width; however, this illusion is caused by the delamination and out-of-plane rotation of the interconnect, causing the steep sides of the interconnect to rotate into view, thereby exposing new area that in the first image was not visible. This also becomes clear from the zoomed images in Figure 17, where the yellow encircled feature moves away from the 'edge' and another feature appears below it.
Since this is a real experiment, it is not possible to determine the accuracy of the correlation by means of error fields, in contrast to the virtual experiments. The accuracy therefore has to be determined using the residual fields. These fields before and after mesh refinement are shown in Figure 19. The buckles clearly show up in the first residual field, as regions with an averaged value that is systematically above (red) or below (blue) zero, which indicates that a finer mesh is desired in these areas to calculate the displacement field accurately. When the mesh refines, in two steps, the final residual field does not have regions with an averaged value systematically Figure 17: Images f and g zoomed in at the same area (blue boxes in insets). It is observed that the pattern features look distinctly different in the two images. For example, the feature in the circles is (almost) unrecognizable. This complicates correlation of the two images Figure 18: The reference image f (left) is shown along with the deformed image g (right). Both images are blurred. The initial mesh and final refined mesh are plotted on top different from zero anymore. Note that refinement in the corners does not improve the residual in this area significantly and occasionally even causes the residual to increase, as for example in the right top corner. This is because the increased number of degrees of freedom also increases the sensitivity to the measurement artefacts, which can lead to poorer correlation.
The out-of-plane displacement field ( Figure 20C) clearly shows the buckles observed in Figures 16 and 18B. The in-plane deformation represents mainly the rigid body rotation of the rectilinear part of the interconnect that occurs upon stretching, as is also observed in these figures. In all, the three-dimensional displacement field in Figure 20 seems to have been measured accurately.
To conclude, since the height profiles resulting from this experiment are difficult to process for any digital height correlation algorithm, a good initial guess was necessary to analyse the data with the adaptive isogeometric DHC method. However, then the method was able to provide accurate results, corresponding to the observed in-plane and out-of-plane displacements in the measured height profiles, and correlated to a relatively low residual field. The mesh was optimized autonomously to be able to describe the complex displacement field accurately, but refined in unnecessary areas due to sensitivity of the small corner shape functions to measurement artefacts. Still, autonomy of the algorithm was partly lost due to the necessity of the preconditioning. As the DHC algorithm seems to work correctly, further improvement should be obtained by application of a height pattern with larger correlation length and more robust measurement of surface height profiles.

Conclusion
A method has been developed which uses an adaptive refinement algorithm to nearly autonomously refine shape functions in a global digital height correlation technique. This method is useful in experiments where the kinematics of the deforming sample is not known in advance. The Figure 19: The residual images are shown after correlation using the corresponding meshes from Figure 18. It can be seen that the buckles cannot be described accurately with the initial mesh, resulting in a high residual in this area (RMS value 0.133 μm). After refinement in this region the residual decreases significantly (RMS value 0.062 μm) Figure 20: The calculated displacement fields uðxÞ in x-direction (a), vðxÞ in y-direction (b) and wðxÞ in z-direction (out-of-plane) (c) are shown on top of undeformed image f. Combining the two in-plane displacement fields yields a rigid body rotation, which can also be seen from the difference between Figures 18(a) and (b), as well as a slight elongation of the structure. In the out-of-plane displacement field, the observed buckles appear clearly, indicating an accurate calculation of the displacement field mesh autonomously adjusts to the displacement field, i.e. optimizing the set of shape functions for capturing the displacement field. An optimized number of shape functions, i.e. number of degrees of freedom, are beneficial in DHC problems. Sufficient DoFs are needed to capture the kinematics of the displacement field, but too many of them make the problem too sensitive to noise. Another advantage of this method is that only little user experience is required to construct a reliable discretization of the problem.
Non-uniform rational B-splines shape functions are used both for the discretization of the DHC problem and the parametrization of the sample geometry. NURBS were originally developed for CAD modelling, and in this work, the CAD representation of the sample is directly used for DHC. NURBS shape functions are geometrically rich and can describe many geometrical shapes and displacement fields. With the use of a CAD program for constructing the mesh, nearly any specimen geometry can be meshed. Moreover, one is not restricted to a particular polynomial order of the shape functions, as the order can simply be selected in the CAD program. In this work, it was chosen to use second-order NURBS, since for buckled samples, the curvature of the buckles is of particular interest and to this end second-order derivatives of the calculated displacement field are desired. A hierarchical approach has been implemented for the adaptive refinement of the shape functions. This way refinement is executed in an efficient way and stays local, in contrast to knot insertion where an entire column and row of shape functions are refined due to the tensor product structure. The adaptive refinement is based on the residual field.
A proof of concept of the novel method is given with a virtual experiment with an out-of-plane displacement field with two sinusoidal peaks, representing a localized buckle pattern. The algorithm works adequately and refines in the expected area, yielding an accurate result. The method is also applied to a virtual experiment where a more complex sample shape is used: a stretchable electronics interconnect. It was shown that the method succeeds in capturing both in-plane and out-of-plane deformation fields accurately and refines the mesh in the expected areas.
Finally the adaptive iso-DHC method is tested in an experimental setting, where an in situ tensile experiment is performed on a stretchable interconnect in a profilometer. This experiment formed a challenge, since the structure is of such small dimensions that reproducible measurement of the microscopic height pattern forms a significant challenge. Successful correlation of the resulting images would be difficult for any DIC algorithm. The problems were overcome by applying blurring of the images and supplying a good initial guess. This is a compromise, since the autonomy of the algorithm is decreased; however, the mesh is still autonomously refined in the necessary regions yielding accurate threedimensional deformation fields. Further improvement should be sought in more robust measurement of surface height profiles.
Concluding, the novel-adaptive isogeometric DHC algorithm is a versatile technique for correlating displacement fields using the height profiles of many different types of experiment, including samples of which the deformation mechanism is unknown in advance (e.g. stretchable electronics interconnects). The shape functions used for discretizing the displacement field adjust autonomously to enrich the kinematics in regions where this is needed, obviating the need for detailed mechanical knowledge in advance. Accurate results have been obtained with the method, making this a promising technique for experimental mechanics of solids.