Residual-based error estimation and adaptivity for stabilized immersed isogeometric analysis using truncated hierarchical B-splines

We propose an adaptive mesh refinement strategy for immersed isogeometric analysis, with application to steady heat conduction and viscous flow problems. The proposed strategy is based on residual-based error estimation, which has been tailored to the immersed setting by the incorporation of appropriately scaled stabilization and boundary terms. Element-wise error indicators are elaborated for the Laplace and Stokes problems, and a THB-spline-based local mesh refinement strategy is proposed. The error estimation .and adaptivity procedure is applied to a series of benchmark problems, demonstrating the suitability of the technique for a range of smooth and non-smooth problems. The adaptivity strategy is also integrated in a scan-based analysis workflow, capable of generating reliable, error-controlled, results from scan data, without the need for extensive user interactions or interventions.

For mixed formulations, such as standard weak forms of the Stokes and Navier-Stokes equations, the immersed isogeometric analysis setting imposes an additional challenge.In order to satisfy the inf-sup condition [41,42] in boundary-fitting (isogeometric) analyses, generally use is made of stable pairs of basis functions (e.g., Taylor-Hood [43][44][45][46] or Raviart-Thomas [45,[47][48][49]).Alternatively, stabilization techniques such as GLS [50][51][52], VMS [53][54][55] or projection methods [56,57] can be used.Direct utilization of these elements or stabilization techniques in the immersed setting can lead to non-physical spurious oscillations in the solution, even with relatively large and regular cut element configurations [58,59].One remedy for tackling this issue is to employ a skeleton-stabilized immersed isogeometric technique [59].The fundamental idea of this stabilization technique is to penalize (high-order) pressure derivative jumps over the edges/faces of the background mesh, resulting in stable discretizations using equal-order spline spaces.The technique proposed in Ref. [59] is inspired by the (continuous) interior penalty ((C)IP) and the ghost penalty (GP) methods [27], extending these techniques to the case of high-regularity isogeometric analysis.
An appraised property of immersed methods in general, and immersed isogeometric analysis in particular, is that the discretization resolution can be controlled independently of the geometry parametrization.The immersed analysis concept avoids the need for geometry-induced mesh refinements in the vicinity of geometric details that are irrelevant in relation to the objective of an analysis.This decoupling of the discretization resolution from the geometry makes it natural to consider immersed finite elements in combination with adaptive discretization strategies.In fact, adaptivity in the form of local p-and hp−refinements has always been an integral part of the finite cell method [60][61][62][63].
A posteriori error estimation and adaptivity techniques are well-established in the context of finite element methods; see, e.g., the reviews [64][65][66].A variety of error estimation and adaptivity techniques has been studied in isogeometric analysis, such as residual-based error estimators for T-splines [67] and hierarchical splines [68][69][70], and goal-oriented techniques [71].The contemporary overview [72] is also noteworthy, as is the advanced industrial application considered in Ref. [73].In the context of Nitsche-based finite element methods (see Refs. [74,75] for an overview), studies on a posteriori error estimators have been conducted [4,[76][77][78][79].Local refinement strategies in immersed methods are predominantly feature based, i.e., either based on geometric features such as boundaries, or based on solution features such as sharp gradients in the solution fields; see, e.g., [11,80,81] for examples of local refinement capabilities in finite cell simulations.Goal-oriented error estimation and adaptivity for immersed methods has also been studied [71,[82][83][84].In the context of stabilized immersed finite elements, Ref. [85] considered a posteriori element-wise error estimation and adaptivity to improve boundary approximations.
Although the computational setting of immersed isogeometric analysis enables the use of volumetric spline patches, the standard h, p and k-type refinement strategies in patch-based isogeometric analysis [7] are not suitable because of the non-local propagation of refinements.Various alternative refinement strategies have been proposed over the last decade to construct local spline refinements, the most prominent of which are T-splines [11,16,[86][87][88][89][90], LRB-splines [91,92], U-splines [93], and (Truncated) Hierarchical B-splines [94].
In the context of immersed isogeometric analysis on volumetric domains, hierarchical splines are particularly suitable, as they optimally leverage the advantages offered by the geometrically simple background mesh.
In this contribution we propose a computational strategy for the application of residual-based a-posteriori error estimation and mesh adaptivity to stabilized immersed isogeometric analyses.We study various computational aspects of the framework that are non-standard in comparison to error estimation and adaptivity for boundary-fitting analyses, viz.: (i) In immersed analyses, the discretization basis is constructed over a mesh comprised of all elements in an ambient mesh that intersect with the computational domain.As a direct consequence of this setting, the support of the computational basis in general changes under refinement operations.The same holds for the mesh skeleton, which is a key ingredient of the considered stabilization methods.The considered computational strategy preserves the geometry of the computational domain under local mesh refinements, despite the change of the background mesh; (ii) Weak formulations in stabilized immersed isogeometric analysis generally involve operators with an explicit dependence on the mesh size.While this mesh size is unambiguously defined in the case of a uniform background mesh, the local mesh refinements considered in the adaptive setting warrant careful consideration of the scaling of the stabilization terms.We herein propose and study a scaling of the stabilization terms based on the local element sizes.
We demonstrate the performance of the proposed computational strategy using a series of test cases for steady heat conduction problems (Poisson problem) and steady viscous flow problems (Stokes problem).We consider the application of the proposed adaptivity technique in a scan-based isogeometric analysis setting, and demonstrate that a robust automatic simulation workflow is realized when the methodology presented herein is combined with the topology-preserving image segmentation algorithm presented in Ref. [22].
This paper is outlined as follows.Section 2 introduces the immersed isogeometric analysis framework, along with a detailed stability analysis for the considered model problems.This analysis focuses particularly on the scaling relations for the stabilization terms.In Section 3 the residual-based error estimator is introduced, and a mesh-refinement strategy is proposed.Benchmark simulation results are then presented in Section 4 for both the steady heat conduction problem and the viscous flow problem, after which the developed framework is applied in a scan-based setting in Section 5. Conclusions are finally drawn in Section 6.

Stabilized immersogeometric analysis with local mesh refinements
In this section we introduce the stabilized immersed isogeometric analysis formulations for the steady heat conduction (Laplace) problem and steady viscous flow (Stokes) problem.We commence with presenting the general setting of the problems in Section 2.1, after which the stabilized formulations are presented in Section 2.2.In preparation of the a posteriori error estimation concept discussed in Section 3, in Section 2.3 we study the stability of the considered formulations.

The finite cell setting
We consider a physical domain Ω ∈ R d (with d ∈ {2, 3}) with boundary ∂Ω, as illustrated in Figure 1.The boundary is composed of a Neumann part, ∂Ω N , and a Dirichlet part, ∂Ω D , such that ∂Ω N ∪∂Ω D = ∂Ω and ∂Ω N ∩ ∂Ω D = ∅.The outward-pointing unit normal vector to the boundary is denoted by n.
The physical domain is immersed in a geometrically simple ambient domain, i.e., A ⊃ Ω, on which a locally refined ambient mesh T A with elements K is defined.In this work, the ambient domain is chosen to be rectangular or cuboid, to facilitate simple, tensor-product, spline discretizations.The locally-refined meshes are constructed by sequential bisectioning of (selections of) elements in the mesh, starting from a Cartesian mesh.Truncated hierarchical B-splines can be formed on such meshes, as will be elaborated in Section 2.2.
Figure 1: (a) A physical domain Ω, with boundary ∂Ω, is embedded in an ambient domain A. The background mesh T , which consists of all elements that intersect the physical domain, is constructed by locally refining the ambient domain mesh T A .The zoom illustrates the employed bisectioning procedure to capture the immersed boundaries.The integration subcells are marked in blue, whereas the background cells are marked in black.The skeleton mesh, F skeleton , and ghost mesh, F ghost , are shown in panels (b) and (c), respectively.
Elements that do not intersect with the physical domain can be omitted from the ambient mesh, resulting in the locally refined (active) background mesh In the remainder, with the abuse of notation, we will use T (and other meshes) to denote both the set of elements in the mesh and the geometry obtained from the union of these elements.The local mesh size of the locally refined background mesh is denoted by By cutting the elements that are intersected by the immersed boundary ∂Ω, a mesh that conforms to the physical domain Ω is obtained: The collection of elements in the background mesh that are crossed by the immersed boundary ∂Ω is defined as In immersed methods, the geometry of the physical domain is captured by the integration procedure on the cut elements, i.e., elements that are intersected by the immersed boundary ∂Ω.We herein employ an octree integration procedure [22,82], which we close at the lowest level of bisectioning with a tessellation procedure.The considered integration procedure is illustrated in Figure 1 (in blue) for a typical cut element; see Ref. [22] for further details.The employed tessellation provides an explicit parametrization of a polygonal approximation of the immersed boundary ∂Ω through the set of boundary faces All The formulations considered in the remainder of this work incorporate stabilization terms formulated on the edges of the background mesh (see Section 2.2), which we refer to as the skeleton mesh Note that the boundary of the background mesh is not part of the skeleton mesh.In addition to the skeleton mesh, we define the ghost mesh as the subset of the skeleton mesh that contain a face of an element intersected by the domain boundary As will be detailed in Section 2.3, the stabilization terms formed on the skeleton and ghost mesh account for stability and ill-conditioning effects related to unfavorably cut elements, as well as for preventing pressure oscillations in equal-order discretizations of the Stokes problem.

Immersogeometric analysis
We consider the immersogeometric analysis of a single-field steady heat-conduction problem and of a two-field viscous flow problem.Both problems are represented by the abstract Galerkin problem Find u h ∈ U h such that: with mesh-dependent bilinear and linear forms, a h and b h , respectively.Note that the superscript h is used to indicate mesh-dependence.The finite dimensional trial and test spaces, U h and V h , are spanned by truncated hierarchical B-spline (THB-spline) [69,94] basis functions of degree k and regularity α constructed over the locally-refined background mesh, viz.
with P k (K) the set of d-variate polynomials on the element K constructed by the tensor-product of univariate polynomials of order k.Truncated hierarchical B-splines, which are illustrated in Figure 2, form a partition of unity and have a reduced support compared to their non-truncated counterpart, which is advantageous from the perspective of system matrix sparsity.Our implementation is based on the open source finite element library Nutils [95].
Since the imposition of strong Dirichlet boundary conditions over the immersed boundary ∂Ω is intractable in the immersogeometric analysis setting, such boundary conditions are imposed weakly through Nitsche's method; see, e.g., Ref. [26].A mesh-dependent consistent stabilization term is introduced in order to ensure the well-posedness of the Galerkin problem (8).

Steady heat conduction
Steady heat conduction is governed by the Poisson problem, which, in dimensionless form, can be formulated as where u is the scalar temperature field, f is a heat source term, q represents the prescribed heat flux on the Neumann boundary, and g is the prescribed temperature on the Dirichlet boundary.The normal gradient is defined as The discretized solution to the strong formulation (10) with the Dirichlet conditions enforced by Nitsche's method is denoted by u h ∈ U h = S k α (T ) ⊂ H 1 (T ), with the corresponding test functions given by v h ∈ V h = U h .We herein consider maximum regularity B-splines, i.e., α = k − 1.The bilinear and linear forms in equation (8) are where β is the Nitsche stabilization parameter.This parameter should be selected and scaled (with the mesh size) appropriately, being large enough to ensure stability, while not being too large to cause a reduction in accuracy (see, e.g., Refs.[31,37]).The ghost-penalty operator in (11a) controls the k th -order normal derivative jumps, indicated by • , over the interfaces of the elements which are intersected by the domain boundary ∂Ω.Since in this contribution B-splines of degree k with C k−1 -continuity are considered, only the k th normal derivative is non-vanishing at the ghost mesh.As will be discussed in detail in Section 2.3, upon approriate selection and scaling (with the mesh size) of γg , a Nitsche stabilization parameter, β, can be selected in such a way that stability of the formulation can be assured independent of the cut-cell configurations.To avoid loss of accuracy, the ghost-penalty parameter γg should also not be too large [96].
Figure 2: Illustration of truncated hierarchical B-splines [69,94] in the immersogeometric analysis setting.The left column shows the hierarchical levels of the mesh T in Figure 1, while the right column illustrates the concept for a one-dimensional immersed domain Ω.The background mesh at the level = 0, • • • , max (with max = 3 in this illustration) is defined as where T A is a regular mesh with mesh size parameter 2 − h.Note that the meshes are nested, in the sense that the domain covered by the physical mesh at level , T , is completely inside that of level − 1, T −1 , i.e., T ⊆ T −1 .The THB-spline basis, H(T ), is constructed by selection and truncation of the basis functions in the B-spline basis At the most refined level, i.e., at = max, all basis functions that are completely inside T max are selected: At the coarser levels, i.e., < max, the functions that are completely inside the domain T but not completely inside the refined domain T +1 are selected and truncated: The truncation operation reduces the support of the B-spline functions by projecting away basis functions retained from the refined levels.The THB-spline basis then follows as H(T ) = ∪ max =0 H(T ).The reader is referred to Ref. [69] for details of THB-spline basis and Ref. [94] for THB-spline basis construction.

Steady viscous flow
Steady viscous flow can be modeled by the Stokes equations, with velocity u, pressure p, constant viscosity µ, body force f , Dirichlet data g and Neumann data t.
By consideration of the solution in the abstract Galerkin problem (8) as a velocity-pressure pair, i.e., , the aggregate bilinear and linear forms corresponding to (12) follow as where For the selection of the Nitsche parameter, β, and ghost stabilization constant, γg , the same arguments apply as for the steady heat conduction problem discussed above.A discussion on the selection and scaling of these parameters for the Stokes problem is presented in Section 2.3.2.An additional stability issue is encountered for the immersed Stokes flow problem (13) on account of the selected equal-order optimal regularity spline spaces of degree k.In the conforming setting, inf-sup stability is achieved by adopting a suitable velocity-pressure pair, e.g., Taylor-Hood [43][44][45][46] or Raviart-Thomas [45,[47][48][49].In the immersed setting, such pairs can still lead to pressure oscillations in the vicinity of cut elements [58].To resolve these pressure oscillations, the immersogeometric skeleton stabilization technique developed in Ref. [59] is applied.This stabilization technique can be regarded as the higher-order continuous version of the method proposed in Ref. [97], which has also been applied in the conforming isogeometric analysis setting [58].
From equation (14c) it is seen that the skeleton stabilization term penalizes jumps in higher-order pressure gradients, where the parameter γs should be selected such that oscillations are suppressed, while the influence of the additional term on the accuracy of the solution remains limited.The purpose of the skeleton stabilization method is to avoid pressure oscillations induced by inf-sup stability problems, allowing for the utilization of identical spaces for the velocity components and the pressure.Since the inf-sup stability problem is not restricted to the immersed boundary, the skeleton stabilization pertains to all interfaces of the background mesh.The appropriate selection and scaling of the skeleton stability parameter is discussed in detail in Section 2.3.2.

Selection of the stabilization parameters: continuity and coercivity of the formulation
Before considering a-posteriori error estimation in Section 3, we first study the continuity and coercivity of the immersed formulations introduced above.We commence with the introduction of the following inequalities: • Using Young's inequality, it follows that for any constant ε > 0 it holds that In combination with the Cauchy-Schwarz inequality, this inequality can be applied to obtain • For any background element K crossed by the boundary ∂Ω, with E = K ∩∂Ω, under the assumption of shape regularity (i.e., provided with an upper bound on the length of the intersection of the boundary within one single element meas (K ∩ ∂Ω)), it holds that (see, e.g., Ref. [98,Lemma 4.2]) where it is noted that this inequality holds for the finite-dimensional space P k of tensor-product polynomials of order k (not for functions in H 1 in general).The constant C T > 0, referred to as the trace inequality constant, is independent of the size of the element, but dependent on the order k.
Note that the right part of the inequality contains the norm over the full background element K, and not just its intersection with the physical domain.
Using inequality (17), the following bound for the normal gradient of u h on the immersed boundary is obtained: with h T defined in Eq. ( 2) and where, with abuse of notation, the constant C T is used to both represent the local trace inequality constant (second line) and its global maximum (third line).
• Norms of functions over the entire background domain T can be bounded by norms over the physical domain Ω and the ghost penalty.Using the ghost-penalty, the gradients on the background mesh are bounded by those in the physical domain.To demonstrate this bound, we split the norm over the background mesh as To show the last inequality, we consider an element K ∈ G which shares the interface F with an element K / ∈ G that completely lies inside Ω, such that the volume integral over the background element K is included in the norm over Ω.We will first demonstrate that the gradients on K are controlled by the ghost penalty and the norms on the physical domain.Later on, elements in G that do not share an interface with an element in T \ G will be considered by means of recursion.To demonstrate that the gradients on K are bound by those in the physical domain, we define the polynomial extension of u h K as the global polynomial ūh K ∈ P k (see Figure 3).Using this extension, the spline function u h on the element K can be decomposed as Let us consider x F as a projection of x on the straight or flat interface F , such that x can be written as x F +x n n F , where x n = (x−x F )•n F .Here, the interface coordinate x F ∈ F is interpreted to be on the side of the element K, and related to the coordinate x ∈ K.The function ũh K has no support on K and has vanishing normal derivatives up to order k at the interface F .By Taylor-series expansion one can infer This splitting is very natural through the use of maximum regularity splines (i.e., ũh K contains all degrees of freedom of K that are independent of K ).

For the polynomial extension ūh
K it holds that where the constant C Q is independent of the mesh size, but dependent on the order of the approximation and the ratio of the size of the elements at either side of the interface.The order-dependence of this constant is illustrated in Figure 3c.The presented results have been computed by solving the generalized eigenvalue problem corresponding to Eq. ( 22).
From the definition of the expansion ũh K in equation ( 21) it follows that with h F the size of K in the direction normal to the interface and where ∇ F denotes the surface gradient in the interface F and where use has been made of the polynomial inequality . The dependence of the constant C F in the inequality (23) on the order is illustrated in Figure 3d.This constant is independent of the mesh size.

Substituting the decomposition (20) in equation (19) yields
Using the inequalities ( 22) and ( 23), and noting that since where To obtain this result, the inequality is first applied to the layer of elements in G that share an interface with the interior mesh T \ G.With control over the gradients in this first layer, the inequality is then applied to a second layer of elements.This recursive application is repeated until all elements in G have been considered.As a result of this recursive application of the ghost inequality, the constant C G depends on the number of layers, which in turn depends on the mesh size.
The boundedness of the gradients on the background mesh finally follows by substitution of ( 25) in (19): • For any p h ∈ U h p , there exists a with a h 2 according to (14b) and ||| • ||| u the velocity energy norm (see Section 2.3.2), for certain positive constants C 1 , C 3 > 0 and a non-negative constant C 2 ≥ 0. Existence of a velocity field w h in accordance with (27) is established in [100,Lemma 3.11] for piecewise linear (k = 1) polynomials.However, this result generalizes to polynomial orders k ≥ 1 and increased continuity of the pressure and velocity spaces.Proof of (27) in the general case is however extensive, and beyond the scope of the present manuscript.

Steady heat conduction
Continuity of the bilinear form (11a) cannot be shown in the H 1 -norm on account of the immersed boundary terms, and coercivity cannot be shown on the infinite-dimensional space.However, with an appropriate selection of the stabilization parameters, continuity and coercivity can be established with respect to the mesh-dependent norm which we refer to as the energy norm.
The bilinear form (11a) is continuous on U h × V h if there exists a constant, C > 0, independent of the mesh size, such that Using the Cauchy-Schwarz inequality, for all Since each of the norms in this expression is bounded from above by the energy norm (28), it follows that a h (u h , v h ) ≤ 5 u h v h .Hence the bilinear form is continuous.The bilinear form (11a) is uniformly (i.e., independent of h) coercive on U h if there exists a constant, c > 0, such that To demonstrate that this is indeed the case, we apply the inequalities ( 16) and ( 26) to obtain Application of the trace inequality (18) and collecting terms then yields for arbitrary ϕ > 0. By selecting element-wise constants 0 < ε ≤ ϕ and 0 < ϕ ≤ for all elements K and interfaces F , where the positive constants β ≥ C T C G and γ g ≥ C F C −1 G are independent of the mesh size.The interface length scale is defined as h F = max (h K , h K ) with K and K being the elements on either side of the interface F .The rational behind this choice is that the ghost stabilization term scales with h 2k−1 F (k ≥ 1) and that hence the larger element size ensures that the stability constant is sufficiently large.

Steady viscous flow
Recalling that for the Stokes problem u h = (u h , p h ), we define the mesh-dependent energy norm as with Continuity of the bilinear form (13a) with respect to this energy norm in the sense of (29) follows directly by application of the Cauchy-Schwarz inequality to all terms in (13a).
With an appropriate selection of the stability parameters for the Stokes problem, it holds that the bilinear form (13a) is inf-sup stable in accordance with sup where c > 0 is referred to as the inf-sup stability constant.To demonstrate this stability property, we recall the splitting of the bilinear form a h according to (13a) and (14).We now take a function where w h depends on p h in accordance with (27), and with some constant α > 0, such that sup Following Section 2.3.1, a h 1 (u h , u h ) is coercive (with constant c u ) and a h 1 (u h , w h ) is continuous (with constant C u ) with respect to the velocity energy norm (33a) in accordance with Eqs. ( 29) and (30), respectively.Hence, where use has been made of From the inequalities (27) it follows that sup which, using Young's inequality (15) with ε = 1, can be reformulated as sup Inf-sup stability as in (34) then holds, provided that We note that the skeleton penalty has two purposes: i) It extends the stability of the pressure field to the background grid T as in (27), in the same way as for the ghost penalty discussed in Section 2.3.1.Since stability is here defined with respect to the L 2 -norm of the pressure field, the skeleton stability constant γs scales with h 2k+1 , following the same reasoning as in Eq. ( 26); ii) It ensures the inf-sup stability for equal-order discretizations, essentially meaning that pressure oscillations in the interior are penalized.This is the reason why this term is applied over the complete skeleton and not only the ghost interfaces.

Error estimation and adaptivity
We study a posteriori error estimation and adaptivity for immersogeometric analysis.In Section 3.1 we first introduce a residual-based error indicator, and elaborate it for the heat conduction problem and viscous flow problem introduced in the previous section.In Section 3.2 the refinement strategy is discussed.

Residual-based error estimation
We propose an error estimator pertaining to the background mesh, T , of the form where the element-wise error indicators, η K , will serve to guide an adaptive refinement procedure.The derivations of the indicators for the heat conduction problem and viscous flow problem will be elaborated in the following sections.
From an abstract perspective, the element-wise error indicators are defined in such a way that the estimator (41) bounds the residual from above as In this expression, the residual and its (dual) norm are defined as The function space V h ⊃ V h corresponds to a suitable extension of V h in such a manner that V h contains an approximation of the solution (possibly the solution itself) that is sufficiently accurate to estimate the error in the approximation u h ∈ V h .An example of such an extended space is an order elevated approximation space on the same mesh and with the same regularity as the space V h , or an approximation space with the same order and the same regularity on a hierarchically refined mesh.The Galerkin approximation problem in The bilinear form a h : V h × V h → R is an extension of the original bilinear form a h according to where the auxiliary symmetric bilinear form s h contains additional stabilization terms.For instance, for an order-elevated space the bilinear form s h contains jumps of higher-order normal derivatives (cf.(11a) and see [101,102]), and for a hierarchically refined mesh s h contains the stabilization terms on the supplementary faces of the mesh.The additional stabilization terms vanish on the original approximation space V h , i.e.
We equip V h with the extended energy norm according to With a suitable choice of the stabilization parameters in a h and s h , the bilinear form a h is weakly coercive and continuous.It is to be noted that this may require that the stabilization parameters in a h are larger than would be required for weak coercivity of a h on V h ×V h .By virtue of ( 44)-( 46) and the weak coercivity and linearity of a h , the following chain of inequalities holds: with the error in the ultimate expression according to e := u h − u h .The chain of inequalities in (48) implies that the error estimator E controls the error in the extended energy norm, |||e||| V h .The reason for defining the residual as a map from V h to V h * is that the stabilization terms in the residual are generally unbounded in the ambient space of the continuum problem, viz.H 1 (Ω) for the steady heat equation and H 1 (Ω, R d ) × L 2 (Ω) for the steady viscous-flow equation.As we will elaborate in Sections 3.1.1and 3.1.2,the refined approximation u h is not required for the calculation of the residual-based estimator (41).The extended space V h merely serves to establish the error-control relation (48).

Steady heat conduction
To derive the error indicators for the steady heat conduction problem introduced in Section 2.2.1, it is first noted that because of Galerkin orthogonality where ṽ = v h − Π h v h ∈ V h and Π h : V h → V h is an interpolation operator [103,104].Note that, for notational convenience, we will drop the diacritic and superscript from v h ∈ V h in the remainder of this section, i.e., v h = v.
Using the definition of the residual (43b) in combination with the definitions of the bilinear and linear forms (11a) and (11b), (reverse) integration by parts yields where The factor 1 2 in the jump and ghost terms accounts for the presence of the associated interfaces in two elements.Using the Cauchy-Schwarz inequality it then follows that Using standard interpolation inequalities [44,105] and the definition of the norm (28), and noting that we consider the functions ṽ and v to be piecewise polynomials, it follows that ṽ where K is the support extension [44] of the element K and K is the element that shares the interface F with element K.The residual can then be bounded as which, using the discrete Cauchy-Schwarz inequality can be rewritten as Using the definition of the residual norm (43b) it follows that with the element error indicators defined as This error indicator reflects that the total element error for all elements that do not intersect the boundary of the domain is composed of the interior residual and the residual term for the jump in the solution normal derivative across the element interfaces.It is noted that for higher-order continuous discretizations, i.e., α > 0, the jump contribution vanishes.For elements that intersect the Neumann boundary, additional error contributions are obtained from the Neumann residual and the ghost penalty residual, while additional Nitsche-related contributions appear for elements intersecting the Dirichlet boundary.

Steady viscous flow
For the Stokes problem introduced in Section 2.2.2, using (reverse) integration by parts, the error indicators in equation ( 41) are obtained by considering the residual (43b) as where ṽ = v − Π h v = (ṽ, q) and Application of the Cauchy-Schwarz inequality gives which, using the inequalities (53) and can be rewritten as Note that the factor 3 in front of the Nitsche residual results from the fact that both terms 2µ (∇ s ṽ) n L 2 (K∩∂Ω D ) and q L 2 (K∩∂Ω D ) are bound by the same norm.Following the same steps as for the heat conduction problem we then obtain the element error indicators as Compared to the error indicators for the heat conduction problem, we here get one additional term to represent the error in the balance of mass, i.e., r h int,p L 2 (K∩Ω) , and one term related to the skeletonstabilization, i.e., r h skeleton L 2 (∂K∩F ) .Moreover, note that the mass and momentum balance terms are scaled with µ − 1 2 and µ 1 2 , respectively, in order to be dimensionally-consistent with the energy norm (32).

Adaptive solution procedure
We employ the residual-based error estimator introduced above in an iterative mesh refinement procedure.In each iteration, for the given mesh we solve the Galerkin problem (8) and subsequently compute the element-wise error indicators (41) (and the corresponding estimator).Based on the indicators, certain elements are then refined, after which the procedure is repeated on the refined mesh.These iterations are continued until a stopping criterion is satisfied.
We consider Dörfler marking [106] to select the elements to be refined.In this marking strategy, the marked set, M, is defined as a minimal set of elements such that with λ a selected fraction of the error estimator.For the considered (truncated) hierarchical spline meshes, refining elements does not necessarily result in a refinement of the approximation space [71,94].To ensure that the approximation space is refined, an additional step is required in which a refinement mask M ⊃ M is defined.To determine the refinement mask, for each element K in the marked set M we determine the support extension and then refine the elements in each support extension which are not smaller than the element K, i.e., During the element refinement procedure the geometry approximation is not altered, as illustrated in Figure 4.In our implementation, the bisectioning depth used to determine the integration subcells is lowered under refinement, resulting in the preservation of the integration subcells under refinement.This ensures that the boundary of the segmented geometry is invariant under mesh refinement.A consequence of this choice is that an element can only be refined up to the level of the integration subcells.Elements requiring refinement beyond the level of the integration subcells are discarded from the refinement list, and the adaptive refinement procedure is stopped if there are no more elements that can be refined.

Cut elements Uncut elements Void elements
Figure 4: Illustration of the refinement procedure for cut elements.The original element is subdivided in integration subcells (blue borders) using the recursive bisectioning procedure detailed in Ref. [22].At the lowest level of bisectioning, a triangulation procedure is employed.After one refinement of the original element, the original element is split into 4 elements, of which one is now an uncut element and the other three are cut elements.The bisectioning depth for the determination of the integration subcells is reduced by one level compared to the original element, so that the subcells remain identical under the element refinement operation.After one further refinement step, each of the four elements in the first refinement is further refined, resulting now also in elements that are void and are hence discarded from the background mesh.

Benchmark simulations
In this section we assess the developed residual-based adaptive refinement technique on a range of numerical experiments.For both the heat conduction problem (Section 4.1) and the viscous flow problem (Section 4.2), both singular and non-singular test cases are considered.For all simulations exact reference solutions are available, allowing for a rigorous study of the stability and accuracy of the developed adaptive immersed isogeometric analysis framework.For all simulations the octree subdivision depth is set equal to the desired maximum number of refinements (see Section 3.2) and the refinement threshold is set to λ = 0.8.Throughout this section, the problems are considered to be in dimensionless form.

Steady heat conduction
We consider the two-dimensional heat conduction problem on a unit square and on a star-shaped domain with a smooth exact solution, and on a domain with a re-entrant corner, for which the exact solution has a reduced regularity (Section 4.1.3).The problems are discretized with linear (k = 1) and quadratic (k = 2) (TH)B-splines using both uniform and adaptive refinement.All examples consider a non-conforming ambient mesh positioned at an angle of 20 degrees (see Figure 5a and Figure 9a), unless specified otherwise.The empirically selected Nitsche and ghost penalty parameters are set to β = 50 and γ g = 10 −(k+2) , respectively.

Unit square
Let Ω = [− 1 2 , 1 2 ] 2 be a unit square with Dirichlet boundary ∂Ω D (see Figure 5a).We define the exact solution of the problem (10) as u(x, y) = sin(πx) + sin(πy), which is shown in Figure 5b.The heat source f corresponding to this exact solution is equal to zero, and the Dirichlet data is set to g = u| ∂Ω D , matching the exact solution.
Figure 6 shows error-analysis results using both uniform and adaptive refinements for the linear case (Figure 6a) and for the quadratic case (Figure 6b).Both refinement procedures start from an initial mesh consisting of 8 × 8 elements covering the ambient domain [−1, 1] 2 .Optimal convergence rates are obtained for both the error in the L 2 -norm (i.e., O(n − 1 2 (k+1) )) and in the H 1 -norm (i.e., O(n − 1 2 k )), with n denoting the number of degrees of freedom.Moreover, as the number of refinement steps increases, the energy norm and H 1 -norm of the error coincide, indicating that the error is dominated by the H 1 -semi-norm contribution in Eq. (28).The estimator (41) is observed to converge at the same rate as the energy norm, bounding the energy norm from above, consistent with Eq. (48).Because of the smooth solution (67), the refinement pattern following from the adaptive refinement procedure closely resembles the uniform refinements, as observed from the close correspondence between the error results for the uniform and adaptive simulations in Figure 6.

Star-shaped domain
To study the sensitivity of the adaptive simulation framework to the cut-cell configurations, we consider the star-shaped domain shown in Figure 7a for various orientation angles ϑ.The star-shaped domain is constructed using the level set function with R 1 = 0.6, R 2 = 0.2 and n fold = 5 [107].On the boundary of the domain, the Dirichlet data is set equal to the same exact solution (67) as in the previous example.For all orientations, an initial mesh of 10 × 10 elements covering the ambient domain [−1, 1] 2 is considered, after which local refinements using second-order THB-splines are performed until the smallest elements have been refined six times.
Figures 7b-7f show the error u − u h after completion of the refinement procedure.These figures convey that both the error and the refinement pattern are similar for all orientations.This is corroborated by the results in Figure 8, which indicates that both the number of degrees of freedom and the errors (in various norms) are insensitive to the orientation angle.

Re-entrant corner
To study the behavior of the adaptive simulation strategy for problems with (weakly) singular solutions, we consider a domain with a re-entrant corner, as shown in Figure 9a.The data on the Dirichlet and Neumann boundaries, u| ∂Ω D = g = 0 and ∂ n u| ∂Ω N = q, is set to match the exact solution [62,71] u(x, y) = (x 2 + y 2 ) The convergence behavior of the L 2 -error, H 1 -error, energy norm error (28) and the residual-based estimator (57) is studied for uniform refinement and residual-based adaptive refinement.Both refinement procedures start from an initial mesh of 10 × 10 elements formed on the ambient domain [− 3 2 , 3 2 ] 2 .The convergence results for first and second order B-splines are shown in Figure 10a and Figure 10b, respectively.
Under uniform refinement, the convergence rates are impeded by the weak singularity at the re-entrant corner.For the L 2 -error and H 1 -error, suboptimal rates of O(n − 2 3 ) and O(n − 1 3 ) are observed, which is in agreement with the expected rates [108].These rates are independent of the order of the approximation, as the regularity of the exact solution limits the rate already for the linear case.As for the cases considered above, the energy error and estimator follow the convergence of the H 1 -error.Using the adaptive refinement strategy with linear basis functions, the optimal rates of O(n −1 ) and O(n − 1 2 ) are recovered for the L 2 -error and H 1 -error, respectively.For the quadratic case, rates that are substantially higher than the theoretical rates are observed.We attribute this to pre-asymptotic behavior, in which the refinement pattern as shown in Figure 10 is strongly focused on the re-entrant corner singularity.After the first two steps, the errors become dominated by the singularity at the re-entrant corner, which results in the further refinement of the few elements in the vicinity of the corner.These refinements do reduce the error, while they only introduce a limited number of additional degrees of freedom.The observed flattening in the rate of the L 2 -error in the quadratic case is caused by the refinement reaching the maximum level in the elements in the corner, which causes the marking strategy to tag elements that do not carry the largest error contributions.

Steady viscous flow
We regard the two-dimensional Stokes flow problem on a quarter annulus ring domain with a smooth solution and on the above-introduced re-entrant corner domain with a singular solution.We consider equalorder discretizations for the velocity and pressure fields using optimal regularity (TH)B-splines of degree k = 1 and k = 2.For the Nitsche and ghost-penalty parameter the same settings are used as for the Laplace problem considered above, i.e., β = 50 and γ g = 10 −(k+2) .In addition, a skeleton-penalty parameter of γ s = 10 −(k+1) is used for all simulations.

Quarter annulus ring
We consider an annulus ring domain Ω = {(x, y) ∈ R 2 >0 : R 2 1 < x 2 + y 2 < R 2 2 } with inner radius R 1 = 1, outer radius R 2 = 4, Dirichlet boundary ∂Ω D and Neumann boundary ∂Ω N , as shown in Figure 12a.The Dirichlet data g and Neumann data t are prescribed in accordance with the divergence-free manufactured solution [58] The body force f in the Stokes problem (12) is determined based on this manufactured solution, with the viscosity set to µ = 1.).The error in the energy norm ( 32) is observed to converge with the same rate as the H 1 -norm velocity error and L 2 -norm pressure error, which is in agreement with the definition of the energy norm.As expected, the error estimator bounds the error in the energy norm from above.
Although optimal convergence rates are obtained using uniform refinements, the adaptive refinement procedure is observed to substantially improve the error for a fixed number of degrees of freedom.This behavior is explained by the observed refinement patterns, as shown in Figure 14.Although the exact solution ( 69) is smooth, in particular the steep gradients in the velocity solution lead to local refinements.This effectively reduces the error when compared to a uniform refinement with a similar number of degrees of freedom.

Re-entrant corner
As a final benchmark problem we consider the Stokes problem (12) on the re-entrant corner domain with mixed Dirichlet and Neumann boundaries introduced above, as shown in Figure 9a.The weakly singular exact solution is taken from Ref. [109] as with constants α = 856399/1572864 and ω = 3 2 π, and with The exact pressure and velocity fields are illustrated in Figure 15.The corresponding Stokes problem (12) is considered with the viscosity set to µ = 1, no body force, f = 0, a no slip condition on Γ D , such that u D = 0, and the Neumann data g on Γ N matching the exact solution.
Figure 16 displays the error convergence results obtained using uniform and adaptive refinements, for both linear and quadratic (TH)B-splines.As for the Laplace case, the weak singularity in the exact solution (70) limits the convergence rate when uniform refinements are considered.Using adaptive mesh refinement  results in a recovery of the optimal rates in the case of linear basis functions, with even higher rates observed for the quadratic splines on account of the highly-focussed refinements resulting from the residual-based error estimator as observed in Figure 17.

Scan-based simulations
In this section we apply the developed adaptive immersed isogeometric analysis framework in the context of scan-based analysis.We consider the viscous flow problem on a two-dimensional image domain and on a three-dimensional patient-specific problem based on a µCT-scan of a carotid artery, represented by grayscale voxels.The primary purpose of the two-dimensional setting is to test the scan-based analysis framework.For all simulations, the octree subdivision depth is set equal to 8 in two dimensions and 3 in three dimensions.The refinement threshold related to the Dörfler marking is set to λ = 0.8.
Our scan-based analysis workflow is illustrated in Figure 18.The first step in this workflow is to smoothen the original grayscale voxel data using a convolution operation on a B-spline basis formed on the voxel grid [82].Since this smoothing operator behaves as a Gaussian filter, geometric features that are similar in size to the voxels can be lost [110].To avoid this loss of features, the topology-preservation procedure proposed in Ref. [110] is employed.This procedure locally refines the convolution basis to retain small geometric features in the smoothing procedure.Once the smooth level set representation has been obtained, the octree segmentation procedure with mid-point tessellation of Ref. [22] is used to obtain the immersed geometry represented on an ambient domain mesh.It is important to note that this ambient domain mesh, on which the solution to the flow problem is computed, can be chosen independently of the voxel size, and hence it is independent of the mesh on which the level set function is constructed.
The considered computational domain is illustrated in Figure 19a.Neumann conditions are imposed on the inflow and outflow boundaries, with the traction on the inflow boundary acting in the normal direction with a traction data, t = −pn, where p is the pressure magnitude.Homogeneous Dirichlet conditions are imposed along the immersed boundaries in accordance with the no slip condition.It is to be noted that a Neumann condition at an inflow boundary generally leads to an ill-posed boundary value problem for the Navier-Stokes equations, but the Stokes problem is well-posed.In all simulations we consider second-order (k = 2) (TH)B-splines and set the stabilization parameters to β = 100, γ g = 10 −(k+2) and γ s = 10 −(k+1) , which have been determined empirically.

Two-dimensional prototypical geometry
To test the developed methodology in the scan-based setting, we first consider the prototypical twodimensional geometry shown in Figure 19a, which is constructed from 32 × 32 grayscale voxel data.The ambient domain, which matches the scan window, is taken as a unit square (L = 1) which is covered by an 8 × 8 elements ambient mesh.The viscosity is set equal to µ = 1 and the pressure to p = 1.
Various steps in the adaptive refinement procedure are depicted in Figure 20.In the first step virtually all elements covering the flow domain are refined, indicating that the initial mesh of only 8×8 elements is too coarse to resolve the solution globally.After the first refinement step, the refinement strategy starts to focus on the regions where the errors are largest, i.e., near boundaries and narrow sections, as also illustrated in Figure 19b.Under further refinement, the procedure resolves prominent solution details, most importantly the (Poiseuille-like) profile in the carotid part of the artery and the velocity profiles at the inflow and outflow boundaries.
Further results of the viscous flow problem solved using uniform and adaptive refinements are shown in Figure 21 in the form of the flux through the left and right outflow channels.The minor difference in results on the initial mesh (left-most points) are caused by a different selection of the octree-depth for the uniform and adaptive simulations.Both methods are observed to converge to the same fluxes under refinement, but an excellent approximation of the reference solution (computed on a uniform overkill refinement, consistent with the result reported in Ref. [110]) is obtained by means of the adaptive mesh refinement procedure using substantially fewer degrees of freedom than for uniform refinements.This is consistent with the observations on the velocity field discussed above, where in particular the ability of the adaptive refinement procedure to resolve the flow in the carotid part is essential.

Three-dimensional patient-specific geometry
To demonstrate the residual-based adaptivity procedure in a real scan-based setting, we consider the patient-specific carotid artery used in Ref. [110].The geometry of the carotid artery is obtainted from CT-scan data containing 80 slices of 85 × 70 voxels.The size of each voxel is 300 × 300 µm 2 and the distance between the slices is 400 µm.The total size of the scan domain is 25.6×21.1×32.0mm 3 .We set the viscosity to 4mPa s and pressure to 17.3kPa (130 mm of Hg).
Simulation results for this problem are shown in Figure 22.Note that for the considered scan data, the application of the topology-preservation algorithm in Ref. [110] is essential, as otherwise the narrow channel section in the right artery would disappear.The simulation results are based on a 24 × 24 × 24 ambient domain mesh of 25.6 × 21.1 × 32.0 mm 3 and an octree depth of three.In this setting, after two refinements,      an element is of a similar size as the voxels.The need to substantially refine beyond the voxel size is, from a practical perspective, questionable, as the dominant error in the analysis will then be related to the scan resolution and the segmentation procedure.In this sense, the constraint of not being able to refine beyond the octree depth is not a crucial problem in the considered simulations.Different steps in the adaptive refinement procedure are illustrated in Figures 22 and 23.In all the refinement steps, the refinement strategy starts to focus on the regions where the errors are largest, i.e., near the stenosed section (i.e., the narrow region at the right artery) and at the outflow section of the left artery.Under local mesh refinement, the procedure resolves prominent solution details, most importantly the velocity field in the left artery and near the stenotic part of the right artery.
The flux at the outlet of the arteries is shown in Figure 24, which is computed with the velocity field obtained by solving the flow problem using adaptive refinements.The solution of the flux in the left artery is observed to gradually converge toward a value of just over 5100 [mm 3 /s].For the right artery, the maximum refinement depth is reached after the second refinement step.As a result, the flux in the right artery does then not substantially change anymore.At this point, the element sizes in the vicinity of the stenotic artery are similar in size to the voxels.The error then becomes dominated by the geometry reconstruction procedure, which also explains why the observed flux in the right artery deviates from the uniform mesh results in Ref. [110], viz.max = 2 instead of the presently applied max = 3.It is observed that the adaptive procedure terminates after 4 refinement steps, because of reaching the maximum refinement level in all the elements tagged for refinement.At this point, the adaptive simulation uses 12, 816 DOFs, which is substantially lower than the number of DOFs required using uniform refinements [110], which amounts to approximately 10 5 .

Concluding remarks
In the immersed (isogeometric) analysis framework, the geometry representation is decoupled from the discretization.This enables the consideration of spline basis functions on complex volumetric domains, for which boundary-fitting discretizations cannot easily be obtained.Moreover, the decoupling of the geometry and the discretization allows one to have a globally accurate representation of the geometry, but only to refine the mesh in places where the errors are large.Such local mesh refinements have the potential to provide a significant efficiency gain compared to uniform meshes.The adaptive simulation strategy proposed in this work automatically refines the elements in places that significantly contribute to the error in the energy norm.The developed error estimation and adaptivity strategy is based on residual-based error estimation, which is well-established in traditional finite elements and has been successfully applied in boundary-fitting isogeometric analysis.In the considered immersed setting, the residual-based error estimation and adaptivity framework requires the incorporation of the stabilization terms for the weakly imposed Dirichlet boundary conditions, and, in the case of the (mixed) Stokes flow problem, for the treatment of equal-order discretizations of the velocity-pressure pair.Adequate scaling of the stabilization constants with the mesh size is essential for the adaptive procedure to be effective.In particular the order dependence of the stabilization constants and the definition of the local element sizes must be treated adequately.
In contrast to residual-based error estimation for boundary-fitting finite elements and isogeometric analysis, in the stabilized immersed setting it is not evident that the residual-based error estimator bounds the error in the energy norm from above.This is a consequence of the absence of an h-independent weak formulation.In this work, it is reasoned, however, that under the assumption of sufficient smoothness, the residual is expected to be useful in the setting of an adaptive refinement strategy.For all numerical simulations considered, including simulations with reduced regularity, it is observed that the error estimator does provide an upper bound to the error in the energy norm.A rigorous study regarding the relation between the residual and the actual error is warranted.
It is demonstrated that the developed adaptive simulation strategy is particularly useful in a scan-based analysis setting, where manual selection of refinement zones is impractical.When used in combination with advanced image segmentation procedures to obtain a smooth geometry representation while preserving small geometric features, the developed adaptive refinement strategy optimally leverages the advantageous approximation properties of splines for geometrically and topologically complex domains.The adaptivity strategy results in a simulation workflow that is capable of obtaining reliable, error-controlled, results with limited user interaction.
The developed adaptive solution strategy is elaborated for the Laplace problem and the Stokes problem.For other problems, such as, for example, Navier-Stokes or Cahn-Hilliard problems, the starting point of the derivation of the error-estimator remains the same.The estimators are problem-specific, however, and hence need to be elaborated for such problems.The same holds for the consideration of additional or alternative stabilization techniques, specifically when these alter the Galerkin form of the problem.
faces E ⊂ ∂Ω N (respectively E ⊂ ∂Ω D ) are assigned to a set of Neumann faces T ∂Ω N (respectively Dirichlet faces T ∂Ω D ).In general, a single polygon face can overlap with both the Neumann and the Dirichlet boundary.Let us note that in an adaptive refinement procedure, the refinements can serve to provide an increasingly accurate approximation of the transition between the Neumann and Dirichlet boundary.

Figure 3 :
Figure 3: Panel (a) is an illustration of a second order B-spline on an element K ∈ G and its adjacent element K with an interface F .Panel (b) is its second order gradient in the direction normal to the interface (with en the unit vector in the normal direction).Panels (c) and (d) show the dependence of the constants in (22) and (23) on the order k.

Figure 5 :
Figure 5: (a) Problem setup, and (b) the exact solution u(x, y), Eq. (67), for the Laplace problem on the unit square domain.

Figure 6 :
Figure 6: Error convergence results for the Laplace problem on the unit square domain under residual-based adaptive refinement (solid) and uniform refinement (dashed) for linear (k = 1) and quadratic (k = 2) basis functions.

Figure 7 :
Figure 7: (a) Problem setup, and (b)-(f) contour plots of the error, u − u h , for the Laplace problem on the star shaped domain at the end of 6 adaptive refinement steps for different angles of mesh rotation ϑ.

Figure 8 :Figure 9 :
Figure 8: (a) Degrees of freedom, and (b) error norms for the Laplace problem on the star shaped domain after 6 adaptive refinement steps for different angles of mesh rotation ϑ.

Figure 10 :
Figure 10: Error convergence results for the Laplace problem on the re-entrant corner domain under residual-based adaptive refinement (solid) and uniform refinement (dashed) for linear (k = 1) and quadratic (k = 2) basis functions.

Figure 11 :
Figure11: Evolution of the mesh using the adaptive refinement procedure for the Laplace problem on the re-entrant corner domain using k = 2.

Figure 13 :
Figure 13: Error convergence results for the Stokes problem on the quarter annulus ring domain under residual-based adaptive refinement (solid) and uniform refinement (dashed) for linear (k = 1) and quadratic (k = 2) basis functions.

Figure 13
Figure13displays the convergence results for the annulus ring problem.Both the uniform refinement results and the adaptive refinement results are obtained starting from a 9 × 9 uniform mesh on the ambient domain [0, R 2 ] 2 = [0, 4] 2 .A good resemblance with the optimal rates of O(n − 1 2 k ) in the velocity H 1 -norm and pressure L 2 -norm is observed, and, as expected, the rate of the velocity L 2 -error is O(n − 1 2 (k+1) ).The error in the energy norm (32) is observed to converge with the same rate as the H 1 -norm velocity error and L 2 -norm pressure error, which is in agreement with the definition of the energy norm.As expected, the error estimator bounds the error in the energy norm from above.Although optimal convergence rates are obtained using uniform refinements, the adaptive refinement procedure is observed to substantially improve the error for a fixed number of degrees of freedom.This behavior is explained by the observed refinement patterns, as shown in Figure14.Although the exact solution (69) is smooth, in particular the steep gradients in the velocity solution lead to local refinements.This effectively reduces the error when compared to a uniform refinement with a similar number of degrees of freedom.

Figure 14 :
Figure14: Evolution of the mesh using the adaptive refinement procedure for the Stokes problem on the quarter annulus ring domain using k = 2.

Figure 15 :
Figure 15: (a) Velocity magnitude and streamlines, and (b) pressure for the exact solution (70) to the Stokes problem on the re-entrant corner domain.Because of the singular solution, the pressure color bar is truncated to the range −10 and 10.

Figure 16 :
Figure 16: Error convergence results for the Stokes problem on the re-entrant corner domain under residual-based adaptive refinement (solid) and uniform refinement (dashed) for linear (k = 1) and quadratic (k = 2) basis functions.

Figure 17 :
Figure 17: Evolution of the mesh using the adaptive refinement procedure for the Stokes problem on the re-entrant corner domain using k = 2.

Figure 18 :Figure 19 :
Figure 18: Illustration of the scan-based analysis workflow.The original grayscale image in panel (a) is converted to a level set function, shown in panel (b), which is constructed using the topology-preserving segmentation algorithm of Ref.[110].The trimmed geometry, shown in panel (c), is then extracted using the recursive bi-sectioning strategy with mid-point tessellation of Ref.[22].

Figure 20 :
Figure 20: Evolution of the mesh and (magnitude of the) velocity field during the adaptive refinement process for the viscous flow in two dimensions.

Figure 21 :
Figure 21: Mesh convergence of the outflow flux at the (a) left and (b) right channel of the domain in Figure 19a using adaptive (solid) and uniform (dashed) mesh refinements.

Figure 23 :
Figure 23: Evolution of the mesh during the adaptive refinement process for the patient-specific viscous flow problem.

Figure 24 :
Figure 24: Mesh convergence of the flux at the left and right outflow boundary using adaptive mesh refinements for the patient-specific viscous flow problem.