# Interior penalty discontinuous Galerkin methods for the nearly
incompressible elasticity eigenvalue problem with heterogeneous media

Arbaz Khan  Department of Mathematics, Indian Institute of Technology Roorkee,
Roorkee 247667, India<EMAIL_ADDRESS>,  Felipe Lepe  GIMNAP-Departamento
de Matemática, Universidad del Bío - Bío, Casilla 5-C, Concepción, Chile
<EMAIL_ADDRESS>and  Jesus Vellojin  GIMNAP-Departamento de Matemática,
Universidad del Bío - Bío, Casilla 5-C, Concepción, Chile
<EMAIL_ADDRESS>
###### Abstract.

This paper studies the family of interior penalty discontinuous Galerkin
methods for solving the Herrmann formulation of the linear elasticity
eigenvalue problem in heterogeneous media. By employing a weighted Lamé
coefficient norm within the framework of non-compact operators theory, we
prove convergence of both continuous and discrete eigenvalue problems as the
mesh size approaches zero, independently of the Lamé constants. Additionally,
we conduct an a posteriori analysis and propose a reliable and efficient
estimator. The theoretical findings are supported by numerical experiments.

###### Key words and phrases:

Eigenvalue problems, Discontinuous Galerkin method, Error estimates, A
posteriori analysis

###### 2020 Mathematics Subject Classification:

35Q35, 65N15, 65N25, 65N30, 76D07

The first author was partially supported by the Sponsored Research &
Industrial Consultancy (SRIC), Indian Institute of Technology Roorkee, India
through the faculty initiation grant MTD/FIG/100878; by SERB MATRICS grant
MTR/2020/000303; by SERB Core research grant CRG/2021/002569.

The second author was partially supported by DIUBB through project 2120173
GI/C Universidad del Bío-Bío.

The third author was partially supported by the National Agency for Research
and Development, ANID-Chile through FONDECYT Postdoctorado project 3230302

##  1\. Introduction

Recently in [15] the nearly incompressible elasticity eigenvalue problem has
been analyzed. In this reference the authors have proved that the elasticity
spectrum converges to the Stokes eigenvalues when the Poisson ratio tends to
$1/2$. From this fundamental observation, and with the aid of inf-sup stable
families of finite elements for Stokes, the authors approximated the spectrum
accurately. Motivated by these findings and in alignment with our research
trajectory, we now turn our attention to the application of the interior
penalty discontinuous Galerkin method.

Interior Penalty Discontinuous Galerkin Finite Element Method (IPDG)
represents a sophisticated approach to tackling eigenvalue problems,
leveraging the advantages of both discontinuous Galerkin methods and interior
penalty techniques. This method excels in handling problems associated with
singularities, discontinuities, or complex geometries by incorporating penalty
terms within the weak formulation. The inclusion of these penalty terms
facilitates the imposition of continuity constraints across element
boundaries, effectively minimizing numerical oscillations and ensuring
stability. IPDG proves particularly adept at capturing eigenvalues in problems
with variable coefficients or irregular domains. By introducing penalty terms
to penalize jumps in the solution and its gradient across element interfaces,
this method strikes a balance between accuracy and stability, making it a
reliable choice for eigenvalue problems characterized by challenging features
and intricate geometries.

The application of the IPDG methods for eigenvalue problems is a current
subject of study and we refer to [2, 4, 17, 19, 18, 20] as our main references
which explore applications encompassing the Laplace operator, the Maxwell’s
equations, Stokes, elasticity and acoustics. These studies collectively
demonstrate the broad applicability and efficacy of IPDG techniques in diverse
scientific domains. Here the main advantage lies, in one hand, on the
flexibility of the method where hanging nodes are allowed, and on the other,
the easy computational implementation where with no major difficulties is
possible to consider high order approximations for the method. This inherent
adaptability not only enables accurate approximation of the spectrum but also
highlights the method’s suitability for tackling complex eigenvalue problems.
However, there is a price to pay when using this method and is related to the
stabilization parameter that precisely penalizes the jumps between elements.
While this parameter is crucial for ensuring stability and convergence, its
determination can pose significant computational hurdles. Theoretical
investigations [17, 19, 18, 20] have underscored the significance of the
stabilization parameter in determining the stability of IPDG methods. The
appropriate selection of this parameter depends on several factors, such as
boundary conditions, geometrical assumptions and physical parameters, among
others. The correct choice of the stabilization parameter is fundamental to
achieve an accurate approximation to the physical spectrum, without spurious
eigenvalues. The literature suggests that the stabilization parameters leading
to accurate results are proportional to the square of the polynomial degree,
although the proportionality constant varies depending on the specific
physical parameters involved. This empirical observation highlights the
intricate interplay between numerical methods and physical phenomena, and
emphasizes the need for judicious parameter selection to ensure the fidelity
of the computational results.

###  1.1. Notations

Throughout this work, $\Omega$ is a generic Lipschitz bounded domain of
$\mathbb{R}^{d}$, $d\in\\{2,3\\}$. For $s\geq 0$,
$\lVert\cdot\rVert_{s,\Omega}$ stands indistinctly for the norm of the
Hilbertian Sobolev spaces $\mathrm{H}^{s}(\Omega)$ or
$[\mathrm{H}^{s}(\Omega)]^{2}$ with the convention
$\mathrm{H}^{0}(\Omega):=\mathrm{L}^{2}(\Omega)$. If ${X}$ and ${Y}$ are
normed vector spaces, we write ${X}\hookrightarrow{Y}$ to denote that ${X}$ is
continuously embedded in ${Y}$. We denote by ${X}^{\prime}$ and
$\|\cdot\|_{{X}}$ the dual and the norm of ${X}$, respectively. Finally, we
employ $\mathbf{0}$ to denote a generic null vector and the relation
$\texttt{a}\lesssim\texttt{b}$ indicates that $\texttt{a}\leq C\texttt{b}$,
with a positive constant $C$ which is independent of a, b, and the size of the
elements in the mesh. The value of $C$ might change at each occurrence. We
remark that we will write the constant $C$ and the estimates dependencies only
when is needed.

Let us recall the following standard notations for any tensor field
$\boldsymbol{\tau}=\left(\tau_{ij}\right)_{1\leq i,j\leq d}$, any vector field
$\boldsymbol{v}=\left(v_{i}\right)_{1\leq i\leq d}$ and any scalar field $v$:

$\begin{gathered}\nabla\cdot\boldsymbol{v}:=\sum_{i=1}^{d}\partial_{i}v_{i},\quad\nabla
v:=\left(\partial_{i}v\right)_{1\leq i\leq
d},\quad\boldsymbol{\nabla}\cdot\boldsymbol{\tau}:=\left(\sum_{j=1}^{d}\partial_{j}\tau_{ij}\right)_{1\leq
i\leq d},\\\ \nabla\boldsymbol{v}:=\left(\partial_{i}v_{j}\right)_{1\leq
i,j\leq
d},\quad\boldsymbol{\tau}^{\mathrm{t}}:=\left(\tau_{ji}\right),\quad\operatorname{tr}(\boldsymbol{\tau}):=\sum_{i=1}^{d}\tau_{ii},\quad\boldsymbol{\tau}^{a}:\boldsymbol{\tau}^{b}:=\sum_{i,j=1}^{d}\tau_{ij}^{a}\tau_{ij}^{b}.\end{gathered}$

###  1.2. Model problem

Let $\Omega\subset\mathbb{R}^{d}$ with $d\in\\{2,3\\}$ be an open and bounded
domain with Lipschitz boundary $\partial\Omega$. Let us assume that the
boundary is divided into two parts $\partial\Omega:=\Gamma_{D}\cup\Gamma_{N}$,
where $|\Gamma_{D}|>0$. In particular, we assume that the structure is fixed
on $\Gamma_{D}$ and free of stress on $\Gamma_{N}$.

The eigenvalue problem of interest is given by: Find
$(\widehat{\kappa},\boldsymbol{u})$ such that

(1.1)
$\displaystyle-\boldsymbol{\nabla}\cdot\left[2\widehat{\mu}\boldsymbol{\varepsilon}(\boldsymbol{u})+\widehat{\lambda}\mathop{\mathrm{tr}}\nolimits(\boldsymbol{\varepsilon}(\boldsymbol{u}))\boldsymbol{I}\right]$
$\displaystyle=\widehat{\kappa}\rho\boldsymbol{u}$ $\displaystyle\text{ in
}\Omega,$ $\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{0}$
$\displaystyle\mbox{ on }\Gamma_{D},$
$\displaystyle\left[2\widehat{\mu}\boldsymbol{\varepsilon}(\boldsymbol{u})+\widehat{\lambda}\mathop{\mathrm{tr}}\nolimits(\boldsymbol{\varepsilon}(\boldsymbol{u}))\boldsymbol{I}\right]\boldsymbol{n}$
$\displaystyle={\boldsymbol{0}}$ $\displaystyle\mbox{ on }\Gamma_{N},$

where $\boldsymbol{I}\in\mathbb{R}^{d\times d}$ represents the identity
matrix, $\sqrt{\widehat{\kappa}}$ is the natural frequency, $\boldsymbol{u}$
is the displacement and $\rho$ is the material density. The Lamé parameters
$\widehat{\lambda}$ and $\widehat{\mu}$ are defined by

$\widehat{\lambda}:=\frac{E(x)\nu}{(1+\nu)(1-2\nu)}\quad\text{and}\quad\widehat{\mu}:=\frac{E(x)}{2(1+\nu)},$

where $E(x)$ represents the space-variable Young’s modulus, $\nu$ represents
the Poisson ratio which we assume as a constant. The strain tensor is
represented by $\boldsymbol{\varepsilon}(\cdot)$ and defined as
$\boldsymbol{\varepsilon}(\boldsymbol{v}):=\frac{1}{2}(\nabla\boldsymbol{v}+(\nabla\boldsymbol{v})^{\texttt{t}})$,
where t represents the transpose operator. Here and thereafter we assume the
existence of two positive constants $E_{\min}$ and $E_{\max}$ such that
$E_{\min}\leq E(x)\leq E_{\max}.$

For the analysis of (1.1) we scale the system with the quantity $(1+\nu)$
leading to the following problem. Find $(\kappa,\boldsymbol{u})$ with
$\boldsymbol{u}\neq\boldsymbol{0}$ such that

(1.2)
$\displaystyle-\boldsymbol{\nabla}\cdot\left[2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})+\lambda\mathop{\mathrm{tr}}\nolimits(\boldsymbol{\varepsilon}(\boldsymbol{u}))\boldsymbol{I}\right]$
$\displaystyle=\kappa\rho\boldsymbol{u}$ $\displaystyle\text{ in }\Omega,$
$\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{0}$
$\displaystyle\mbox{ on }\Gamma_{D},$
$\displaystyle\left[2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})+\lambda\mathop{\mathrm{tr}}\nolimits(\boldsymbol{\varepsilon}(\boldsymbol{u}))\boldsymbol{I}\right]\boldsymbol{n}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ on }\Gamma_{N},$

where the scaled eigenvalue and the Lamé coefficients are given by

(1.3)
$\kappa:=(1+\nu)\widehat{\kappa},\qquad\lambda:=\frac{E(x)\nu}{1-2\nu},\qquad\mu:=\frac{E(x)}{2}.$

Note that for a given $E(x)$, the eigenvalue problem (1.2) has the advantage
of having $\mu$ independent of $\nu$. The incompressibility of the material is
explicitly given only by the behavior of the constant $\nu$ through the Lamé
parameter $\lambda$.

With the scaled parameters specified in (1.3) we rewrite system (1.1) as the
following system: We seek the eigenvalue $\kappa\in\mathbb{R}^{+}$, a vector
field $\boldsymbol{u}\neq\boldsymbol{0}$, and a scalar field $p\neq 0$ such
that

$\displaystyle-\boldsymbol{\nabla}\cdot(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u}))+\nabla
p$ $\displaystyle=\kappa\rho\boldsymbol{u}$ $\displaystyle\mbox{ in }\Omega,$
$\displaystyle\nabla\cdot\boldsymbol{u}+\frac{1}{\lambda}p$ $\displaystyle=0$
$\displaystyle\mbox{ in }\Omega,$ (1.4) $\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ on }\Gamma_{D},$
$\displaystyle\left[2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})-p\boldsymbol{I}\right]\boldsymbol{n}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\mbox{ on }\Gamma_{N}.$

We note from (1.3) that if the Poisson ratio approaches $1/2$, then
$\lambda\uparrow\infty$. In such case, we recover the Stokes system. This is
an important matter in view of deriving robust stability and error bounds.
From now and on, and for simplicity, in the analysis below we consider
$\rho=1$.

###  1.3. Main contribution

To the best of our knowledge, the realm of IPDG methods remains uncharted
territory in the context of the model problem (1.2). This paper marks the
inaugural exploration in this direction. The primary focus is on the family of
IPDG methods employed for solving the eigenvalue problem in nearly
incompressible elasticity within heterogeneous media. A cornerstone of this
research lies in the comprehensive examination of both a priori and a
posteriori error analyses, shedding light on the method’s performance.
Notably, our proposed estimator exhibits remarkable robustness with respect to
the Poisson ratio, underpinning its suitability for addressing complex, real-
world scenarios. Furthermore, it is established that the spectrum of the
nearly incompressible elasticity eigenvalue problem within heterogeneous media
converges to the spectrum of the Stokes problems as the parameter $\lambda$
tends towards infinity

###  1.4. Outline

The subsequent sections of the paper are structured as follows: Section 2
delves into the weak formulation and presents some initial findings. Section 3
explores the family of IPDG methods concerning the model problem (1.2).
Convergence results and error estimates are thoroughly examined in Section 4.
A comprehensive a posteriori error analysis is detailed in Section 5. Lastly,
Section 6 showcases numerical results, demonstrating the effectiveness of the
IPDG methods.

##  2\. Variational formulation and preliminary results

Let us begin by defining the spaces
$\mathrm{H}:=\left\\{\boldsymbol{v}\in\mathrm{H}^{1}(\Omega)^{d}\;:\;\boldsymbol{v}=\boldsymbol{0}\text{
on }\Gamma_{D}\right\\}$ and $\mathrm{Q}:=\mathrm{L}^{2}(\Omega)$ where we
will seek the displacement and the pressure, respectively. Now, by testing
system (1.2) with adequate functions and imposing the boundary conditions, we
end up with the following saddle point variational formulation: Find
$(\kappa,(\boldsymbol{u},p))\in\mathbb{R}^{+}\times\mathrm{H}\times\mathrm{Q}$
with $(\boldsymbol{u},p)\neq(\boldsymbol{0},0)$ such that

