Discrete Stokes Polygonal Complex

• Discrete Stokes polygonal complex is a compatible discretization that reproduces the continuous Stokes sequence on polygonal meshes using exact discrete curl and divergence operators.
• It leverages methods like the Virtual Element Method, Morley-type, and Scott–Vogelius-type constructions to achieve robust divergence-free velocity approximations and optimal convergence.
• The complex underpins diverse frameworks for pressure robustness and divergence-preserving reconstructions, enhancing the efficiency and accuracy of fluid dynamics simulations in complex geometries.

Discrete Stokes polygonal complex most commonly denotes a compatible discretization of the two-dimensional Stokes sequence on polygonal meshes in which a scalar stream-function space, a velocity space, and a pressure space are linked by discrete curl and divergence operators in an exact sequence. The clearest polygonal realization in the cited literature is the Virtual Element Method sequence

$0 \xrightarrow[]{\,i\,} \Phi_h \xrightarrow[]{\,\curl\,} V_h \xrightarrow[]{\,\div\,} Q_h \xrightarrow[]{\,0\,} 0,$

which reproduces the continuous Stokes complex on simply connected polygonal domains and yields exactly divergence-free discrete velocities (Veiga et al., 2018). Closely related constructions appear on polygonal or polyhedral cell complexes through generalized finite element systems, Morley-type and Scott–Vogelius-type virtual spaces, weak Galerkin bases, compatible discrete operator schemes, and quadrilateral or rectangular nonconforming complexes (Christiansen et al., 2016).

1. Continuous sequence and the meaning of exactness

On a simply connected polygonal domain $\Omega\subset \mathbb R^2$, the continuous Stokes complex used throughout this literature is

$0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$

With the convention

$\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$

every divergence-free $H^1$-velocity is the curl of a scalar stream function, and every zero-mean pressure is the divergence of an $H_0^1$-field. Exactness therefore identifies the divergence-free kernel with the image of the scalar potential space rather than only with a weak constraint (Veiga et al., 2018).

This exactness has several concrete consequences in the cited works. It yields pointwise or exact discrete incompressibility, supports equivalent stream-function/curl formulations of Stokes or Navier–Stokes, and separates velocity approximation from spurious pressure contamination. In the Morley-type stream-function formulation for Navier–Stokes, the same continuous identities are used in the form

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$

so that the Stokes sequence becomes the structural device for recovering velocity, vorticity, and pressure from a scalar stream approximation (Adak et al., 2022).

2. Exact polygonal complexes in the Virtual Element Method

The most explicit polygonal discrete Stokes complex in the supplied literature is the two-dimensional Virtual Element construction of Beirão da Veiga, Russo, and Vacca. On general polygonal meshes satisfying the standard assumptions that each element $E$ is star-shaped with respect to a ball of radius comparable to $h_E$ and that the distance between any two vertices is bounded below by a multiple of $h_E$, the pressure space is

$\Omega\subset \mathbb R^2$0

The velocity space $\Omega\subset \mathbb R^2$1 is an $\Omega\subset \mathbb R^2$2-conforming vector VEM space with polynomial degree $\Omega\subset \mathbb R^2$3, and the new scalar stream-function space $\Omega\subset \mathbb R^2$4 is a global $\Omega\subset \mathbb R^2$5-conforming, $\Omega\subset \mathbb R^2$6-virtual space on polygons. The local and global exact sequences are

$\Omega\subset \mathbb R^2$7

and

$\Omega\subset \mathbb R^2$8

The key identities are

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

so every discrete divergence-free velocity in $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$0 is exactly the curl of a discrete stream function in $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$1, while

$0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$2

The same paper also gives a reduced exact complex with piecewise constant pressure (Veiga et al., 2018).

The local velocity space is defined through edge traces, divergence moments, and moments against $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$3, while the local stream-function space is defined through polynomial traces, polynomial gradient traces, and a biharmonic residual in $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$4. The degrees of freedom are chosen so that the canonical projections

$0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$5

are computable for $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$6, and the stream-function degrees of freedom determine the velocity degrees of freedom of $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$7. This computability is essential because virtual functions are not known in closed form inside each polygon.

The exact sequence is also used to rewrite the discrete Navier–Stokes problem on the smaller scalar space. Since

$0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$8

