The investigation of effective adaptive refinement procedures has recently become an active area of research in the context of fast and efficient solvers for isogeometric analysis (IgA) [24, 25]. The adaptivity scheme is naturally linked with reliable and quantitatively efficient a posteriori error estimation tools. The latter ones are expected to identify the parts of a considered computational domain with relatively high discretisation errors and provide a fully automated refinement strategy in order to reach desired accuracy levels for an approximate solution.
Due to a tensor-product setting of IgA splines, mesh refinement has global effects, which include a large percentage of superfluous control points in data analysis, unwanted ripples on the surface, etc. These issues produce certain challenges at the design stage as well as complications in handling big amounts of data, and, therefore, naturally trigger the development of local refinement strategies for IgA. At the moment, four different IgA approaches for adaptive mesh refinement are known, i.e., T-splines, hierarchical splines, PHT-splines, and LR splines.
The localised splines of the first type, T-splines, were introduced in [66, 65] and analysed in [1, 4, 63, 64]. They are based on the T-junctions that allow eliminating redundant control points from NURBS model. The thorough study confirmed that this approach generates an efficient local refinement algorithm for analysis-suitable T-splines  and avoids the excessive propagation of control points. In [3, 9], it was proposed to characterise such splines as dual-compatible T-splines, and in  a refinement strategy with linear complexity was described for the bivariate case.
The alternative approach that implies the local control of refinement is based on hierarchical B-splines (HB-splines), such that a selected refinement region basis functions are replaced with the finer ones of the same type. The procedure of designing a basis for the hierarchical spline space was suggested in [17, 29, 23] and extended in [70, 16, 62]. Such construction guarantees the linear independence of the basis and provides nested approximation spaces. However, since the partition of unity is not preserved for these splines, truncated hierarchical B-splines (THB-splines) have been developed (see ). In addition to good stability and approximation properties inherited from HB-splines [19, 67], THB-splines form a convex partition of unity. Therefore, they are suitable for the application in CAD. Various usage of THB-spline for arbitrary topologies can be found, e.g., in [72, 76, 77].
The locally defined splines of the third type, namely, polynomial splines over hierarchical T-meshes, are constructed for the entire space of piecewise polynomials with given smoothness on the subdivision of considered domain. The corresponding application can be found in [46, 71]. However, in this case, one must assume the reduced regularity of basis  or fulfil a certain constraint on admissible mesh configuration . In , the adaptive procedure in based on the recovery-based error estimator, in which discontinuous enrichment functions are added to the IgA approximation.
Finally, locally refined splines (LR-splines) rely on the idea of splitting basis functions. This technique achieves localisation but creates difficulties with linear independence , which has been studied in [6, 7]. The application of such type of splines has been thoroughly investigated in . In , one can find the summary of a detailed comparison of (T)HB-splines and LR splines with respect to sparsity and condition numbers. The study concludes that even though LR splines have smaller support than THB-splines, the numerical experiments did not reveal any significant advantages of the first ones with respect to the sparsity patterns or condition numbers of mass and stiffness matrices.
The refinement tools of IgA mentioned above were combined with various a posteriori error estimation techniques. For instance, the a posteriori error estimates based on hierarchical splines were investigated in [14, 70]. In [27, 71, 8, 30]
, the authors used the residual-based a posteriori error estimates and their modifications in order to construct mesh refinement algorithms. The latter ones, in particular, require the computation of constants related to the Clement-type interpolation operators, which are mesh-dependent and often difficult to compute for general element shapes. Moreover, these constants must be re-evaluated every time a new mesh is generated. The goal-oriented error estimators, which are rather naturally adapted to practical applications, have lately been introduced for IgA approximations and can be found in[69, 10, 31, 32].
In the current work, the terms error estimate and error indicator distinguish from each other. The first one is considered as the total upper (or lower) bound of true energy error. These are very important characteristics related to the approximate solution since they can be used to judge whether obtained data are reliable or not. In order to locate the areas of the discretised domain that have the highest error in the approximation, a quantitively sharp error indicator
is required. The methods of a posteriori error estimation listed above are rather error indicators in this terminology and indeed were successfully used for mimicking the approximation error distribution. However, their use in the error control, i.e., a reliable estimation of the accuracy of obtained data, is rather heuristic in nature.
Below we investigate a different functional method providing fully guaranteed error estimates, the upper (and lower) bounds of the exact error in the various weighted norms equivalent to the global energy norm. These estimates include only global constants (independent of the mesh characteristic ) and are valid for any approximation from admissible functional space. One of the most advantageous properties of functional error estimates is their independence of the numerical method used for calculating approximate solutions. The strongest assumption about approximations is that they are conforming in the sense that they belong to a certain natural Sobolev space suited for the problem. It is important to emphasise that this is still a rather weak assumption and that no further restrictions, such as Galerkin orthogonality, are needed.
Functional error estimates were initially introduced in [59, 60] and later applied to different mathematical models summarised in monographs [45, 52, 37]. They provide guaranteed, sharp, and fully computable upper and lower bounds of errors. A pioneering study on the combination of functional type error estimates with the IgA approximations generated by tensor-product splines is presented in  for elliptic boundary value problems (BVP). The extensive numerical tests presented in this work confirmed that majorant produces not only good upper bounds of the error but also a quantitatively sharp error indicator. Moreover, the authors suggest the heuristic algorithm that allows using the smoothness of B-splines for a rather efficient calculation of true error upper bound.
The current work further extends the ideas used in  for B-splines (NURBS) and combines the functional approach to the error control with THB-splines. Moreover, our focus is concentrated not only on the qualitative and quantitative performance of error estimates but also on the required computation time for their reconstruction. The systematic analysis of majorant’s numerical properties is based on a collection of extensive tests performed on the problems of different complexity. For the error control implemented with the help of tensor-structured B-splines (NURBS) and THB-splines, we manage to obtain an impressive speed-up in majorant reconstruction by exploiting high smoothness of B-splines to our advantage. However, for the problems with sharp local changes or various singularities in the solution, the THB-splines implementation in G+Smo restricts the performance speed-up when it comes to solving the optimal system for the error majorant as well as for its element-wise evaluation. We restrict this study only to the domains modelled by a single patch, which provides at least -continuity of the approximate solutions inside the patch. However, the application of studied majorants can be extended to a multi-patch domain, since the error estimates for stationary problems are flexible enough to handle fully non-conforming approximations (this issue has been in details addressed in [35, 68, 61]).
The error control for the problems defined on domains of complicated shapes induces another issue related to the estimation of Friedrichs’ constant used by functional error estimates not only as the weight but also as the geometrical characteristics of the considered problem. When such domains are concerned, one can perform their decomposition into a collection of non-overlapping convex sub-domains, such that the global constant can be replaced by constants in local embedding inequalities (Poincaré and Poincaré-type inequalities [50, 51]). The reliable estimates of these local constants can be found in [48, 2, 44, 42]. The derivation of functional error estimates exploiting these ideas is discussed in [52, 54] for the elliptic BVP and in [40, 39, 41] for the parabolic initial boundary value problem (I-BVP). In order to use this method, one needs to impose a crucial restriction on the multi-patch configuration, namely, each patch must be a convex sub-domain. Since in the IgA framework patches are treated as mappings from the reference domain , the estimation of local constants is reduced to the analysis of the IgA mapping and calculating the corresponding constant for .
The paper proceeds with the following structure. Section 2 formulates the general statement of the considered problem and recalls the definition of functional error estimates and their main properties in the context of reliable energy error estimation and efficient error-distribution indication. The next section serves as an overview of IgA techniques used in the current work, i.e., B-splines, NURBS, and THB-splines. In Section 4, we focus on the algorithms and details of the functional error estimates integration into the IgA framework. Last but not least, Section 5 presents the systematic selection of most relevant numerical examples and obtained results that illustrate numerical properties of studied error estimates and indicators.
2 Functional approach to the error control
In this section, we present a model problem, recall the well-posedness results for linear parabolic PDEs, which have been thoroughly studied in [33, 75, 73]. We also introduce functional a posteriori error estimates for the stated model and discuss its crucial properties.
Let , , be a bounded domain with Lipschitz boundary . The elliptic BVP is formulated as
We assume that the operator is symmetric and satisfies the condition of uniform ellipticity for almost all (a.a.) , which reads
with . Throughout the paper, the following notation for the norms is used:
where stands for a weighted scalar-product for all . After multiplying (1) by the test function
We consider the functional error estimates, which provide a guaranteed two-sided bound of the distance between the generalised solution of (7) and any function . It is important to emphasise that the suggested functional approach to error estimates’ derivation is universal for any numerical method used to discretise bilinear form (7). This fact makes it rather unique in comparison with alternative approaches, which are always tailored to the discretised version of the identity . Later on, the considered is generated numerically, and the distance to is evaluated in terms of the total energy norm
as well as its element-wise contributions such that
Here, represents the elements of the mesh introduced on . Hence, besides providing the guaranteed upper bound of total error (8), the majorant yields a quantitatively sharp indicator of the local error distribution.
It has been proved that the integrand of the majorant tends to the distribution of the true error in the sense of the measure, if the majorant globally converges (from above) to the true value of the error. In other words, the measure of the set, where the local difference between majorant and true error is bigger than a given epsilon, tends to zero. For the purposes of the error indication, it is enough ([52, Section 3]).
To derive the upper bound, we need to transform (7) by subtracting from left- (LHS) and right-hand side (RHS) and by setting . Thus, one obtains the error identity
The main idea of functional approach is the introduction of an auxiliary vector-valued variable
In further calculations, the above-introduced variable allows an additional degree of freedom for the majorant (additional optimisation step), whereas, for instance, the residual error estimates do not have this flexibility in improving its values. Next, we add the identity (10) to the RHS of (9), which yields
(a) For any functions and , we have the estimate
(b) For , the variational problem
has a solution (with the corresponding zero-value for the functional), and its minimum is attained if and only if and .
Proof: For the detailed proof of this theorem we refer the reader to [53, Section 3.2].
For the real-life problems, the exact solution is usually not known, therefore evaluation of becomes impossible. For the efficiency verification of , the lower bound of the error (further also referred as minorant) has been derived using variational arguments (see Theorem 2).
Proof: For the detailed proof of this theorem, we refer the reader to [53, Section 4.1].
Remarks below summarise several essential properties of the error estimate derived in Theorem 1.
Each term on the RHS of (13) serves as the penalty of the error that might occur in equations (3) and (4). The positive weight can be selected optimally in order to get the best value of the majorant. The constant acts as a geometric characteristic for the considered domain (unlike, for instance, in the least-square methods, where the weights are selected after the terms and were calculated). From the author’s point of view, this constant is essential and cannot be excluded since it scales proportionally to the diameter of the considered . Moreover, in order to guarantee the reliability of , the constant must be estimated from above in a reliable way. Since in practice the term is rather small compared to the dominating term , the Friedrichs constant can be replaced by some penalty constant (even though it might affect the ratio of the majorant to the error). In what follows, to characterise the efficiency of (13), we use the quantity that measures the gap between and .
The functional generates the upper bound of the error for any auxiliary and , therefore the choice of might vary. The first and most straightforward way to select this variable is to set , where is a certain gradient-averaging operator. For this case, using the IgA framework becomes quite advantageous since for splines of the degree an obtained is a -continuous function and is already in . Therefore, no additional post-processing is needed. On the other hand, due to the quadratic structure of the majorant, it is rather obvious that the optimal error estimate value is achieved at , i.e.,
One of the numerical methods providing an efficient reconstruction of both dual and primal variables is a mixed method. It generates an efficient approximation of the pair that can be straightforwardly substituted into the majorant . Moreover, if the error is measured in terms of the combined norm, i.e., including the norm of the error in primal and in dual variables
it is controlled by the residuals of the majorant as follows:
We note that the ratio between the majorant (that does not include any constants) and the error is controlled by , which proves the robustness of such error estimate. The series of works (see, e.g., [56, 58, 57]) has confirmed the efficiency of combination of mixed methods and functional error estimates.
The alternative approach providing an accurate -reconstruction follows from the minimisation problem
The latter one is equivalent to the variational formulation for the optimal , i.e.,
where the optimal is given by with
In this work, using the IgA approximation schemes’ setting, we apply the second method of the efficient -reconstruction described in detail in Section 4. To compare the performance of with alternative error estimates we use the standard residual error estimator (applied, e.g., in )
where denotes the diameter of cell , and stands for approximation reconstructed by the IgA scheme. The term measuring the jumps across the element edges, which is usually included into residual error estimates, vanishes in (18) due to the properties of produced by the IgA schemes. It is provided in the G+Smo package and can be accessed by using the class available from the G+Smo library [38, stable/src/gsErrEstPoissonResidual.h].
3 IgA overview: B-splines, NURBS, and THB-splines
For the consistency of exposition, we first give an overview of the general IgA framework, the definitions of B-splines, NURBS, and THB-splines, their use in the geometrical representation of the computational domain and in the construction of IgA discretisation spaces.
Let denote the degree of polynomials used for the IgA approximations, and let be the number of basis functions used to construct a -spline curve. The knot-vector in is a non-decreasing set of coordinates in the parameter domain, written as , , where and . The knots can be repeated, and the multiplicity of the -th knot is indicated by . Throughout the paper, we consider only open knot vectors, i.e., . For the one-dimensional parametric domain , denotes a locally quasi-uniform mesh, where each element is constructed by distinct neighbouring knots. The global size of is denoted by
Henceforth, we assume locally quasi-uniform meshes, i.e., the ratio of two neighbouring elements and satisfies the inequality
The univariate B-spline basis functions are defined by means of the Cox-de Boor recursion formula
where a division by zero is defined to be zero. The B-splines are -times continuously differentiable across the -th knot with multiplicity . Hence, if for inner knots, the B-splines of the degree e.o.c. are continuous across them.
The multivariate B-splines on the parameter domain , , are defined as tensor products of the corresponding univariate ones. In the multidimensional case, we define a knot-vector dependent on the coordinate direction , , where indicates the direction (in space or time). Furthermore, we introduce a set of multi-indices
and a multi-index indicating the order of polynomials. The tensor-product of univariate B-spline basis functions generates multivariate B-spline basis functions
The univariate and multivariate NURBS basis functions are defined in a parametric domain by means of B-spline basis functions, i.e., for a given and any , NURBS basis functions are defined as
where . To recall basic definitions related to THB-splines, we follow the structure outlined in  and consider a finite sequence of nested -variate tensor-product spline spaces defined on the axis aligned box-domain . To each space we assign a tensor-product B-spline basis of degree
where is a set of multi-indices for each level, and denotes the number of univariate B-spline basis functions in the -th coordinate direction. After assuming that has a fixed ordering and rewriting the basis as it can be considered as a column-vector of basis functions. Then, a spline function is defined by and a coefficient matrix , i.e.,
where are row-coefficients of .
Since , the basis can be represented by the linear combination of , namely,
where is a refinement matrix. Its entries can be obtained from B-splines refinement rules (see ). Along with nested space, a corresponding sequence of nested domains is considered
where each is covered by a collection of cells with respect to the tensor-product grid of level . In this work, we focus on dyadic cell refinement for the bi- and trivariate cases with uniform degrees for all levels and coordinate directions, therefore, in further exposition.
Let the characteristic matrix of w.r.t. domains and is defined as
Next, for each level , the set of the indices of active functions can be defined with
To store the indices of all active functions at all hierarchical levels, we define anindex set
The initial hierarchical data structure defined by the tensor-product mesh (see Figure (a)a). In particular, we illustrate the knot lines of the spaces (where levels increase form the left to the right). By means of the insertion operation, new subdomains can be added to obtain new representations . The sets of active basis functions as well as the characteristic matrices for all levels are extracted simultaneously with the new box insertion and initialisation of the basis. Figure illustrates meshes at the refinement levels .
The THB-spline basis related to the hierarchical domains is defined as
where the truncation of any function w.r.t. level is defined by
denotes an identity matrixof size the multiplication of by represents w.r.t. to the level , and additional multiplication by performs the truncation operation. For the detailed discussion of truncation operation, we refer the reader to [21, 22, 20]. We illustrate an effect of the truncation for the case of univariate quadratic spline basis functions (see Figure 8). Figure (a)a presents the HB-splines for the hierarchical levels , whereas Figure (b)b shows the same levels for THB-splines. On the last raw, basis functions influenced by the truncation on coarser levels are exposed, i.e., THB-splines before and after being truncated are denoted by the grey and black marker, respectively.
The physical domain is defined by the geometrical mapping of the parametric domain :
where are control points, and stands for either B-splines, NURBS, or THB-basis functions . The mesh discretising consists of elements that are the images of , i.e.,
The global mesh-size is denoted by
Moreover, we assume that is a quasi-uniform mesh, i.e., there exists a positive constant independent of , such that
4 Functional error estimates within the IgA framework
In this section, we present the algorithms used for general reliable computations and functional-type error estimates reconstruction. Then we proceed with commenting on the implementation of these error estimates in G+Smo and their integration into the library’s structure. Finally, we present a series of examples demonstrating numerical properties of derived error majorants.
4.1 Reliable reconstruction of IgA approximations. Algorithms
Let the approximation
Here, is generated with NURBS of degree , i.e., . Due to the one-patch setting and restriction on the knots’ multiplicity of , the smoothness is automatically provided. Since no numerical algorithms specific to the hierarchical levels of the localised splines will be discussed below, we use the same notation for spaces generated by THB-splines. Therefore, the constructed approximation can be written as
where is a vector of degrees of freedom (d.o.f.) defined by the system
The majorant corresponding to the problem (25) reads as
where and are defined in (17), and . Here, the approximation space for
is generated by the push-forward of a corresponding space in the parametric domain
Here, is a space of NURBS with the degree for each of components of The details of the numerical reconstruction of (27) were thoroughly studied in . The best estimate follows from the optimisation of w.r.t. function
The basis functions generate the space , whereas is a vector of d.o.f of defined by a system
According to the numerical results obtained in , the most efficient majorant reconstruction (with uniform refinement) is obtained when is set substantially higher than . Let us assume that , . At the same time, when is reconstructed on the mesh , we use a coarser one , in order to recover . For the reader’s convenience, all used notation is summarised in Table 1. The initial mesh and the basis functions defined on it are assumed to be given via the geometry representation of the computational domain. The exact representation of geometry on the initial (the coarsest) level is preserved in the process of mesh refinement.
For the reconstruction of , let the approximation
Here, is approximation space generated with NURBS of degree on the parameter domain, i.e., . Then, the auxiliary approximation can be written as
where is a vector of degrees of freedom (d.o.f.) defined by the system
Analogously to the selection of the for the space , we let , . At the same time, we use a coarser mesh , for the approximation.
|degree of the splines used for approximation|
|degree of the splines used for approximation|
|degree of the splines used for approximation|
|()||approximation space for the scalar-functions generated by splines|
|approximation space for the -dimensional vector-functions generated by splines|
|approximation space for the two-dimensional vector-functions generated by splines|
|coarsening ratio of the global size of the mesh for approximation to the global size of the mesh for|