---
abstract: 'We propose a Nitsche-based fictitious domain method for the three field Stokes problem in which the boundary of the domain is allowed to cross through the elements of a fixed background mesh.  The dependent variables of velocity, pressure and extra-stress tensor are discretised on the background mesh using linear finite elements.  This equal order approximation is stabilized using a continuous interior penalty (CIP) method.  On the unfitted domain boundary, Dirichlet boundary conditions are weakly enforced using Nitsche’s method.  We add CIP-like ghost penalties in the boundary region and prove that our scheme is inf-sup stable and that it has optimal convergence properties independent of how the domain boundary intersects the mesh.  Additionally, we demonstrate that the condition number of the system matrix is bounded independently of the boundary location.  We corroborate our theoretical findings numerically.' 
author:
- Erik Burman
- Susanne Claus
- André Massing
title: '**A Stabilized Cut Finite Element Method for the Three Field Stokes Problem**'
---

Three field Stokes, continuous interior penalty, fictitious domain, cut finite element method, ghost penalty, Nitsche’s method, viscoelasticity

65N12, 65N15, 65N30, 76D07

Introduction
============

In this article, we develop a stabilized finite element method for the so-called three field Stokes system [@BonvinPicassoStenberg2001; @Bonito2006; @Bonito2008].  In the three field Stokes equation, the extra-stress tensor is considered as a separate variable, additionally to velocity and pressure.  This description of the Stokes system is of particular interest in viscoelastic fluid mechanics, where the extra-stress tensor is related to the rate of deformation tensor through a non-linear constitutive equation and therefore can no longer be easily substituted into the momentum equation [@Owens2002].\
In particular, we develop a finite element scheme in which the surface of the fluid can cut elements in the computational mesh in an arbitrary manner.  Such cut finite element methods are especially beneficial for applications in which the use of interface tracking techniques such as arbitrary Lagrangian Eulerian methods [@DoneaALE], where the mesh is fitted and moved with the interface, involve frequent re-meshing and sophisticated mesh moving algorithms that can be prohibitively expensive.  One such application of high interest is the simulation of viscoelastic free surface flows in which the fluid surface undergoes large deformations and in which drop detachments may occur.  This type of free surface flow of viscoelastic liquids plays a key role in a wide range of industrial applications such as surface coating for molten plastics, filtration operations of engine oils or inkjet printing.\
Our cut finite element method is based on an earlier formulation for fitted meshes presented in [@Bonito2008; @Bonito2006], where equal order approximation spaces for all variables are combined with a continuous interior penalty method to obtain a stable and optimally convergent method for the three field Stokes equation.  We employ Nitsche’s method [@Nitsche1971] to weakly enforce the boundary conditions on the unfitted boundary domain.  Nitsche-type methods for unfitted interface problems and fictitious domain methods have previously been developed in [@HansboHansbo2002; @HansboHansboLarson2003; @DoHa09; @HaDo10] for elliptic problems, in [@HaHa04; @BeckerBurmanHansbo2009] for elasticity problems and in [@MassingLarsonLoggEtAl2013; @MassingLarsonLoggEtAl2013a; @BurmanHansbo2013] for Stokes problems.  A particular complication here is the potential occurrence of elements which are only partially covered by the physical domain $\Omega$. In such a case it has been demonstrated in [@BurmanHansbo2012; @Burman2010] that the sole application of Nitsche’s method results in suboptimal schemes, where the discretization error and the condition number of the discrete system are highly dependent on the position of the boundary with respect to the mesh.  In our contribution, we therefore apply so-called ghost penalties [@BurmanHansbo2012; @Burman2010] in the vicinity of the boundary to extend the solution of velocity, pressure and extra-stress from the cut part of the element partially covered by the physical domain to the whole element in the interface zone resulting in a fictitious domain approach.  The ghost penalties consist of penalties on the gradient jumps of the velocity, pressure and extra-stress and are applied to all edges in the interface cell layer.  Due to their similarity to the continuous interior penalty terms, our scheme allows to cope with both inf-sup and fictitious domain related instabilities in a unified way.  The resulting formulation is weakly consistent and we prove that our scheme satisfies a uniform inf-sup condition and exhibits optimal convergence order independent of how the boundary cuts the mesh.  Numerical experiments demonstrate that the resulting discrete system is insensitive to the position of the boundary within the computational domain. 

The remainder of this paper is organised as follows.  Firstly, we briefly review the strong and weak formulation of the three field Stokes system in Section \[sec:tfstokes\].  In Section \[sec:fem-formulation\], the novel cut finite element method for the three field Stokes problem is introduced.  Section \[sec:approximation-properties\] summarizes some useful inequalities and provides certain interpolation estimates which are necessary to proceed with the stability analysis in Section \[sec:stability\] and to establish the a priori estimates in Section \[sec:apriori-analysis\].  Finally, we present numerical results corroborating the theoretical findings along with numerical investigations of the conditioning of the discrete system.  Moreover, we demonstrate the applicability of our discretization method to complex 3D geometries. 

