Neural Approximated Virtual Element Method

• NAVEM is a hybrid numerical method that replaces traditional virtual basis functions with neural network approximations on polygonal elements.
• The method simplifies elasticity formulations by eliminating projection and stabilization operators while achieving comparable convergence rates with lower error constants.
• Continuous extensions like B-NAVEM and P-NAVEM ensure exact inter-element continuity and robust performance in large-deformation and nonlinear applications.

Searching arXiv for NAVEM and closely related papers to ground the article in the cited literature. Neural Approximated Virtual Element Method (NAVEM) denotes a family of polygonal discretization methods in which virtual basis functions are replaced by neural-network approximations inside an explicit ansatz space. In the elasticity setting, NAVEM is presented as a hybrid technique combining classical concepts from the Finite Element Method and the Virtual Element Method with recent advances in deep neural networks; specifically, it is a polygonal method in which the virtual basis functions are element-wise approximated by a neural network, eliminating the need for stabilization or projection operators typical of the standard virtual element method (Berrone et al., 8 Jul 2025). The approach was first introduced for lowest-order elliptic problems, then extended to general polygonal elements and hanging-node configurations, and later complemented by globally continuous variants termed B-NAVEM and P-NAVEM (Berrone et al., 2023, Berrone et al., 2024, Berrone et al., 14 Jan 2026).

1. Standard VEM background and research lineage

In the elasticity formulation, NAVEM is best understood as a modification of the lowest-order Virtual Element Method on a polygon $E$ with vertices $\{v_j\}_{j=1}^{N_V}$. The standard displacement space is

$$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$

with degrees of freedom given by the nodal values $v(v_j)$. Standard VEM introduces the projection

$\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$

through gradient orthogonality and a boundary-average constraint, together with an $L^2$ projector $\Pi^{0,E}_0$ into constants. Its local discrete bilinear form is

$$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$

where $S_E$ is a user-chosen stabilization and $\alpha_E(w_h)$ is a material- and state-dependent scaling. The need to choose $\{v_j\}_{j=1}^{N_V}$0 and $\{v_j\}_{j=1}^{N_V}$1-operators is identified as a major drawback, especially in non-linear settings (Berrone et al., 8 Jul 2025).

The published NAVEM line proceeds in stages. The 2023 formulation addressed lowest-order elliptic problems and showed that the projector and stabilization can be removed by learning approximations of local VEM basis functions in a harmonic polynomial space (Berrone et al., 2023). The 2024 polygonal extension introduced different training strategies, including treatments for general polygonal meshes and triangular meshes with hanging nodes, and provided a Strang-lemma-based justification for lowest-order convergence (Berrone et al., 2024). The 2025 elasticity formulation specialized the method to linear and non-linear elasticity, with a discrete nonlinear problem and Newton–Raphson tangent stiffness built directly on neural basis functions (Berrone et al., 8 Jul 2025). The 2026 work introduced two globally continuous extensions, B-NAVEM and P-NAVEM, in which exact edge continuity is enforced by construction (Berrone et al., 14 Jan 2026).

Paper	Problem class	Main contribution
"The lowest-order Neural Approximated Virtual Element Method"	Elliptic problems	First NAVEM formulation without projection or stabilization
"The lowest-order Neural Approximated Virtual Element Method on polygonal elements"	General polygonal meshes	Training strategies, polygonal generalization, hanging-node emphasis
"The Neural Approximated Virtual Element Method for Elasticity Problems"	Linear and non-linear elasticity	Discrete elasticity formulation and nonlinear tests
"Two continuous extensions of the Neural Approximated Virtual Element Method"	Continuous neural variants	Exact continuity via B-NAVEM and P-NAVEM

2. Elasticity model and variational setting

The elasticity version considers an elastic body $\{v_j\}_{j=1}^{N_V}$2, clamped on $\{v_j\}_{j=1}^{N_V}$3 and loaded by a body force $\{v_j\}_{j=1}^{N_V}$4. The displacement field $\{v_j\}_{j=1}^{N_V}$5 solves

$\{v_j\}_{j=1}^{N_V}$6

$\{v_j\}_{j=1}^{N_V}$7

$\{v_j\}_{j=1}^{N_V}$8

With

$\{v_j\}_{j=1}^{N_V}$9

the weak problem is: find $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$0 such that

$$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$1

This formulation leaves the constitutive stress $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$2 general, and the numerical experiments instantiate both linear and non-linear elasticity laws (Berrone et al., 8 Jul 2025).

The formulation is significant because it isolates the approximation mechanism from the constitutive model. In standard VEM, nonlinearity enters together with projection choices, stabilization choices, and sometimes state-dependent scaling. In NAVEM, the neural basis functions are inserted directly into the elemental weak form, so the constitutive response is evaluated against explicit approximations of the displacement field and its gradient. This is the structural reason the method is presented as particularly advantageous in handling non-linearities.

