Adaptive Crouzeix-Raviart FEM

• Adaptive Crouzeix-Raviart FEM is a nonconforming method using piecewise polynomial functions with edge-average continuity to discretize PDEs efficiently.
• It employs residual-based a posteriori error estimators and modified marking strategies to achieve instance optimal convergence for Poisson, Stokes, and eigenvalue problems.
• The method overcomes non-nestedness by leveraging quasi-orthogonality and transfer operators, broadening its scope to convection-diffusion and optimal control applications.

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 [(Kreuzer et al., 2014); (Hu et al., 2019); (Du et al., 2016); (Shaikh et al., 2023); (Li et al., 4 Aug 2025); (Bohne et al., 18 Feb 2026)].

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 $\mathrm{CR}_k$0 implicit. The mesh family $\mathrm{CR}_k$1 consists of conforming triangulations obtained from an initial triangulation $\mathrm{CR}_k$2 by newest vertex bisection, and the resulting meshes remain uniformly shape-regular while enjoying partial ordering, overlay, and closure properties (Kreuzer et al., 2014).

On a triangulation $\mathrm{CR}_k$3, the lowest-order CR space is

$\mathrm{CR}_k$4

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

$\mathrm{CR}_k$5

so the broken gradient is piecewise constant on each triangle. For Stokes, the discrete velocity and pressure spaces are

$\mathrm{CR}_k$6

with discrete divergence

$\mathrm{CR}_k$7

which ensures exact local divergence-freeness. The CR interpolation operator $\mathrm{CR}_k$8 is defined by edge averages,

$\mathrm{CR}_k$9

and satisfies the projection property that $\Omega \subset \mathbb{R}^2$0 is the $\Omega \subset \mathbb{R}^2$1-projection of $\Omega \subset \mathbb{R}^2$2 into piecewise constants, together with the approximation and stability estimate recorded as equation (3.6) in the paper (Kreuzer et al., 2014).

The same nonconforming design recurs in other problem classes. For elliptic eigenvalue problems, the CR space $\Omega \subset \mathbb{R}^2$3 is defined on a shape-regular triangulation by piecewise $\Omega \subset \mathbb{R}^2$4 functions that are continuous in average across interior edges and have zero edge average on the boundary; the broken bilinear form is

$\Omega \subset \mathbb{R}^2$5

The commuting property

$\Omega \subset \mathbb{R}^2$6

for each $\Omega \subset \mathbb{R}^2$7 is essential for eigenvalue error analysis (Hu et al., 2019).

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 (Kreuzer et al., 2014) is side-based rather than element-based: $\Omega \subset \mathbb{R}^2$8 with global quantity

$\Omega \subset \mathbb{R}^2$9

The oscillation term is

$u \in H^1_0(\Omega)$0

and reliability and efficiency are expressed by

$u \in H^1_0(\Omega)$1

For Stokes, the side indicator has the same structure, now applied componentwise to the velocity and coupled to the pressure error: $u \in H^1_0(\Omega)$2 (Kreuzer et al., 2014).

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

$u \in H^1_0(\Omega)$3

while the paper also develops two asymptotically exact global estimators $u \in H^1_0(\Omega)$4 and $u \in H^1_0(\Omega)$5 for the eigenvalue error, both driven by recovered gradients and computable approximations of the consistency error term (Hu et al., 2019). For convection-diffusion eigenvalue problems, the estimator includes an element residual

$u \in H^1_0(\Omega)$6

together with normal and tangential jump terms on edges, and analogous quantities for the dual eigenproblem (Du et al., 2016). 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 (Shaikh et al., 2023). For the first eigenpair of the $u \in H^1_0(\Omega)$7-Laplacian, the estimator separates a volume contribution

$u \in H^1_0(\Omega)$8

from a nonconformity contribution

$u \in H^1_0(\Omega)$9

and for high-order odd-degree $$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).$$0 the residual-type estimator contains both normal and tangential projections of gradient jumps (Li et al., 4 Aug 2025, Bohne et al., 18 Feb 2026).

3. Adaptive loop and marking strategies

The canonical adaptive cycle is

$$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).$$1

In the instance-optimal analysis for Poisson and Stokes, the decisive point is the modified maximum marking strategy. For a marking parameter $$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).$$2, the contribution attached to a current side $$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).$$3 is not simply $$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).$$4, but

$$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).$$5