The three field Stokes problem {#sec:tfstokes}
==============================

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ ($d=2$ or $3$) with a Lipschitz boundary $\Gamma = \partial \Omega$. The three field Stokes problem reads: Find the extra-stress tensor ${\sigma}: \Omega
\rightarrow
  \mathbb{R}^{d \times d}$, the velocity ${{u}}: \Omega \rightarrow
  \mathbb{R}^{d}$ and the pressure $p: \Omega \rightarrow
  \mathbb{R}$ such that $$\left\{
  \begin{aligned}
      {\sigma}- 2 \eta {{\epsilon}}({{u}}) &= 0 \mbox{ in } \Omega,
    \\
      - \nabla \cdot {\sigma}+ \nabla p &= {{f}}\mbox{ in } \Omega,
    \\
      \nabla \cdot {{u}}&= 0 \mbox{ in } \Omega, 
    \\
      {{u}}&= {{g}}\mbox{ on } \partial \Omega. 
  \end{aligned}
    \label{eq:tfstokes-strong}
    \right.$$ Here, ${{\epsilon}}({{u}}) = \frac{1}{2} \left(\nabla {{u}}+ \nabla
  {{u}}^{\top}\right)$ is the rate of deformation tensor, ${{f}}:\Omega \rightarrow \mathbb{R}^d$ is the body force, $\eta$ is the fluid viscosity and $ {{g}}$ is the Dirichlet boundary value.  To be compatible with the divergence constraint in , the boundary data is supposed to satisfy $\int_{\partial \Omega} {{n}}\cdot {{g}}{\,\mathrm{d}s}= 0$ where ${{n}}$ denotes the outward pointing boundary normal.  Note that the only difference between the Stokes equation and the three field Stokes equation is that the extra-stress tensor is kept as a separate variable.  This type of equation system is of particular interest in viscoelastic fluid mechanics, where the stress tensor depends on the rate of deformation tensor through a non-linear constitutive equation.  Hence, the extra-stress tensor can no longer be substituted into the momentum equation.\
The weak formulation of the three field Stokes system is obtained by multiplying with test functions $({{{ \tau}}},{{v}},p) \in
\left[L^2(\Omega)\right]^{d \times d} \times \left[H^1_0(\Omega)\right]^{d} \times L^2_0(\Omega)$, by integrating over $\Omega$, and by integrating by parts.  The resulting weak formulation reads: Find $({\sigma}, {{u}}, p) \in \left[L^2(\Omega)\right]^{d \times d }
\times\left[H^1_0(\Omega)\right]^{d} \times
L^2_0(\Omega) $ such that $$\begin{aligned}
\frac{1}{2\eta}
\left({\sigma}, \, {{{ \tau}}}\right)_{\Omega} - \left( {{\epsilon}}({{u}}), \,
{{{ \tau}}}\right)_{\Omega} &=& 0 \nonumber,
  \\
\left({\sigma}, \, {{\epsilon}}({{v}}) \right)_{\Omega} - \left(p, \, \nabla
\cdot {{v}}\right)_{\Omega} &=& \left({{f}}, \, {{v}}\right)_{\Omega}
  \nonumber, \\
\left(\nabla \cdot {{u}}, \, q \right)_{\Omega} &=&  0 \nonumber
  \\
\label{eq:tfstokes-weak}\end{aligned}$$ for all $({{{ \tau}}}, {{v}}, q) \in \left[L^2(\Omega)\right]^{d \times d}
\times \left[H^1_g(\Omega)\right]^{d} \times L^2_0(\Omega)$. Here and throughout this work, we use the notation $H^s(U)$ and $[H^s(U)]^d$ for the standard Sobolev space of order $s \in {\mathbb{R}}$ and their ${\mathbb{R}}^d$-valued equivalents defined on the (possibly lower-dimensional) domain $U \subseteq {\mathbb{R}}^d$. The associated inner products and norms are written as $(\cdot, \cdot)_{s, U}$ and $\| \cdot \|_{s, U}$. If $s = 0$, we usually drop the index $s$ if no ambiguities occur.  For $s > 1/2$, we let $[H_g^s(\Omega)]^d$ consist of all functions in $[H(\Omega)]^d$ whose boundary traces are equal to $g$. Finally, $L^2_0(\Omega)$ denotes the functions in $L^2(\Omega)$ with zero average over $\Omega$.

Stabilized Cut Finite Element Formulation {#sec:fem-formulation}
=========================================

Let $\Omega$ be an open and bounded domain in ${\mathbb{R}}^d$ ($d = 2, 3$) with Lipschitz boundary $\Gamma = \partial \Omega$ and let ${\mathcal{T}_h}$ be a quasi-uniform tesselation that covers the domain $\Omega$. We do not assume that the mesh ${\mathcal{T}_h}$ is fitted to the boundary of $\Omega$, but we require that $T \cap \Omega \neq \emptyset \,, \, 
{\forall\,}T \in {\mathcal{T}_h}$. Typically, the mesh ${\mathcal{T}_h}$ can be thought of as a suitable sub-mesh of a larger and easy to generate mesh, see Figure \[fig:computational-domain\].  The domain ${\Omega^{\ast}}$ consisting of the union of all elements $T \in {\mathcal{T}_h}$ is called the *fictitious domain*.  For mesh faces in ${\mathcal{T}_h}$, i.e. edges of elements in two dimensions and faces in three dimensions, we distinguish between *exterior faces*, $\mathcal{F}_e$, which are faces that belong to one element only and are thus part of the boundary $\partial {\Omega^{\ast}}$ and *interior faces*, $\mathcal{F}_i$, which are faces that are shared by two elements in ${\mathcal{T}_h}$.

Next, let $G_h$ be the subset of elements in ${\mathcal{T}_h}$ that intersect the boundary ${\Gamma}$ $$\label{eq:define-cutting-cell-mesh}
G_h = \{T \in {\mathcal{T}_h}: T \cap {\Gamma}\neq \emptyset \}$$ and let us introduce the notation ${\mathcal{F}}_{\Gamma}$ for the set of all interior faces belonging to elements intersected by the boundary ${\Gamma}$ $$\label{eq:define-cutting-face-mesh}
  {\mathcal{F}}_{\Gamma} = \{ F \in {\mathcal{F}}_i :\;
  T^+_F \cap {\Gamma}\neq  \emptyset
  \vee
  T^-_F \cap {\Gamma}\neq  \emptyset
  \}.$$ Here, $T^+_F$ and $T^-_F$ are the two elements sharing the interior face $F \in \mathcal{F}_i$. Figure \[subfig: face notation\] illustrates these notations.  To ensure that the boundary $\Gamma$ is reasonably resolved by ${\mathcal{T}_h}$, we make the following assumptions:

-   G1: The intersection between $\Gamma$ and a face $F \in
      {\mathcal{F}}_i$ is simply connected; that is, $\Gamma$ does not cross an interior face multiple times. 

-   G2: For each element $T$ intersected by $\Gamma$, there exists a plane $S_T$ and a piecewise smooth parametrization $\Phi: S_T \cap T
      \rightarrow \Gamma \cap T$.

-   G3: We assume that there is an integer $N>0$ such that for each element $T \in G_h$ there exists an element $T' \in
      {\mathcal{T}_h}\setminus G_h$ and at most $N$ elements $\{T\}_{i=1}^N$ such that $T_1 = T,\,T_N = T'$ and $T_i \cap T_{i+1}
      \in {\mathcal{F}}_i,\; i = 1,\ldots N-1$. In other words, the number of faces to be crossed in order to “walk” from a cut element $T$ to a non-cut element $T' \subset \Omega$ is uniformly bounded. 

Similar assumptions were made in [@HansboHansbo2002; @BurmanHansbo2012; @MassingLarsonLoggEtAl2013a].  Next, we introduce the continuous linear finite element spaces $$\begin{aligned}
  \mathcal{V}_h &= \left\{ v_h \in C^0(\Omega): { \left.  v_h \right|_{T}} \in
    \mathcal{P}_1(T)\, {\forall\,}T \in {\mathcal{T}_h}\right\}\end{aligned}$$ for the pressure, $$\begin{aligned}
  \mathcal{V}_h^d &= \left\{{{v}}_h \in C^0(\Omega):
    { \left. {{ v}}_h \right|_{T}} \in [\mathcal{P}_1(T)]^d\, {\forall\,}T \in
    {\mathcal{T}_h}\right\}\end{aligned}$$ for the velocity and $$\begin{aligned}
  \mathcal{V}_h^{d \times d} &= \left\{ {\sigma}_h \in C^0(\Omega):
    { \left.  \mathbf{{\sigma}}_h \right|_{T}} \in [\mathcal{P}_1(T)]^{d \times d}\,
    {\forall\,}T \in {\mathcal{T}_h}\right\}\end{aligned}$$ for the extra-stress tensor.  We combine these spaces in the mixed finite element space $$\label{eq:total-approx-space}
  {\mathbb{V}}_h = \mathcal{V}_h^{d \times d}
  \times \mathcal{V}_h^d \times \mathcal{V}_h.$$ When we assume that the computational mesh matches the domain, the use of equal order approximation spaces for pressure, velocity and extra-stress tensor does not result in a stable discretization in the sense of Babuška–Brezzi.  However, a stable discretization can be obtained by augmenting the variational formulation  with suitable stabilization terms, see for instance [@BonvinPicassoStenberg2001], where Galerkin least-square techniques were employed, and [@Bonito2008], where interior penalty operators were used.  Remarkably, the same interior penalty techniques have been proved beneficial in [@BurmanHansbo2010; @BurmanHansbo2012; @BurmanHansbo2013; @MassingLarsonLoggEtAl2013a] to devise fictitious domain methods which are robust irrespective of how the boundary cuts the mesh.  In this work, we therefore employ interior penalty operators to cope with both inf-sup and cut geometry related instabilities in a unified way. 

For a quantity $x$ representing either a scalar, vector or tensor-valued function on ${\Omega^{\ast}}$, the interior penalty operators are defined by $$\label{eq:stabilization_operator}
  s_{k}(x,y) = \sum\limits_{F \in \mathcal{F}_i} h^k  
  \left( {\left\llbracket \nabla x\right\rrbracket}_n {\left\llbracket \nabla y\right\rrbracket}_n \right)_F.$$ Here, ${\left\llbracket \nabla x\right\rrbracket}_n$ denotes the normal jump of the quantity $x$ over the face, $F$, defined as ${\left\llbracket \nabla x\right\rrbracket}_n = \left.  \nabla x
\right|_{T_F^+} n_F  - \left.  \nabla x
\right|_{T_F^-} n_F$, where $n_F$ denotes a unit normal to the face $F$ with fixed but arbitrary orientation.  With this notation, the stabilization operators for pressure-velocity and velocity-extra-stress coupling employed in [@Bonito2008] are then given by $$\begin{aligned}
  \label{eq:stabilization-operator-p}
  s_p(p_h,q_h) &= \frac{\gamma_p}{2 \eta}
  s_{3}(p_h,q_h),
  \\
  \label{eq:stabilization-operator-u}
  s_u({{v}}_h,{{u}}_h) &= 2 \eta \gamma_u
  s_{1}({{v}}_h,{{u}}_h),\end{aligned}$$ with $\gamma_u$, $\gamma_p$ being positive penalty parameters to be determined later. 

We would like to emphasize that the face-wise contributions in the stabilization operator  are always computed on the *entire* face $F$, even for cut faces where $F
  \cap \overline{\Omega} \not = F$. Following the nomenclature in [@Burman2010; @BeckerBurmanHansbo2009], we refer to such penalties evaluated outside the physical domain as ghost-penalties.  Such ghost-penalties have been proven to be crucial in extending classical finite elements with weakly imposed Dirichlet boundary conditions to the fictitious domain case.  In that sense, the stabilization serves two purposes: to stabilize the discretisation scheme when using equal order interpolation spaces and to make the scheme insensitive to the boundary position with respect to the mesh. 

In addition, we define a stabilization operator for the extra-stress variable $$\label{eq:stabilization-operators}
  s_{{\sigma}}({\sigma}_h,{{{ \tau}}}_h) = \frac{\gamma_{\sigma}}{2\eta}
  \sum_{F \in {\mathcal{F}_{\Gamma}}} h^3
  ({\left\llbracket \nabla {\sigma}_h\right\rrbracket}_n, {\left\llbracket \nabla {{{ \tau}}}_h\right\rrbracket}_n)_F,$$ which acts *only* on ${\mathcal{F}_{\Gamma}}$ as opposed to the stabilization operators $s_u(\cdot,\cdot)$ and $s_p(\cdot,\cdot)$ which act on the set ${\mathcal{F}}_i$ of all interior faces.  The sole purpose of the additional stabilization $s_{\sigma}(\cdot,\cdot)$ is to ensure the stability of the discretization with respect to how the boundary cuts the mesh.  Here, $\gamma_{\sigma}$ denotes a positive penalty parameter. 

We are now in the position to formulate our stabilized cut finite element method for the three field Stokes problem.  Introducing the bilinear forms $$\begin{aligned}
  \label{eq:a_h-definition}
  a_h\left({\sigma}_h, {{v}}_h \right) &=  \left({\sigma}_h, \,
{{\epsilon}}({{v}}_h) \right)_{\Omega} - 
{\left(  
{\sigma}_h \cdot {{n}}, \, {{v}}_h  \right)}_{\partial \Omega}, 
\\
\label{eq:b_h-definition}
b_h\left(p_h, {{v}}_h \right) &= - \left(p_h, \, \nabla \cdot
    {{v}}_h \right)_{\Omega} + 
    {\left( p_h {{n}}, \, {{v}}_h  \right)}_{\partial \Omega},
  \end{aligned}$$ the proposed discretization scheme reads: Find $U_h:=\left({\sigma}_h,{{u}}_h,p_h\right) \in {\mathbb{V}}_h$ such that for all $V_h:=({{{ \tau}}}_h,{{v}}_h,q_h) \in {\mathbb{V}}_h$ $$\begin{aligned}
    A_h(U_h,V_h) + S_h(U_h,V_h) = L(V_h), 
\label{prob:tfstokes-cutfem}\end{aligned}$$ where $$\begin{aligned}
  \label{eq:A_h-definition}
    A_h(U_h,V_h) =& 
    \frac{1}{2\eta}  \left({\sigma}_h, \, {{{ \tau}}}_h \right)_{\Omega} 
    + \frac{\gamma_b \eta}{h} 
    {\left( 
    u_h,\,v_h \right)}_{\partial \Omega} 
    \nonumber
    + a_h\left({\sigma}_h,  {{v}}_h \right) 
    - a_h\left({{{ \tau}}}_h, {{u}}_h \right) 
    \\
    &\phantom{=}
    + b_h\left(p_h, {{v}}_h \right) 
    - b_h\left(q_h, {{u}}_h \right), 
    \\
  \label{eq:L_h-definition}
  L_h(V_h) =&  ({{f}}, {{v}}_h)_{\Omega} + 
  {\left(  {{{ \tau}}}_h \cdot {{n}}, \, {{g}}\right)}_{\partial \Omega} 
  - {\left(  q_h {{n}}, \, {{g}}\right)}_{\partial \Omega} + \frac{\gamma_b \eta}{h} 
  {\left( 
    {{g}},\,{{v}}_h
   \right)}_{\partial \Omega},
  \end{aligned}$$ with the stabilization operator $$\begin{aligned}
  \label{eq:S_h-definition}
    S_h(U_h,V_h) =& 
    s_{{\sigma}}({\sigma}_h, {{{ \tau}}}_h)
    + s_u({{u}}_h,{{v}}_h) 
    + s_p(p_h,q_h). 
  \end{aligned}$$ Here, the positive penalty constant $\gamma_b$ arises from the weak enforcement of Dirichlet boundary conditions through Nitsche’s method and $h=\max_{T\in {\mathcal{T}_h}} h_T$ is the mesh size, where $h_T$ denotes the diameter of $T$.

For the sake of keeping the technical details presented in this work at a moderate level, we assume for the following a priori error analysis that the contributions from the cut elements $\{ T \cap \Omega : T \in G_h\}$ and boundary parts $\{ T \cap
  \Gamma : T \in G_h\}$ can be computed exactly.  For a thorough treatment of variational crimes arising from the discretization of a curved boundary in the context of cut finite element methods, we refer to [@BurmanHansboLarsonEtAl2014]. 

Approximation properties {#sec:approximation-properties}
========================

In this section, we summarize some useful inequalities and interpolation estimates which are necessary to proceed with the stability and a priori analysis in Section \[sec:stability\] and \[sec:apriori-analysis\]. 

Norms
-----

Recalling the notation from Section \[sec:fem-formulation\], we introduce the triple norm $$\begin{aligned}
  {|\mspace{-1mu}|\mspace{-1mu}|}U {|\mspace{-1mu}|\mspace{-1mu}|}^2 
  &= \dfrac{1}{2\eta} \| {\sigma}\|^2_{\Omega}
  + 2\eta \| {{\epsilon}}({{u}}) \|^2_{\Omega}
  + \dfrac{1}{2\eta} \| p \|^2_{\Omega}
  + 2 \eta \gamma_b \|h^{-1/2} {{u}}\|^2_{\Gamma}
  \label{eq:triple-norm}\end{aligned}$$ for the three field Stokes problem  and its discrete counterpart $$\begin{aligned}
  {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h^2 
  &= {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}^2  
  + S_h(U_h, U_h)
  \label{eq:triple-norm-disc}\end{aligned}$$ which will be instrumental in studying the stability and convergence properties of problem  in the following sections. 

A key point in the definition of the discrete norm ${|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h$ is that the inclusion of the stabilization terms allows to reconstruct natural norms for the discrete function $U_h$ defined on the *entire* fictitious domain ${\Omega^{\ast}}$. More specifically, @BurmanHansbo2012 [@MassingLarsonLoggEtAl2013a] proved the following lemma:

\[lem:control-norm-with-ghost-penalties\] Let $\Omega$, ${\Omega^{\ast}}$ and ${\mathcal{F}_{\Gamma}}$ be defined as in Section \[sec:fem-formulation\].  Then for all $v_h \in {\mathcal{V}}_h$ it holds $$\| v_h \|_{{\Omega^{\ast}}}^2 
  \lesssim
  \| v_h \|_{\Omega}^2
  + \sum_{F \in  {\mathcal{F}_{\Gamma}}} h_F^3
  ({\left\llbracket \nabla v_h\right\rrbracket}_n, {\left\llbracket \nabla v_h\right\rrbracket}_n)_F
  \lesssim
  \| v_h \|_{{\Omega^{\ast}}}^2,
  \label{eq:u-fd-norm-est1}$$ $$\| \nabla v_h \|_{{\Omega^{\ast}}}^2
  \lesssim
  \| \nabla v_h
  \|_{\Omega}^2
  + \sum_{F \in {\mathcal{F}_{\Gamma}}} h_F
  ({\left\llbracket \nabla v_h\right\rrbracket}_{n}, {\left\llbracket \nabla v_h\right\rrbracket}_n)_F
  \lesssim
  \| \nabla v_h \|_{{\Omega^{\ast}}}^2. 
  \label{eq:u-fd-norm-est}$$

Here and throughout, we use the notation $ a \lesssim b$ for $ a
\leqslant C b$ for some generic constant $C$ which varies with the context but is always independent of the mesh size $h$.

Useful inequalities
-------------------

We recall the following trace inequality (for a proof, see e.g. [@DipietroErn2011]) for $v \in H^1({\Omega^{\ast}})$ $${\left\| v\right\|}_{\partial T} 
    \lesssim
    h_T^{-1/2} {\left\| v\right\|}_{T} +
    h_T^{1/2} {\left\| \nabla v\right\|}_{T}
    \quad {\forall\,}T \in {\mathcal{T}_h}. 
    \label{eq:trace-inequality}$$ If the intersection $\Gamma \cap T$ does not coincide with a boundary edge of the mesh and if the intersection is subject to conditions G1)–G3), then the corresponding inequality $$\begin{aligned}
    {\left\| v\right\|}_{T \cap \partial \Omega} 
    \lesssim
      h_T^{-1/2} {\left\| v\right\|}_{T} 
      + h_T^{1/2} {\left\| \nabla v\right\|}_{T} 
    \label{eq:trace-inequality-unfitted}
  \end{aligned}$$ holds (see @HansboHansbo2002).  We will also use the following well-known inverse estimates for $v_h \in {\mathcal{V}}_h$ $$\begin{aligned}
{3}
    \label{eq:inverse-estimate-I}
    {\left\| \nabla v_h\right\|}_{T} 
    & \lesssim
    h_T^{-1} 
    {\left\| v_h\right\|}_{T} & & \quad {\forall\,}T \in {\mathcal{T}_h},
    \\ 
    \label{eq:inverse-estimate-II}
    {\left\| \nabla v_h \cdot n\right\|}_{\partial T} 
    & \lesssim
    h_T^{-1/2} 
    {\left\| \nabla v_h\right\|}_{T} & & \quad {\forall\,}T \in {\mathcal{T}_h}. 
  \end{aligned}$$ For elements $T$ intersected by the boundary $\Gamma$, we have $$\begin{aligned}
{3}
    {\left\| \nabla v_h \cdot n\right\|}_{T \cap \Gamma} 
    & \lesssim
    h_T^{-1/2} {\left\| \nabla v_h\right\|}_{T} & & \quad {\forall\,}T \in
    {\mathcal{T}_h},
    \label{eq:inverse-estimate-III}
  \end{aligned}$$ which is proven in [@HansboHansbo2002].  Finally, we recall the well-known Korn inequalities, stating that $$\begin{gathered}
    \| \nabla {{v}}\|_{\Omega} 
    \lesssim 
    \| {{\epsilon}}({{v}}) \|_{\Omega} \quad {\forall\,}{{v}}\in
    [H_0^1(\Omega)]^d,
    \label{eq:korn-I}
    \\
    \| {{v}}\|_{1,\Omega} 
    \lesssim 
    \| {{\epsilon}}({{v}}) \|_{\Omega} 
    + \| {{v}}\|_{\Omega}
    \quad {\forall\,}{{v}}\in
    [H^1(\Omega)]^d. 
    \label{eq:korn-II}
  \end{gathered}$$ Since the norm  incorporates the boundary data of $u$, it will be most convenient to work with the following variant of Korn’s inequality :

For ${{v}}\in [H^1(\Omega)]^d$, it holds that $$\begin{aligned}
      \| {{v}}\|_{1,\Omega} 
      \lesssim 
      \| {{\epsilon}}({{v}}) \|_{\Omega} 
      + \| {{v}}\|_{\partial \Omega}. 
      \label{eq:korn-III}
    \end{aligned}$$ \[lem:korn-III\]

For the sake of completeness, we provide a short proof here, which is established by contradiction.  Assuming that  does not hold, we can construct a sequence $\{{{v}}_n\}_n$ such that $\| {{v}}_n \|_{1,\Omega} = 1$ and $$\begin{aligned}
    \| {{\epsilon}}({{v}}_n) \|_{\Omega} 
    + \| {{v}}_n \|_{\partial \Omega}
    \leqslant \dfrac{1}{n}. 
    \label{eq:korn-contradiction}
    \end{aligned}$$ The compact embedding $H^1(\Omega) \subset\subset L^2(\Omega)$ implies that there is a subsequence $\{{{v}}_{n'}\}_{n'}$ which converges in the $\| \cdot \|_{0, \Omega}$-norm.  Since by construction, $\| {{\epsilon}}({{v}}_n - {{v}}_m) \|_{\Omega} \leqslant \tfrac{1}{n} +
    \tfrac{1}{m}$, we conclude using Korn’s inequality  that $\{{{v}}_{n'}\}_{n'}$ is also a Cauchy sequence in $[H^1(\Omega)]^d$ with ${{v}}_{n'}
    \overset{\|\cdot\|_{1,\Omega}}{\to} {{v}}'$ for some ${{v}}' \in [H^1(\Omega)]^d$. Due to the boundedness of the trace operator $T:[H^1(\Omega)]^d \to [L^2(\Omega)]^d$ and  , we observe that ${{v}}' \in [H^1_0(\Omega)]^d$. Now Poincaré’s inequality together with Korn’s inequality – gives the contradiction $$\begin{aligned}
      1 = \| {{v}}' \|_{1,\Omega} \lesssim \|{{\epsilon}}({{v}}')\|_{\Omega}
      + \| {{v}}'
      \|_{\Omega} \lesssim \| \nabla {{v}}' \|_{\Omega} 
      \lesssim  \|{{\epsilon}}({{v}}')\|_{\Omega} = 0. 
    \end{aligned}$$

Interpolation and projection operators
--------------------------------------

Before we construct various interpolation operators $L^2(\Omega) \to
{\mathcal{V}}_h$, we recall that for a Lipschitz-domain $\Omega$, an extension operator $$E: H^s({\Omega}) \rightarrow H^s({\Omega^{\ast}})$$ can be defined which is bounded $$\label{eq:extension-operator-boundedness}
  {\left\| E v\right\|}_{s,{\Omega^{\ast}}} \lesssim {\left\| v\right\|}_{s,\Omega}, \quad s=0,1,2,$$ see [@Stein1970] for a proof.  Occasionally, we write $v^* = Ev$. Then, for any interpolation operator $I_h: H^s({\Omega^{\ast}}) \to {\mathcal{V}}_h$, we can define its “fictitious domain” variant $I_h^*: H^s(\Omega) \to {\mathcal{V}}_h$ by simply requiring that $$I_h^*u = I_h(u^*)$$ for $u \in H^s(\Omega)$. In particular, we will choose $I_h$ to be the Clément and Oswald interpolation operators, which we denote by ${\mathcal{C}}_h$ and ${\mathcal{O}}_h$, respectively (see for instance [@ErnGuermond2004]).  Recalling that the standard interpolation estimates for the Clément interpolant $$\begin{aligned}
{3}
\| v - {\mathcal{C}}_h v \|_{r,T} 
& \lesssim
  h^{s-r}| v |_{s,\omega(T)},
&\quad 0\leqslant r \leqslant s \leqslant 2 \quad &{\forall\,}T\in {\mathcal{T}_h},
  \label{eq:interpest0}
  \\
\| v - {\mathcal{C}}_h v \|_{r,F} &\lesssim h^{s-r-1/2}| v |_{s,\omega(T)},
&\quad 0\leqslant r \leqslant s \leqslant 2 \quad &{\forall\,}F\in
  {\mathcal{F}}_i 
  \label{eq:interpest1}\end{aligned}$$ hold if $v \in H^s({\Omega^{\ast}})$, we observe that the extended Clément interpolant ${\mathcal{C}}_h^*$ satisfies $$\begin{aligned}
{3}
  \| v^* - {\mathcal{C}}_h^* v \|_{r,{\Omega^{\ast}}} 
  & \lesssim
  h^{s-r}\| v \|_{s,\Omega},
  & &\quad 0\leqslant r \leqslant s \leqslant 2,
  \label{eq:interpest0-ast}
  \\
  \sum_{F\in {\mathcal{F}}} 
  \| v^* - {\mathcal{C}}_h^* v \|_{r,F} 
  &\lesssim h^{s-r-1/2} \| v \|_{s,\Omega},
  & &\quad 0\leqslant r \leqslant s \leqslant 2
  \label{eq:interpest1-ast}\end{aligned}$$ due to the boundedness of the extension operator .  Here, $\omega(T)$ is the set of elements in ${\mathcal{T}_h}$ sharing at least one vertex with $T$ (for ) and sharing at least one vertex with $F \in \mathcal{F}_i$ (for ), respectively.  The Oswald interpolation operator ${\mathcal{O}}_h: H^2({\mathcal{T}_h}) \to {\mathcal{V}}_h$ is of particular use in the context of continuous interior penalty methods, as it allows to control the fluctuation $\nabla v_h -
{\mathcal{O}}_h(\nabla v_h)$ in terms of the stabilization operator .  More precisely, @BurmanFernandezHansbo2006 proved the following lemma:

\[lem:cip-fluctuation-control\] For all $v_h \in {\mathcal{V}}_h$ $$\label{eq:cip-fluctuation-control}
    {\left\| h \left( \nabla v_h - {\mathcal{O}}_h\left(\nabla v_h
    \right)\right)\right\|}_{{\Omega^{\ast}}}^2 
    \lesssim
    s_{3}(v_h,v_h).$$

To exploit this control given by the stabilization operators $s_{3}(\cdot, \cdot)$ in the stability analysis of our fictitious domain method, we define the *stabilized* approximate $L^2$-projection $$\Pi^*_h : L^2(\Omega) \rightarrow {\mathcal{V}}_h$$ by $$\left(\Pi_h^* u, v_h \right)_{\Omega}
  + s_{3}\left(\Pi_h^* u, v_h \right)
  =
  \left( u, v_h \right)_{\Omega}. 
  \label{eq:stab-projection-def}$$ We conclude this section by proving certain approximation properties of the stabilized $L^2$-projection.  We start with the following. 

\[lem:l2-stability\] For $u\in L^2(\Omega)$ it holds $$\begin{aligned}
    {\left\| \Pi_h^* u\right\|}_{\Omega} \leq {\left\| 
    u\right\|}_{\Omega}, \\
    {\left\| \Pi_h^* u\right\|}_{\Omega^*} \lesssim {\left\| 
    u\right\|}_{\Omega}. 
    \label{eq:l2-stability}
  \end{aligned}$$

Using the property of the stabilized $L^2$-projection, a Cauchy-Schwarz inequality and Lemma \[lem:control-norm-with-ghost-penalties\], we obtain $$\begin{aligned}
    {\left\| \Pi_h^* u\right\|}_{\Omega}^2 
    &\leq  {\left\|  \Pi_h^* u \right\|}_{\Omega}^2  + s_3(\Pi_h^* u, \Pi_h^* u) \nonumber \\
    &= \left(u, \Pi_h^* u \right)_{\Omega}  \nonumber \\
    &\leq  {\left\| u\right\|}_{\Omega}   {\left\| \Pi_h^* u\right\|}_{\Omega}.  
  \end{aligned}$$ $$\begin{aligned}
{\left\| \Pi_h^*
  u \right\|}_{{\Omega^{\ast}}}^2 
    &\lesssim  {\left\|  \Pi_h^* u \right\|}_{\Omega}^2  + s_3(\Pi_h^* u, \Pi_h^* u). 
  \end{aligned}$$

\[prop:interpolation-estimate-stabilized-projection\] Assuming a quasi-uniform triangulation, the stabilized projection operator $\Pi_h^*$ satisfies the following approximation property for $u \in H^s(\Omega)$ $$h^{1/2} {\left\| u - \Pi_h^* u\right\|}_{\partial \Omega} + {\left\| u - \Pi_h^*
      u\right\|}_{\Omega} + h {\left\| \nabla \left(u - \Pi_h^* u \right)\right\|}_{\Omega} \leq C
    h^s {\left| u\right|}_{s}. 
    \label{eq:interpolation-estimate-stabilized-projection}$$

We begin by writing equation  as $I + II +
  III$. We first proof the $L^2$-error estimate $$II = {\left\| \Pi_h^* u - u\right\|}_{\Omega}  \leq C h^s
    {\left| u\right|}_{H^s(\Omega)}.$$ Using Lemma \[lem:control-norm-with-ghost-penalties\] and the triangle inequality, we obtain $$\begin{aligned}
     {\left\| \Pi_h^*u - u\right\|}_{\Omega}  \lesssim  {\left\| \Pi_h^*u - u^*\right\|}_{{\Omega^{\ast}}} &\leq
    \| \underbrace{\Pi_h^*u - \mathcal{C}_h^*u}_{\xi_h} \|_{{\Omega^{\ast}}}
    + {\left\| \mathcal{C}_h^*u - u^*\right\|}_{{\Omega^{\ast}}}. 
  \end{aligned}$$ Next, using the approximation properties of the Clément operator ${\mathcal{C}}_h^{\ast}$ , it is enough to estimate $\xi_h$ as follows $$\begin{aligned}
 {\left\| \xi_h\right\|}_{{\Omega^{\ast}}}^2 &\lesssim \left( \xi_h,\, \xi_h\right)_{\Omega} + s_3 \left( \xi_h ,\, \xi_h\right) \nonumber \\
    &= \left( u - \mathcal{C}_h^* u ,\, \xi_h\right)_{\Omega} +
    s_3 \left( \mathcal{C}_h^* u, \, \xi_h\right) \nonumber \\
    &\lesssim {\left\| u - \mathcal{C}_h^* u\right\|}_{{\Omega^{\ast}}}
    {\left\| \xi_h\right\|}_{\Omega}
    + s_3 \left( \mathcal{C}_h^* u, \,  \mathcal{C}_h^* u \right)^{1/2}  s_3 \left( \xi_h, \,  \xi_h \right)^{1/2}  \nonumber \\
    &\lesssim \left( {\left\| u - \mathcal{C}_h^* u\right\|}_{{\Omega^{\ast}}}
      +  s_3 \left( \mathcal{C}_h^* u, \,  \mathcal{C}_h^* u
    \right)^{1/2} \right) {\left\| \xi_h\right\|}_{{\Omega^{\ast}}}, \nonumber 
  \end{aligned}$$ where in the first and fourth line, we used Lemma \[lem:control-norm-with-ghost-penalties\] and in the second line, we used the property   to pass from $\Pi_h^* u$ to $u$. Consequently, $$\begin{aligned}
   {\left\| \xi_h\right\|}_{{\Omega^{\ast}}}  &\lesssim h^s {\left| u^*\right|}_{H^s({\Omega^{\ast}})} 
      \nonumber 
      \\
    &\lesssim h^s {\left| u\right|}_{H^s(\Omega)}.\end{aligned}$$ Next, we will prove $$III =     h {\left\| \nabla(u - \Pi_h^* u)\right\|}_{\Omega} \lesssim   h {\left\| \nabla(u^* - \Pi_h^* u)\right\|}_{\Omega^*}  \leq C h^s
    {\left| u\right|}_{H^s(\Omega)}.$$ Using Lemma \[lem:control-norm-with-ghost-penalties\], the approximation properties of the Clément interpolant, the triangle inequality, an inverse inequality in combination with the $L^2$-approximation property of the stabilized projection, the third term can be estimated as follows: $$\begin{aligned}
  h {\left\| \nabla(u^* - \Pi_h^* u)\right\|}_{\Omega^*} &\lesssim  h( {\left\| \nabla(u^* -
  {\mathcal{C}}_h^{\ast} u)\right\|}_{\Omega^*} + {\left\| \nabla({\mathcal{C}}_h^{\ast} u - \Pi_h^*
      u)\right\|}_{\Omega^*})
    \\
    & \lesssim  h ({\left\| \nabla(u^* - {\mathcal{C}}_h^{\ast} u)\right\|}_{{\Omega^{\ast}}} +
    h^{-1}{\left\| {\mathcal{C}}_h^{\ast} u^* - \Pi_h^* u\right\|}_{{\Omega^{\ast}}})
    \\
    & \lesssim  h ({\left\| \nabla(u^* - {\mathcal{C}}_h^{\ast} u)\right\|}_{{\Omega^{\ast}}} +
    h^{-1}({\left\| {\mathcal{C}}_h^{\ast} u - u\right\|}_{{\Omega^{\ast}}} + {\left\| u^* - \Pi_h^*
      u\right\|}_{{\Omega^{\ast}}} ))
    \\
    & \lesssim  h^s {\left| u\right|}_{s,\Omega}. 
  \end{aligned}$$ We conclude the proof by bounding the first term via the trace inequality  $$\begin{aligned}
    h^{1/2} {\left\| u - \Pi_h^* u\right\|}_{\Gamma} \lesssim {\left\|  u^* - \Pi_h^*
      u\right\|}_{{\Omega^{\ast}}} + h {\left\| \nabla(u^* - \Pi_h^* u)\right\|}_{{\Omega^{\ast}}}
  \end{aligned}$$ and using the estimates for $II$ and $III$.

\[lem:h1-stability\] For $u\in H^1(\Omega)$ it holds $${\left\| \nabla \Pi_h^* u\right\|}_{\Omega} \lesssim {\left\| \nabla
    u\right\|}_{\Omega}. 
    \label{eq:h1-stability}$$

The desired estimate follows from Proposition \[prop:interpolation-estimate-stabilized-projection\]: $$\begin{aligned}
    {\left\| \nabla \Pi_h^* u\right\|}_{\Omega} &\leq  {\left\| \nabla (\Pi_h^* u - u)\right\|}_{\Omega} + {\left\| \nabla
    u\right\|}_{\Omega} \nonumber \\
&\leq C {\left\| \nabla
    u\right\|}_{\Omega}  + {\left\| \nabla
    u\right\|}_{\Omega}.  \end{aligned}$$

Stability estimates {#sec:stability}
===================

In this section, we prove that the stabilized cut finite element formulation  for the three field Stokes problem  fulfills an inf-sup condition in the Babuška–Brezzi sense.  As a first step, we demonstrate that the pressure stabilization given in  allows to formulate a weakened inf-sup condition for the pressure-velocity coupling when equal-order interpolation spaces are employed.  Similar estimates have previously been stated in [@BurmanHansbo2006a] for the classical matching mesh case and in [@BurmanHansbo2013; @MassingLarsonLoggEtAl2013a] for a Nitsche-based fictitious domain formulation for the Stokes problem.  Introducing the discrete velocity norm $$\begin{aligned}
  \| v_h \|_{1, h}^2 = \| v_h \|_{1, {\Omega^{\ast}}}^2 + \gamma_b \|h^{-1/2}
  v_h \|_{\Gamma}^2
  \label{eq:discrete-vel-norm}\end{aligned}$$ for $v_h \in [{\mathcal{V}}_h]^d$, we can state the following

\[prop:mod-inf-sup-bh\] Let $p_h \in {\mathcal{V}}_h$, then there is a constant $c > 0$ such that $$\begin{aligned}
    \sup_{{{v}}_h \in {\mathcal{V}}_h^d\setminus \{0\}} \dfrac{b_h(p_h,{{v}}_h)}{\| {{v}}_h
    \|_{1,h}}
    \gtrsim
    \| p_h \|_{\Omega} -  c s_{3}(p_h,p_h)^{1/2}.  
    \label{eq:mod-inf-sup-bh}
  \end{aligned}$$

Due to the surjectivity of the divergence operator $\nabla \cdot
  : [H_0^1(\Omega)]^d \to L^2(\Omega)$, there exists a $v^p \in
  [H_0^1(\Omega)]^d$ such that $ \nabla \cdot v^p = p_h$ and $ \| v^p
  \|_{1,\Omega} \lesssim \| p_h \|_{\Omega}$. Setting $v_h =
  \Pi^{\ast}_h v^p$ and using the $H^1$-stability of the stabilized $L^2$-projection stated in Lemma \[lem:h1-stability\], we thus obtain $$\begin{aligned}
    \nonumber
    b_h(p_h, v_h) 
    &= b_h(p_h, v^p) + b_h(p_h, \Pi_h^{\ast} v^p - v^p) 
    \\
    &\gtrsim \| p_h \|_{\Omega} \| v^p \|_{1,\Omega} + 
    b_h(p_h, \Pi_h^{\ast} v^p - v^p). 
    \label{eq:b_h-split}
  \end{aligned}$$ Next, we estimate the remaining term in .  Recalling definition  of $b_h(\cdot, \cdot)$ and integrating by parts gives $$\begin{aligned}
    b_h(p_h, \Pi_h^{\ast} v^p - v^p)
    &= \left( \nabla p_h, \Pi^{\ast}_h v^p - v^p \right)_{\Omega}. 
    \label{eq:b_h-ibp}
  \end{aligned}$$ We now exploit the (almost) orthogonality of the stabilized $L^2$-projection $\Pi_h^{\ast} : [H^1_0(\Omega)]^d \to {\mathcal{V}}^{d}$ by inserting ${\mathcal{O}}_h(\nabla p_h) \in {\mathcal{V}}^{d}$ into , yielding $$\begin{aligned}
    \nonumber
    b_h(p_h, \Pi_h^{\ast} v^p - v^p)
    &= \left( \nabla p_h - {\mathcal{O}}_h(\nabla p_h), \Pi^{\ast}_h v^p - v^p \right)_{\Omega}
    - s_{3}({\mathcal{O}}_h(\nabla p_h),  \Pi^{\ast}_h v^p) = I + II. 
  \end{aligned}$$ Combining Lemma \[lem:cip-fluctuation-control\] with the stability and approximation properties of $\Pi_h^{\ast}$, cf.   and , the first term can be bounded as follows: $$\begin{aligned}
    I
    &\gtrsim 
    -s_{3}(p_h, p_h)^{1/2} \|h^{-1} (\Pi^{\ast}_h v^p - v^p)
    \|_{\Omega}
    \gtrsim 
    - s_{3}(p_h, p_h)^{1/2} \| v^p \|_{1, \Omega}
    \label{eq:b_h-I-estimate}
  \end{aligned}$$ To estimate $II$, recall the definition of ${{v}}_h$ and apply successively a Cauchy-Schwarz inequality and Lemma \[lem:control-norm-with-ghost-penalties\] to obtain $$\begin{aligned}
    \nonumber
   II
   \gtrsim - s_{1}(v_h,v_h)^{1/2} s_{5}\left({\mathcal{O}}_h(\nabla
   p_h),{\mathcal{O}}_h(\nabla p_h)\right)^{1/2}
   &\gtrsim 
   - \| v_h \|_{1,{\Omega^{\ast}}} 
   s_{5}\left({\mathcal{O}}_h(\nabla p_h),{\mathcal{O}}_h(\nabla p_h)\right)^{1/2}. 
  \end{aligned}$$ Using successively the discrete trace inequality , the inverse inequality  and Lemma \[lem:cip-fluctuation-control\], the last term can be bounded in the following way: $$\begin{aligned}
    \nonumber
   s_{5}\left({\mathcal{O}}_h(\nabla p_h),{\mathcal{O}}_h(\nabla p_h)\right)
   &= 
   s_{5}\left({\mathcal{O}}_h(\nabla p_h) - \nabla p_h,
   {\mathcal{O}}_h(\nabla p_h) - \nabla p_h \right)
   \\
    \nonumber
   & \lesssim
   \sum_{T \in {\mathcal{T}_h}} h^4 \| \nabla ({\mathcal{O}}_h(\nabla p_h) - \nabla p_h )
   \|_T
   \\
   \label{eq:b_h-inconsistency-estimate}
   & \lesssim
   \sum_{T \in {\mathcal{T}_h}} h^2 \| {\mathcal{O}}_h(\nabla p_h) - \nabla p_h \|_T
   \lesssim
   s_{3}\left(p_h, p_h \right). 
  \end{aligned}$$ Consequently, $$\begin{aligned}
  II \gtrsim - \| v_h \|_{1,{\Omega^{\ast}}} s_{3}(p_h,p_h)^{1/2}. 
   \label{eq:b_h-II-estimate}
  \end{aligned}$$ Combining  , and , we find that for some constants $c_1$ and $c_2$ $$\begin{aligned}
    b_h(p_h, v_h) 
    \gtrsim \left(
      \|p_h\|_{\Omega} 
      - c_1 s_{3}(p_h,p_h)^{1/2}
    \right) \| v^p \|_{1,\Omega} 
    - c_2 s_{3}(p_h,p_h)^{1/2} \| v_h \|_{1,{\Omega^{\ast}}}.  
    \label{eq:b_h-estimate-final}
  \end{aligned}$$ To conclude the proof, we note that since $v_p \in [H^1_0(\Omega)]^d$, we have $$\begin{aligned}
    \|v_h \|_{1, h}^2 
    &= \| v_h \|_{1, {\Omega^{\ast}}}^2 + h^{-1}\| v_h - v_p \|_{\Gamma}^2
    \lesssim
    \| v_p \|_{1, \Omega}^2,
  \end{aligned}$$ thanks to the stability  of the operator $\Pi^{\ast}_h$, the interpolation estimate  applied for $s = 1$ and our choice of $v_h$. As a result, $$\begin{aligned}
    \dfrac{b_h(p_h, v_h)}
    {\| v_h \|_{1,h}}
    \gtrsim 
    \dfrac{b_h(p_h, v_h)}
    {\| v^p \|_{1,\Omega}}
    \gtrsim 
      \|p_h\|_{\Omega} 
      - c s_{3}(p_h,p_h)^{1/2},
      \label{eq:b_h-estimate-final-II}
  \end{aligned}$$ if $b_h(p_h, v_h) \geqslant 0$, otherwise we can simply use $\widetilde{v}_h = - v_h$ to arrive at  with $v_h$ replaced by $\widetilde{v}_h$.

As a second step, we state and prove a weakened inf-sup condition for the coupling between the velocity and the extra-stress.  Here, the “defect” of the inf-sup condition is quantified in terms of the velocity stabilization form in  and the boundary penalization in . 

\[prop:mod-inf-sup-ah\] Let ${{u}}_h \in {\mathcal{V}}_h^d$, then there is a constant $c > 0$ such that $$\begin{aligned}
    \sup_{{{{ \tau}}}_h \in {\mathcal{V}}_h^{d\times d}\setminus \{0\}}
    \dfrac{a_h({{{ \tau}}}_h,\,{{u}}_h)}{\| {{{ \tau}}}_h \|_{{\Omega^{\ast}}}}
     \gtrsim     \| {{u}}_h \|_{1,\Omega} -  c 
    \left( s_{1}({{u}}_h,{{u}}_h)^{1/2} +
    {\left( h^{-1} {{u}}_h, {{u}}_h \right)}^{1/2}_{\Gamma}
    \right). 
    \label{eq:mod-inf-sup-ah}
  \end{aligned}$$

Choose ${{{ \tau}}}_h = \Pi_h^{\ast} {{\epsilon}}({{u}}_h)$. Then, by adding and subtracting ${{\epsilon}}({{u}}_h)$ and then ${\left( h^{-1} {{u}}_h, {{u}}_h \right)}_{\Gamma}^{1/2}$, we obtain $$\begin{aligned}
      \nonumber
        a_h( {{{ \tau}}}_h ,\, {{u}}_h) 
        &= (\Pi^{\ast}_h {{\epsilon}}({{u}}_h),\, {{\epsilon}}({{u}}_h))_{\Omega}
        - {\left( \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}},
        {{u}}_h \right)}_{\Gamma}
        \\
      \nonumber
        &= \| {{\epsilon}}({{u}}_h) \|_{\Omega}^2
            +\left(\Pi^{\ast}_h {{\epsilon}}({{u}}_h) -
            {{\epsilon}}({{u}}_h),\,{{\epsilon}}({{u}}_h)\right)_{\Omega}
          - {\left( \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}},
          {{u}}_h \right)}_{\Gamma}
        \\
      \nonumber
        &\gtrsim 
        \bigl(
          \| {{\epsilon}}({{u}}_h) \|_{\Omega}
          + {\left( h^{-1} {{u}}_h, {{u}}_h \right)}_{\Gamma}^{1/2}
        \bigr)
        \| {{\epsilon}}({{u}}_h) \|_{\Omega}
        - {\left( h^{-1} {{u}}_h, {{u}}_h \right)}_{\Gamma}^{1/2} 
        \| {{\epsilon}}({{u}}_h) \|_{\Omega}
        \\
      \nonumber
        & \phantom{\gtrsim}
        +
            \left(\Pi^{\ast}_h {{\epsilon}}({{u}}_h) -
            {{\epsilon}}({{u}}_h),\,{{\epsilon}}({{u}}_h)\right)_{\Omega}
            - 
          {\left( \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}},
          {{u}}_h \right)}_{\Gamma}
    \\
      \nonumber
        &\gtrsim 
    \| {{u}}_h \|_{1,\Omega}
        \| {{\epsilon}}({{u}}_h) \|_{\Omega}
        - {\left( h^{-1} {{u}}_h, {{u}}_h \right)}_{\Gamma}^{1/2} 
        \| {{\epsilon}}({{u}}_h) \|_{\Omega}
        \\
      \nonumber
        & \phantom{\gtrsim}
        + \underbrace{
            \left(\Pi^{\ast}_h {{\epsilon}}({{u}}_h) -
            {{\epsilon}}({{u}}_h),{{\epsilon}}({{u}}_h)\right)_{\Omega}}_{\text{I}}
            - \underbrace{
          {\left( \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}},
          {{u}}_h \right)}_{\Gamma}
    }_{\text{II}},
    \end{aligned}$$ where we used the $L^2$-stability of the stabilized $L^2$-projection and variant  of Korn’s inequality in the last two steps.  We proceed by estimating the terms $I$ and $II$ separately.\
