remarkRemark hypothesisHypothesis claimClaim Analysis of a DG Method for the
BP problemI. P. A. Papadopoulos

# Numerical analysis of a discontinuous Galerkin method for the
Borrvall–Petersson topology optimization problem ††thanks: Submitted DATE. The
author is supported by the EPSRC Centre for Doctoral Training in Partial
Differential Equations: Analysis and Applications [grant number EP/L015811/1]
and The MathWorks, Inc.

Ioannis P. A. Papadopoulos Mathematical Institute, University of Oxford,
Oxford, UK<EMAIL_ADDRESS>
###### Abstract

Divergence-free discontinuous Galerkin (DG) finite element methods offer a
suitable discretization for the pointwise divergence-free numerical solution
of Borrvall and Petersson’s model for the topology optimization of fluids in
Stokes flow [Topology optimization of fluids in Stokes flow, International
Journal for Numerical Methods in Fluids 41 (1) (2003) 77–107]. The convergence
results currently found in literature only consider $H^{1}$-conforming
discretizations for the velocity. In this work, we extend the numerical
analysis of Papadopoulos and Süli to divergence-free DG methods with an
interior penalty [I. P. A. Papadopoulos and E. Süli, Numerical analysis of a
topology optimization problem for Stokes flow, arXiv preprint
arXiv:2102.10408, (2021)]. We show that, given an isolated minimizer of the
infinite-dimensional problem, there exists a sequence of DG finite element
solutions, satisfying necessary first-order optimality conditions, that
strongly converges to the minimizer.

###### keywords:

topology optimization, nonconvex variational problem, multiple solutions,
finite element method, discontinuous Galerkin method

35Q35, 65K10, 65N30, 90C26

##  1 Introduction

The Borrvall–Petersson problem is the first model in literature for the
topology optimization of fluid flow [13]. Given a restriction on the
proportion of a design domain that a Stokes fluid can occupy, the problem
attempts to find the optimal channels from the inlets to the outlets that
minimize the power dissipation of the flow. No prior assumptions are required
for the shape or topology of the solution resulting in a flexible optimization
process [4, 11]. This flexibility comes at the cost of a nonconvex
optimization problem with PDE, box, and inequality constraints that often
supports multiple minimizers.

The topology of the solution is encoded in the material distribution,
$\rho:\Omega\to[0,1]$; a function that maps from the coordinates of the domain
to the unit interval. Regions where $\\{\rho=1\\}\subset\Omega$ are physically
interpreted as the channels in which the fluid flows through. Whereas the
regions $\\{\rho=0\\}\subset\Omega$ are highly impermeable and thus there is
negligible fluid flow in these areas. Intermediate values
$\\{0<\rho<1\\}\subset\Omega$ are difficult to interpret. However, these
regions are penalized in the Borrvall–Petersson model via an inverse
permeability function $\alpha$.

Since Borrvall and Petersson’s seminal work, there have been numerous
extensions. Evgrafov [22], Olesen et al. [39], and Gersborg-Hansen et al. [31]
extended the model to fluids satisfying the steady-state Navier–Stokes flow.
Kreissl et al. [36] and Deng et al. [20] were the first to consider unsteady
Navier–Stokes flow and Deng et al. [19] later included body forces. Alonso et
al. extended the model to rotating bodies in cylindrical coordinates [5, 6].
For a detailed review on the literature of the topology optimization of
fluids, we refer to the work of Alexandersen and Andreasen [3] .

Due to the nonlinear nature of the problem, solutions of the infinite-
dimensional problem are often difficult to find. Hence, the problem is
typically discretized with the finite element method and the minimizers are
computed numerically. The first result concerning finite element convergence
can be found in the original paper by Borrvall and Petersson [13]. They showed
that a conforming and inf-sup stable finite element discretization for the
velocity-pressure pair together with a piecewise constant discretization for
the material distribution weakly(-*) converges to an unspecified minimizer of
the problem. Thore [50] proved a similar result for a low-order finite element
approximation for the velocity-pressure pair that utilizes a penalty on the
jumps of the pressure to overcome the violation of the inf-sup stability.
Recently, the original Borrvall–Petersson result was improved and extended by
Papadopoulos and Süli [42] for conforming and inf-sup stable finite element
discretizations for the velocity-pressure pair and conforming discretizations
for the material distribution. They showed that for _every_ isolated minimizer
of the problem, there exists a sequence of finite element solutions, to the
first-order optimality conditions, that _strongly_ converges to the minimizer.
Their analysis resolved a number of outstanding issues; namely the lack of
checkerboarding in material distribution approximation as the mesh size tends
to zero and whether every minimizer of the problem could be well approximated
by a finite element method. The difficulty in the analysis is primarily due to
the nonconvexity of problem. In order to account for the possibility of
multiple minimizers, the authors fix a minimizer of the infinite-dimensional
problem and construct a modified optimization problem with the fixed minimizer
as its unique minimizer. Then, strong convergence of the finite element
minimizers to the modified optimization problem is proven. The modified
optimization problem is then related back to the original optimization problem
by showing that a subsequence of the finite element minimizers also satisfy
the first-order optimality conditions of the original optimization problem.

The analysis by Borrvall and Petersson [13], Thore [50], and Papadopoulos and
Süli [42] heavily relied on the use of an $H^{1}$-conforming finite element
discretization for the velocity. In the past couple of decades, discontinuous
Galerkin (DG) methods for fluid flow have become increasingly popular [18, 17,
30, 34, 35]. This is in part due to the existence of divergence-free DG finite
element methods. Some stable finite element methods for fluid flow, such as
the Taylor–Hood finite element pair, do not satisfy the incompressibility
constraint, $\mathrm{div}(\boldsymbol{u})=0$, pointwise. This manifests as a
dependence of the error in the velocity on the best approximation error in the
pressure. In some problems, pointwise violation of the incompressibility
constraint has been observed to support instabilities that result in
nonphysical solutions [33, 37]. In divergence-free finite element methods, the
incompressibility constraint is satisfied pointwise which is useful for
ensuring pressure robustness [49] and deriving error bounds on the velocity
that are independent of the error of the pressure.

In Borrvall–Petersson topology optimization problems, a natural mesh
refinement to obtain sharper solutions is in regions where $0<\rho<1$ a.e. It
can be empirically checked that mesh refinement in these regions does little
to improve the error in the pressure. If the convergence for the velocity and
material distribution rely heavily on the convergence of the pressure, then
only doing mesh refinement in those areas caps the improvement in the errors
for the velocity and material distribution. This motivates the need for
discretizations that decouple the dependence of errors of the velocity and
material distribution with the approximation error of the pressure.
Divergence-free finite element discretizations also allow for an easier
characterization of the kernel of the discretized grad-div term. This
characterization has applications in preconditioners for systems arising in
incompressible fluid flow [26, 27, 32, 40, 48].

$H^{1}$-conforming divergence-free finite element methods exist, for example
the
Scott–Vogelius finite element [49]. To ensure inf-sup stability for a general
mesh in a $k$-th order Scott–Vogelius finite element method, the polynomial
order for the velocity space must be $k\geq 2d$ where $d\in\\{2,3\\}$ is the
dimension of the problem [49, 51]. The expense of the high order method is
normally justified by the accompanying high convergence rate. However, the
material distribution $\rho$ is often discretized with piecewise constant or
continuous piecewise linear finite elements due to the box constraints on the
material distribution. The box constraints not only cause algorithmic
restrictions but also reduce the regularity of $\rho$. The relatively low
order approximation of the material distribution then caps the order of
convergence of the velocity and pressure [42, Sec. 5.1] which negates the
advantage of the high order method. Inf-sup stability can be achieved for
$k\geq d$ if the mesh is barycentrically refined [44]. This was successfully
implemented for the double-pipe problem [41, Sec. 4.1], an example of the
Borrvall–Petersson problem with a local and global minimizer. However,
barycentrically refined meshes can be difficult to align with jumps in the
material distribution that solves the infinite-dimensional problem, which can
lead to poorly resolved solutions. Moreover, barycentrically refined meshes
complicate the generation of a mesh hierarchy for robust multigrid cycles
[26]. In contrast, there exist low-order divergence-free DG finite element
methods that are inf-sup stable on general meshes.

In this paper, we extend the results of Papadopoulos and Süli [42] to
divergence-free DG finite element methods. A standard interior penalty method
is used to control the jumps across the facets. Our main result is to show
that for every isolated minimizer of the Borrvall–Petersson optimization
problem, there exists a sequence of DG finite element solutions to the
discretized first-order optimality conditions that strongly converges in the
appropriate norms to the minimizer. In particular, if $(\boldsymbol{u},\rho)$
is an isolated velocity-material distribution minimizer to the
Borrvall–Petersson problem, and $p$ is the associated pressure, then there
exists a sequence of strongly converging finite element solutions
$(\boldsymbol{u}_{h},p_{h},\rho_{h})$ such that
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$, $\|p-p_{h}\|_{L^{2}(\Omega)}\to 0$, and
$\|\rho-\rho_{h}\|_{L^{s}(\Omega)}\to 0$, $s\in[1,\infty)$, where
$\|\cdot\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}$ is the broken
$H^{1}_{\boldsymbol{g}}$-norm as defined in Eq. 18. This analysis ensures that
_every_ isolated minimizer is well approximated by the divergence-free DG
finite element method as the mesh size tends to zero.

##  2 Topology optimization of Stokes flow

Given a volume constraint on a fluid in a fixed bounded Lipschitz domain
$\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, the Borrvall–Petersson model
attempts to minimize the energy lost by the flow due to viscous dissipation,
whilst maximizing the flow velocities at the applied body force. More
precisely, the objective is to find $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ that
minimizes

(BP) $\displaystyle
J(\boldsymbol{u},\rho)\coloneqq\frac{1}{2}\int_{\Omega}\left(\alpha(\rho)|\boldsymbol{u}|^{2}+\nu|\nabla\boldsymbol{u}|^{2}-2\boldsymbol{f}\cdot\boldsymbol{u}\right)\text{d}x,$

where $\boldsymbol{u}:\Omega\to\mathbb{R}^{d}$ denotes the velocity of the
fluid, $\rho:\Omega\to[0,1]$ is the material distribution of the fluid,
$|\cdot|$ denotes the Euclidean norm for $\boldsymbol{u}$ and the Frobenius
norm for $\nabla\boldsymbol{u}$, and

$\displaystyle H^{1}_{\boldsymbol{g}}(\Omega)^{d}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in
H^{1}(\Omega)^{d}:\boldsymbol{v}=\boldsymbol{g}\;\text{on}\;\partial\Omega\\},$
$\displaystyle H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in
H^{1}_{\boldsymbol{g}}(\Omega)^{d}:\mathrm{div}(\boldsymbol{v})=0\;\text{a.e.\
in}\;\Omega\\},$ $\displaystyle C_{\gamma}$
$\displaystyle\coloneqq\left\\{\eta\in L^{\infty}(\Omega):0\leq\eta\leq
1\;\text{a.e.},\;\;\int_{\Omega}\eta\;\text{d}x\leq\gamma|\Omega|\right\\}.$

Here, $H^{s}(\Omega)$, $0<s<\infty$, and $L^{q}(\Omega)$, $0<q\leq\infty$,
denote the standard Sobolev ($W^{s,2}(\Omega)$) and Lebesgue spaces,
respectively [2]. Furthermore, $\boldsymbol{f}\in L^{2}(\Omega)^{d}$ is a body
force, $\nu>0$ is the (constant) viscosity, and $\gamma\in(0,1)$ is the volume
fraction. The equality $\boldsymbol{u}=\boldsymbol{g}$ on $\partial\Omega$ is
to be understood in the boundary trace sense [21, Sec. 5.5]. Moreover, the
boundary data $\boldsymbol{g}\in H^{1/2}(\partial\Omega)^{d}$ and
$\boldsymbol{g}=\boldsymbol{0}$ on $\Gamma\subset\partial\Omega$, with
$\mathcal{H}^{d-1}(\Gamma)>0$, i.e. $\Gamma$ has nonzero Hausdorff measure on
the boundary. Borrvall and Petersson introduced the inverse permeability term,
$\alpha$, which models the influence of the material distribution on the flow.
For values of $\rho$ close to one, $\alpha(\rho)$ is small, permitting fluid
flow; for values of $\rho$ close to zero, $\alpha(\rho)$ is very large,
restricting fluid flow. The function $\alpha$ satisfies the following
properties:

  1. (A1)

$\alpha:[0,1]\to[\underline{\alpha},\overline{\alpha}]$ with
$0\leq\underline{\alpha}<\overline{\alpha}<\infty$;

  2. (A2)

$\alpha$ is strongly convex and monotonically decreasing;

  3. (A3)

$\alpha(0)=\overline{\alpha}$ and $\alpha(1)=\underline{\alpha}$;

  4. (A4)

$\alpha$ is twice continuously differentiable,

generating an operator also denoted $\alpha:C_{\gamma}\to
L^{\infty}(\Omega;[\underline{\alpha},\overline{\alpha}])$. Typically, in the
literature $\alpha$ takes the form [13, 23]

(1)
$\displaystyle\alpha(\rho)=\bar{\alpha}\left(1-\frac{\rho(q+1)}{\rho+q}\right),$

where $q>0$ is a penalty parameter, so that
$\lim_{q\to\infty}\alpha(\rho)=\bar{\alpha}(1-\rho)$. Borrvall and Petersson
[13, Sec. 3.2] remark that as $q\to\infty$ the material distribution tends to
a 0-1 solution.

We define the following spaces that will be used throughout this work:

(2) $\displaystyle L^{2}_{0}(\Omega)$ $\displaystyle\coloneqq\left\\{q\in
L^{2}(\Omega):\int_{\Omega}q\;\mathrm{d}x=0\right\\},$ (3) $\displaystyle
H^{1}_{0}(\Omega)^{d}$ $\displaystyle\coloneqq\\{\boldsymbol{v}\in
H^{1}(\Omega)^{d}:\boldsymbol{v}=\boldsymbol{0}\;\text{on}\;\partial\Omega\\},$
(4) $\displaystyle\boldsymbol{H}(\mathrm{div};\Omega)$
$\displaystyle\coloneqq\left\\{\boldsymbol{v}\in
L^{2}(\Omega)^{d}:\mathrm{div}(\boldsymbol{v})\in L^{2}(\Omega)\right\\},$ (5)
$\displaystyle\boldsymbol{H}_{0}(\mathrm{div};\Omega)$
$\displaystyle\coloneqq\left\\{\boldsymbol{v}\in\boldsymbol{H}(\mathrm{div};\Omega):\boldsymbol{v}\cdot\boldsymbol{n}=0\;\text{on}\;\partial\Omega\right\\},$
(6) $\displaystyle\boldsymbol{H}_{\boldsymbol{g}}(\mathrm{div};\Omega)$
$\displaystyle\coloneqq\left\\{\boldsymbol{v}\in\boldsymbol{H}(\mathrm{div};\Omega):(\boldsymbol{v}-\boldsymbol{g})\cdot\boldsymbol{n}=0\;\text{on}\;\partial\Omega\right\\},$
(7)
$\displaystyle\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)$
$\displaystyle\coloneqq\left\\{\boldsymbol{v}\in\boldsymbol{H}_{\boldsymbol{g}}(\mathrm{div};\Omega):\mathrm{div}(\boldsymbol{v})=0\;\text{a.e.~{}in}\;\Omega\right\\}.$

We note that $\boldsymbol{v}\cdot\boldsymbol{n}$ on $\partial\Omega$ is well-
defined for all $\boldsymbol{v}\in\boldsymbol{H}(\mathrm{div};\Omega)$ [8, Th.
3.12]. Moreover, $\boldsymbol{H}(\mathrm{div};\Omega)$ is a Hilbert space when
equipped with the inner product

(8)
$\displaystyle(\boldsymbol{u},\boldsymbol{v})_{\boldsymbol{H}(\mathrm{div};\Omega)}\coloneqq\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}+\mathrm{div}(\boldsymbol{u})\,\mathrm{div}(\boldsymbol{v})\,\mathrm{d}x.$

We extend the definition of $J$ in Eq. BP to functions $\boldsymbol{v}\not\in
H^{1}_{\boldsymbol{g}}(\Omega)^{d}$ by

(9) $\displaystyle J(\boldsymbol{v},\eta)=+\infty\;\;\text{for
all}\;\;\boldsymbol{v}\not\in H^{1}_{\boldsymbol{g}}(\Omega)^{d},\,\eta\in
C_{\gamma}.$

The following existence theorem is due to Borrvall and Petersson [13, Th.
3.1].

######  Theorem 2.1.

