Divergence-Free Mixed VEM

• Divergence-free mixed VEM are mixed virtual element discretizations that ensure exact, pointwise incompressibility on polygonal and polyhedral meshes.
• They utilize local polynomial projectors, tailored stabilization, and exact sequence structures to obtain pressure-robust and reduced formulations.
• These methods extend to nonlinear, three-dimensional, and coupled problems like magnetohydrodynamics and quad-curl formulations while preserving key solenoidal properties.

Divergence-free mixed virtual element methods are mixed virtual element discretizations for incompressible or solenoidal partial differential equations in which the discrete field is constrained by construction to lie in an exactly divergence-free kernel, often pointwise on each element. In the Stokes setting, the foundational idea is to choose local virtual velocity spaces so that $\operatorname{div} V_h^K \subset P_{k-1}(K)$ and, in the classical conforming construction, $\operatorname{div} V_h^K=P_{k-1}(K)$, which makes the discrete continuity equation exact and yields pointwise divergence-free discrete velocities on polygonal meshes (Veiga et al., 2015). The same principle has since been extended to polygonal and polyhedral discretizations of Navier–Stokes, nonconforming Stokes, non-Newtonian flow, magnetohydrodynamics, phase-field models, and quad-curl problems, usually with computable polynomial projectors, stabilization on non-polynomial components, and exact-sequence or Hodge-theoretic structure (Veiga et al., 2017, Veiga et al., 2019, Kwak et al., 2021, Antonietti et al., 2024, Alvarez et al., 2020, Veiga et al., 2022, Silgado et al., 26 Jan 2026, Brenner et al., 20 Apr 2026).

1. Foundational formulation and defining property

The canonical prototype is the Stokes problem on a polygonal domain $\Omega\subset\mathbb R^2$, written in mixed form on $V=[H_0^1(\Omega)]^2$ and $Q=L_0^2(\Omega)$. The divergence-free mixed VEM replaces these by discrete spaces

$$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$

with local velocity space

$$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$

Because $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$, the discrete divergence bilinear form is exact on the chosen pressure space; in the original divergence-free Stokes construction, $b(u_h,q_h)=0$ for all $q_h\in Q_h$ forces $\operatorname{div} V_h^K=P_{k-1}(K)$0 to vanish on each element as a polynomial of degree at most $\operatorname{div} V_h^K=P_{k-1}(K)$1, so the discrete velocity is pointwise divergence-free (Veiga et al., 2015).

The same structural property underlies the 2D virtual element discretization of Navier–Stokes. There the global kernel

$\operatorname{div} V_h^K=P_{k-1}(K)$2

satisfies $\operatorname{div} V_h^K=P_{k-1}(K)$3, where $\operatorname{div} V_h^K=P_{k-1}(K)$4 is the continuous divergence-free space. Consequently, any $\operatorname{div} V_h^K=P_{k-1}(K)$5 is pointwise divergence-free on each polygonal element, and the convective term can be discretized on a solenoidal discrete velocity field rather than on a merely weakly incompressible one (Veiga et al., 2017).

This exactness is the central distinguishing feature of the divergence-free mixed VEM. It separates the framework from methods in which incompressibility is imposed only after projection or only in a cell-average sense.

2. Local virtual spaces, degrees of freedom, projectors, and stabilization

The virtual element philosophy is to define approximation spaces through trace conditions, internal differential constraints, and moment conditions, while retaining only degrees of freedom that make the necessary polynomial projections computable. In the conforming Stokes family, the local velocity degrees of freedom are vertex values, values at $\operatorname{div} V_h^K=P_{k-1}(K)$6 distinct points on each edge, internal moments against $\operatorname{div} V_h^K=P_{k-1}(K)$7, and moments of $\operatorname{div} V_h^K=P_{k-1}(K)$8 against $\operatorname{div} V_h^K=P_{k-1}(K)$9; the local pressure degrees of freedom are the moments against $\Omega\subset\mathbb R^2$0 (Veiga et al., 2015).

The nonconforming formulation makes this structure more explicit. On each element $\Omega\subset\mathbb R^2$1 one introduces

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

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

and the enhanced virtual space

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

A unisolvent set of local degrees of freedom is provided by edge normal moments, edge tangential moments, and cell moments. The actual nonconforming space $\Omega\subset\mathbb R^2$5 is selected from $\Omega\subset\mathbb R^2$6 through orthogonality conditions involving the energy projector $\Omega\subset\mathbb R^2$7, and one obtains the global broken space by imposing jump conditions on the moments across interior edges. In this construction $\Omega\subset\mathbb R^2$8 on each element, hence $\Omega\subset\mathbb R^2$9 (Kwak et al., 2021).

Computability is organized around local polynomial projectors. The most recurrent are the $V=[H_0^1(\Omega)]^2$0-seminorm or energy projector $V=[H_0^1(\Omega)]^2$1 and the $V=[H_0^1(\Omega)]^2$2 projector $V=[H_0^1(\Omega)]^2$3. Local bilinear forms are then split into a polynomially consistent part and a stabilization on the kernel of the projector,

$V=[H_0^1(\Omega)]^2$4

