The symmetric
C0 interior penalty method
is one of the most popular
discontinuous Galerkin methods
for the biharmonic equation.
This paper introduces an automatic selection
of the involved stability parameter
in terms of the geometry
of the underlying triangulation.
It avoids any extra input variables
and guarantees a stable discretization
with stability $1/2$.
The practical realization
of this parameter-free scheme
in less than 120 lines of MATLAB
is described and explained in detail.
Numerical experiments show
stability constants smaller than 1
and so over-stabilization is excluded.
The approach is documented
for the lowest-order scheme
in 2D for triangles,
but the methodology behind can be generalized
for automatic local parameter decisions.

C0 interior penalty method,
discontinuous Galerkin method,
biharmonic equation,

65N12, 65N15, 65N30, 65N50, 65Y20

§ INTRODUCTION

Let $\Omega \subset \R^2$ be a bounded Lipschitz domain
with polygonal boundary $\partial \Omega$.
For a given source term $f \in L^2(\Omega)$,
suppose $u \in V \coloneqq H^2_0(\Omega)$ solves
the biharmonic equation $\Delta^2 u = f$.
The weak formulation seeks \(u \in V\) with

\begin{align}\label{eq:weakform}
\coloneqq
\int_\Omega \D^2 u : \D^2 v \dx
\int_\Omega f v \dx
\quad \text{for all } v \in V.
\end{align}

The Riesz representation theorem
applies in the Hilbert space $(V,a)$ and proves
well-posedness of the formulation (<ref>).
Elliptic regularity theory verifies that
$f \in L^2(\Omega)$ implies $u \in H^{2+\alpha}(\Omega) \cap H^2_0(\Omega)$
[1, 5, 19, 21].
The pure Dirichlet boundary conditions
in the model example
lead to $\alpha > 1/2$ and
allow for a control
of the traces $\D^2 u$ in the jump terms.
The conforming finite element method for the plate problem
originates from the work of Argyris [2]
with a quintic polynomial ansatz space $P_5(T)$
and $21$ degrees of freedom on each triangle $T$.
The practical application is less prominent
due to the higher computational efforts,
classical nonconforming FEMs
[11, 13, 14, 18]
and discontinuous Galerkin schemes [3, 10, 11, 16, 24]
appear more popular in the literature.
Discontinuous Galerkin methods allow for discontinuities
across triangle interfaces leading to more flexible schemes
compared to their classical nonconforming pendants.

The key idea of discontinuous Galerkin methods is
the replacement of the space \(V\)
in the weak formulation (<ref>)
by the piecewise Sobolev space

\begin{equation}
\label{eq:piecewise_Sobolev}
\coloneqq
\big\{
v \in L^2(\Omega)
\forall T \in \T,\:
\in H^2(T) \coloneqq H^2(\Int(T))
\big\}
\end{equation}

with respect to a shape-regular triangulation \(\T\)
of the domain $\Omega$ into closed triangles.

For \(v \in H^2(\T)\),
the piecewise application of the distributional derivatives
leads to the Hessian \(\D^2_\PW v \in L^2(\Omega; \mathbb{S})\),
\((\D^2_\PW v)\vert_T \coloneqq \D^2 (v\vert_T)\)
with values in the space \(\mathbb{S} \subset \R^{2 \times 2}\)
of symmetric 2-by-2 matrices.

The piecewise Sobolev functions \(v \in H^2(\T)\),
allow for the evaluation of
averages \(\langle v\rangle_E\) and
jumps \([v]_E\) across an edge \(E \in \E\).
Each interior edge \(E \in \E(\Omega)\)
of length \(h_E \coloneqq \vert E \vert\)
is the common edge of exactly two triangles \(T_+, T_- \in T\), written
\(E = \partial T_+ \cap \partial T_-\).

\begin{equation}
\label{eq:jump}
\coloneqq
(v \vert_{T_+})\vert_E - (v \vert_{T_-})\vert_E
\quad\text{and}\quad
\langle v \rangle_E
\coloneqq
\frac12 (v \vert_{T_+})\vert_E + \frac12 (v \vert_{T_-})\vert_E.
\end{equation}

The unit normal vector \(\nu_E\) is oriented such that
\(\nu_E \cdot \nu_{T_\pm}\vert_E = \pm 1\)
for the outward unit normal vectors \(\nu_{T_\pm}\) of \(T_\pm\). For each boundary edge \(E \in \E(\partial\Omega)\),
let \(T_+ \in \T\) denote the unique triangle with edge
\(E \in \E(T_+)\) and set

\(
\coloneqq
\langle v \rangle_E
\coloneqq
(v \vert_{T_+})\vert_E.
\)

Analogous definitions apply for vector- or matrix-valued polynomials.
For further details on the notation,
the reader is referred to the following Section <ref>.

The symmetric C0 interior penalty method
seeks a function
\(
u_\IP \in V_h
\coloneqq
\coloneqq
P_2(\T) \cap H^1_0(\Omega)
\)
in the quadratic Lagrange finite element space with

\begin{equation}
\label{eq:discrete_form}
A_h(u_\IP, v_\IP)
\int_\Omega f v_\IP \dx
\quad \text{for all } v_\IP \in V_h.
\end{equation}

The three contributions to
the bilinear form $A_h: V_h \times V_h \to \R$,

\begin{equation}
\label{eq:definition_Ah}
A_h(u_{\IP}, v_{\IP})
\coloneqq
- (\mathcal{J}(u_\IP,v_\IP)
+ \mathcal{J}(v_\IP, u_\IP))
+ c_\IP(u_\IP, v_\IP)
\end{equation}

for all $u_\IP, v_\IP \in V_h$ originate from the triangle-wise application
of the integration by parts
in the derivation of the weak formulation (<ref>)
This leads to
a piecewise version of the energy scalar product
\(a_\PW: V_h \times V_h \to \R\)
and a jump term \(\mathcal{J}: V_h \times V_h \to \R\),

\begin{equation}
\label{eq:stiffness_and_jump}
\begin{split}
\sum_{T \in \mathcal{T}}
\int_T \D^2 u_\IP: \D^2 v_\IP \dx, \\
\mathcal{J}(u_\IP,v_\IP)
\sum_{E \in \mathcal{E}}
\int_E
\langle \D^2_\PW u_\IP \nu_E \rangle_E
\cdot [\nabla v_\IP]_E \dx
\end{split}
\end{equation}

for all \(u_\IP, v_\IP \in V_h\).

The \(H^2\) conformity is imposed weakly by the
additional penalty term $c_\IP: V_h \times V_h \to \R$,

\begin{equation}
\label{eq:penalty_term}
\coloneqq
\sum_{E \in \mathcal{E}}
\frac{\sigma_{\IP,E}}{h_E}
\int_E
[\nabla u_\IP \cdot \nu_E]_E
[\nabla v _\IP \cdot \nu_E]_E
\ds \text{ for all } u_\IP, v_\IP \in V_h,
\end{equation}

with some edge-dependent parameter
\(\sigma_{\IP,E} > 0\) for each \(E \in \E\).

While the boundedness of \(A_h\)
follows from standard arguments,
its coercivity is typically subject to the assumption
of a sufficiently large penalty parameter
$\sigma_{\IP,E}$ [10, 12].
The contribution of this paper is twofold.
First, a novel stability analysis
introduces an automated mesh-dependent selection
of \(\sigma_{\IP, E}\) with some \(a > 1\) by

\begin{align}
\label{eq:definition_sigma}
\sigma_{\IP,E}
\coloneqq
\begin{cases}
\displaystyle
\frac{3a h_E^2}{4}
\bigg(
\frac{1}{\vert T_+\vert} + \frac{1}{\vert T_- \vert}
\bigg)
\text{if }
E = \partial T_+ \cap \partial T_- \in
\E(\Omega),\medskip\\
\displaystyle
\frac{3a h_E^2}{\vert T_+\vert}
\text{if }
E \in \E(T_+) \cap \E(\partial\Omega).
\end{cases}
\end{align}

Every choice of \(a > 1\) leads to guaranteed stability
in the main result Theorem <ref>
in Section <ref> below.

Second, the numerical realization of the
C0 interior penalty method
is provided in great detail.
The implementation bases on the evaluation of derivatives
of discrete functions \(v_h \in V_h\) presented in
Section <ref> below.
The Sections <ref>–<ref>
cover the computation of the local contributions
to the bilinear forms \(a_\PW\), \(\mathcal{J}\),
and \(c_\IP\) from (<ref>)–(<ref>)
and Section <ref>
describes the assembling of the linear system of equations.

Particular focus is put on the procedure for
evaluating the jumps \([\nabla v_h \cdot \nu_E]_E\)
and the averages \(\langle \D^2_\PW v_h \nu_E \rangle_E\)
in \(\mathcal{J}\) and \(c_\IP\)
in Section <ref> below.

The automatic computation of the penalty parameter \(\sigma_{\IP, E}\)
solely requires four lines of MATLAB and avoids
any manual selection by the user.

The presentation aims to be easily accessible
and serves as a documentation
of the complete implementation
given in the Appendix <ref>.

Numerical experiments with the implementation at hand
investigate the stability of the scheme with
the suggested automatic penalty selection
in Section <ref>.
The coercivity constants are computed numerically
and compared with the theoretically established value.

the performance of an adaptive mesh-refinement algorithm
is examined. The appendix includes the complete implementation
as well as a realization of an a posteriori
error estimator.

§ NOTATION

Standard notation of Lebesgue and Sobolev spaces, their norms,
and $L^2$ scalar products applies throughout the paper.

Let $\T$ denote a shape regular triangulation
of the polygonal Lipschitz domain $\Omega$ into closed triangles
and \(\V\) (resp. $\V(\Omega)$ or $\V(\partial \Omega)$)
the set of all (resp. interior or boundary) vertices [7, 9].
Let \(\E\) (resp. $\E(\Omega)$ or $\E(\partial \Omega)$)
be the set of all (resp. interior or boundary) edges.

For each triangle \(T \in \T\) of area \(\vert T \vert\),
let \(\V(T)\) denote the set of its three vertices and
\(\E(T)\) the set of its three edges.
Abbreviate the edge patch
by \(\omega_E \coloneqq \Int(T_+ \cup T_-) \subseteq \Omega\)
for an interior edge
\(E = \partial T_+ \cap \partial T_- \in \E(\Omega)\)
and by
\(\omega_E \coloneqq \Int(T_+) \subseteq \Omega\)
for a boundary edge \(E \in \E(T_+) \cap \E(\partial\Omega)\).

Define the space of piecewise polynomials

\[
\coloneqq
\{
v \in L^2(\Omega)
\;:\;
v\vert_{T} \in P_k(T) \ \text{for all} \ T \in \T
\}
\]

of total degree at most $k \in \N_0$.

The notation $A \lesssim B$ abbreviates $A \le C B$ for a positive,
generic constant $C$, solely depending on the domain $\Omega$,
the polynomial degree, and the shape regularity of the triangulation $\T$;
$A \approx B$ abbreviates $A \lesssim B \lesssim A$.

For two indices \(j, k \in \N\),
let \(\delta_{jk} \in \{0, 1\}\) denote the Kronecker symbol
defined by \(\delta_{jk} \coloneqq 1\) if and only if \(j = k\).
The enclosing single bars \(\vert \bullet \vert\)
apply context-sensitively and
denote not only the modulus of real numbers,
the Euclidian norm of vectors in \(\R^2\),
but also the cardinality of finite sets,
the area of two-dimensional Lebesgue sets,
and the length of edges.

§ STABILITY

This section develops a novel stability analysis
with a penalty parameter $\sigma_{\IP,E}$
in the bilinear form \(c_\IP\) from (<ref>).

Let the discrete space \(V_h\)
be equipped with the mesh-dependent norm
$\Vert \bullet \Vert_h$ defined,
for \(v_\IP \in V_h\), by

\[
\Vert v_\IP \Vert_h^2
\coloneqq
\vvvert v_\IP \vvvert_{\PW}^2 + c_\IP(v_\IP, v_\IP).
\]