Suppose that $\Omega\subset\mathbb{R}^{d}$ is a Lipschitz domain, with
$d\in\\{2,3\\}$ and $\alpha$ satisfies properties (A1)–(A3). Then,
$J(\boldsymbol{v},\eta)$ is weak$\times$weak-* lower semicontinuous on
$H^{1}(\Omega)^{d}\times L^{\infty}(\Omega)$. Moreover, there exists a pair
$(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ that
minimizes $J$ (as defined in Eq. BP).

Although we are guaranteed the existence of the minimizer from Theorem 2.1,
the lack of convexity in the functional $J$ means that Eq. BP can support
multiple minimizers.

We now define the following forms:

(10) $\displaystyle a(\boldsymbol{u},\boldsymbol{v};\rho)$
$\displaystyle\coloneqq\int_{\Omega}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{v}+\nu\nabla\boldsymbol{u}:\nabla\boldsymbol{v}\;\mathrm{d}x,\quad
l(\boldsymbol{v})\coloneqq\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\;\mathrm{d}x,$
(11) $\displaystyle c(\rho,\eta;\boldsymbol{u})$
$\displaystyle\coloneqq\frac{1}{2}\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}\eta\;\mathrm{d}x,$
(12) $\displaystyle b(\boldsymbol{u},p)$
$\displaystyle\coloneqq-\int_{\Omega}p\,\mathrm{div}(\boldsymbol{u})\;\mathrm{d}x.$

######  Definition 2.2 (Strict minimizer).

Let $Z$ be a Banach space and suppose that $z_{0}\in Z$ is a local or global
minimizer of the functional $J:Z\to\mathbb{R}$. We say that $z_{0}$ is a
strict minimizer if there exists an open neighborhood $E\subset Z$ of $z_{0}$
such that $J(z_{0})<J(z)$ for all $z\neq z_{0}$, $z\in E$.

######  Definition 2.3 (Isolated minimizer).

Let $Z$ be a Banach space and suppose that $z_{0}\in Z$ is a local or global
minimizer of the functional $J:Z\to\mathbb{R}$. We say that $z_{0}$ is
isolated if there exists an open neighborhood $E\subset Z$ of $z_{0}$ such
that there are no other minimizers contained in $E$.

######  Remark 2.4.

If $z$ is an isolated minimizer, then it is also a strict minimizer.

The proof of the following proposition on first-order optimality conditions
for isolated minimizers of Eq. BP can be found in Papadopoulos and Süli [42,
Prop. 2].

######  Proposition 2.5.

Suppose that $\Omega\subset\mathbb{R}^{d}$ is a Lipschitz domain, with
$d\in\\{2,3\\}$ and $\alpha$ satisfies properties (A1)–(A4). Consider a
minimizer $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ whose
existence is guaranteed by Theorem 2.1. Then, there exists a unique Lagrange
multiplier $p\in L^{2}_{0}(\Omega)$ such that the following necessary first-
order optimality conditions hold:

(FOC1) $\displaystyle
a(\boldsymbol{u},\boldsymbol{v};\rho)+b(\boldsymbol{v},p)$
$\displaystyle=l(\boldsymbol{v})\;\;$ $\displaystyle\text{for
all}\;\;\boldsymbol{v}\in H^{1}_{0}(\Omega)^{d},$ (FOC2) $\displaystyle
b(\boldsymbol{u},q)$ $\displaystyle=0\;\;$ $\displaystyle\text{for
all}\;\;q\in L^{2}_{0}(\Omega),$ (FOC3) $\displaystyle
c(\rho,\eta-\rho;\boldsymbol{u})$ $\displaystyle\geq 0\;\;$
$\displaystyle\text{for all}\;\;\eta\in C_{\gamma}.$

##  3 Discretization

In this section we fix our choice of discretization. We denote the finite
element spaces for the velocity, pressure, and material distribution by
$\boldsymbol{X}_{h}\subset\boldsymbol{H}(\mathrm{div};\Omega)$, $M_{h}\subset
L^{2}_{0}(\Omega)$ and $C_{\gamma,h}\subset C_{\gamma}$, respectively. In the
following we introduce the notation that is required to define the discretized
optimization problem.

Let $(\mathcal{T}_{h})_{h\in(0,1]}$ denote a family of triangulations of the
domain $\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, characterized by the
mesh size $h\coloneqq\max_{K\in\mathcal{T}_{h}}h_{K}$, where $h_{K}$ is the
diameter of the element $K\in\mathcal{T}_{h}$. We assume that
$(\mathcal{T}_{h})_{h\in(0,1]}$ satisfies:

  1. (M1)

(Shape regularity). There exists constants $c_{1},c_{2}>0$ such that

$\displaystyle c_{1}h^{d}_{K}\leq|K|\leq c_{2}h^{d}_{K}\;\;\text{for
all}\;K\in\mathcal{T}_{h}.$

Moreover, we assume that $(\mathcal{T}_{h})_{h\in(0,1]}$ satisfies a submesh
condition as found in Buffa and Ortner [16, Assumption 2.1]. For a given $h$,
let the set $\mathcal{F}_{h}$ denote the set of all facets of the
triangulation $\mathcal{T}_{h}$ and $h_{F}$ represent the diameter of each
facet $F\in\mathcal{F}_{h}$. We split the set of facets into the union
$\mathcal{F}_{h}=\mathcal{F}^{i}_{h}\cup\mathcal{F}^{\partial}_{h}$ where
$\mathcal{F}^{i}_{h}$ is the subset of interior facets and
$\mathcal{F}^{\partial}_{h}$ collects all Dirichlet boundary facets
$F\subset\partial\Omega$. We note that we only consider Dirichlet boundary
conditions in this work, and, hence, all the facets on the boundary are
elements of the set $\mathcal{F}^{\partial}_{h}$. We assume the following:

  1. (M2)

(Contact regularity). There exists a constant $c_{1}>0$ such that

$\displaystyle c_{1}h^{d-1}_{K}\leq\mathcal{H}^{d-1}(F)\;\;\text{for
all}\;F\in\mathcal{F}_{h},\,K\in\mathcal{T}_{h}\;\;\text{such
that}\;\;F\subset\bar{K};$

  2. (M3)

(Boundary regularity). There exists a constant $c_{1}>0$ such that:

$\displaystyle c_{1}h\leq h_{F}\;\;\text{for
all}\;F\in\mathcal{F}^{\partial}_{h}.$

If $F\in\mathcal{F}^{i}_{h}$, then $F=\overline{\partial
K^{+}}\cap\overline{\partial K^{-}}$ for two elements
$K^{-},K^{+}\in\mathcal{T}_{h}$. We write $\boldsymbol{n}^{+}_{F}$ and
$\boldsymbol{n}^{-}_{F}$ to denote the outward normal unit vectors to the
boundaries $\partial K^{+}$ and $\partial K^{-}$, respectively. If
$F\in\mathcal{F}^{\partial}_{h}$, then $\boldsymbol{n}_{F}$ is the outer unit
normal vector $\boldsymbol{n}$. We denote the space of discontinuous finite
element functions with degree no higher than $k$ by

(13) $\displaystyle X_{\mathrm{DG}_{k}}\coloneqq\\{v\in
L^{1}(\Omega):v|_{K}\in\mathcal{P}_{k}\;\text{for
all}\;K\in\mathcal{T}_{h}\\},$

where $\mathcal{P}_{k}$ denotes the set of polynomials of order no higher than
$k$. Let $\boldsymbol{\phi}\in(X_{\mathrm{DG}_{k}})^{d}$ and
$\boldsymbol{\mathsf{\Phi}}\in(X_{\mathrm{DG}_{k}})^{d\times d}$ be any
piecewise vector or matrix-valued function, respectively with traces from
within the interior of $K^{\pm}$ denoted by $\boldsymbol{\phi}^{\pm}$ and
$\boldsymbol{\mathsf{\Phi}}^{\pm}$, respectively. We define the jump
$[\\![\cdot]\\!]_{F}$ and the average $\\{\\!\\!\\{\cdot\\}\\!\\!\\}_{F}$
operators across interior facets $F\in\mathcal{F}^{i}_{h}$ by

(14)
$\displaystyle[\\![\boldsymbol{\phi}]\\!]_{F}=\boldsymbol{\phi}^{+}\otimes\boldsymbol{n}_{F}^{+}+\boldsymbol{\phi}^{-}\otimes\boldsymbol{n}_{F}^{-}\quad\text{and}\quad\\{\\!\\!\\{\boldsymbol{\mathsf{\Phi}}\\}\\!\\!\\}_{F}=\frac{1}{2}\left(\boldsymbol{\mathsf{\Phi}}^{+}+\boldsymbol{\mathsf{\Phi}}^{-}\right).$

Here, $\otimes$ denotes the tensor product, i.e. for two vectors
$\boldsymbol{\mathbf{a}},\boldsymbol{\mathbf{b}}\in\mathbb{R}^{d}$,
$\boldsymbol{\mathbf{a}}\otimes\boldsymbol{\mathbf{b}}\in\mathbb{R}^{d\times
d}$,
$[\boldsymbol{\mathbf{a}}\otimes\boldsymbol{\mathbf{b}}]_{ij}=\boldsymbol{\mathbf{a}}_{i}\boldsymbol{\mathbf{b}}_{j}$.
If $F\in\mathcal{F}^{\partial}_{h}$, we set
$[\\![\boldsymbol{\phi}]\\!]_{F}=\boldsymbol{\phi}\otimes\boldsymbol{n}_{F}$
and
$\\{\\!\\!\\{\boldsymbol{\mathsf{\Phi}}\\}\\!\\!\\}_{F}=\boldsymbol{\mathsf{\Phi}}$.
For any $F\in\mathcal{F}^{\partial}_{h}$, we define
$\int_{F}|[\\![\boldsymbol{\phi}]\\!]_{F}|^{2}\,\mathrm{d}s=\int_{F}|\boldsymbol{\phi}|^{2}\,\mathrm{d}s$.

In general functions $v\in X_{\mathrm{DG}_{k}}$ do not live in $H^{1}(\Omega)$
due to the jumps across the facets of the elements. However, on each element,
$K\in\mathcal{T}_{h}$, $v|_{K}$ is a polynomial and, therefore, $v\in
H^{1}(K)$. We define the broken Sobolev space $H^{1}(\mathcal{T}_{h})$ as:

(15) $\displaystyle H^{1}(\mathcal{T}_{h})\coloneqq\\{v\in L^{1}(\Omega):v\in
H^{1}(K)\;\text{for all}\;K\in\mathcal{T}_{h}\\}.$

Moreover, for a function $\boldsymbol{v}\in H^{1}(\mathcal{T}_{h})^{d}$, we
define the broken $H^{1}$-seminorm and norms as:

(16) $\displaystyle|\boldsymbol{v}|^{2}_{H^{1}(\mathcal{T}_{h})}$
$\displaystyle\coloneqq\sum_{K\in\mathcal{T}_{h}}\|\nabla\boldsymbol{v}\|^{2}_{L^{2}(K)}+\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{v}]\\!]_{F}|^{2}\mathrm{d}s,$
(17) $\displaystyle\|\boldsymbol{v}\|^{2}_{H^{1}(\mathcal{T}_{h})}$
$\displaystyle\coloneqq\|\boldsymbol{v}\|^{2}_{L^{2}(\Omega)}+|\boldsymbol{v}|^{2}_{H^{1}(\mathcal{T}_{h})},$
(18)
$\displaystyle\|\boldsymbol{v}\|^{2}_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}$
$\displaystyle\coloneqq\|\boldsymbol{v}\|^{2}_{H^{1}(\mathcal{T}_{h})}+\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{v}-\boldsymbol{g}]\\!]_{F}|^{2}\mathrm{d}s.$

The two families of DG finite elements of interest for the velocity are the
Brezzi–Douglas–Marini (BDM) finite element [14, 15] and the Raviart–Thomas
(RT) finite element [46, 38]. The $k$-th order BDM finite element is defined
for $d=2$ in [15, Sec. 2] and for $d=3$ in [14, Sec. 2]. Similarly the $k$-th
order RT finite element is defined in [46, Sec. 3] and [38, Sec. 2] for $d=2$
and $d=3$, respectively. The finite element spaces induced by the $k$-th order
BDM and RT finite elements are denoted by $\boldsymbol{X}_{\mathrm{BDM}_{k}}$
and $\boldsymbol{X}_{\mathrm{RT}_{k}}$, respectively. We note that
$\boldsymbol{X}_{\mathrm{RT}_{k}}\subset\boldsymbol{X}_{\mathrm{BDM}_{k}}\subset\boldsymbol{Z}_{h}$
where, for a given $k\geq 1$,

(19) $\displaystyle\boldsymbol{Z}_{h}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in(X_{\mathrm{DG}_{k}})^{d}:\,\mathrm{div}(\boldsymbol{v})\in
X_{\mathrm{DG}_{k-1}}\cap L^{2}(\Omega)\\}.$

We note that $\boldsymbol{Z}_{h}\subset\boldsymbol{H}(\mathrm{div};\Omega)$
and $\boldsymbol{Z}_{h}\subset H^{1}(\mathcal{T}_{h})^{d}$.

We define the following subspaces of $\boldsymbol{X}_{h}$:

(20) $\displaystyle\boldsymbol{X}_{h,0}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in\boldsymbol{X}_{h}:\boldsymbol{v}\cdot\boldsymbol{n}=0\;\text{on}\;\partial\Omega\\},$
(21) $\displaystyle\boldsymbol{X}_{h,\boldsymbol{g}}$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in\boldsymbol{X}_{h}:(\boldsymbol{v}-\boldsymbol{g})\cdot\boldsymbol{n}=0\;\text{on}\;\partial\Omega\\}.$

In general the boundary data $\boldsymbol{g}$ cannot be represented in the
finite element space. Hence we instead approximate the boundary data with a
finite element function $\boldsymbol{g}_{h}$ (which can be represented) and
assume that

  1. (F1)

$h^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}\to 0$ as
$h\to 0$.

We also assume that:

  1. (F2)

$\boldsymbol{X}_{h,0}$ and $M_{h}$ satisfy the following inf-sup condition for
some $c_{b}>0$, independent of $h$,

(22) $\displaystyle c_{b}\leq\inf_{q_{h}\in
M_{h}\backslash\\{0\\}}\sup_{\boldsymbol{v}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{b(\boldsymbol{v}_{h},q_{h})}{\|\boldsymbol{v}_{h}\|_{H^{1}(\mathcal{T}_{h})}\|q_{h}\|_{L^{2}(\Omega)}}.$

  2. (F3)

The finite element spaces are dense in their respective function spaces, i.e.,
for any $(\boldsymbol{v},\eta,q)\in H^{1}(\mathcal{T}_{h})^{d}\times
C_{\gamma}\times L^{2}_{0}(\Omega)$,

$\displaystyle\lim_{h\to
0}\inf_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h}}\|\boldsymbol{v}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}$
$\displaystyle=\lim_{h\to 0}\inf_{\zeta_{h}\in
C_{\gamma,h}}\|\eta-\zeta_{h}\|_{L^{2}(\Omega)}$
$\displaystyle\indent=\lim_{h\to 0}\inf_{r_{h}\in
M_{h}}\|q-r_{h}\|_{L^{2}(\Omega)}=0.$

######  Remark 3.1.

The inf-sup (F2) and density (F3) conditions are satisfied by either
$\boldsymbol{X}_{h}=\boldsymbol{X}_{\mathrm{BDM}_{k}}$ or
$\boldsymbol{X}_{h}=\boldsymbol{X}_{\mathrm{RT}_{k}}$ with the choice of the
pressure finite element space $M_{h}=X_{\mathrm{DG}_{k-1}}$ [17].

######  Remark 3.2.

Although $C_{\gamma}$ is not separable, we will only require the density of
$C_{\gamma,h}$ in $C_{\gamma}$ with respect to the $L^{2}$-norm.

We now define the discrete Borrvall–Petersson power dissipation functional for
a DG finite element discretization. Consider the functions
$\boldsymbol{u}_{h}\in\boldsymbol{X}_{h}\subset\boldsymbol{Z}_{h}$ and
$\rho_{h}\in C_{\gamma,h}$. We note that $J(\boldsymbol{u}_{h},\rho_{h})$ is
ill-defined as, in general, $\nabla\boldsymbol{u}_{h}\not\in
L^{2}(\Omega)^{d}$. Hence, the term
$\int_{\Omega}|\nabla\boldsymbol{u}_{h}|^{2}\,\mathrm{d}x$, as found in
$J(\boldsymbol{u}_{h},\rho_{h})$, might not be finite. Given a penalization
parameter $\sigma>0$, we define the discrete analogue $J_{h}$ as

(23)
$\displaystyle\begin{split}J_{h}(\boldsymbol{u}_{h},\rho_{h})&\coloneqq\frac{1}{2}\int_{\Omega}\left(\alpha(\rho_{h})|\boldsymbol{u}_{h}|^{2}-2\boldsymbol{f}\cdot\boldsymbol{u}_{h}\right)\text{d}x+\frac{\nu}{2}\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla\boldsymbol{u}_{h}|^{2}\mathrm{d}x\\\
&+\frac{\nu}{2}\sum_{F\in\mathcal{F}^{i}_{h}}\sigma
h_{F}^{-1}\int_{F}|[\\![\boldsymbol{u}_{h}]\\!]_{F}|^{2}\mathrm{d}s-\nu\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}\\{\\!\\!\\{\nabla\boldsymbol{u}_{h}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{u}_{h}]\\!]_{F}\,\mathrm{d}s\\\
&+\frac{\nu}{2}\sum_{F\in\mathcal{F}^{\partial}_{h}}\sigma
h_{F}^{-1}\int_{F}|[\\![\boldsymbol{u}_{h}-\boldsymbol{g}_{h}]\\!]_{F}|^{2}\mathrm{d}s\\\
&-\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}\\{\\!\\!\\{\nabla\boldsymbol{u}_{h}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{u}_{h}-\boldsymbol{g}_{h}]\\!]_{F}\,\mathrm{d}s.\end{split}$