with $V=[H_0^1(\Omega)]^2$5 chosen symmetric positive definite and scaled like the continuous energy on $V=[H_0^1(\Omega)]^2$6 (Veiga et al., 2015).

In nonlinear settings the same pattern persists, but the stabilization is adapted to the constitutive law. For steady non-Newtonian incompressible flow with Carreau–Yasuda stress,

$V=[H_0^1(\Omega)]^2$7

the local enhanced velocity space is chosen so that $V=[H_0^1(\Omega)]^2$8 and the resulting discrete sequence

$V=[H_0^1(\Omega)]^2$9

is exact. The discrete nonlinear operator uses projected symmetric gradients together with a “dofi-dofi” stabilization tailored to mimic the continuous operator’s monotonicity and boundedness (Antonietti et al., 2024).

3. Exact kernels, divergence-free bases, and discrete complexes

Once $Q=L_0^2(\Omega)$0 is available, the discrete kernel is not merely algebraic but geometric: it is a subspace of genuinely solenoidal functions. This viewpoint is particularly sharp in the nonconforming Stokes construction, where

$Q=L_0^2(\Omega)$1

is characterized by local moment conditions and then endowed with an explicit basis (Kwak et al., 2021).

The basis construction for $Q=L_0^2(\Omega)$2 proceeds in four mutually linearly independent families: interior-vertex modes, edge-tangential modes, edge-normal enriched modes, and cell-interior modes. The total number of these functions matches

$Q=L_0^2(\Omega)$3

This yields a direct sum decomposition of the divergence-free subspace. For $Q=L_0^2(\Omega)$4 on a triangular mesh, only the vertex and edge-tangential families remain, and the construction recovers exactly the classical Crouzeix–Raviart divergence-free basis of Brenner–Thomasset (Kwak et al., 2021).

The basis perspective is one manifestation of a broader exact-sequence structure. In three dimensions, the Stokes complex for virtual elements establishes an exact discrete complex on general polyhedral partitions,

$Q=L_0^2(\Omega)$5

Here $Q=L_0^2(\Omega)$6 is the 3D divergence-free mixed VEM velocity space, $Q=L_0^2(\Omega)$7 is the discontinuous polynomial pressure space, and the intermediate spaces are designed so that the discrete gradient, curl, and divergence operators reproduce the topological structure of the continuous Stokes complex (Veiga et al., 2019).

This exact-complex viewpoint is not an ancillary algebraic refinement. It is the mechanism that explains why divergence constraints remain exact on general polytopal meshes, why compatible magnetic or vorticity-like fields can be coupled without spurious modes, and why elimination strategies can be formulated on the kernel itself.

4. Pressure elimination, reduced formulations, and pressure-robustness

A major practical consequence of exact divergence-free kernels is that the mixed saddle-point problem can often be reduced without altering the discrete velocity. In the original divergence-free VEM for Stokes, most internal velocity moments of divergence type and all pressure degrees of freedom except one constant per element can be condensed out locally. The reduced local spaces are

$Q=L_0^2(\Omega)$8

$Q=L_0^2(\Omega)$9

and the reduced global problem has a velocity solution coinciding with that of the full problem, while the reduced pressure is $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$0 (Veiga et al., 2015).

The nonconforming basis construction pushes this idea to an explicit kernel formulation. If $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$1 and $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$2 denote bases of $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$3 and $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$4, the mixed matrix system is

$$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$5

If $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$6 is the incidence matrix whose columns are the coordinate vectors of the divergence-free basis in the $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$7 basis, then $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$8 and one may write $$V_h:=\{v\in[H_0^1(\Omega)]^2:\; v|_K\in V_h^K\ \forall K\},\qquad Q_h:=\{q\in L_0^2(\Omega):\; q|_K\in P_{k-1}(K)\ \forall K\},$$9, obtaining the reduced system

$$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$0

which is symmetric positive definite on the divergence-free subspace (Kwak et al., 2021).

A recurrent misconception is that exact divergence-free discretization automatically implies pressure-robustness. That claim is explicitly rejected in the analysis of pressure-robust VEM for Stokes: divergence-free VEM on polygonal meshes is not really pressure-robust as long as the right-hand side is not discretized in a careful manner. The standard load approximation based on an $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$1-best approximation of the virtual test function destroys the divergence and therefore destroys the orthogonality between divergence-free test functions and gradient forces. To repair this, a divergence-preserving reconstruction $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$2 is built on Raviart–Thomas spaces over local subtriangulations of each polygon; for $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$3 it preserves the discrete divergence exactly and leads to the modified right-hand side

$$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$4

The hydrostatic tests in that work show locking for classical right-hand-side discretizations as $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$5, whereas the pressure-robust variants remain at machine precision (Frerichs et al., 2020).

5. Nonlinear flow models, three dimensions, and conservative extensions

The divergence-free mixed VEM framework extends beyond linear Stokes flow without abandoning exact incompressibility. For the 2D Navier–Stokes equations on polygonal meshes, the viscous form is discretized by the same projector-plus-stabilization mechanism, while the convective term is assembled in both straight and skew-symmetric versions. The discrete kernel remains pointwise divergence-free, the discrete inf–sup condition is uniform in $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$6, and the method is stable and optimally convergent under the small-data condition $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$7 (Veiga et al., 2017).

