CutIGA with Basis Function Removal

We consider a cut isogeometric method, where the boundary of the domain is allowed to cut through the background mesh in an arbitrary fashion for a second order elliptic model problem. In order to stabilize the method on the cut boundary we remove basis functions which have small intersection with the computational domain. We determine criteria on the intersection which guarantee that the order of convergence in the energy norm is not affected by the removal. The higher order regularity of the B-spline basis functions leads to improved bounds compared to standard Lagrange elements.


Introduction
Background and Earlier Work. CutFEM and CutIGA, are methods where the geometry of the domain is allowed to cut through the background mesh in an arbitrary fashion, which manufactures so called cut elements at the boundary. This approach typically leads to some loss of stability and ill conditioning of the resulting stiffness matrix that must be handled in some way. Several approaches have been proposed: • Gradient jump penalties or some related stabilization term, see [3] and [4].
• Adding a small amount of extra stiffness to each active element as is done in the finite cell method, see [7] and [12].
• Element merging where small elements are associated with a neighbor which has a large intersection. For DG methods see [11] and for CG methods see [1].
• Basis function removal where basis functions with support that has a small intersection with the domain are removed. For the case of isogeometric spline spaces see [8].
For a general introduction to CutFEM we refer to the overview paper [4] and for an introduction to isogeometric analysis we refer to [6].
New Contributions. We investigate the basis function removal approach based on simply eliminating basis functions that has a small intersection with the domain in the context of isogeometric analysis, more precisely we employ B-spline spaces of order p with maximal regularity C p−1 . To this end we need to make the meaning of small intersection precise and our guideline will be that we should not lose order in a given norm. In particular, we consider the the error in the energy norm and show that we may remove basis functions with sufficiently small energy norm and still retain optimal order convergence. We also quantify the meaning of a basis function with sufficiently small energy norm in terms of the size of the intersection between the support of the basis function and the domain. In order to measure the size of the intersection we consider a corner inside the domain and we let δ i , i = 1, . . . , d with d the dimension, be the distance from the corner to the intersection of edge E i with the boundary. If there is no intersection δ i = h. We then identify a condition on δ i in terms of the mesh parameter h which guarantees that we have optimal order convergence in the energy norm. The energy norm of the basis functions may be approximated by the diagonal element of the stiffness matrix and we propose a convenient selection procedure based on the diagonal elements in the stiffness matrix which is easy to implement.
We also derive the condition on δ i corresponding to the W 1 ∞ norm, which will be tighter since the norm is stronger and here we also need the continuity of the derivative of the basis functions. We discuss the approach in the context of standard Lagrange basis function where we note that we get much a tighter condition on δ i in the energy norm and in the W 1 ∞ norm we find that it is not possible to remove basis functions. We impose Dirichlet conditions weakly using nonsymmetric Nitsche, which is coercive by definition. Since the energy norm used in the nonsymmetric Nitsche method does not control the normal gradient on the Dirichlet boundary we do however need to add a standard least squares stabilization term on the elements in the vicinity of the boundary. Note that this term is element wise in contrast to the stabilization terms usually used in CutFEM.
When symmetric Nitsche is used to enforce Dirichlet boundary conditions stabilization appears to be necessary to guarantee that a certain inverse estimate holds. This bound is not improved by the higher regularity of the splines and will not be enforced in a satisfactory manner by basis function removal.
Outline: In Section 2 we introduce the model problem and the method, in Chapter 3 we derive properties of the bilinear form, define the interpolation operator, define the criteria for basis function removal, derive error bounds, and quantify δ in terms of h for various norms, and finally in Section 4 we present some illustrating numerical examples.

Model Problem
Let Ω be a domain in R d with smooth boundary ∂Ω and consider the problem: find u : Ω → R such that For sufficiently regular data there exists a unique solution to this problem and we will be interested in higher order methods and therefore we will always assume that the solution satisfies the regularity estimate for s ≥ 2. Here and below a b means that there is a positive constant C such that a ≤ Cb.