######  Remark 3.3.

This particular choice $J_{h}$ as the discrete analogue of $J$ is motivated by
an interior penalty approach for DG formulations. In Proposition 4.11 we prove
that velocity minimizers of $J_{h}$ satisfy a fluid momentum equation
featuring terms that arise in the interior penalty DG discretization of the
Stokes equations [18, 30].

######  Remark 3.4.

For any $\boldsymbol{v}_{h}\in\boldsymbol{Z}_{h}$ the terms
$\\{\\!\\!\\{\nabla\boldsymbol{v}_{h}\\}\\!\\!\\}_{F}$ and
$[\\![\boldsymbol{v}_{h}]\\!]_{F}$ as they appear in $J_{h}$ are well-defined
[9, Sec. 3.1].

######  Proposition 3.5 (Consistency of $J_{h}$).

Consider any $(\boldsymbol{v},\eta)\in
H^{1}_{\boldsymbol{g}}(\Omega)^{d}\times C_{\gamma}$ such that
$\boldsymbol{v}\in H^{r}(\Omega)^{d}$ for some $r>3/2$. Then, $J_{h}$, $h>0$
is consistent, i.e.

(24) $\displaystyle J_{h}(\boldsymbol{v},\eta)=J(\boldsymbol{v},\eta).$

######  Proof 3.6.

Since $\boldsymbol{v}\in H^{r}(\Omega)^{d}$, for $r>3/2$, we note that
$J_{h}(\boldsymbol{v},\eta)$ is well-defined and there can be no jumps in
$\boldsymbol{v}$ across elements. Hence, for all $F\in\mathcal{F}^{i}_{h}$,
integrals involving $[\\![\boldsymbol{v}]\\!]_{F}$ are equal to zero. Moreover
$\boldsymbol{v}=\boldsymbol{g}$ on $\partial\Omega$ and, therefore,
$\int_{F}|[\\![\boldsymbol{v}-\boldsymbol{g}]\\!]_{F}|^{2}\,\mathrm{d}s=0$ for
all $F\in\mathcal{F}^{\partial}_{h}$. Hence,

(25) $\displaystyle
J_{h}(\boldsymbol{v},\eta)=\frac{1}{2}\int_{\Omega}\left(\alpha(\eta)|\boldsymbol{v}|^{2}-2\boldsymbol{f}\cdot\boldsymbol{v}\right)\text{d}x+\frac{\nu}{2}\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla\boldsymbol{v}|^{2}\mathrm{d}x=J(\boldsymbol{v},\eta).$

For a sufficiently large penalization parameter $\sigma>0$, we define the
broken form $a_{h}(\boldsymbol{u},\boldsymbol{v};\rho)$ by

(26)
$\displaystyle\begin{split}a_{h}(\boldsymbol{u},\boldsymbol{v};\rho)&\coloneqq\sum_{K\in\mathcal{T}_{h}}\int_{K}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{v}+\nabla\boldsymbol{u}:\nabla\boldsymbol{v}\;\mathrm{d}x+\nu\sum_{F\in\mathcal{F}_{h}}\sigma
h_{F}^{-1}\int_{F}[\\![\boldsymbol{u}]\\!]_{F}:[\\![\boldsymbol{v}]\\!]_{F}\mathrm{d}s\\\
&-\nu\sum_{F\in\mathcal{F}_{h}}\int_{F}\\{\\!\\!\\{\nabla\boldsymbol{u}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{v}]\\!]_{F}\mathrm{d}s-\nu\sum_{F\in\mathcal{F}_{h}}\int_{F}[\\![\boldsymbol{u}]\\!]_{F}:\\{\\!\\!\\{\nabla\boldsymbol{v}\\}\\!\\!\\}_{F}\mathrm{d}s,\end{split}$

and the linear functional $l_{h}(\boldsymbol{v};\boldsymbol{g})$ as

(27)
$\displaystyle\begin{split}l_{h}(\boldsymbol{v};\boldsymbol{g})&\coloneqq\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\;\mathrm{d}x\\\
&\indent+\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\sigma
h_{F}^{-1}\int_{F}[\\![\boldsymbol{g}]\\!]_{F}:[\\![\boldsymbol{v}]\\!]_{F}\;\mathrm{d}s-\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}[\\![\boldsymbol{g}]\\!]_{F}:\\{\\!\\!\\{\nabla\boldsymbol{v}\\}\\!\\!\\}_{F}\,\mathrm{d}s.\end{split}$

######  Proposition 3.7 (Consistency of $a_{h}$).

Suppose that $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ is a
minimizer of Eq. BP and let $p\in L^{2}_{0}(\Omega)$ denote the Lagrange
multiplier such that $(\boldsymbol{u},\rho,p)$ satisfy Eq. FOC1–Eq. FOC3.
Moreover, assume that $\boldsymbol{u}\in H^{r}(\Omega)^{d}$ for some $r>3/2$.
Then, for all $\boldsymbol{v}_{h}\in
H^{1}(\mathcal{T}_{h})^{d}\cap\boldsymbol{H}_{0}(\mathrm{div};\Omega)$, we
have that

(28) $\displaystyle
a_{h}(\boldsymbol{u},\boldsymbol{v}_{h};\rho)+b(\boldsymbol{v}_{h},p)=l_{h}(\boldsymbol{v}_{h};\boldsymbol{g}).$

######  Proof 3.8.

By Eq. FOC1, we have that, for all $\boldsymbol{\phi}\in
H^{1}_{0}(\Omega)^{d}$, and therefore, for all $\boldsymbol{\phi}\in
C_{c}^{\infty}(\Omega)^{d}$,

(29) $\displaystyle
a(\boldsymbol{u},\boldsymbol{\phi};\rho)+b(\boldsymbol{\phi},p)=l(\boldsymbol{\phi}),$

where $C_{c}^{\infty}(\Omega)$ denotes set of smooth and compactly supported
functions in $\Omega$. By an integration by parts, we see that

(30)
$\displaystyle\int_{\Omega}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{\phi}-\nu\Delta\boldsymbol{u}\cdot\boldsymbol{\phi}+\nabla
p\cdot\boldsymbol{\phi}\,\mathrm{d}x=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{\phi}\,\mathrm{d}x.$

We note that
$\int_{\Omega}\boldsymbol{\psi}\cdot\boldsymbol{\phi}\,\mathrm{d}x$ is well-
defined for any $\boldsymbol{\psi}\in(C_{c}^{\infty}(\Omega)^{d})^{*}$ (the
space of distributions). Hence, Eq. 30 is well-defined since
$\alpha(\rho)\boldsymbol{u}$, $\Delta\boldsymbol{u}$, $\nabla p$,
$\boldsymbol{f}\in(C_{c}^{\infty}(\Omega)^{d})^{*}$.

As the set of smooth functions is dense in $L^{2}(\Omega)$, we can test Eq. 30
against any $\boldsymbol{v}_{h}\in
H^{1}(\mathcal{T}_{h})^{d}\cap\boldsymbol{H}_{0}(\mathrm{div};\Omega)\subset
L^{2}(\Omega)^{d}$. Thus, by performing a second integration by parts, we have
that

(31)
$\displaystyle\begin{split}&\sum_{K\in\mathcal{T}_{h}}\int_{K}\alpha(\rho)\boldsymbol{u}\cdot\boldsymbol{v}_{h}+\nu\nabla\boldsymbol{u}:\nabla\boldsymbol{v}_{h}-p\,\mathrm{div}(\boldsymbol{v}_{h})\,\mathrm{d}x\\\
&\indent-\nu\sum_{F\in\mathcal{F}_{h}}\int_{F}\\{\\!\\!\\{\nabla\boldsymbol{u}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{v}_{h}]\\!]_{F}\,\mathrm{d}s=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}_{h}\,\mathrm{d}x.\end{split}$

The element-wise surface integrals arising by the integration by parts of the
$\nabla p$ term drop out due to the continuity of
$\boldsymbol{v}_{h}\cdot\boldsymbol{n}$ across elements for all
$\boldsymbol{v}_{h}\in\boldsymbol{H}_{0}(\mathrm{div};\Omega)$. Similarly the
boundary surface integrals drop out since
$\boldsymbol{v}_{h}\cdot\boldsymbol{n}=0$ on $\partial\Omega$. As
$\boldsymbol{u}\in H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}$, for all
$F\in\mathcal{F}^{i}_{h}$, we have that $[\\![\boldsymbol{u}]\\!]_{F}=0$ and
for all $F\in\mathcal{F}^{\partial}_{h}$,
$[\\![\boldsymbol{u}]\\!]_{F}=[\\![\boldsymbol{g}]\\!]_{F}$. As
$\boldsymbol{u}\in H^{r}(\Omega)^{d}$, for some $r>3/2$, the traces of
$\nabla\boldsymbol{u}$ on $F\in\mathcal{F}_{h}$ are well-defined. We conclude
that Eq. 28 holds.

######  Proposition 3.9 (Coercivity and boundedness of $a_{h}$).

There exists a $\sigma_{0}>0$, such that for all $\sigma\geq\sigma_{0}$,
$\boldsymbol{w}_{h},\boldsymbol{u}_{h}\in\boldsymbol{X}_{h}$ and $\eta\in
C_{\gamma}$, there exists constants $c_{a},C_{a}>0$ such that

(32) $\displaystyle c_{a}\|\boldsymbol{w}_{h}\|^{2}_{H^{1}(\mathcal{T}_{h})}$
$\displaystyle\leq a_{h}(\boldsymbol{w}_{h},\boldsymbol{w}_{h};\eta),$ (33)
$\displaystyle a_{h}(\boldsymbol{w}_{h},\boldsymbol{u}_{h};\eta)$
$\displaystyle\leq
C_{a}\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}\|\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}.$

######  Proof 3.10.

We note that, by assumption (A1), $0\leq\alpha(\eta)\leq\bar{\alpha}$ for all
$\eta\in C_{\gamma}$. Hence, the result follows from classical coercivity and
boundedness results for DG discretizations for interior penalty methods [9,
Sec. 4.1–4.2].

######  Definition 3.11.

We define the spaces $\boldsymbol{U}_{h,\boldsymbol{g}_{h}}$ and
$\boldsymbol{U}_{h,0}$ as:

(34) $\displaystyle\boldsymbol{U}_{h,\boldsymbol{g}_{h}}$
$\displaystyle\coloneqq\\{\boldsymbol{u}\in\boldsymbol{X}_{h,\boldsymbol{g}_{h}}:b(\boldsymbol{u}_{h},q_{h})=0\;\text{for
all}\;q_{h}\in X_{\mathrm{DG}_{k-1}}\\}$ (35)
$\displaystyle\boldsymbol{U}_{h,0}$
$\displaystyle\coloneqq\\{\boldsymbol{u}\in\boldsymbol{X}_{h,0}:b(\boldsymbol{u}_{h},q_{h})=0\;\text{for
all}\;q_{h}\in X_{\mathrm{DG}_{k-1}}\\}.$

In the following lemma we provide the proof that functions
$\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,0}$ and
$\boldsymbol{v}_{h}\in\boldsymbol{U}_{\boldsymbol{g}_{h},0}$ are pointwise
divergence-free.

######  Lemma 3.12 (Pointwise divergence-free).

Suppose that $\boldsymbol{X}_{h}\subset\boldsymbol{Z}_{h}$. Consider a
function $\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,0}$ or
$\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}$. Then,
$\mathrm{div}(\boldsymbol{v}_{h})=0$ a.e. in $\Omega$.

######  Proof 3.13.

Since
$\boldsymbol{U}_{h,0},\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\subset\boldsymbol{X}_{h}\subset\boldsymbol{Z}_{h}$
then, by definition, $\mathrm{div}(\boldsymbol{v}_{h})\in
X_{\mathrm{DG}_{k-1}}$. Hence, there exists a $q_{h}\in X_{\mathrm{DG}_{k-1}}$
such that $q_{h}=\mathrm{div}(\boldsymbol{v}_{h})$. Therefore,

(36) $\displaystyle
b(\boldsymbol{v}_{h},q_{h})=-\|\mathrm{div}(\boldsymbol{v}_{h})\|^{2}_{L^{2}(\Omega)}=0,$

which implies that $\mathrm{div}(\boldsymbol{v}_{h})=0$ a.e. in $\Omega$.

To prove the convergence of a DG finite element method, we require the
existence of sequences in $\boldsymbol{U}_{h,\boldsymbol{g}_{h}}$ that
converge strongly to $\boldsymbol{u}$.

######  Lemma 3.14 (Strongly converging sequences).

Suppose that (F1)–(F3) hold and $\boldsymbol{X}_{h}\subset\boldsymbol{Z}_{h}$.
Consider any minimizer $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ of Eq. BP.
Then, there exists a sequence of functions
$(\boldsymbol{\tilde{u}}_{h},\tilde{p}_{h})\in\boldsymbol{U}_{\boldsymbol{g}_{h},h}\times
M_{h}$ such that
$\|\boldsymbol{u}-\boldsymbol{\tilde{u}}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$ and $\|p-\tilde{p}_{h}\|_{L^{2}(\Omega)}\to 0$.

######  Proof 3.15.

For sufficiently large $\sigma>0$ and fixed $\rho\in C_{\gamma}$, consider the
problem, find
$(\boldsymbol{\tilde{u}}_{h},\tilde{p}_{h})\in\boldsymbol{U}_{\boldsymbol{g}_{h},h}\times
M_{h}$ that satisfies

(37) $\displaystyle
a_{h}(\boldsymbol{\tilde{u}}_{h},\boldsymbol{v}_{h};\rho)+b(\boldsymbol{v}_{h},\tilde{p}_{h})$
$\displaystyle=l_{h}(\boldsymbol{v}_{h})$ $\displaystyle\text{for
all}\;\boldsymbol{v}_{h}\in\boldsymbol{X}_{h,0},$ (38) $\displaystyle
b(\boldsymbol{\tilde{u}}_{h},q_{h})$ $\displaystyle=0$ $\displaystyle\text{for
all}\;q_{h}\in M_{h}.$

Then, under assumptions (F1)–(F3), by standard results for
$\boldsymbol{H}(\mathrm{div};\Omega)$ finite element discretizations of the
Stokes and Stokes–Brinkman equations with an interior penalty [18, 34, 35],
the pair $(\boldsymbol{\tilde{u}}_{h},\tilde{p}_{h})$ exists, is unique, and
$\|\boldsymbol{u}-\boldsymbol{\tilde{u}}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$, $\|p-\tilde{p}_{h}\|_{L^{2}(\Omega)}\to 0$ as $h\to 0$.

##  4 Convergence of a DG finite element method

In their original paper Borrvall and Petersson [13, Sec. 3.3] considered a
piecewise constant finite element approximation of the material distribution
coupled with an inf-sup stable quadrilateral finite element approximation of
the velocity and the pressure. In particular their velocity finite element
space was $H^{1}$-conforming. They showed that if the domain $\Omega$ is a
polygonal domain in two dimensions or a polyhedral Lipschitz domain in three
dimensions, such approximations of the velocity and material distribution
$(\boldsymbol{u}_{h},\rho_{h})$ that minimize $J(\boldsymbol{u}_{h},\rho_{h})$
converge to an unspecified solution $(\boldsymbol{u},\rho)$ of Eq. BP in the
following sense [13, Th. 3.2]:

$\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\rightharpoonup\boldsymbol{u}\;\;\text{weakly
in}\;\;H^{1}(\Omega)^{d},$ $\displaystyle\rho_{h}$
$\displaystyle\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho\;\;\text{weakly-*
in}\;\;L^{\infty}(\Omega),$ $\displaystyle\rho_{h}$
$\displaystyle\to\rho\;\;\text{strongly
in}\;\;L^{s}(\Omega_{b}),\;\;s\in[1,\infty),$

where $\Omega_{b}$ is any measurable subset of $\Omega$ where $\rho$ is equal
to zero or one a.e. Their analysis suggests that a finite element method is a
suitable discretization, but it left a number of open problems:

  1. (P1)

It is not clear which minimizer the sequence is converging to as the
nonconvexity of the problem provides multiple candidates for the limits;

  2. (P2)

The convergence is weak-* in the material distribution in regions where
$\\{0<\rho<1\\}\subset\Omega$ which permits the presence of checkerboard
patterns as $h\to 0$;

  3. (P3)

There are no convergence results for the finite element approximation of the
pressure, $p$.

In general (P1) means that their result does not imply that there exists a
sequence of finite element solutions that converges to the global minimizer.

