1 Introduction
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 tensorproduct 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., Tsplines, hierarchical splines, PHTsplines, and LR splines.
The localised splines of the first type, Tsplines, were introduced in [66, 65] and analysed in [1, 4, 63, 64]. They are based on the Tjunctions that allow eliminating redundant control points from NURBS model. The thorough study confirmed that this approach generates an efficient local refinement algorithm for analysissuitable Tsplines [36] and avoids the excessive propagation of control points. In [3, 9], it was proposed to characterise such splines as dualcompatible Tsplines, and in [43] 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 Bsplines (HBsplines), 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 Bsplines (THBsplines) have been developed (see [21]). In addition to good stability and approximation properties inherited from HBsplines [19, 67], THBsplines form a convex partition of unity. Therefore, they are suitable for the application in CAD. Various usage of THBspline 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 Tmeshes, 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 [11] or fulfil a certain constraint on admissible mesh configuration [74]. In [47], the adaptive procedure in based on the recoverybased error estimator, in which discontinuous enrichment functions are added to the IgA approximation.
Finally, locally refined splines (LRsplines) rely on the idea of splitting basis functions. This technique achieves localisation but creates difficulties with linear independence [13], which has been studied in [6, 7]. The application of such type of splines has been thoroughly investigated in [13]. In [26], one can find the summary of a detailed comparison of (T)HBsplines and LR splines with respect to sparsity and condition numbers. The study concludes that even though LR splines have smaller support than THBsplines, 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 residualbased 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 Clementtype interpolation operators, which are meshdependent and often difficult to compute for general element shapes. Moreover, these constants must be reevaluated every time a new mesh is generated. The goaloriented 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 tensorproduct splines is presented in [28] 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 Bsplines for a rather efficient calculation of true error upper bound.
The current work further extends the ideas used in [28] for Bsplines (NURBS) and combines the functional approach to the error control with THBsplines. 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 tensorstructured Bsplines (NURBS) and THBsplines, we manage to obtain an impressive speedup in majorant reconstruction by exploiting high smoothness of Bsplines to our advantage. However, for the problems with sharp local changes or various singularities in the solution, the THBsplines implementation in G+Smo restricts the performance speedup when it comes to solving the optimal system for the error majorant as well as for its elementwise 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 multipatch domain, since the error estimates for stationary problems are flexible enough to handle fully nonconforming 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 nonoverlapping convex subdomains, 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 (IBVP). In order to use this method, one needs to impose a crucial restriction on the multipatch configuration, namely, each patch must be a convex subdomain. 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 errordistribution indication. The next section serves as an overview of IgA techniques used in the current work, i.e., Bsplines, NURBS, and THBsplines. 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 wellposedness 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
(1)  
(2) 
where supposed to be in . Alternatively, the problem (1)–(2) can be viewed a system with two (primal and dual) unknowns
(3)  
(4)  
(5) 
We assume that the operator is symmetric and satisfies the condition of uniform ellipticity for almost all (a.a.) , which reads
(6) 
with . Throughout the paper, the following notation for the norms is used:
where stands for a weighted scalarproduct for all . After multiplying (1) by the test function
we arrive at the standard generalised formulation of (1)–(2): find satisfying the integral identity
(7) 
According to [33], the generalised problem (7) has a unique solution in provided that and condition (6) holds.
We consider the functional error estimates, which provide a guaranteed twosided 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
(8) 
as well as its elementwise 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.
Remark 1
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 righthand side (RHS) and by setting . Thus, one obtains the error identity
(9) 
The main idea of functional approach is the introduction of an auxiliary vectorvalued variable
satisfying
(10) 
In further calculations, the aboveintroduced 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(11) 
The equilibrated and dual residualfunctionals obtained in the RHS of (11) mimic equations (3) and (4), respectively, and are denoted by
(12) 
Theorem 1
(a) For any functions and , we have the estimate
(13) 
where the residuals and are defined in (12), is a positive parameter, and is the constant in the Friedrichs inequality [18]
(b) For , the variational problem
has a solution (with the corresponding zerovalue 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 reallife 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).
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.
Remark 2
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 leastsquare 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 .
Remark 3
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 gradientaveraging 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 postprocessing 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.,
(15) 
From (15), it is easy to see that if the auxiliary is chosen optimally and is set to zero (in the RHS of (15)), there is no gap between and .
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
(16) 
The latter one is equivalent to the variational formulation for the optimal , i.e.,
where the optimal is given by with
(17) 
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 [20])
(18) 
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: Bsplines, NURBS, and THBsplines
For the consistency of exposition, we first give an overview of the general IgA framework, the definitions of Bsplines, NURBS, and THBsplines, 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 knotvector in is a nondecreasing 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 onedimensional parametric domain , denotes a locally quasiuniform mesh, where each element is constructed by distinct neighbouring knots. The global size of is denoted by
Henceforth, we assume locally quasiuniform meshes, i.e., the ratio of two neighbouring elements and satisfies the inequality
The univariate Bspline basis functions are defined by means of the Coxde Boor recursion formula
(19) 
where a division by zero is defined to be zero. The Bsplines are times continuously differentiable across the th knot with multiplicity . Hence, if for inner knots, the Bsplines of the degree e.o.c. are continuous across them.
The multivariate Bsplines on the parameter domain , , are defined as tensor products of the corresponding univariate ones. In the multidimensional case, we define a knotvector dependent on the coordinate direction , , where indicates the direction (in space or time). Furthermore, we introduce a set of multiindices
and a multiindex indicating the order of polynomials. The tensorproduct of univariate Bspline basis functions generates multivariate Bspline basis functions
(20) 
The univariate and multivariate NURBS basis functions are defined in a parametric domain by means of Bspline basis functions, i.e., for a given and any , NURBS basis functions are defined as
(21) 
where . To recall basic definitions related to THBsplines, we follow the structure outlined in [20] and consider a finite sequence of nested variate tensorproduct spline spaces defined on the axis aligned boxdomain . To each space we assign a tensorproduct Bspline basis of degree
where is a set of multiindices for each level, and denotes the number of univariate Bspline 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 columnvector of basis functions. Then, a spline function is defined by and a coefficient matrix , i.e.,
where are rowcoefficients of .
Since , the basis can be represented by the linear combination of , namely,
where is a refinement matrix. Its entries can be obtained from Bsplines refinement rules (see [49]). Along with nested space, a corresponding sequence of nested domains is considered
(22) 
where each is covered by a collection of cells with respect to the tensorproduct 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 an
index setThe initial hierarchical data structure defined by the tensorproduct 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 THBspline basis related to the hierarchical domains is defined as
where the truncation of any function w.r.t. level is defined by
Here,
denotes an identity matrix
of 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 HBsplines for the hierarchical levels , whereas Figure (b)b shows the same levels for THBsplines. On the last raw, basis functions influenced by the truncation on coarser levels are exposed, i.e., THBsplines 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 :
(23) 
where are control points, and stands for either Bsplines, NURBS, or THBbasis functions . The mesh discretising consists of elements that are the images of , i.e.,
The global meshsize is denoted by
(24) 
Moreover, we assume that is a quasiuniform 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 functionaltype 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
In order to keep the presentation concise, we restrict (3)–(5) to the Dirichlet–Poisson problem
(25) 
Let the approximation
Here, is generated with NURBS of degree , i.e., . Due to the onepatch 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 THBsplines. Therefore, the constructed approximation can be written as
where is a vector of degrees of freedom (d.o.f.) defined by the system
(26) 
The majorant corresponding to the problem (25) reads as
(27) 
where and are defined in (17), and . Here, the approximation space for
is generated by the pushforward 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 [28]. 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
(28) 
where
According to the numerical results obtained in [28], 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
(29) 
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 scalarfunctions generated by splines  
approximation space for the dimensional vectorfunctions generated by splines  
approximation space for the twodimensional vectorfunctions generated by splines  
coarsening ratio of the global size of the mesh for approximation to the global size of the mesh for 
Comments
There are no comments yet.