---
title: Adaptive Crouzeix-Raviart FEM
url: https://www.emergentmind.com/topics/adaptive-crouzeix-raviart-finite-element-method
type: topic
---

# Adaptive Crouzeix-Raviart FEM

Searching arXiv for the cited papers and closely related adaptive Crouzeix–Raviart FEM work.
Adaptive Crouzeix-Raviart finite element methods are adaptive finite element methods built on the nonconforming Crouzeix-Raviart discretization, in which trial and test functions are piecewise polynomial on each simplex and satisfy continuity only in the sense of vanishing integral jumps across faces or edges. In the form analyzed for the Poisson and incompressible Stokes problems in 2D, the method combines newest vertex bisection, side-based a posteriori estimators, and a modified maximum marking strategy, and yields instance optimality. Subsequent work develops closely related adaptive CR frameworks for elliptic eigenvalue problems, convection-diffusion eigenvalue problems, distributed optimal control governed by Stokes equations, the first eigenpair of the \(p\)-Laplacian, and high-order odd-degree \(\mathrm{CR}_k\) discretizations [1404.3065; 1905.08243; 1606.03171; 2307.00555; 2508.02077; 2602.16588].

## 1. Problem classes, meshes, and nonconforming spaces

For the Poisson model on a bounded polygonal domain \(\Omega \subset \mathbb{R}^2\), the continuous problem is: find \(u \in H^1_0(\Omega)\) such that
\[
a(u,v) := \int_\Omega \nabla u \cdot \nabla v \, dx = (f,v) := \int_\Omega f v \, dx, \quad \forall v \in H^1_0(\Omega).
\]
For the incompressible Stokes problem in 2D, one seeks \((u,p)\) with \(u \in [H^1_0(\Omega)]^2\) and \(p \in L^2_0(\Omega)\) such that
\[
\int_\Omega \nabla u : \nabla v \, dx + \int_\Omega p \operatorname{div} v\, dx + \int_\Omega q \operatorname{div} u\, dx = \int_\Omega f \cdot v\, dx,
\]
for all \((v,q) \in [H^1_0(\Omega)]^2 \times L^2_0(\Omega)\), with homogeneous Dirichlet boundary conditions and \(\nu=1\) implicit. The mesh family \(\mathbb{T}\) consists of conforming triangulations obtained from an initial triangulation \(\mathcal{T}_\perp\) by newest vertex bisection, and the resulting meshes remain uniformly shape-regular while enjoying partial ordering, overlay, and closure properties [1404.3065].

On a triangulation \(\mathcal{T}\), the lowest-order CR space is
\[
\mathrm{CR}(\mathcal{T}) := \left\{ v \in L^2(\Omega) : v|_T \in P_1(T)\ \forall T\in\mathcal{T},\ \int_S \llbracket v \rrbracket\, ds = 0\ \forall S \in \mathcal{S}(\mathcal{T}) \right\}.
\]
Its degrees of freedom are mean values over edges, and boundary conditions are incorporated via edge averages on boundary edges. The nonconforming gradient is defined elementwise by
\[
(\nabla_{\mathrm{NC}} v)|_{T} := \nabla (v|_{T}), \quad \forall T \in \mathcal{T},
\]
so the broken gradient is piecewise constant on each triangle. For Stokes, the discrete velocity and pressure spaces are
\[
V_{\mathrm{CR}}(\mathcal{T}) := \mathrm{CR}(\mathcal{T}) \times \mathrm{CR}(\mathcal{T}), \qquad
Q_0(\mathcal{T}) := \left\{ q \in L^2_0(\Omega) : q|_T \in P_0(T),\ \forall T\in \mathcal{T} \right\},
\]
with discrete divergence
\[
\operatorname{div}_{\mathrm{NC}} v := \operatorname{tr} \nabla_{\mathrm{NC}} v \in Q_0(\mathcal{T}),
\]
which ensures exact local divergence-freeness. The CR interpolation operator \(I_{\mathcal{T}}\) is defined by edge averages,
\[
\int_S (I_{\mathcal{T}} v)\, ds = \int_S v\, ds,\quad \forall S\in\mathcal{S}(\mathcal{T}),
\]
and satisfies the projection property that \(\nabla_{\mathrm{NC}}(I_{\mathcal{T}} v)\) is the \(L^2\)-projection of \(\nabla v\) into piecewise constants, together with the approximation and stability estimate recorded as equation (3.6) in the paper [1404.3065].