Recently Borrvall and Petersson’s result was extended and refined by
Papadopoulos and Süli [42, Th. 4]. They considered conforming discretizations
of the material distribution and conforming inf-sup stable finite element
discretizations for the velocity-pressure pair. Once again, the velocity
finite element space $\boldsymbol{X}_{h,\boldsymbol{g}_{h}}$ was assumed to be
$H^{1}$-conforming. They showed that, for any isolated minimizer
$(\boldsymbol{u},\rho)$ of Eq. BP, there exists a sequence of solutions
$(\boldsymbol{u}_{h},\rho_{h},p_{h})\in\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\times
C_{\gamma,h}\times M_{h}$ satisfying the discretized first-order optimality
condititions, for all
$(\boldsymbol{v}_{h},\eta_{h},q_{h})\in\boldsymbol{X}_{h,0}\times
C_{\gamma,h}\times M_{h}$:

(39) $\displaystyle
a(\boldsymbol{u}_{h},\boldsymbol{v}_{h};\rho_{h})+b(\boldsymbol{v}_{h},p_{h})$
$\displaystyle=l(\boldsymbol{v}_{h}),$ (40) $\displaystyle
b(\boldsymbol{u}_{h},q_{h})$ $\displaystyle=0,$ (41) $\displaystyle
c(\rho_{h},\eta_{h}-\rho_{h};\boldsymbol{u}_{h})$ $\displaystyle\geq 0,$

such that

$\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\to\boldsymbol{u}\;\;\text{strongly in}\;\;H^{1}(\Omega)^{d},$
$\displaystyle\rho_{h}$ $\displaystyle\to\rho\;\;\text{strongly
in}\;\;L^{s}(\Omega),\;s\in[1,\infty),$ $\displaystyle p_{h}$
$\displaystyle\to p\;\;\text{strongly in}\;\;L^{2}(\Omega).$

Their analysis resolved the open issues (P1)–(P3). The assumption that
$\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\subset H^{1}(\Omega)^{d}$ was crucial
for the compactness results utilized in their proof.

Our goal in this section is to prove a similar result for when
$\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\not\subset H^{1}(\Omega)^{d}$ but
$\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\subset\boldsymbol{Z}_{h}\subset\boldsymbol{H}(\mathrm{div};\Omega)$.
The following theorem is the main result of this work.

######  Theorem 4.1 (Convergence of the DG finite element method).

Let $\Omega\subset\mathbb{R}^{d}$ be a polygonal domain in two dimensions or a
polyhedral Lipschitz domain in three dimensions. Suppose that the inverse
permeability $\alpha$ satisfies (A1)–(A4) and there exists an isolated
minimizer $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ of Eq. BP
that has the additional regularity $\boldsymbol{u}\in H^{r}(\Omega)^{d}$ for
some $r>3/2$. Moreover, assume that, for $\theta>0$, $U_{\theta}$ is the
subset of $\Omega$ where $|\boldsymbol{u}|^{2}\geq\theta$ a.e. in $U_{\theta}$
and suppose that there exists a $\theta^{\prime}>0$ such that $U_{\theta}$ is
closed and has non-empty interior for all $\theta\leq\theta^{\prime}$. Let $p$
denote the unique Lagrange multiplier associated with $(\boldsymbol{u},\rho)$
such that $(\boldsymbol{u},\rho,p)$ satisfy the first-order optimality
conditions Eq. FOC1–Eq. FOC3.

Consider the finite element spaces
$\boldsymbol{X}_{h}\subset\boldsymbol{Z}_{h}$, $C_{\gamma,h}\subset
C_{\gamma}$, and $M_{h}\subset L^{2}_{0}(\Omega)$ and suppose that the
assumptions (F1)–(F3) hold.

Then, there exists an $\bar{h}>0$ such that, for $h\leq\bar{h}$, $h\to 0$,
there is a family of solutions
$(\boldsymbol{u}_{h},\rho_{h},p_{h})\in\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\times
C_{\gamma,{h}}\times M_{h}$ to the following discretized first-order
optimality conditions

(FOC1-h) $\displaystyle
a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h};\rho_{h})+b(\boldsymbol{v}_{h},p_{h})$
$\displaystyle=l_{h}(\boldsymbol{v}_{h};\boldsymbol{g}_{h})$
$\displaystyle\text{for all}\;\boldsymbol{v}_{h}\in\boldsymbol{X}_{h,0},$
(FOC2-h) $\displaystyle b(\boldsymbol{u}_{h},q_{h})$ $\displaystyle=0$
$\displaystyle\text{for all}\;q_{h}\in M_{h},$ (FOC3-h) $\displaystyle
c(\rho_{h},\eta_{h}-\rho_{h};\boldsymbol{u}_{h})$ $\displaystyle\geq 0$
$\displaystyle\text{for all}\;\eta_{h}\in C_{\gamma,{h}},$

such that,
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$, $\rho_{h}\to\rho$ strongly in $L^{s}(\Omega)$, $s\in[1,\infty)$, and
$p_{h}\to p$ strongly in $L^{2}(\Omega)$ as $h\to 0$.

######  Remark 4.2.

Convergence of $\boldsymbol{u}_{h}$ in the norm
$\|\cdot\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}$ ensures that the
Dirichlet boundary condition is correctly satisfied in the limit. This cannot
be immediately deduced if $\boldsymbol{u}_{h}$ only strongly converges in the
norm $\|\cdot\|_{H^{1}(\mathcal{T}_{h})}$.

We first introduce some auxiliary propositions to facilitate the proof of
Theorem 4.1. In all the following propositions and corollaries, we assume that
the conditions in Theorem 4.1 hold and we fix an isolated minimizer
$(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ of Eq. BP.
We also fix an $r>0$ such that $(\boldsymbol{u},\rho)$ is the unique local
minimizer of Eq. BP in $B_{r,\boldsymbol{H}(\mathrm{div};\Omega)\times
L^{2}(\Omega)}(\boldsymbol{u},\rho)\cap(\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)\times
C_{\gamma})$, where

(42)
$\displaystyle\begin{split}&B_{r,\boldsymbol{H}(\mathrm{div};\Omega)\times
L^{2}(\Omega)}(\boldsymbol{u},\rho)\\\
&\coloneqq\\{\boldsymbol{v}\in\boldsymbol{H}(\mathrm{div};\Omega),\;\eta\in
C_{\gamma}:\|\boldsymbol{u}-\boldsymbol{v}\|_{\boldsymbol{H}(\mathrm{div};\Omega)}+\|\rho-\eta\|_{L^{2}(\Omega)}\leq
r\\}.\end{split}$

Such an $r$ is guaranteed to exist by the definition of an isolated minimizer
and the extension of $J$ in Eq. 9 to functions
$\boldsymbol{v}\in\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)$
such that $\boldsymbol{v}\not\in H^{1}_{\boldsymbol{g}}(\Omega)^{d}$. We also
define $B_{r,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$ and
$B_{r,L^{2}(\Omega)}(\rho)$ by

(43) $\displaystyle B_{r,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$
$\displaystyle\coloneqq\\{\boldsymbol{v}\in\boldsymbol{H}(\mathrm{div};\Omega):\|\boldsymbol{u}-\boldsymbol{v}\|_{\boldsymbol{H}(\mathrm{div};\Omega)}\leq
r\\},$ (44) $\displaystyle B_{r,L^{2}(\Omega)}(\rho)$
$\displaystyle\coloneqq\\{\eta\in C_{\gamma}:\|\rho-\eta\|_{L^{2}(\Omega)}\leq
r\\}.$

We note that

$\displaystyle(\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$ $\displaystyle\indent\subset
B_{r,\boldsymbol{H}(\mathrm{div};\Omega)\times
L^{2}(\Omega)}(\boldsymbol{u},\rho)\cap(\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)\times
C_{\gamma})$

and hence $(\boldsymbol{u},\rho)$ is also the unique minimizer in
$(\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$.

######  Remark 4.3.

The extension of $J$ in Eq. 9 to functions
$\boldsymbol{v}\in\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)$
such that $\boldsymbol{v}\not\in H^{1}_{\boldsymbol{g}}(\Omega)^{d}$ means
that functions $\boldsymbol{v}\not\in H^{1}_{\boldsymbol{g}}(\Omega)^{d}$
cannot be minimizers.

######  Remark 4.4.

We make the assumption that $\rho\in C_{\gamma}$ is isolated with respect to
the $L^{2}$-norm (as opposed to the $L^{\infty}$-norm). This is a stronger
isolation assumption as discussed in [42, Rem. 7]. However, it is equivalent
to being isolated with respect to any $L^{s}$-norm for $s\in[1,\infty)$ as
discussed in [42, Rem. 8]. As far as we are aware, the $L^{2}$-isolation
assumption is valid for all practical problems found in the literature, in
particular it holds for the example found in Section 5.

######  Remark 4.5.

We are required to make the assumption that $\boldsymbol{u}$ is isolated with
respect to the $\boldsymbol{H}(\mathrm{div};\Omega)$-norm (as opposed to the
$H^{1}$-norm), i.e. for all
$\boldsymbol{v}\in\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)$
such that $\boldsymbol{v}\in
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$, then
$\boldsymbol{v}$ cannot be part of a minimizing pair. This is necessary to
construct the discretized problem Eq. BP-h below. This is a stronger isolation
assumption than isolation with respect to the $H^{1}$-norm. However, in
examples found in the literature and, in particular, the example we consider
in Section 5, the velocity minimizers are always isolated with respect to the
$\boldsymbol{H}(\mathrm{div};\Omega)$-norm.

######  Proposition 4.6 (Weak convergence of ($\boldsymbol{u}_{h},\rho_{h})$
in $\boldsymbol{H}(\mathrm{div};\Omega)\times L^{2}(\Omega)$).

For a given $h>0$, consider the finite-dimensional optimization problem: find
$(\boldsymbol{u}_{h},\rho_{h})\in(\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$ that minimizes

(BP-h) $\displaystyle J_{h}(\boldsymbol{u}_{h},\rho_{h}).$

Then, a global minimizer $(\boldsymbol{u}_{h},\rho_{h})$ of Eq. BP-h exists
and there exist subsequences (up to relabeling) such that as $h\to 0$:

(45) $\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\rightharpoonup\boldsymbol{u}\;\text{weakly
in}\;\boldsymbol{H}(\mathrm{div};\Omega),$ (46)
$\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\to\boldsymbol{u}\;\;\text{strongly in}\;\;L^{q}(\Omega)^{d},$
(47) $\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\to\boldsymbol{u}\;\;\text{strongly
in}\;\;L^{r}(\partial\Omega)^{d},$ (48) $\displaystyle\rho_{h}$
$\displaystyle\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho\;\text{weakly-*
in}\;L^{\infty}(\Omega),$ (49) $\displaystyle\rho_{h}$
$\displaystyle\rightharpoonup\rho\;\text{weakly
in}\;L^{s}(\Omega),\;s\in[1,\infty),$

where $1\leq q,r<\infty$ in two dimensions and $1\leq q<6$, $1\leq r<4$ in
three dimensions.

######  Proof 4.7.

The functional $J_{h}$ is continuous and

(50) $\displaystyle(\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$

is a finite-dimensional, closed and bounded set and, for sufficiently small
$h$, non-empty, therefore, sequentially compact by the Heine–Borel theorem
[28, Th. 11.18]. Hence, $J_{h}$ attains its infimum in
$(\boldsymbol{U}_{\boldsymbol{g}_{h},h}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$ and, therefore, a global minimizer
$(\boldsymbol{u}_{h},\rho_{h})$ exists.

By a corollary of Kakutani’s Theorem [29, Th. A.65], if a Banach space is
reflexive then every norm-closed, bounded and convex subset of the Banach
space is weakly compact and thus, by the Eberlein–Šmulian theorem [29, Th.
A.62], sequentially weakly compact. It can be checked that
$\boldsymbol{H}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$ and
$C_{\gamma}\cap B_{r/2,L^{2}(\Omega)}(\rho)$ are norm-closed, bounded and
convex subsets of the reflexive Banach spaces
$\boldsymbol{H}(\mathrm{div};\Omega)$ and $L^{2}(\Omega)$, respectively.
Therefore, $\boldsymbol{H}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$ is weakly
sequentially compact in $\boldsymbol{H}(\mathrm{div};\Omega)$ and
$C_{\gamma}\cap B_{r/2,L^{2}(\Omega)}(\rho)$ is weakly sequentially compact in
$L^{2}(\Omega)$.

Hence we extract subsequences (up to relabeling), $(\boldsymbol{u}_{h})$ and
$(\rho_{h})$ of the sequence generated by the global minimizers of Eq. BP-h
such that

(51) $\displaystyle\boldsymbol{u}_{h}$
$\displaystyle\rightharpoonup\hat{\boldsymbol{u}}\in\boldsymbol{H}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})\;\text{weakly
in}\;\boldsymbol{H}(\mathrm{div};\Omega),$ (52) $\displaystyle\rho_{h}$
$\displaystyle\rightharpoonup\hat{\rho}\in C_{\gamma}\cap
B_{r/2,L^{2}(\Omega)}(\rho)\;\text{weakly in}\;L^{2}(\Omega).$

By assumption (F3), there exists a sequence of finite element functions
$\tilde{\rho}_{h}\in C_{\gamma,h}$ that strongly converges to $\rho$ in
$L^{2}(\Omega)$. Moreover, Lemma 3.14 implies the existence of a sequence
$(\boldsymbol{\tilde{u}}_{h})\in\boldsymbol{U}_{\boldsymbol{g}_{h},h}$ that
satisfies
$\|\boldsymbol{u}-\boldsymbol{\tilde{u}}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$.

We now wish to identify the limits $\hat{\boldsymbol{u}}$ and $\hat{\rho}$.
Consider the following bound:

(53)
$\displaystyle\begin{split}&2|J_{h}(\boldsymbol{\tilde{u}}_{h},\tilde{\rho}_{h})-J(\boldsymbol{u},\rho)|\\\
&\indent\leq\int_{\Omega}|(\alpha(\rho)-\alpha(\tilde{\rho}_{h}))|\boldsymbol{u}|^{2}|+|\alpha(\tilde{\rho}_{h})(|\boldsymbol{u}|^{2}-|\boldsymbol{\tilde{u}}_{h}|^{2})|+2|\boldsymbol{f}\cdot(\boldsymbol{u}-\boldsymbol{\tilde{u}}_{h})|\;\mathrm{d}x\\\
&\indent\indent+\nu\sum_{K\in\mathcal{T}_{h}}\int_{K}\left||\nabla\boldsymbol{u}|^{2}-|\nabla\boldsymbol{\tilde{u}}_{h}|^{2}|\right|\;\mathrm{d}x\\\
&\indent\indent+\nu\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}\sigma
h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}]\\!]_{F}|^{2}\,\mathrm{d}s+\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}\sigma
h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}-\boldsymbol{g}_{h}]\\!]_{F}|^{2}\,\mathrm{d}s\\\
&\indent\indent+2\nu\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}|\\{\\!\\!\\{\nabla\boldsymbol{\tilde{u}}_{h}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{\tilde{u}}_{h}]\\!]_{F}|\,\mathrm{d}s\\\
&\indent\indent+2\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}|\\{\\!\\!\\{\nabla\boldsymbol{\tilde{u}}_{h}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{\tilde{u}}_{h}-\boldsymbol{g}_{h}]\\!]_{F}|\,\mathrm{d}s.\end{split}$

For all $\boldsymbol{v}\in H^{1}(\mathcal{T}_{h})^{d}$,
$\boldsymbol{\mathsf{\Phi}}\in(X_{\mathrm{DG}_{k}})^{d\times d}$, $h>0$, we
have that [16, Lem. 7]

(54)
$\displaystyle\begin{split}&\sum_{F\in\mathcal{F}_{h}}\int_{F}|\\{\\!\\!\\{\boldsymbol{\mathsf{\Phi}}\\}\\!\\!\\}_{F}:[\\![\boldsymbol{v}]\\!]_{F}|\,\mathrm{d}s\\\
&\leq
C\left(\sum_{F\in\mathcal{F}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{v}]\\!]_{F}|^{2}\,\mathrm{d}s\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}}\|\boldsymbol{\mathsf{\Phi}}\|^{2}_{L^{2}(K)}\right)^{1/2},\end{split}$

for a constant $C$ that only depends on the mesh quality. Hence, we see that