$\displaystyle
2\int_{\Omega}\mu(x)\boldsymbol{\varepsilon}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v})\,dx-\int_{\Omega}p\nabla\cdot\boldsymbol{v}\,dx$
$\displaystyle=\kappa\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,dx$
$\displaystyle\quad\forall\boldsymbol{v}\in\mathrm{H},$
$\displaystyle-\int_{\Omega}q\nabla\cdot\boldsymbol{u}\,dx-\int_{\Omega}\frac{1}{\lambda(x)}pq\,dx$
$\displaystyle=\,0$ $\displaystyle\quad\forall q\in\mathrm{Q}.$

This variational problem can be rewritten as follows: Find
$(\kappa,(\boldsymbol{u},p))\in\mathbb{R}\times\mathrm{H}\times\mathrm{Q}$
such that

(2.1)
$\left\\{\begin{array}[]{rcll}a(\boldsymbol{u},\boldsymbol{v})+b(\boldsymbol{v},p)&=&\kappa
d(\boldsymbol{u},\boldsymbol{v})&\forall\boldsymbol{v}\in\mathrm{H},\\\
b(\boldsymbol{u},q)-c(p,q)&=&\;0&\forall q\in\mathrm{Q},\end{array}\right.$

where the bilinear forms $a:\mathrm{H}\times\mathrm{H}\to\mathbb{R}$,
$b:\mathrm{H}\times\mathrm{Q}\to\mathbb{R}$,
$c:\mathrm{Q}\times\mathrm{Q}\to\mathbb{R}$, and
$d:\mathrm{H}\times\mathrm{H}\to\mathbb{R}$ are defined by

$\displaystyle a(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle:=2\int_{\Omega}\mu(x)\boldsymbol{\varepsilon}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v})\,dx,\quad
b(\boldsymbol{v},q):=-\int_{\Omega}q\nabla\cdot\boldsymbol{v}\,dx,$
$\displaystyle d(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle:=\int_{\Omega}\boldsymbol{u}\cdot\boldsymbol{v}\,dx,\quad
c(p,q):=\int_{\Omega}\frac{1}{\lambda(x)}p\,q\,dx,$

for all $\boldsymbol{u},\boldsymbol{v}\in\mathrm{H}$, and $p,q\in\mathrm{Q}$.

To perform the analysis, we need a suitable norm which in particular depends
on the parameters of the problem. With this in mind, and for all
$(\boldsymbol{v},q)\in\mathrm{H}\times Q$, we define the following weighted
norm,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}\times
Q}^{2}:=\|\mu(x)^{1/2}\nabla\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mu(x)^{-1/2}q\|_{0,\Omega}^{2}+\|\lambda(x)^{-1/2}q\|_{0,\Omega}^{2}.$

Moreover, in what follows, we assume that
$\mu(x),1/\lambda(x)\in\mathrm{W}^{1,\infty}(\Omega)$.

Let us introduce and define the kernel of $b(\cdot,\cdot)$ by

$\mathcal{K}:=\\{\boldsymbol{v}\in\mathrm{H}\,:\,b(\boldsymbol{v},q)=0,\quad\forall\,q\in
Q\\}\,=\,\\{\boldsymbol{v}\in\mathrm{H}\,:\,\nabla\cdot\boldsymbol{v}\,=\,0\,\,\,\,{\rm
in}\,\,\Omega\\},$

and let us recall that the bilinear form $b(\cdot,\cdot)$ satisfies the inf-
sup condition:

$\sup_{\stackrel{{\scriptstyle\scriptstyle\boldsymbol{v}\in\mathrm{H}}}{{\boldsymbol{v}\neq
0}}}\frac{|b(\boldsymbol{v},q)|}{\|\mu(x)^{1/2}\nabla\boldsymbol{v}\|_{0,\Omega}}\geq\beta_{2}\|\mu(x)^{-1/2}q\|_{0,\Omega}\quad\forall
q\in\mathrm{Q},$

with an inf-sup constant $\beta_{2}>0$ only depending on $\Omega$; see [16],
for instance.

Now, we are in a position to introduce the solution operator which we denote
by
$\boldsymbol{T}_{\nu}:\mathrm{L}^{2}(\Omega)^{d}\rightarrow\mathrm{L}^{2}(\Omega)^{d}$
and is such that
$\boldsymbol{f}\mapsto\boldsymbol{T}_{\nu}\boldsymbol{f}:=\widehat{\boldsymbol{u}}$,
where the pair
$(\widehat{\boldsymbol{u}},\widehat{p})\in\mathrm{H}\times\mathrm{Q}$ is the
solution of the following well posed source problem: Given
$\boldsymbol{f}\in\mathrm{L}^{2}(\Omega)^{d}$ find
$(\widehat{\boldsymbol{u}},\widehat{p})\in\mathrm{H}\times\mathrm{Q}$ such
that

(2.2)
$\left\\{\begin{array}[]{rcll}a(\widehat{\boldsymbol{u}},\boldsymbol{v})+b(\boldsymbol{v},\widehat{p})&=&d(\boldsymbol{f},\boldsymbol{v})&\forall\boldsymbol{v}\in\mathrm{H},\\\
b(\widehat{\boldsymbol{u}},q)-c(\widehat{p},q)&=&\;0&\forall
q\in\mathrm{Q}.\end{array}\right.$

Thanks to the Babuška-Brezzi theory we have that $\boldsymbol{T}_{\nu}$ is
well defined and the following estimate holds
$\|\boldsymbol{T}_{\nu}\boldsymbol{f}\|_{0,\Omega}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{T}_{\nu}\boldsymbol{f},\widehat{p})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\lesssim\|\boldsymbol{f}\|_{0,\Omega}$
with a constant depending on the inverse of $\mu$ and the Poincaré constant.
It is easy to check that $\boldsymbol{T}_{\nu}$ is a selfadjoint operator with
respect to the $\mathrm{L}^{2}$ inner product. Also, let $\chi$ be a real
number such that $\chi\neq 0$. Notice that
$(\chi,\boldsymbol{u})\in\mathbb{R}\times\mathrm{H}$ is an eigenpair of
$\boldsymbol{T}_{\nu}$ if and only if there exists $p\in\mathrm{Q}$ such that,
$(\kappa,(\boldsymbol{u},p))$ solves problem (2.1) with $\chi:=1/\kappa$.

The following result states an additional regularity result for problem (2.2).

######  Lemma 2.1.

Let $(\widehat{\boldsymbol{u}},\widehat{p})$ be the unique solution of (2.2),
then there exists $s\in(0,1]$ such that
$(\widehat{\boldsymbol{u}},\widehat{p})\in\mathrm{H}^{1+s}(\Omega)^{d}\times\mathrm{H}^{s}(\Omega)$
and the following estimate holds (see for instance [11] or [21])

$\|\widehat{\boldsymbol{u}}\|_{1+s,\Omega}+\|\widehat{p}\|_{s,\Omega}\lesssim\|\boldsymbol{f}\|_{0,\Omega},$

where the hidden constant depends on $\Omega$ and the Lamé coefficients.

Let us remark that Lemma 2.1 holds true even for $\lambda\uparrow\infty$
(Further comments on this fact can be seen in [18]).

As a consequence of this result, we conclude the compactness of
$\boldsymbol{T}_{\nu}$. In fact, since
$\mathrm{H}^{1+s}(\Omega)^{d}\hookrightarrow\mathrm{L}^{2}(\Omega)^{d}$ we
have that $\boldsymbol{T}_{\nu}$ is a compact operator. Hence, we are in
position to establish the spectral characterization of $\boldsymbol{T}_{\nu}$.

######  Theorem 2.2.

The spectrum of $\boldsymbol{T}_{\nu}$ satisfies
$\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T}_{\nu})=\\{0\\}\cup\\{\chi_{k}\\}_{k\in\mathbb{N}}$,
where $\\{\chi_{k}\\}_{k\in\mathbb{N}}$ is a sequence of positive eigenvalues.

On the other hand, from [8] we have the following regularity result for the
eigenfunctions of the eigenproblem (2.1).

######  Theorem 2.3.

If $(\kappa,(\boldsymbol{u},p))\in\mathbb{R}\times\mathrm{H}\times\mathrm{Q}$
solves (2.1), then, there exists $r>0$ such that
$\boldsymbol{u}\in\mathrm{H}^{1+r}(\Omega)^{d}$ and
$p\in\mathrm{H}^{r}(\Omega)$. Moreover, the following estimate holds

$\|\boldsymbol{u}\|_{1+r,\Omega}+\|p\|_{r,\Omega}\lesssim\|\boldsymbol{u}\|_{0,\Omega}.$

Now, Let us end this section introducing the bilinear form
$\mathcal{A}:(\mathrm{H}\times\mathrm{Q})\times(\mathrm{H}\times\mathrm{Q})\rightarrow\mathbb{R}$
defined by

$\mathcal{A}((\boldsymbol{u},p),(\boldsymbol{v},q)):=a(\boldsymbol{u},\boldsymbol{v})+b(\boldsymbol{v},p)+b(\boldsymbol{u},q)+c(p,q),$

which allows us to rewrite the eigenvalue problem as follows. Find
$\kappa\in\mathbb{R}$ and
$(\boldsymbol{0},0)\neq(\boldsymbol{u},p)\in\mathrm{H}\times\mathrm{Q}$ such
that

(2.3) $\mathcal{A}((\boldsymbol{u},p),(\boldsymbol{v},q))=\kappa
d(\boldsymbol{u},\boldsymbol{v})\quad\forall(\boldsymbol{v},q)\in\mathrm{H}\times
Q.$

Let us remark that $\mathcal{A}(\cdot,\cdot)$ is a bounded bilinear form, in
the sense that there exists a positive constant
$\widehat{C}:=\max\\{2,\mu_{min}^{-1/2}\mu_{max}^{1/2}\\}$ such that

$\mathcal{A}((\boldsymbol{u},p),(\boldsymbol{v},q))\leq\widehat{C}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u},p)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}\times
Q}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}\times
Q}\quad\forall(\boldsymbol{u},p),(\boldsymbol{v},q)\in\mathrm{H}\times\mathrm{Q}.$

Next we state the stability result for the bilinear form
$\mathcal{A}(\cdot,\cdot)$.

######  Lemma 2.4.

For any $(\boldsymbol{u},p)\in\mathrm{H}\times\mathrm{Q}$, there exists
$(\boldsymbol{v},q)\in\mathrm{H}\times\mathrm{Q}$ with
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u},p)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
such that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u},p)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}\lesssim\mathcal{A}((\boldsymbol{u},p),(\boldsymbol{v},q)).$

##  3\. The DG method

Now our aim is to introduce the DG methods. In order to do this, we need to
set some notations and definitions, inherent for these type of methods, as DG
spaces, jumps, averages, and discrete bilinear forms. For the sake of
completeness, in this section we will rigorously monitor the constants in
order to provide bounds that explicitly reflect the importance of the physical
parameters in the stability of the method.

###  3.1. Preliminaries

Let $\mathcal{T}_{h}$ be a shape regular family of meshes which subdivide the
domain $\bar{\Omega}$ into triangles/tetrahedra that we denote by $K$. Let us
denote by $h_{K}$ the diameter of any element $K\in\mathcal{T}_{h}$ and let
$h$ be the maximum of the diameters of all the elements of the mesh, i.e.
$h:=\max_{K\in\mathcal{T}_{h}}\\{h_{K}\\}$.

Let $F$ be a closed set. We say that $F\subset\overline{\Omega}$ is an
interior edge/face if $F$ has a positive $(n-1)$-dimensional measure and if
there are distinct elements $K$ and $K^{\prime}$ such that
$F=\bar{K}\cap\bar{K}^{\prime}$. A closed subset $F\subset\overline{\Omega}$
is a boundary edge/face if there exists $K\in\mathcal{T}_{h}$ such that $F$ is
an edge/face of $K$ and $F=\bar{K}\cap\Gamma$. Let $\mathcal{F}_{h}^{0}$ and
$\mathcal{F}_{h}^{\partial}$ be the sets of interior edges/faces and boundary
edges/faces, respectively. We assume that the boundary mesh
$\mathcal{F}_{h}^{\partial}$ is compatible with the partition
$\Gamma=\Gamma_{D}\cup\Gamma_{N}$, namely,

$\bigcup_{F\in\mathcal{F}_{h}^{D}}F=\Gamma_{D}\qquad\text{and}\qquad\bigcup_{F\in\mathcal{F}_{h}^{N}}F=\Gamma_{N},$

where $\mathcal{F}_{h}^{D}:=\\{F\in\mathcal{F}_{h}^{\partial};\quad
F\subset\Gamma_{D}\\}$ and
$\mathcal{F}_{h}^{N}:=\\{F\in\mathcal{F}_{h}^{\partial}:\quad
F\subset\Gamma_{N}\\}$. Also we denote
$\mathcal{F}_{h}:=\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{\partial}$ and
$\mathcal{F}^{*}_{h}:=\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{D}$. Also, for
any element $K\in\mathcal{T}_{h}$, we introduce the set
$\mathcal{F}(K):=\\{F\in\mathcal{F}_{h}:\,\,F\subset\partial K\\}$ of
edges/faces composing the boundary of $K$.

For any $t\geq 0$, we define the following broken Sobolev space

$\mathrm{H}^{t}(\mathcal{T}_{h})^{n}:=\\{\boldsymbol{v}\in[\mathrm{L}^{2}(\Omega)]^{n}:\quad\boldsymbol{v}|_{K}\in[\mathrm{H}^{t}(K)]^{n}\quad\forall
K\in\mathcal{T}_{h}\\}.$

Also, the space of the skeletons of the triangulations $\mathcal{T}_{h}$ is
defined by
$\mathrm{L}^{2}(\mathcal{F}_{h}):=\prod_{F\in\mathcal{F}_{h}}\mathrm{L}^{2}(F).$

In the forthcoming analysis,
$h_{\mathcal{F}}\in\mathrm{L}^{2}(\mathcal{F}_{h})$ will represent the
piecewise constant function defined by $h_{\mathcal{F}}|_{F}:=h_{F}$ for all
$F\in\mathcal{F}_{h}$, where $h_{F}$ denotes the diameter of an edge/face $F$.

Let $\mathcal{P}_{m}(\mathcal{T}_{h})$ be the space of piecewise polynomials
respect with to $\mathcal{T}_{h}$ of degree at most $m\geq 0$; namely,

$\mathcal{P}_{m}(\mathcal{T}_{h}):=\left\\{\boldsymbol{v}\in[\mathrm{L}^{2}(\Omega)]^{d}:\,v|_{K}\in[\mathcal{P}_{m}(K)]^{d},\forall
K\in\mathcal{T}_{h}\right\\}.$

Let $k\geq 1$. To approximate the displacement we introduce the space

$\mathrm{H}_{h}:=\\{\boldsymbol{v}_{h}\in[\mathrm{L}^{2}(\Omega)]^{d}\,:\quad\boldsymbol{v}_{h}|_{K}\in[\mathcal{P}_{k}(K)]^{d},\,\,\,\forall
K\in{\mathcal{T}}_{h}\\},$