where $$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).$$6 denotes the set of sides that must be bisected with $$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).$$7 to maintain conformity. One marks all sides 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).$$8

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 (Kreuzer et al., 2014).

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 (Kreuzer et al., 2014) 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 $$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).$$9, or all $(u,p)$0 such that $(u,p)$1, with $(u,p)$2 given as an example (Hu et al., 2019). For convection-diffusion eigenvalue problems, the marked subset $(u,p)$3 satisfies

$(u,p)$4

so both primal and dual errors drive refinement (Du et al., 2016). For distributed Stokes control and for high-order $(u,p)$5, Dörfler marking is likewise used, whereas the $(u,p)$6-Laplacian paper adopts maximum marking and notes that Dörfler’s bulk marking with parameter $(u,p)$7 is possible (Shaikh et al., 2023, Bohne et al., 18 Feb 2026, Li et al., 4 Aug 2025).

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: $(u,p)$8 As a consequence, standard Galerkin orthogonality is absent. The analysis therefore replaces it by quasi-orthogonality, transfer operators, and quasi-best approximation results (Kreuzer et al., 2014).

In the Poisson and Stokes theory of (Kreuzer et al., 2014), the framework is an adaptation of Diening-Kreuzer-Stevenson to the nonconforming setting. Its basic ingredients are a generalized energy $(u,p)$9, 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

$u \in [H^1_0(\Omega)]^2$0

and the main instance-optimal estimate is

$u \in [H^1_0(\Omega)]^2$1

For Stokes, the pressure term is added: $u \in [H^1_0(\Omega)]^2$2 Here $u \in [H^1_0(\Omega)]^2$3 denotes the number of elements, and the constants are independent of the data but may depend on $u \in [H^1_0(\Omega)]^2$4 (Kreuzer et al., 2014).

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 $u \in [H^1_0(\Omega)]^2$5, 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

$u \in [H^1_0(\Omega)]^2$6

on a first-layer patch of refined elements (Shaikh et al., 2023, Bohne et al., 18 Feb 2026).

5. Specialized adaptive CR frameworks

For elliptic eigenvalue problems, the nonconformity creates a consistency error term in the exact identity

$u \in [H^1_0(\Omega)]^2$7

The paper (Hu et al., 2019) addresses this by combining high-accuracy gradient recovery with computable approximations of the consistency term. It proves asymptotic exactness for two eigenvalue estimators: $u \in [H^1_0(\Omega)]^2$8 and introduces a postprocessed eigenvalue

$u \in [H^1_0(\Omega)]^2$9

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 $p \in L^2_0(\Omega)$0, 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 $p \in L^2_0(\Omega)$1 because no upwinding or extra stabilization is used (Du et al., 2016).

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

$p \in L^2_0(\Omega)$2

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 (Shaikh et al., 2023).

For the first Dirichlet eigenpair of the $p \in L^2_0(\Omega)$3-Laplacian, the adaptive CR method works on $p \in L^2_0(\Omega)$4 polyhedral domains and proves that the estimator has vanishing limit,

$p \in L^2_0(\Omega)$5

that the discrete eigenvalues converge,

$p \in L^2_0(\Omega)$6

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 (Li et al., 4 Aug 2025).

For arbitrary odd degree $p \in L^2_0(\Omega)$7, high-order $p \in L^2_0(\Omega)$8 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 $p \in L^2_0(\Omega)$9 because for even $$\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,$$0 the degrees of freedom are not strictly local and require global constraints. Numerical experiments reported in (Bohne et al., 18 Feb 2026) illustrate optimal convergence rates for all considered values of $$\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,$$1.

6. Scope, limitations, and recurring distinctions

The classical instance-optimal theory in (Kreuzer et al., 2014) is proved in 2D only, for polygonal domains and right-hand sides $$\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,$$2 for Poisson and $$\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,$$3 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 (Kreuzer et al., 2014).

Several broader limitations recur across the literature. Elliptic eigenvalue asymptotic exactness in (Hu et al., 2019) assumes $$\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,$$4 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 (Hu et al., 2019, Du et al., 2016). The $$\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,$$5-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 (Li et al., 4 Aug 2025). The high-order $$\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,$$6 theory is restricted to 2D Poisson and odd polynomial degrees (Bohne et al., 18 Feb 2026).

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 [(Kreuzer et al., 2014); (Du et al., 2016); (Bohne et al., 18 Feb 2026)].