(55)
$\displaystyle\begin{split}&2|J_{h}(\boldsymbol{\tilde{u}}_{h},\tilde{\rho}_{h})-J(\boldsymbol{u},\rho)|\\\
&\indent\leq
L_{\alpha}\|\boldsymbol{u}\|^{2}_{L^{4}(\Omega)}\|\tilde{\rho}_{h}-\rho\|_{L^{2}(\Omega)}+2\|\boldsymbol{f}\|_{L^{2}(\Omega)}\|\boldsymbol{\tilde{u}}_{h}-\boldsymbol{u}\|_{L^{2}(\Omega)}\\\
&\indent\indent+\bar{\alpha}\|\boldsymbol{\tilde{u}}_{h}-\boldsymbol{u}\|_{L^{2}(\Omega)}(\|\boldsymbol{\tilde{u}}_{h}-\boldsymbol{u}\|_{L^{2}(\Omega)}+2\|\boldsymbol{u}\|_{L^{2}(\Omega)})\\\
&\indent\indent+\nu\sum_{K\in\mathcal{T}_{h}}\|\nabla\boldsymbol{\tilde{u}}_{h}-\nabla\boldsymbol{u}\|_{L^{2}(K)}(\|\nabla\boldsymbol{\tilde{u}}_{h}-\nabla\boldsymbol{u}\|_{L^{2}(K)}+2\|\nabla\boldsymbol{u}\|_{L^{2}(K)})\\\
&\indent\indent+\nu\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}\sigma
h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}]\\!]_{F}|^{2}\mathrm{d}s+\nu\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}\sigma
h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}-\boldsymbol{g}_{h}]\\!]_{F}|^{2}\mathrm{d}s\\\
&\indent\indent+C\nu\left(\sum_{F\in\mathcal{F}^{i}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}]\\!]_{F}|^{2}\mathrm{d}s\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}}\|\nabla\boldsymbol{\tilde{u}}_{h}\|^{2}_{L^{2}(K)}\right)^{1/2}\\\
&\indent\indent+C\nu\left(\sum_{F\in\mathcal{F}^{\partial}_{h}}\int_{F}h_{F}^{-1}|[\\![\boldsymbol{\tilde{u}}_{h}-\boldsymbol{g}_{h}]\\!]_{F}|^{2}\mathrm{d}s\right)^{1/2}\left(\sum_{K\in\mathcal{T}_{h}}\|\nabla\boldsymbol{\tilde{u}}_{h}\|^{2}_{L^{2}(K)}\right)^{1/2},\end{split}$

where $L_{\alpha}$ denotes the Lipschitz constant for $\alpha$. Thanks to the
strong convergence of $\boldsymbol{\tilde{u}}_{h}$ in the broken
$H^{1}_{\boldsymbol{g}}$-norm to $\boldsymbol{u}$ and by assumption (F1), from
Eq. 55 we deduce that

$\displaystyle J_{h}(\boldsymbol{\tilde{u}}_{h},\tilde{\rho}_{h})\to
J(\boldsymbol{u},\rho)\;\;\text{as}\;\;h\to 0.$

Furthermore, for sufficiently small $h>0$ we note that

$\displaystyle(\boldsymbol{\tilde{u}}_{h},\tilde{\rho}_{h})\in(\boldsymbol{U}_{\boldsymbol{g}_{h},h}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho)).$

Therefore, since $(\boldsymbol{u}_{h},\rho_{h})$ is a global minimizer of Eq.
BP-h in
$B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho))$,

(56) $\displaystyle J_{h}(\boldsymbol{u}_{h},\rho_{h})\leq
J_{h}(\boldsymbol{\tilde{u}}_{h},\tilde{\rho}_{h}).$

By taking the limit as $h\to 0$ and utilizing the strong convergence of
$\boldsymbol{\tilde{u}}_{h}$ and $\tilde{\rho}_{h}$ to $\boldsymbol{u}$ and
$\rho$, respectively, we see that

(57) $\displaystyle\lim_{h\to 0}J_{h}(\boldsymbol{u}_{h},\rho_{h})\leq
J(\boldsymbol{u},\rho).$

By assumption (F3), for every $q\in L^{2}_{0}(\Omega)$, there exists a
sequence of $\tilde{q}_{h}\in M_{h}$ such that $\tilde{q}_{h}\to q$ strongly
in $L^{2}(\Omega)$. Since
$\boldsymbol{u}_{h}\rightharpoonup\hat{\boldsymbol{u}}$ weakly in
$\boldsymbol{H}(\mathrm{div};\Omega)$ and
$\boldsymbol{u}_{h}\in\boldsymbol{U}_{\boldsymbol{g}_{h},h}$, we see that

(58) $\displaystyle b(\hat{\boldsymbol{u}},q)=\lim_{h\to
0}b(\boldsymbol{u}_{h},\tilde{q}_{h})+\lim_{h\to
0}b(\boldsymbol{u}_{h},q-\tilde{q}_{h})=0\;\;\text{for all}\;\;q\in
L^{2}_{0}(\Omega).$

Hence, $\hat{\boldsymbol{u}}$ is pointwise divergence-free. The final step to
identify $\hat{\boldsymbol{u}}$ as $\boldsymbol{u}$ is to show that
$\hat{\boldsymbol{u}}\in H^{1}_{\boldsymbol{g}}(\Omega)^{d}$. Now, the
sequence $(\boldsymbol{u}_{h})$ also defines a bounded sequence in
$H^{1}(\mathcal{T}_{h})^{d}$ such that

$\displaystyle\sup_{h>0}\left[\|\boldsymbol{u}_{h}\|_{L^{1}(\Omega)}+|\boldsymbol{u}_{h}|_{H^{1}(\mathcal{T}_{h})}\right]<+\infty.$

Hence, by the compact embedding lemma, as found in Buffa and Ortner [16, Lem.
8], there exists a subsequence (up to relabeling) and a limit
$\hat{\boldsymbol{w}}\in H^{1}(\Omega)^{d}$ such that

(59) $\displaystyle\boldsymbol{u}_{h}\to\hat{\boldsymbol{w}}\;\;\text{strongly
in}\;\;L^{q}(\Omega)^{d},$

where $1\leq q<\infty$ in two dimensions and $1\leq q<6$ in three dimensions.
By the uniqueness of limits $\hat{\boldsymbol{w}}=\hat{\boldsymbol{u}}$ a.e.
in $\Omega$ and thus $\hat{\boldsymbol{u}}\in H^{1}(\Omega)^{d}$. Moreover,
the same compact embedding lemma implies that

(60) $\displaystyle\boldsymbol{u}_{h}\to\hat{\boldsymbol{u}}\;\;\text{strongly
in}\;\;L^{r}(\partial\Omega)^{d},$

where $1\leq r<\infty$ in two dimensions and $1\leq r<4$ in three dimensions.
If $\|\boldsymbol{u}_{h}-\boldsymbol{g}\|_{L^{2}(\partial\Omega)}\not\to 0$,
then $J_{h}(\boldsymbol{u}_{h},\rho_{h})\to+\infty$. Since
$(\boldsymbol{u}_{h})$ is a bounded sequence, we must have that
$\|\boldsymbol{u}_{h}-\boldsymbol{g}\|_{L^{2}(\partial\Omega)}\to 0$. Hence,

(61)
$\displaystyle\|\hat{\boldsymbol{u}}-\boldsymbol{g}\|_{L^{2}(\partial\Omega)}$
$\displaystyle\leq\|\hat{\boldsymbol{u}}-\boldsymbol{u}_{h}\|_{L^{2}(\partial\Omega)}+\|\boldsymbol{u}_{h}-\boldsymbol{g}\|_{L^{2}(\partial\Omega)}\to
0.$

Thus, Eq. 58, Eq. 59, and Eq. 61 imply that $\hat{\boldsymbol{u}}\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$.

In order to identify the weak limit $(\boldsymbol{\hat{u}},\hat{\rho})$ with
the isolated minimizer $(\boldsymbol{u},\rho)$, we require a weak lower
semicontinuity result. Consider the following decomposition of the functional
$J_{h}(\boldsymbol{u}_{h},\rho_{h})=J_{1,h}(\boldsymbol{u}_{h},\rho_{h})+J_{2,h}(\boldsymbol{u}_{h})$
where

(62) $\displaystyle
J_{1,h}(\boldsymbol{u}_{h},\rho_{h})=\frac{1}{2}\int_{\Omega}\alpha(\rho_{h})|\boldsymbol{u}_{h}|^{2}-2\boldsymbol{f}\cdot\boldsymbol{u}_{h}\,\mathrm{d}x,$

and
$J_{2,h}(\boldsymbol{u}_{h})=J_{h}(\boldsymbol{u}_{h},\rho_{h})-J_{1,h}(\boldsymbol{u}_{h},\rho_{h})$.
It follows from Eq. 59 and a small modification to the proof in [13, Th. 3.1]
that

(63) $\displaystyle J_{1,h}(\boldsymbol{\hat{u}},\hat{\rho})\leq\liminf_{h\to
0}J_{1,h}(\boldsymbol{u}_{h},\rho_{h}).$

Moreover, it follows from a convergence result in Buffa and Ortner [16, Th.
6.1] that

(64) $\displaystyle
J_{2,h}(\boldsymbol{u}_{h})\to\frac{\nu}{2}\int_{\Omega}|\nabla\boldsymbol{\hat{u}}|^{2}\,\mathrm{d}x.$

Hence, we have the following weak lower semicontinuity result:

(65) $\displaystyle J(\boldsymbol{\hat{u}},\hat{\rho})\leq\liminf_{h\to
0}J_{h}(\boldsymbol{u}_{h},\rho_{h}).$

Since $(\boldsymbol{u},\rho)$ is the unique minimizer of Eq. BP in

$\displaystyle(\boldsymbol{H}_{\boldsymbol{g},\mathrm{div}}(\mathrm{div};\Omega)\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u}))\times(C_{\gamma}\cap
B_{r/2,L^{2}(\Omega)}(\rho)),$

we see that $J(\boldsymbol{u},\rho)\leq J(\boldsymbol{\hat{u}},\hat{\rho})$.
Hence, from Eq. 57 and Eq. 65, it follows that

(66) $\displaystyle
J(\boldsymbol{\hat{u}},\hat{\rho})=J(\boldsymbol{u},\rho).$

Since $(\boldsymbol{u},\rho)$ is the unique minimizer in the spaces we
consider, we identify the limits $\hat{\boldsymbol{u}}$ and $\hat{\rho}$ as
$\boldsymbol{u}$ and $\rho$, respectively, and state that
$\boldsymbol{u}_{h}\rightharpoonup\boldsymbol{u}$ weakly in
$\boldsymbol{H}(\mathrm{div};\Omega)$, $\boldsymbol{u}_{h}\to\boldsymbol{u}$
strongly in $L^{q}(\Omega)^{d}$, $\boldsymbol{u}_{h}\to\boldsymbol{u}$
strongly in $L^{r}(\partial\Omega)^{d}$ and $\rho_{h}\rightharpoonup\rho$
weakly in $L^{2}(\Omega)$, where $1\leq q,r<\infty$ in two dimensions and
$1\leq q<6$, $1\leq r<4$ in three dimensions.

We note that by the Banach–Alaoglu theorem [29, Th. A.52], the closed unit
ball of the dual space of a normed vector space, (for example
$L^{1}(\Omega)$), is compact in the weak-* topology. Hence we also find a
subsequence such that
$\rho_{h}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\hat{\rho}\in
C_{\gamma}\cap\\{\eta:\|\rho-\eta\|_{L^{\infty}(\Omega)}\leq r/2\\}$ weakly-*
in $L^{\infty}(\Omega)$. By the uniqueness of the weak limit, we identify
$\hat{\rho}=\rho$ a.e. in $\Omega$ and, thus, we deduce that
$\rho_{h}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho$ weakly-* in
$L^{\infty}(\Omega)$. Consequently, $\rho_{h}\rightharpoonup\rho$ weakly in
$L^{s}(\Omega)$ for all $s\in[1,\infty)$.

######  Proposition 4.8 (Strong convergence of $\rho_{h}$ in $L^{s}(\Omega)$,
$s\in[1,\infty)$).

There exists a subsequence of minimizers, $(\rho_{h})$, of Eq. BP-h such that

(67) $\displaystyle\rho_{h}\to\rho\;\text{strongly
in}\;L^{s}(\Omega),\;\;s\in[1,\infty).$

The proof of Proposition 4.8, with some small modifications, can be found in
[42, Prop. 5]. For the convenience of the reader, we reproduce the proof in
Appendix A with the necessary changes.

######  Proposition 4.9 (Strong convergence of $\boldsymbol{u}_{h}$ in the
${H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}$-norm).

There exists a subsequence of minimizers, $(\boldsymbol{u}_{h})$, of Eq. BP-h
such that

(68)
$\displaystyle\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0.$

######  Proof 4.10.

We note that $\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$ is a convex set,
and hence for any
$\boldsymbol{w}_{h}\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$, $t\in[0,1]$, we
have that
$\boldsymbol{u}_{h}+t(\boldsymbol{w}_{h}-\boldsymbol{u}_{h})\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$. Since
$(\boldsymbol{u}_{h},\rho_{h})$ is a global minimizer of Eq. BP-h, we note
that

(69)
$\displaystyle\frac{1}{t}\left[J_{h}(\boldsymbol{u}_{h}+t(\boldsymbol{w}_{h}-\boldsymbol{u}_{h}),\rho_{h})-J_{h}(\boldsymbol{u}_{h},\rho_{h})\right]\geq
0.$

By taking the limit $t\to 0$, a calculation shows that, for all
$\boldsymbol{w}_{h}\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$,

(70) $\displaystyle
a_{h}(\boldsymbol{u}_{h},\boldsymbol{w}_{h}-\boldsymbol{u}_{h};\rho_{h})\geq
l_{h}(\boldsymbol{w}_{h}-\boldsymbol{u}_{h};\boldsymbol{g}_{h}).$

We note that $\boldsymbol{w}_{h}-\boldsymbol{u}_{h}\in\boldsymbol{U}_{h,0}$.
Hence, from Proposition 3.7 and Lemma 3.12, we deduce that

(71) $\displaystyle
a_{h}(\boldsymbol{u},\boldsymbol{w}_{h}-\boldsymbol{u}_{h};\rho)=l_{h}(\boldsymbol{w}_{h}-\boldsymbol{u}_{h};\boldsymbol{g}).$

Therefore, from Eq. 70 and Eq. 71, we see that

(72)
$\displaystyle\begin{split}a_{h}(\boldsymbol{u}_{h},\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\rho_{h})&\leq
a_{h}(\boldsymbol{u},\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\rho)\\\
&\indent+l_{h}(\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\boldsymbol{g}_{h})-l_{h}(\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\boldsymbol{g}).\end{split}$

Hence, by subtracting
$a_{h}(\boldsymbol{w}_{h},\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\rho_{h})$
from both sides of Eq. 72, and utilizing the coercivity of
$a_{h}(\cdot,\cdot;\cdot)$ as stated in Proposition 3.9, we have that

(73)
$\displaystyle\begin{split}c_{a}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|^{2}_{H^{1}(\mathcal{T}_{h})}&\leq
a_{h}(\boldsymbol{u},\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\rho)-a_{h}(\boldsymbol{w}_{h},\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\rho_{h})\\\
&\indent+l_{h}(\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\boldsymbol{g}_{h})-l_{h}(\boldsymbol{u}_{h}-\boldsymbol{w}_{h};\boldsymbol{g}).\end{split}$

Now by assumption (M3), for all $F\in\mathcal{F}^{\partial}_{h}$ there exists
a $c>0$ such that $h^{-1}_{F}\leq ch^{-1}$, where $c$ depends on the mesh
regularity. By taking the absolute value of the right-hand side of Eq. 73,
collecting terms, utilizing the inequality Eq. 54, and the boundedness of
$a_{h}$ by $C_{a}$ as stated in Proposition 3.9, we have that

(74)
$\displaystyle\begin{split}c_{a}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|^{2}_{H^{1}(\mathcal{T}_{h})}&\leq\bar{\alpha}\|\boldsymbol{u}-\boldsymbol{w}_{h}\|_{L^{2}(\Omega)}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{L^{2}(\Omega)}\\\
&\indent+\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{L^{2}(\Omega)}\\\
&\indent+C_{a}\|\boldsymbol{u}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}\\\
&\indent+Ch^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{L^{2}(\partial\Omega)}\\\
&\indent+Ch^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})},\end{split}$

for some constant $C=C(\sigma)$ that also depends on the mesh regularity.

We note that
$\|\boldsymbol{u}-\boldsymbol{w}_{h}\|_{L^{2}(\Omega)}\leq\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}$
by definition. Moreover, by the broken trace theorem as found in Buffa and
Ortner [16, Th. 4.4], there exists a constant $C_{\mathrm{BT}}$ such that, for
all $\boldsymbol{v}\in H^{1}(\mathcal{T}_{h})^{d}$, $d\in\\{2,3\\}$ we have

(75) $\displaystyle\|\boldsymbol{v}\|_{L^{2}(\partial\Omega)}\leq
C_{\mathrm{BT}}\|\boldsymbol{v}\|_{H^{1}(\mathcal{T}_{h})}.$

Therefore, by bounding the $L^{2}(\Omega)$ and $L^{2}(\partial\Omega)$-norms
of $\boldsymbol{u}_{h}-\boldsymbol{w}_{h}$ above by the broken $H^{1}$-norm,
and dividing through by
$c_{a}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}$ we
see that

(76)
$\displaystyle\begin{split}\|\boldsymbol{u}_{h}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}&\leq
C\|\boldsymbol{u}-\boldsymbol{w}_{h}\|_{L^{2}(\Omega)}+C\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}\\\
&\indent+C\|\boldsymbol{u}-\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}+Ch^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)},\\\
\end{split}$

for some constant $C$ that depends on
$c_{a},C_{a},C_{\mathrm{BT}},\bar{\alpha},\sigma$ and the mesh regularity.