and for the pressure, the finite element space is
$\mathrm{Q}_{h}:=\\{\boldsymbol{v}_{h}\in\mathrm{L}^{2}(\Omega)\,:\,\boldsymbol{v}_{h}|_{K}\in\mathcal{P}_{k-1}(K),\,\,\,\forall
K\in{\mathcal{T}}_{h}\\}$.

Let us assume $\mathcal{F}(K)\cap\mathcal{F}(K^{\prime})\neq\emptyset$. For a
scalar field $q$, we define the average
$\\{q\\}\in\mathrm{L}^{2}(\mathcal{F}_{h})$ and the jump $\llbracket
q\rrbracket\in\mathrm{L}^{2}(\mathcal{F}_{h})$ by
$\\{q\\}:=(q_{K}+q_{K^{\prime}})/2$ and $\llbracket
q\rrbracket:=q_{K}\boldsymbol{n}_{K}+q_{K^{\prime}}\boldsymbol{n}_{K^{\prime}}$,
respectively, where $\boldsymbol{n}_{K}$ is the outward unit normal vector to
$\partial K$ and $q_{K}$ represents the restriction $q|_{K}$. Similarly, if
$\boldsymbol{v}$ is a vector field, we define the average
$\\{\boldsymbol{v}\\}\in[\mathrm{L}^{2}(\mathcal{F}_{h})]^{n}$ by
$\\{\boldsymbol{v}\\}:=(\boldsymbol{v}_{K}+\boldsymbol{v}_{K^{\prime}})/2$,
the scalar jump
$\llbracket\underline{\boldsymbol{v}}\rrbracket\in\mathrm{L}^{2}(\mathcal{F}_{h})$
by
$\llbracket\underline{\boldsymbol{v}}\rrbracket:=\boldsymbol{v}_{K}\cdot\boldsymbol{n}_{K}+\boldsymbol{v}_{K^{\prime}}\cdot\boldsymbol{n}_{K^{\prime}}$,
and the tensor or total jump
$\llbracket\boldsymbol{v}\rrbracket\in[\mathrm{L}^{2}(\mathcal{F}_{h})]^{d\times
d}$ as
$\llbracket\boldsymbol{v}\rrbracket:=\boldsymbol{v}_{K}\otimes\boldsymbol{n}_{K}+\boldsymbol{v}_{K^{\prime}}\otimes\boldsymbol{n}_{K^{\prime}}$.
Finally, if $\boldsymbol{\tau}$ is a tensor field, we define the corresponding
average and jump as
$\\{\boldsymbol{\tau}\\}:=(\boldsymbol{\tau}_{K}+\boldsymbol{\tau}_{K^{\prime}})/2\in[\mathrm{L}^{2}(\mathcal{F}_{h})]^{d\times
d}$ and
$\llbracket\boldsymbol{\tau}\boldsymbol{n}\rrbracket:=\boldsymbol{\tau}_{K}\boldsymbol{n}_{K}+\boldsymbol{\tau}_{K^{\prime}}\boldsymbol{n}_{K^{\prime}}\in[\mathrm{L}^{2}(\mathcal{F}_{h})]^{d}$,
respectively. If we have an element $K$ such that a facet $\mathcal{F}$
satisfies $\mathcal{F}\in\mathcal{F}_{h}^{\partial}$, we can obtain the
definition of average and jump in the domain boundary by taking $K=K^{\prime}$
and $K^{\prime}=0$ in the above definitions, respectively.

Inspired by the analysis of [2], let us define
$\mathrm{H}(h):=\mathrm{H}+\mathrm{H}_{h}$ which we endow with the following
norm

$\|\boldsymbol{v}\|_{\mathrm{H}(h)}^{2}=\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}\rrbracket\|_{0,\mathcal{F}_{h}}^{2},$

which coincides with the natural norm of $\mathrm{H}$. On the other hand, we
also introduce the following norm for $\mathrm{Q}_{h}$

$\|q_{h}\|_{\mathrm{Q}}^{2}:=\|\mu^{-1/2}(x)q_{h}\|_{0,\Omega}^{2}+\|\lambda^{-1/2}(x)q_{h}\|_{0,\Omega}^{2}\quad\forall
q_{h}\in\mathrm{Q}_{h},$

where the natural product norm for the discrete spaces is defined by
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times\mathrm{Q}}^{2}:=\|\boldsymbol{v}_{h}\|_{\mathrm{H}(h)}^{2}+\|q\|_{\mathrm{Q}}^{2}$.
Finally we introduce the following trace inequality [20]

(3.1) $\|h_{\mathcal{F}}^{1/2}\\{\xi(x)v\\}\|_{0,\mathcal{F}}\leq
C_{\texttt{tr}}\|v\|_{0,\Omega}\qquad\forall
v\in\mathcal{P}_{k}(\mathcal{T}_{h}),$

where $\xi(x)\in\mathrm{L}^{\infty}(\Omega)$ and $C_{\texttt{tr}}>0$ is a
constant depending on $\|\xi(x)\|_{0,\infty}$.

###  3.2. Symmetric and nonsymmetric DG schemes

With the discrete spaces defined previously at hand, we introduce the discrete
counterpart of (2.1) as follows: Find $\kappa_{h}\in\mathbb{C}$ and
$(\boldsymbol{0},0)\neq(\boldsymbol{u}_{h},p_{h})\in\mathrm{H}_{h}\times\mathrm{Q}_{h}$
such that

(3.2)
$\begin{array}[]{rcll}a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})+b_{h}(\boldsymbol{v}_{h},p_{h})&=&\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})&\forall\
\boldsymbol{v}_{h}\in\mathrm{H}_{h},\\\ \displaystyle
b_{h}(\boldsymbol{u}_{h},q_{h})-c(p_{h},q_{h})&=&0&\forall\
q_{h}\in\mathrm{Q}_{h},\end{array}$

where the sesquilinear form
$a_{h}:\mathrm{H}_{h}\times\mathrm{H}_{h}\rightarrow\mathbb{C}$ is defined by

(3.3)
$a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}):=\int_{{\mathcal{T}}_{h}}2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}):\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})+\int_{\mathcal{F}^{*}_{h}}\frac{\texttt{a}_{S}}{h_{\mathcal{F}}}2\mu(x)\llbracket\boldsymbol{u}_{h}\rrbracket\cdot\llbracket\boldsymbol{v}_{h}\rrbracket\\\
-\int_{\mathcal{F}^{*}_{h}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h})\\}\cdot\llbracket\boldsymbol{v}_{h}\rrbracket-\epsilon\int_{\mathcal{F}^{*}_{h}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})\\}\cdot\llbracket\boldsymbol{u}_{h}\rrbracket,$

where
$\boldsymbol{\varepsilon}_{h}:=\frac{1}{2}(\nabla_{h}(\cdot)+(\nabla_{h}(\cdot))^{\texttt{t}})$
and $\texttt{a}_{S}>0$ is a positive constant which we will refer as the
stabilization parameter and $\varepsilon\in\\{-1,0,1\\}$. On the other hand,
we define the bounded bilinear form
$b_{h}:\mathrm{H}_{h}\times\mathrm{Q}_{h}\rightarrow\mathbb{R}$ by

$b_{h}(\boldsymbol{v}_{h},q_{h}):=-\int_{{\mathcal{T}}_{h}}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}_{h}q_{h}+\int_{\mathcal{F}_{h}^{*}}\\{q_{h}\\}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket,\quad\forall\boldsymbol{v}\in\mathrm{H}_{h},\,\forall
q_{h}\in Q_{h},$

and $c:\mathrm{Q}_{h}\times\mathrm{Q}_{h}\rightarrow\mathbb{R}$ is defined by
$\displaystyle c(p_{h},q_{h}):=\int_{\Omega}\frac{p_{h}q_{h}}{\lambda(x)}$.

We have that $a_{h}(\cdot,\cdot)$ is bounded on the broken weighted norm.
Indeed, for $\boldsymbol{u}_{h},\boldsymbol{v}_{h}\in\mathrm{H}_{h}$ and
applying (3.1) we have

$|a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})|\leq\left|\int_{\Omega}2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}):\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})\right|+\left|\int_{\mathcal{F}_{h}^{*}}\frac{\texttt{a}_{S}}{h_{\mathcal{F}_{h}^{*}}}2\mu(x)\llbracket\boldsymbol{u}_{h}\rrbracket\cdot\llbracket\boldsymbol{v}_{h}\rrbracket\right|\\\
+\left|\int_{\mathcal{F}_{h}^{*}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h})\\}\cdot\llbracket\boldsymbol{v}_{h}\rrbracket\right|+\left|\epsilon\int_{\mathcal{F}_{h}^{*}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})\\}\cdot\llbracket\boldsymbol{u}_{h}\rrbracket\right|\\\
\leq
2C\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{u}_{h}\|_{0,\Omega}\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}\\\
+\texttt{a}_{S}\|\mu^{1/2}h_{\mathcal{F}^{*}}^{-1/2}\llbracket\boldsymbol{u}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\|h_{\mathcal{F}_{h}}^{-1/2}h_{\mathcal{F}_{h}^{*}}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
+2\|h_{\mathcal{F}}^{1/2}\\{\mu(x)\nabla\boldsymbol{u}_{h}\\}\|_{0,\mathcal{F}_{h}^{*}}\|h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}+2\epsilon\|h_{\mathcal{F}}^{1/2}\\{\mu(x)\nabla\boldsymbol{v}_{h}\\}\|_{0,\mathcal{F}_{h}^{*}}\|h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{u}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
\leq
2C\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{u}_{h}\|_{0,\Omega}\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}\\\
+\texttt{a}_{S}C_{E_{\min}}\|\mu^{1/2}(x)h_{\mathcal{F}^{*}}^{-1/2}\llbracket\boldsymbol{u}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\|\mu^{1/2}(x)h_{\mathcal{F}}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
+2C_{\texttt{tr}}C_{E_{\min}}^{2}\|\mu^{1/2}(x)\nabla\boldsymbol{u}_{h}\|_{0,\Omega}\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
+2\epsilon
C_{\texttt{tr}}C_{E_{\min}}^{2}\|\mu^{1/2}(x)\nabla\boldsymbol{v}_{h}\|_{0,\Omega}\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{u}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
\leq
C_{1}\left(\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{u}_{h}\|_{0,\Omega}^{2}+\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{u}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\right)^{1/2}\\\
\times\left(\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}^{2}+\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\right)^{1/2}=C_{1}\|\boldsymbol{u}_{h}\|_{\mathrm{H}(h)}\|\boldsymbol{v}_{h}\|_{\mathrm{H}(h)},$

where
$C_{1}:=2\max\\{C,C_{E_{\min}}\texttt{a}_{S},C_{\texttt{tr}}C_{E_{\min}}^{2},\epsilon
C_{\texttt{tr}}C_{E_{\min}}^{2},1\\}$.

For $b_{h}(\cdot,\cdot)$ the continuity is as follows

$|b_{h}(\boldsymbol{v}_{h},q_{h})|\leq\left|\int_{\Omega}\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}_{h}q_{h}\right|+\left|\int_{\mathcal{F}_{h}^{*}}\\{q_{h}\\}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket\right|\\\
\leq\|\mathop{\mathrm{div}}\nolimits_{h}\boldsymbol{v}_{h}\|_{0,\Omega}\|q_{h}\|_{0,\Omega}+\|h^{1/2}_{\mathcal{F}}\\{q_{h}\\}\|_{0,\mathcal{F}_{h}^{*}}\|h_{\mathcal{F}}^{-1/2}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
\leq
C\|\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}\|q_{h}\|_{0,\Omega}+C_{\texttt{tr}}C_{E_{\min}}\|q_{h}\|_{0,\Omega}\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
\leq
CC_{E_{\min}}C_{E_{\max}}\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}\|\lambda^{-1/2}(x)q_{h}\|_{0,\Omega}\\\
+C_{\texttt{tr}}C_{E_{\max}}C_{E_{\min}}\|\mu^{-1/2}(x)q_{h}\|_{0,\Omega}\|\mu^{1/2}h_{\mathcal{F}}^{-1/2}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}\\\
\leq
C_{2}\left(\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}^{2}+\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\underline{\boldsymbol{v}_{h}}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\right)^{1/2}\\\
\times(\|\mu^{-1/2}(x)q_{h}\|_{0,\Omega}^{2}+\|\lambda^{-1/2}(x)q_{h}\|_{0,\Omega}^{2})^{1/2}=C_{2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v}_{h},q_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q},$

where
$C_{2}:=\max\\{CC_{E_{\min}}C_{E_{\max}},C_{\texttt{tr}}C_{E_{\max}}C_{E_{\min}}\\}$.

######  Remark 3.1.

Observe that the parameter $\epsilon\in\\{-1,0,1\\}$ dictates if the IPDG
method is symmetric or non-symmetric. More precisely, if $\epsilon=1$ we
obtain the classic symmetric interior penalty method (SIP), if $\epsilon=-1$
we obtain the non-symmetric interior penalty method (NIP) and if $\epsilon=0$
the incomplete interior penalty method (IIP). Clearly when non-symmetric
methods are considered, complex computed eigenvalues are expected when the
spectrum is approximated. Moreover, the SIP attains optimal order of
convergence for the approximation of the eigenvalues and eigenfunctions,
whereas the NIP and IIP methods deliver suboptimal orders.

Finally, the IPDG discretization of (2.3) reads as follows: Find
$\kappa_{h}\in\mathbb{C}$ and
$(\boldsymbol{0},0)\neq(\boldsymbol{u}_{h},p_{h})\in\mathrm{H}_{h}\times\mathrm{Q}_{h}$
such that

(3.4)
$\mathcal{A}_{h}((\boldsymbol{u}_{h},p_{h}),(\boldsymbol{v}_{h},q_{h}))=\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})\quad\forall(\boldsymbol{v}_{h},q_{h})\in\mathrm{H}_{h}\times\mathrm{Q}_{h}.$

To establish the well posedness of our discrete problem (3.2), we have from
[12, Proposition 10] the following discrete inf-sup condition

$\displaystyle\sup_{\boldsymbol{\tau}_{h}\in\mathrm{H}_{h}}\frac{b_{h}(\boldsymbol{\tau}_{h},q_{h})}{\|\boldsymbol{\tau}_{h}\|_{\mathrm{H}(h)}}\geq\beta\|\mu^{-1/2}(x)\,q_{h}\|_{0,\Omega}\quad\forall
q_{h}\in\mathrm{Q}_{h},$

where the constant $\beta>0$ is independent of $h$. On the other hand, let us
define the discrete kernel $\mathcal{K}_{h}$ of $b_{h}(\cdot,\cdot)$ as
follows

$\mathcal{K}_{h}:=\\{\boldsymbol{\tau}_{h}\in\mathrm{H}_{h}\,:\,b_{h}(\boldsymbol{\tau}_{h},q_{h})=0\,\,\,\,\forall
q_{h}\in\mathrm{Q}_{h}\\}.$

With this space at hand, we prove the following coercivity result for
$a_{h}(\cdot,\cdot)$.