the mixed formulation and the curl formulation produce the same discrete velocity $0 \xrightarrow[]{\,i\,} H_0^2(\Omega) \xrightarrow[]{\,\curl\,} [H_0^1(\Omega)]^2 \xrightarrow[]{\,\div\,} L_0^2(\Omega) \xrightarrow[]{\,0\,} 0.$9. Numerically, the curl formulation uses $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$0 fewer global unknowns than the mixed formulation, but its condition number behaves like $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$1, versus $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$2 for the velocity-pressure formulation. The same study reports that the predicted convergence rates are confirmed and that the three trilinear form choices behave almost identically on the tested problems (Veiga et al., 2018).

3. Polygonal VEM variants: divergence-free, Morley-type, and Scott–Vogelius-type constructions

A second major polygonal line is the divergence-free VEM for the Stokes problem on polygonal meshes. There the local velocity space is again a Stokes-type virtual space, but the decisive structural property is stated directly as

$\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$3

Because $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$4 on each element and the discrete incompressibility equation is imposed against all of $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$5, the final discrete velocity is pointwise divergence-free. The same paper proves that the full mixed problem is immediately equivalent to a reduced problem in which the velocity is unchanged and the reduced pressure is the elementwise constant projection of the full pressure; the saving is $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$6 degrees of freedom (Veiga et al., 2015).

The Morley-type virtual element method for Navier–Stokes uses a different polygonal pairing. Its stream-function space $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$7 is a nonconforming $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$8-type polygonal analogue of the Morley element, its companion vector space $\curl \psi := \left(\frac{\partial \psi}{\partial y}, -\frac{\partial \psi}{\partial x}\right),$9 is a Crouzeix–Raviart-type virtual element space, and the discrete Stokes-complex identity actually used is

$H^1$0

where $H^1$1 is the subspace of $H^1$2 with vanishing discrete divergence pairing against the piecewise constant pressure space $H^1$3. In this setting the complex is the mechanism for pressure recovery: once $H^1$4 is known, the recovered pressure $H^1$5 is obtained from an auxiliary mixed problem on $H^1$6, and the paper proves optimal-order estimates for recovered velocity, vorticity, and pressure on general polygonal meshes, with domains not necessarily convex (Adak et al., 2022).

A different polygonal extension is the VEM generalization of the Scott–Vogelius method. Here the velocity space is obtained by applying the standard scalar conforming VEM space componentwise and pairing it with discontinuous elementwise polynomial pressures,

$H^1$7

This construction is explicitly presented as a polygonal generalization of Scott–Vogelius, but not as a full discrete Stokes complex. Its compatibility is expressed instead by the exact identity

$H^1$8

and by the projected divergence-free property

$H^1$9

For meshes satisfying the standard VEM assumptions, the paper proves a discrete inf-sup condition for $H_0^1$0, while the lowest-order case on triangular and square meshes reproduces the classical unstable situations (Manzini et al., 2021).

4. Generalized and compatible discrete frameworks on cell complexes

Beyond VEM, the supplied literature contains two broader structural frameworks. The first is the generalized finite element systems framework for smooth differential forms. It extends finite element systems by allowing pullback, full trace, and double-trace restrictions, and it is formulated on a cellular complex. In dimension two, it produces low-order and high-order exact complexes with regularity $H_0^1$1 and $H_0^1$2. For example, on a Clough–Tocher split $H_0^1$3 of a triangle, one exact sequence is

$H_0^1$4

and another is

$H_0^1$5

The last two spaces in these complexes provide conforming Stokes pairs with continuous or discontinuous pressure. The framework also guarantees commuting interpolators and identifies some of the resulting spaces as minimal. The paper is explicit that the abstract framework applies to a cellular complex, but its concrete local constructions are composite polynomial spaces on refined simplices rather than arbitrary single polygonal cells (Christiansen et al., 2016).

The second framework is the compatible discrete operator approach on three-dimensional polyhedral meshes. There the discrete unknowns are organized as cochains on primal and dual complexes, the topological laws are discretized exactly by incidence matrices, and the constitutive relations are approximated by discrete Hodge operators. The discrete gradient, curl, and divergence operators satisfy the cochain identities

$H_0^1$6

and the schemes preserve local mass and momentum conservation. Two families are analyzed: one with pressure degrees of freedom at mesh vertices and one with pressure degrees of freedom at cells. Discrete stability is proved by new discrete Poincaré inequalities, and first-order error estimates are derived by commutator arguments for smooth solutions. This is a polyhedral, curl-based realization of the same compatible philosophy, but expressed in mimetic and cochain language rather than in virtual element exact sequences (Bonelle et al., 2014).