For sufficiently small $h$, we note that
$\boldsymbol{\tilde{u}}_{h}\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$ (where
$\boldsymbol{\tilde{u}}_{h}$ is defined in Lemma 3.14) and
$\boldsymbol{\tilde{u}}_{h}\to\boldsymbol{u}$ strongly in
$H^{1}(\mathcal{T}_{h})^{d}$. Fix
$\boldsymbol{w}_{h}=\boldsymbol{\tilde{u}}_{h}$.

From Proposition 4.8, we know that there exists a subsequence (not indicated)
such that $\rho_{h}\to\rho$ strongly in $L^{4}(\Omega)$. We now observe that

(77)
$\displaystyle\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}\leq
L_{\alpha}\|\rho-\rho_{h}\|_{L^{4}(\Omega)}\|\boldsymbol{u}\|_{L^{4}(\Omega)},$

where $L_{\alpha}$ is the Lipschitz constant for $\alpha$.
$\|\boldsymbol{u}\|_{L^{4}(\Omega)}$ is bounded for $d\in\\{2,3\\}$ thanks to
the Sobolev embedding theorem. Hence, by taking the limit as $h\to 0$ in Eq.
76, from (F1), Lemma 3.14, and Eq. 77, we deduce that
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\to 0$ as $h\to
0$. In Proposition 4.6, we showed that $\boldsymbol{u}_{h}\to\boldsymbol{u}$
strongly in $L^{2}(\partial\Omega)^{d}$. Hence, we conclude that
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$ as $h\to 0$.

######  Proposition 4.11 (Discretized first-order optimality conditions).

There exists an $\bar{h}>0$ such that for all $h<\bar{h}$, there exists a
unique Lagrange multiplier $p_{h}\in M_{h}$ such that the functions
$(\boldsymbol{u}_{h},\rho_{h})$, that minimize Eq. BP-h, satisfy the first-
order optimality conditions Eq. FOC1-h–Eq. FOC3-h.

######  Proof 4.12.

From Proposition 4.9, we know that $\boldsymbol{u}_{h}\to\boldsymbol{u}$
strongly in $H^{1}(\mathcal{T}_{h})^{d}$. Hence by definition of strong
convergence, there exists an $\bar{h}_{1}>0$ such that, for all
$h\leq\bar{h}_{1}$,
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\leq r/4$.
Moreover, since $\boldsymbol{u}\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}$, we have that
$\mathrm{div}(\boldsymbol{u})=0$ a.e. in $\Omega$ and by Lemma 3.12, we have
that $\mathrm{div}(\boldsymbol{u}_{h})=0$ a.e. in $\Omega$. Therefore,

(78)
$\displaystyle\begin{split}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}_{\boldsymbol{H}(\mathrm{div};\Omega)}&=\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}_{L^{2}(\Omega)}+\|\mathrm{div}(\boldsymbol{u}-\boldsymbol{u}_{h})\|^{2}_{L^{2}(\Omega)}\\\
&=\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}_{L^{2}(\Omega)}\leq\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}_{H^{1}(\mathcal{T}_{h})}\leq
r^{2}/16.\end{split}$

Hence, for each $\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,0}$, if
$|t|<r/(4\|\boldsymbol{v}_{h}\|_{H^{1}(\mathcal{T}_{h})})$ then
$\boldsymbol{u}_{h}+t\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,\boldsymbol{g}_{h}}\cap
B_{r/2,\boldsymbol{H}(\mathrm{div};\Omega)}(\boldsymbol{u})$. From Proposition
4.8 we have that $\rho_{h}\to\rho$ strongly in $L^{2}(\Omega)$. Hence, there
exists an $\bar{h}_{2}>0$ such that, for all $h\leq\bar{h}_{2}$,
$\|\rho-\rho_{h}\|_{L^{2}(\Omega)}\leq r/4$. Therefore, for each $\eta_{h}\in
C_{\gamma,h}$, if $0<t<r/(4\|\eta_{h}-\rho_{h}\|_{L^{2}(\Omega)})$ then
$\rho_{h}+t(\eta_{h}-\rho_{h})\in C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho)$. Let $\bar{h}=\min\\{\bar{h}_{1},\bar{h}_{2}\\}$
and consider $h\leq\bar{h}$.

Since $(\boldsymbol{u}_{h},\rho_{h})$ is a global minimizer of Eq. BP-h, then,
for all $\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,0}$, if
$|t|<r/(4\|\boldsymbol{v}_{h}\|_{H^{1}(\mathcal{T}_{h})})$ we have

(79)
$\displaystyle\frac{1}{t}\left[J_{h}(\boldsymbol{u}_{h}+t\boldsymbol{v}_{h},\rho_{h})-J_{h}(\boldsymbol{u}_{h},\rho_{h})\right]\geq
0.$

By considering the limits for $t\to 0_{+}$ and $t\to 0_{-}$, we have that, for
all $\boldsymbol{v}_{h}\in\boldsymbol{U}_{h,0}$,

(80) $\displaystyle
a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h};\rho_{h})=l_{h}(\boldsymbol{v}_{h};\boldsymbol{g}_{h}).$

From Eq. 80, the existence of a unique $p_{h}\in M_{h}$ such that
$(\boldsymbol{u}_{h},\rho_{h},p_{h})$ satisfy Eq. FOC1-h–Eq. FOC2-h follows
from the inf-sup condition (F2) and the argument can be found in [42, Prop.
2].

Similarly, since $(\boldsymbol{u}_{h},\rho_{h})$ is a global minimizer of Eq.
BP-h, then, for all $\eta_{h}\in C_{\gamma,h}$, if
$0<t<r/(4\|\eta_{h}-\rho_{h}\|_{L^{2}(\Omega)})$ we have

(81)
$\displaystyle\frac{1}{t}\left[J_{h}(\boldsymbol{u}_{h},\rho_{h}+t(\eta_{h}-\rho_{h}))-J_{h}(\boldsymbol{u}_{h},\rho_{h})\right]\geq
0.$

By taking the limit as $t\to 0$, we deduce that Eq. FOC3-h holds.

######  Proposition 4.13 (Strong convergence of $p_{h}$ in $L^{2}(\Omega)$).

There is a subsequence of the unique $p_{h}\in M_{h}$ defined in Proposition
4.11 that converges strongly in $L^{2}(\Omega)$ to the $p\in
L^{2}_{0}(\Omega)$ that solves Eq. FOC1–Eq. FOC3 for the given isolated
minimizer $(\boldsymbol{u},\rho)$.

######  Proof 4.14.

The inf-sup condition (F2) for $M_{h}$ and $\boldsymbol{X}_{h,0}$ implies
that, for any $q_{h}\in M_{h}$,

$\displaystyle c_{b}\|q_{h}-p_{h}\|_{L^{2}(\Omega)}$
$\displaystyle\leq\sup_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{b(\boldsymbol{w}_{h},q_{h}-p_{h})}{\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}}$
$\displaystyle=\sup_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{b(\boldsymbol{w}_{h},p-p_{h})+b(\boldsymbol{w}_{h},q_{h}-p)}{\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}}$
$\displaystyle\leq\sup_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{|b(\boldsymbol{w}_{h},p-p_{h})|+|b(\boldsymbol{w}_{h},q_{h}-p)|}{\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}}.$

From Proposition 3.7 and Proposition 4.11, it follows that

(82) $\displaystyle
b(\boldsymbol{w}_{h},p-p_{h})=a_{h}(\boldsymbol{u}_{h},\boldsymbol{w}_{h};\rho_{h})-a_{h}(\boldsymbol{u},\boldsymbol{w}_{h};\rho)+l_{h}(\boldsymbol{w}_{h};\boldsymbol{g})-l_{h}(\boldsymbol{w}_{h};\boldsymbol{g}_{h}).$

Therefore,

(83) $\displaystyle\begin{split}&c_{b}\|q_{h}-p_{h}\|_{L^{2}(\Omega)}\\\
&\leq\sup_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{|a_{h}(\boldsymbol{u}_{h},\boldsymbol{w}_{h};\rho_{h})-a_{h}(\boldsymbol{u},\boldsymbol{w}_{h};\rho)|}{\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}}\\\
&\indent+\sup_{\boldsymbol{w}_{h}\in\boldsymbol{X}_{h,0}\backslash\\{0\\}}\frac{|l_{h}(\boldsymbol{w}_{h};\boldsymbol{g})-l_{h}(\boldsymbol{w}_{h};\boldsymbol{g}_{h})|+|b(\boldsymbol{w}_{h},q_{h}-p)|}{\|\boldsymbol{w}_{h}\|_{H^{1}(\mathcal{T}_{h})}}.\end{split}$

By using the same argument we used to bound Eq. 74 by Eq. 76, we see that Eq.
83 implies that

(84) $\displaystyle\begin{split}c_{b}\|q_{h}-p_{h}\|_{L^{2}(\Omega)}&\leq
C\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}+C\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\\\
&\indent+Ch^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}+C_{b}\|p-q_{h}\|_{L^{2}(\Omega)}.\end{split}$

where $C_{b}$ is the boundedness constant for $b(\cdot,\cdot)$ and $C$ is
dependent on $C_{a}$, $C_{\mathrm{BT}}$, $\bar{\alpha}$, $L_{\alpha}$, $\nu$,
$\sigma$ and the mesh regularity. Hence, by an application on the
Cauchy–Schwarz inequality,

(85) $\displaystyle\begin{split}\|p-p_{h}\|_{L^{2}(\Omega)}&\leq
C\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}+C\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\\\
&\indent+Ch^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}+C\|p-q_{h}\|_{L^{2}(\Omega)},\end{split}$

where $C$ is dependent on $c_{b}$, $C_{b}$, $C_{a}$, $C_{\mathrm{BT}}$,
$\bar{\alpha}$, $L_{\alpha}$, $\nu$, $\sigma$ and the mesh regularity.

By assumption (F3), there exists a sequence of finite element functions,
$\tilde{p}_{h}\in M_{h}$ that satisfies $\tilde{p}_{h}\to p$ strongly in
$L^{2}(\Omega)$. Let $q_{h}=\tilde{p}_{h}$. We have already shown that
$\boldsymbol{u}_{h}\to\boldsymbol{u}$ strongly in $H^{1}(\mathcal{T}_{h})^{d}$
in Proposition 4.9. Similarly, in the proof of Proposition 4.9 we also showed
that $\|(\alpha(\rho)-\alpha(\rho_{h}))\boldsymbol{u}\|_{L^{2}(\Omega)}\to 0$.
By assumption (F1),
$h^{-1}\|\boldsymbol{g}-\boldsymbol{g}_{h}\|_{L^{2}(\partial\Omega)}\to 0$.
Hence, we conclude that $p_{h}\to p$ strongly in $L^{2}(\Omega)$.

We now have the required results to prove Theorem 4.1.

######  Proof 4.15 (Proof of Theorem 4.1).

Fix an isolated minimizer $(\boldsymbol{u},\rho)$ of Eq. BP and its unique
associated Lagrange multiplier $p$. By the results of Propositions 4.6, 4.8,
4.9, and 4.11, there exists a mesh size $\bar{h}$ such that for, $h<\bar{h}$,
there exists a sequence of finite element solutions
$(\boldsymbol{u}_{h},\rho_{h},p_{h})\in\boldsymbol{X}_{h,\boldsymbol{g}_{h}}\times
C_{\gamma,h}\times M_{h}$ satisfying Eq. FOC1-h–Eq. FOC3-h that converges to
$(\boldsymbol{u},\rho,p)$. By taking a subsequence if necessary (not
indicated), Proposition 4.8, implies that $\rho_{h}\to\rho$ strongly in
$L^{s}(\Omega)$, $s\in[1,\infty)$, Proposition 4.9 implies that
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}_{\boldsymbol{g}}(\mathcal{T}_{h})}\to
0$, and Proposition 4.13 implies that $p_{h}\to p$ strongly in
$L^{2}(\Omega)$.

##  5 Numerical results

In this section we consider an example of a Borrvall–Petersson topology
optimization problem that supports two minimizers. Our goal is to discretize
the velocity-pressure pair with a Brezzi–Douglas–Marini
$\mathrm{BDM}_{1}\times\mathrm{DG}_{0}$ discretization, and the material
distribution with a piecewise constant $\mathrm{DG}_{0}$ discretization, and
numerically verify the existence of two sequences of solutions to Eq.
FOC1-h–Eq. FOC3-h that converge to the two different minimizers. Moreover, we
compare the violation of the incompressibility constraint, measured by
$\|\mathrm{div}(\boldsymbol{u}_{h})\|_{L^{2}(\Omega)}$, with the equivalent
solution as computed via a Taylor–Hood
$(\mathrm{CG}_{2})^{2}\times\mathrm{CG}_{1}$ discretization for the velocity-
pressure pair on the same meshes.

We note that $\mathrm{DG}_{0}\not\subset C_{\gamma}$ and
$\mathrm{DG}_{0}\not\subset L^{2}_{0}(\Omega)$, as, in general,
$\eta\in\mathrm{DG}_{0}$ does not satisfy $0\leq\eta\leq 1$ a.e. and
$\int_{\Omega}\eta\,\mathrm{d}x\leq\gamma|\Omega|$, and, in general,
$q\in\mathrm{DG}_{0}$ does not satisfy $\int_{\Omega}q\,\mathrm{d}x=0$.
However, the choice of optimization algorithm (described below) only finds
solutions $\rho_{h}\in C_{\gamma}$ and $p_{h}\in L^{2}_{0}(\Omega)$. Hence, we
are effectively working with the correct conforming finite element spaces for
the pressure and material distribution and, hence, our discretization
restricted by the optimization strategy satisfies the conditions of Theorem
4.1.

The problem is implemented with the finite element software Firedrake [45] and
the computational domains are triangulated with simplices. The solutions are
computed using the deflated barrier method [41]. The deflated barrier method
reformulates Eq. FOC1-h–Eq. FOC3-h into a mixed complementarity problem and
solves the nonlinear system with a primal-dual active set solver that enforces
the box constraints on the material distribution [12]. The volume constraint
is enforced via a one-dimensional Lagrange multiplier. The zero mean value
constraint on the pressure is either enforced via a one-dimensional Lagrange
multiplier or by orthogonalizing against the nullspace of constants. The
global nonlinear convergence is aided by the continuation of barrier terms. A
key feature of the deflated barrier method is that it can systematically
discover multiple solutions of topology optimization problems by utilizing the
deflation technique [24, 25]. In the BDM discretization, the linear systems
arising in the deflated barrier method are solved with FGMRES [47]
preconditioned with block preconditioning techniques and the Schur complements
are controlled with an augmented Lagrangian term [40]. The individual blocks
are solved by a sparse LU factorization with MUMPS [7] and PETSc [10]. In the
Taylor–Hood discretization, we invert the entire linear system directly with
MUMPS. There are no known solutions of the infinite-dimensional problem for
choices of the inverse permeability, $\alpha$, used in practice. Hence, the
errors are measured with respect to the most heavily-refined finite element
solution where $h=1.41\times 10^{-3}$ resulting in 16,389,121 degrees of
freedom for the first-order BDM discretization.

Although, for each isolated minimizer, Theorem 4.1 guarantees the existence of
a converging sequence of solutions, it does not guarantee that the sequence is
unique. In principle, there can be infinitely many different subsequences of
finite element solutions that strongly converge to the same minimizer of the
infinite-dimensional problem at different convergence rates. Separate
subsequences can appear in oscillations in the error and cause difficulty in
the interpretation of the convergence plots [42, Sec. 5]. Here, we attempt to
find solutions along the same sequence by first computing the solutions on a
coarse mesh, uniformly refining the mesh, and successively interpolating the
solutions onto the finer mesh as initial guesses for the deflated barrier
method.

We now describe the double-pipe problem of Borrvall and Petersson [13, Sec.
4.5] with a modification to the boundary conditions that improves the
guaranteed regularity of the solutions. The design domain is a rectangle
$\Omega=(0,3/2)\times(0,1)$ with two inlets and two outlets. The volume
fraction is chosen to be $\gamma=1/3$ and the inverse permeability $\alpha$ is
given by Eq. 1 with $q=1/10$ and $\bar{\alpha}=2.5\times 10^{4}$. The boundary
conditions are given by:

(86)
$\displaystyle\boldsymbol{g}(x,y)=\begin{cases}\left(\exp(1-\frac{1}{1-(12y-9)^{2}}),0\right)^{\top}&\text{if}\;\;2/3\leq
y\leq 5/6,\,x=0\;\text{or}\;3/2,\\\
\left(\exp(1-\frac{1}{1-(12y-3)^{2}}),0\right)^{\top}&\text{if}\;\;1/6\leq
y\leq 1/3,\,x=0\;\text{or}\;3/2,\\\ (0,0)^{\top}&\text{elsewhere
on}\;\partial\Omega.\end{cases}$

This problem supports exactly two isolated minimizers in form of a straight
channels solution and double-ended wrench solution as shown in Fig. 1.