[Estimate (I)]{}: The use of the stabilized $L^2$-projection $\Pi^{\ast}_h$ allows to insert the Oswald interpolant of ${{\epsilon}}({{u}}_h)$ yielding $$\begin{aligned}
      \nonumber
        I &= 
      ( \Pi^{\ast}_h
      {{\epsilon}}({{u}}_h) - {{\epsilon}}({{u}}_h),\, {{\epsilon}}({{u}}_h) - {\mathcal{O}}_h{{\epsilon}}({{u}}_h))_{\Omega} 
      - s_{3}\left({\mathcal{O}}_h {{\epsilon}}({{u}}_h), \Pi^{\ast}_h {{\epsilon}}({{u}}_h)\right)
      \nonumber
        \\ 
      \nonumber
        &\gtrsim -s_{1}({{u}}_h,{{u}}_h)^{1/2} \| \Pi^{\ast}_h {{\epsilon}}( {{u}}_h) -{{\epsilon}}(
        {{u}}_h) \|_{\Omega}
      - s_{3}\left({\mathcal{O}}_h {{\epsilon}}({{u}}_h), \Pi^{\ast}_h {{\epsilon}}({{u}}_h)\right)
      \nonumber
        \\
      \nonumber
        &\gtrsim 
      - s_{1}({{u}}_h,{{u}}_h)^{1/2} \|{{\epsilon}}( {{u}}_h) \|_{\Omega}
      - s_{3}\left({\mathcal{O}}_h {{\epsilon}}({{u}}_h), {\mathcal{O}}_h {{\epsilon}}({{u}}_h)\right)^{1/2}
        \| \Pi^{\ast}_h {{\epsilon}}( {{u}}_h) \|_{{\Omega^{\ast}}},
    \end{aligned}$$ where we successively applied Lemma \[lem:control-norm-with-ghost-penalties\], the $L^2$-boundedness of $\Pi^{\ast}_{h}$, and finally a Cauchy-Schwarz inequality.  By an argument similar to  in the previous Proposition \[prop:mod-inf-sup-bh\], we can show that $$\begin{aligned}
      \nonumber
        s_{3}\left({\mathcal{O}}_h {{\epsilon}}({{u}}_h), {\mathcal{O}}_h {{\epsilon}}({{u}}_h)\right)
        \lesssim s_{1}({{u}}_h,{{u}}_h)
      \end{aligned}$$ and hence we arrive at $$\begin{aligned}
      \nonumber
        I 
        \gtrsim  
      - s_{1}({{u}}_h,{{u}}_h)^{1/2} \|{{\epsilon}}( {{u}}_h) \|_{\Omega}
      - s_{1}\left( {{u}}_h, {{u}}_h \right)^{1/2}
        \|  \Pi^{\ast}_h {{\epsilon}}( {{u}}_h)  \|_{{\Omega^{\ast}}}. 
      \end{aligned}$$\