The Finite Element Method
The B-Spline Spaces.
• Let T h , h ∈ (0, h 0 ], be a family of uniform tensor product meshes in R d with mesh parameter h. • Let V h = C p−1 Q p (R d ) be the space of C p−1 tensor product B-splines of order p defined on T h . Let B = {ϕ i } i∈ I be the standard basis in V h , where I is an index set.
• Let B = {ϕ ∈ B : supp(ϕ) ∩ Ω = ∅} be the set of basis functions with support that intersects Ω. Let I be an index set for B. Let V h = span{B} and let T h = {T ∈ T h : T ⊂ ∪ ϕ∈B supp(ϕ)}. An illustration of the basis functions in 1D is given in Figure 1.
• Let B = B a ∪ B r be a partition into a set B a of active basis functions which we keep and a set B r of basis functions which we remove. Let I = I a ∪ I r be the corresponding partition of the index set. Let V h,a = span{B a } be the active finite element space.

Remark 2.1
To construct the basis functions in V h we start with the one dimensional line R and define a uniform partition, with nodes x i = ih, i ∈ Z, where h is the mesh parameter, and elements The basis functions ϕ i,p are then defined by the Cox-de Boor recursion formula we note that these basis functions are C p−1 and supported on [x i , x i+p+1 ] which corresponds to p + 1 elements, see Figure 1. We then define tensor product basis functions in R d of the form The Nonsymmetric Method. Find u h,a ∈ V h,a such that The forms are defined by with positive parameters β and τ . Furthermore, we used the notation and with δ ∼ h and B δ (x) the open ball with center x and radius δ. We note that it follows from (2.14) that U δ (∂Ω D ) ⊂ T h,D .

Galerkin Orthogonality. It holds
In practice, T h,D may be taken as the set of all elements that intersect the Dirichlet boundary ∂Ω D and their neighbors, i.e.

Remark 2.3 (The Symmetric
Method) The symmetric version of (2.8) takes the form: find u h,a ∈ V h,a such that The forms are defined by where β and γ are positive parameters, F h,D is the set of interior faces which belong to an element in T h (∂Ω D ), and D n F = n F · ∇ is the directional derivative normal to the face F . The stabilization term s h,sym provides the control where we note that we indeed obtain control on the full elements T ∈ T h (∂Ω D ). The control (2.21) is employed in the proof of the coercivity of A h in the symmetric case. More precisely, (2.21) is used as follows where we used an inverse inequality in the first estimate In the symmetric formulation we stabilize to ensure that coercivity holds and this stabilization also implies that the resulting linear system of equations is well conditioned. Therefore, in the symmetric case, we do not employ basis function removal on the Dirichlet boundary.

Basic Properties of A h
Energy norm. Define the norms where V = H 2 (Ω). This result follows directly from the definition and the fact that the parameters τ ≥ 0 and β > 0.

Continuity. The form A h is continuous
Proof. First we note that We proceed with an estimate of the first term on the right hand side. To that end let χ : Ω → [0, 1] be a smooth function such that where U δ (∂Ω δ ) is defined in (2.15). Splitting the term (∇v, ∇w) Ω using χ and then applying Green's formula for the term in the vicinity of ∂Ω D followed by some obvious bounds give Next using the bound see [5], we conclude that Finally, choosing δ ∼ h and using the fact that which in combination with (3.6) concludes the proof.