######  Lemma 3.2 (ellipticity of $a_{h}(\cdot,\cdot)$).

For any $\varepsilon\in\\{-1,0,1\\}$, there exists a positive parameter
$\texttt{a}^{*}$ such that for all $\texttt{a}_{S}\geq\texttt{a}^{*}$ there
holds

$a_{h}(\boldsymbol{v}_{h},\boldsymbol{v}_{h})\geq\alpha\|\boldsymbol{v}_{h}\|_{\mathrm{H}(h)}^{2}\quad\forall\boldsymbol{v}_{h}\in\mathcal{K}_{h},$

where $\alpha>0$ is independent of $h$.

###### Proof.

Let $\boldsymbol{v}_{h}\in\mathcal{K}_{h}$. Hence, applying Cauchy-Schwarz and
Korn’s inequalities and the the well known identity $\displaystyle
ab\leq\frac{a^{2}}{2\eta}+\frac{\eta b^{2}}{2}$ for $a,b\in\mathbb{R}$ and
$\eta>0$, we have

$a_{h}(\boldsymbol{v}_{h},\boldsymbol{v}_{h})=\int_{\Omega}2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h}):\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})+\int_{\mathcal{F}_{h}^{*}}\frac{\texttt{a}_{S}}{h_{\mathcal{F}}}2\mu(x)\llbracket\boldsymbol{v}_{h}\rrbracket\cdot\llbracket\boldsymbol{v}_{h}\rrbracket\\\
-2(1+\epsilon)\int_{\mathcal{F}_{h}^{*}}\\{\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h})\\}\cdot\llbracket\boldsymbol{v}_{h}\rrbracket\\\
\geq
2C_{\texttt{Korn}}E_{\min}\|\widetilde{\mu}^{1/2}\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}^{2}+2\texttt{a}_{S}\|\mu^{1/2}h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\\\
+2(1+\epsilon)E_{\max}\left(-\frac{\|h_{\mathcal{F}}^{1/2}\widetilde{\mu}^{1/2}\boldsymbol{\varepsilon}(\boldsymbol{v}_{h})\|_{0,\mathcal{F}_{h}^{*}}^{2}}{2\eta}-\eta\frac{\|h_{\mathcal{F}}^{-1/2}\llbracket\widetilde{\mu}^{1/2}\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}}{2}\right)\\\
\geq\frac{1}{E_{\max}}\underbrace{\left(2C_{\texttt{Korn}}E_{\min}-\frac{2(1+\epsilon)C_{E_{\max}}}{2\eta}\right)}_{:=\widehat{C}_{1}}\|\mu^{1/2}(x)\nabla_{h}\boldsymbol{v}_{h}\|_{0,\Omega}^{2}\\\
+\underbrace{\left(\texttt{a}_{S}-\frac{E_{\max}(1+\epsilon)\eta}{2E_{\min}}\right)}_{:=\widehat{C}_{2}}\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\boldsymbol{v}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}^{2}\geq\min\\{\widehat{C}_{1},\widehat{C}_{2}\\}\|\boldsymbol{v}_{h}\|_{\mathrm{H}(h)}^{2}.$

We observe that to conclude the proof, the parameter $\eta$ must be chosen as

$\eta>\frac{(1+\epsilon)C_{E_{\max}}}{2C_{\texttt{Korn}}E_{\min}}\quad\text{and
hence}\quad\texttt{a}_{S}>\texttt{a}^{*}:=\frac{(1+\epsilon)^{2}C_{E_{\max}}E_{\max}}{4E_{\min}^{2}C_{\texttt{Korn}}}.$

This leads to the positiveness of $\widehat{C}_{1}$ and $\widehat{C}_{2}$ and
hence, the proof. ∎

Let us introduce the discrete solution operator defined by

$\boldsymbol{T}_{h}:\mathrm{H}\rightarrow\mathrm{H}_{h},\qquad\boldsymbol{f}\mapsto\boldsymbol{T}_{h}\boldsymbol{f}:=\widehat{\boldsymbol{u}}_{h},$

where given $\boldsymbol{f}\in\mathrm{H}$, the pair
$(\widehat{\boldsymbol{u}}_{h},\widehat{p}_{h})\in\mathrm{H}_{h}\times\mathrm{Q}_{h}$
solves the following discrete source problem

$\mathcal{A}_{h}((\widehat{\boldsymbol{u}}_{h},\widehat{p}_{h}),(\boldsymbol{v}_{h},q_{h}))=(\boldsymbol{f},\boldsymbol{v}_{h})\quad\forall(\boldsymbol{v}_{h},q_{h})\in\mathrm{H}_{h}\times\mathrm{Q}_{h}.$

It is easy to check that there exists a positive constant $\widehat{\alpha}$
depending on the constant provided by Lemma 3.2 such that

(3.5)
$\mathcal{A}_{h}((\widehat{\boldsymbol{u}}_{h},\widehat{p}_{h}),(\widehat{\boldsymbol{u}}_{h},\widehat{p}_{h}))\geq\widehat{\alpha}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}_{h},p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\mathrm{H}(h)\times\mathrm{Q}_{h}}\quad\forall(\boldsymbol{u}_{h},p_{h})\in\mathrm{H}(h)\times\mathrm{Q}_{h}.$

##  4\. Convergence and error estimates

The aim of this section is to derive convergence results and error estimates
for our DG methods. Despite to the fact that $\boldsymbol{T}$ is compact, the
classic theory of compact operators is not enough to conclude the convergence
in norm between the continuous and discrete solution operators, since the
numerical method of our interest is non conforming. Is this reason, and
following the spirit of [2], why we resort to the theory of non compact
operators of [6, 7].

In what follows, we will denote by
$\lVert\cdot\rVert_{\mathcal{L}(\mathrm{H}(h),\mathrm{H}(h))}$ the
corresponding norm acting from $\mathrm{H}(h)$ into the same space. In
addition, we will denote by
$\lVert\cdot\rVert_{\mathcal{L}(\mathrm{H}_{h},\mathrm{H}(h))}$ the norm of an
operator restricted to the discrete subspace $\mathrm{H}_{h}$; namely, if
$\boldsymbol{L}:\mathrm{H}(h)\to\mathrm{H}(h)$, then

$\lVert\boldsymbol{L}\rVert_{\mathcal{L}(\mathrm{H}_{h},\mathrm{H}(h))}:=\sup_{\mathbf{0}\neq\boldsymbol{\tau}_{h}\in\mathrm{H}_{h}}\frac{\lVert\boldsymbol{L}\boldsymbol{\tau}_{h}\rVert_{\mathrm{H}(h)}}{\lVert\boldsymbol{\tau}_{h}\rVert_{\mathrm{H}(h)}}.$

Our first task is to prove the following properties

  * •

P1.
$\lVert\boldsymbol{T}-\boldsymbol{T}_{h}\rVert_{\mathcal{L}(\mathrm{H}_{h},\mathrm{H}(h))}\to
0$ as $h\to 0$.

  * •

P2. $\forall\boldsymbol{\tau}\in\mathrm{H}$, there holds
$\inf_{\boldsymbol{\tau}_{h}\in\mathrm{H}_{h}}\lVert\boldsymbol{\tau}-\boldsymbol{\tau}_{h}\rVert_{\mathrm{H}(h)}\,\,\to
0\quad\text{as}\quad h\to 0$.

We need this properties in order to establish spectral correctness (see [6])
for all the discrete methods (symmetric or nonsymmetric). Property P2 is
immediate as a consequence of the density of continuous piecewise degree $k$
polynomial functions in $\mathrm{H}_{h}$. The task now is to prove P1.

The following convergence result holds for the continuous and discrete
solution operators.

######  Lemma 4.1.

For all $\boldsymbol{f}\in\mathrm{H}$, the following estimate holds

$\|(\boldsymbol{T}-\boldsymbol{T}_{h})\boldsymbol{f}\|_{\mathrm{H}(h)}\lesssim
h^{s}\|\boldsymbol{f}\|_{0,\Omega},$

where $s>0$ and the hidden constant is independent of $h$.

###### Proof.

For $\boldsymbol{f}\in\mathrm{H}$ we have
$\widehat{\boldsymbol{u}}:=\boldsymbol{T}\boldsymbol{f}$ and
$\widehat{\boldsymbol{u}}_{h}:=\boldsymbol{T}_{h}\boldsymbol{f}$. Then, there
exists a constant $\widetilde{C}>0$ such that the following Céa estimate holds

$\|(\boldsymbol{T}-\boldsymbol{T}_{h})\boldsymbol{f}\|_{\mathrm{H}(h)}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\widehat{\boldsymbol{u}}-\widehat{\boldsymbol{u}}_{h},\widehat{p}-\widehat{p}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times\mathrm{Q}_{h}}\\\
\leq\widetilde{C}\inf_{(\widehat{\boldsymbol{v}}_{h},\widehat{q}_{h})\in\mathrm{H}(h)\times
Q_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\widehat{\boldsymbol{u}}-\widehat{\boldsymbol{v}}_{h},\widehat{p}-\widehat{q}_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times\mathrm{Q}_{h}}\\\
=\widetilde{C}\left(\inf_{\widehat{\boldsymbol{v}}_{h}\in\mathrm{H}(h)}\|\widehat{\boldsymbol{u}}_{h}-\widehat{\boldsymbol{v}}_{h}\|_{\mathrm{H}(h)}+\inf_{\widehat{q}_{h}\in\mathrm{Q}_{h}}\|\widehat{p}-\widehat{q}_{h}\|_{\mathrm{Q}_{h}}\right).$

For the first infimum on the identity above, if $\Pi_{h}$ represents the
Lagrange interpolation operator, we have

$\displaystyle\inf_{\widehat{\boldsymbol{v}}_{h}\in\mathrm{H}(h)}$
$\displaystyle\|\widehat{\boldsymbol{u}}_{h}-\widehat{\boldsymbol{v}}_{h}\|_{\mathrm{H}(h)}\leq\|\widehat{\boldsymbol{u}}_{h}-\Pi_{h}\widehat{\boldsymbol{u}}_{h}\|_{\mathrm{H}(h)}$
$\displaystyle=\|\mu^{1/2}(x)\nabla(\widehat{\boldsymbol{u}}_{h}-\Pi_{h}\widehat{\boldsymbol{u}}_{h})\|_{0,\Omega}+\|\mu^{1/2}(x)h_{\mathcal{F}}^{-1/2}\llbracket\widehat{\boldsymbol{u}}_{h}-\Pi_{h}\widehat{\boldsymbol{u}}_{h}\rrbracket\|_{0,\mathcal{F}_{h}^{*}}$
(4.1) $\displaystyle\leq
C_{E_{\max}}\|\widehat{\boldsymbol{u}}_{h}-\Pi_{h}\widehat{\boldsymbol{u}}_{h}\|_{1,\Omega}\leq
C_{E_{\max}}C\|\widehat{\boldsymbol{u}}\|_{1+s}\leq
C_{E_{\max}}C\|\boldsymbol{f}\|_{0,\Omega}.$

For the second infimum, if $\textrm{P}_{h}$ represents the
$\mathrm{L}^{2}$-projection operator, we have

$\displaystyle\inf_{\widehat{q}_{h}\in\mathrm{Q}_{h}}\|\widehat{p}-\widehat{q}_{h}\|_{\mathrm{Q}}$
$\displaystyle\leq\|\widehat{p}-\textrm{P}_{h}\widehat{p}\|_{\mathrm{Q}}=\|\mu^{-1/2}(x)(\widehat{p}-\textrm{P}_{h}\widehat{p})\|_{0,\Omega}+\|\lambda^{-1/2}(x)(\widehat{p}-\textrm{P}_{h}\widehat{p})\|_{0,\Omega}$
(4.2)
$\displaystyle\leq(E_{\min}+E_{\min})\|\widehat{p}-\textrm{P}_{h}\widehat{p}\|_{0,\Omega}\leq
C(2E_{\min})\|\widehat{p}\|_{s,\Omega}\leq
CE_{\min}\|\boldsymbol{f}\|_{0,\Omega}.$

Then, combining (4) and (4) we conclude the proof. ∎

Now we prove the analogous of the previous lemma, but considering discrete
sources. Since the proof is analogous to the above lemma, we do not include
the details and present the result as the following corollary.

######  Corollary 4.2.

For all $\boldsymbol{f}_{h}\in\mathrm{H}_{h}$, the following estimate holds

$\|(\boldsymbol{T}-\boldsymbol{T}_{h})\boldsymbol{f}_{h}\|_{\mathrm{H}(h)}\lesssim
h^{s}\|\boldsymbol{f}_{h}\|_{\mathrm{H}(h)},$

where $s$ is as in Theorem 2.3 and the hidden constant is independent of $h$.

Now we are in position to establish P1 and its proof is identical to the one
in [17, Lemma 3].

######  Lemma 4.3.

There following estimate holds

$\lVert\boldsymbol{T}-\boldsymbol{T}_{h}\rVert_{\mathcal{L}(\mathrm{H}_{h},\mathrm{H}(h))}\lesssim
h^{s},$

where the hidden constant is independent of $h$.

From now on, $\mathcal{D}$ denotes the unitary disk defined in the complex
plane by $\mathcal{D}:=\\{z\in\mathbb{C}\,:\,|z|\leq 1\\}$. The following
result proves that the continuous resolvent is bounded in the $\mathrm{H}(h)$
norm.

######  Lemma 4.4.

There exists a constant $C>0$ independent of $h$ such that for all
$z\in\mathcal{D}\setminus\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})$ there
holds

$\lVert(z\boldsymbol{I}-\boldsymbol{T})\boldsymbol{f}\rVert_{\mathrm{H}(h)}\geq
C|z|\,\lVert\boldsymbol{f}\rVert_{\mathrm{H}(h)}\quad\forall\boldsymbol{f}\in\mathrm{H}(h).$

######  Remark 4.5.

Lemma 4.4 implies that the resolvent of $\boldsymbol{T}$ is bounded. In other
words, if $J$ is a compact subset of
$\mathcal{D}\setminus\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})$, then
there exists a positive constant $C$ independent of $h$, satisfying the
estimate

$\displaystyle\lVert(z\boldsymbol{I}-\boldsymbol{T})^{-1}\rVert_{\mathcal{L}(\mathrm{H}(h),\mathrm{H}(h))}\leq
C\qquad\forall z\in J.$

Our next goal is to derive the boundedness of the discrete resolvent, when $h$
is small enough.

######  Lemma 4.6.

If $z\in\mathcal{D}\setminus\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})$,
there exists $h_{0}>0$ such that for all $h\leq h_{0}$,

$\lVert(z\boldsymbol{I}-\boldsymbol{T}_{h})\boldsymbol{f}\rVert_{\mathrm{H}(h)}\geq
C\,\lVert\boldsymbol{f}\rVert_{\mathrm{H}(h)}\quad\forall\boldsymbol{f}\in\mathrm{H}(h),$

with $C>0$ independent of $h$ but depending on $|z|$.