[ Estimate (II)]{}: Here we use inverse estimate  and the Nitsche penalty to control the boundary contribution: $$\begin{aligned}
      \nonumber
      {\left( \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}}, {{u}}_h \right)}_{\Gamma}
            &= 
            {\left(  h^{1/2} \Pi^{\ast}_h {{\epsilon}}({{u}}_h) \cdot {{n}},
            h^{-1/2} {{u}}_h \right)}_{\Gamma}
            \\
            \nonumber
            &\lesssim \| \Pi^{\ast}_h  {{\epsilon}}({{u}}_h) \|_{{\Omega^{\ast}}}
            {\left( h^{-1} {{u}}_h, {{u}}_h \right)}^{1/2}_{\Gamma}
            \\
            \nonumber
            &\lesssim 
            \| {{\epsilon}}({{u}}_h) \|_{\Omega}
            {\left( h^{-1} {{u}}_h, {{u}}_h \right)}^{1/2}_{\Gamma}. 
    \end{aligned}$$ Collecting the estimates for $I$ and $II$ gives $$\begin{aligned}
      \nonumber
      a_h({{u}}_h, {{{ \tau}}}_h)
      &\gtrsim
      \| {{u}}_h\|_{1,\Omega}
      \| {{\epsilon}}({{u}}_h) \|_{\Omega}
      - s_{1}({{u}}_h,{{u}}_h)^{1/2} \|{{\epsilon}}( {{u}}_h) \|_{\Omega}
      \\
      &\quad - s_{1}\left( {{u}}_h, {{u}}_h \right)^{1/2}
        \| {{{ \tau}}}_h \|_{{\Omega^{\ast}}}
      +  \| {{\epsilon}}({{u}}_h) \|_{\Omega}
      {\left( h^{-1} {{u}}_h, {{u}}_h \right)}^{1/2}_{\Gamma}. 
      \label{eq:a_h-final-estimate}
    \end{aligned}$$ Finally, we divide  by $\| {{\epsilon}}({{u}}_h) \|_{\Omega}$ and recall $\| {{{ \tau}}}_h \|_{{\Omega^{\ast}}} \lesssim \| {{\epsilon}}({{u}}_h)
    \|_{\Omega}$ to find that for some $c > 0$ $$\begin{aligned}
      \nonumber
    \dfrac{a_h({{{ \tau}}}_h, \, {{u}}_h)}{\| {{{ \tau}}}_h \|_{{\Omega^{\ast}}}}
 \gtrsim 
    \dfrac{a_h({{{ \tau}}}_h, \, {{u}}_h)}{\| {{\epsilon}}({{u}}_h) \|_{\Omega}}
 \gtrsim
    \| {{u}}_h \|_{1,\Omega} -  c ( s_{1}({{u}}_h,{{u}}_h)^{1/2} +
    {\left( h^{-1} {{u}}_h, {{u}}_h \right)}^{1/2}_{\Gamma})
    \label{eq:mod-inf-sup-ah}
  \end{aligned}$$ if $a_h({{{ \tau}}}_h, \, {{u}}_h) \geqslant 0$, otherwise we proceed as in the previous proof. 

