Virtual Element Method Overview

• Virtual Element Method (VEM) is a Galerkin-type numerical method that discretizes PDEs on arbitrary polyhedral meshes without explicit shape functions.
• It employs projection operators to decompose displacement fields and compute strain energies accurately for elasticity problems.
• The method ensures optimal convergence and robustness on irregular meshes through effective stabilization and geometric integration techniques.

The Virtual Element Method (VEM) is a Galerkin-type numerical framework for the discretization of partial differential equations, distinguished by its capacity to operate on general polytopal (polygonal or polyhedral) meshes. Unlike classical finite element methods, VEM eliminates the need for explicit shape functions within each element by instead utilizing locally defined degrees of freedom (DOFs) coupled with projection and stabilization operations. This allows for the accurate and stable approximation of solutions to a wide class of boundary value problems, particularly in elasticity and related fields, on meshes with highly irregular or nonconvex elements.

1. Foundational Concepts: Space Decomposition and Projection Operators

VEM is grounded in the rigorous kinematic decomposition of element-level discrete function spaces. For three-dimensional elasticity on a polyhedral element $E$, any displacement field $v$ is decomposed as

$v = T_R v + T_C v + (v - T_P v),$

where:

• $T_R v$ projects onto the rigid body motion space $R$ (i.e., all translations and infinitesimal rotations; see equations (28) and (30)).
• $T_C v$ projects onto the space $C$ of constant (uniform) strain, i.e., fields with constant symmetric gradient (see equation (32)).
• $T_P = T_R + T_C$ projects onto the full space $P$ of linear displacements, and $(v - T_P v)$ is the nonpolynomial higher-order component.

The variational formulation is constructed such that the energy contribution from the polynomial (constant strain and rigid body) subspaces is computed exactly by geometric projection. The "virtual" or "nonpolynomial" component is controlled by a stabilization term defined exclusively via the DOFs. These projection maps rely on the geometry (volume, centroid, face normals, etc.) and are implemented through explicit formulas—enabling precise control of consistency, patch test satisfaction, and energy decomposition.

2. Mathematical Structure and Stiffness Assembly

The elementwise bilinear form for elasticity adopts the structure: $$a_h(u, v) = a_F(T_C u, T_C v) + s_E(u - T_P u, v - T_P v) \quad \text{(see equation (41))},$$ where:

• $a_F(\cdot, \cdot)$ is the continuous elasticity bilinear form (exactly integrated on constant strain projections).
• $s_E$ is the stabilization term acting exclusively on the kernel of the projection (i.e., on nonpolynomial modes).

The assembled element stiffness matrix is given by: $$K_E = |E| (W_C D W_C^T) + (I - P_P) S_E (I - P_P) \quad \text{(see equation (70))},$$ where $D$ is a 6×6 constitutive-geometry matrix, $W_C$ encodes the canonical basis for constant-strain fields, and $P_P$ represents the nodal matrix instantiating $T_P$. All these matrices are assembled from mesh geometry only, rendering the scheme independent of interior explicit basis construction and replacing quadrature with geometric summations (e.g., via divergence theorem reductions, as in equations (31), (35)).

3. Algorithmic Workflow and Implementation Details

A standard first-order VEM implementation for 3D elasticity proceeds as follows:

1. Mesh Generation: Partition the computational domain $\Omega$ into arbitrary polyhedral elements; extract geometry (vertex positions, face connectivity, centroids, normals).
2. Local Space Definition: Set up nodal DOFs at vertices (typically, physical displacement values); enforce $C^0$ continuity through face-level barycentric or harmonic coordinates.
3. Local Projections: Compute projected averages and strains ($\bar{v}$, $\bar{\varepsilon}$) via boundary quadrature; assemble $T_R$, $T_C$ by explicit geometric formulas.
4. Stiffness Construction: Evaluate exact energy for rigid body and linear strain via geometric quantities; define and apply a stabilization matrix $s_E$ to the nonpolynomial (hourglass) modes; assemble as in (41)/(70).
5. Global Assembly and Solve: Build the global system by finite element-type assembly; enforce physical boundary conditions; solve the linear system.