with the piecewise semi-norm
\(
\vvvert \bullet \vvvert_{\PW}
\coloneqq
\vert \bullet \vert_{H^2(\T)}
\coloneqq
\Vert \D^2_\PW \bullet \Vert_{L^2(\Omega)}
\).

The parameter \(\sigma_{\IP,E}\) from (<ref>) solely depends
on the underlying triangulation and
allows for a guaranteed stability of the bilinear from
\(A_h\) from (<ref>) with respect to the
mesh-dependent norm \(\Vert \bullet \Vert_h\)
in the following theorem.

Given any $a > 1$,
define the penalty parameter $\sigma_{\IP,E} > 0$ as in (<ref>).
For every choice of \(a > 1\),
the constant \(\kappa \coloneqq 1 - 1/\sqrt{a} > 0\)
and all discrete functions $v_\IP \in V_h$ satisfy
the stability estimate

\begin{equation}
\label{eq:stability}
A_h(v_\IP, v_\IP) \geq \kappa\: \Vert v_\IP \Vert_h^2.
\end{equation}

Two remarks are in order
before the proof of the Theorem <ref>
concludes this section.

The penalty parameter \(\sigma_{\IP,E}\) in (<ref>)
contains a prefactor $a > 1$.
The striking point of Theorem <ref> is
that every choice of \(a > 1\) leads
to guaranteed stability.
In the case of very large $a \to \infty$, i.e.,
$\sigma_{\IP, E} \to \infty$,
the stability constant $\kappa$ tends to 1.
Hence, the theoretical  stability constant
can be improved by a large penalization.

This is of particular interest,
because optimal convergence rates of the
adaptive discontinuous Galerkin method in [6]
rely on a sufficiently large penalty parameter.

Although the presentation restricts to the quadratic case,
the concept can be generalized for guaranteed stability
in the case of polynomial degree $k \geq 3$.
Even the assumption of uniform polynomial degree can be weakened,
allowing for a different polynomial degree $k_T$
on every triangle $T \in \T$.

The point of departure is the difference,
for every $\kappa > 0$,

\begin{equation}
\label{eq:stability_kappa}
\begin{split}
A_h(v_\IP, v_\IP) - \kappa \Vert v_\IP \Vert_h^2
(1-\kappa) \vvvert v_\IP \vvvert_\PW^2
- 2 \mathcal{J}(v_\IP, v_\IP) \\
&\hphantom{{}={}}+ (1 - \kappa) \sum_{E \in \mathcal{E}}
\frac{\sigma_{\IP,E}}{h_E}
\Vert [\nabla v_\IP \cdot \nu_E]_E \Vert_{L^2(E)}^2.
\end{split}
\end{equation}

The first step of the proof
bounds the jump term $\mathcal{J}$
with the weighted Young inequality,
for every $\varepsilon > 0$,

\[
2 \vert \mathcal{J}(v_\IP, v_\IP) \vert
\quad\leq
\sum \limits_{E \in \mathcal{E}}
\left(
\frac{\varepsilon \sigma_{\IP,E}}{h_E}
\Vert [\nabla v_\IP]_E \Vert_{L^2(E)}^2
+ \frac{h_E}{\varepsilon \sigma_{\IP,E}}
\Vert \langle
\D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\right).
\]

Inserting this estimate into (<ref>)
confirms the lower bound

\begin{equation}
\label{eq:lower_bound_Ah}
\begin{split}
A_h(v_\IP, v_\IP) - \kappa \Vert v_\IP \Vert_h^2\\
(1-\kappa) \vvvert v_\IP \vvvert_\PW^2
+ (1-\kappa) \sum_{E \in \mathcal{E}}
\frac{\sigma_{\IP,E}}{h_E}
\Vert [\nabla v_\IP \cdot \nu_E]_E \Vert_{L^2(E)}^2\\
- \sum_{E \in \mathcal{E}}
\left(
\frac{\varepsilon \sigma_{\IP,E}}{h_E}
\Vert [\nabla v_\IP]_E \Vert_{L^2(E)}^2
+ \frac{h_E}{\varepsilon \sigma_{\IP,E}}
\Vert \langle
\D^2_\PW v_\IP
\rangle_E \Vert_{L^2(E)}^2
\right).
\end{split}
\end{equation}

The second step
estimates the average $\langle \D^2_\PW v_\IP \nu_E \rangle_E$.
Since the piecewise second derivative $\D^2_\PW v_{\IP}$
of any function $v_\IP \in V_h$
is constant on each triangle \(T\),
it holds that

\begin{equation}
\label{eq:Hess_identity}
\Vert \D^2_\PW v_\IP \vert_{T} \Vert_{L^2(E)}^2
\frac{h_E}{\vert T \vert} \Vert \D^2_\PW v_\IP \Vert_{L^2(T)}^2.
\end{equation}

For an interior edge
\(E=\partial T_+ \cap \partial T_- \in \E(\Omega)\)
with the neighboring triangles \(T_+\) and \(T_-\),
the definition of the average from (<ref>)
and \(\vert \nu_E \vert = 1\) verify

\[
\Vert \langle
\D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\leq
\frac{1}{4}
\Vert
\D^2 v_\IP \vert_{T_+} + \D^2 v_\IP \vert_{T_-}
\Vert_{L^2(E)}^2.
\]

The combination with a triangle inequality,
another weighted Young inequality with \(\alpha > 0\),
and (<ref>)
results in

\[
\Vert \langle \D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\leq
\frac{h_E}{4}
\bigg(
\frac{1 + \alpha}{\vert T_+ \vert}
\Vert \D^2 v_\IP \Vert_{L^2(T_+)}^2
\frac{1 + 1/\alpha}{\vert T_- \vert}
\Vert \D^2 v_\IP \Vert_{L^2(T_-)}^2
\bigg).
\]

the optimal value \(\alpha = \vert T_+ \vert / \vert T_- \vert\)
leads to

\[
\Vert \langle \D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\leq
\frac{h_E}{4}
\bigg(
\frac{1}{\vert T_+ \vert}
+ \frac{1}{\vert T_- \vert}
\bigg)
\Vert \D^2_\PW v_\IP \Vert_{L^2(\omega_E)}^2.
\]

The multiplication with
\(h_E / (\varepsilon \sigma_{\IP, E})\)
and the definition of \(\sigma_{\IP, E}\) in (<ref>)

\begin{equation}
\label{eq:average_estimate}
\frac{h_E}{\varepsilon \sigma_{\IP,E}}
\Vert \langle \D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\leq
\frac{1}{3 a \varepsilon}
\Vert \D^2_\PW v_\IP \Vert_{L^2(\omega_E)}^2.
\end{equation}

For a boundary edge \(E \in \E(\partial \Omega)\) with
adjacent triangle \(T_+ \in \T\),
the definition of the average,
the equality \(\vert \nu_E \vert = 1\),
and (<ref>) verify

\[
\Vert \langle
\D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\Vert
\D^2_\PW v_\IP \nu_E
\Vert_{L^2(E)}^2
\leq
\Vert
\D^2_\PW v_\IP
\Vert_{L^2(E)}^2
\frac{h_E}{\vert T_+ \vert}
\Vert \D^2 v_\IP \Vert_{L^2(T_+)}^2.
\]

the multiplication with
\(h_E / (\varepsilon \sigma_{\IP, E})\)
and the definition of the penalty parameter \(\sigma_{\IP, E}\)
result in the analogous estimate (<ref>)
for \(E \in \E(\partial\Omega)\).

The sum of (<ref>) over all edges \(E \in \E\)
and the finite overlap of the edge patches
$(\omega(E): E \in \mathcal{E})$ lead to

\begin{equation}
\label{eq:bound_hessian_normal}
\sum_{E \in \mathcal{E}}
\frac{h_E}{\varepsilon\sigma_{\IP,E}}
\Vert \langle
\D^2_\PW v_\IP \nu_E
\rangle_E \Vert_{L^2(E)}^2
\leq
\frac{1}{a \varepsilon}
\vvvert v_\IP \vvvert_\PW^2.
\end{equation}

The third step investigates the jump
\([\nabla v_\IP]_E\).
The analysis of the jump term employs
the orthogonal split into the tangential and normal component
On each edge $E \in \mathcal{E}$,
the jump $[\nabla v_\IP]_E$ fulfills

\[
[\nabla v_\IP]_E
\big(
[\nabla v_\IP]_E \cdot \tau_E
\big)\tau_E
\big(
[\nabla v_\IP]_E \cdot \nu_E
\big)\nu_E.
\]

For any interior edge
$E = \partial T_+ \cap \partial T_- \in \E(\Omega)$
with the adjacent triangles $T_+$ and $T_-$,
the function $v_\IP$ is differentiable in $T_+$ and $T_-$.
The Hadamard jump condition for the globally continuous function
$v_\IP$ assures that the tangential derivatives
from the adjacent triangles coincide
This means that the term $[\nabla v_\IP]_E \cdot \tau_E$
vanishes and

\begin{equation}
\label{eq:jump_condition}
[\nabla v_\IP]_E
\big(
[\nabla v_\IP]_E \cdot \nu_E
\big)\nu_E.
\end{equation}

For any boundary edge $E \in \mathcal{E}(\partial \Omega)$,
the homogeneous boundary values of the function $v_\IP \in V_h$
guarantee that the tangential derivative vanishes.
This and the jump with zero in the definition (<ref>)

\[
[\nabla v_\IP]_E
\nabla (v_\IP \vert_{T_+}) \vert_E
\big(
\nabla (v_\IP \vert_{T_+}) \vert_E \cdot \nu_E
\big) \nu_E
\big(
[\nabla v_\IP]_E \cdot \nu_E
\big) \nu_E.
\]

Hence, for all \(E \in \E\),

\begin{align}
\label{eq:normal_jump}
\Vert [\nabla v_\IP]_E \Vert_{L^2(E)}^2
\Vert [\nabla v_\IP]_E \cdot \nu_E \Vert_{L^2(E)}^2.
\end{align}

The fourth step concludes the proof with
a combination of the lower bound (<ref>),
the estimate (<ref>), and
the identity (<ref>) for

\begin{equation}
\label{eq:final_lower_bound}
\begin{split}
A_h(v_\IP, v_\IP) - \kappa \Vert v_\IP \Vert_h^2
\geq
&\left(1-\kappa - \frac{1}{a \varepsilon}\right)
\vvvert v_\IP \vvvert_\PW^2 \\
&+ (1-\kappa - \varepsilon)
\sum_{E \in \mathcal{E}}
\frac{\sigma_{\IP,E}}{h_E}
\Vert [\nabla v_\IP \cdot \nu_E]_E \Vert_{L^2(E)}^2.
\end{split}
\end{equation}

Every choice of
\(0 < \kappa \leq \min\{1 - 1/(a \varepsilon), 1-\varepsilon\}\)
leads to a nonnegative lower bound in (<ref>), and proves
the claim (<ref>).
The optimal value \(\varepsilon = 1/ \sqrt{a}\)
allows for \(\kappa = 1 - 1/\sqrt{a}\).

§ DISCRETE FUNCTIONS AND THEIR DERIVATIVES

The presentation of the MATLAB realization
of the C0 interior penalty method
departs with the definition of the quadratic
shape functions.
This section illustrates the evaluation of
their derivatives of first and second order
and provides the foundation for the computation
of the local contributions to the bilinear forms
in \(A_h\) in the
Sections <ref>–<ref>

Two key data structures represent the underlying
triangulation \(\T\).
The rows of the array
\(\in \R^{\vert \V \vert \times 2}\)
contain the two coordinates of each vertex in \(\V\).
The order of the vertices in the array determines
their enumeration.
Three vertices form a triangle and their numbers
are stored in the rows of the array
\(\in \N^{\vert \T \vert \times 3}\)
in counter-clockwise order.

Let \(T \in \T\) denote a triangle
of global number
\(\in \{1, \dots, \vert \T \vert\}\), i.e,
the global indices of its three vertices
\(\V(T) = \{P_1,  P_2, P_3\}\)
belong to the -th row
of .
The enumeration of the midpoints \(M_1, M_2, M_3\)
of the three edges \(\E(T) = \{E_1,  E_2, E_3\}\)
follows the convention that
the edge \(E_j\) is opposite to the node \(P_j\)
for \(j = 1, 2, 3\)
as depicted in Figure <ref>.