The previous lemma states that if we consider a compact subset $E$ of the
complex plane such that
$E\cap\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})=\emptyset$ for $h$ small
enough and for all $z\in E$, operator $z\boldsymbol{I}-\boldsymbol{T}_{h}$ is
invertible. Moreover, there exists a positive constant $C$ independent of $h$
such that
$\lVert(z\boldsymbol{I}-\boldsymbol{T}_{h})^{-1}\rVert_{\mathcal{L}(\mathrm{H}(h),\mathrm{H}(h))}\leq
C$ for all $z\in E$. This fact is important since it determines that the
numerical method is spurious free for $h$ small enough. This is summarized in
the following result proved in [6].

######  Theorem 4.7.

Let $E\subset\mathbb{C}$ be a compact subset not intersecting
$\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})$. Then, there exists $h_{0}>0$
such that, if $h\leq h_{0}$, then
$E\cap\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T}_{h})=\emptyset.$

We recall the definition of the resolvent operator of $\boldsymbol{T}$ and
$\boldsymbol{T}_{h}$ respectively:

$\displaystyle(z\boldsymbol{I}-\boldsymbol{T})^{-1}\,:\,\mathrm{H}\to\mathrm{H}\,,\quad
z\in\mathbb{C}\setminus\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T}),$
$\displaystyle(z\boldsymbol{I}-\boldsymbol{T}_{h})^{-1}\,:\,\mathrm{H}_{h}\to\mathrm{H}_{h}\,,\quad
z\in\mathbb{C}\setminus\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T}_{h}).$

We introduce the definition of the gap $\widehat{\delta}$ between two closed
subspaces $\mathcal{X}$ and $\mathcal{Y}$ of $\mathrm{L}^{2}(\Omega)$:

$\widehat{\delta}(\mathcal{X},\mathcal{Y}):=\max\big{\\{}\delta(\mathcal{X},\mathcal{Y}),\delta(\mathcal{Y},\mathcal{X})\big{\\}},$

where

$\delta(\mathcal{X},\mathcal{Y}):=\sup_{x\in\mathcal{X}:\
\left\|x\right\|_{0,\Omega}=1}\left(\inf_{y\in\mathcal{Y}}\left\|x-y\right\|_{0,\Omega}\right).$

Let $\kappa$ be an isolated eigenvalue of $\boldsymbol{T}$ and let $D$ be an
open disk in the complex plane with boundary $\gamma$ such that $\kappa$ is
the only eigenvalue of $\boldsymbol{T}$ lying in $D$ and
$\gamma\cap\mathop{\mathrm{sp}}\nolimits(\boldsymbol{T})=\emptyset$. We
introduce the spectral projector corresponding to the continuous and discrete
solution operators $\boldsymbol{T}$ and $\boldsymbol{T}_{h}$, respectively

$\displaystyle\boldsymbol{\mathfrak{E}}:=\frac{1}{2\pi
i}\int_{\gamma}\left(z\boldsymbol{I}-\boldsymbol{T}\right)^{-1}\,dz:\mathrm{H}(h)\longrightarrow\mathrm{H}(h),$
$\displaystyle\boldsymbol{\mathfrak{E}}_{h}:=\frac{1}{2\pi
i}\int_{\gamma}\left(z\boldsymbol{I}-\boldsymbol{T}_{h}\right)^{-1}\,dz:\mathrm{H}(h)\longrightarrow\mathrm{H}(h).$

The following approximation result for the spectral projections holds.

######  Lemma 4.8.

There holds

$\lim_{h\rightarrow
0}\|\boldsymbol{\mathfrak{E}}-\boldsymbol{\mathfrak{E}}_{h}\|_{\mathcal{L}(\mathrm{H}_{h},\mathrm{H}(h))}=0.$

###### Proof.

See [2, Theorem 5.1]. ∎

The next result provides, for the proposed IPDG methods, an error estimate for
the eigenfunctions and a double order of convergence for the eigenvalues. More
precisely, the result presents estimates for symmetric and nonsymmetric
methods, where optimal and suboptimal order of convergence are attained,
respectively.

######  Theorem 4.9.

There exists a strictly positive constant $h_{0}$ such that, for
$\varepsilon\in\\{-1,0,1\\}$ and $h<h_{0}$ there holds

$\widehat{\delta}_{h}(\boldsymbol{\mathfrak{E}}(\mathrm{H},\boldsymbol{\mathfrak{E}}_{h}(\mathrm{H}_{h}))\lesssim
h^{\min\\{r,k\\}},\quad\text{ and }\quad\max_{1\leq i\leq
m}|\kappa-\kappa_{i,h}|\lesssim\,h^{\sigma\min\\{r,k\\}},$

where $\sigma:=\frac{1}{2}\left[(\varepsilon+2)^{\varepsilon}-1\right]+1$ and
the hidden constants are independent of $h$.

###### Proof.

The proof for the eigenfunction estimate follows the same arguments of [17,
Lemma 7]. Let
$(\kappa_{i,h},(\boldsymbol{u}_{h},p_{h}))\in\mathbb{C}\times\mathrm{H}_{h}\times\mathrm{Q}_{h}$
be the solution of (3.2) with $i=1,\ldots,m$ and let
$(\kappa,(\boldsymbol{u},p))\in\mathbb{R}\times\mathrm{H}\times\mathrm{Q}$.
Let us consider the following well known algebraic identity

(4.3)
$\mathcal{A}_{h}((\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}),(\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}))-\kappa(\boldsymbol{u}-\boldsymbol{u}_{h},\boldsymbol{u}-\boldsymbol{u}_{h})=\left(\kappa_{i,h}-\kappa\right)(\boldsymbol{u}_{h},\boldsymbol{u}_{h}).$

Hence, taking modulus to $\mathcal{A}(\cdot,\cdot)$ we have

$\displaystyle|\mathcal{A}_{h}((\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}),(\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}))|$
$\displaystyle\leq|a_{h}(\boldsymbol{u}-\boldsymbol{u}_{h},,\boldsymbol{u}-\boldsymbol{u}_{h},)|+|c(p-p_{h},p-p_{h})|$
(4.4)
$\displaystyle\leq\|\boldsymbol{u}-\boldsymbol{u}_{h}\|^{2}+\|p-p_{h}\|^{2}$
$\displaystyle\leq\widehat{\delta}_{h}(\boldsymbol{\mathfrak{E}}(\mathrm{H},\boldsymbol{\mathfrak{E}}_{h}(\mathrm{H}_{h}))\lesssim
h^{2\min\\{r,k\\}}.$

On the other hand

(4.5)
$\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,\Omega}^{2}\leq\widehat{\delta}_{h}(\boldsymbol{\mathfrak{E}}(\mathrm{H},\boldsymbol{\mathfrak{E}}_{h}(\mathrm{H}_{h}))\lesssim
h^{2\min\\{r,k\\}}.$

Also, invoking (3.5) we have

(4.6)
$\displaystyle(\boldsymbol{u}_{h},\boldsymbol{u}_{h})=\frac{\mathcal{A}_{h}((\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}),(\boldsymbol{u}-\boldsymbol{u}_{h},p-p_{h}))}{|\kappa_{i,h}|}\geq\frac{\widehat{\alpha}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}_{h},p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\mathrm{H}(h)\times\mathrm{Q}_{h}}}{|\kappa_{i,h}|}\geq\widehat{C}>0.$

Hence, plugin (4), (4.5), and (4.6) in (4.3) we conclude the proof. ∎

##  5\. A posteriori error analysis

This section aims to present a viable residual-based error estimator for the
Stokes/elasticity eigenvalue problem. The primary objective is to demonstrate
the equivalence between the proposed estimator and the actual error.
Furthermore, in the subsequent analysis, our focus will be exclusively on
eigenvalues with a multiplicity of 1. Having established the a priori error
estimates accurately in terms of the stabilization constant, the bounding
constants will not be followed up during this section.

###  5.1. The local and global indicators

In this section, we introduce an a posteriori error estimator based on
residuals for the elasticity eigenvalue problem. The proposed estimator also
works with the limit case $\lambda\rightarrow\infty$. For non-symmetric IPDG
methods, the eigenvalue can be complex. However, from [19] it is observed that
the imaginary part is negligible. Hence, we provide the analysis for the real
part of the spectrum, which works for symmetric and non-symmetric IPDG.

Consider an eigenpair approximation
$(\kappa_{h},(\boldsymbol{u}_{h},p_{h}))\in\mathbb{R}^{+}\times\mathrm{H}_{h}\times\mathrm{Q}_{h}$.
For every element $K\in\mathcal{T}_{h}$, the interior residual estimator
$\eta_{R_{K}}$ is defined as follows:

$\displaystyle\eta_{R_{K}}^{2}:=$
$\displaystyle\|(2\mu_{h})^{-1/2}h_{K}(\kappa_{h}\boldsymbol{u}_{h}+\nabla\cdot(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}))-\nabla
p_{h})\|_{0,K}^{2}$
$\displaystyle+\|((2\mu_{h})^{-1}+\lambda^{-1})^{-1/2}[\nabla\cdot\boldsymbol{u}_{h}+(1/\lambda)\,p_{h}]\|_{0,K}^{2},$

and the facet residual estimator $\eta_{\mathcal{F}_{K}}$ by

$\displaystyle\eta_{\mathcal{F}_{K}}^{2}$
$\displaystyle:=\sum_{\mathcal{F}\in\partial
K\cap\mathcal{F}_{h}^{0}}\|(2\mu_{h})^{-1/2}h_{\mathcal{F}}^{1/2}[\\![(p_{h}\boldsymbol{I}-2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}))\textbf{n}]\\!]\|^{2}_{0,\mathcal{F}}$
$\displaystyle\hskip 113.81102pt+\sum_{\mathcal{F}\in\partial
K\cap\mathcal{F}_{h}^{N}}\|(2\mu_{h})^{-1/2}h_{\mathcal{F}}^{1/2}(p_{h}\boldsymbol{I}-2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}))\textbf{n}\|^{2}_{0,\mathcal{F}},$

where $\boldsymbol{I}$ represents the identity matrix. Moving forward, we
present the estimator $\eta_{J_{K}}$, designed to quantify the jump in the
approximate solution $\boldsymbol{u}_{h}$.

$\eta_{J_{K}}^{2}:=\sum_{\mathcal{F}\in\partial
K\cap\mathcal{F}_{h}^{0}}\|(2\mu_{h})^{1/2}\gamma_{h}^{1/2}[\\![\boldsymbol{u}_{h}]\\!]\|^{2}_{0,\mathcal{F}}+\sum_{\mathcal{F}\in\partial
K\cap\mathcal{F}_{h}^{D}}\|(2\mu_{h})^{1/2}\gamma_{h}^{1/2}\boldsymbol{u}_{h}\otimes\boldsymbol{n}\|^{2}_{0,\mathcal{F}},$

with $\gamma_{h}=\texttt{a}_{S}/h_{\mathcal{F}}$, where $\texttt{a}_{S}$
denotes the penalty parameter as discussed in the section 3. The local error
indicator, obtained by summing the three aforementioned terms, is defined as:

$\eta_{K}^{2}:=\eta_{R_{K}}^{2}+\eta_{\mathcal{F}_{K}}^{2}+\eta_{J_{K}}^{2}.$

Lastly, we present the (global) a posteriori error estimator

(5.1) $\eta:=\left(\sum_{K\in\mathcal{T}_{h}}\eta_{K}^{2}\right)^{1/2},$

and the data oscillation term is given by

$\Theta:=\left(\sum_{K\in\mathcal{T}_{h}}\Theta_{K}^{2}\right)^{1/2},$

where
$\Theta_{K}:=\|(2\mu_{h})^{-1/2}(\mu-\mu_{h})\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h})\|_{0,K}$.

###  5.2. Reliability

Following the approach in [13, 14], we define
$\mathrm{H}_{h}^{c}=\mathrm{H}_{h}\cap\mathrm{H}$. The orthogonal complement
of $\mathrm{H}_{h}^{c}$ in $\mathrm{H}_{h}$ with respect to the norm
${\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}$
is denoted as $\mathrm{H}_{h}^{\perp}$. Consequently, we have the
decomposition $\mathrm{H}_{h}=\mathrm{H}_{h}^{c}\oplus\mathrm{H}_{h}^{\perp}$,
allowing us to uniquely decompose the DG velocity approximation into
$\boldsymbol{u}_{h}=\boldsymbol{u}_{h}^{c}+\boldsymbol{u}_{h}^{r}$, where
$\boldsymbol{u}_{h}^{c}\in\mathrm{H}_{h}^{c}$ and
$\boldsymbol{u}_{h}^{r}\in\mathrm{H}_{h}^{\perp}$. Applying the triangle
inequality, we can express:

${\left|\kern-1.07639pt\left|\boldsymbol{u}-\boldsymbol{u}_{h}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}\leq{\left|\kern-1.07639pt\left|\boldsymbol{u}-\boldsymbol{u}_{h}^{c}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}+{\left|\kern-1.07639pt\left|\boldsymbol{u}_{h}^{r}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$

Referring to [13, Proposition 4.1], we derive the upper bound for the second
term.

(5.2)
${\left|\kern-1.07639pt\left|\boldsymbol{u}_{h}^{r}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}\lesssim\left(\sum_{K\in\mathcal{T}_{h}}\eta_{J_{K}}^{2}\right)^{1/2}.$

It is important to remark that the DG form $a_{h}(\cdot,\cdot)$ lacks a well-
defined definition for functions $\boldsymbol{u},\boldsymbol{v}$ within
$\mathrm{H}$. This challenge can be addressed by employing an appropriate
lifting operator, as outlined in [5, 14]. However, we explore an alternative
approach here, wherein the DG form $a_{h}(\cdot,\cdot)$ is decomposed into
multiple components.

$a_{h}(\boldsymbol{u},\boldsymbol{v})=(2\mu\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}),\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}))+C_{h}(\boldsymbol{u},\boldsymbol{v})+J_{h}(\boldsymbol{u},\boldsymbol{v}),$

where $C_{h}$ and $J_{h}$ are defined by

$\displaystyle C_{h}(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle:=-\int_{\mathcal{F}^{*}_{h}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{u})\\}\cdot\llbracket\boldsymbol{v}\rrbracket-\epsilon\int_{\mathcal{F}^{*}_{h}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v})\\}\cdot\llbracket\boldsymbol{u}\rrbracket,$
$\displaystyle J_{h}(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle:=\int_{\mathcal{F}^{*}_{h}}\frac{\texttt{a}_{S}}{h_{\mathcal{F}}}{2\mu(x)}\llbracket\boldsymbol{u}\rrbracket\cdot\llbracket\boldsymbol{v}\rrbracket.$

The following result provides an estimate for the term $C_{h}(\cdot,\cdot)$ in
terms of the indicator $\eta_{J_{K}}$.

######  Lemma 5.1.

Suppose $\boldsymbol{u}_{h}\in\mathrm{H}_{h}$ and
$\boldsymbol{v}_{h}^{c}\in\mathrm{H}_{h}^{c}$; then, the following relation
holds:

$C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})\lesssim|\epsilon|\texttt{a}_{S}^{-1/2}\left(\sum_{K\in\mathcal{T}_{h}}\eta_{J_{K}}^{2}\right)^{1/2}{\left|\kern-1.07639pt\left|\boldsymbol{v}^{c}_{h}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$

###### Proof.

Since $\boldsymbol{v}_{h}^{c}\in\mathrm{H}_{h}^{c}$, it follows that

$C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})=-\epsilon\int_{\mathcal{F}^{*}_{h}}\\{2\mu(x)\boldsymbol{\varepsilon}_{h}(\boldsymbol{v})\\}\cdot\llbracket\boldsymbol{u}\rrbracket.$