Combining the modified inf-sup conditions  and  enable us to prove an inf-sup condition for the discrete variational problem  with respect to the total approximation space ${\mathbb{V}}_h$.

It holds $$\sup_{V_h \in {\mathbb{V}}_h\setminus \{0\}} 
    \dfrac{A_h(U_h,V_h) + S_h(U_h,V_h)}{{|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h}
    \gtrsim
    {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h, \quad {\forall\,}U_h \in {\mathbb{V}}_h. 
    \label{eq:inf-sup-Ah-Sh}$$ \[thm:inf-sup-Ah-Sh\]

Given $U_h = ({\sigma}_h, {{u}}_h, p_h)$, we construct a proper test function $V_h$ in four steps.\
[**I**]{}: Choosing $V^1_h = U_h$, we obtain $$\begin{aligned}
    \label{eq:inf-sup-A1-term}
    A_h(U_h,V_h^1)  + S_h(U_h, V_h^1)
    = \dfrac{1}{2 \eta} \| {\sigma}_h \|_{\Omega}^2
    + 2 \eta \gamma_b \|h^{-1/2}{{u}}_h \|_{\Gamma}^2
    + S_h(U_h, U_h). 
  \end{aligned}$$\
[**II**]{}: Now we choose $V_h^2 = (0, {{v}}_h^p,0)$ , where ${{v}}_h^p$ attains the supremum in  for given $p_h$ and is rescaled such that $\eta \| {{v}}_h^p \|_{1, h}^2 = \tfrac{1}{\eta} \| p_h
  \|_{\Omega}^2$. Then writing $A_2 = A_h(U_h,V_h^2) + S_h(U_h,V_h^2)$ and applying Cauchy-Schwarz and the modified inf-sup condition , we obtain $$\begin{aligned}
    A_2 
    &= a_h(\sigma_h, {{v}}_h^p) + b_h(p_h, {{v}}_h^p) + s_u({{u}}_h,{{v}}_h^p)  + \frac{\gamma_b \eta}{h} 
    {\left( 
    u_h,\,v_h^p \right)}_{\Gamma} 
    \nonumber
    \\
    &\gtrsim
    - \| {\sigma}_h \|_{\Omega} \| {{\epsilon}}({{v}}_h^p) \|_{\Omega} 
    - \| h^{1/2} {\sigma}_h \|_{\Gamma} \| h^{-1/2} {{v}}_h^p \cdot n \|_{\Gamma}
    \nonumber
    \\
    &\quad 
    + \| p_h \|_{\Omega}  \| {{v}}_h^p \|_{1,h} 
    - s_p(p_h, p_h)^{1/2} \| {{v}}_h^p \|_{1,h} 
    - \dfrac{\eta \gamma_b}{h} \|{{u}}_h \|_{\Gamma} \| {{v}}_h \|_{\Gamma}
    + s_u({{u}}_h, {{v}}_h^p)
    \nonumber
    \\
    &\gtrsim
    - \dfrac{\delta^{-1}}{\eta} \| {\sigma}_h \|_{\Omega}^2
    -        \delta \eta \| {{\epsilon}}({{v}}_h^p) \|_{\Omega}^2 
    - \dfrac{\delta^{-1}}{\eta} \| h^{1/2} {\sigma}_h \|_{\Gamma}^2
    -        \delta \eta \| h^{-1/2} {{v}}_h^p \cdot n \|_{\Gamma}^2
    \nonumber
    \\
    &\quad 
    + \dfrac{1}{\eta}\| p_h \|_{\Omega}^2 
    - \dfrac{\delta^{-1}}{\eta} s_p(p_h, p_h) - \delta \eta  \| {{v}}_h^p
    \|_{1,h}^2
    \nonumber
  - \delta^{-1} \eta \gamma_b \|h^{-1/2} {{u}}_h \|_{\Gamma}^2
  - \delta \eta \gamma_b \|h^{-1/2} {{v}}_h^p \|_{\Gamma}^2
    \nonumber
    \\
    &\quad 
    - \dfrac{\delta^{-1}}{\eta} s_u({{u}}_h, {{u}}_h)
    - \delta \eta s_u({{v}}_h^p, {{v}}_h^p),
  \end{aligned}$$ where a $\delta$-weighted arithmetic-geometric inequality was used in the last step.  Due to the scaling choice $\eta \| {{v}}_h^p \|_{1,h}^2 = \tfrac{1}{\eta} \| p_h
  \|_{\Omega}^2$, all ${{v}}_h^p$-related terms can be absorbed into $\| p_h \|_{\Omega}$ by choosing $\delta$ small enough.  If we then combine an inverse estimate and Lemma \[lem:control-norm-with-ghost-penalties\] to estimate the boundary term $\| h^{1/2} {\sigma}_h \|_{\Gamma}$ by $$\begin{aligned}
    \dfrac{1}{\eta}\| h^{1/2} {\sigma}_h \|_{\Gamma}^2 
    \lesssim 
    \dfrac{1}{ \eta}
    \| 
    {\sigma}_h
    \|_{{\Omega^{\ast}}}^2
    \lesssim  
    \dfrac{1}{ \eta}
    \| {\sigma}_h \|_{\Omega}^2 + 
    s_{{\sigma}}({\sigma}_h,
    {\sigma}_h),
  \end{aligned}$$ we see that there exists constants $C_p$ such that $$\begin{aligned}
    A_2 
    &\gtrsim
    \dfrac{1}{\eta}\| p_h \|_{{\Omega^{\ast}}}^2 
    - C_p 
    \left(
    \dfrac{1}{\eta} \| {\sigma}_h \|_{\Omega}^2
    + \eta \gamma_b \|h^{-1/2}{{u}}_h \|_{\Gamma}^2
    + S_h(U_h, U_h)
    \right),
    \label{eq:inf-sup-A2-term}
  \end{aligned}$$ where we also applied Lemma \[lem:control-norm-with-ghost-penalties\] to pass from $\| p_h \|_{\Omega}$ to $\| p_h \|_{{\Omega^{\ast}}}$ via the term $s_p(p_h, p_h)$.

[**III**]{}: Next, we pick $V_h^3 = ({{{ \tau}}}_h^u, 0 ,0)$, where ${{{ \tau}}}^u_h$ attains the supremum in  for the given ${{u}}_h$ and is rescaled such that $ \dfrac{1}{\eta} \| {{{ \tau}}}_h^u
  \|_{{\Omega^{\ast}}}^2 =  \eta \| {{u}}_h \|_{1,\Omega}^2$. Introducing $A_3 =  A_h(U_h,V_h^3) +
  S_h(U_h,V_h^3)$, we can bound $A_3$ along the same lines as in the previous step: $$\begin{aligned}
    A_3 
    &= \dfrac{1}{2 \eta} ({\sigma}_h, {{{ \tau}}}_h^p)_{\Omega}
    - a_h({{{ \tau}}}_h^p, {{u}}_h)_{\Omega}
    + s_{{\sigma}}({\sigma}_h, {{{ \tau}}}_h^u)
    \nonumber
    \\
    & \gtrsim
    - \dfrac{\delta^{-1}}{\eta} \| {\sigma}_h \|_{\Omega}^2
    - \dfrac{\delta}{ \eta} \| {{{ \tau}}}_h^u \|_{\Omega}^2
    + \eta \| {{u}}_h \|_{1,\Omega}^2 
    - \delta^{-1} s_u({{u}}_h, {{u}}_h)
    - \delta^{-1} \eta \| h^{-1/2} {{u}}_h \|_{\Gamma}^2
    \nonumber
    \\
    &\quad 
    - \delta^{-1} s_{\sigma}({\sigma}_h, {\sigma}_h)
    - \delta s_{\sigma}({{{ \tau}}}_h^u, {{{ \tau}}}_h^u)
    \nonumber
    \\
    & \gtrsim
    \eta \| {{u}}_h \|_{1,{\Omega^{\ast}}}^2 -
     C_u 
     \left(
    \dfrac{1}{\eta} \| {\sigma}_h \|_{\Omega}^2
    + \eta \gamma_b \|h^{-1/2}{{u}}_h \|_{\Gamma}^2
    + S_h(U_h, U_h)
     \right). 
    \label{eq:inf-sup-A3-term}
  \end{aligned}$$ [**IV**]{}: Finally, we define $V_h = V_h^1 + \alpha V_h^2 + \beta
  V_h^3$. Combining the estimates , and , we observe that by choosing $\alpha$ and $\beta$ small enough, it holds that $$\begin{aligned}
    A_h(U_h,V_h) + S_h(U_h,V_h)
    \gtrsim {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h^2. 
    \nonumber
  \end{aligned}$$ To conclude the proof, we note that by our choices of $V_h^i$, $i =
  1,2,3$, $$\begin{aligned}
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h^1 {|\mspace{-1mu}|\mspace{-1mu}|}_h^2 &= {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h^2,
    \\
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h^2 {|\mspace{-1mu}|\mspace{-1mu}|}_h^2 &= 2 \eta \| \epsilon(v_h^p) \|_{\Omega}^2 + 2 \eta \gamma_b
    \|h^{-1/2} v_h^p \|_{\Gamma}^2 + s_u(v_h^p, v_h^p)
    \lesssim 
    \eta \| v_h^p \|_{1,h}^2 = \dfrac{1}{\eta}
    \| p_h \|^2_{\Omega},
    \\
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h^3 {|\mspace{-1mu}|\mspace{-1mu}|}_h^2 &= \dfrac{1}{2 \eta} {|\mspace{-1mu}|\mspace{-1mu}|}\tau_h^u {|\mspace{-1mu}|\mspace{-1mu}|}_{\Omega}^2 +
    s_{\sigma}(\tau_h^u, \tau_h^u)
    \lesssim \dfrac{1}{\eta} \| \tau_h^u \|_{{\Omega^{\ast}}}
    = \eta \| u_h \|^2_{1,\Omega},
  \end{aligned}$$ and thus ${|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h \lesssim {|\mspace{-1mu}|\mspace{-1mu}|}U_h {|\mspace{-1mu}|\mspace{-1mu}|}_h$ which proves the desired estimate. 

A priori estimates {#sec:apriori-analysis}
==================

In this section, we state and prove the a priori estimate for the error in the discrete solution, defined by problem .  Before we present the main result, we state two lemmas which quantify the effect of the consistency error introduced by the stabilization term $S_h$. The first lemma ensures that a weakened form of the Galerkin orthogonality holds:

Let $(\sigma_h, u_h, p_h) \in {\mathbb{V}}_h$ be the finite element approximation defined by  and assume that the weak solution $({\sigma}, {{u}}, p)$ of the three field Stokes problem  is in $[H^1(\Omega)]^{d\times d} 
  \times [H^2_0(\Omega)]^d \times H^1(\Omega)$. Then $$\begin{aligned}
    A_h(U - U_h, V_h) =  S(U_h, V_h). 
  \end{aligned}$$

The proof follows immediately from the definition of the weak variational problem  and the easily verified fact that the continuous solution $U$ satisfies $A_h(U,V_h) = L_h(V_h)$.

The second lemma ensures that the consistency error does not make the convergence rate deteriorate. 

Suppose that $U = ({\sigma}, {{u}}, p) \in [H^1(\Omega)]^{d \times d} \times
  [H^2(\Omega)]^d \times H^1(\Omega)$, then it holds that $$\begin{aligned}
    \label{eq:consistency-error}
    | S_h({\mathcal{C}}_h^* U, V_h) | \lesssim h
    \left(
    \eta^{1/2}\|{{u}}\|_{2,\Omega} 
    +\frac{1}{\eta^{1/2}}\|p\|_{1,\Omega}
    +\frac{1}{\eta^{1/2}}\|{\sigma}\|_{1,\Omega}
    \right)
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h. 
  \end{aligned}$$

By definition, $$\begin{aligned}
    S_h({\mathcal{C}}_h^* U, V_h) 
    = s_{\sigma}({\mathcal{C}}_h^* {\sigma}, {{{ \tau}}}_h) 
    + s_u({\mathcal{C}}_h^* u, v_h) 
    + s_p({\mathcal{C}}_h^* p, q_h). 
  \end{aligned}$$ We start with the velocity related terms.  Since we assume that $u \in H^2(\Omega) \cap H^1_0(\Omega)$, we have $s_u({{u}}^{\ast},{{v}}_h) = 0$ for its extension $u^{\ast} =
  Eu$ to ${\Omega^{\ast}}$. Exploiting this fact together with the trace inequality , the inverse estimate , the interpolation estimate  and the stability of the interpolation operator ${\mathcal{C}}_h^*$, we might estimate the velocity part of the consistency error as follows: $$\begin{aligned}
    | s_u({\mathcal{C}}_h^* {{u}}, {{v}}_h) |
    &=  
    | s_u({\mathcal{C}}_h^* {{u}}- {{u}}^{\ast}, {{v}}_h) |
    \lesssim
    \eta \sum_{F\in {\mathcal{F}}_i} 
    h^{1/2}\| {\partial_{{{n}}}}({\mathcal{C}}_h^* {{u}}- {{u}}^{\ast}) \|_{F}
    \,
    h^{1/2}\| {\partial_{{{n}}}}{{v}}_h \|_{F}
    \\
    & \lesssim
    \eta^{1/2}  \left(\sum_{T\in {\mathcal{T}_h}} 
    \left(h
    \| \nabla ({\mathcal{C}}_h^* {{u}}- {{u}}^{\ast}) \|^2_{T}
    + \|  ({\mathcal{C}}_h^* {{u}}- {{u}}^{\ast}) \|^2_{T}
  \right) \right)^{\frac12}
    \eta^{1/2} \| \nabla {{v}}_h \|_{{\Omega^{\ast}}}
    \\
    &\lesssim h \eta^{1/2} \| {{u}}^{\ast} \|_{2,{\Omega^{\ast}}} \| {{v}}_h \|_{1,{\Omega^{\ast}}}
    \lesssim  h \eta^{1/2} \| {{u}}\|_{2,\Omega} {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h. 
  \end{aligned}$$

For the pressure, applying the inverse inequality  and the boundedness of the interpolation operator  gives $$\begin{aligned}
    | s_p({\mathcal{C}}_h^* p, q_h ) |
    \lesssim
    h \eta^{-1} \| \nabla {\mathcal{C}}_h^* p \|_{{\Omega^{\ast}}}  
    h \| \nabla q_h \|_{{\Omega^{\ast}}}
    \lesssim
    h \eta^{-1} \| p \|_{1,\Omega}  
    \| q_h \|_{{\Omega^{\ast}}}
    \lesssim
    h \eta^{-1/2} \| p \|_{1,\Omega}  
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h. 
  \end{aligned}$$ Finally, we observe that the consistency error in ${\sigma}_h$ might be bounded by applying the same steps as for the pressure related terms. 

We are now in the position to state our main result. 

Let $U = ({\sigma},{{u}},p) \in [H^1(\Omega)]^{d \times d} \times
  [H^2(\Omega)]^d \times H^1(\Omega)$ be the solution of the three field Stokes problem  and let $U_h = ({\sigma}_h, {{u}}_h, p_h)$ be the solution to the discrete problem .  Then the following error estimate holds: $${\interleave U-U_h\interleave}\lesssim h 
  \left(
    \eta^{1/2}\|{{u}}\|_{2,\Omega} +
    \frac{1}{\eta^{1/2}}\|p\|_{1,\Omega} + \frac{1}{\eta^{1/2}}
    \|{\sigma}\|_{1,\Omega}
  \right),$$ where the hidden constant is independent of how the boundary cuts the mesh. 

Using the triangle inequality ${|\mspace{-1mu}|\mspace{-1mu}|}U - U_h {|\mspace{-1mu}|\mspace{-1mu}|}\lesssim {|\mspace{-1mu}|\mspace{-1mu}|}U -
  {\mathcal{C}}_h^* U {|\mspace{-1mu}|\mspace{-1mu}|}+ {|\mspace{-1mu}|\mspace{-1mu}|}U_h - {\mathcal{C}}_h^* U {|\mspace{-1mu}|\mspace{-1mu}|}_h$ and the standard interpolation estimates , we can see that the error ${|\mspace{-1mu}|\mspace{-1mu}|}U - {\mathcal{C}}_h^* U {|\mspace{-1mu}|\mspace{-1mu}|}$ satisfies the desired estimate.  By inf-sup condition  and the weak Galerkin orthogonality, there exists a $V_h$ such that $$\begin{aligned}
    {|\mspace{-1mu}|\mspace{-1mu}|}U_h -  {\mathcal{C}}_h^* U {|\mspace{-1mu}|\mspace{-1mu}|}_h
    &\lesssim 
    \dfrac{A_h(U_h -  {\mathcal{C}}_h^* U, V_h) + S_h(U_h - {\mathcal{C}}_h^* U, V_h) }{{|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h}
    \\
    &= 
    \dfrac{A_h(U -  {\mathcal{C}}_h^* U, V_h) - S_h({\mathcal{C}}_h^* U, V_h) }{{|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h}
    = A + S.
    \label{}
  \end{aligned}$$ Recalling the bound for the consistency error , it suffices to estimate $$\begin{aligned}
    \nonumber
    A &= 
    \dfrac{1}{2 \eta} ({\sigma}- {\mathcal{C}}_h^* {\sigma}, {{{ \tau}}}_h)_{\Omega}
    + 2 \eta \gamma_b {\left( h^{-1} ({{u}}-  {\mathcal{C}}_h^* {{u}}),
    {{v}}_h \right)}_{\Gamma}
    + a_h({\sigma}- {\mathcal{C}}_h^* {\sigma}, {{v}}_h)
    - a_h({{{ \tau}}}_h, {{u}}- {\mathcal{C}}_h^* {{u}}) 
    \\
    \nonumber
    &\quad
    + b_h(p - {\mathcal{C}}_h^* p, {{v}}_h)
    - b_h(q_h, {{u}}- {\mathcal{C}}_h^* {{u}}).  
  \end{aligned}$$ For the first term, we simply have $$\begin{aligned}
   |\dfrac{1}{2 \eta} ({\sigma}- {\mathcal{C}}_h^* {\sigma},
   {{{ \tau}}}_h)_{\Omega}|
   &\lesssim
   \dfrac{h}{\eta^{1/2}} \| {\sigma}\|_{1,\Omega} {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h,
  \end{aligned}$$ while for the second term, combining the trace inequality  with the interpolation estimate  yields $$\begin{aligned}
    2 \eta \gamma_b {\left( h^{-1} ({{u}}-  {\mathcal{C}}_h^* {{u}}),
    {{v}}_h \right)}_{\Gamma}
    \lesssim (2 \eta \gamma_b)^{1/2} h \|{{u}}\|_{2,\Omega} {|\mspace{-1mu}|\mspace{-1mu}|}V_h
    {|\mspace{-1mu}|\mspace{-1mu}|}_h. 
  \end{aligned}$$ Next, the third term can be estimated by $$\begin{aligned}
   |a_h({\sigma}- {\mathcal{C}}_h^* {\sigma}, {{v}}_h)|
   &\lesssim
   \dfrac{1}{\eta^{1/2}} 
   ( \| {\sigma}- {\mathcal{C}}_h^* {\sigma}\|_{\Omega}
   + \|h^{1/2}({\sigma}- {\mathcal{C}}_h^* {\sigma})\|_{\Gamma}
   )
   \cdot \eta^{1/2}
   (\| {{\epsilon}}({{v}}_h) \|_{\Omega}
    + \| h^{-1/2} {{v}}_h \|_{\Gamma}
    )
    \\
    &\lesssim
    \dfrac{h}{\eta^{1/2}} \| {\sigma}\|_{1,\Omega}
    {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h. 
  \end{aligned}$$ Similarly, the fourth term can be bounded $$\begin{aligned}
    | a_h({{{ \tau}}}_h, {{u}}- {\mathcal{C}}_h^* {{u}}) |
    & \lesssim
    \dfrac{1}{\eta^{1/2}}
    (
    \| {{{ \tau}}}_h \|_{\Omega} + \| h^{1/2} {{{ \tau}}}_h \|_{\Gamma}
    )
    \; \eta^{1/2}
    (
    \| {{\epsilon}}({{u}}) - {\mathcal{C}}_h^*{{u}}\|_{\Omega}
    + \| h^{-1/2}({{u}}- {\mathcal{C}}_h^*{{u}}) \|_{\Gamma}
    )
    \nonumber
    \\
    &\lesssim
    h {|\mspace{-1mu}|\mspace{-1mu}|}V_h {|\mspace{-1mu}|\mspace{-1mu}|}_h \eta^{1/2} \| {{u}}\|_{2,\Omega}.  
    \label{eq:apriori-step-3}
  \end{aligned}$$ Here, we estimated the boundary term in  by successively applying the trace inequality , standard interpolation estimates and the boundedness of the extension operator $E: H^2(\Omega) \to H^2({\Omega^{\ast}})$, cf.  , which yields $$\begin{aligned}
    \| h^{-1/2}({{u}}- {\mathcal{C}}_h^*{{u}}) \|_{\Gamma}
    \lesssim
    h^{-1} \| {{u}}^{\ast} - {\mathcal{C}}_h^*{{u}}\|_{{\Omega^{\ast}}}
    + 
    h \| {{u}}^{\ast} - {\mathcal{C}}_h^* {{u}}\|_{1,{\Omega^{\ast}}}
    \lesssim h \| {{u}}\|_{2,\Omega}. 
  \end{aligned}$$ The estimates for the remaining terms involving $b_h(\cdot,\cdot)$ are completely analogous, which concludes the proof. 

To reduce the system matrix stencil, one may use the element-based penalty terms $$\begin{aligned}
s_p(p_h,q_h) &= \frac{\gamma_p}{2 \eta} \sum_{T \in {\mathcal{T}_h}} h^2(\nabla p_h, \nabla q_h)_{T}, \\ 
  s_{\sigma}(\sigma_h,\tau_h) &= \frac{\gamma_{\sigma}}{2 \eta} \sum_{T \in {\mathcal{T}_h}} h^2(\nabla \sigma_h, \nabla \tau_h)_{T} \end{aligned}$$ for the pressure and stress, instead of the face-based penalty terms , over gradient jumps, if linear finite element spaces are chosen for velocity, pressure and stress.  Note that both the face and element-based penalty terms are weakly consistent for ${\sigma}\in [H^1(\Omega)]^{d \times d}$ and $p \in
H^1(\Omega)$. However, the face-based penalty term  for $u \in [H^2(\Omega)]^d$ is strongly consistent and thus strictly necessary as the analogous element-based penalty term leads to a consistency error which deteriorates the overall convergence order. 

Numerical results {#sec: numerical results}
=================

In this section, we will demonstrate that the theoretical estimates of Section \[sec:stability\] and Section \[sec:apriori-analysis\] hold.  In particular, we will show that the finite element solution of velocity, pressure and extra-stress tensor converge with optimal order to a $\sin$-$\cos$ reference solution of the three field Stokes system and we will demonstrate that the ghost penalties yield independence of the quality of the solution on the boundary location.  All numerical simulations have been performed using our software package libCutFEM which will be made available soon at <http://www.cutfem.org>.  LibCutFEM is an open source library which extends the finite element library DOLFIN [@LoggWells2010a] and the FEniCS framework [@LoggMardalEtAl2011] for automated computing of finite element variational problems with cut finite element capabilities.  The inner workings of libCutFEM are described as part of the review article [@BurmanClausHansboEtAl2014]. 

Convergence study for reference solution
----------------------------------------

To evaluate the accuracy of our scheme, we investigate the rate of convergence of the numerical solution to the following reference solution $$\begin{aligned}
{{u}}_{ex} = \left[\begin{matrix} -\sin(\pi y)\cos(\pi x) \\ \sin(\pi x)\cos(\pi y) \end{matrix}\right], \nonumber \\
p_{ex} = -2\eta\cos(\pi x)\sin(\pi y), \nonumber \\
{\sigma}_{ex} = \left[\begin{matrix}2.0 \pi \eta \sin{\left (\pi x \right )} \sin{\left (\pi y \right )} & 0\\0 & - 2.0 \pi \eta \sin{\left (\pi x \right )} \sin{\left (\pi y \right )}\end{matrix}\right], \nonumber \\
{{f}}= \left[\begin{matrix}2 \pi \eta \sin{\left (\pi x \right )} \sin{\left (\pi y \right )} - 2 \pi^{2} \eta \sin{\left (\pi y \right )} \cos{\left (\pi x \right )}\\2 \pi^{2} \eta \sin{\left (\pi x \right )} \cos{\left (\pi y \right )} - 2 \pi \eta \cos{\left (\pi x \right )} \cos{\left (\pi y \right )}\end{matrix}\right] 
\label{equ: 2d sincos solution}\end{aligned}$$ of the three field Stokes system.  Here, we choose $\gamma_u= 0.01$, $\gamma_p= 0.1$, $\gamma_{{\sigma}}=0.1$, $\gamma_b=15.0$ and $\eta= 0.5$ and compute the velocity, pressure and extra-stress in a unit circle embedded in a fixed background mesh.  We set $
{{u}}= {{u}}_{ex} \mbox{ at } \partial \Omega. 
$ For the velocity and the extra-stress tensor, the sum of the error of the components is evaluated.\
The rate of convergence for the $L^2$-error, $|| U_h - U_{ex}||_0$, and for the $H^1$ error of $|| u_h - u_{ex}||_1$ are displayed in Figure \[fig: convergence rate 2d\].  We obtain a convergence order of 1.05 for the velocity in the $H^1$-norm, which is what we expect from our error analysis.  In the $L^2$ norm, the velocity converges with order 2.18.  We obtain a convergence order of 1.77 and 1.99 for the extra-stress and pressure which is better than expected.  However, this can be explained by the smoothness of the solution. 

Stabilizing effect of ghost penalty terms
-----------------------------------------

In this section, we investigate the quality of the solution with respect to how the interface cuts the mesh.  A boundary location in which only very small parts of the elements in the interface zone are covered by the physical domain can lead to an ill-conditioned system matrix and to an unbounded Nitsche boundary penalty parameter $\gamma_b$. To demonstrate that this ill-conditioning and the unboundedness of the penalty parameter can be alleviated using ghost penalties, we investigate the quality of the solution in terms of a boundary location parameter $0<\epsilon<1$. This parameter $\epsilon$ indicates the relative height of thin fluid stripes in a boundary cut parallel to an element edge, see Figure \[subfig: sliverness def\].  We call this type of cut configuration, the sliver case and $\epsilon$ indicates the sliver size. 

### Dependence of the quality of the solution on the sliver size

Consider the reference solution in a square domain $\Omega=[-1,1]^2$ embedded in a dilated background mesh of size $${\Omega^{\ast}}=[-1-l,1+l] \mbox{ with } l=\frac{2(1-\epsilon)} {N-2(1-\epsilon)},$$ where $N$ is the number of elements in the $x$ and in the $y$ direction.  Figure \[subfig: sliver test mesh\] shows the approximated interface location of the quadratic domain in dilated background meshes for $\epsilon=0.5$ and $\epsilon=0.1$. We investigate the effect of the ghost penalty parameter $\gamma_{{\sigma}}$ on the quality of the sin-cos reference solution for $\epsilon=\left\{0.5,0.1,0.02,0.004\right\}$. Throughout this section, we set $\gamma_b=15.0$, $\gamma_{{{u}}}=0.1$ and $\gamma_p=0.1$. Figure \[subfig: conv sliver\] shows that for $\epsilon=0.02$ and $\gamma_{\sigma}=0.1$, the extra-stress, the velocity and the pressure converge with the optimal order of convergence as predicted by the analysis in Section \[sec:apriori-analysis\].  Setting the ghost penalty parameter to $\gamma_{{\sigma}}=0.0$ for $\epsilon=0.02$ causes an upward shift of the error for the extra-stress tensor as shown in Figure \[subfig: conv sliver2\].  Figure \[subfig: stress wo gammas\] shows this increase of the error in the extra-stress tensor with decreasing sliver size $\epsilon$ for the unstabilized extra-stress tensor variable.  Using the ghost penalty stabilization ($\gamma_{{\sigma}}=0.1$), this increase in error can be alleviated and the solution becomes independent of the boundary location (see Figure \[subfig: stress w gammas\]).  The cause for the error increase for the unstabilized extra-stress variable can be observed in Figure \[fig: gammas ghost stab\] for the extra-stress tensor component $\sigma_{xx}$. Without the ghost penalty stabilization, we have huge spikes appearing at the corner of the domain in the solution and the solution shows large oscillations along the boundary.  Even though these spikes and oscillations decrease with mesh refinement the solution of the extra-stress tensor component is polluted by the poor solution in the boundary region.  Setting $\gamma_{{\sigma}}=0.1$ alleviates this problem and the solution does not undergo any large spikes or oscillations in the boundary region. 

\

\

### Condition number

In this section, we investigate the condition number of the system matrix $A$   in dependence to the boundary location for a ghost stabilization parameter of $\gamma_{{\sigma}}=\left\{ 0.0, 0.001,0.1,1.0\right\}$. Here, we consider a fixed fictitious domain ${\Omega^{\ast}}= [-1,1]^2$ with a fixed mesh size $h$ and a shrinking physical domain $\Omega=[-1+(1-\epsilon)\Delta x, 1-(1-\epsilon)\Delta x]^2$, where $\Delta x$ is the edge length of the elements in $x$-direction and $y$-direction.  We choose $\gamma_b=15.0, \gamma_u = 0.1, \gamma_p =
0.1$. Figure \[fig: condition number\] shows the condition number with respect to the sliver size parameter $\epsilon$. We observe that for $\gamma_{{\sigma}}=0.0$, the condition number is unbounded while for $\gamma_{{\sigma}}=\left\{ 0.001,0.1,1.0\right\}$ the condition number is bounded.  Hence, even for very small ghost penalty stabilization parameters the ill-conditioning dependence on the boundary location is alleviated. 

\

Three field Stokes in an aneurysm {#p2:ssec:flow-in-complex-geometry}
---------------------------------

As a final numerical example, we present the computation of a fluid flow governed by the three field Stokes problem in a three-dimensional domain with a complex boundary geometry.  The boundary geometry is taken from a part of an arterial network known as the Circle of Willis which is located close to the human brain.  It is known that the network is prone to develop aneurysms and therefore the computer-assisted study of the blood flow in the Circle of Willis has been a recent subject of interest, see for instance @Steinman2003 [@IsaksenBazilevsKvamsdalEtAl2008; @Valen-Sendstad2011].\
However, the purpose of this example is not to perform a realistic study of the blood flow dynamics.  Rather, we would like to demonstrate the principal applicability of the developed method to simulation scenarios where complex three-dimensional geometries are involved.  The blood vessel geometry is embedded in a structured background mesh as illustrated in Figure \[fig:aneurysm-box-mesh\].\
The velocity is prescribed on the entire boundary $\Gamma$, where we set ${{u}}= 0$ on the arterial walls and ${{u}}= 1200\,
\mathrm{mm/s}$ on the inlet boundary.  The two outflow velocities are set such that the total flux is balanced.  We choose $\eta= 1.0$, $\gamma_{{{u}}}= 0.1$, $\gamma_p= 0.1$, $\gamma_{{\sigma}}=0.1$ and $\gamma_b=10.0$.\
Figure \[fig: aneurysm\] displays the pressure, velocity and extra-stress profiles in the aneurysm geometry.  The extra-stress tensor is displayed in terms of the van Mises stress measure [@Mises1913] given by $$\sigma_{v}^2 = \frac{1}{2} \left[ \left( \sigma_{xx} - \sigma_{yy}\right)^2 + \left( \sigma_{yy} - \sigma_{zz}\right)^2 + \left( \sigma_{zz} - \sigma_{xx}\right)^2  + 6 \left( \sigma_{xy}^2 + \sigma_{yz}^2 + \sigma_{xz}^2\right)\right].$$ This stress measure provides an indication of the strength of normal stress differences and shear stresses in the fluid.  Although the fictitious domain mesh ${\mathcal{T}_h}$ provides only a coarse resolution of the aneurysm geometry, the values of the velocity approximation clearly conforms to the required boundary values on the actual surface geometry. 

\

Conclusions
===========

In this article, we have developed a novel fictitious domain method for the three field Stokes equation.  We have demonstrated theoretically and numerically that our scheme is inf-sup stable and possesses optimal convergence order properties independent of the boundary location.  We have approximated the velocity, pressure and extra stress tensor with linear finite elements and we have stabilized our scheme using a continuous interior penalty approach combined with ghost penalty terms in the boundary region.  We have demonstrated that the ghost penalties in the boundary region guarantee a stable and accurate solution independent of how the boundary intersects the mesh.  Additionally, we have demonstrated numerically that the ghost penalty stabilization yields a bounded condition number independent of the boundary location.\
In a future contribution, we will extend the scheme developed in this paper to multi-phase three field Stokes problems. 

Acknowledgements {#acknowledgements .unnumbered}
================

The work for this article was supported by the EPSRC grant EP/J002313/2 on “Computational methods for multiphysics interface problems” and a Center of Excellence grant from the Research Council of Norway to the Center for Biomedical Computing at Simula Research Laboratory.  The authors wish to thank Sebastian Warmbrunn for providing the surface geometry used in Section \[p2:ssec:flow-in-complex-geometry\] and the anonymous referees for helpful comments.