Enumeration of nodes

node style/.style=
minimum width=.8cm,
inner sep=0pt

regular polygon,
regular polygon sides=3,
minimum size=4cm
(T) at (0,0) \(T\);

(T.corner 1) node[node style] \(P_1\);
(T.corner 2) node[node style] \(P_2\);
(T.corner 3) node[node style] \(P_3\);
(T.side 1) node[node style] \(M_3\);
(T.side 2) node[node style] \(M_1\);
(T.side 3) node[node style] \(M_2\);

Enumeration of degrees of freedom

node style/.style=
minimum width=.5cm,
inner sep=0pt

regular polygon,
regular polygon sides=3,
minimum size=4cm
(T) at (0,0) \(T\);

(T.corner 1) node[node style] \(1\);
(T.corner 2) node[node style] \(2\);
(T.corner 3) node[node style] \(3\);
(T.side 1) node[node style] \(6\);
(T.side 2) node[node style] \(4\);
(T.side 3) node[node style] \(5\);

Local enumeration of nodes associated
to the piecewise quadratic basis functions in the
the triangle \(T\).

Let \(\lambda_j \in P_1(T)\) denote
the barycentric coordinate
associated to the vertex \(P_j\)
of the triangle \(T\) for \(j = 1, 2, 3\).
\(
\lambda_1(x) P_1 +
\lambda_2(x) P_2 +
\lambda_3(x) P_3
\)
and \(1 = \lambda_1(x) + \lambda_2(x) + \lambda_3(x)\)
for all \(x \in T\),
the gradients \(\nabla\lambda_1,  \nabla \lambda_2, \nabla \lambda_3\)
are determined by the linear system of equations

\begin{equation}
\label{eq:gradient_barycentric_coordinates}
\begin{bmatrix}
1 & 1 & 1\\
P_1 & P_2 & P_3
\end{bmatrix}
\begin{bmatrix}
\nabla \lambda_1^\top \\
\nabla \lambda_2^\top \\
\nabla \lambda_3^\top
\end{bmatrix}
\begin{bmatrix}
0 & 0 \\
1 & 0 \\
0 & 1
\end{bmatrix}
\end{equation}

Let the three-dimensional array
\(\in \R^{2 \times 3 \times \vert \T \vert}\)
contain the matrix
\([P_1\: P_2\: P_3] \in \R^{2 \times 3}\)
for each triangle \(T \in \T\).
The following line of MATLAB computes
the gradients \(\nabla\lambda_1,  \nabla \lambda_2, \nabla \lambda_3\),
according to (<ref>).

GradP1 = [ones(1,3); c4e(:,:,e)]  [zeros(1,2); eye(2)];

The backslash operator in MATLAB calls
the function  and solves the
linear system (<ref>)
up to machine precision.
The three rows of the solution array
\(\in \R^{3 \times 2}\)
contain the two components of the constant gradients.

The shape functions
\(\phi_1, \dots, \phi_6\) of \(P_2(T)\)
can be represented using products
of the barycentric coordinates.
Throughout this section,
let \(\{j, k, \ell\} = \{1, 2, 3\}\)
denote any pairwise distinct indices.
The quadratic shape function
\(\phi_j \coloneqq \lambda_j (2 \lambda_j - 1)\)
is associated to the vertex \(P_j\)
and \(\phi_{3 + j} \coloneqq 4 \lambda_k \lambda_\ell\)
to the midpoint \(M_j\) of \(E_j\).
They satisfy the duality

\[
\phi_j(P_n)
\delta_{jn}
\phi_{3 + j}(M_n)
\quad
\text{and}
\quad
\phi_j(M_n)
\phi_{3 + j}(P_n) \quad \text{for all } j, n = 1, 2, 3.
\]

The enumeration of \(\phi_1, \dots, \phi_6\) is illustrated in
Figure <ref>.
The representation via barycentric coordinates
and the product rule
lead to the following formulas
for the derivatives of first and second order.

For pairwise distinct indices
\(\{j, k, \ell\} = \{1, 2, 3\}\),

\begin{align}
\label{eq:gradient}
\nabla \phi_j
(4 \lambda_j - 1) \nabla \lambda_j,
\nabla \phi_{3 + j}
4 (\lambda_k \nabla \lambda_\ell
+ \lambda_\ell \nabla \lambda_k),\\
\label{eq:Hessian}
\D^2 \phi_j
4 \nabla \lambda_j \otimes \nabla \lambda_j,
\D^2 \phi_{3 + j}
\nabla \lambda_k \otimes \nabla \lambda_\ell
\nabla \lambda_\ell \otimes \nabla \lambda_k
\end{align}

Hence, the evaluation of the gradients
\(\nabla \phi_1, \dots, \nabla \phi_6\)
in the three nodes \(P_k\), \(M_j\), and \(P_\ell\) reads,
for pairwise distinct indices
\(\{j, k, \ell\} = \{1, 2, 3\}\),

\begin{equation}
\label{eq:gradient_evaluation}
\begin{bmatrix}
\nabla \phi_j(P_j) &
\nabla \phi_j(P_k) \\
\nabla \phi_j(M_j) &
\nabla \phi_j(M_k) \\
\nabla \phi_{3 + j}(P_j) &
\nabla \phi_{3 + j}(P_k) \\
\nabla \phi_{3 + j}(M_j) &
\nabla \phi_{3 + j}(M_k)
\end{bmatrix}
\begin{bmatrix}
3 \nabla \lambda_j &
- \nabla \lambda_j \\
- \nabla \lambda_j &
\nabla \lambda_j \\
0 &
4 \nabla \lambda_\ell \\
2 (\nabla \lambda_k + \nabla \lambda_\ell) &
2 \nabla \lambda_k
\end{bmatrix}.
\end{equation}

The following lines of MATLAB realize the evaluation
of the gradients from (<ref>).
Gradient = cat(3, [3 -1 -1 -1 1 1] .* GradP1(1,:)', ...
[-1 3 -1 1 -1 1] .* GradP1(2,:)', ...
[-1 -1 3 1 1 -1] .* GradP1(3,:)', ...
[0 0 4 2 2 0] .* GradP1(2,:)' ...
+ [0 4 0 2 0 2] .* GradP1(3,:)', ...
[0 0 4 2 2 0] .* GradP1(1,:)' ...
+ [4 0 0 0 2 2] .* GradP1(3,:)', ...
[0 4 0 2 0 2] .* GradP1(1,:)' ...
+ [4 0 0 0 2 2] .* GradP1(2,:)');

Linear combinations of the gradients
of the barycentric coordinates from
provide the formulas in (<ref>).
They are rearranged and stored in the
three-dimensional array
\(\in \R^{2 \times 6 \times 6}\).

The following lines of MATLAB evaluate the piecewise constant
Hessian matrices from (<ref>).
H = permute(GradP1, [2 3 1 4]) .* permute(GradP1, [3 2 4 1]);
Hessian = 4 * cat(3, H(:,:,1,1), H(:,:,2,2), H(:,:,3,3), ...
H(:,:,2,3) + H(:,:,3,2), ...
H(:,:,1,3) + H(:,:,3,1), ...
H(:,:,1,2) + H(:,:,2,1));

All pairwise tensor products of the gradients
of the barycentric coordinates
from Lemma <ref>
result in nine 2-by-2 matrices.
They are stored in the four-dimensional array
\(\in \R^{2 \times 2 \times 3 \times 3}\).
The linear combinations of these matrices
provide the Hessian matrices
\(\D^2 \phi_1, \dots, \D^2 \phi_6\)
in the array
\(\in \R^{2 \times 2 \times 6}\).

The Hessian matrices
\(\in \R^{2 \times 2 \times 6}\)
from Lemma <ref>
and the gradient evaluations
\(\in \R^{2 \times 6 \times 6}\)
are computed elementwise and
stored in the global variables

\(\in \R^{6 \times 2 \times 2 \times \vert \T \vert}\)

\(\in \R^{6 \times 2 \times 6 \times \vert \T \vert}\)
Hess4e(:,:,:,e) = permute(Hessian, [3 1 2]);
Grad4e(:,:,:,e) = permute(Gradient, [3 1 2]);

The evaluations of the derivatives are supposed to be used
in a parallelized computation over all edges
for the local contributions in the
Sections <ref>–<ref> below.
In order to comply with the requirements of MATLAB's parallel toolbox,
the data structures are rearranged in an edge-oriented way

\[
\mathtt{Hess4s}
\in \R^{6 \times 2 \times 2 \times 2 \times \vert \E \vert}
\quad\text{and}\quad
\mathtt{Grad4s}
\in \R^{6 \times 2 \times 6 \times 2 \times \vert \E \vert}.
\]

For each edge, they contain the corresponding information
on the adjacent triangles \(T_+, T_- \in \T\)
specified by the index in the fourth component
(index 1 for \(T_+\) and 2
for \(T_-\)).
The rearrangement employs a data structure
whose rows contain the global indices
of \(T_+\) and \(T_-\) for interior edges
or the global index of \(T_+\) and \(0\) for boundary edges.
For the realization, see the lines 70–75 of the program
in Appendix <ref>.

§ STIFFNESS MATRIX AND RIGHT-HAND SIDE

The contribution to the bilinear form \(a_\PW\)
from the triangle \(T \in \T\)
leads to the local stiffness matrix
\(A(T) \in \R^{6 \times 6}\) with,
for \(\alpha, \beta = 1,\dots, 6\),

\begin{equation}
\label{eq:local_stiffness}
\coloneqq
\int_T \D^2 \phi_\alpha : \D^2 \phi_\beta \;\mathrm{d}x
\vert T \vert\, (\D^2 \phi_\alpha : \D^2 \phi_\beta).
\end{equation}

Let the array
\(\in \R^{\vert \T \vert}\)
contain the area \(\vert T \vert\)
of each triangle \(T \in \T\).
The following lines of MATLAB
compute (<ref>).
Integrand = sum(Hessian .* permute(Hessian, [1 2 4 3]), [1 2]);
A4e(:,:,e) = area4e(e) * squeeze(Integrand);

The pairwise Frobenius scalar products
of the Hessian matrices
\(\D^2\phi_\alpha\), \(\alpha = 1,\dots,6\)
from the array
provide the array
\(\in \R^{6 \times 6}\).
The multiplication with the area
completes the computation of (<ref>).

The local contribution \(b(T) \in \R^6\)
from the triangle \(T \in \T\) to the right-hand side
consists of approximations to the integral
\(\int_T f \phi_\alpha \,\mathrm{d} x\) for
\(\alpha = 1, \dots, 6\).
Its computation employs
the edge-oriented midpoint quadrature rule

\begin{equation}
\label{eq:local_RHS}
\coloneqq
\frac{\vert T \vert}{3}
\sum_{\ell = 1}^3
f(M_\ell) \phi_\alpha(M_\ell)
\begin{cases}
0 & \text{for } \alpha = 1, 2, 3,\\
\vert T \vert\, f(M_{\alpha - 3})
& \text{for } \alpha = 4, 5, 6.
\end{cases}
\end{equation}

Let the rows of the array

\(\in \N^{\vert \E \vert \times 2}\)
specify the global numbers
of the two endpoints of each edge.
Its order determines an enumeration of the edges
of the triangulation and the numbers
of the three edges of each triangle
form the rows of an array

\(\in \N^{\vert \T \vert \times 3}\).
Given a function handle
describing the right-hand side \(f: \overline\Omega \to \R\),
the evaluations of \(f\) required for (<ref>)
f4mid = f((c4n(n4s(:,1),:) + c4n(n4s(:,2),:)) / 2);
valF4e = f4mid(s4e);

On the triangle \(T\) with index ,
the following lines of MATLAB compute (<ref>).
b4e(:,e) = [zeros(1, 3), area4e(e) / 3 * valF4e(e,:)];

The right-hand side \(f\) is evaluated
in the midpoints of each edge first.
The array  allows
to select the three point evaluations
of \(f\) on each triangle into the array
\(\in \R^{\vert \T \vert \times 3}\).
On the triangle \(T\),
the multiplication with the quadrature weigths
completes the computation of (<ref>).

§ JUMP

This section presents the MATLAB realization
of averages and jumps of the derivatives
from Section <ref>
as well as the local contributions
to the bilinear form \(\mathcal{J}\).
Let \(E \in \E\) denote the edge
of global number
\(\in \{1, \dots, \vert \E \vert \}\), i.e,
the global indices of its two vertices
belong to the -th row
of the array .

If \(E \in \E(\Omega)\) is an interior edge,
two adjacent triangles \(T_+, T_- \in \T\)
share the edge  \(\partial T_+ \cap \partial T_- = E\).
The local enumeration of the interpolation nodes
\(P^\pm_1, P^\pm_2, P^\pm_3\), \(M^\pm_1, M^\pm_2, M^\pm_3\)
in \(T_\pm\) from Figure <ref>
determine two indices \(q, r \in \{1, 2, 3\}\)
fixed for each edge \(E\) such that
\(\Mid(E) = M^+_q = M^-_r\).
The remaining nodes are numbered accordingly
in counter-clockwise order
as displayed in Figure <ref>
in a cyclic notation with index $0$ resp. $1$ identified
with $3$ resp. $4$.

Enumeration in the adjacent triangles

line cap=round,
line join=round,
node style/.style=
minimum width=.8cm,
inner sep=0pt
regular polygon,
regular polygon sides=3,
minimum size=4cm,
] (Tp) at (0, 0) ;
(Tp) node \(T_+\);