Utilizing the Cauchy-Schwarz inequality yields:

$\displaystyle C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})$
$\displaystyle\lesssim|\epsilon|\left(\sum_{\mathcal{F}\in\mathcal{F}_{h}^{*}}\gamma_{h}^{-1}\|[\\![(2\mu)^{1/2}\nabla\boldsymbol{v}_{h}^{c}]\\!]\|_{0,\mathcal{F}}^{2}\right)^{1/2}\left(\sum_{\mathcal{F}\in\mathcal{F}_{h}^{*}}\gamma_{h}\|(2\mu)^{1/2}[\\![\boldsymbol{u}_{h}]\\!]\|_{0,\mathcal{F}}^{2}\right)^{1/2}.$

Applying a trace estimate in conjunction with a discrete inverse inequality
for an edge/face $\mathcal{F}\in\mathcal{F}_{h}^{*}$, where
$\mathcal{F}=K_{1}\cap K_{2}$ if $\mathcal{F}_{h}^{0}$ and
$\mathcal{F}=K_{1}$, $K_{2}=\emptyset$ if $\mathcal{F}\subset\partial\Omega$,
results in:

$\|\nabla\boldsymbol{v}_{h}^{c}\|_{0,E}\lesssim
h_{K}^{-1/2}\|\nabla\boldsymbol{v}_{h}^{c}\|_{0,K_{1}\cup K_{2}}.$

Thus we have

$C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})\lesssim|\epsilon|\texttt{a}_{S}^{-1/2}\left(\sum_{K\in\mathcal{T}_{h}}\|(2\mu)^{1/2}\nabla\boldsymbol{v}_{h}^{c}\|^{2}_{0,K}\right)^{1/2}\,\left(\sum_{K\in\mathcal{T}_{h}}\eta_{J_{K}}^{2}\right)^{1/2}.$

This concludes the proof. ∎

Let us denote by ${\Pi}_{h}:\mathrm{H}\to\mathrm{H}_{h}^{c}$ the Scott-Zhang
interpolation operator [22] (see also [14, eq. (4.5) and eq. (4.6)]), which is
stable
$\|\nabla({\Pi}_{h}\boldsymbol{v})\|_{0}\lesssim\|\nabla\boldsymbol{v}\|_{0}$
and satisfies the following interpolation property

(5.3)
$\sum_{K\in\mathcal{T}_{h}}h_{K}^{-2}\|\boldsymbol{v}-{\Pi}_{h}\boldsymbol{v}\|^{2}_{0,K}+\sum_{\mathcal{F}\in\mathcal{F}_{h}^{*}}h_{\mathcal{F}}^{-1}\|\boldsymbol{v}-\Pi_{h}\boldsymbol{v}\|^{2}_{0,E}\lesssim\|\nabla\boldsymbol{v}\|_{0}^{2},$

for any $\boldsymbol{v}\in\mathrm{H}$.

######  Lemma 5.2.

Let $\boldsymbol{v}_{h}^{c}=\Pi_{h}\boldsymbol{v}\in\mathrm{H}_{h}^{c}$, where
$\Pi_{h}$ denotes the Scott-Zhang interpolation of
$\boldsymbol{v}\in\mathrm{H}$. For any $\boldsymbol{u}_{h}\in\mathrm{H}_{h}$,
$p_{h}\in Q_{h}$, and $\kappa_{h}\in\mathbb{R}^{+}$, the following estimate
holds:

$\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}-\boldsymbol{v}_{h}^{c})-(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}-\boldsymbol{v}_{h}^{c}))+(p_{h},\nabla\cdot(\boldsymbol{v}-\boldsymbol{v}_{h}^{c}))\lesssim\eta{\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$

###### Proof.

Employing integration by parts on each element $K\in\mathcal{T}_{h}$, we
obtain:

$\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}-\boldsymbol{v}_{h}^{c})-(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}-\boldsymbol{v}_{h}^{c}))+(p_{h},\nabla\cdot(\boldsymbol{v}-\boldsymbol{v}_{h}^{c}))\\\
\qquad=\sum_{K\in\mathcal{T}_{h}}\int_{K}(\kappa_{h}\boldsymbol{u}_{h}+\nabla\cdot(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}))-\nabla
p_{h})(\boldsymbol{v}-\boldsymbol{v}_{h}^{c})d\boldsymbol{x}\\\
\qquad\qquad+\sum_{K\in\mathcal{T}_{h}}\int_{\partial
K}(p_{h}\boldsymbol{I}-2\mu_{h}\boldsymbol{\varepsilon}_{h}({\boldsymbol{u}}_{h})\textbf{n}_{K}\cdot(\boldsymbol{v}-\boldsymbol{v}_{h}^{c})d\boldsymbol{s}=T_{1}+T_{2}.$

Now the task is to estimate the terms $T_{1}$ and $T_{2}$. Let us begin with
$T_{1}$. By applying the Cauchy-Schwarz inequality and (5.3), we arrive at:

$T_{1}\lesssim\left(\sum_{K\in\mathcal{T}_{h}}\|h_{K}(2\mu_{h})^{-1/2}(\kappa_{h}\boldsymbol{u}_{h}+\nabla\cdot
2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h})-\nabla
p_{h})\|_{0,K}^{2}\right)^{1/2}\\\
\times\left(\sum_{K\in\mathcal{T}_{h}}h_{K}^{-2}\|(2\mu_{h})^{1/2}(\boldsymbol{v}-\boldsymbol{v}_{h}^{c})\|^{2}_{0,K}\right)^{1/2}\lesssim\left(\sum_{K\in\mathcal{T}_{h}}\eta_{R_{K}}^{2}\right)^{1/2}{\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$

Given that
$(\boldsymbol{v}-\boldsymbol{v}_{h}^{c})|_{\partial\Omega}=\boldsymbol{0}$, we
can express $T_{2}$ in terms of a sum over interior facets:

$T_{2}=\sum_{\mathcal{F}_{h}^{0}}\int_{E}[\\![(p_{h}\boldsymbol{I}-2\mu_{h}\boldsymbol{\varepsilon}\boldsymbol{u}_{h})\boldsymbol{n}]\\!](\boldsymbol{v}-\boldsymbol{v}_{h}^{c})d\boldsymbol{s}.$

Once more, the application of the Cauchy-Schwarz inequality and (5.3) results
in:

$T_{2}\lesssim\left(\sum_{\mathcal{F}_{h}^{0}}h_{\mathcal{F}}\|(2\mu_{h})^{-1/2}[\\![(p_{h}\boldsymbol{I}-2\mu_{h}\boldsymbol{\varepsilon}\boldsymbol{u}_{h})\boldsymbol{n}]\\!]\|^{2}_{0,\mathcal{F}}\right)^{1/2}\\\
\times\left(\sum_{\mathcal{F}_{h}^{0}}h_{\mathcal{F}}^{-1}\|(2\mu_{h})^{1/2}(\boldsymbol{v}-\boldsymbol{v}_{h}^{c})\|^{2}_{0,\mathcal{F}}\right)^{1/2}\lesssim\left(\sum_{K\in\mathcal{T}_{h}}\eta_{\mathcal{F}_{K}}^{2}\right)^{1/2}{\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$

Consolidating the aforementioned estimates establishes the desired result. ∎∎

######  Lemma 5.3.

Let
$(\kappa,(\boldsymbol{u},p))\in\mathbb{R}^{+}\times\mathrm{H}\times\mathrm{Q}$
be the solution of the continuous eigenvalue problem (2.1) and let
$(\kappa_{h},(\boldsymbol{u}_{h},p_{h}))\in\mathbb{R}^{+}\times\mathrm{H}_{h}\times\mathrm{Q}_{h}$
be the solution of the IPDG discretization given by (3.2). Then we have the
following upper bound for the conforming velocity and pressure errors

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}\lesssim\eta+\|(2\mu_{h})^{-1/2}(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h})\|_{0}+\Theta.$

###### Proof.

Using [15, Lemma 2.2], there exists a pair
$(\boldsymbol{v},q)\in\mathrm{H}\setminus\\{0\\}\times\mathrm{Q}$ such that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}^{2}\lesssim{\mathcal{A}}(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h};\boldsymbol{v},q),$

and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}$. Since
$\boldsymbol{u},\boldsymbol{u}_{h}^{c},\boldsymbol{v}\in\mathrm{H}$, we have

${\mathcal{A}}(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h};\boldsymbol{v},q)=(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u}-\boldsymbol{u}_{h}^{c}),\boldsymbol{\varepsilon}(\boldsymbol{v}))-(p-p_{h},\nabla\cdot\boldsymbol{v})\\\
-(q,\nabla\cdot(\boldsymbol{u}-\boldsymbol{u}_{h}^{c}))-(q,1/\lambda(p-p_{h})).$

Based on (2.3), we derive

$\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h};\boldsymbol{v},q)=\kappa(\boldsymbol{u},\boldsymbol{v})-(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u}_{h}^{c}),\boldsymbol{\varepsilon}(\boldsymbol{v}))+(p_{h},\nabla\cdot\boldsymbol{v})+(q,\nabla\cdot\boldsymbol{u}_{h}^{c})+(q,1/\lambda
p_{h})\\\
=\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v})+(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h},\boldsymbol{v})-(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{v}))+(p_{h},\nabla\cdot\boldsymbol{v})\\\
+(q,\nabla\cdot\boldsymbol{u}_{h})+(q,1/\lambda
p_{h})-(q,\nabla\cdot\boldsymbol{u}_{h}^{r})+(2\mu\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}^{r}),\boldsymbol{\varepsilon}(\boldsymbol{v}))-(2(\mu-\mu_{h})\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{v})).$

Consider $\boldsymbol{v}_{h}^{c}=\Pi_{h}\boldsymbol{v}\in\mathrm{H}_{h}^{c}$
as the Scott-Zhang interpolation of $\boldsymbol{v}$. The utilization of (3.4)
results in:

$0=\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})-(2\mu\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}_{h}(\boldsymbol{v}_{h}^{c}))-C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c})+(p_{h},\nabla\cdot\boldsymbol{v}_{h}^{c}).$

Then we have
$\mathcal{A}(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h};\boldsymbol{v},q)=T_{1}+T_{2}+T_{3}+T_{4}$,
where the terms $T_{i}$, $i=1,2,3,4$ are defined by

$\displaystyle T_{1}$
$\displaystyle:=\kappa_{h}(\boldsymbol{u}_{h},\boldsymbol{v}-\boldsymbol{v}_{h}^{c})-(2\mu_{h}\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{v}-\boldsymbol{v}_{h}^{c}))+(p_{h},\nabla\cdot\boldsymbol{v}-\boldsymbol{v}_{h}^{c})$
$\displaystyle+(q,\nabla\cdot\boldsymbol{u}_{h})+(q,1/\lambda p_{h}),$
$\displaystyle T_{2}$
$\displaystyle:=(q,\nabla\cdot\boldsymbol{u}_{h}^{r})+(2\mu\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}^{r}),\boldsymbol{\varepsilon}(\boldsymbol{v})),$
$\displaystyle T_{3}$
$\displaystyle:=C_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h}^{c}),$
$\displaystyle T_{4}$
$\displaystyle:=(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h},\boldsymbol{v})-(2(\mu-\mu_{h})\boldsymbol{\varepsilon}_{h}(\boldsymbol{u}_{h}),\boldsymbol{\varepsilon}(\boldsymbol{v})).$

By employing Lemma 5.2, we obtain
$T_{1}\lesssim\eta{\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}.$
The Cauchy-Schwarz inequality, along with (5.2), demonstrates:

$T_{2}\lesssim{\left|\kern-1.07639pt\left|\boldsymbol{u}_{h}^{r}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}\lesssim\eta{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}.$

By applying Lemma 5.1 to establish the bound of $T_{3}$, we obtain
$T_{3}\lesssim\texttt{a}_{S}^{-1/2}\eta{\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}$.
The application of Cauchy-Schwarz and the Poincaré inequality results in:

$T_{4}\lesssim\left(\|(2\mu_{h})^{-1/2}(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h})\|_{0}+\Theta\right){\left|\kern-1.07639pt\left|\boldsymbol{v}\right|\kern-1.07639pt\right|_{\mathrm{H}(h)}}$

Finally, with the combination of the aforementioned information with the
estimate
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{v},q)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}$, we attain the desired outcome. ∎

######  Theorem 5.4.

Let
$(\kappa,(\boldsymbol{u},p))\in\mathbb{R}^{+}\times\mathrm{H}\times\mathrm{Q}$
be the solution of the continuous eigenvalue problem (2.1) and let
$(\kappa_{h},(\boldsymbol{u}_{h},p_{h}))\in\mathbb{R}^{+}\times\mathrm{H}_{h}\times\mathrm{Q}_{h}$
be the solution of the IPDG discretization given by (3.2). Then the following
a posteriori error bound is derived:

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}\lesssim\eta+\|(2\mu_{h})^{-1/2}(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h})\|_{0}+\Theta,$

where the hidden constant in the expression is independent of the penalty
parameter $\texttt{a}_{S}$, provided that $\texttt{a}_{S}$ is sufficiently
large and satisfies $\texttt{a}_{S}\geq 1$.

###### Proof.

The assertion is readily established by jointly applying Lemma 5.3 and (5.2).
∎

Applying the algebraic identity in conjunction with Theorem 5.4 results in the
following theorem.

######  Theorem 5.5.

If $\boldsymbol{u}\in\mathrm{H}^{2}(\mathcal{T}_{h})$ and
$p\in\mathrm{H}^{1}(\mathcal{T}_{h})$, then the eigenvalue error satisfies

$|\kappa-\kappa_{h}|\lesssim\eta^{2}+\|(2\mu_{h})^{-1/2}(\kappa\boldsymbol{u}-\kappa_{h}\boldsymbol{u}_{h})\|_{0}^{2}+\Theta^{2}.$

###  5.3. Efficiency

The efficiency analysis adheres to conventional arguments grou-nded in bubble
functions [1, 24, 9, 15]. Consequently, we solely provide the theorem
statement for efficiency without presenting its proof.

######  Theorem 5.6.