Figure 1: The material distribution of the straight channels (left) and
double-ended wrench (right) solutions of the double-pipe optimization problem.
The black regions represent where $\rho=0$ and the white regions are where
$\rho=1$. The mesh size is $h=1.41\times 10^{-3}$ and a first-order BDM
discretization for the velocity-pressure and a discontinuous piecewise
constant discretization for the material distribution are used. This results
in 16,389,121 degrees of freedom.

It can be checked that the boundary datum $\boldsymbol{g}$ can be expressed as
the trace of a function $\hat{\boldsymbol{g}}\in H^{2}(\Omega)^{d}$. Hence,
since the domain is convex and the forcing term is smooth, we have that, by
regularity results for isolated minimizers of the Borrvall–Petersson problem
[42, Lem. 5], $(\boldsymbol{u},p)\in H^{2}(\Omega)^{d}\times H^{1}(\Omega)$
for both minimizers. Hence, the trace of $\nabla\boldsymbol{u}$ is well-
defined on the faces of each element and the consistency result in Proposition
3.7 holds. All the conditions of Theorem 4.1 are satisfied and hence there
exists a sequence of solutions to Eq. FOC1-h–Eq. FOC3-h that converges
strongly to the straight channel solution and a different sequence of
solutions that converges to the double-ended wrench. The existence of these
sequences are numerically verified in Fig. 2 for a
$\mathrm{DG}_{0}\times\mathrm{BDM}_{1}\times\mathrm{DG}_{0}$ discretization
for $(\rho_{h},\boldsymbol{u}_{h},p_{h})$.

Figure 2: The convergence of $\boldsymbol{u}_{h}$, $\rho_{h}$, and $p_{h}$ for
the double-pipe problem for both the straight channels and double-ended wrench
solutions on a sequence of uniformly refined meshes with a
$\mathrm{DG}_{0}\times\mathrm{BDM}_{1}\times\mathrm{DG}_{0}$ discretization
for $(\rho_{h},\boldsymbol{u}_{h},p_{h})$.

We report the values of $\|\mathrm{div}(\boldsymbol{u}_{h})\|_{L^{2}(\Omega)}$
in Table 1 for the BDM discretization alongside the equivalent solutions
computed with a Taylor–Hood $(\mathrm{CG}_{2})^{2}\times\mathrm{CG}_{1}$
discretization for the velocity-pressure pair and a $\mathrm{DG}_{0}$
discretization for the material distribution on the same meshes. Even on
coarse meshes, the $L^{2}$-norm of the divergence of the velocity in the BDM
discretization is small with values in the range of $10^{-6}\sim 10^{-9}$ for
both minimizers. Many of the BDM discretization values are roughly the square
root of Float64 machine precision, denoted eps. Suppose that Eq. FOC2-h is
satisfied up to machine precision, then we have that
$|b(\boldsymbol{u}_{h},q_{h})|\leq\mathrm{\texttt{eps}}$. Now, by choosing
$q_{h}=\mathrm{div}(\boldsymbol{u}_{h})$, we note that
$\|\mathrm{div}(\boldsymbol{u}_{h})\|_{L^{2}(\Omega)}\leq\sqrt{\mathrm{\texttt{eps}}}$.
By contrast, the pointwise violation of the incompressibility constraint for
the Taylor–Hood discretization manifests as relatively large values of
$\|\mathrm{div}(\boldsymbol{u}_{h})\|_{L^{2}(\Omega)}$. Even on the finest
mesh, where $h=2.82\times 10^{-3}$ resulting in 4,512,004 degrees of freedom,
the $L^{2}$-norm is still $\mathcal{O}(10^{-3})$, 4 orders of magnitude larger
than the equivalent BDM discretization.

| Straight channels | Double-ended wrench
---|---|---
$h$ | BDM | Taylor–Hood | BDM | Taylor–Hood
$4.51\times 10^{-2}$ | $1.00\times 10^{-8}$ | $2.49\times 10^{-1}$ | $2.69\times 10^{-6}$ | $3.25\times 10^{-1}$
$2.25\times 10^{-2}$ | $6.35\times 10^{-9}$ | $1.09\times 10^{-1}$ | $2.75\times 10^{-8}$ | $1.35\times 10^{-1}$
$1.13\times 10^{-2}$ | $1.59\times 10^{-7}$ | $3.95\times 10^{-2}$ | $2.62\times 10^{-8}$ | $4.66\times 10^{-2}$
$5.63\times 10^{-3}$ | $4.19\times 10^{-8}$ | $1.19\times 10^{-2}$ | $1.48\times 10^{-7}$ | $1.36\times 10^{-2}$
$2.82\times 10^{-3}$ | $4.97\times 10^{-7}$ | $3.17\times 10^{-3}$ | $2.98\times 10^{-7}$ | $3.58\times 10^{-3}$
Table 1: Reported values for
$\|\mathrm{div}(\boldsymbol{u}_{h})\|_{L^{2}(\Omega)}$ in a BDM and
Taylor–Hood discretization for the double-pipe problem as measured on five
meshes in a uniformly refined mesh hierarchy.

Code availability: For reproducibility, the implementation of the deflated
barrier method used in this work, as well as scripts to generate the
convergence plots and solutions can be found at
https://github.com/ioannisPApapadopoulos/fir3dab/. The version of the software
used in this paper is archived on Zenodo [1].

##  6 Conclusions

In this work we studied the convergence of a divergence-free discontinuous
Galerkin finite element discretization of the fluid topology optimization
model of Borrvall and Petersson [13]. Our approach extends the techniques used
by Papadopoulos and Süli [42] for $H^{1}$-conforming finite element
approximations of the velocity. The nonconvexity of the optimization problem
was handled by fixing any isolated minimizer and introducing a modified
optimization problem with the chosen isolated minimizer as its unique
solution. We then showed that there exists a sequence of discretized solutions
that converges to the minimizer in the appropriate norms. In particular,
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\to 0$,
$\|\rho-\rho_{h}\|_{L^{s}(\Omega)}\to 0$, $s\in[1,\infty)$, and
$\|p-p_{h}\|_{L^{2}(\Omega)}\to 0$. The modified optimization problem was
related back to the original optimization problem by showing that a
subsequence of the strongly converging minimizers also satisfy the first-order
optimality conditions of the original problem. Moreover, these first-order
optimality conditions can be solved numerically. Finally, we numerically
verified that these sequences exist and compared the improvement in the
$L^{2}$-norm of the incompressibility constraint in the discretized solutions.
Future work could include adapting these methods to other topology
optimization problems including extensions of the Borrvall–Petersson problem
with a more sophisticated fluid flow, as well as topology optimization
formulations for cantilevers and MBB beams that utilize linear elasticity.

##  Appendix A Proof of Proposition 4.8

We first quote a couple of propositions that are required in the proof of
Proposition 4.8. The proof of the following proposition concerning the support
of $\rho$ can be found in Papadopoulos and Süli [42, Prop. 3].

######  Proposition A.1 (Support of $\rho$).

Suppose that $\Omega\subset\mathbb{R}^{d}$ is a Lipschitz domain,
$d\in\\{2,3\\}$, and $\alpha$ satisfies properties (A1)–(A4). Further assume
that the minimizer $(\boldsymbol{u},\rho)\in
H^{1}_{\boldsymbol{g},\mathrm{div}}(\Omega)^{d}\times C_{\gamma}$ of (BP) is a
strict minimizer. Then, $\mathrm{supp}(\rho)\subseteq U$, where
$U\coloneqq\mathrm{supp}(\boldsymbol{u})$.

######  Remark A.2.

Proposition A.1 is proved by contradiction. Suppose that there exists a strict
minimizer $(\boldsymbol{u},\rho)$ such that $\mathrm{supp}(\rho)\not\subseteq
U$. Then, it is possible to construct a pair $(\boldsymbol{u},\tilde{\rho})$
arbitrarily close to $(\boldsymbol{u},\rho)$ such that
$J(\boldsymbol{u},\rho)=J(\boldsymbol{u},\tilde{\rho})$. Hence,
$(\boldsymbol{u},\rho)$ cannot be a strict minimizer.

The following proposition concerning strong convergence of $\rho_{h}$ in sets
where $\rho=0$ or $\rho=1$ a.e. is thanks to Petersson [43, Cor. 3.2] and can
be found, as stated, in Papadopoulos and Süli [42, Cor. 1].

######  Proposition A.3 (Strong convergence of $\rho_{h}$ in
$L^{s}(\Omega_{b})$).

Fix an isolated minimizer $(\boldsymbol{u},\rho)$ of Eq. BP and suppose that
the conditions of Theorem 4.1 hold. Let $\Omega_{b}$ be any measurable subset
of $\Omega$ of positive measure on which $\rho$ is equal to zero or one a.e.
(if such a set exists). Suppose that there exists a sequence of finite element
minimizers, $\rho_{h}$, of Eq. BP-h such that
$\rho_{h}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho$ weakly-* in
$L^{\infty}(\Omega)$. Then, $\rho_{h}\to\rho$ strongly in $L^{s}(\Omega_{b})$,
where $s\in[1,\infty)$.

######  Proof A.4.

If $\rho_{h}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho$ weakly-* in
$L^{\infty}(\Omega)$, then by definition,
$\int_{\Omega}\rho_{h}\phi\,\mathrm{d}x\to\int_{\Omega}\rho\phi\,\mathrm{d}x$
for all $\phi\in L^{1}(\Omega)$. Consider
$\Omega_{b}=\Omega_{0}\cup\Omega_{1}$ where $\Omega_{0}\coloneqq\\{\rho=0\\}$
and $\Omega_{1}\coloneqq\\{\rho=1\\}$. Suppose that $|\Omega_{0}|>0$. Choose
$\phi=\chi_{\Omega_{0}}\in L^{1}(\Omega)$ (the characteristic function for
$\Omega_{0}$). Then,

(87)
$\displaystyle\int_{\Omega_{0}}\rho_{h}\,\mathrm{d}x\to\int_{\Omega_{0}}\rho\,\mathrm{d}x=0,$

as $\rho=0$ a.e. in $\Omega_{0}$. Hence,
$\int_{\Omega_{0}}|\rho_{h}-\rho|\mathrm{d}x\to 0$ as $0\leq\rho_{h}$.
Similarly, if $|\Omega_{1}|>0$, by choosing $\phi=\chi_{\Omega_{1}}\in
L^{1}(\Omega)$ and utilizing that $\rho_{h}\leq 1$, we find that
$\int_{\Omega_{1}}|\rho-\rho_{h}|\mathrm{d}x\to 0$. Therefore,
$\|\rho-\rho_{h}\|_{L^{1}(\Omega_{b})}\to 0$. Consider any $s\in(1,\infty)$.
Then,

(88)
$\displaystyle\int_{\Omega_{b}}|\rho-\rho_{h}|^{s}\mathrm{d}x\leq\int_{\Omega_{b}}|\rho-\rho_{h}|^{s-1}|\rho-\rho_{h}|\mathrm{d}x\leq\|\rho-\rho_{h}\|_{L^{1}(\Omega_{b})}\to
0,$

where the second inequality holds since $0\leq\rho,\rho_{h}\leq 1$. Therefore,
$\|\rho-\rho_{h}\|_{L^{s}(\Omega_{b})}\to 0$ for any $s\in[1,\infty)$.

We now reproduce the proof of Proposition 4.8 as found in [42, Prop. 5], with
some small modifications.

######  Proof A.5 (Proof of Proposition 4.8).

We note that $C_{\gamma,h}\cap B_{r/2,L^{2}(\Omega)}(\rho)$ is a convex set,
and hence for any $\eta_{h}\in C_{\gamma,h}\cap B_{r/2,L^{2}(\Omega)}(\rho)$,
$t\in[0,1]$, we have that $\rho_{h}+t(\eta_{h}-\rho_{h})\in C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho)$. Since $(\boldsymbol{u}_{h},\rho_{h})$ is a
global minimizer of Eq. BP-h, we note that

(89)
$\displaystyle\frac{1}{t}\left[J_{h}(\boldsymbol{u}_{h},\rho_{h}+t(\eta_{h}-\rho_{h}))-J_{h}(\boldsymbol{u}_{h},\rho_{h})\right]\geq
0.$

By taking the limit $t\to 0$ and noting that by assumption (A4), $\alpha$ is
continuously differentiable, we deduce that

(90)
$\displaystyle\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\eta_{h}-\rho_{h})\mathrm{d}x\geq
0\;\;\text{for all}\;\;\eta_{h}\in C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho).$

Hence, Eq. FOC3 and Eq. 90 imply that for all $\eta\in C_{\gamma}$ and
$\eta_{h}\in C_{\gamma,h}\cap B_{r/2,L^{2}(\Omega)}(\rho)$ we have that

(91)
$\displaystyle\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}\rho\;\mathrm{d}x$
$\displaystyle\leq\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}\eta\;\mathrm{d}x,$
(92)
$\displaystyle\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}\rho_{h}\;\mathrm{d}x$
$\displaystyle\leq\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}\eta_{h}\;\mathrm{d}x.\;\;$

By subtracting
$\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}\rho_{h}\mathrm{d}x$
from Eq. 91 and
$\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}\rho\,\mathrm{d}x$
from Eq. 92, we see that

(93)
$\displaystyle\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}(\rho-\rho_{h})\;\mathrm{d}x$
$\displaystyle\leq\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}(\eta-\rho_{h})\mathrm{d}x,$
(94)
$\displaystyle\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\rho_{h}-\rho)\;\mathrm{d}x$
$\displaystyle\leq\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\eta_{h}-\rho)\;\mathrm{d}x.$

Summing Eq. 93 and Eq. 94 and rearranging the left-hand side, we see that

(95)
$\displaystyle\begin{split}&\int_{\Omega}(\alpha^{\prime}(\rho)-\alpha^{\prime}(\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})\mathrm{d}x+\int_{\Omega}\alpha^{\prime}(\rho_{h})(|\boldsymbol{u}|^{2}-|\boldsymbol{u}_{h}|^{2})(\rho-\rho_{h})\mathrm{d}x\\\
&\indent\leq\int_{\Omega}\alpha^{\prime}(\rho)|\boldsymbol{u}|^{2}(\eta-\rho_{h})\mathrm{d}x+\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\eta_{h}-\rho)\mathrm{d}x.\end{split}$

By fixing $\eta=\rho_{h}\in C_{\gamma}$ and subtracting the second term on the
left-hand side of Eq. 95 from both sides we deduce that

(96)
$\displaystyle\begin{split}&\int_{\Omega}(\alpha^{\prime}(\rho)-\alpha^{\prime}(\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})\mathrm{d}x\\\
&\indent\leq\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\eta_{h}-\rho)\mathrm{d}x+\int_{\Omega}\alpha^{\prime}(\rho_{h})(|\boldsymbol{u}_{h}|^{2}-|\boldsymbol{u}|^{2})(\rho-\rho_{h})\mathrm{d}x.\end{split}$

By an application of the mean value theorem, we note that there exists a
$c\in(0,1)$ such that

(97)
$\displaystyle\begin{split}&\int_{\Omega}(\alpha^{\prime}(\rho)-\alpha^{\prime}(\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})\mathrm{d}x\\\
&\indent=\int_{\Omega}\alpha^{\prime\prime}(\rho_{h}+c(\rho-\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})^{2}\mathrm{d}x.\end{split}$

Since by assumption (A2), $\alpha$ is strongly convex and by (A4) it is twice
continuously differentiable, there exists a constant
$\alpha^{\prime\prime}_{\mathrm{min}}>0$ such that

(98)
$\displaystyle\alpha^{\prime\prime}_{\mathrm{min}}\leq\alpha^{\prime\prime}(y)\;\;\text{for
all}\;\;y\in[0,1].$

Therefore by Eq. 98 and the definition of $U_{\theta}$ (given in Theorem 4.1)
we bound Eq. 97 from below:

(99)
$\displaystyle\begin{split}&\int_{\Omega}\alpha^{\prime\prime}(\rho_{h}+c(\rho-\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})^{2}\mathrm{d}x\\\
&\indent\geq\int_{U_{\theta}}\alpha^{\prime\prime}(\rho_{h}+c(\rho-\rho_{h}))|\boldsymbol{u}|^{2}(\rho-\rho_{h})^{2}\mathrm{d}x\geq\alpha^{\prime\prime}_{\mathrm{min}}\theta\|\rho-\rho_{h}\|^{2}_{L^{2}(U_{\theta})}.\end{split}$

Now we bound the right-hand side of Eq. 96 as follows,

(100)
$\displaystyle\begin{split}&\int_{\Omega}\alpha^{\prime}(\rho_{h})|\boldsymbol{u}_{h}|^{2}(\eta_{h}-\rho)\mathrm{d}x+\int_{\Omega}\alpha^{\prime}(\rho_{h})(|\boldsymbol{u}_{h}|^{2}-|\boldsymbol{u}|^{2})(\rho-\rho_{h})\mathrm{d}x\\\
&\indent\leq
2\alpha^{\prime}_{\text{max}}(\|\boldsymbol{u}\|^{2}_{L^{4}(\Omega)}+\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}_{L^{4}(\Omega)})\|\rho-\eta_{h}\|_{L^{2}(\Omega)}\\\
&\indent\indent+\alpha^{\prime}_{\text{max}}\|\rho-\rho_{h}\|_{L^{q}(\Omega)}\|\boldsymbol{u}+\boldsymbol{u}_{h}\|_{L^{q^{\prime}}(\Omega)}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)},\end{split}$