(Tp.corner 1) node[node style] \(P^+_{q-1}\);
(Tp.corner 2) node[node style] \(P^+_q\);
(Tp.corner 3) node[node style] \(P^+_{q+1}\);
(Tp.side 1) node[node style] \(M^+_{q+1}\);
(Tp.side 2) node[node style] \(M^+_{q-1}\);
(Tp.side 3) node[node style] \(M^+_q\);

regular polygon,
regular polygon sides=3,
minimum size=4cm,
] (Tm) at (2.8, 1.6) ;
(Tm) node \(T_-\);

(Tm.corner 1) node[node style] \(P^-_{r-1}\);
(Tm.corner 2) node[node style] \(P^-_r\);
(Tm.corner 3) node[node style] \(P^-_{r+1}\);
(Tm.side 1) node[node style] \(M^-_{r+1}\);
(Tm.side 2) node[node style] \(M^-_{r-1}\);
(Tm.side 3) node[node style] \(M^-_r\);

[very thick] ($(Tp.corner 3)!0.5!(Tm.corner 1)$)
node[left] \(E\)
– ($(Tp.corner 1)!0.5!(Tm.corner 3)$);

Enumeration in the edge patch

line cap=round,
line join=round,
node style/.style=
minimum width=.7cm,
inner sep=0pt

regular polygon,
regular polygon sides=3,
minimum size=4cm
] (T) at (0, 0) ;
(T) node \(T_+\);
($(T)-(T.corner 2)$) node \(T_-\);

($2*(T.side 3)-(T.corner 2)$) node[coordinate] (N) \(8\);

(T.corner 3) – node[node style] \(Q_7\) (N)
node[node style] \(Q_8\)
– node[node style] \(Q_9\)
(T.corner 1);

(T.corner 1) node[node style] \(Q_3\);
(T.corner 2) node[node style] \(Q_5\);
(T.corner 3) node[node style] \(Q_1\);
(T.side 1) node[node style] \(Q_4\);
(T.side 2) node[node style] \(Q_6\);
(T.side 3) node[node style] \(Q_2\);

Local enumeration of nodes
associated to the piecewise quadratic basis functions
in the adjacent triangles \(T_+\) and \(T_-\)
and in the whole patch \(\omega_E\) of the edge \(E\).
The local enumeration in each triangle
from Figure <ref>
induces the choice two indices \(q, r \in \{1,  2, 3\}\)
such that \(\Mid(E) = M^+_q = M^-_r\).

This information is stored in an array
\(\in \N^{6 \times 2}\)
containing the corresponding permutations of
\(\{1, \dots, 6\}\) as a column vector
or \(T_+\) and \(T_-\).

Abbreviate the nodes belonging to
the edge patch \(\omega_E\) for \(E \in \E(\Omega)\) by

\begin{equation}
\label{eq:Simpson_nodes}
\begin{split}
\big[
Q_1, Q_2, Q_3
\big]
\big[
\big]
\big[
\big],\\
\big[
Q_4, \dots, Q_9
\big]
\big[
\big].
\end{split}
\end{equation}

The nine global basis functions \(\phi_1, \dots, \phi_9 \in S^2(\T)\)
of the edge patch \(\omega_E\)
are associated to the nodes \(Q_1, \dots, Q_9\)
as in Figure <ref>.
This means that
\(\phi_1, \phi_2, \phi_3\) belong to the edge \(E\),
\(\phi_4, \phi_5, \phi_6\) exclusively to \(T_+\),
and \(\phi_7, \phi_8, \phi_9\) exclusively to \(T_-\).
For \(\alpha = 1, \dots, 9\),
abbreviate the restrictions
\(\phi^\pm_\alpha \coloneqq \phi_\alpha\vert_{T_\pm} \in P_2(T_\pm)\) of $\phi_\alpha$ to $T_\pm$.

In the case of a boundary edge \(E \in \E(\partial\Omega)\),
the nodes \(Q_7, Q_8, Q_9\) corresponding to \(T_-\) are omitted
and the basis functions \(\phi_\alpha \coloneqq 0\)
vanish for \(\alpha = 7, 8, 9\).

The evaluations of the normal jumps
of the gradients \(\nabla \phi_1, \dots, \nabla \phi_9\)
in the nodes \(Q_1, Q_2, Q_3\) on the edge \(E\)
for \(\alpha = 1, \dots, 9\)
and \(\ell = 1, 2, 3\),

\begin{equation}
\label{eq:evaluation_normal_jump}
[\nabla \phi_\alpha (Q_\ell) \cdot \nu_E]_E
\begin{cases}
\big(
\nabla \phi^+_\alpha (Q_\ell)
- \nabla \phi^-_\alpha (Q_\ell)
\big) \cdot \nu_E
&\quad\text{if } E \in \E(\Omega),\\
\nabla \phi^+_\alpha (Q_\ell) \cdot \nu_E
&\quad\text{if } E \in \E(\partial\Omega),
\end{cases}
\end{equation}

In particular,
the vanishing tangential component of the gradient jumps
from (<ref>)

\begin{equation}
\label{eq:evaluation_jump}
[\nabla \phi_\beta (Q_2)]_E
[\nabla \phi_\beta (Q_2) \cdot \nu_E]_E \nu_E.
\end{equation}

Let the rows of the array
\(\in \R^{\vert \E \vert \times 2}\)
contain the two components of the unit normal
vector \(\nu_E \in \R^2\) for each edge \(E\).
Assume that its sign is chosen such that
\(\nu_E = \nu_{T_+}\) for the outward unit normal vector
\(\nu_{T_+}\) of the triangle \(T_+\).
Using the five-dimensional array
from Remark <ref>,
the following lines of MATLAB compute
the jumps (<ref>)–
GradNormal = squeeze(sum(Grad4s(:,:,:,:,s) .* normal4s(s,:), 2));
JumpNormal = [GradNormal(Ind(1:3,1),Ind(1:3,1),1) ...
- GradNormal(Ind(3:-1:1,2),Ind(3:-1:1,2),2); ...
GradNormal(Ind(4:6,1),Ind(1:3,1),1); ...
- GradNormal(Ind(4:6,2),Ind(3:-1:1,2),2)];
JumpGrad = JumpNormal(:,2) * normal4s(s,:);

The array
\(\in \R^{6 \times 6 \times 2}\)
contains the evaluations of the normal derivatives
\(\nabla \phi_\alpha^\pm \cdot \nu_E\) for \(\alpha = 1,\dots, 6\)
at the nodes \(P_1^\pm, P_2^\pm, P_3^\pm\),
\(M_1^\pm, M_2^\pm, M_3^\pm\) on \(T_\pm\).

Their linear combination lead to the jumps
along the edge in the array
\(\in \R^{9 \times 3}\)
for the basis functions \(\phi_\alpha\)
associated to the nodes \(Q_\alpha\)
for \(\alpha = 1,\dots,9\).
The first three basis functions \(\phi_1, \phi_2, \phi_3\)
belong to the edge \(E\) and the jump
consists of the difference of contributions
from \(T_+\) and \(T_-\).
The array  selects the corresponding
basis functions on the adjacent triangles
as well as the corresponding nodes \(Q_1, Q_2, Q_3\)
belonging to the edge \(E\).
Note that  vanishes on \(T_-\)
for boundary edges.
The basis functions \(\phi_4, \phi_5, \phi_6\)
vanish on \(T_-\) and
the jump of their normal derivatives from \(T_+\)
is computed with 0.
Vice versa, the basis functions \(\phi_7, \phi_8, \phi_9\)
lead to the jump of 0 with the normal derivatives from \(T_-\).

The selection of the evaluation in \(Q_2\)
and the multiplication with the unit normal vector
\(\nu_E\) in
realize the computation of (<ref>)
for  \(\in \R^{9 \times 2}\).

The averages of the piecewise Hessian matrices
from (<ref>) read,
for \(\alpha = 1, \dots, 9\),

\begin{equation}
\label{eq:Hessian_average}
\langle \D^2_\PW \phi_\alpha \nu_E \rangle_E
\begin{cases}
\frac12
(\D^2 \phi^+_\alpha + \D^2 \phi^-_\alpha) \nu_E
&\quad\text{if } E \in \E(\Omega),\\
\D^2 \phi^+_\alpha \, \nu_E
&\quad\text{if } E \in \E(\partial\Omega).
\end{cases}
\end{equation}

Let the boolean array

\(\in \{0, 1\}^{\vert \E \vert}\)
determine whether an edge belongs to the
interior of the domain.
The following lines of MATLAB
compute the average (<ref>).

HessNormal = squeeze(sum(Hess4s(:,:,:,:,s) .* normal4s(s,:), 2));
MeanHessNormal = ...
[HessNormal(Ind(1:3,1),:,1) + HessNormal(Ind(3:-1:1,2),:,2);...
HessNormal(Ind(4:6,1),:,1); HessNormal(Ind(4:6,2),:,2)];
if isInteriorSide(s); MeanHessNormal = MeanHessNormal / 2; end

The array
\(\in \R^{6 \times 2 \times 2}\)
contains the product
\(
(\nu_E^\top \D^2 \phi^\pm_\alpha)^\top
\D^2 \phi^\pm_\alpha \, \nu_E
\)
of each Hessian matrix
\(\D^2 \phi^\pm_\alpha\) for \(\alpha = 1, \dots, 6\)
on \(T_+\) and \(T_-\)
with the unit normal vector \(\nu_E\).

The linear combination for the averages
\(\langle \D^2_\PW \phi_\alpha \nu_E \rangle_E\)
along \(E\) for \(\alpha = 1, \dots, 9\)
is computed analogously to
in the proof of Lemma <ref>
with constant Hessian matrices replacing
the gradient evaluations.
The two components of the vector-valued averages are
stored in the array
\(\in \R^{9 \times 2}\).

Finally, the prefactor \(1/2\) is applied for
interior edges only according to (<ref>).

The contribution to the bilinear form \(\mathcal{J}\)
from the edge \(E \in \E\)
consists of the local matrix \(J(E) \in \R^{9 \times 9}\) with,
for \(\alpha, \beta = 1,\dots, 9\),

\begin{equation}
\label{eq:local_jump_term}
\coloneqq
\int_E
\langle \D^2_\PW \phi_\alpha \nu_E \rangle_E
\cdot
[\nabla \phi_\beta]_E
\;\mathrm{d}s
\vert E \vert\,
\langle \D^2_\PW \phi_\alpha \nu_E \rangle_E
\cdot
[\nabla \phi_\beta (Q_2)]_E.
\end{equation}

If \(E \in \E(\partial\Omega)\),
\(\phi_7, \phi_8, \phi_9 \equiv 0\)
implies that
\(J_{\alpha\beta}(E) = 0\) vanishes for all \(\alpha, \beta = 1,\dots,9\)
with \(7 \leq \alpha\) or \(7 \leq \beta\).