The same nonconforming design recurs in other problem classes. For elliptic eigenvalue problems, the CR space \(V_h^{CR}\) is defined on a shape-regular triangulation by piecewise \(P_1\) functions that are continuous in average across interior edges and have zero edge average on the boundary; the broken bilinear form is
\[
a_h(w_h, v_h) := \sum_{K \in \mathcal{T}_h} \int_K \nabla_h w_h \cdot \nabla_h v_h.
\]
The commuting property
\[
\int_K \nabla(w - \Pi_h w)\cdot \nabla v_h = 0
\]
for each \(K\) is essential for eigenvalue error analysis [1905.08243].

## 2. Residual-based a posteriori estimation

A defining feature of adaptive CR methods is that the estimators are tailored to nonconformity. For the Poisson problem, the estimator in [1404.3065] is side-based rather than element-based:
\[
\mathcal{E}_{\mathcal{T}}(S) := \| h_{\mathcal{T}} f \|_{L^2(\omega_S)}^2 + h_S\ \| \llbracket \nabla_{\mathrm{NC}} u_{\mathcal{T}} \cdot t \rrbracket \|_{L^2(S)}^2,
\]
with global quantity
\[
\mathcal{E}_{\mathcal{T}}^2(\mathcal{S}) := \sum_{S \in \mathcal{S}} \mathcal{E}_{\mathcal{T}}(S).
\]
The oscillation term is
\[
\operatorname{osc}(\mathcal{T})^2 := \sum_{T \in \mathcal{T}} \int_T h_T^2 |f - f_T|^2\, dx,
\qquad
f_T = \frac{1}{|T|} \int_T f\, dx,
\]
and reliability and efficiency are expressed by
\[
\| \nabla_{\mathrm{NC}} (u - u_{\mathcal{T}}) \|_{L^2(\Omega)}^2 + \operatorname{osc}(\mathcal{T})^2 \lesssim \mathcal{E}_{\mathcal{T}}^2(\mathcal{S}(\mathcal{T})) \lesssim \| \nabla_{\mathrm{NC}} (u - u_{\mathcal{T}}) \|_{L^2(\Omega)}^2 + \operatorname{osc}(\mathcal{T})^2.
\]
For Stokes, the side indicator has the same structure, now applied componentwise to the velocity and coupled to the pressure error:
\[
\| \nabla_{\mathrm{NC}} (u - u_{\mathcal{T}}) \|_{L^2}^2 + \| p - p_{\mathcal{T}} \|_{L^2}^2 + \operatorname{osc}(\mathcal{T})^2 \lesssim \mathcal{E}_{\mathcal{T}}^2(\mathcal{S}(\mathcal{T})) \lesssim \| \nabla_{\mathrm{NC}} (u - u_{\mathcal{T}}) \|_{L^2(\Omega)}^2 + \| p - p_{\mathcal{T}} \|_{L^2(\Omega)}^2 + \operatorname{osc}(\mathcal{T})^2
\]
[1404.3065].

Later CR estimators preserve this residual-plus-jump pattern but adapt it to the underlying PDE. For elliptic eigenvalue problems, local indicators are given by
\[
\eta_K^2 = h_K^2 \| \Delta_h u_h \|^2_{0, K} + \sum_{e\subset \partial K} h_e \left\| \Big[ \frac{\partial u_h}{\partial n} \Big] \right\|^2_{0, e} + h_e \left\| \Big[ \frac{\partial u_h}{\partial t} \Big] \right\|^2_{0,e},
\]
while the paper also develops two asymptotically exact global estimators \(F^{CR}_{CR,1}\) and \(F^{CR}_{CR,2}\) for the eigenvalue error, both driven by recovered gradients and computable approximations of the consistency error term [1905.08243]. For convection-diffusion eigenvalue problems, the estimator includes an element residual
\[
\eta_{h,K} := h_K \left\| \lambda_h u_h + \Delta_h u_h - \mathbf{b}\cdot\nabla_h u_h \right\|_{L^2(K)}
\]
together with normal and tangential jump terms on edges, and analogous quantities for the dual eigenproblem [1606.03171]. For distributed optimal control governed by Stokes equations, the global estimator combines state and adjoint residuals, jump terms, and oscillation, and the estimator is both reliable and efficient up to oscillation for the coupled state-adjoint-pressure-control error [2307.00555]. For the first eigenpair of the \(p\)-Laplacian, the estimator separates a volume contribution
\[
\eta_{h,1}(\mu, v, K) := \mu^q h_K^q \| v \|_{L^p(K)}^p
\]
from a nonconformity contribution
\[
\eta_{h,2}(v, K) := \sum_{E \subset \partial K \cap \Omega} \frac{1}{2} h_E \| [v] \|_{L^p(E)}^p + \sum_{E \subset \partial K \cap \partial\Omega} h_E \| [v] \|_{L^p(E)}^p,
\]
and for high-order odd-degree \(\mathrm{CR}_k\) the residual-type estimator contains both normal and tangential projections of gradient jumps [2508.02077; 2602.16588].