3. Neural approximation of local virtual bases

The central construction is local and element-wise. On a reference square $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$3, one defines a “rich” ansatz space

$$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$4

containing all harmonic polynomials up to some degree $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$5 together with a small number of Laplace-solved “hanging-node” functions mimicking VEM boundary behaviour. In the elasticity paper,

$$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$6

The true scalar VEM basis $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$7 is approximated by

$$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$8

where the coefficient vector $$V_h^1(E):=\{v\in[H^1(E)]^2:\Delta v=0 \text{ in }E,\; v|_e\in[P_1(e)]^2\ \forall \text{ edge } e\subset\partial E\},$$9 is predicted by a small feed-forward network,

$v(v_j)$0

$v(v_j)$1

with activation $v(v_j)$2. To reduce gradient oscillations, a second network predicts coefficients

$v(v_j)$3

for an approximation of $v(v_j)$4,

$v(v_j)$5

The associated training errors are

$v(v_j)$6

and the losses are

$v(v_j)$7

$v(v_j)$8

with $v(v_j)$9 (Berrone et al., 8 Jul 2025).

Several design features recur across the NAVEM literature. The 2023 formulation encoded element geometry by recentering and rescaling the polygon through an affine inertial map, then rotating so that the $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$0-th vertex sits at $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$1, after which the network outputs coefficients in an orthonormal basis of harmonic polynomials on a reference square (Berrone et al., 2023). The 2024 polygonal version enlarged the approximation space by adding three specially shifted “singular” harmonic functions to capture corner behavior and established that matching tangential boundary derivatives suffices to guarantee small interior gradient error, via Proposition 3.2 (Berrone et al., 2024). Across these variants, the common principle is explicit recovery of basis functions and gradients, rather than indirect computability through projectors.

4. Discrete formulation, training protocol, and assembly

In the elasticity formulation, the online phase on each element $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$2 and vertex $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$3 is:

• S.1 classify and pick the network for $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$4,
• S.2 encode $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$5,
• S.3 predict $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$6 and $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$7, then build $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$8 and $\Pi^{\nabla,E}_1:V_h^1(E)\to [P_1(E)]^2$9.

The resulting local spaces are

$L^2$0

The discrete nonlinear problem is: find $L^2$1 such that

$L^2$2

For Newton–Raphson linearization, the tangent stiffness is

$L^2$3

No extra $L^2$4 or $L^2$5-operators appear. The assembly-and-solve stage is:

• S.4 assemble $L^2$6 and $L^2$7 from the tangent form and residual equation,
• S.5 solve $L^2$8, then set $L^2$9 (Berrone et al., 8 Jul 2025).

The reported training protocol uses Python and TensorFlow, with first $\Pi^{0,E}_0$0 epochs of ADAM and then $\Pi^{0,E}_0$1 updates of self-scaled BFGS. Two datasets are employed: RDQM (random quadrilaterals) for $\Pi^{0,E}_0$2, and VM (Voronoi cells) for $\Pi^{0,E}_0$3. Training sets contain $\Pi^{0,E}_0$4 polygons per $\Pi^{0,E}_0$5. Each network uses $\Pi^{0,E}_0$6 layers $\Pi^{0,E}_0$7 $\Pi^{0,E}_0$8 neurons, $\Pi^{0,E}_0$9, hence $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$0. The offline cost is described as $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$1minutes to hours depending on dataset size, but reused for all subsequent solves. The online cost is $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$2 network evaluations, and for fine meshes the neural-network cost is negligible compared to element-wise quadratures (Berrone et al., 8 Jul 2025).

This discrete structure is closely aligned with the original 2023 formulation, where the learned space $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$3 was assembled with the exact energy form

$$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$4

yielding a purely “FEM-like” assembly of stiffness matrices on polygonal meshes (Berrone et al., 2023). The elasticity version extends that principle from linear scalar elliptic operators to nonlinear vector-valued constitutive laws.

5. Numerical behaviour in linear and non-linear elasticity

The elasticity paper reports three numerical tests. The first is a linear elasticity convergence study on $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$5 with $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$6, $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$7, and exact displacement chosen so that

$$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$8

The meshes are distorted $$a_h^E(u_h,v_h;w_h) =\int_E \sigma(x,\Pi^{0,E}_0\nabla u_h):\Pi^{0,E}_0\nabla v_h +\alpha_E(w_h)\,S_E\bigl((I-\Pi^{\nabla,E}_1)u_h,(I-\Pi^{\nabla,E}_1)v_h\bigr),$$9 quadrilaterals. Errors are measured through

$S_E$0

with $S_E$1 and $S_E$2 defined via $S_E$3 and $S_E$4. Both methods achieve $S_E$5 in $S_E$6 and $S_E$7 in $S_E$8, and the NAVEM constants are $S_E$9 VEM (Berrone et al., 8 Jul 2025).