Let the array  \(\in \R^{\vert \E \vert}\)
contain the length of each edge.
The following line of MATLAB
computes (<ref>).
J4s(:,:,s) = length4s(s) * (MeanHessNormal * JumpGrad');

The pairwise scalar products for
\(\alpha, \beta = 1,\dots, 9\)
of the evaluations from the
Lemmas <ref>–<ref>
followed by the multiplication
with the length of the edge \(E\)
realize to the computation of
the integral in (<ref>).

§ PENALIZATION

Recall the enumeration of nodes and basis functions
in the edge patch $\omega_E$ from Section <ref>.

Given a parameter
\( \coloneqq a > 1\),
the following lines of MATLAB compute
the automatic penalization parameter \(\sigma_{\IP,E}\)
from (<ref>).
sigma4s = 3 * a * (length4s.^2) ./ area4e(e4s(:,1));
sigma4s(isInteriorSide) = sigma4s(isInteriorSide) / 4 ...
+ 3/4 * a * length4s(isInteriorSide).^2 ...
./ area4e(e4s(isInteriorSide,2));

The parameter \(\sigma_{\IP,E}\)
consists of the sum of one contribution
from each adjacent triangle \(T_+\) and \(T_-\).
The value for all edges \(E \in \E\) is initialized
by \(3a\,  h_E^2 / \vert T_+ \vert\).
In case of an interior edge,
the contribution from \(T_+\)
is multiplied with the prefactor \(1/4\)
and the term \(3a\, h_E^2 / (4\vert T_- \vert)\)
for \(T_-\) is added.

The contribution to the penalty term \(c_\IP\)
from the edge \(E \in \E\)
leads to the local matrix \(C(E) \in \R^{9 \times 9}\) with,
for \(\alpha, \beta = 1, \dots, 9\),

\[
\coloneqq
\frac{\sigma_{\IP,E}}{h_E}
\int_E
[\nabla \phi_\alpha \cdot \nu_E]_E\,
[\nabla \phi_\beta \cdot \nu_E]_E
\;\mathrm{d}s.
\]

The three-point one-dimensional Simpson quadrature rule
with quadrature nodes \(Q_1,  Q_2, Q_3\)
from (<ref>)
and weights \(w_1 = w_3 = 1/6\) and \(w_2 = 2/3\)
integrates the product of the jump terms exactly and

\begin{equation}
\label{eq:local_penalty_term}
\sigma_{\IP,E}
\sum_{j = 1}^3
[\nabla \phi_\alpha(Q_j) \cdot \nu_E]_E \,
[\nabla \phi_\beta(Q_j) \cdot \nu_E]_E.
\end{equation}

If \(E \in \E(\partial\Omega)\),
\(\phi_7, \phi_8, \phi_9 \equiv 0\)
implies that
\(C_{\alpha\beta}(E) = 0\) vanishes for all \(\alpha, \beta = 1,\dots,9\)
with \(7 \leq \alpha\) or \(7 \leq \beta\).

The following lines of MATLAB
compute (<ref>).
Integrand = JumpNormal .* permute(JumpNormal, [3 2 1]);
C4s(:,:,s) = sigma4s(s)/6 * squeeze(sum([1 4 1] .* Integrand, 2));

The pairwise products of the evaluations
of the normal jumps
\([\nabla \phi_\alpha \cdot \nu_E]_E\),
\(\alpha = 1, \dots, 9\),
from Lemma <ref>
form the array
\(\in \R^{9 \times 3 \times 9}\).
The multiplication with the penalty parameter
\(\sigma_{\IP,E}\) and the quadrature weights
followed by the sum over the three terms
of the quadrature rule
complete the computation of (<ref>).

§ ASSEMBLING

This section describes the assembling of the global
system matrix using the local contributions
from the Sections <ref>–

To this end,
let \(\V = \{z_1, \dots, z_K\}\) denote
a global enumeration of
the \(K \coloneqq \vert \V \vert \in \N\) vertices
and \(\E = \{E_1, \dots, E_L\}\) of
the \(L \coloneqq \vert \E \vert \in \N\) edges
of the triangulation \(\T\).
The edge midpoints
\(\Mid(\E) = \{m_\ell \coloneqq \Mid(E_\ell) \;:\; \ell = 1,\dots, L\}\)
are numbered accordingly.
This induces a global enumeration of the
\(\mathtt{N} \coloneqq N \coloneqq K + L\) basis functions
\(\varphi_1, \dots, \varphi_N \in S^2(\T)\) via,
for all \(j,k = 1, \dots, K\) and \(\ell,n = 1, \dots, L\),

\[
\varphi_j(z_k) = \delta_{jk},
\quad
\varphi_{K+\ell}(m_n) = \delta_{\ell n},
\quad\text{and}\quad
\varphi_j(m_n) = 0 = \varphi_{K+\ell}(z_k).
\]

The bilinear form \(a_\PW\)
leads to the global stiffness matrix \(A \in \R^{N \times N}\)

\[
\coloneqq
a_\PW(\varphi_j, \varphi_k)
\sum_{T \in \T}
\int_T
\D^2 \varphi_j : \D^2 \varphi_k
\;\mathrm{d}x \quad \text{for all } j,k = 1, \dots, N.
\]

The integration of the right-hand side
\(\int_\Omega f \, \varphi_j \dx\)
from (<ref>)
employs the edge-oriented midpoint quadrature rule.
The resulting right-hand side vector \(b \in \R^N\) reads,
for all \(j = 1, \dots, K\) and \(\ell = 1,\dots, L\),

\[
\coloneqq
\quad\text{and}\quad
b_{K + \ell}
\coloneqq
\sum_{T \in \T, E_\ell \subset T}
\frac{\vert T \vert}{3}
\]

The assembling of \(A\) and \(b\)
with the local contributions \(A(T)\) and \(b(T)\)
from Section <ref>
requires the translation of
local indices of shape functions on the triangle \(T \in \T\)
into global indices of basis functions.
The index set

\begin{align*}
\big\{
(\alpha, j) \in \{1,2,3\} \times \{1, \dots, K\}
\;:\;
P_\alpha = z_j
\big\}\\
\cup
\big\{
(3 + \alpha, K + \ell)
\;:\;
(\alpha, \ell) \in \{1,2,3\} \times \{1, \dots, L\}
\text{ with }
M_\alpha = m_\ell
\big\}.
\end{align*}

is stored in the variable

in that the index \(j\) has position \(\alpha\)
in the column associated to \(T\)
and analogously for \(K+\ell\) and \(3 + \alpha\).
The set \(I(T)\) leads to the assembling of
\(A \in \R^{N \times N}\)
and \(b \in \R^N\)

\begin{align}
\label{eq:assembling_A}
\sum_{T \in \T}
\sum_{(\alpha, j) \in I(T)}
\sum_{(\beta, k) \in I(T)}
A_{\alpha\beta}(T) \, e_j \otimes e_k,
\\
\label{eq:assembling_b}
\sum_{T \in \T}
\sum_{(\alpha, j) \in I(T)}
b_\alpha(T) \, e_j.
\end{align}

The following lines of MATLAB compute \(A\) and \(b\)
from (<ref>)–(<ref>).
RowIndices4e = permute(repmat(dof4e, [1 1 size(dof4e, 2)]), [1 3 2]);
ColIndices4e = permute(RowIndices4e, [2 1 3]);
A = sparse(RowIndices4e(:), ColIndices4e(:), A4e(:), N, N);
b = accumarray(dof4e(:), b4e(:), [N 1]);

The indices of \(I(T)\) in
are rearranged to match the requirements
of MATLAB's  command.
This results in the two arrays
\(\in \N^{6 \times 6 \times \vert \T \vert}\)
of the same size as .
For each entry in ,
the corresponding entries in
provide the row and column indices of the global
position in .
The  function realizes the
same procedure for dense one-dimensional

The bilinear forms \(\mathcal{J}\) and \(c_\IP\)
lead to the global matrices
\(J \in \R^{N \times N}\)
and \(C \in \R^{N \times N}\)
with, for all \(j,k = 1, \dots, N\),

\begin{align*}
\mathcal{J}(\varphi_j, \varphi_k)
\sum_{E \in \E}
\int_E
\langle \D^2_\PW \varphi_j \nu_E \rangle_E
\cdot
[\nabla \varphi_k]_E
\;\mathrm{d}s,\\
c_\IP(\varphi_j, \varphi_k)
\sum_{E \in \E}
\frac{\sigma_{\IP,E}}{h_E}
\int_E
[\nabla \varphi_j \cdot \nu_E]_E\,
[\nabla \varphi_k \cdot \nu_E]_E
\;\mathrm{d}s.
\end{align*}

Recall the enumeration of the nodes \(Q_1, \dots, Q_9\)
from (<ref>)
as displayed in Figure <ref>.

The translation of
local indices of shape functions on the edge patch \(\omega_E\)
into global indices of basis functions
employs the index sets
\(I(E) \coloneqq I^+(E) \cup I^-(E)\)
defined by

\begin{align*}
\big\{
(\alpha, j)
\in \{1,3,5\} \times \{1, \dots, K\}
\;:\;
Q_\alpha = z_j
\big\}
\\
\cup
\big\{
(\alpha, M + \ell)
\;:\;
\alpha \in \{2,4,6\}, \ell \in \{1, \dots, L\},
\:
Q_\alpha = m_\ell
\big\}
\end{align*}

and \(I^-(E)\) according to the location of \(E\).
For an interior edge \(E \in \E(\Omega)\),
the indices from \(T_-\) read

\begin{align*}
\big\{
(8, j) \;:\; j \in \{1, \dots, K\},\: Q_8 = z_j
\big\}
\\
\cup
\big\{
(\alpha, K + \ell)
\;:\;
\alpha \in \{7, 9\},
\ell \in \{1, \dots, K\},\: Q_\alpha = m_\ell
\big\}.
\end{align*}

For a boundary edge \(E \in \E(\partial\Omega)\),
recall that the local contributions
\(J_{\alpha\beta}(E) = C_{\alpha\beta}(E) = 0\)
from the Sections <ref> and <ref>
vanish for all \(\alpha, \beta = 1,\dots,9\)
with \(7 \leq \alpha\) or \(7 \leq \beta\).
The vanishing entries do not influence the assembling
of the global matrix and
the corresponding index pairs can be chosen arbitrarily,
\(
\coloneqq
\{(7,1), (8,1), (9,1)\}
\).

The set \(I(E)\) is stored in a variable
\(\in \N^{\vert \E \vert \times 9}\)
similarly to  for \(I(T)\).
For the detailed realization, the reader is referred
to the lines 23–41 of the code
in the Appendix <ref>.

The assembling of \(J \in \R^{N \times N}\)
with the local contributions \(J(E)\)
from Section <ref> reads

\begin{equation}
\label{eq:edge_assembling}
\sum_{E \in \E}
\sum_{(\alpha, j) \in I(E)}
\sum_{(\beta, k) \in I(E)}
J_{\alpha\beta}(E) \, e_j \otimes e_k,
\end{equation}

and analogously for
\(C \in \R^{N \times N}\)
with \(C_{\alpha\beta}(E)\)
from Section <ref>
replacing \(J_{\alpha\beta}(E)\).

The MATLAB realization of (<ref>)
is verbatim to Lemma <ref>
with  replacing .

The sum of the three contributions results in the
system matrix

\[
\coloneqq
A - J - J^\top + C
\in
\R^{N \times N}.
\]

Since the discrete solution \(u_\IP \in S^2_0(\T)\)
to (<ref>)
vanishes along the boundary,
the indices of the degrees of freedom in the set

\[
\coloneqq
\big\{
j \in \{1, \dots, M\}
\;:\;
z_j \in \V(\Omega)
\big\}
\cup
\big\{
M + \ell
\;:\;
\ell \in \{1, \dots, L\},\:
E_\ell \in \E(\Omega)
\big\}
\]

restrict to all interpolation nodes
in the interior of the domain.

The discrete formulation (<ref>)
seeks the coefficient vector
\(x \in \R^{M + L}\) of

\[
\sum_{j = 1}^{M + L}
x_j \varphi_j
\in
\]

with \(x_j = 0\) for all boundary indices
\(j \in \{1, \dots, M + L\} \setminus I(\Omega)\)
the linear system of equations

\[
(B_{jk})_{j,k \in I(\Omega)} (x_k)_{k \in I(\Omega)}
(b_j)_{j \in I(\Omega)}.
\]

The solution of this system concludes the realization
of the C0 interior penalty
method from (<ref>).

§ NUMERICAL EXPERIMENTS

This section investigates the implementation from Appendix <ref>
for four computational benchmark examples with domains
in Figure <ref>.

[L-shaped domain]

eta/.style =
mark = *,
every mark/.append style=
scale = 1.1,
fill  = colEta!60!white

uniform/.style =
every mark/.append style = fill = black!20!white
adaptive/.style = solid

axis equal,
width            = 0.3,
ymajorgrids      = true,
font             = ,
legend style     =
legend columns = 1,
legend pos     = outer north east,
font           =
[HUblue!50!white] coordinates (1,-1)(0,0)(0,-1);

[HUblue!50!white] coordinates (-1,-1)(0,0)(-1,0);

[HUblue!50!white] coordinates (-1,1)(0,0);

[black,thick] coordinates (-1,-1)(-1,1)(0,1)(0,0)(1,0)(1,-1)(-1,-1);

[1/8 cusp domain]

eta/.style =
mark = *,
every mark/.append style=
scale = 1.1,
fill  = colEta!60!white

uniform/.style =
every mark/.append style = fill = black!20!white
adaptive/.style = solid

axis equal,
width            = 0.3,
ymajorgrids      = true,
font             = ,
legend style     =
legend columns = 1,
legend pos     = outer north east,
font           =
[HUblue!50!white] coordinates (1,-1)(0,0)(0,-1);

[HUblue!50!white] coordinates (-1,-1)(0,0)(-1,0);

[HUblue!50!white] coordinates (-1,1)(0,0)(0,1);

[black,thick] coordinates (-1,-1)(-1,1)(1,1)(0,0)(1,0)(1,-1)(-1,-1);

[Dumbbell-slit domain]

eta/.style =
mark = *,
every mark/.append style=
scale = 1.1,
fill  = colEta!60!white

uniform/.style =
every mark/.append style = fill = black!20!white
adaptive/.style = solid

width            = 0.55,
ymajorgrids      = true,
font             = ,
axis equal image,
legend style     =
legend columns = 1,
legend pos     = outer north east,
font           =
[HUblue!50!white] coordinates (-1,-1)(0,0)(0,-1)(1,0)(0,0)(1,1);
[HUblue!50!white] coordinates (-1,1)(0,0)(0,1);
[HUblue!50!white] coordinates (0.5,-1)(0.5,-0.5)(1,-0.5)(0.5,-1);
[HUblue!50!white] coordinates (0.75,-1)(0.75,-0.75)(1,-0.75);

[HUblue!50!white] coordinates (0,0)(1,-1)(1,-0.75)(1.25,-1)(1.25,-0.75)(1.5,-1)(1.5,-0.75)(1.75,-1)(1.75,-0.75)(2,-1)(2,-0.75);
[HUblue!50!white] coordinates (2,-1)(2.25,-0.75)(2.25,-1)(2.5,-0.75)(2.5,-1)(2.75,-0.75)(2.75,-1)(3,-0.75)(3,-1);

[HUblue!50!white] coordinates (5,-1)(4,0)(4,-1)(3,0)(4,0)(3,1);
[HUblue!50!white] coordinates (3,-1)(5,1);
[HUblue!50!white] coordinates (3,-0.75)(3.25,-0.75)(3.25,-1);
[HUblue!50!white] coordinates (3.5,-1)(3.5,-0.5)(3,-0.5)(3.5,-1);
[HUblue!50!white] coordinates (5,0)(4,0)(4,1);

[black,thick] coordinates (-1,-1)(5,-1)(5,1)(3,1)(3,-0.75)(1,-0.75)(1,1)(-1,1)(-1,-1);
[black,thick] coordinates (-1,0)(0,0);

[Four-slit domain]

eta/.style =
mark = *,
every mark/.append style=
scale = 1.1,
fill  = colEta!60!white

uniform/.style =
every mark/.append style = fill = black!20!white
adaptive/.style = solid

axis equal,
width            = 0.3,
ymajorgrids      = true,
font             = ,
legend style     =
legend columns = 1,
legend pos     = outer north east,
font           =
[HUblue!50!white] coordinates (-1,-0.5)(1,-0.5);
[HUblue!50!white] coordinates (-1,0.5)(1,0.5);

[HUblue!50!white] coordinates (-0.5,-1)(-0.5,1);
[HUblue!50!white] coordinates (0.5,-1)(0.5,1);

[HUblue!50!white] coordinates (-0.5,0)(0.5,0);
[HUblue!50!white] coordinates (0,-0.5)(0,0.5);

[HUblue!50!white] coordinates (-0.5,-1)(1,0.5)(0.5,1)(-1,-0.5)(-0.5,-1);
[HUblue!50!white] coordinates (0.5,-1)(1,-0.5)(-0.5,1)(-1,0.5)(0.5,-1);

[black,thick] coordinates (-1,-1)(1,-1)(1,1)(-1,1)(-1,-1);

[black,thick] coordinates (0,-1)(0,-0.5);
[black,thick] coordinates (0,0.5)(0,1);
[black,thick] coordinates (-1,0)(-0.5,0);
[black,thick] coordinates (0.5,0)(1,0);

Visualization of the initial triangulations $\mathcal{T}_0$
for the four benchmark problems.

The experiments focus on a comparison
of uniform and adaptive mesh-refinement
driven by the a posteriori error estimator \(\eta\)
from Appendix <ref> below.
This employs the Dörfler marking strategy [17] and
newest-vertex bisection according to [23]
in the adaptive Algorithm <ref>.
The marking strategy with minimal cardinality of
the set \(\mathcal{M}_\ell\) of marked triangles
can be realized in linear computational complexity

regular triangulation $\T_0$ and bulk parameter $0 < \theta \leq 1$.
$\ell = 0,1,2,\dots$
Solve (<ref>) with respect to triangulation
$\T_\ell$ for solution $u_\ell \in S^2_0(\T_\ell)$.
Compute refinement indicator $\eta^2(\T_\ell, T)$
from (<ref>)
for all \(T \in \T_\ell\).
Mark a minimal subset
$\mathcal{M}_\ell \subseteq \T_\ell$
by the Dörfler criterion

\[
\theta\: \eta^2(\T_\ell)
\leq
\sum_{T \in \mathcal{M}_\ell} \eta^2(\T_\ell, T).
\]

Refine $\T_\ell$ to $\T_{\ell + 1}$
by newest-vertex bisection
such that $\mathcal{M}_\ell\subseteq \T_\ell \setminus \T_{\ell+1}$.
sequence of triangulations $\T_\ell$ with $u_\ell$
and $\eta(\T_\ell)$ for $\ell \in \mathbb{N}_0$.

Adaptive C0 interior penalty method

Another main aspect of the numerical investigation
concerns the nonparametric choice of the stabilization
introduced in Section <ref>.

To this end,
the principal eigenvalue $\lambda_1$
to the following eigenvalue problem is computed:
Seek $(\Phi_\IP, \mu) \in V_h \times \R$ such that
$\Phi_\IP \neq 0$ and

\begin{equation}
\label{eq:EVP}
\mu (a_\PW(\Phi_\IP,v_\IP) + c_\IP(\Phi_\IP,v_\IP))
\quad
\text{for all }
v_\IP \in V_h.
\end{equation}

The practical implementation employs
MATLAB's  function
in one line.

[ , lambda] = eigs(B(dof,dof), A(dof,dof) + C(dof,dof), 1, 0);

The Rayleigh-Ritz principle shows that
the principal eigenvalue $\lambda_1$
is an upper bound for the stability constant
\(\kappa = 1/2\) in the case of $a=4$
in Theorem <ref>.

§.§ Singular solution on the L-shaped domain

This benchmark considers
the L-shaped domain $\Omega = (-1,1)^2 \setminus [0,1)^2$ with
interior angle $\omega = 3\pi/2$ at the reentrant corner.
This determines the noncharacteristic solution
$\alpha \coloneqq 0.5444837$ to
$\sin^2(\alpha \omega) = \alpha^2 \sin^2(\omega)$ in

\begin{equation}
\label{eq:angle_component}
\begin{split}
\bigg(
\frac{\sin((\alpha-1)\omega)}{\alpha-1}
- \frac{\sin((\alpha+1)\omega)}{\alpha+1}
\bigg)
\big(
\cos((\alpha-1)\varphi) - \cos((\alpha+1)\varphi)
\big) \\
\bigg(
\frac{\sin((\alpha-1)\varphi)}{\alpha-1}
- \frac{\sin((\alpha+1)\varphi)}{\alpha-1}
\bigg)
\big(
\cos((\alpha-1)\omega) - \cos((\alpha+1)\omega)
\big).
\end{split}
\end{equation}

The exact singular solution in polar coordinates
from <cit.> reads

\begin{align}
\label{eq:singularsolutionLshape}
(1 - r^2 \cos^2\varphi)^2 (1 - r^2 \sin^2\varphi)^2
r^{1+\alpha} g_{\alpha, \omega}\Big(\varphi-\frac{\pi}{2}\Big).
\end{align}

The right-hand side \(f \coloneqq \Delta^2 u\) is computed accordingly.
The reduced regularity of the exact solution $u$
leads to the empirical convergence rate $\alpha/2$
for uniform mesh-refinement displayed by dashed lines in
Figure <ref>.

The adaptive refinement strategy
results in local mesh-refining at the reentrant corner
and recovers the optimal convergence rate $1/2$
with respect to the number of degrees of freedom
(ndof) displayed by the solid lines in
Figure <ref>.

axis equal image,
width = 8cm,
xmin = -1.1,
xmax = 1.1,
ymin = -1.1,
ymax = 1.1

graphics [xmin=-1, xmax=1, ymin=-1, ymax=1]

axis equal image,
width = 0.36,
xmin = -1.15,
xmax = 1.15,
ymin = -1.15,
ymax = 1.15

graphics [xmin=-1, xmax=1, ymin=-1, ymax=1]

Adaptively refined mesh (\(\theta = 0.5\))
with $10406$ degrees of freedom
for the L-shaped domain
from Section <ref>.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.5,
xlabel           = ndof,
ylabel           = error,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos  = outer north east

[col sep=comma]

[col sep=comma]

\(\eta(\mathcal{T}_\ell)\);
\(|| u - u_\ell ||_{\textup{h}}\);

empty legend;

\(\theta = 1\)\(.5\);
\(\theta = 0.5\);

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;
+ [error, uniform, forget plot]
table [x=ndof4lvl, y=err4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;
+ [error, adaptive, forget plot]
table [x=ndof4lvl, y=err4lvl] ;

Convergence history plot
for the benchmark problem on the L-shaped domain
from Section <ref>.

§.§ Singular solution on the 1/8 cusp domain

This benchmark problem investigates the \(1/8\) cusp domain
$\Omega \coloneqq (-1,1)^2 \setminus \conv\{(0,0),(1,-1),(1,0)\}$
with interior angle $\omega = 7\pi / 4$
as depicted in Figure <ref>.
The right-hand side \(f \coloneqq \Delta^2 u\) is given by
the exact singular solution in polar coordinates

\begin{align}
\label{eq:singular_solution_cusp}
(1 - r^2 \cos^2\varphi)^2 (1 - r^2 \sin^2\varphi)^2
r^{1+\alpha} g_{\alpha, \omega}\Big(\varphi - \frac\pi4\Big)
\end{align}

analogously to (<ref>)
with \(g_{\alpha,\omega}\) from (<ref>)
for the parameter $\alpha = 0.50500969$.

The singularity causes an empirical convergence rate $\alpha/2$
for uniform mesh-refinement in Figure <ref>,
while the adaptive algorithm recovers the optimal rate $1/2$.
Undisplayed triangulation plots
confirm the increased adaptive refinement
towards the reentrant corner.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.5,
xlabel           = ndof,
ylabel           = error,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]

[col sep=comma]

\(\eta(\mathcal{T}_\ell)\);
\(|| u - u_\ell ||_{\textup{h}}\);

empty legend;

\(\theta = 1\)\(.0\);
\(\theta = 0.5\);

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;
+ [error, uniform, forget plot]
table [x=ndof4lvl, y=err4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;
+ [error, adaptive, forget plot]
table [x=ndof4lvl, y=err4lvl] ;

Convergence history plot
for the benchmark problem on the
\(1/8\) cusp domain
from Section <ref>.

§.§ Uniform load on the dumbbell-slit domain

The benchmark problem applies the uniform force $f \equiv 1$
on the dumbbell-slit domain

\[
\Omega
\Big(
(-1,1) \times (-1,5)
\setminus
[1,3] \times [-0.75,1)
\Big)
\setminus
(-1,0] \times \{0\}
\]

depicted in Figure <ref>.
Although \(\Omega\) is no Lipschitz domain,
Figures <ref>
and <ref>
confirm the optimal convergence rate with adaptive mesh-refinement and
the stability result from Theorem <ref>.
Undisplayed plots of the adaptive triangulations
show an increased refinement towards the slit
as well as the two reentrant corners at
\((1, -0.75)\) and \((3, -0.75)\).

§.§ Uniform load on the four-slit domain

\[
\Omega
\Big(
(-1,1)^2 \setminus
\big(
[-1,-0.5) \cup [0.5,1)
\big)
\times \{0\}
\Big)
\setminus \{0\} \times ([-1,-0.5) \cup [0.5,1))
\]

denote the square domain with one slit
at each edge as depicted in Figure <ref>.
This benchmark also considers the
uniform load $f \equiv 1$.

Undisplayed triangulation plots indicate that the
adaptive algorithm detects the singularities
of the unknown exact solution
at tips of the four slits.

§.§ Conclusions

While [6] provides a proof of the quasi-optimality
of adaptive discontinuous Galerkin schemes for the diffusion problem,
a rigorous proof of optimal convergence rates
remains open for the biharmonic problem.
Nevertheless, the numerical experiments
in Figure <ref> exhibit
rate-optimal convergence of the
adaptive C0 interior penalty method.
In particular,
even moderate bulk parameters $\theta \leq 0.75$
in the Dörfler marking in Algorithm <ref>
improve the suboptimal convergence rates
for uniform refinement
and lead to the optimal rate of \(1/2\)
in all four benchmark problems.

L-shaped domain from Subsection <ref>.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =
adaptive2/.style =
mark = otimes
adaptive3/.style =
mark = o
adaptive4/.style =
mark = oplus*
adaptive5/.style =
mark = Mercedes star

width            = 0.36,
xlabel           = ndof,
ylabel           = error estimator,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]

[col sep=comma]

[col sep=comma]

[col sep=comma]

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive2, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive3, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

\(1/8\) cusp domain from Subsection <ref>.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =
adaptive2/.style =
mark = otimes
adaptive3/.style =
mark = o
adaptive4/.style =
mark = oplus*
adaptive5/.style =
mark = Mercedes star

width            = 0.36,
xlabel           = ndof,

ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]