Interpolation Error Estimates
Definition of the Interpolant. There is an extension operator E : see [9]. Define the interpolant by where π Cl,h is a Clement type interpolation operator onto the spline space. We have the expansion where (π h (Ev)) i is the coefficient corresponding to basis function ϕ i . We define the interpolant on the active and removed finite element spaces by We then have π h (Ev) = π h,a (Ev) + π h,r (Ev) (3.28) Below we simplify the notation and write v = Ev and π h (Ev) = π h v. Proof. Using the identity π h v = π h,a v + π h,r v and the triangle inequality by standard spline interpolation results [2]. To estimate the second term on the right hand side we introduce the scalar product v, w h, = (∇v, ∇w) Ω + h(n · ∇v, n · ∇w) associated with the norm ||| · ||| h, . Expanding π h,r v in the basis B r we get and we have the bound j∈Ir δ ij ≤ (2p + 1) d (3.41) • We used the L ∞ (N h (Ω)) stability of the interpolant π h and then the L ∞ stability of the extension operator and finally the Sobolev embedding theorem

Error Estimate
We have the following error estimate.
Remark 3.1 Note that if we take τ = 0, i.e. we use the method without least squares stabilization in the vicinity of the Dirichlet boundary. We may still derive an error estimate as follows

50)
Now we note that and thus we obtain the bound where the second term on the right hand side is a residual term involving the computed solution u h . The resulting bound is thus of a priori -a posteriori type. One may estimate the residual term on elements in the interior of Ω but for elements which are cut we do not have access to the required inverse estimate.

Bounds in Terms of the Geometry of the Cut Elements
In this section we derive a criterion in terms of the geometry of the cut support of the basis function which implies (3.29). This criterion will in general not be used in practice but it provides insight into the effect of the higher order regularity of the B-splines.
Assuming that there are h −(d−1) such elements we have the estimate up to constants and in local coordinates with origo X 0 , and