Assuming $(\kappa,(\boldsymbol{u},p))$ is the solution in the functional
spaces $\mathbb{R}_{+}\times\mathrm{H}\times\mathrm{Q}$ for the eigenvalue
problem(2.1), and $(\kappa_{h},(\boldsymbol{u}_{h},p_{h}))$ is its DG-
approximation obtained through (3.2), with $\eta$ the a posteriori error
estimator defined in (5.1), the following efficiency bound holds

$\eta\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{u}-\boldsymbol{u}_{h}^{c},p-p_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mathrm{H}(h)\times
Q}+\Theta+h.o.t.$

The hidden constant in the expression is independent of the penalty parameter
$\texttt{a}_{S}$, provided that $\texttt{a}_{S}$ is sufficiently large and
satisfies $\texttt{a}_{S}\geq 1$.

##  6\. Numerical experiments

In this section we report some numerical tests in order to assess the
performance of the proposed numerical method in the computation of the
eigenvalues of problem (3.2). These results have been obtained using the
FEniCSx project [3, 23], while the meshes are constructed using GMSH [10].The
convergence rates of the eigenvalues have been obtained with a standard least
square fitting and highly refined meshes.

We denote by $N$ the mesh refinement level, whereas dof denotes the number of
degrees of freedom. We denote by $\kappa_{h_{i}}$ the i-th discrete
eigenfrequency, and by $\widehat{\kappa}_{h_{i}}:=\kappa_{h_{i}}/(1+\nu)$ the
unscaled version of $\kappa_{h_{i}}$.

Hence, we denote the error on the $i$-th eigenvalue by
$\texttt{err}(\widehat{\kappa}_{i})$ with

$\texttt{err}(\widehat{\kappa}_{i}):=|\widehat{\kappa}_{h_{i}}-\widehat{\kappa}_{i}|.$

We remark that computing $\sqrt{\widehat{\kappa}_{h_{i}}}$ or
$\sqrt{\widehat{\kappa}_{i}}$ gives the corresponding discrete and exact
eigenfrequency.

With the aim of assessing the performance of our estimator, we consider
domains with singularities in two and three dimensions in order to observe the
improvement of the convergence rate. On each adaptive iteration, we use the
blue-green marking strategy to refine each $T^{\prime}\in{\mathcal{T}}_{h}$
whose indicator $\eta_{T^{\prime}}$ satisfies

$\eta_{T^{\prime}}\geq 0.5\max\\{\eta_{T}\,:\,T\in{\mathcal{T}}_{h}\\}.$

We define the effectivity indexes with respect to $\eta$ and the eigenvalue
$\widehat{\kappa_{i}}$ by
$\texttt{eff}(\widehat{\kappa_{i}}):=\texttt{err}(\widehat{\kappa_{i}})/\eta^{2}$.

We note that other spectral analysis using DG methods introduce spurious
eigenvalues when the stabilization parameter is not correctly chosen (see for
instance [18, 19]). In the first part of this section, we study the safe range
where the methods is spurious-free by comparing with references values from
the literature. We will analyze the influence stabilization parameter
$\texttt{a}_{S}$. More precisely, the stabilization parameter $\texttt{a}_{S}$
in the bilinear form $a_{h}(\cdot,\cdot)$ (cf. (3.3)) will be chosen
proportionally to the square of the polynomial degree $k$ as
$\texttt{a}_{S}=\texttt{a}k^{2}$ with $\texttt{a}>0$. Also, taking into
account the theoretical convergence rates presented Theorem 4.9, the results
in Table 2 below and the results from the a posteriori analysis, the
experiments will be focused on the SIP method in order to observe the recovery
of the optimal double convergence rate.

###  6.1. Square with bottom clamped boundary condition

We borrow this test from [19]. The computational domain is the unit square
domain $\Omega:=(0,1)^{2}$, with
$\Gamma_{D}:=\\{(x,y)\in\partial\Omega\,:\,y=0\\}$ and
$\Gamma_{N}:=\partial\Omega\backslash\Gamma_{D}$. The Poisson ratio is taken
to be $\nu=0.35$ and $E=1$. The boundary conditions are
$\boldsymbol{u}=\boldsymbol{0}$ on $\Gamma_{D}$ and
$\left(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})-p\boldsymbol{I}\right)\boldsymbol{n}=\boldsymbol{0}$
on $\Gamma_{N}$. The test consider a uniform mesh with $N=8$. The spurious
study for this case have been previously studied in [19] in the context of a
mixed discontinuous Galerkin scheme.

In Table 1 we report the results of computing the first 10 lowest eigenvalues
for different values of the stabilization parameter and different polynomial
orders. The spurious eigenvalues are highlighted in bold numbers. It notes
that the higher the penalization value is, the less spurious eigenvalues
appear. Also, for $k=3$ the spurious appearance is considerably low. For the
NIP and IIP methods we observed a similar behavior as those of [19], so we
refer to that study for further details. The experiment ends with Table 2, in
which we study the convergence of the first two eigenvalues. Since we have a
mixed boundary condition setting, then there is a singularity in the corners
where Dirichlet changes to Neumann boundary. This is reflected in the
convergence rate, where we observe a $\mathcal{O}(h^{r}),\,r\geq 1.2$ for all
the values of $\nu$ selected.

Table 1. Test 6.1. Computed eigenfrequencies $\sqrt{\kappa_{h,i}}$ by using the SIP method ($\varepsilon=1$) for polynomial degrees $k=1,2,3$ and mesh level $N=8$ on the unit square domain. The Poisson ratio for this case is $\nu=0.35$. $k$ |  Ref. [19] |  $\texttt{a}=$1/4 |  $\texttt{a}=$1/2 |  $\texttt{a}=$1 |  $\texttt{a}=$2 |  $\texttt{a}=$4 |  $\texttt{a}=$8
---|---|---|---|---|---|---|---
1 | 0.6808 | 0.6555 | 0.6612 | 0.6631 | 0.6738 | 0.6809 | 0.6881
1.6993 | $\bf{1.0131}$ | 1.6781 | 1.6914 | $\bf{0.7873}$ | 1.6991 | 1.7065
1.8222 | 1.6719 | 1.7973 | 1.7771 | 1.6887 | 1.8310 | 1.8426
2.9477 | 1.7880 | $\bf{2.0082}$ | $\bf{1.8554}$ | $\bf{1.8162}$ | 2.9702 | 3.0085
3.0181 | $\bf{2.0217}$ | $\bf{2.0735}$ | $\bf{2.7417}$ | 1.8440 | 3.0434 | 3.0720
3.4433 | 2.8190 | 2.8669 | 2.9291 | $\bf{1.9966}$ | 3.4740 | 3.5052
4.1418 | 2.9139 | 2.9480 | 2.9524 | 2.9651 | 4.2061 | 4.2830
4.6312 | $\bf{3.2805}$ | 3.3931 | $\bf{3.3762}$ | 3.0285 | 4.6810 | 4.7681
4.7616 | 3.3557 | 3.8363 | 3.4530 | 3.4592 | 4.8488 | 4.9367
| 4.7887 | $\bf{3.4927}$ | 4.3341 | $\bf{3.5020}$ | $\bf{3.5864}$ | 4.8720 | 4.9775
2 | 0.6808 | 0.6776 | 0.6768 | 0.6701 | 0.6803 | 0.6816 | 0.6825
1.6993 | 1.6963 | $\bf{0.7360}$ | 1.6839 | 1.6985 | 1.7002 | 1.7013
1.8222 | 1.8213 | 1.6974 | 1.8217 | 1.8225 | 1.8226 | 1.8229
2.9477 | $\bf{2.7325}$ | 1.8208 | $\bf{2.5678}$ | 2.9481 | 2.9485 | 2.9489
3.0181 | 2.9482 | $\bf{2.2148}$ | 2.9309 | 3.0176 | 3.0203 | 3.0225
3.4433 | 3.0332 | 2.9466 | 3.0554 | 3.4435 | 3.4444 | 3.4450
4.1418 | 3.4410 | 3.0229 | $\bf{3.0947}$ | 4.1436 | 4.1459 | 4.1479
4.6312 | $\bf{3.8448}$ | 3.4425 | 3.4309 | 4.6336 | 4.6364 | 4.6395
4.7616 | $\bf{4.0832}$ | 4.1361 | $\bf{3.9605}$ | 4.7635 | 4.7665 | 4.7686
| 4.7887 | 4.1458 | 4.6241 | 4.1301 | 4.7899 | 4.7954 | 4.7992
3 | 0.6808 | 0.6760 | 0.6724 | 0.6808 | 0.6808 | 0.6814 | 0.6817
1.6993 | 1.6972 | 1.6929 | 1.6989 | 1.6993 | 1.6999 | 1.7003
1.8222 | 1.8220 | 1.8221 | 1.8222 | 1.8222 | 1.8223 | 1.8224
2.9477 | $\bf{2.0530}$ | 2.9477 | 2.9477 | 2.9477 | 2.9477 | 2.9477
3.0181 | 2.9476 | 2.9990 | 3.0179 | 3.0181 | 3.0191 | 3.0198
3.4433 | 3.0176 | 3.4420 | $\bf{3.2003}$ | 3.4433 | 3.4435 | 3.4436
4.1418 | 3.4427 | 4.1407 | 3.4434 | 4.1419 | 4.1422 | 4.1424
4.6312 | 4.1412 | 4.6258 | 4.1418 | 4.6313 | 4.6320 | 4.6325
4.7616 | 4.6289 | $\bf{4.6365}$ | 4.6312 | 4.7617 | 4.7619 | 4.7619
| 4.7887 | 4.7612 | 4.7614 | 4.7618 | 4.7887 | 4.7898 | 4.7904
Table 2. Test 6.1. Convergence behavior of the first three lowest computed eigenvalues by using the SIP method ($\varepsilon=1$) on the unit square domain . The stabilization parameter is set to be $\texttt{a}=10$. $\nu$ | $N=20$ | $N=30$ | $N=40$ | $N=50$ | Order | $\widehat{\kappa}_{\text{extr}}$
---|---|---|---|---|---|---
$k=1$
0.35 | 0.6832 | 0.6821 | 0.6817 | 0.6815 | 1.50 | 0.6809
1.7018 | 1.7008 | 1.7003 | 1.7000 | 1.41 | 1.6993
0.49 | 0.7033 | 0.7018 | 0.7011 | 0.7008 | 1.27 | 0.6996
1.8439 | 1.8412 | 1.8399 | 1.8392 | 1.33 | 1.8373
0.50 | 0.7056 | 0.7040 | 0.7033 | 0.7029 | 1.26 | 0.7016
1.8557 | 1.8528 | 1.8515 | 1.8508 | 1.33 | 1.8487

###  6.2. Clamped square with different materials

In this test we study the behavior of the numerical scheme when we consider an
heterogeneous media along the domain. For this, we consider two different
materials on the unit square domain. The boundary conditions for this
experiments are

$\boldsymbol{u}=\boldsymbol{0}\text{ on
}\Gamma_{D}:=(0,y)\cup(1,y),\,y\in[0,1],$

while
$\left(2\mu\boldsymbol{\varepsilon}(\boldsymbol{u})-p\boldsymbol{I}\right)\boldsymbol{n}=\boldsymbol{0}$
is assumed on the rest of the boundary. The computational domain is decomposed
such that $\Omega:=\Omega_{1}\cup\Omega_{2}$, where
$\Omega_{1}:=(0,0)\times(0,1/2)$ and $\Omega_{2}:=(0,1/2)\times(1,1)$. A
sample mesh for this decomposition is portrayed in Figure 1. The Young modulus
and density parameters are then given by