[col sep=comma]

[col sep=comma]

[col sep=comma]

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive2, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive3, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

Dumbbell-slit domain from Subsection <ref>.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =
adaptive2/.style =
mark = otimes
adaptive3/.style =
mark = o

width            = 0.36,
xlabel           = ndof,
ylabel           = error estimator,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]

[col sep=comma]

[col sep=comma]

[col sep=comma]

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive2, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive3, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

Four-slit domain from Subsection <ref>.

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =
adaptive2/.style =
mark = otimes
adaptive3/.style =
mark = o

width            = 0.36,
xlabel           = ndof,

ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]ndof4lvl,eta4lvl,err4lvl,lambda4lvl

[col sep=comma]

[col sep=comma]

[col sep=comma]

+ [eta, uniform, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive2, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

+ [eta, adaptive3, forget plot]
table [x=ndof4lvl, y=eta4lvl] ;

eta/.style =
mark = *,
mark size = 3pt,
every mark/.append style =
fill = colEta!60!white

error/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
fill = colError!60!white

eta2/.style =
mark = triangle*,
mark size = 3pt,
every mark/.append style =
fill = colEta2!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =
adaptive2/.style =
mark = otimes
adaptive3/.style =
mark = o
adaptive4/.style =
mark = oplus*
adaptive5/.style =
mark = Mercedes star

hide axis,
legend columns=-1,
legend style=
column sep=1.5ex,

\(\eta(\mathcal{T}_\ell)\);

empty legend;

\(\theta = 1\)\(.75\);
\(\theta = 0.75\);
\(\theta = 0.5\)\(5\);
\(\theta = 0.25\);

Convergence history plot of the error estimator $\eta$
from (<ref>)
for several choices of the bulk parameters \(\theta\).

Figure <ref>
displays the stability constant $\lambda_1$
from (<ref>)
for adaptive and uniform mesh-refinement
in all four benchmark problems.
The computed values for $\lambda_1$ are slightly larger than the guaranteed lower bound \(1/2\) from
Theorem <ref> for \(a = 4\).

The graphs tend towards values between \(0.6\) and \(0.7\)
indicating a small over-stabilization
of the automated choice of $\sigma_{\IP,E}$.

An interesting empirical observation is that adaptive mesh-refinement leads to slightly
more stable discretization with larger stability constants.

The plots in Figure <ref> exhibit no significant difference
between the four domains.
This supports the advantage of the automated choice
of the penalization with $\sigma_{\IP,E}$
for guaranteed stability. The proposed automatic parameter selection appears to be compatible with adaptive mesh-refining. An over-stabilization as in [6] is not necessary in the computational benchmarks.

L-shaped domain from Subsection <ref>.

stability/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
rotate = 45,
fill = colEfficiency!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.38,
xlabel           = ndof,
ylabel           = stability constant,
ymin             = 0.6,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     =  north east

[col sep=comma]

[col sep=comma]

\(\lambda_{1,\ell}\);

empty legend;

\(\theta=1\)\(.5\);
\(\theta = 0.5\);

+ [stability, uniform, forget plot]
table [x=ndof4lvl, y=lambda4lvl];
+ [stability, adaptive, forget plot]
table [x=ndof4lvl, y=lambda4lvl];

\(1/8\) cusp domain from Subsection <ref>.

stability/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
rotate = 45,
fill = colEfficiency!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.38,
xlabel           = ndof,

ymajorgrids      = true,
ymin             = 0.6,
grid style       =
densely dotted,
legend style     =
legend pos  = north east

[col sep=comma]

[col sep=comma]

+ [stability, uniform, forget plot]
table [x=ndof4lvl, y=lambda4lvl];
+ [stability, adaptive, forget plot]
table [x=ndof4lvl, y=lambda4lvl];

Dumbbell-slit domain from Subsection <ref>.

stability/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
rotate = 45,
fill = colEfficiency!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.38,
xlabel           = ndof,
ylabel           = stability constant,
ymin             = 0.6,
ymajorgrids      = true,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]