The second test is a benchmark nonlinear model with

$\alpha_E(w_h)$0

on $\alpha_E(w_h)$1, clamped on the full boundary, with moderate-deformation cases

$\alpha_E(w_h)$2

Structured squares and distorted Voronoi meshes are used. In the large-deformation case, VEM requires incremental force and careful stabilization, either norm-based or trace-based, to avoid loss of convergence or spurious oscillations. NAVEM works robustly without increments or stabilization scaling, with lower error constants and fewer Newton steps (Berrone et al., 8 Jul 2025).

The third test replaces $\alpha_E(w_h)$3 by the First Piola–Kirchhoff stress of neo-Hookean hyperelasticity,

$\alpha_E(w_h)$4

on $\alpha_E(w_h)$5, with $\alpha_E(w_h)$6 clamped on $\alpha_E(w_h)$7 and $\alpha_E(w_h)$8. On an unstructured polygonal mesh of $\alpha_E(w_h)$9 elements, VEM must choose a projection of $\{v_j\}_{j=1}^{N_V}$00 and a stabilization; the reported VEM solution exhibits visible spurious oscillations in the deformed displacement. NAVEM evaluates $\{v_j\}_{j=1}^{N_V}$01 pointwise, requires no stabilization, and yields a smooth, oscillation-free solution comparable to a fine FEM reference (Berrone et al., 8 Jul 2025).

Taken together, these tests support three empirical statements emphasized in the elasticity work: the absence of a user-tuned stabilization term greatly simplifies nonlinear iterations; NAVEM attains the same convergence rates as standard VEM/FEM but with smaller error constants; and NAVEM is more robust in large-deformation regimes, with no need for incremental loading or specialized scaling of stabilization.

6. Theoretical extensions, continuity, and open issues

The broader NAVEM literature clarifies both what the method is and what it is not. In the 2024 polygonal formulation, the global NAVEM trial/test space is assembled from learned local surrogates $\{v_j\}_{j=1}^{N_V}$02, and Strang’s lemma for broken spaces yields

$\{v_j\}_{j=1}^{N_V}$03

where

$\{v_j\}_{j=1}^{N_V}$04

Under decay of the boundary-fit error, one recovers the standard lowest-order rates

$\{v_j\}_{j=1}^{N_V}$05

The same paper also reports that NAVEM is free of delicate stabilization tuning in a strongly anisotropic diffusion test and converges in fewer Newton iterations in a nonlinear diffusion problem (Berrone et al., 2024).

A frequent misconception is to treat NAVEM as merely a learned stabilization of VEM. The published formulations describe a stronger modification: the abstract virtual spaces and associated projectors and stabilizations are replaced by explicit neural-network-based basis functions in a known approximation space. Another misconception concerns continuity. Standard NAVEM only enforces boundary linearity approximately, through losses on values and tangential derivatives along edges. By contrast, the continuous 2026 extensions introduce a boundary operator

$\{v_j\}_{j=1}^{N_V}$06

where $\{v_j\}_{j=1}^{N_V}$07 is a bubble vanishing on $\{v_j\}_{j=1}^{N_V}$08 and $\{v_j\}_{j=1}^{N_V}$09 is a transfinite interpolant equal to the VEM basis on $\{v_j\}_{j=1}^{N_V}$10. This construction enforces exact Dirichlet data on $\{v_j\}_{j=1}^{N_V}$11 and exact continuity across adjacent elements (Berrone et al., 14 Jan 2026).

The two continuous variants pursue different principles. B-NAVEM uses a PINN approximation of the local Laplace problem

$\{v_j\}_{j=1}^{N_V}$12

so it remains closer in spirit to the classical harmonic virtual basis. P-NAVEM abandons harmonicity and instead enforces partition of unity and linear reproduction,

$\{v_j\}_{j=1}^{N_V}$13

The reported trade-off is explicit: B-NAVEM preserves approximate harmonicity inside but is $\{v_j\}_{j=1}^{N_V}$14 slower in wall-time during training because automatic differentiation of Laplacians is expensive, whereas P-NAVEM is the most accurate on concave-polygon meshes and does not approximate the true VEM space (Berrone et al., 14 Jan 2026).

The limitations and open directions are also stated plainly in the polygonal literature. Offline training requires representative polygon samples; extension to higher polynomial orders needs new approximation spaces; and rigorous a priori bounds on $\{v_j\}_{j=1}^{N_V}$15 versus network size remain to be studied. Possible extensions include multi-patch adaptivity, unsteady PDEs, and coupling NAVEM with domain decomposition or multigrid (Berrone et al., 2024). For elasticity in particular, the existing evidence indicates that the most distinctive advantage of NAVEM lies not in changing the weak form, but in simplifying the discrete mechanics of polygonal non-linear analysis by removing the projection and stabilization machinery that standard VEM requires.