CutIGA with basis function removal
Daniel Elfverson1,
Mats G. Larson1 &
Karl Larsson ORCID: orcid.org/0000-0001-7838-13071
Advanced Modeling and Simulation in Engineering Sciences volume 5, Article number: 6 (2018) Cite this article
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.
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].
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 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 \(\delta _i\), \(i=1,\ldots ,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 \(\delta _i = h\). We then identify a condition on \(\delta _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 \(\delta _i\) corresponding to the \(W^1_\infty \) 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 \(\delta _i\) in the energy norm and in the \(W^1_\infty \) 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.
In "The model problem and method" section 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 \(\delta \) in terms of h for various norms, and finally in "Numerical results" section we present some illustrating numerical examples.
The model problem and method
Model problem
Let \(\Omega \) be a domain in \(\mathbb {R}^d\) with smooth boundary \(\partial \Omega \) and consider the problem: find \(u: \Omega \rightarrow \mathbb {R}\) such that
$$\begin{aligned} -\Delta u =f \qquad {\text {in}}\, {\Omega } \end{aligned}$$
$$\begin{aligned} u =g_D \qquad {\text {on}}\, {\partial \Omega _D} \end{aligned}$$
$$\begin{aligned} n\cdot \nabla u =g_N \qquad {\text {on}}\, {\partial \Omega _N} \end{aligned}$$
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
$$\begin{aligned} \Vert u \Vert _{H^{s}(\Omega )} \lesssim 1 \end{aligned}$$
for \(s\ge 2\). Here and below \(a \lesssim b\) means that there is a positive constant C such that \(a \le C b\).
B-spline basis functions in one dimension. The set B of basis functions with non-empty support in \(\Omega \) are indicated in deep purple. Note that basis functions crossing the boundary of \(\Omega \) are defined analogously to interior basis functions. a \(C^1 Q^2(\mathbb {R})\). b \(C^2 Q^3(\mathbb {R})\)
The B-spline spaces
Let \(\widetilde{\mathcal {T}}_{h}\), \(h \in (0,h_0]\), be a family of uniform tensor product meshes in \(\mathbb {R}^d\) with mesh parameter h.
Let \(\widetilde{V}_{h}=C^{p-1}Q^p(\mathbb {R}^d)\) be the space of \(C^{p-1}\) tensor product B-splines of order p defined on \(\widetilde{\mathcal {T}}_{h}\). Let \(\widetilde{B} = \{\varphi _i\}_{ i\in \widetilde{I}}\) be the standard basis in \(\widetilde{V}_h\), where \(\widetilde{I}\) is an index set.
Let \(B = \{ \varphi \in \widetilde{B} \, : \, {{\mathrm{supp}}}(\varphi ) \cap \Omega \ne \emptyset \}\) be the set of basis functions with support that intersects \(\Omega \). Let I be an index set for B. Let \(V_h = {{\mathrm{span}}}\{B\}\) and let \(\mathcal {T}_h = \{ T \in \widetilde{\mathcal {T}}_h : T \subset \cup _{\varphi \in B } {{\mathrm{supp}}}(\varphi ) \}\). An illustration of the basis functions in 1D is given in Fig. 1.
Let \(B = B_a \cup 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 \cup I_r\) be the corresponding partition of the index set. Let \(V_{h,a} = {{\mathrm{span}}}\{B_a\}\) be the active finite element space.
To construct the basis functions in \(\widetilde{V}_h\) we start with the one dimensional line \(\mathbb {R}\) and define a uniform partition, with nodes \(x_i = i h\), \(i \in \mathbb {Z}\), where h is the mesh parameter, and elements \(I_i = [x_{i-1},x_i)\). We define
$$\begin{aligned} \varphi _{i,0}(x) = {\left\{ \begin{array}{ll} 1 &{} x \in I_i \\ 0 &{} x \in \mathbb {R}\setminus I_i \end{array}\right. } \end{aligned}$$
The basis functions \(\varphi _{i,p}\) are then defined by the Cox-de Boor recursion formula
$$\begin{aligned} \varphi _{i,p} = \frac{x - x_i}{x_{i+p} - x_i} \varphi _{i,p-1}(x) + \frac{x_{i+p+1} - x}{x_{i+p+1} - x_{i+1}} \varphi _{i+1,p-1}(x) \end{aligned}$$
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 Fig. 1. We then define tensor product basis functions in \(\mathbb {R}^d\) of the form
$$\begin{aligned} \varphi _{i_1,\dots ,i_d} (x) = \prod _{k=1}^d \varphi _{i_k}(x_k) \end{aligned}$$
The nonsymmetric method
Find \(u_{h,a} \in V_{h,a}\) such that
$$\begin{aligned} A_{h}(u_{h,a},v) = L_h(v)\qquad v \in V_{h,a} \end{aligned}$$
The forms are defined by
$$\begin{aligned} A_h(v,w)&= a_h(v,w) + \tau h^2(\Delta v,\Delta w)_{\mathcal {T}_{h,D} \cap \Omega } \end{aligned}$$
$$\begin{aligned} L_h(v)&= l_h(v) + \tau h^2 (f,\Delta v)_{\mathcal {T}_{h,D}\cap \Omega } \end{aligned}$$
$$\begin{aligned} a_h(v,w)&= (\nabla v,\nabla w)_\Omega - (n\cdot \nabla v,w)_{\partial \Omega _D} + (v, n\cdot \nabla w)_{\partial \Omega _D} + \beta h^{-1} (v,w)_{\partial \Omega _D} \end{aligned}$$
$$\begin{aligned} l_h(v)&= (f, v)_\Omega + (g_N,v)_{\partial \Omega _N} + (g_D,n\cdot \nabla v)_{\partial \Omega _D} + \beta h^{-1}(g_D,v)_{\partial \Omega _D} \end{aligned}$$
with positive parameters \(\beta \) and \(\tau \). Furthermore, we used the notation
$$\begin{aligned} (v,w)_{\mathcal {T}_{h,D} \cap \Omega } = \sum _{T \in \mathcal {T}_{h,D}} (v,w)_{T\cap \Omega } \end{aligned}$$
where \(\mathcal {T}_{h,D} \subset \mathcal {T}_h\) is defined by
$$\begin{aligned} \mathcal {T}_{h,D} = \mathcal {T}_h\left( U_\delta (\partial \Omega _D)\right) = \{T \in \mathcal {T}_h : T \cap U_\delta (\partial \Omega _D) \ne \emptyset \} \end{aligned}$$
$$\begin{aligned} U_\delta (\partial \Omega _D) = \left( \bigcup _{x\in \partial \Omega _D} B_\delta (x)\right) \cap \Omega \end{aligned}$$
with \(\delta \sim h\) and \(B_\delta (x)\) the open ball with center x and radius \(\delta \). We note that it follows from (14) that \(U_\delta (\partial \Omega _D) \subset \mathcal {T}_{h,D}\).
Galerkin orthogonality It holds
$$\begin{aligned} A_h( u - u_h, v) = 0 \qquad \forall v \in V_h \end{aligned}$$
In practice, \(\mathcal {T}_{h,D}\) may be taken as the set of all elements that intersect the Dirichlet boundary \(\partial \Omega _D\) and their neighbors, i.e. \(\mathcal {T}_{h,D} = \mathcal {N}_h(\mathcal {T}_h(\partial \Omega _D))\).
(The Symmetric Method) The symmetric version of (8) takes the form: find \(u_{h,a} \in V_{h,a}\) such that
$$\begin{aligned} a_{h,\mathrm {sym}}(u_{h,a},v) + s_{h,\mathrm {sym}}(u_{h,a},v) = l_{h,\mathrm {sym}}(v)\qquad v \in V_{h,a} \end{aligned}$$
$$\begin{aligned} a_{h,\mathrm {sym}}(v,w)&= (\nabla v,\nabla w)_\Omega - (n\cdot \nabla v,w)_{\partial \Omega _D} - (v, n\cdot \nabla w)_{\partial \Omega _D} + \beta h^{-1} (v,w)_{\partial \Omega _D} \end{aligned}$$
$$\begin{aligned} s_{h,\mathrm {sym}}(v,w)&= \gamma h^{2p - 1} \left( [D^p_{n_F} v ], [D^p_{n_F} w ]\right) _{\mathcal {F}_{D,h}} \end{aligned}$$
$$\begin{aligned} l_{h,\mathrm {sym}}(v)&= (f,v)_\Omega + (g_N,v)_{\partial \Omega _N} - (g_D,n\cdot \nabla v)_{\partial \Omega _D} + \beta h^{-1}(g_D,v)_{\partial \Omega _D} \end{aligned}$$
where \(\beta \) and \(\gamma \) are positive parameters, \(\mathcal {F}_{h,D}\) is the set of interior faces which belong to an element in \(\mathcal {T}_h(\partial \Omega _D)\), and \(D_{n_F} = n_F \cdot \nabla \) is the directional derivative normal to the face F.
The stabilization term \(s_{h,\mathrm {sym}}\) provides the control
$$\begin{aligned} \Vert \nabla v \Vert ^2_{\mathcal {T}_h(\partial \Omega _D)} \lesssim \Vert \nabla v \Vert ^2_{\Omega } + \Vert v \Vert ^2_{s_{h,\mathrm {sym}}} \qquad v \in V_{h} \end{aligned}$$
where we note that we indeed obtain control on the full elements \(T \in \mathcal {T}_{h}(\partial \Omega _D)\). The control (21) is employed in the proof of the coercivity of \(A_h\) in the symmetric case. More precisely, (21) is used as follows
$$\begin{aligned} h \Vert n \cdot \nabla v \Vert ^2_{\partial \Omega _D} \lesssim \Vert \nabla v \Vert ^2_{\mathcal {T}_h(\partial \Omega _D)} \lesssim \Vert \nabla v \Vert ^2_{\Omega } + \Vert v \Vert ^2_{s_{h,\mathrm {sym}}} \end{aligned}$$
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.
Error estimates
Basic properties of \(\varvec{A}_{\varvec{h}}\)
Energy norm Define the norms
$$\begin{aligned} |||v |||_h^2&= \Vert \nabla v \Vert ^2_{\Omega } + h^{-1}\Vert v \Vert ^2_{\partial \Omega _D} + \tau h^2 \Vert \Delta v \Vert ^2_{\mathcal {T}_{h,D} \cap \Omega } \end{aligned}$$
$$\begin{aligned} |||v |||_{h,\bigstar }^2&= \Vert \nabla v \Vert ^2_{\Omega } + h^{-1}\Vert v \Vert ^2_{\partial \Omega _D} + \tau h^2 \Vert \Delta v \Vert ^2_{\mathcal {T}_{h,D} \cap \Omega } + h \Vert n \cdot \nabla v \Vert ^2_{\partial \Omega _D} \end{aligned}$$
Coercivity For \(\beta >0\) the form \(A_h\) is coercive
$$\begin{aligned} |||v |||_{h}^2 \lesssim A_h(v,v) \qquad v \in V + V_h \end{aligned}$$
where \(V = H^2(\Omega )\). This result follows directly from the definition and the fact that the parameters \(\tau \ge 0\) and \(\beta >0\).
Continuity The form \(A_h\) is continuous
$$\begin{aligned} A_h(v,w) \lesssim |||v |||_h |||w |||_{h,\bigstar } \qquad v,w \in V + V_h \end{aligned}$$
First we note that
$$\begin{aligned} A_h(v,w)&= (\nabla v, \nabla w)_\Omega - (n\cdot \nabla v, w)_{\partial \Omega _D} \nonumber \\&\qquad + (v, n \cdot \nabla w )_{\partial \Omega _D} + \beta h^{-1} (v, w )_{\partial \Omega _D} + \tau h^2 (\Delta v, \Delta w )_{\mathcal {T}_{h,D} \cap \Omega } \end{aligned}$$
$$\begin{aligned}&\lesssim |(\nabla v, \nabla w)_\Omega - (n\cdot \nabla v, w)_{\partial \Omega _D}| + |||v|||_h |||w |||_{h,\bigstar } \end{aligned}$$
We proceed with an estimate of the first term on the right hand side. To that end let \(\chi :\Omega \rightarrow [0,1]\) be a smooth function such that
$$\begin{aligned} {\left\{ \begin{array}{ll} \chi = 1 \quad \text {on } \partial \Omega _D\\ {{\mathrm{supp}}}(\chi ) \subset U_\delta (\partial \Omega _D )\\ \Vert \nabla \chi \Vert _{L^\infty (U_\delta (\partial \Omega _D ))} \lesssim \delta ^{-1} \end{array}\right. } \end{aligned}$$
where \(U_\delta (\partial \Omega _\delta )\) is defined in (15). Splitting the term \((\nabla v, \nabla w)_\Omega \) using \(\chi \) and then applying Green's formula for the term in the vicinity of \(\partial \Omega _D\) followed by some obvious bounds give
$$\begin{aligned}&(\nabla v, \nabla w )_{\Omega } - (n\cdot \nabla v, w)_{\partial \Omega _D} \nonumber \\&\qquad = \left( \nabla v, (1 - \chi ) \nabla w \right) _{\Omega } + (\nabla v, \chi \nabla w )_{\Omega } - (n\cdot \nabla v, \chi w)_{\partial \Omega _D} \end{aligned}$$
$$\begin{aligned}&\qquad = \left( \nabla v, (1 - \chi ) \nabla w \right) _{\Omega } - (\nabla \cdot (\chi \nabla v) , w )_{\Omega }\end{aligned}$$
$$\begin{aligned}&\qquad = \left( \nabla v, (1 - \chi ) \nabla w \right) _{\Omega } - (\nabla \chi \cdot \nabla v , w )_{\Omega } - (\chi \Delta v , w )_{\Omega }\end{aligned}$$
$$\begin{aligned}&\qquad \lesssim \Vert \nabla v\Vert _\Omega \Vert \nabla w \Vert _{\Omega } + \delta ^{-1} \Vert \nabla v\Vert _{U_\delta (\partial \Omega _D)} \Vert w \Vert _{U_\delta (\partial \Omega _D)} + \Vert \Delta v\Vert _{U_\delta (\partial \Omega _D)} \Vert w \Vert _{U_\delta (\partial \Omega _D)}\qquad \end{aligned}$$
where we in (31) assume \(v \in H^2(\Omega )\) which holds if \(V_h\) is a space of \(C^{p-1}\) tensor product B-splines of order \(p \ge 2\).
Next using the bound
$$\begin{aligned} \Vert w \Vert ^2_{U_\delta (\partial \Omega _D)} \lesssim \delta \Vert w \Vert ^2_{\partial \Omega _D} + \delta ^2 \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} \end{aligned}$$
see [5], we conclude that
$$\begin{aligned}&(\nabla v, \nabla w )_{\Omega } - (n\cdot \nabla v, w)_{\partial \Omega _D}\end{aligned}$$
$$\begin{aligned}&\qquad \lesssim \Vert \nabla v\Vert _\Omega \Vert \nabla w \Vert _{\Omega } +\delta ^{-1} \Vert \nabla v\Vert _{U_\delta (\partial \Omega _D)} ( \delta \Vert w \Vert ^2_{\partial \Omega _D} + \delta ^2 \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} )^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \qquad + \Vert \Delta v\Vert _{U_\delta (\partial \Omega _D)} ( \delta \Vert w \Vert _{\partial \Omega _D} + \delta ^2 \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} )^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \lesssim \Vert \nabla v\Vert _\Omega \Vert \nabla w \Vert _{\Omega } +\Vert \nabla v\Vert _{U_\delta (\partial \Omega _D)} ( \delta ^{-1}\Vert w \Vert ^2_{\partial \Omega _D} + \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} )^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \qquad + \delta \Vert \Delta v\Vert _{U_\delta (\partial \Omega _D)} ( \delta ^{-1} \Vert w \Vert _{\partial \Omega _D} + \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} )^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \lesssim (\Vert \nabla v\Vert ^2_\Omega + \Vert \nabla v\Vert ^2_{U_\delta (\partial \Omega _D)} + \delta ^2 \Vert \Delta v\Vert ^2_{U_\delta (\partial \Omega _D)})^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \qquad \times (\Vert \nabla w \Vert ^2_{\Omega } + \delta ^{-1} \Vert w \Vert ^2_{\partial \Omega _D} + \Vert \nabla w \Vert ^2_{U_\delta (\partial \Omega _D)} )^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \lesssim (\Vert \nabla v\Vert ^2_\Omega + \delta ^2 \Vert \Delta v\Vert ^2_{U_\delta (\partial \Omega _D)})^{1/2}\end{aligned}$$
$$\begin{aligned}&\qquad \qquad \times (\Vert \nabla w \Vert ^2_{\Omega } + \delta ^{-1} \Vert w \Vert ^2_{\partial \Omega _D} )^{1/2} \end{aligned}$$
Finally, choosing \(\delta \sim h\) and using the fact that \(U_\delta (\partial \Omega _D) \subset \mathcal {T}_{h,D}\) we obtain
$$\begin{aligned} (\nabla v, \nabla w )_{\Omega } - (n\cdot \nabla v, w)_{\partial \Omega _D} \lesssim |||v |||_h |||w |||_{h,\bigstar } \end{aligned}$$
which in combination with (28) concludes the proof. \(\square \)
Interpolation error estimates
Definition of the interpolant There is an extension operator \(E:W^k_q(\Omega ) \rightarrow W^k_q(\mathbb {R}^d)\), \(k\ge 0\) and \(q\ge 1\), such that
$$\begin{aligned} \Vert E v \Vert _{W^k_q(\mathbb {R}^d)} \lesssim \Vert v \Vert _{W^k_q(\Omega )} \end{aligned}$$
see [9]. Define the interpolant by
$$\begin{aligned} \pi _{h} : H^s(\Omega ) \ni u \mapsto \pi _{Cl,h} ( E u ) \in V_h \end{aligned}$$
where \(\pi _{Cl,h}\) is a Clement type interpolation operator onto the spline space. We have the expansion
$$\begin{aligned} \pi _h ( E v ) = \sum _{\varphi _i \in I} (\pi _h ( E v ))_i \varphi _i \end{aligned}$$
where \((\pi _h ( E v ))_i\) is the coefficient corresponding to basis function \(\varphi _i\). We define the interpolant on the active and removed finite element spaces by
$$\begin{aligned} \pi _{h,a} v = \sum _{\varphi _i \in I_a} (\pi _h (E v) )_i \varphi _i \end{aligned}$$
$$\begin{aligned} \pi _{h,r} v = \sum _{\varphi _i \in I_r} (\pi _h (E v) )_i \varphi _i \end{aligned}$$
We then have
$$\begin{aligned} \pi _h (E v) = \pi _{h,a} ( E v ) + \pi _{h,r} (E v ) \end{aligned}$$
Below we simplify the notation and write \(v = Ev\) and \(\pi _h ( E v ) = \pi _h v\).
Basis function removal condition Let \(B_r\), with corresponding index set \(I_r\), be such that
$$\begin{aligned} \sum _{i \in I_r} |||\varphi _i |||^2_{h,\bigstar } \lesssim tol^2 \end{aligned}$$
Selection procedure To determine \(B_r\) we may thus compute \(|||\varphi _i |||_{h,\bigstar }\), \(i \in I\), sort the basis functions in increasing order and then simply add functions to \(I_r\) as long as (51) is satisfied. If we wish to avoid computing \(|||\varphi _i |||_{h,\bigstar }\) we may use the directly available diagonal values \(A_h(\varphi _i,\varphi _i)\) of the stiffness matrix as approximations.
Lemma 1
(Interpolation error estimate) Let \(\pi _{h,a}\) be defined by (48) with \(B = B_a \cup B_r\) such that \(B_r\) satisfies (51), then
$$\begin{aligned} |||v - \pi _{h,a} v |||_{h,\bigstar }&\lesssim (h^{p} + tol) \Vert v \Vert _{H^{p+1}(\Omega )} \end{aligned}$$
Using the identity \(\pi _h v = \pi _{h,a} v + \pi _{h,r} v\) and the triangle inequality
$$\begin{aligned} |||v - \pi _{h,a} v |||^2_{h,\bigstar }&\lesssim |||v - \pi _{h} v |||^2_{h,\bigstar } + |||\pi _{h,r} v |||^2_{h,\bigstar } \end{aligned}$$
$$\begin{aligned}&\lesssim h^{2p} \Vert v \Vert ^2_{H^{p+1}(\Omega )} + |||\pi _{h,r} v |||^2_{h,\bigstar } \end{aligned}$$
by standard spline interpolation results [2]. To estimate the second term on the right hand side we introduce the scalar product
$$\begin{aligned} \langle v,w \rangle _{h,\bigstar } = (\nabla v, \nabla w)_\Omega \,+\, h (n\cdot \nabla v, n\cdot \nabla w)_{\partial \Omega _D} \,+\, h^{-1} (v, w)_{\partial \Omega _D}\,+ \,h^2(\Delta v,\Delta w)_{\mathcal {T}_{h,D} \cap \Omega } \end{aligned}$$
associated with the norm \(|||\cdot |||_{h,\bigstar }\). Expanding \(\pi _{h,r} v\) in the basis \(B_r\) we get
$$\begin{aligned} |||\pi _{h,r} v |||^2_{h,\bigstar }&=\sum _{i,j\in I_r} (\pi _h v)_i (\pi _h v)_j \langle \varphi _i,\varphi _j \rangle _{h,\bigstar } \end{aligned}$$
$$\begin{aligned}&\le \sum _{i \in I_r} \sum _{j\in I_r} \delta _{ij} |(\pi _h v)_i|\, |(\pi _h v)_j| \, |||\varphi _i |||_{h,\bigstar } |||\varphi _j |||_{h,\bigstar }\end{aligned}$$
$$\begin{aligned}&\le \sum _{i \in I_r} \sum _{j\in I_r} \frac{\delta _{ij}}{2} |(\pi _h v)_i|^2 |||\varphi _i |||_{h,\bigstar } + \frac{\delta _{ij}}{2} |(\pi _h v)_j|^2 |||\varphi _j |||_{h,\bigstar }\end{aligned}$$
$$\begin{aligned}&= \sum _{i \in I_r} \left( \sum _{j\in I_r} \delta _{ij} \right) |(\pi _h v)_i|^2 |||\varphi _i |||^2_{h,\bigstar }\end{aligned}$$
$$\begin{aligned}&\lesssim \Vert \pi _h v \Vert ^2_{L^{\infty }(\mathcal {N}_h(\Omega ))} \left( \sum _{i \in I_r} |||\varphi _i |||^2_{h,\bigstar } \right) \end{aligned}$$
$$\begin{aligned}&\lesssim \Vert v \Vert ^2_{H^{p+1}(\Omega )} tol^2 \end{aligned}$$
We defined
$$\begin{aligned} \delta _{ij} = {\left\{ \begin{array}{ll} 1 &{} \text {if }{{\mathrm{supp}}}(\varphi _i) \cap {{\mathrm{supp}}}(\varphi _j) \ne \emptyset \\ 0 &{} \text {if } {{\mathrm{supp}}}(\varphi _i) \cap {{\mathrm{supp}}}(\varphi _j) = \emptyset \end{array}\right. } \end{aligned}$$
and we have the bound
$$\begin{aligned} \sum _{j \in I_r} \delta _{ij} \le (2p+1)^d \end{aligned}$$
We used the \(L^\infty (\mathcal {N}_h(\Omega ))\) stability of the interpolant \(\pi _h\) and then the \(L^\infty \) stability of the extension operator and finally the Sobolev embedding theorem
$$\begin{aligned} \Vert \pi _h v \Vert _{L^{\infty }(\mathcal {N}_h(\Omega ))} \lesssim \Vert v \Vert _{L^{\infty }(\mathcal {N}_h(\Omega ))} \lesssim \Vert v \Vert _{L^{\infty }(\Omega )} \lesssim \Vert v \Vert _{H^{p+1}(\Omega )} \end{aligned}$$
\(\square \)
Error estimate
We have the following error estimate.
Theorem 1
Let \(u_{h,a}\) be the solution to (8) with \(V_{h,a} = {{\mathrm{span}}}\{B_a\}\) the active spline space, \(V_h = {{\mathrm{span}}}\{B\}\) the full spline space, and \(B = B_a \cup B_r\), where \(B_r\) satisfies (51) with \(tol\sim h^p\), then
$$\begin{aligned} |||u - u_{h,a} |||_h \lesssim h^{p} \Vert u \Vert _{H^{p+1}(\Omega )} \end{aligned}$$
Using coercivity (25), Galerkin orthogonality (16), and continuity (26), we obtain
$$\begin{aligned} |||u - u_{h,a} |||^2_h&\lesssim A_h( u - u_{h,a}, u - u_{h,a} ) \end{aligned}$$
$$\begin{aligned}&= A_h( u - u_{h,a}, u - \pi _{h,a} u ) \end{aligned}$$
$$\begin{aligned}&\lesssim |||u - u_{h,a} |||_{h} |||u - \pi _{h,a} u |||_{h,\bigstar } \end{aligned}$$
Thus we arrive at
$$\begin{aligned} |||u - u_{h,a} |||_h \lesssim |||u - \pi _{h,a} u |||_{h,\bigstar } \end{aligned}$$
which together with the interpolation error estimate (52) completes the proof of (65). \(\square \)
Note that if we take \(\tau =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
$$\begin{aligned} \Vert \nabla ( u - u_{h,a}) \Vert ^2_\Omega + \Vert u - u_h \Vert ^2_{\partial \Omega _D}&\lesssim A_h(u - u_{h,a}, u - u_{h,a} ) \end{aligned}$$
$$\begin{aligned}&= A_h( u - u_{h,a}, u - \pi _{h,a} u )\end{aligned}$$
Now we note that
$$\begin{aligned} |||u - u_{h,a} |||^2_{h}&= \Vert \nabla ( u - u_{h,a}) \Vert ^2_\Omega + \Vert u - u_h \Vert ^2_{\partial \Omega _D} + h^2 \Vert \Delta ( u - u_{h,a} ) \Vert ^2_{\mathcal {T}_{h,D}\cap \Omega }\end{aligned}$$
$$\begin{aligned}&=\Vert \nabla ( u - u_{h,a}) \Vert ^2_\Omega + \Vert u - u_h \Vert ^2_{\partial \Omega _D} + h^2 \Vert f - \Delta u_{h,a} \Vert ^2_{\mathcal {T}_{h,D}\cap \Omega } \end{aligned}$$
and thus we obtain the bound
$$\begin{aligned} \Vert \nabla ( u - u_{h,a}) \Vert ^2_\Omega + \Vert u - u_h \Vert ^2_{\partial \Omega _D}&\lesssim h^{2p} \Vert u \Vert ^2_{H^{p+1}(\Omega )} + h^2 \Vert f - \Delta u_{h,a} \Vert ^2_{\mathcal {T}_{h,D}\cap \Omega } \end{aligned}$$
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 \(\Omega \) 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 (51). 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
$$\begin{aligned} \sum _{i \in I_r} |||\varphi _i |||^2_{h,\bigstar } \lesssim h^{-(d-1)} \max _{i \in I_r } |||\varphi _i |||^2_{h,\bigstar } \end{aligned}$$
and setting \(tol\sim h^p\) we get
$$\begin{aligned} \max _{i \in I_r } |||\varphi _i |||^2_{h,\bigstar } \lesssim h^{d-1} tol \lesssim h^{2p + d -1} \end{aligned}$$
and we may define \(B_r\) as all basis functions \(\varphi \in B\) such that
$$\begin{aligned} |||\varphi |||^2_{h,\bigstar } \lesssim h^{d-1} tol \lesssim h^{2p + d -1} \end{aligned}$$
Let us for simplicity consider a basis function \(\varphi \) such that \({{\mathrm{supp}}}(\varphi ) \subset \partial \Omega _D = \emptyset \), i.e. a basis function that reside on the Neumann part of the boundary. In this case \( |||\varphi |||_{h,\bigstar } = \Vert \nabla \varphi \Vert _{{{\mathrm{supp}}}(\varphi ) \cap \Omega } \) and thus \(\varphi \in B_r\) if
$$\begin{aligned} \Vert \nabla \varphi \Vert ^2_{{{\mathrm{supp}}}(\varphi ) \cap \Omega } \lesssim h^{2p + d -1} \end{aligned}$$
The 1D case: energy norm Let \(\Omega = [0,1]\) and consider a basis function \(\varphi \) with support \([X_0, X_1]\) such that \(X_0 \in [0,1]\) and \({{\mathrm{supp}}}(\varphi ) \cap [0,1] = [X_0,1]\) is an interval of length \(\delta \). Then for \(\delta \) small enough we have
$$\begin{aligned} \varphi (x) = \left( \frac{x}{h} \right) ^p, \qquad |D \varphi (x)|^2 = \frac{p^2}{h^2}\left( \frac{x}{h} \right) ^{2(p-1)} \end{aligned}$$
up to constants and in local coordinates with origo \(X_0\), and
$$\begin{aligned} \Vert D \varphi \Vert ^2 = \int _0^\delta \frac{p^2}{h^2}\left( \frac{x}{h} \right) ^{2(p-1)} = \frac{p}{h} \frac{p}{2p-1}\left( \frac{\delta }{h}\right) ^{2p - 1} \end{aligned}$$
Condition (79) thus takes the form
$$\begin{aligned} \frac{p}{h} \frac{p}{2p-1} \left( \frac{\delta }{h}\right) ^{2p - 1} \lesssim h^{2p + d - 1} \quad \Longrightarrow \quad \frac{\delta }{h} \lesssim h^{\frac{2p + 1}{2p - 1}} \end{aligned}$$
For Lagrange basis functions we instead have \(|D\varphi (x)| \sim h^{-1}\) and we therefore obtain the condition
$$\begin{aligned} \delta h^{-2} \lesssim h^{2p + d - 1} \quad \Longrightarrow \quad \frac{x}{\delta } \lesssim h^{2p + 1} \end{aligned}$$
An illustration of both B-spline and Lagrange basis functions in this setting is given in Fig. 2. Comparing (82) and (83) we note that the condition is much stronger for the Lagrange functions and higher order p.
B-spline (top row) and Lagrange (bottom row) basis functions of order \(p=2,3\) in a 1D element intersecting \(\Omega \). Note that gradient of the blue B-spline basis functions is \(O(h^{-1}(\frac{\delta }{h})^{p-1})\) within \(\Omega \) while the gradient of Lagrange basis functions is \(O(h^{-1})\) regardless of p. a \(C^1 Q^2\) basis. b \(C^1 Q^2\) gradient. c \(C^2 Q^3\) basis. d \(C^2 Q^3\) gradient. e \(Q^2\) basis. f \(Q^2\) gradient. g \(Q^3\) basis. h \(Q^3\) gradient
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
$$\begin{aligned} \Vert D \varphi \Vert _{L^\infty ({{\mathrm{supp}}}(\varphi ) \cap \Omega )} \lesssim h^{-1} \left( \frac{\delta }{h} \right) ^{p-1} \end{aligned}$$
while for Lagrange basis functions
$$\begin{aligned} \Vert D \varphi \Vert _{L^\infty ({{\mathrm{supp}}}(\varphi ) \cap \Omega )} \lesssim h^{-1} \end{aligned}$$
which in the latter case can not be controlled by decreasing \(\delta \), see Fig. 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
$$\begin{aligned} \frac{\delta }{h} \lesssim h^{\frac{p+1}{p-1}} \end{aligned}$$
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 \({{\mathrm{supp}}}(\varphi )\) which reside in the interior of \(\Omega \). Let \(\{e_i\}_{i=1}^d\) 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\}_{i=1}^d\) of \({{\mathrm{supp}}}(\varphi )\) which originates at \(X_0\), see Fig. 3. Using the local coordinates in the vicinity of \(X_0\) we have the expansions
$$\begin{aligned} \varphi (x_1,x_2)= & {} \left( \frac{x_1}{h} \right) ^p \left( \frac{x_2}{h}\right) ^p \end{aligned}$$
$$\begin{aligned} |\nabla \varphi (x_1,x_2)|^2= & {} \frac{1}{h^2} \left( \frac{x_1}{h} \right) ^{2p-2} \left( \frac{x_2}{h}\right) ^{2p} + \frac{1}{h^2} \left( \frac{x_1}{h} \right) ^{2p} \left( \frac{x_2}{h}\right) ^{2p-2} \end{aligned}$$
Let \(\delta _i = \Vert X_i - X_0 \Vert _{\mathbb {R}^d}\) be the distance from the vertex \(X_0\) to the intersection \(X_i\) of edge \(E_i\) with the boundary \(\partial \Omega \). Assume that
$$\begin{aligned} {{\mathrm{supp}}}(\varphi ) \cap \Omega \subset [0,\delta _1] \times [0,\delta _2] \end{aligned}$$
Integrating over \([0,\delta _1] \times [0,\delta _2]\) we obtain
$$\begin{aligned} \int _0^{\delta _1} \int _0^{\delta _2} |\nabla \varphi |^2&\lesssim \left( \frac{\delta _1}{h} \right) ^{2p-1} \left( \frac{\delta _2}{h}\right) ^{2p+1} + \left( \frac{\delta _1}{h} \right) ^{2p+1} \left( \frac{\delta _2}{h}\right) ^{2p-1} \end{aligned}$$
$$\begin{aligned} \left( \frac{\delta _1}{h} \right) ^{2p-1} \left( \frac{\delta _2}{h}\right) ^{2p+1} + \left( \frac{\delta _1}{h} \right) ^{2p+1} \left( \frac{\delta _2}{h}\right) ^{2p-1} \lesssim h^{2p + d - 1} \end{aligned}$$
which implies
$$\begin{aligned} \frac{\delta _1}{h} \lesssim h \left( \frac{\delta _2}{h} \right) ^{-\frac{2p-1}{2p+1}} \qquad \text {and} \qquad \frac{\delta _2}{h} \lesssim h \left( \frac{\delta _1}{h} \right) ^{-\frac{2p-1}{2p+1}} \end{aligned}$$
See Fig. 4 for an illustration of this condition.
Illustration of the geometric quantities used in intersection conditions (51) in energy norm and (96) in max norm
Illustrations of the basis function intersection condition (51) in energy norm and (96) in max norm for splines of polynomial order \(p=1,2,\dots ,5\). a Energy norm, \(h=0.1\). b Energy norm, \(h=0.05\). c Max norm, \(h=0.1\). d Max norm, \(h=0.05\)
The 2D case: max norm Starting from the expansion (88) and observing that for small enough \(\delta \) parameters \(|\nabla \varphi |^2\) is increasing when we move out from the vertex. Using assumption (89) we thus conclude that
$$\begin{aligned} \Vert \nabla \varphi \Vert _{L^\infty ( {{\mathrm{supp}}}(\varphi ) \cap \Omega } \lesssim |\nabla \varphi (\delta _1,\delta _2) | \end{aligned}$$
$$\begin{aligned} \nabla \varphi (x_1,x_2) = \left[ \frac{p}{h} \left( \frac{x_1}{h} \right) ^{p-1} \left( \frac{x_2}{h}\right) ^p, \; \frac{p}{h} \left( \frac{x_1}{h} \right) ^{p} \left( \frac{x_2}{h}\right) ^{p-1} \right] \end{aligned}$$
Setting \(x_1 = \delta _1\) and \(x_2 = \delta _2\) we get the conditions
$$\begin{aligned} \frac{p}{h} \left( \frac{\delta _1}{h} \right) ^{p-1} \left( \frac{\delta _2}{h}\right) ^p \lesssim h^p \quad \text {and} \quad \frac{p}{h} \left( \frac{\delta _1}{h} \right) ^{p} \left( \frac{\delta _2}{h}\right) ^{p-1} \lesssim h^p \end{aligned}$$
which we may write in the form
$$\begin{aligned} \frac{\delta _1}{h} \lesssim \frac{1}{p} h^{\frac{p+1}{p}} \left( \frac{\delta _2}{h} \right) ^{-\frac{p-1}{p}} \quad \text {and} \quad \frac{\delta _2}{h}\lesssim \frac{1}{p} h^{\frac{p+1}{p}} \left( \frac{\delta _1}{h} \right) ^{-\frac{p-1}{p}} \end{aligned}$$
Numerical results
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:\Omega \rightarrow \mathbb {R}^d\) such that
$$\begin{aligned} -\sigma (u)\cdot \nabla = f \qquad {\text {in}}\,{\Omega } \end{aligned}$$
$$\begin{aligned} \sigma (u) \cdot n = g_N \qquad {\text {on}}\,{\partial \Omega _N} \end{aligned}$$
$$\begin{aligned} u = g_D \qquad {\text {on}}\,{\partial \Omega _D} \end{aligned}$$
where the stress and strain tensors are defined by
$$\begin{aligned} \sigma (u) = 2\mu \epsilon (u) + \lambda \text {tr}(\epsilon (u)), \qquad \epsilon (u) = \frac{1}{2}\Big ( u\otimes \nabla + \nabla \otimes u \Big ) \end{aligned}$$
with Lamé parameters \(\lambda \) and \(\mu \); f, \(g_N\), \(g_D\) are given data; \(a \otimes b\) is the tensor product of vectors a and b with elements \((a \otimes b)_{ij} = a_i b_j\).
The nonsymmetric method for linear elasticty
Find \(u_{h,a} \in [V_{h,a}]^d\) such that
$$\begin{aligned} A_{h}(u_{h,a},v) = L_h(v)\qquad v \in [V_{h,a}]^d \end{aligned}$$
$$\begin{aligned} A_h(v,w)&= a_h(v,w) + \tau h^2(\epsilon (v)\cdot \nabla ,\epsilon (w)\cdot \nabla )_{\mathcal {T}_{h,D} \cap \Omega }\end{aligned}$$
$$\begin{aligned} L_h(v)&= l_h(v) + \tau h^2 (f,\epsilon (v)\cdot \nabla )_{\mathcal {T}_{h,D}\cap \Omega } \end{aligned}$$
$$\begin{aligned} a_h(v,w)&= (\sigma (v),\epsilon (w))_\Omega - (\sigma (v) \cdot n ,w)_{\partial \Omega _D} + (v, \sigma (w)\cdot n )_{\partial \Omega _D} + \beta h^{-1} (v,w)_{\partial \Omega _D}\nonumber \\\end{aligned}$$
$$\begin{aligned} l_h(v)&= (f,v)_\Omega + (g_N,v)_{\partial \Omega _N} + (g_D,\sigma (v)\cdot n)_{\partial \Omega _D} + \beta h^{-1}(g_D,v)_{\partial \Omega _D} \end{aligned}$$
with positive parameters \(\beta \) and \(\tau \). Furthermore, the energy norm is defined
$$\begin{aligned} |||v |||_h^2 = ( \sigma (v) , \epsilon (v) )_\Omega + h^{-1} \Vert v \Vert _{\partial \Omega _D}^2 + \tau h^2 \Vert \epsilon (v)\cdot \nabla \Vert ^2_{\mathcal {T}_h(\partial \Omega _D) \cap \Omega } \end{aligned}$$
Geometries in the two model problems. Boundaries with non-homogeneous Neumann conditions are indicated in blue and Dirichlet boundaries are indicated in red. a Neumann problem. b Manufactured problem
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 Fig. 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 \(\nu =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
$$\begin{aligned}&\Omega = [0,1]^2 , \ \partial \Omega _D = \{x\in [0,1],y=0\} , \ \partial \Omega _N = \partial \Omega \backslash \partial \Omega _D \end{aligned}$$
$$\begin{aligned}&u(x,y) = [ -\cos (\pi x)\sin (\pi y), \, \sin (\pi x/7)\sin (\pi y/3) ]/10 \end{aligned}$$
see Fig. 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 \(\partial \Omega _D\) we in our current implementation neglect the least squares term in the vicinity of \(\partial \Omega _D\), i.e. we choose \(\tau =0\).
Illustration of the selection procedure
We utilize the selection procedure based on the stiffness matrix proposed in "Interpolation error estimates" section. Some realizations of this selection are visualized in Fig. 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 = c h^p\). This can also be seen in Fig. 7 where we investigate how the choice of this constant effects the number of removed basis functions. In Fig. 8 we note that the use of basis removal is quite effective and also gives better quality stresses along the boundary.
Four realizations of removed basis functions using on the stiffness matrix based selection procedure described in "Interpolation error estimates" section; all using the same constant \(c=0.01\) for the tolerance \(tol = ch^{p} \times \sqrt{E}\) in (51). 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. a \(C^1Q^2\), \(h=0.4\). b \(C^1Q^2\), \(h=0.2\). c \(C^1Q^2\), \(h=0.1\). d \(C^2Q^3\), \(h=0.2\)
Studies of how the choice of constant c for the tolerance \(tol=c h^p \times \sqrt{E}\) in (51) relates to the number of removed basis functions. The set-up here is the same as in Fig. 6. a a \(h=0.2\). b \(h=0.1\)
Displacements and von-Mises stresses from numerical solutions with and without basis removal in the Neumann problem using \(C^1 Q^2\)-splines and mesh size \(h=0.1\). In the detailed view we note poor quality of the stresses on the boundary in the standard solution which is remedied when removing the problematic basis function. a Standard solution. b Detail in standard solution. c Basis removal solution. d Detail in basis removal solution
To estimate the convergence we use the manufactured problem described in "Linear elasticity" section and the cut situations are induced by rotating the background grid \(\pi /7\) radians as illustrated by the mesh with removed basis functions in Fig. 9 together with the corresponding numerical solution. In Fig. 10a we present convergence studies in energy norm for various choices of the constant c in the tolerance \(tol = c h^{p} \times \sqrt{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 Fig. 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})\).
Example of numerical solution using \(C^1Q^2\) splines and mesh size \(h=0.1\). The mesh is rotated \(\pi /7\) radians to induce cut situations and the removed basis functions are selected using the tolerance \(tol=0.01 h^2 \times \sqrt{E}\). a Mesh and removed spline basis functions. b Numerical solution
Convergence in \(|||\cdot |||_h\) norm and condition numbers for the manufactured problem using basis removal with \(C^1Q^2\)-spline basis. The tolerances used in the selection procedure is \(tol = c h^p \times \sqrt{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\)). a Energy norm convergence. b Condition number
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.
Badia S, Verdugo F, Martín AF. The aggregated unfitted finite element method for elliptic problems. Sept: ArXiv e-prints; 2017.
Bazilevs Y, Beirão da Veiga L, Cottrell JA, Hughes TJR, Sangalli G. Isogeometric analysis: approximation, stability and error estimates for \(h\)-refined meshes. Math Models Methods Appl Sci. 2006;16(7):1031–90.
Burman E. Ghost penalty. C R Math Acad Sci Paris. 2010;348(21–22):1217–20.
Burman E, Claus S, Hansbo P, Larson MG, Massing A. CutFEM: discretizing geometry and partial differential equations. Int J Numer Methods Eng. 2015;104(7):472–501.
Burman E, Hansbo P, Larson MG. A cut finite element method with boundary value correction. Math Comput. 2018;87:633–57.
Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric anlysis: toward integration of CAD and FEA. Chichester: John Wiley & Sons, Ltd.; 2009.
Dauge M, Düster A, Rank E. Theoretical and numerical investigation of the finite cell method. J Sci Comput. 2015;65(3):1039–64.
Embar A, Dolbow J, Harari I. Imposing Dirichlet boundary conditions with Nitsche's method and spline-based finite elements. Int J Numer Methods Eng. 2010;83(7):877–98.
Folland GB. Introduction to partial differential equations. 2nd ed. Princeton: Princeton University Press; 1995.
Hansbo P, Larson MG, Larsson K. Cut finite element methods for linear elasticity problems. In: Bordas S, Burman E, Larson M, Olshanskii M, (eds), In: Proceedings of the UCL Workshop 2016: geometrically unfitted finite element methods and applications. Berlin: Springer; 2018. To be published.
Johansson A, Larson MG. A high order discontinuous Galerkin Nitsche method for elliptic problems with fictitious boundary. Numer Math. 2013;123(4):607–28.
Parvizian J, Düster A, Rank E. Finite cell method: \(h\)- and \(p\)-extension for embedded domain problems in solid mechanics. Comput Mech. 2007;41(1):121–33.
All authors have prepared the manuscript. All authors read and approved the final manuscript.
This research was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants Nos. 2013-4708, 2017-03911, and the Swedish Research Programme Essence.
Department of Mathematics and Mathematical Statistics, Umeå University, 90 187, Umeå, Sweden
Daniel Elfverson
, Mats G. Larson
& Karl Larsson
Search for Daniel Elfverson in:
Search for Mats G. Larson in:
Search for Karl Larsson in:
Correspondence to Karl Larsson.
Elfverson, D., Larson, M.G. & Larsson, K. CutIGA with basis function removal. Adv. Model. and Simul. in Eng. Sci. 5, 6 (2018) doi:10.1186/s40323-018-0099-2
Unfitted techniques in computational mechanics