## 3. Adaptive loop and marking strategies

The canonical adaptive cycle is
\[
\texttt{SOLVE} \longrightarrow \texttt{ESTIMATE} \longrightarrow \texttt{MARK} \longrightarrow \texttt{REFINE}.
\]
In the instance-optimal analysis for Poisson and Stokes, the decisive point is the modified maximum marking strategy. For a marking parameter \(\theta \in (0,1]\), the contribution attached to a current side \(S \in \mathcal{S}(\mathcal{T}_k)\) is not simply \(\mathcal{E}_{\mathcal{T}_k}(S)\), but
\[
\sum_{S' \in \operatorname{ref}'d(\mathcal{T}_k; S)} \mathcal{E}_{\mathcal{T}_k}(S'),
\]
where \(\operatorname{ref}'d(\mathcal{T}_k; S)\) denotes the set of sides that must be bisected with \(S\) to maintain conformity. One marks all sides such that
\[
\sum_{S' \in \operatorname{ref}'d(\mathcal{T}_k; S)} \mathcal{E}_{\mathcal{T}_k}(S') \geq \theta \cdot \max_{S'' \in \mathcal{S}(\mathcal{T}_k)} \sum_{S''' \in \operatorname{ref}'d(\mathcal{T}_k; S'')} \mathcal{E}_{\mathcal{T}_k}(S'''),
\]
and then refines by newest vertex bisection. The modification is necessary because standard maximum marking on nonconforming CR spaces does not ensure the estimator is reduced enough when conformity constraints force additional refinements [1404.3065].

A recurrent source of confusion is the role of conformity in the marking stage. In CR AFEM, the object marked for reduction may be a side, an element, or an aggregate attached to a refinement dependency, but the analytical reason is the same: local error reduction must be compatible with the refinement closure generated by NVB. This is why the marking criterion in [1404.3065] accumulates all sides that must be marked together.

In later CR adaptive schemes, Dörfler or bulk-chasing strategies are also used. For elliptic eigenvalue problems, the algorithm marks elements with largest \(\eta_K\), or all \(K\) such that \(\eta_K > \theta \max_{K'} \eta_{K'}\), with \(\theta = 0.3\) given as an example [1905.08243]. For convection-diffusion eigenvalue problems, the marked subset \(\widehat{\mathcal{T}_{h_l}}\) satisfies
\[
\sum_{K\in\widehat{\mathcal{T}_{h_l}} \left( \eta_{h_l}(K)^2 + \eta_{h_l}^*(K)^2 \right) \geq \theta \sum_{K\in\mathcal{T}_{h_l}} \left( \eta_{h_l}(K)^2 + \eta_{h_l}^*(K)^2 \right),
\]
so both primal and dual errors drive refinement [1606.03171]. For distributed Stokes control and for high-order \(\mathrm{CR}_k\), Dörfler marking is likewise used, whereas the \(p\)-Laplacian paper adopts maximum marking and notes that Dörfler’s bulk marking with parameter \(\theta \in (0,1)\) is possible [2307.00555; 2602.16588; 2508.02077].

## 4. Analytical structure: non-nestedness, quasi-orthogonality, and optimality

The central analytical difficulty of adaptive CR methods is that the discrete spaces are not nested under refinement:
\[
\mathrm{CR}(\mathcal{T}) \not\subset \mathrm{CR}(\mathcal{T}_+).
\]
As a consequence, standard Galerkin orthogonality is absent. The analysis therefore replaces it by quasi-orthogonality, transfer operators, and quasi-best approximation results [1404.3065].

In the Poisson and Stokes theory of [1404.3065], the framework is an adaptation of Diening-Kreuzer-Stevenson to the nonconforming setting. Its basic ingredients are a generalized energy \(\mathcal{G}(\mathcal{T})\), monotonic under refinement; localized energy differences equivalent to quasi-error; a lower diamond estimate supporting overlay arguments; quasi-best approximation; estimator reduction and stability; and linear complexity. For Poisson, the energy difference satisfies
\[
\mathcal{G}(\mathcal{T}) - \mathcal{G}(\mathcal{T}_*) \simeq \| \nabla_{\mathrm{NC}} (u_{\mathcal{T}_*} - u_{\mathcal{T}}) \|_{L^2(\Omega)}^2 + \| h_{\mathcal{T}} f \|_{L^2(\Omega \setminus \mathcal{T}_*)}^2,
\]
and the main instance-optimal estimate is
\[
\| \nabla_{\mathrm{NC}}(u-u_{\mathcal{T}_k}) \|_{L^2(\Omega)}^2 + \operatorname{osc}(\mathcal{T}_k)^2
\le C\, \inf_{\mathcal{T}\,:\,\#\mathcal{T}\le \#\mathcal{T}_k}
\Big( \| \nabla_{\mathrm{NC}}(u-u_{\mathcal{T}}) \|_{L^2(\Omega)}^2 + \operatorname{osc}(\mathcal{T})^2 \Big).
\]
For Stokes, the pressure term is added:
\[
\| \nabla_{\mathrm{NC}}(u-u_{\mathcal{T}_k}) \|_{L^2}^2+ \|p - p_{\mathcal{T}_k}\|_{L^2}^2 + \operatorname{osc}(\mathcal{T}_k)^2
\le C\, \inf_{\mathcal{T}\,:\, \#\mathcal{T}\le \#\mathcal{T}_k}
\Big( \|\nabla_{\mathrm{NC}}(u-u_{\mathcal{T}})\|_{L^2}^2 + \|p - p_{\mathcal{T}}\|_{L^2}^2 + \operatorname{osc}(\mathcal{T})^2 \Big).
\]
Here \(\#\mathcal{T}\) denotes the number of elements, and the constants are independent of the data but may depend on \(\theta\) [1404.3065].

A broader abstract viewpoint appears in later work through the axioms of adaptivity: stability, reduction, discrete reliability, and quasi-orthogonality. For distributed optimal control governed by Stokes equations, these axioms yield quasi-optimal convergence rates of the adaptive algorithm. For high-order odd-degree \(\mathrm{CR}_k\), the principal new theoretical result is the proof of Axiom 3, discrete reliability, using new local quasi-interpolation operators, companion operators, and an intersection operator, with
\[
\|\nabla_{\widehat{\mathcal{T}}} \widehat{u}_h - \nabla_{\mathcal{T}} u_h \|_{L^2(\Omega)} \lesssim \eta(u_h; \mathcal{R}^1, f)
\]
on a first-layer patch of refined elements [2307.00555; 2602.16588].

## 5. Specialized adaptive CR frameworks

For elliptic eigenvalue problems, the nonconformity creates a consistency error term in the exact identity
\[
\lambda - \lambda_h = |u - u_h|_{1,h}^2 + 2\left(a_h(u, u_h) - \lambda_h (u, u_h)\right) - \lambda_h \|u - u_h\|^2_{0,\Omega}.
\]
The paper [1905.08243] addresses this by combining high-accuracy gradient recovery with computable approximations of the consistency term. It proves asymptotic exactness for two eigenvalue estimators:
\[
\lim_{h\to 0} \frac{ |\lambda - \lambda_h|}{|F^{CR}_{CR,i}| } = 1, \qquad i=1,2,
\]
and introduces a postprocessed eigenvalue
\[
\lambda_{pp} = \frac{F_1}{F_1+F_2}\, \lambda^{(2)} + \frac{F_2}{F_1+F_2}\, \lambda^{(1)},
\]
whose weights are entirely computed from a posteriori estimators and which requires only one eigenvalue problem to be solved.

For convection-diffusion eigenvalue problems, the CR method is applied to the non-selfadjoint operator \(\mathcal{L}u := -\Delta u + \mathbf{b}\cdot\nabla u\), so both right and left eigenproblems are included in the adaptive loop. Reliability and efficiency are established for the primal and dual estimators, and the marking uses their sum. Numerical results reported in the source show that the estimator successfully concentrates the mesh in layer regions created by strong convection, while the paper also notes deterioration for \(\mathbf{b}=(10,0)^T\) because no upwinding or extra stabilization is used [1606.03171].

For distributed optimal control governed by the Stokes equations, the nonconforming lowest order Crouzeix-Raviart element and piecewise constant spaces discretize the velocity and pressure, while the control is handled by variational discretization through
\[
u_h = -\alpha^{-1} z_h,
\]
or by the projected form in the constrained case. The resulting estimator combines state and adjoint residuals, edge jumps, and oscillation, and the paper establishes error equivalence at both continuous and discrete levels together with quasi-optimal convergence rates under the general axiomatic framework [2307.00555].

For the first Dirichlet eigenpair of the \(p\)-Laplacian, the adaptive CR method works on \(d=2,3\) polyhedral domains and proves that the estimator has vanishing limit,
\[
\lim_{k \to \infty} \eta_{k,1}(\mu_k, u_k) + \eta_{k,2}(u_k) = 0,
\]
that the discrete eigenvalues converge,
\[
\lim_{k \to \infty} \mu_k = \lambda_1,
\]
and that the discrete eigenfunctions converge to the eigenspace in the mesh-dependent broken norm. The key novelty is a discrete compactness result for CR finite elements on adaptively generated meshes, driven by the vanishing of the jump term [2508.02077].

For arbitrary odd degree \(k\geq 1\), high-order \(\mathrm{CR}_k\) adaptive methods for the 2D Poisson problem retain the SOLVE–ESTIMATE–MARK–REFINE structure and a residual-type estimator. The analysis is confined to odd \(k\) because for even \(k\) the degrees of freedom are not strictly local and require global constraints. Numerical experiments reported in [2602.16588] illustrate optimal convergence rates for all considered values of \(k\).

## 6. Scope, limitations, and recurring distinctions

The classical instance-optimal theory in [1404.3065] is proved in 2D only, for polygonal domains and right-hand sides \(f \in L^2(\Omega)\) for Poisson and \(f \in [L^2(\Omega)]^2\) for Stokes. The paper explicitly states that extension to 3D is open because the population tree structure and combinatorial ingredients are not established in higher dimensions. It also states that the analysis is robust for piecewise constant coefficients with modifications, while smooth non-piecewise constant coefficients require additional oscillation terms and may not admit instance optimality with current methods [1404.3065].

Several broader limitations recur across the literature. Elliptic eigenvalue asymptotic exactness in [1905.08243] assumes \(u \in H^{7/2}(\Omega)\) for the sharpest results, and the second type estimator is noted to be more robust for nonsmooth eigenfunctions. The convection-diffusion eigenvalue paper does not discuss clustered eigenvalues or complexity optimality and emphasizes that no upwind stabilization is used [1905.08243; 1606.03171]. The \(p\)-Laplacian analysis requires a sufficiently fine initial mesh so that adaptive refinement targets the correct first eigenpair, and full rate and complexity theory remain open [2508.02077]. The high-order \(\mathrm{CR}_k\) theory is restricted to 2D Poisson and odd polynomial degrees [2602.16588].

Two distinctions are especially important. First, adaptive CR methods are nonconforming but not arbitrary: the loss of nestedness forces the use of quasi-orthogonality, transfer operators, or compactness arguments instead of standard Galerkin orthogonality. Second, CR estimators are not confined to normal derivative jumps. In the nonconforming setting, tangential jumps play a structural role, and in several formulations both normal and tangential components enter the estimator. A plausible implication is that the analytic and algorithmic identity of adaptive CR methods lies less in a single estimator formula than in the combination of nonconforming spaces, refinement-closure-aware marking, and proof techniques that compensate for non-nestedness [1404.3065; 1606.03171; 2602.16588].

Source: https://www.emergentmind.com/topics/adaptive-crouzeix-raviart-finite-element-method