--- title: Divergence-Free Mixed VEM url: https://www.emergentmind.com/topics/divergence-free-mixed-virtual-element-method type: topic --- # Divergence-Free Mixed VEM 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 [1510.01655]. 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 [1703.00437, 1905.01579, 2108.09967, 2403.03886, 2004.11467, 2201.04417, 2601.18758, 2604.18406]. ## 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}u_h$ to vanish on each element as a polynomial of degree at most $k-1$, so the discrete velocity is pointwise divergence-free [1510.01655]. The same structural property underlies the 2D virtual element discretization of Navier–Stokes. There the global kernel \[ Z_h:=\{v_h\in V_h:\; b_h(v_h,q_h)=0\ \forall q_h\in Q_h\} \] satisfies $Z_h\subset Z$, where $Z$ is the continuous divergence-free space. Consequently, any $v_h\in Z_h$ 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 [1703.00437]. 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 $k-1$ distinct points on each edge, internal moments against $\mathcal G_{k-2}(K)^\perp$, and moments of $\operatorname{div}v$ against $P_{k-1}(K)\setminus\mathbb R$; the local pressure degrees of freedom are the moments against $P_{k-1}(K)$ [1510.01655]. The nonconforming formulation makes this structure more explicit. On each element $K$ one introduces \[ W_h^1(K)=\Bigl\{v\in H^1(K)^2:\;\operatorname{div}v\in P_{k-1}(K),\;\operatorname{rot}v=0,\; v\cdot n_K|_e\in P_k(e)\ \forall e\subset\partial K\Bigr\}, \] \[ \Phi_h(K)=\Bigl\{\phi\in H^2(K):\;\Delta^2\phi\in P_{k-3}(K),\; \phi|_e=0,\;(\Delta\phi)|_e\in P_{k-1}(e)\ \forall e\subset\partial K\Bigr\}, \] and the enhanced virtual space \[ \widetilde V_h(K)=W_h^1(K)\oplus\operatorname{curl}\Phi_h(K). \] A unisolvent set of local degrees of freedom is provided by edge normal moments, edge tangential moments, and cell moments. The actual nonconforming space $V_h(K)$ is selected from $\widetilde V_h(K)$ through orthogonality conditions involving the energy projector $\Pi_K^\nabla$, and one obtains the global broken space by imposing jump conditions on the moments across interior edges. In this construction $\operatorname{div}V_h(K)=P_{k-1}(K)$ on each element, hence $\operatorname{div}_hV_{h,0}=Q_h$ [2108.09967]. Computability is organized around local polynomial projectors. The most recurrent are the $H^1$-seminorm or energy projector $\Pi_k^{\nabla,K}$ and the $L^2$ projector $\Pi_\ell^{0,K}$. Local bilinear forms are then split into a polynomially consistent part and a stabilization on the kernel of the projector, \[ a_h^K(u_h,v_h)=a^K(\Pi_k^{\nabla,K}u_h,\Pi_k^{\nabla,K}v_h) +S^K((I-\Pi_k^{\nabla,K})u_h,(I-\Pi_k^{\nabla,K})v_h), \] with $S^K$ chosen symmetric positive definite and scaled like the continuous energy on $\ker\Pi_k^{\nabla,K}$ [1510.01655]. 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, \[ \sigma(x,\tau)=\mu(x)\bigl(\delta^\alpha+|\tau|^\alpha\bigr)^{(r-2)/\alpha}\tau, \] the local enhanced velocity space is chosen so that $\operatorname{div}v_h\in P_{k-1}(E)$ and the resulting discrete sequence \[ [P_{k+1}] \xrightarrow{\nabla} V_h \xrightarrow{\operatorname{div}} Q_h \to 0 \] 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 [2403.03886]. ## 3. Exact kernels, divergence-free bases, and discrete complexes Once $\operatorname{div}V_h=Q_h$ 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 \[ Z_{h,0}=\{v\in V_{h,0}:\;\operatorname{div}_h v=0\} \] is characterized by local moment conditions and then endowed with an explicit basis [2108.09967]. The basis construction for $Z_{h,0}$ 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 \[ \dim Z_{h,0} =\dim V_{h,0}-\dim Q_h =N_{V,i}+kN_{E,i}+(k-1)N_{E,i}+\frac{(k-1)(k-2)}2\,N_P. \] This yields a direct sum decomposition of the divergence-free subspace. For $k=1$ 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 [2108.09967]. 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, \[ \mathbb R \xrightarrow{i} W_h \xrightarrow{\nabla} \Sigma_h \xrightarrow{\operatorname{curl}} V_h \xrightarrow{\operatorname{div}} Q_h \to 0. \] Here $V_h$ is the 3D divergence-free mixed VEM velocity space, $Q_h$ 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 [1905.01579]. 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 \[ \widehat V_h^K := \{v\in[H^1(K)]^2:\; v|_{\partial K}\in[B_k(\partial K)]^2,\; \operatorname{div}v\in P_0(K),\; -\nu\Delta v-\nabla s\in\mathcal G_{k-2}(K)^\perp\}, \] \[ \widehat Q_h^K:=P_0(K), \] and the reduced global problem has a velocity solution coinciding with that of the full problem, while the reduced pressure is $\Pi_0p_h$ [1510.01655]. The nonconforming basis construction pushes this idea to an explicit kernel formulation. If $[\Phi]$ and $[\Psi]$ denote bases of $V_{h,0}$ and $Q_h$, the mixed matrix system is \[ \begin{pmatrix}A & B^T\\ B & 0\end{pmatrix} \begin{pmatrix}U\\ P\end{pmatrix} = \begin{pmatrix}F\\ 0\end{pmatrix}. \] If $Z$ is the incidence matrix whose columns are the coordinate vectors of the divergence-free basis in the $V_{h,0}$ basis, then $BZ=0$ and one may write $U=Z\alpha$, obtaining the reduced system \[ Z^T A Z\,\alpha = Z^T F, \] which is symmetric positive definite on the divergence-free subspace [2108.09967]. 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 $L^2$-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 $I_{RT_m}$ is built on Raviart–Thomas spaces over local subtriangulations of each polygon; for $m=k-1$ it preserves the discrete divergence exactly and leads to the modified right-hand side \[ F_h^{PR}(v_h):=\int_\Omega f\cdot I_{RT_{k-1}}(v_h)\,dx. \] The hydrostatic tests in that work show locking for classical right-hand-side discretizations as $\nu\to0$, whereas the pressure-robust variants remain at machine precision [2002.01830]. ## 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 $h$, and the method is stable and optimally convergent under the small-data condition $\gamma_h<1$ [1703.00437]. For steady non-Newtonian incompressible flow, the method is formulated in $W^{1,r}$ spaces with $10$ one even observes $O(h^k)$ velocity-pressure rates [2403.03886]. 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 $O(h^k)$ rates for the $H^1$-error of the velocity and the $L^2$-error of the pressure on Cartesian, tetrahedral, centroidal Voronoi, and random Voronoi meshes [1905.01579]. 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 $h^k$ energy-norm estimates for velocity and magnetic field together with $h^{k+1}$ $L^2$ estimates under dual regularity [2410.18376]. ## 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 $H(\operatorname{curl})$-conforming electric field and an $H(\operatorname{div})$-conforming magnetic flux. The VEM spaces reproduce the mixed curl-div structure on general polygonal meshes, and the discrete evolution preserves \[ \nabla\cdot B_h(\cdot,t)=0 \] exactly on every element for all $t$ provided the initial field is discrete divergence-free. Numerical tests on triangular, perturbed quadrilateral, and Voronoi meshes report $\|\operatorname{div}(B-B_h)\|_{L^2}=0$ at machine-zero level [2004.11467]. The 3D resistive MHD model uses a four-field VEM on general polyhedral meshes with velocity space $W_h\subset[H_0^1(\Omega)]^3$, edge space $E_h\subset H_0(\operatorname{curl};\Omega)$, face space $F_h\subset H_0(\operatorname{div};\Omega)$, and piecewise constant pressure. The associated commuting diagram ensures that no spurious divergence is created by the discrete $\nabla$ and $\nabla\times$ operators, and the computed fields satisfy $\|\nabla\cdot B_h\|_{L^2}$ and $\|\nabla\cdot u_h\|_{L^2}$ at the level of $10^{-13}$–$10^{-11}$ through refinements [2201.04417]. 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 \[ 0\to V_h \xrightarrow{\operatorname{rot}} E_h \xrightarrow{\operatorname{div}} P_h \to 0. \] 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 [2104.04096]. The divergence-free paradigm has also been coupled with high-order $C^1$ 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 $C^1$-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 [2601.18758]. 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 \[ u_h=\operatorname{curl}\Pi_h^1\phi_h+\sum_{j=1}^m c_{j,h}\nabla\Pi_h^1\varphi_{j,h}. \] Each term is exactly divergence-free at the discrete level, so no Lagrange multiplier or penalty is needed to enforce $\operatorname{div}u_h=0$. This places divergence-free VEM within a wider Hodge-decomposition program for structure-preserving discretization on polygonal meshes [2604.18406]. 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 [2207.01361]. Source: https://www.emergentmind.com/topics/divergence-free-mixed-virtual-element-method