[col sep=comma]

+ [stability, uniform, forget plot]
table [x=ndof4lvl, y=lambda4lvl];
+ [stability, adaptive, forget plot]
table [x=ndof4lvl, y=lambda4lvl];

Four-slit domain from Subsection <ref>.

stability/.style =
mark = square*,
mark size = 3pt,
every mark/.append style =
rotate = 45,
fill = colEfficiency!60!white

uniform/.style =
every mark/.append style =
fill = black!20!white

adaptive/.style =

width            = 0.38,
xlabel           = ndof,

ymajorgrids      = true,
ymin             = 0.6,
grid style       =
densely dotted,
legend style     =
legend pos     = outer north east

[col sep=comma]ndof4lvl,eta4lvl,err4lvl,lambda4lvl

[col sep=comma]

+ [stability, uniform, forget plot]
table [x=ndof4lvl, y=lambda4lvl];
+ [stability, adaptive, forget plot]
table [x=ndof4lvl, y=lambda4lvl];

Plot of stability constant $\lambda_{1, \ell}$.

§.§ Data availability

Data sharing not applicable as to this article as no datasets were generated or analyzed during the current study.

§ DECLARATIONS

The authors declare that they have no conflict of interests.

[1]
S. Agmon, Lectures on elliptic boundary value problems, AMS Chelsea
Publishing, Providence, RI, 2010.
Prepared for publication by B. Frank Jones, Jr. with the assistance
of George W. Batten, Jr., Revised edition of the 1965 original.

[2]
J. H. Argyris, I. Fried, and D. W. Scharpf, The tuba family of plate
elements for the matrix displacement method, Aeronautical Journal, 72
(1968), pp. 701–709.

[3]
G. A. Baker, Finite element methods for elliptic equations using
nonconforming elements, Math. Comp., 31 (1977), pp. 45–59,
<https://doi.org/10.2307/2005779>, <https://doi.org/10.2307/2005779>.

[4]
J. M. Ball and C. Carstensen, Interaction of martensitic
microstructures in adjacent grains, in Proceedings of the International
Conference on Martensitic Transformations: Chicago, A. P. Stebner and G. B.
Olson, eds., Cham, 2018, Springer International Publishing, pp. 29–33.

[5]
H. Blum and R. Rannacher, On the boundary value problem of the
biharmonic operator on domains with angular corners, Math. Methods Appl.
Sci., 2 (1980), pp. 556–581.

[6]
A. Bonito and R. H. Nochetto, Quasi-optimal convergence rate of an
adaptive discontinuous Galerkin method, SIAM J. Numer. Anal., 48 (2010),
pp. 734–771.

[7]
D. Braess, Finite elements, Cambridge University Press, Cambridge,
third ed., 2007.
Theory, fast solvers, and applications in elasticity theory,
Translated from the German by Larry L. Schumaker.

[8]
S. C. Brenner, $C^0$ interior penalty methods, in Frontiers in
numerical analysis—Durham 2010, vol. 85 of Lect. Notes Comput. Sci. Eng.,
Springer, Heidelberg, 2012, pp. 79–147.

[9]
S. C. Brenner and L. R. Scott, The mathematical theory of finite
element methods, vol. 15 of Texts in Applied Mathematics, Springer, New
York, third ed., 2008.

[10]
S. C. Brenner and L.-Y. Sung, $C^0$ interior penalty methods for
fourth order elliptic boundary value problems on polygonal domains, J. Sci.
Comput., 22/23 (2005), pp. 83–118.

[11]
S. C. Brenner, L.-y. Sung, H. Zhang, and Y. Zhang, A Morley finite
element method for the displacement obstacle problem of clamped Kirchhoff
plates, J. Comput. Appl. Math., 254 (2013), pp. 31–42.

[12]
C. Carstensen, G. Mallik, and N. Nataraj, A priori and a posteriori error control of discontinuous Galerkin finite element
methods for the von Kármán equations, IMA J. Numer. Anal., 39
(2019), pp. 167–200.

[13]
C. Carstensen and N. Nataraj, Adaptive Morley FEM for the von
Kármán equations with optimal convergence rates, SIAM J. Numer.
Anal., 59 (2021), pp. 696–719.

[14]
C. Carstensen and N. Nataraj, A priori and a posteriori error
analysis of the Crouzeix-Raviart and Morley FEM with original and
modified right-hand sides, Comput. Methods Appl. Math., 21 (2021),
pp. 289–315.

[15]
C. Carstensen and N. Nataraj, Lowest-order equivalent nonstandard
finite element methods for biharmonic plates, ESAIM Math. Model. Numer.
Anal., 56 (2022), pp. 41–78.

[16]
D. A. Di Pietro and A. Ern, Mathematical aspects of discontinuous
Galerkin methods, vol. 69, Springer, 2011.

[17]
W. Dörfler, A convergent adaptive algorithm for Poisson's
equation, SIAM J. Numer. Anal., 33 (1996), pp. 1106–1124.

[18]
D. Gallistl, Morley finite element method for the eigenvalues of the
biharmonic operator, IMA J. Numer. Anal., 35 (2015), pp. 1779–1811.