Condition (3.57) thus takes the form
For Lagrange basis functions we instead have |Dϕ(x)| ∼ h −1 and we therefore obtain the condition An illustration of both B-spline and Lagrange basis functions in this setting is given in Figure 2. Comparing (3.60) and (3.61) we note that the condition is much stronger for the Lagrange functions and higher order p.
The 1D Case: Max Norm. The difference between the B-splines and Lagrange basis functions is even more drastic if we consider instead evaluating the max norm of the derivative. Then for B-splines we have while for Lagrange basis functions which in the latter case can not be controlled by decreasing δ, see Figure 2. Thus for Lagrange basis functions we get a pointwise error of order h −1 if we remove a basis function while for quadratic and higher order B-splines we may retain optimal order local accuracy by choosing The 2D Case: Energy Norm. We now extend our calculation to the 2D case. The higher dimensional case can be handled using a similar approach. Let X 0 be a vertex of supp(ϕ) which reside in the interior of Ω. Let {e i } d i=1 be an orthonormal coordinate system centered at X 0 and with basis vectors e i , and coordinates x i , aligned with the edges {E i } d i=1 of supp(ϕ) which originates at X 0 , see Figure 3. Using the local coordinates in the vicinity of X 0 we have the expansions Let δ i = X i − X 0 R d be the distance from the vertex X 0 to the intersection X i of edge E i with the boundary ∂Ω. Assume that Integrating over [0, Condition (3.57) thus takes the form See Figure 4 for an illustration of this condition.
The 2D Case: Max Norm. Starting from the expansion (3.66) and observing that for small enough δ parameters |∇ϕ| 2 is increasing when we move out from the vertex. Using assumption (3.67) we thus conclude that which we may write in the form See Figure 4 for an illustration of this condition.

Linear Elasticity
While we for simplicity use the Poisson model problem in the above analysis the same analysis holds also for other second order elliptic problems which may be of more practical interest. We therefore in the numerical results apply our findings to the linear elasticity problem: find the displacement u : where the stress and strain tensors are defined by σ(u) = 2µ (u) + λtr( (u)), with Lamé parameters λ and µ; f , g N , g D are given data; a ⊗ b is the tensor product of vectors a and b with elements (a ⊗ b) ij = a i b j .
The Nonsymmetric Method for Linear Elasticty. Find u h,a ∈ [V h,a ] d such that The forms are defined by where with positive parameters β and τ . Furthermore, the energy norm is defined  A Neumann Problem. To illustrate the selection of spline basis functions to remove we first consider a pure Neumann problem with the geometry presented in Figure 5a. The domain is symmetrically pulled from the left and the right using a unitary traction load. We assume a linear isotropic material with an E-modulus of E = 100 and a Poisson ratio of ν = 0.3. To ensure the discretized problem is well posed we seek solutions orthogonal to the rigid body modes by using Lagrange multipliers.
A Manufactured Problem. To numerically estimate convergence rates we use the following manufactured problem from [10]. The geometry and the solution is given by see Figure 5b. Assuming a linear isotropic material with the material parameters of steel we deduce expressions for the input data f , g N and g D . Note that while this problem does include a Dirichlet boundary ∂Ω D we in our current implementation neglect the least squares term in the vicinity of ∂Ω D , i.e. we choose τ = 0.

Illustration of the Selection Procedure
We utilize the selection procedure based on the stiffness matrix proposed in Section 3.2. Some realizations of this selection are visualized in Figure 6 where we note that the selection becomes more restrictive as the mesh size decreases. This is a natural effect as the selection procedure is developed to ensure optimal approximation properties of the active spline space V h,a . We also note that when increasing spline order more basis functions are removed when using the same constant in the tolerance tol = ch p . This can also be seen in Figure 7 where we investigate how the choice of this constant effects the number of removed basis functions. In Figure 8 we note that the use of basis removal is quite effective and also give better quality stresses along the boundary. Figure 6: Four realizations of removed basis functions using on the stiffness matrix based selection procedure described in Section 3.2; all using the same constant c = 0.01 for the tolerance tol = ch p × √ E in (3.29). Each cross marks a removed basis function and the domain of its support is visualized in pink. In (a)-(c) we note that the selection becomes more restrictive with smaller mesh size h. Comparing (b) and (d) we also note that more basis functions typically may be removed as the spline order increases.

Convergence
To estimate the convergence we use the manufactured problem described in Section 4.1 and the cut situations are induced by rotating the background grid π/7 radians as illustrated by the mesh with removed basis functions in Figure 9 together with the corresponding numerical solution. In Figure 10a we present convergence studies in energy norm for various choices of the constant c in the tolerance tol = ch p × √ E used in the selection procedure. As can be seen, a larger constant naturally means a larger error, but the convergence rates remain optimal. The stiffness matrix condition numbers corresponding to these convergence studies is presented in Figure 10b. It can be noted that while basis removal greatly reduce the size of the condition numbers, basis removal alone does not yield an optimal scaling of O(h −2 ).

Conclusion
We have shown that: • Basis function removal can be done in a rigorous way which guarantees optimal order of convergence and that the resulting linear system is not arbitrarily close to singular. These results critically depend on the smoothness of the B-spline spaces.
• Basis function removal is easy to implement and efficient since there is no fill-in in the stiffness matrix as is the case in for instance face based stabilization. Furthermore, basis function removal is consistent in contrast to the finite cell method.
We note however that even though the stiffness matrix is not arbitrarily close to singular the resulting condition number will in general be worse than O(h −2 ), which is the optimal scaling for standard finite element approximation of second order elliptic problems and therefore a direct solver or preconditioning in combination with an iterative solver is necessary in practice.   : Example of numerical solution using C 1 Q 2 splines and mesh size h = 0.1. The mesh is rotated π/7 radians to induce cut situations and the removed basis functions are selected using the tolerance tol = 0.01h 2 × √ E. (b) Condition number Figure 10: Convergence in ||| · ||| h norm and condition numbers for the manufactured problem using basis removal with C 1 Q 2 -spline basis. The tolerances used in the selection procedure is tol = ch p × √ E and we note that the choices c = 10 −1 and c = 10 −2 give no visible difference in the error compared to using the full approximation space (c = 0).