(6.1) $E(x):=\left\\{\begin{aligned} &7.72\cdot 10^{10}\text{Pa}\text{ in
}\Omega_{1},\\\ &1.10\cdot 10^{11}\text{Pa}\text{ in
}\Omega_{2},\end{aligned},\right.\qquad\rho(x)=\left\\{\begin{aligned}
&19300\,\text{kg/m}^{3}\text{ in }\Omega_{1},\\\ &8850\,\text{kg/m}^{3}\text{
in }\Omega_{2},\end{aligned}\right.$

which consist of gold in $\Omega_{1}$, and copper on $\Omega_{2}$. By the
results from the previous example, it is enough to consider $k=1$ to study the
convergence on the a priori case.

The convergence with uniform meshes is presented in Table 3. Here, we observe
a clear suboptimal rate of convergence for the four lowest computed
eigenvalues. By observing Figure 1 we note that the heterogeneity of the
medium yields to a considerable difference in the displacement field for the
first eigenmode with $\nu=0.35$ and $\nu=0.5$. Also, the pressure contour
shows, along with the corner singularities, two zones of high gradients near
the intersection of the material discontinuity and the boundary domain.

When studying the a posteriori refinements for this case, we consider $k=1,2$
and a total of $15$ adaptive iterations. Figure 2 depicts several intermediate
meshes for the chosen polynomials degrees. We note that the refinements
concentrate along the singularities and the number of marked elements is lower
for the higher degree. The adaptive refinement also makes evident the
difference between the materials on the domain. Those results are followed by
the error curves, presented in Figure 3, where we observe an experimental rate
of $\mathcal{O}(h^{k+1})$. The reliability and the efficiency of the estimator
are evidenced on Figure 4, where a proper boundedness of $\texttt{eff}(\cdot)$
and a decay behaving like $\mathcal{O}(h^{k+1})$ for $\eta^{2}$ are observed.

Figure 1. Test 6.2. Left: Sample mesh for the subdivided computational domain with two materials. The next two figures correspond to the warped domain using the lowest mode $\boldsymbol{u}_{1,h}$ for $\nu=0.35$ and $\nu=0.5$, respectively. The right figure correspond to the first pressure mode, which is similar for $\nu=0.35$ and $\nu=0.5$.  Table 3. Test 6.2. Convergence behavior of the first fourth lowest computed eigenvalues by using the SIP method ($\varepsilon=1$) for $k=1$ on the unit square domain with different materials. The stabilization parameter is set to be $\texttt{a}=10$. $\nu$ | $N=8$ | $N=16$ | $N=32$ | $N=64$ | Order | $\widehat{\kappa}_{\text{extr}}$
---|---|---|---|---|---|---
0.35 | 4506.7208 | 4457.1445 | 4439.2239 | 4433.2168 | 1.50 | 4429.6821
7546.2296 | 7454.7702 | 7421.5220 | 7410.1874 | 1.49 | 7403.5352
| 8043.1998 | 7864.1790 | 7812.3865 | 7798.2104 | 1.82 | 7792.2188
| 10533.2314 | 10301.4824 | 10223.0867 | 10200.2366 | 1.62 | 10187.2085
0.49 | 4462.9780 | 4419.4198 | 4401.2061 | 4394.2151 | 1.30 | 4389.1357
7623.5227 | 7454.1769 | 7396.1249 | 7376.9013 | 1.57 | 7366.7357
| 8692.8982 | 8521.8508 | 8468.1473 | 8451.5251 | 1.69 | 8443.9533
| 11173.4976 | 10698.5980 | 10551.1393 | 10507.2968 | 1.74 | 10488.3190
0.50 | 4464.8940 | 4421.5136 | 4403.1190 | 4395.9793 | 1.28 | 4390.6467
7611.7557 | 7442.4339 | 7384.3312 | 7365.0352 | 1.57 | 7354.9169
| 8771.4492 | 8591.9091 | 8535.1211 | 8517.4054 | 1.68 | 8509.2642
| 11166.2453 | 10690.8606 | 10543.1916 | 10499.1463 | 1.73 | 10479.2863

$\nu=0.35,k=1,iter=10$

$\nu=0.35,k=2,iter=10$

$\nu=0.35,k=1,iter=15$

$\nu=0.35,k=2,iter=15$

$\nu=0.50,k=1,iter=10$

$\nu=0.50,k=2,iter=10$

$\nu=0.50,k=1,iter=15$

$\nu=0.50,k=2,iter=15$

Figure 2. Test 6.2. Intermediate meshes of the square domain with different
materials obtained with the adaptive algorithm and different values of $\nu$,
with $k=1,2$.  Figure 3. Test 6.2. Error curves obtained from the adaptive
algorithm in the square domain compared with the lines
$\mathcal{O}(\texttt{dof}^{-1})$ and $\mathcal{O}(\texttt{dof}^{-2})$.

Figure 4. Test 6.2. Estimator and efficiency curves obtained from the adaptive
algorithm in the square domain with different values of $\nu$ and $k=1,2$.

###  6.3. The 3D L-shape domain

We now extend the results of our numerical scheme to three dimensions. In this
experiment we test the estimator in the case of a domain with a dihedral
singularity. The domain under consideration is the L-shaped domain defined by

$\Omega:=(-1,1)\times(-1,1)\times(-1,0)\backslash\bigg{(}(-1,0)\times(-1,0)\times(-1,0)\bigg{)}.$

For simplicity, the lowest order $k=1$ and $8$ adaptive iterations are
considered in this experiment. Fully fixed boundary conditions are assumed
across $\partial\Omega$. The domain is divided in two subdomains $\Omega_{1}$
and $\Omega_{2}$ with different materials, where

$\displaystyle\Omega_{1}:=(-1,1)\times(-1,1)\times(-1,-0.5)\backslash\bigg{(}(-1,0)\times(-1,0)\times(-1,-0.5)\bigg{)},$
$\displaystyle\Omega_{2}:=(-1,1)\times(-1,1)\times(-0.5,0)\backslash\bigg{(}(-1,0)\times(-1,0)\times(-0.5,0)\bigg{)}.$

The physical parameters are similar to those of Section 6.2, namely, the Young
modulus and density are taken similar to (6.1). A sample of the computational
domain is depicted in Figure 5. Here, the initial mesh is such that $h\approx
1/2$. Because of the singularities presented in this domain, the first
eigenvalue is singular for each value of $\nu$. The extrapolated values,
computed using highly refined meshes are given in Table 4.

We present the first eigenmodes on Figure 5 consisting on the displacement
field streamlines. We observe that high values of the displacement magnitude
are concentrated on $\Omega_{1}$ (bottom domain) because of the difference in
the density of the materials. Different adaptive meshes for this domain are
presented in Figure 6, where we observe that the estimator is capable of
detecting the singular behavior of the pressure for each Poisson ratio and
refine near that zone. Also, the algorithm marks more elements near the line
singularity belonging to $\overline{\Omega}_{1}$, as expected.

Let us finish this experiment with Figure 3. The figure represent the error
curves and estimator values obtained with adaptive refinements, where an
optimal rate is observed. For the 3D case, a rate of
$\mathcal{O}(h^{2})\approx\mathcal{O}(\texttt{dof}^{-0.66})$ was expected.

Figure 5. Test 6.3. Initial mesh for the subdivided computational domain $\Omega:=\Omega_{1}\cup\Omega_{2}$ (left) together with a comparison of the first displacement eigenmodes for $\nu=0.35$ (middle) and $\nu=0.5$ (right) at the last adaptive iteration. Table 4. Test 6.3. Reference lowest computed eigenvalues in the 3D Lshape domain. $\nu$ | $\widehat{\kappa}_{1}$
---|---
0.35 | $7187.76090678178$
0.49 | $8591.22953946104$
0.5 | $8600.14945418864$

$\nu=0.35,iter=5$

$\nu=0.35,iter=7$

$\nu=0.35,iter=8$

$\nu=0.50,iter=5$

$\nu=0.50,iter=7$

$\nu=0.50,iter=8$

Figure 6. Test 6.3. Intermediate meshes of the square domain with
heterogeneous media obtained with the adaptive algorithm and different values
of $\nu$. Figure 7. Test 6.3. Error curves and estimator values obtained from
the adaptive algorithm in the Lshape domain compared with the line
$\mathcal{O}(\texttt{dof}^{-0.66})$.

###  6.4. Robustness in 2D and 3D geometries

This final experiment aims to assess the robustness of the proposed DG
methods. For the two-dimensional case, we consider again the unit square
domain $\Omega:=(0,1)^{2}$ and the unit cube domain $\Omega:=(0,1)^{3}$ is
considered for the three-dimensional case. The square and the cube are
considered to be clamped at the bottom and free of stress on the rest of the
facets. In the cube domain, we consider the bottom as the set of facets such
that $y=0$. The extrapolated values for this experiment are calculated from a
reference value for each value of $\nu$. More precisely, in each test we have
the extrapolated values in Table 5 for $E=10^{j},j=1,2,3...$.

Table 5. Test 6.4. Reference eigenvalues for different values of $\nu$ and $E=10^{j},j=1,2,3...$ . | Square domain | Cube domain
---|---|---
$\nu$ | $\widehat{\kappa}_{1}$ | $\widehat{\kappa}_{1}$
0.35 | $0.46355423498481496\cdot E$ | $0.444317882233217\cdot E$
0.49 | $0.48938358373431\cdot E$ | $0.44833605908518\cdot E$
0.5 | $0.492273855811713\cdot E$ | $0.44915941661888\cdot E$

A stabilization parameter $\texttt{a}=10$ is chosen for all the choices of $E$
and $\nu$. From the results displayed in Table 6 we observe that a almost a
constant value for $\texttt{eff}(\cdot)$ is attained for the selected values
of $E$ in two and three dimensions. We recall that similar results were
obtained when considering Young modulus $E\in\\{10^{-1},10^{-2},10^{-4}\\}$.
Hence, the proposed estimator is robust with respect to the physical
parameters.

Table 6. Test 6.4. Efficiency indexes for different values of $\nu$ and $E$ on the unit square and unit cube domain when using the SIP method ($\varepsilon=1$). | $\nu=0.35$ | $\nu=0.49$ | $\nu=0.50$
---|---|---|---
Unit square domain $\Omega=(0,1)^{2}$
dof | $\texttt{eff}(\widehat{\kappa}_{1})$ | $\texttt{eff}(\widehat{\kappa}_{1})$ | $\texttt{eff}(\widehat{\kappa}_{1})$
126 | 9.01e-02 | 9.01e-02 | 9.01e-02 | 8.17e-02 | 8.19e-02 | 8.19e-02 | 8.13e-02 | 8.15e-02 | 8.15e-02
350 | 9.60e-02 | 9.60e-02 | 9.60e-02 | 9.16e-02 | 9.18e-02 | 9.18e-02 | 9.11e-02 | 9.14e-02 | 9.14e-02
1134 | 1.02e-01 | 1.02e-01 | 1.02e-01 | 9.92e-02 | 9.94e-02 | 9.94e-02 | 9.86e-02 | 9.89e-02 | 9.89e-02
4046 | 1.08e-01 | 1.08e-01 | 1.08e-01 | 1.04e-01 | 1.04e-01 | 1.04e-01 | 1.03e-01 | 1.04e-01 | 1.04e-01
15246 | 1.12e-01 | 1.12e-01 | 1.12e-01 | 1.07e-01 | 1.07e-01 | 1.07e-01 | 1.05e-01 | 1.06e-01 | 1.06e-01
Unit cube domain $\Omega=(0,1)^{3}$
dof | $\texttt{eff}(\widehat{\kappa}_{1})$ | $\texttt{eff}(\widehat{\kappa}_{1})$ | $\texttt{eff}(\widehat{\kappa}_{1})$
9750 | 4.53e-01 | 4.54e-01 | 4.54e-01 | 2.21e-01 | 2.21e-01 | 2.21e-01 | 1.80e-01 | 1.80e-01 | 1.80e-01
78000 | 8.05e-01 | 8.06e-01 | 8.06e-01 | 2.33e-01 | 2.33e-01 | 2.33e-01 | 1.43e-01 | 1.43e-01 | 1.43e-01
263250 | 1.13e+00 | 1.13e+00 | 1.13e+00 | 2.60e-01 | 2.60e-01 | 2.60e-01 | 1.27e-01 | 1.28e-01 | 1.28e-01
624000 | 1.43e+00 | 1.44e+00 | 1.44e+00 | 2.93e-01 | 2.94e-01 | 2.94e-01 | 1.19e-01 | 1.19e-01 | 1.19e-01
1218750 | 1.72e+00 | 1.72e+00 | 1.72e+00 | 3.29e-01 | 3.30e-01 | 3.30e-01 | 1.13e-01 | 1.14e-01 | 1.14e-01
$E$ | $10$ | $10^{2}$ | $10^{4}$ | $10$ | $10^{2}$ | $10^{4}$ | $10$ | $10^{2}$ | $10^{4}$

## References

  * [1] Mark Ainsworth and J. Tinsley Oden, _A posteriori error estimation in finite element analysis_ , Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000. MR 1885308
  * [2] Paola F. Antonietti, Annalisa Buffa, and Ilaria Perugia, _Discontinuous Galerkin approximation of the Laplace eigenproblem_ , Comput. Methods Appl. Mech. Engrg. 195 (2006), no. 25-28, 3483–3503. MR 2220929
  * [3] Igor A Barrata, Joseph P Dean, Jωrgen S Dokken, Michal HABERA, Jack HALE, Chris Richardson, Marie E Rognes, Matthew W Scroggs, Nathan Sime, and Garth N Wells, _DOLFINx: The next generation fenics problem solving environment_ , (2023).
  * [4] Annalisa Buffa and Ilaria Perugia, _Discontinuous Galerkin approximation of the Maxwell eigenproblem_ , SIAM J. Numer. Anal. 44 (2006), no. 5, 2198–2226. MR 2263045
  * [5] Bernardo Cockburn, Guido Kanschat, and Dominik Schötzau, _A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations_ , J. Sci. Comput. 31 (2007), no. 1-2, 61–73.
  * [6] Jean Descloux, Nabil Nassif, and Jacques Rappaz, _On spectral approximation. I. The problem of convergence_ , RAIRO Anal. Numér. 12 (1978), no. 2, 97–112, iii. MR 483400
  * [7] by same author, _On spectral approximation. II. Error estimates for the Galerkin method_ , RAIRO Anal. Numér. 12 (1978), no. 2, 113–119, iii. MR 483401
  * [8] E. B. Fabes, C. E. Kenig, and G. C. Verchota, _The Dirichlet problem for the Stokes system on Lipschitz domains_ , Duke Math. J. 57 (1988), no. 3, 769–793. MR 975121
  * [9] Joscha Gedicke and Arbaz Khan, _Divergence-conforming discontinuous Galerkin finite elements for Stokes eigenvalue problems_ , Numer. Math. 144 (2020), no. 3, 585–614. MR 4071826
  * [10] Christophe Geuzaine and Jean-François Remacle, _Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities_ , International journal for numerical methods in engineering 79 (2009), no. 11, 1309–1331.
  * [11] P Grisvard, _Problemes aux limites dans les polygones. mode d’emploi_ , Bulletin de la Direction des Etudes et Recherches Series C Mathematiques, Informatique 1 (1986), 21–59.
  * [12] Peter Hansbo and Mats G. Larson, _Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’s method_ , Comput. Methods Appl. Mech. Engrg. 191 (2002), no. 17-18, 1895–1908. MR 1886000
  * [13] Paul Houston, Dominik Schötzau, and Thomas P. Wihler, _Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Stokes problem_ , J. Sci. Comput. 22/23 (2005), 347–370.
  * [14] Guido Kanschat and Dominik Schötzau, _Energy norm a posteriori error estimation for divergence-free discontinuous Galerkin approximations of the Navier-Stokes equations_ , Internat. J. Numer. Methods Fluids 57 (2008), no. 9, 1093–1113.
  * [15] Arbaz Khan, Felipe Lepe, David Mora, and Jesus Vellojin, _Finite element analysis of the nearly incompressible linear elasticity eigenvalue problem with variable coefficients_ , 2023.
  * [16] Arbaz Khan, Catherine E. Powell, and David J. Silvester, _Robust a posteriori error estimators for mixed approximation of nearly incompressible elasticity_ , Internat. J. Numer. Methods Engrg. 119 (2019), no. 1, 18–37. MR 3973678
  * [17] Felipe Lepe, _Interior penalty discontinuous Galerkin methods for the velocity-pressure formulation of the Stokes spectral problem_ , Adv. Comput. Math. 49 (2023), no. 4, Paper No. 61, 31. MR 4623018
  * [18] Felipe Lepe, Salim Meddahi, David Mora, and Rodolfo Rodríguez, _Mixed discontinuous Galerkin approximation of the elasticity eigenproblem_ , Numer. Math. 142 (2019), no. 3, 749–786. MR 3962898
  * [19] Felipe Lepe and David Mora, _Symmetric and nonsymmetric discontinuous Galerkin methods for a pseudostress formulation of the Stokes spectral problem_ , SIAM J. Sci. Comput. 42 (2020), no. 2, A698–A722. MR 4077220
  * [20] Felipe Lepe, David Mora, and Jesus Vellojin, _Discontinuous Galerkin methods for the acoustic vibration problem_ , J. Comput. Appl. Math. 441 (2024), Paper No. 115700. MR 4673997
  * [21] Felipe Lepe, Gonzalo Rivera, and Jesus Vellojin, _Mixed methods for the velocity-pressure-pseudostress formulation of the Stokes eigenvalue problem_ , SIAM J. Sci. Comput. 44 (2022), no. 3, A1358–A1380. MR 4430561
  * [22] L. Ridgway Scott and Shangyou Zhang, _Finite element interpolation of nonsmooth functions satisfying boundary conditions_ , Math. Comp. 54 (1990), no. 190, 483–493. MR 1011446
  * [23] Matthew W Scroggs, Igor A Baratta, Chris N Richardson, and Garth N Wells, _Basix: a runtime finite element basis evaluation library_ , Journal of Open Source Software 7 (2022), no. 73, 3982.
  * [24] Rüdiger Verfürth, _A posteriori error estimation techniques for finite element methods_ , Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2013. MR 3059294