[19]
D. Gilbarg and N. S. Trudinger, Elliptic partial differential
equations of second order, Classics in Mathematics, Springer-Verlag, Berlin,
Reprint of the 1998 edition.

[20]
P. Grisvard, Singularities in boundary value problems, vol. 22 of
Recherches en Mathématiques Appliquées [Research in Applied
Mathematics], Masson, Paris; Springer-Verlag, Berlin, 1992.

[21]
J. Nečas, Direct methods in the theory of elliptic equations,
Springer Monographs in Mathematics, Springer, Heidelberg, 2012.
Translated from the 1967 French original by Gerard Tronel and Alois
Kufner, Editorial coordination and preface by Šárka Nečasová
and a contribution by Christian G. Simader.

[22]
C.-M. Pfeiler and D. Praetorius, Dörfler marking with minimal
cardinality is a linear complexity problem, Math. Comp., 89 (2020),
pp. 2735–2752.

[23]
R. Stevenson, The completion of locally refined simplicial
partitions created by bisection, Math. Comp., 77 (2008), pp. 227–241.

[24]
E. Süli and I. Mozolevski, hp-version interior penalty
DGFEMs for the biharmonic equation, Comput. math. with appl., 196
(2007), pp. 1851–1863.

§ 116 LINES OF MATLAB

This section describes the realization of the computations
from the Sections <ref>,
<ref>, <ref> and <ref>
as a single MATLAB function.
The code has been implemented and tested with MATLAB version
9.9.0.1718557 (R2020b) Update 6.

function [u, ndof] =  C0IP(c4n, n4e, f, a)

if nargin < 4; a = 4; end

nNodes = size(c4n, 1);
nElem = size(n4e, 1);

[sortedSides, localPosition] = ...
sort([n4e(:,[2 3]); n4e(:,[3 1]); n4e(:,[1 2])], 2);
[ , index, position] = unique(sortedSides, 'rows', 'first');
n4s = sortedSides(sub2ind(size(sortedSides), repmat(index, 1, 2),...
nSides = size(n4s, 1);
s4e = reshape(position, nElem, 3);
ElemNo = repmat((1:nElem)', 3, 1);
e4s = [ElemNo(index), accumarray(position, ElemNo) - ElemNo(index)];
isInteriorSide = (e4s(:,2)  = 0);

tangent4s = c4n(n4s(:,2),:) - c4n(n4s(:,1),:);
length4s = sqrt(sum(tangent4s.^2, 2));
tangent4s = tangent4s ./ length4s;
normal4s = [tangent4s(:,2), - tangent4s(:,1)];

c4e = permute(reshape(c4n(n4e,:), nElem, 3, 2), [3 2 1]);
area4e = squeeze(abs(sum((c4e(2,[2 3 1],:) - c4e(2,[3 1 2],:)) ...
.* c4e(1,:,:), 2))) / 2;

dof4e = [n4e, nNodes + s4e]';

[oppositeIndPlus, ] = ...
find(all(n4e(e4s(:,1),:)  = permute(n4s, [1 3 2]), 3)');
[oppositeIndMinus, ] = ...
find(all(n4e(e4s(isInteriorSide,2),:) ...
= permute(n4s(isInteriorSide,:), [1 3 2]), 3)');

LOCAL = [2 4 3 5 1 6; 3 5 1 6 2 4; 1 6 2 4 3 5]';
localInd4s = ones(6, 2, nSides);
localInd4s(:,1,:) = LOCAL(:,oppositeIndPlus);
localInd4s(:,2,isInteriorSide) = LOCAL(:,oppositeIndMinus);

dof4s = [n4s(:,1), nNodes + (1:nSides)', n4s(:,2), ones(nSides, 6)]';
dof4s(5,:) = n4e(sub2ind(size(n4e), e4s(:,1), oppositeIndPlus));
dof4s([4 6],:) = nNodes + ...
s4e(sub2ind(size(s4e), repmat(e4s(:,1), [1 2]), ...
LOCAL([1 3], oppositeIndPlus)'))';
dof4s(8,isInteriorSide) = ...
n4e(sub2ind(size(n4e), e4s(isInteriorSide,2), oppositeIndMinus));
dof4s([7 9],isInteriorSide) = nNodes + ...
s4e(sub2ind(size(s4e), repmat(e4s(isInteriorSide,2), [1 2]), ...
LOCAL([1 3],oppositeIndMinus)'))';

f4mid = f((c4n(n4s(:,1),:) + c4n(n4s(:,2),:)) / 2);
valF4e = f4mid(s4e);

Grad4e = zeros(6, 2, 6, nElem);

Hess4e = zeros(6, 2, 2, nElem);

A4e = zeros(6, 6, nElem);

b4e = zeros(6, nElem);

parfor e = 1:nElem

GradP1 = [ones(1,3); c4e(:,:,e)]  [zeros(1,2); eye(2)];

H = permute(GradP1, [2 3 1 4]) .* permute(GradP1, [3 2 4 1]);
Hessian = 4 * cat(3, H(:,:,1,1), H(:,:,2,2), H(:,:,3,3), ...
H(:,:,2,3) + H(:,:,3,2), ...
H(:,:,1,3) + H(:,:,3,1), ...
H(:,:,1,2) + H(:,:,2,1));

Integrand = sum(Hessian .* permute(Hessian, [1 2 4 3]), [1 2]);
A4e(:,:,e) = area4e(e) * squeeze(Integrand);

b4e(:,e) = [zeros(1, 3), area4e(e) / 3 * valF4e(j,:)];

Hess4e(:,:,:,e) = permute(Hessian, [3 1 2]);
Gradient = cat(3, [3 -1 -1 -1 1 1] .* GradP1(1,:)', ...
[-1 3 -1 1 -1 1] .* GradP1(2,:)', ...
[-1 -1 3 1 1 -1] .* GradP1(3,:)', ...
[0 0 4 2 2 0] .* GradP1(2,:)' ...
+ [0 4 0 2 0 2] .* GradP1(3,:)', ...
[0 0 4 2 2 0] .* GradP1(1,:)' ...
+ [4 0 0 0 2 2] .* GradP1(3,:)', ...
[0 4 0 2 0 2] .* GradP1(1,:)' ...
+ [4 0 0 0 2 2] .* GradP1(2,:)');
Grad4e(:,:,:,e) = permute(Gradient, [3 1 2]);

Grad4s = zeros(6, 2, 6, 2, nSides);
Grad4s(:,:,:,1,:) = Grad4e(:,:,:,e4s(:,1));
Grad4s(:,:,:,2,isInteriorSide) = Grad4e(:,:,:,e4s(isInteriorSide,2));
Hess4s = zeros(6, 2, 2, 2, nSides);
Hess4s(:,:,:,1,:) = Hess4e(:,:,:,e4s(:,1));
Hess4s(:,:,:,2,isInteriorSide) = Hess4e(:,:,:,e4s(isInteriorSide,2));

sigma4s = 3 * a * (length4s.^2) ./ area4e(e4s(:,1));
sigma4s(isInteriorSide) = sigma4s(isInteriorSide) / 4 ...
+ 3/4 * a * length4s(isInteriorSide).^2 ...
./ area4e(e4s(isInteriorSide,2));

J4s = zeros(9, 9, nSides);

C4s = zeros(9, 9, nSides);

parfor s = 1:nSides

Ind = localInd4s(:,:,s);

GradNormal = squeeze(sum(Grad4s(:,:,:,:,s) .* normal4s(s,:), 2));

JumpNormal = [GradNormal(Ind(1:3,1),Ind(1:3,1),1) ...
- GradNormal(Ind(3:-1:1,2),Ind(3:-1:1,2),2); ...
GradNormal(Ind(4:6,1),Ind(1:3,1),1); ...
- GradNormal(Ind(4:6,2),Ind(3:-1:1,2),2)];

Integrand = JumpNormal .* permute(JumpNormal, [3 2 1]);
C4s(:,:,s) = sigma4s(s) / 6 * squeeze(sum([1 4 1] .* Integrand, 2));

HessNormal = squeeze(sum(Hess4s(:,:,:,:,s) .* normal4s(s,:), 2));

MeanHessNormal = ...
[HessNormal(Ind(1:3,1),:,1) + HessNormal(Ind(3:-1:1,2),:,2);...
HessNormal(Ind(4:6,1),:,1); HessNormal(Ind(4:6,2),:,2)];
if isInteriorSide(s); MeanHessNormal = MeanHessNormal / 2; end

JumpGrad = JumpNormal(:,2) * normal4s(s,:);

J4s(:,:,s) = length4s(s) * (MeanHessNormal * JumpGrad');

RowIndices4e = permute(repmat(dof4e, [1 1 size(dof4e, 2)]), [1 3 2]);
ColIndices4e = permute(RowIndices4e, [2 1 3]);

RowIndices4s = permute(repmat(dof4s, [1 1 size(dof4s, 2)]), [1 3 2]);
ColIndices4s = permute(RowIndices4s, [2 1 3]);

N = nNodes + nSides;
A = sparse(RowIndices4e(:), ColIndices4e(:), A4e(:), N, N);
J = sparse(RowIndices4s(:), ColIndices4s(:), J4s(:), N, N);
C = sparse(RowIndices4s(:), ColIndices4s(:), C4s(:), N, N);
B = A - J - J' + C;
b = accumarray(dof4e(:), b4e(:), [N 1]);

n4sCb = n4s( isInteriorSide,:);
CbNodes = unique(n4sCb(:));
[ , CbSides] = ismember(sort(n4sCb, 2), sort(n4s, 2), 'rows');
dof = setdiff(1:N, [CbNodes, nNodes + CbSides]');
ndof = size(dof, 2);

u = zeros(N, 1);
u(dof) = B(dof,dof)  b(dof);

§ A POSTERIORI ERROR ESTIMATION

The estimator
\(\eta^2(\T) \coloneqq \sum_{T \in \T} \eta^2(T)\)
with, for all \(T \in \T\),

\begin{equation}
\label{eq:estimator}
\begin{split}
\eta^2(T)
\Vert h_T^2 f \Vert_{L^2(T)}^2
\sum_{E \in \E(T)}
\frac{\sigma_{\IP,E}^2}{h_E}
\Vert [ \nabla u_\IP ]_E \cdot \nu_E \Vert_{L^2(E)}^2\\
\sum_{E \in \E(T) \cap \mathcal{E}(\Omega)}
\Vert [\partial_{\nu_E}^2 u_\IP]_E \Vert_{L^2(\Omega)}^2
\end{split}
\end{equation}

from <cit.>
can be evaluated using some of the data structures
from the previous appendix Section <ref>.

vol4e = (area4e.^3) .* sum(valF4e.^2, 2) / 3;
u4e = u(dof4e');
u4s = u(dof4s');
penalty4s = sigma4s .* squeeze(sum(C4s .* permute(u4s, [2 3 1]) ...
.* permute(u4s, [3 2 1]), [1 2]));
tangents4e = (c4e(:,[3 1 2],:) - c4e(:,[2 3 1],:)) ...
./ permute(length4s(s4e), [3 2 1]);
normal4e = [tangents4e(2,:,:); - tangents4e(1,:,:)];
signedNormal4e = ...
permute(reshape(normal4s(s4e(:),:), [nElem 3 2]), [3 2 1]);
HessU4e = squeeze(sum(Hess4e .* permute(u4e, [2 3 4 1]), 1));
BinormalU4e = ...
sum(permute(normal4e, [1 4 3 2]) .* HessU4e ...
.* permute(signedNormal4e, [4 1 3 2]), [1 2]));
jump4s = accumarray(s4e(:), biNormalU4e(:), [nSides 1]);
jump4s = (length4s).^2 .* jump4s.^2;
jump4s(CbSides) = 0;
eta4e = vol4e + sum(jump4s(s4e), 2) + sum(penalty4s(s4e), 2);

§ DOWNLOADS

The MATLAB function
is available online and
can be downloaded from the URL
<https://www.math.hu-berlin.de/ cc/cc_homepage/software/code/2022-PB_CC_JS-Parameter_free_C0IP/C0IP.m>