where $2<q^{\prime}<\infty$ in two dimensions, $2<q^{\prime}\leq 6$ in three
dimensions, and $q=2q^{\prime}/(q^{\prime}-2)$. We note that

(101)
$\displaystyle\begin{split}&\|\boldsymbol{u}+\boldsymbol{u}_{h}\|_{L^{q^{\prime}}(\Omega)}\leq\|\boldsymbol{u}\|_{L^{q^{\prime}}(\Omega)}+\|\boldsymbol{u}_{h}\|_{L^{q^{\prime}}(\Omega)}\\\
&\indent\leq\|\boldsymbol{u}\|_{H^{1}(\Omega)}+\|\boldsymbol{u}_{h}\|_{H^{1}(\mathcal{T}_{h})}\leq\hat{C}<\infty,\end{split}$

where the second inequality holds thanks to the Sobolev embedding theorem and
the broken Friedrichs-type inequality as found in Buffa and Ortner [16, Cor.
4.3]. Combining Eq. 96–Eq. 101 we see that

(102) $\displaystyle\|\rho-\rho_{h}\|^{2}_{L^{2}(U_{\theta})}\leq
C\left(\|\rho-\eta_{h}\|_{L^{2}(\Omega)}+\|\rho-\rho_{h}\|_{L^{q}(\Omega)}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)}\right),$

where
$C=C(\alpha^{\prime}_{\text{max}},\alpha^{\prime\prime}_{\mathrm{min}},\theta,\|\boldsymbol{u}\|_{L^{4}(\Omega)},\hat{C})$.
By assumption (F3), there exists a sequence of finite element functions
$\tilde{\rho}_{h}\in C_{\gamma,h}$ such that $\tilde{\rho}_{h}\to\rho$
strongly in $L^{2}(\Omega)$. Thanks to the strong convergence, we note that
for sufficiently small $h$, $\tilde{\rho}_{h}\in C_{\gamma,h}\cap
B_{r/2,L^{2}(\Omega)}(\rho)$. Hence we fix $\eta_{h}=\tilde{\rho}_{h}$. By
Proposition 4.6, we know that $\boldsymbol{u}_{h}\to\boldsymbol{u}$ strongly
in $L^{2}(\Omega)^{d}$ and since $\rho\in C_{\gamma}$, $\rho_{h}\in
C_{\gamma,h}\subset C_{\gamma}$, then
$\|\rho-\rho_{h}\|_{L^{q}(\Omega)}\leq|\Omega|^{1/q}\|\rho-\rho_{h}\|_{L^{\infty}(\Omega)}\leq|\Omega|^{1/q}$.
Therefore, the right-hand side of Eq. 102 tends to zero as $h\to 0$. Hence, we
deduce that

(103) $\displaystyle\rho_{h}\to\rho\;\;\text{strongly
in}\;\;L^{2}(U_{\theta}),\;\;\theta>0.$

We define $U$ as $U\coloneqq\mathrm{supp}(\boldsymbol{u})$. Now we note that

(104)
$\displaystyle\|\rho-\rho_{h}\|_{L^{2}(\Omega)}=\|\rho-\rho_{h}\|_{L^{2}(U_{\theta})}+\|\rho-\rho_{h}\|_{L^{2}(U\backslash
U_{\theta})}+\|\rho-\rho_{h}\|_{L^{2}(\Omega\backslash U)}.$

If $U\backslash U_{\theta}$ or $\Omega\backslash U$ are empty, we neglect the
corresponding term in Eq. 104 with no loss of generality. Suppose
$\Omega\backslash U$ is non-empty. By definition of $U$,
$\boldsymbol{u}=\boldsymbol{0}$ a.e. in $\Omega\backslash U$. Proposition A.1
implies that $\rho=0$ a.e. in $\Omega\backslash U$. Hence, $\Omega\backslash
U\subseteq\Omega_{b}$ where $\Omega_{b}=\\{\rho=1\\}\cup\\{\rho=0\\}$.
Moreover, in Proposition 4.6 we showed that
$\rho_{h}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}\rho$ weakly-* in
$L^{\infty}(\Omega)$. Therefore, Proposition A.3 implies that

(105) $\displaystyle\rho_{h}\to\rho\;\;\text{strongly
in}\;\;L^{2}(\Omega\backslash U).$

Suppose $U\backslash U_{\theta}$ is non-empty. Since, $\rho,\rho_{h}\in
C_{\gamma}$ we see that

(106) $\displaystyle\|\rho-\rho_{h}\|_{L^{2}(U\backslash
U_{\theta})}\leq|U\backslash U_{\theta}|^{1/2}\to 0\;\;\text{as}\;\;\theta\to
0.$

Therefore, by first taking the limit as $h\to 0$ and then by taking the limit
as $\theta\to 0$, Eq. 103–Eq. 106 imply that $\rho_{h}\to\rho$ strongly in
$L^{2}(\Omega)$.

Since
$\|\rho-\rho_{h}\|_{L^{1}(\Omega)}\leq|\Omega|^{1/2}\|\rho-\rho_{h}\|_{L^{2}(\Omega)}$,
we see that $\rho_{h}\to\rho$ strongly in $L^{1}(\Omega)$. Hence, for any
$s\in[1,\infty)$,

(107)
$\displaystyle\int_{\Omega}|\rho-\rho_{h}|^{s}\mathrm{d}x=\int_{\Omega}|\rho-\rho_{h}|^{s-1}|\rho-\rho_{h}|\mathrm{d}x\leq
1^{s-1}\|\rho-\rho_{h}\|_{L^{1}(\Omega)},$

which implies that $\rho_{h}\to\rho$ strongly in $L^{s}(\Omega)$.

## Acknowledgments

The author would like to thank Pablo Alexei Gazca-Orozco for a useful
discussion on weak compactness results, and Endre Süli and Patrick Farrell for
their comments on the manuscript. The author would also like to thank the
anonymous reviewers for their insightful comments.

## References

  * [1] Software used in ‘Numerical analysis of a discontinuous Galerkin method for the Borrvall–Petersson topology optimization problem’, 2021, https://doi.org/10.5281/zenodo.5146324.
  * [2] R. A. Adams and J. J. Fournier, Sobolev spaces, Elsevier, second ed., 2003.
  * [3] J. Alexandersen and C. S. Andreasen, A review of topology optimisation for fluid-based problems, Fluids, 5 (2020), p. 29, https://doi.org/10.3390/fluids5010029.
  * [4] G. Allaire, Shape optimization by the homogenization method, vol. 146, Springer Science & Business Media, 2012, https://doi.org/10.1007/978-1-4684-9286-6.
  * [5] D. H. Alonso, L. F. N. de Sá, J. S. R. Saenz, and E. C. N. Silva, Topology optimization applied to the design of 2D swirl flow devices, Structural and Multidisciplinary Optimization, 58 (2018), pp. 2341–2364, https://doi.org/10.1007/s00158-018-2078-0.
  * [6] D. H. Alonso, J. S. R. Saenz, and E. C. N. Silva, Non-Newtonian laminar 2D swirl flow design by the topology optimization method, Structural and Multidisciplinary Optimization, (2020), pp. 1–23, https://doi.org/10.1007/s00158-020-02499-2.
  * [7] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, (2001), https://doi.org/10.1137/S0895479899358194.
  * [8] D. N. Arnold, Finite element exterior calculus, SIAM, 2018, https://doi.org/10.1137/1.9781611975543.
  * [9] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis, 39 (2002), pp. 1749–1779, https://doi.org/10.1137/S0036142901384162.
  * [10] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkout, W. Gropp, R. Tran Mills, T. Munson, K. Rupp, P. Sana, B. Smith, S. Zampini, H. Zhang, and H. Zhang, PETSc Users Manual, Tech. Report ANL-95/11 - Revision 3.11, Argonne National Laboratory, 2019, http://www.mcs.anl.gov/petsc.
  * [11] M. P. Bendsøe and O. Sigmund, Topology Optimization, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004, https://doi.org/10.1007/978-3-662-05086-6.
  * [12] S. J. Benson and T. S. Munson, Flexible complementarity solvers for large-scale applications, Optimization Methods and Software, 21 (2003), pp. 155–168, https://doi.org/10.1080/10556780500065382.
  * [13] T. Borrvall and J. Petersson, Topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Fluids, 41 (2003), pp. 77–107, https://doi.org/10.1002/fld.426.
  * [14] F. Brezzi, J. Douglas, R. Durán, and M. Fortin, Mixed finite elements for second order elliptic problems in three variables, Numerische Mathematik, 51 (1987), pp. 237–250, https://doi.org/10.1007/BF01396752.
  * [15] F. Brezzi, J. Douglas, and L. D. Marini, Two families of mixed finite elements for second order elliptic problems, Numerische Mathematik, 47 (1985), pp. 217–235, https://doi.org/10.1007/BF01389710.
  * [16] A. Buffa and C. Ortner, Compact embeddings of broken Sobolev spaces and applications, IMA Journal of Numerical Analysis, 29 (2009), pp. 827–855, https://doi.org/10.1093/imanum/drn038.
  * [17] B. Cockburn, G. Kanschat, and D. Schötzau, A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes equations, Journal of Scientific Computing, 31 (2007), pp. 61–73, https://doi.org/10.1007/s10915-006-9107-7.
  * [18] B. Cockburn, G. Kanschat, D. Schötzau, and C. Schwab, Local discontinuous Galerkin methods for the Stokes system, SIAM Journal on Numerical Analysis, 40 (2002), pp. 319–343, https://doi.org/10.1137/S0036142900380121.
  * [19] Y. Deng, Z. Liu, J. Wu, and Y. Wu, Topology optimization of steady Navier–Stokes flow with body force, Computer Methods in Applied Mechanics and Engineering, 255 (2013), pp. 306–321, https://doi.org/10.1016/j.cma.2012.11.015.
  * [20] Y. Deng, Z. Liu, P. Zhang, Y. Liu, and Y. Wu, Topology optimization of unsteady incompressible Navier–Stokes flows, Journal of Computational Physics, 230 (2011), pp. 6688–6708, https://doi.org/10.1016/j.jcp.2011.05.004.
  * [21] L. C. Evans, Partial Differential Equations, American Mathematical Society, 2 ed., 2010.
  * [22] A. Evgrafov, Topology optimization of slightly compressible fluids, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics, 86 (2006), pp. 46–62, https://doi.org/10.1002/zamm.200410223.
  * [23] A. Evgrafov, State space Newton’s method for topology optimization, Computer Methods in Applied Mechanics and Engineering, 278 (2014), pp. 272–290, https://doi.org/10.1016/j.cma.2014.06.005.
  * [24] P. E. Farrell, Á. Birkisson, and S. W. Funke, Deflation techniques for finding distinct solutions of nonlinear partial differential equations, SIAM Journal on Scientific Computing, 37 (2015), pp. A2026–A2045, https://doi.org/10.1137/140984798.
  * [25] P. E. Farrell, M. Croci, and T. M. Surowiec, Deflation for semismooth equations, Optimization Methods and Software, (2019), pp. 1–24, https://doi.org/10.1080/10556788.2019.1613655.
  * [26] P. E. Farrell, L. Mitchell, L. R. Scott, and F. Wechsung, A Reynolds-robust preconditioner for the Scott-Vogelius discretization of the stationary incompressible Navier-Stokes equations, The SMAI Journal of Computational Mathematics, 7 (2021), pp. 75–96, https://doi.org/10.5802/smai-jcm.72.
  * [27] P. E. Farrell, L. Mitchell, and F. Wechsung, An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier–Stokes equations at high Reynolds number, SIAM Journal on Scientific Computing, 41 (2019), pp. A3073–A3096, https://doi.org/10.1137/18M1219370.
  * [28] P. Fitzpatrick, Advanced calculus, vol. 5, American Mathematical Soc., 2 ed., 2009.
  * [29] I. Fonseca and G. Leoni, Modern Methods in the Calculus of Variations: Lp Spaces, Springer Monographs in Mathematics, Springer New York, New York, NY, 2006, https://doi.org/10.1007/978-0-387-69006-3.
  * [30] N. R. Gauger, A. Linke, and P. W. Schroeder, On high-order pressure-robust space discretisations, their advantages for incompressible high Reynolds number generalised Beltrami flows and beyond, The SMAI Journal of Computational Mathematics, 5 (2019), pp. 89–129, https://doi.org/10.5802/smai-jcm.44.
  * [31] A. Gersborg-Hansen, O. Sigmund, and R. B. Haber, Topology optimization of channel flow problems, Structural and Multidisciplinary Optimization, 30 (2005), pp. 181–192, https://doi.org/10.1007/s00158-004-0508-7.
  * [32] Q. Hong, J. Kraus, J. Xu, and L. Zikatanov, A robust multigrid method for discontinuous Galerkin discretizations of Stokes and linear elasticity equations, Numerische Mathematik, 132 (2016), pp. 23–49, https://doi.org/10.1007/s00211-015-0712-y.
  * [33] V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM review, 59 (2017), pp. 492–544, https://doi.org/10.1137/15M1047696.
  * [34] J. Könnö and R. Stenberg, $H(\mathrm{div})$-conforming finite elements for the Brinkman problem, Mathematical Models and Methods in Applied Sciences, 21 (2011), pp. 2227–2248, https://doi.org/10.1142/S0218202511005726.
  * [35] J. Könnö and R. Stenberg, Numerical computations with $H(\mathrm{div})$-finite elements for the Brinkman problem, Computational Geosciences, 16 (2012), pp. 139–158, https://doi.org/10.1007/s10596-011-9259-x.
  * [36] S. Kreissl, G. Pingen, and K. Maute, Topology optimization for unsteady flow, International Journal for Numerical Methods in Engineering, 87 (2011), pp. 1229–1253, https://doi.org/10.1002/nme.3151.
  * [37] A. Linke and L. G. Rebholz, Pressure-induced locking in mixed methods for time-dependent (Navier–) Stokes equations, Journal of Computational Physics, 388 (2019), pp. 350–356, https://doi.org/10.1016/j.jcp.2019.03.010.
  * [38] J.-C. Nédélec, Mixed finite elements in $\mathbb{R}^{3}$, Numerische Mathematik, 35 (1980), pp. 315–341, https://doi.org/10.1007/BF01396415.
  * [39] L. H. Olesen, F. Okkels, and H. Bruus, A high-level programming-language implementation of topology optimization applied to steady-state Navier–Stokes flow, International Journal for Numerical Methods in Engineering, 65 (2006), pp. 975–1001, https://doi.org/10.1002/nme.1468.
  * [40] I. P. A. Papadopoulos and P. E. Farrell, Preconditioners for computing multiple solutions in three-dimensional fluid topology optimization, arXiv preprint arXiv:2202.08248, (2022).
  * [41] I. P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec, Computing multiple solutions of topology optimization problems, SIAM Journal on Scientific Computing, 43 (2021), pp. A1555–A1582, https://doi.org/10.1137/20M1326209.
  * [42] I. P. A. Papadopoulos and E. Süli, Numerical analysis of a topology optimization problem for Stokes flow, arXiv preprint arXiv:2102.10408, (2021).
  * [43] J. Petersson, A finite element analysis of optimal variable thickness sheets, SIAM Journal on Numerical Analysis, 36 (1999), pp. 1759–1778, https://doi.org/10.1137/S0036142996313968.
  * [44] J. Qin, On the convergence of some low order mixed finite elements for incompressible fluids, PhD thesis, Pennsylvania State University, 1994.
  * [45] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. McRae, G.-T. Bercea, G. R. Markall, and P. H. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions on Mathematical Software (TOMS), 43 (2016), pp. 1–27, https://doi.org/10.1145/2998441.
  * [46] P.-A. Raviart and J.-M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, Springer, 1977, pp. 292–315, https://doi.org/10.1007/BFb0064470.
  * [47] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM Journal on Scientific Computing, 14 (1993), pp. 461–469, https://doi.org/10.1137/0914028.
  * [48] J. Schöberl, Multigrid methods for a parameter dependent problem in primal variables, Numerische Mathematik, 84 (1999), pp. 97–119, https://doi.org/10.1007/s002110050465.
  * [49] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, ESAIM: Mathematical Modelling and Numerical Analysis, 19 (1985), pp. 111–143, https://doi.org/10.1051/m2an/1985190101111.
  * [50] C.-J. Thore, Topology optimization of Stokes flow with traction boundary conditions using low-order finite elements, Computer Methods in Applied Mechanics and Engineering, 386 (2021), p. 114177, https://doi.org/10.1016/j.cma.2021.114177.
  * [51] S. Zhang, Divergence-free finite elements on tetrahedral grids for $k\geq 6$, Mathematics of Computation, 80 (2011), pp. 669–695, https://doi.org/10.1090/S0025-5718-2010-02412-3.