For steady non-Newtonian incompressible flow, the method is formulated in $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$8 spaces with $$V_h^K := \left\{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_{k-1}(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\ \text{for some }s\in L^2(K)\right\}.$$9, using an enhanced velocity space, exact divergence-free enforcement, and a nonlinear stabilization designed to reproduce strong monotonicity and Hölder continuity. Under the regularity assumptions stated in the analysis, the method yields

$\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$0

$\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$1

so that for full regularity one recovers the classical rates $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$2 and $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$3, while for $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$4 one even observes $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$5 velocity-pressure rates (Antonietti et al., 2024).

In three dimensions, the virtual Stokes complex shows that divergence-free mixed VEM is not restricted to planar formulations. On polyhedral meshes the 3D spaces support both Stokes and diffusion-dominated Navier–Stokes approximations, with exact divergence-free discrete velocities, uniform inf–sup stability, and optimal $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$6 rates for the $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$7-error of the velocity and the $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$8-error of the pressure on Cartesian, tetrahedral, centroidal Voronoi, and random Voronoi meshes (Veiga et al., 2019).

A closely related conservative nonconforming formulation has also been developed for stationary incompressible magnetohydrodynamics. There the satisfactory divergence-free property of the virtual velocity field is used to ensure mass conservation, and the resulting method attains optimal $\operatorname{div}V_h^K\subset P_{k-1}(K)=Q_h^K$9 energy-norm estimates for velocity and magnetic field together with $b(u_h,q_h)=0$0 $b(u_h,q_h)=0$1 estimates under dual regularity (Dong et al., 2024).

6. Magnetohydrodynamics, phase-field coupling, quad-curl formulations, and solvers

The same design philosophy has been transported from incompressible flow to coupled field theories with solenoidal constraints. In the 2D resistive MHD electromagnetics subsystem, the discrete unknowns are an $b(u_h,q_h)=0$2-conforming electric field and an $b(u_h,q_h)=0$3-conforming magnetic flux. The VEM spaces reproduce the mixed curl-div structure on general polygonal meshes, and the discrete evolution preserves

$b(u_h,q_h)=0$4

exactly on every element for all $b(u_h,q_h)=0$5 provided the initial field is discrete divergence-free. Numerical tests on triangular, perturbed quadrilateral, and Voronoi meshes report $b(u_h,q_h)=0$6 at machine-zero level (Alvarez et al., 2020).

The 3D resistive MHD model uses a four-field VEM on general polyhedral meshes with velocity space $b(u_h,q_h)=0$7, edge space $b(u_h,q_h)=0$8, face space $b(u_h,q_h)=0$9, and piecewise constant pressure. The associated commuting diagram ensures that no spurious divergence is created by the discrete $q_h\in Q_h$0 and $q_h\in Q_h$1 operators, and the computed fields satisfy $q_h\in Q_h$2 and $q_h\in Q_h$3 at the level of $q_h\in Q_h$4–$q_h\in Q_h$5 through refinements (Veiga et al., 2022).

A broader MHD formulation on polygonal meshes develops two discrete chains, one for electromagnetics and one for fluid flow, and emphasizes the exact discrete de Rham sequence

$q_h\in Q_h$6

Within that structure, a Newton–Krylov linearization produces well-posed saddle-point systems while preserving the divergence-free condition on the magnetic field at each nonlinear iteration (Naranjo-Alvarez et al., 2021).

The divergence-free paradigm has also been coupled with high-order $q_h\in Q_h$7 virtual elements for diffuse-interface flow. In the semi- and fully discrete virtual element methods for the Navier–Stokes–Cahn–Hilliard system, the spatial discretization combines divergence-free velocity spaces with $q_h\in Q_h$8-conforming phase-field spaces, and the skew-symmetric treatment of convection in the Cahn–Hilliard equation yields exact mass conservation together with discrete energy bounds (Silgado et al., 26 Jan 2026).

An allied but structurally revealing example is the quad-curl problem on planar domains. There the discrete field is not approximated directly in a mixed velocity-pressure pair; instead, it is reconstructed as

$q_h\in Q_h$9

Each term is exactly divergence-free at the discrete level, so no Lagrange multiplier or penalty is needed to enforce $\operatorname{div} V_h^K=P_{k-1}(K)$00. This places divergence-free VEM within a wider Hodge-decomposition program for structure-preserving discretization on polygonal meshes (Brenner et al., 20 Apr 2026).

On the linear algebra side, BDDC preconditioners have been extended to the saddle-point systems arising from divergence-free virtual element discretizations of the 2D Stokes equations. Under suitable hypotheses on the choice of primal unknowns, the preconditioned linear system is symmetric positive definite, so the preconditioned conjugate gradient method can be used. The accompanying theory estimates the condition number of the preconditioned system, and numerical experiments show scalability, quasi-optimality, and robustness with respect to polygon shape; a slightly larger coarse space can also accelerate convergence (Bevilacqua et al., 2022).