Reliable implementation on arbitrary polyhedra mandates robust data structures for face and vertex traversal, accurate computation of centroids and normals (especially on nonconvex faces), and careful management of stabilization contributions to ensure spectral stability and eliminate spurious kernel modes.

4. Numerical Performance and Convergence Properties

VEM yields strong numerical performance, validated by rigorous patch and convergence tests:

• Patch Test: On a unit cube, for exactly linear displacement fields, VEM delivers displacement and stress errors near machine precision (order $10^{-16}$) on centroidal Voronoi, random Voronoi, and nonconvex polyhedral meshes.
• Beam and Shear Testing: In cantilever beam benchmarks under shear loading, the scheme attains second-order convergence in displacement and first-order convergence in stress for various mesh types (CVT, random, and uniform hexahedral meshes).
• Mesh Robustness: Optimal convergence rates (matching those of classical FEM) are observed even on highly irregular polyhedral meshes, with negligible accuracy loss compared to uniform hexahedral discretizations.
• Stabilization Scaling: A parameter sweep over stabilization scaling coefficients reveals an optimal regime (typically $y \approx 1$ in stabilization scaling) that remains almost invariant with mesh refinement, indicating robust parameter insensitivity.

Errors plotted against total DOF counts show marginal accuracy advantages for regular hexahedral meshes, but polyhedral VEM maintains competitiveness and is robust against extreme mesh irregularities.

5. Comparison with Classical and Mimetic Discretizations

VEM distinguishes itself by:

• Eliminating Explicit Shape Functions: All computation is boundary- and DOF-based; there is no need for volumetric shape functions or quadrature over nontrivial domains.
• Polyhedral Mesh Compatibility: Arbitrary polyhedral cells are handled directly, relaxing mesh generator constraints and increasing adaptability for complex geometries and evolving domains.
• Guaranteed Patch Test and Consistency: The method exactly reproduces all rigid body and constant strain modes, ensuring it passes patch tests and exhibits both consistency and optimal convergence.
• Geometric Integration Efficiency: All required volume or strain projections reduce to boundary integrals over triangles or polygons, bypassing computationally intensive interior quadrature used in FEM on general polytopes.

Compared to mimetic finite difference (MFD) methods, VEM can be viewed as a conforming Galerkin variant with a virtual (implicit) approximation space, inheriting much of the structural flexibility of MFD while operating in the mathematical framework of projection-based finite elements.

Potential limitations include the burden of designing optimal stabilization terms $s_E$—especially in the presence of highly anisotropic or degenerate cells—and the necessity for careful geometric handling in mesh data structures. However, these are compensated by the method’s consistency, stability, and generality.

6. Theoretical Guarantees and Practical Implications

By construction, VEM's energy decomposition and projection strategy guarantee:

• Exactness on Linear Deformations: All modes in the linear (affine) subspace are captured without error, with strain energies computed exactly from geometric quantities.
• Optimal Order Convergence: Theoretical analysis and experiments show that VEM achieves the optimal rate (first-order in energy, up to second-order in $L^2$) even for nonconvex, star-shaped polyhedral meshes.
• Computational Efficiency: No interior shape functions or high-order quadrature are required; all matrix assembly is boundary-based or geometric, enabling scalable and efficient implementations.
• Consistency and Stability for Complex and Evolving Meshes: Robustness against mesh distortion and cell degeneracy, with stabilization terms ensuring control over spurious and hourglass modes.

These properties, verified by rigorous analysis and numerical experiments, position VEM as an attractive discretization strategy for three-dimensional elasticity, especially where standard finite element methods face meshing or integration barriers.

7. Significance and Research Directions

The formulation and validation of VEM for 3D elasticity on arbitrary polyhedral meshes (Gain et al., 2013) established key architectural features—kinematic decomposition, geometric projection-based assembly, and stabilization—that have subsequently underpinned theoretical and algorithmic advances in VEM for a widening array of PDEs and numerical regimes. This includes extensions to nonlinear elasticity, plasticity, quasilinear problems, coupled multiphysics, and eigenvalue computations. The geometric and algebraic flexibility intrinsic to VEM continues to motivate research in mesh adaptivity, high-order accuracy, and integration with advanced meshing and data-driven methodologies in computational mechanics and engineering analysis.