A third compatible route is the weak Galerkin divergence-free method on polygons and polyhedra. For the lowest-order weak Galerkin element,

$H_0^1$7

the paper defines the discrete divergence-free space

$H_0^1$8

and constructs explicit bases of $H_0^1$9 on general polygonal and polyhedral meshes. In two dimensions the basis decomposes into element-interior functions, edge-tangential functions, and vertex-hull functions; in three dimensions it decomposes into element-interior functions, face-tangential functions, and edge-solid-hull functions. The resulting reduced system on $$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$0 is symmetric and positive definite (Mu et al., 2016).

5. Quadrilateral and rectangular discrete Stokes complexes

A separate line of work develops discrete Stokes complexes on quadrilateral and rectangular meshes using nonconforming finite elements made of piecewise polynomials directly on physical cells. On convex quadrilateral grids, Zhang constructs a stable nonconforming Stokes pair and proves the exact sequence

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$1

Here $$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$2 is a quadrilateral Morley space,

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$3

the scalar velocity space underlying $$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$4 is

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$5

and the pressure space is discontinuous $$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$6 on each quadrilateral. The paper also proves the commuting relations

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$7

and extends the construction to mixed meshes consisting of triangles and convex quadrilaterals. In that precise sense, the quadrilateral complex becomes a polygonal-type complex on meshes built from two polygonal cell types (Zhang, 2013).

On uniform rectangular meshes, a different nonconforming complex is

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$8

with a commuting diagram

$$u=\operatorname{curl}\psi,\qquad \omega=\operatorname{rot}u=-\Delta\psi,$$9

The scalar space $E$0 is a 12-DoF rectangular non-$E$1 plate element, $E$2 is a 12-DoF vector space with weak continuity of normal moments and tangential averages, and $E$3 is piecewise constant pressure. Although the method is nonconforming and only rigorous for uniform rectangular partitions, it provides an exact commuting discrete Stokes complex and unexpectedly high accuracy: $E$4 velocity error in the discrete $E$5-norm, $E$6 velocity error in $E$7 on convex domains, and $E$8 postprocessed pressure accuracy (Zhou et al., 2018).

6. Divergence-preserving reconstruction, pressure robustness, and three-dimensional directions

Several polygonal methods in the supplied literature show that an exact or projected divergence-free pair does not by itself settle the full structure of a discrete Stokes complex; the right-hand side and reconstruction operators must also preserve the relevant orthogonality. The clearest statement appears in the paper on “really pressure-robust” VEM. Starting from a divergence-free polygonal VEM, it argues that standard $E$9-best approximation of virtual test functions destroys the orthogonality to gradient forces, and it introduces a local Raviart–Thomas reconstruction $h_E$0 on a polygonal subtriangulation with the commuting identity

$h_E$1

For $h_E$2, divergence is preserved exactly. The reconstruction also satisfies

$h_E$3

and it restores the pressure-robust annihilation of gradient forces. This provides a computable $h_E$4-conforming transfer from polygonal virtual spaces to Raviart–Thomas spaces on local simplicial refinements, which is a distinctly complex-like ingredient even though no full exact polygonal sequence is written down (Frerichs et al., 2020).

A related but different polygonal method is the lowest-order staggered discontinuous Galerkin discretization on general convex polygonal meshes. Its unknowns live on a primal polygonal partition plus a simplicial star-subdivision, and the right-hand side is modified by a divergence-preserving reconstruction $h_E$5 into a polygonal $h_E$6 space. The key identity

$h_E$7

is the mechanism behind pressure robustness. The reconstructed discrete velocity $h_E$8 is exactly divergence free, and the method achieves optimal first-order convergence together with a superconvergent estimate

$h_E$9

under stronger regularity (Zhao et al., 2020).

Three-dimensional extensions remain materially harder. The Freudenthal-mesh paper states that conforming exactly incompressible discretizations are now well understood in two dimensions but remain poorly understood in three dimensions. It formulates two conjectures for the Scott–Vogelius pair on uniform meshes: inf-sup stability for velocity degree $h_E$0, while the best result available in the literature is for $h_E$1, and the existence of a stable space decomposition of the kernel of the divergence for $h_E$2. It also presents numerical evidence supporting these conjectures (Farrell et al., 2022). This suggests that a fully satisfactory three-dimensional polygonal or polyhedral analogue of the two-dimensional exact Stokes polygonal complex remains an open direction rather than a closed theory.