---
title: Neural Approximated Virtual Element Method
url: https://www.emergentmind.com/topics/neural-approximated-virtual-element-method
type: topic
---

# Neural Approximated Virtual Element Method

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 [2507.05786]. 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 [2311.18534] [2409.15917] [2601.09595].

## 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 \(S_E\) and \(\Pi\)-operators is identified as a major drawback, especially in non-linear settings [2507.05786].

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 [2311.18534]. 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 [2409.15917]. 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 [2507.05786]. The 2026 work introduced two globally continuous extensions, B-NAVEM and P-NAVEM, in which exact edge continuity is enforced by construction [2601.09595].

| 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 \(\Omega\subset\mathbb R^2\), clamped on \(\Gamma_D\subset\partial\Omega\) and loaded by a body force \(f\). The displacement field \(u:\Omega\to\mathbb R^2\) solves
\[
-\nabla\cdot\sigma(x,\nabla u(x))=f(x)\quad\text{in }\Omega,
\]
\[
u=0\quad\text{on }\Gamma_D,
\]
\[
\sigma(x,\nabla u)\cdot n=0\quad\text{on }\Gamma_N:=\partial\Omega\setminus\Gamma_D.
\]
With
\[
V:=H^1_{0,\Gamma_D}(\Omega;\mathbb R^2),
\]
the weak problem is: find \(u\in V\) such that
\[
a(u,v):=\int_\Omega \sigma(x,\nabla u):\nabla v = \int_\Omega f\cdot v
\qquad \forall v\in V.
\]
This formulation leaves the constitutive stress \(\sigma(x,\nabla u)\) general, and the numerical experiments instantiate both linear and non-linear elasticity laws [2507.05786].

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 \(\widehat E\), one defines a “rich” ansatz space
\[
\mathcal A_{j,E}= \operatorname{span}\{h_k\}_{k=1}^{M},
\]
containing all harmonic polynomials up to some degree \(\ell\) together with a small number of Laplace-solved “hanging-node” functions mimicking VEM boundary behaviour. In the elasticity paper,
\[
M=2\ell+4.
\]
The true scalar VEM basis \(\phi_{j,E}\) is approximated by
\[
\widehat \phi_{j,E}(x)=\sum_{k=1}^M c_k^{\phi}(E,j)\,h_k(x),
\]
where the coefficient vector \(c^{\phi}(E,j)\in\mathbb R^M\) is predicted by a small feed-forward network,
\[
x_0=\operatorname{encode}(E,v_j),\qquad
x_\ell=\rho(W_\ell x_{\ell-1}+b_\ell),\ \ \ell=1,\dots,L-1,
\]
\[
c^{\phi}(E,j)=W_Lx_{L-1}+b_L,
\]
with activation \(\rho(t)=\tanh(t)\). To reduce gradient oscillations, a second network predicts coefficients
\[
c^q(E,j)=Net^q(\operatorname{encode}(E,v_j))
\]
for an approximation of \(q_{j,E}:=\nabla\phi_{j,E}\),
\[
\widehat q_{j,E}(x)=\sum_{k=1}^M c_k^q(E,j)\,\nabla h_k(x).
\]
The associated training errors are
\[
\varepsilon_{j,E}^{\phi}=\|\widehat \phi_{j,E}-\phi_{j,E}\|_{H^{1/2}(\partial E)},
\qquad
\varepsilon_{j,E}^{q}=\|(\widehat q_{j,E}-q_{j,E})\cdot t\|_{L^2(\partial E)},
\]
and the losses are
\[
\mathcal L_\phi=
\sqrt{
\frac{1}{|\mathcal T|\,N_V}
\sum_{E\in\mathcal T}\sum_{j=1}^{N_V}(\varepsilon_{j,E}^{\phi})^2
}
+\lambda\|\theta\|^2,
\]
\[
\mathcal L_q=
\sqrt{
\frac{1}{|\mathcal T|\,N_V}
\sum_{E\in\mathcal T}\sum_{j=1}^{N_V}(\varepsilon_{j,E}^{q})^2
}
+\lambda\|\theta\|^2,
\]
with \(\lambda=10^{-8}\) [2507.05786].

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 \(j\)-th vertex sits at \((1,0)\), after which the network outputs coefficients in an orthonormal basis of harmonic polynomials on a reference square [2311.18534]. 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 [2409.15917]. 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 \(E\) and vertex \(j\) is:

- **S.1** classify and pick the network for \(N_V\),
- **S.2** encode \(x_0=\operatorname{encode}(E,v_j)\),
- **S.3** predict \(c^\phi\) and \(c^q\), then build \(\widehat\phi_{j,E}\) and \(\widehat q_{j,E}\).

The resulting local spaces are
\[
\widetilde V_h^1(E)=\operatorname{span}\{\widehat \phi_{j,E}\}_{j=1}^{N_V},
\qquad
\widetilde V_h(E)=[\widetilde V_h^1(E)]^2,
\qquad
\dim \widetilde V_h(E)=2N_V.
\]
The discrete nonlinear problem is: find \(u_h^{NN}\in\widetilde V_h\) such that
\[
\sum_{E\in T_h}\int_E \sigma(x,\nabla u_h^{NN}):\nabla v_h^{NN}
=
\sum_{E\in T_h}\int_E f\cdot v_h^{NN}
\qquad \forall v_h^{NN}\in\widetilde V_h.
\]
For Newton–Raphson linearization, the tangent stiffness is
\[
D\mathcal G_E(u_h^k,v_h^{NN})[\delta_h]
=
\int_E (A(x,\nabla u_h^k):\nabla\delta_h):\nabla v_h^{NN},
\qquad A=\frac{\partial \sigma}{\partial \nabla u}.
\]
No extra \(S_E\) or \(\Pi\)-operators appear. The assembly-and-solve stage is:

- **S.4** assemble \(K^k\) and \(r^k\) from the tangent form and residual equation,
- **S.5** solve \(K^k\delta^k=-r^k\), then set \(u^{k+1}=u^k+\delta^k\) [2507.05786].

The reported training protocol uses Python and TensorFlow, with first \(5\,000\) epochs of ADAM and then \(5\,000\) updates of self-scaled BFGS. Two datasets are employed: RDQM (random quadrilaterals) for \(N_V=4\), and VM (Voronoi cells) for \(N_V\ge 5\). Training sets contain \(O(10^3)\) polygons per \(N_V\). Each network uses \(5\) layers \(\times\) \(50\) neurons, \(\ell=20\), hence \(M=44\). The offline cost is described as \(\sim\)minutes to hours depending on dataset size, but reused for all subsequent solves. The online cost is \(O(N_V\cdot \#\text{elements})\) network evaluations, and for fine meshes the neural-network cost is negligible compared to element-wise quadratures [2507.05786].

This discrete structure is closely aligned with the original 2023 formulation, where the learned space \(\widehat V_h^1(E)\) was assembled with the exact energy form
\[
a_{h,NN}^E(u_h,v_h)=\int_E \nabla u_h\cdot \nabla v_h,
\]
yielding a purely “FEM-like” assembly of stiffness matrices on polygonal meshes [2311.18534]. 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 \(\Omega=(0,1)^2\) with \(\mu=1.5\), \(\lambda=3\), and exact displacement chosen so that
\[
u_1=16x_1(1-x_1)x_2(1-x_2)+1.1,\qquad u_2=5u_1.
\]
The meshes are distorted \(4\times4,\dots,32\times32\) quadrilaterals. Errors are measured through
\[
err_0^{NN}=\sqrt{\sum_E\|u-u_h^{NN}\|^2_{L^2(E)}},
\qquad
err_1^{NN}=\sqrt{\sum_E\|\nabla u-\nabla u_h^{NN}\|^2_{L^2(E)}},
\]
with \(err_0^{VEM}\) and \(err_1^{VEM}\) defined via \(\Pi_0^{0,E}\) and \(\Pi_1^{0,E}\). Both methods achieve \(O(h^2)\) in \(L^2\) and \(O(h)\) in \(H^1\), and the NAVEM constants are \(\lesssim\) VEM [2507.05786].

The second test is a benchmark nonlinear model with
\[
\sigma=(1+\|\varepsilon(u)\|^2)\varepsilon(u)
\]
on \(\Omega=(0,1)^2\), clamped on the full boundary, with moderate-deformation cases
\[
u=C\,x_1(1-x_1)x_2(1-x_2)[1,1]^T,\qquad C=1 \text{ and } C=80.
\]
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 [2507.05786].

The third test replaces \(\sigma\) by the First Piola–Kirchhoff stress of neo-Hookean hyperelasticity,
\[
P(F)=\mu(F-F^{-T})+\lambda\,J\log J\,F^{-T},\qquad J=\det(F),
\]
on \(\Omega=(0,1)^2\), with \(u\) clamped on \(x=0\) and \(f=[100x_2^3,0]\). On an unstructured polygonal mesh of \(\sim 400\) elements, VEM must choose a projection of \(J\) and a stabilization; the reported VEM solution exhibits visible spurious oscillations in the deformed displacement. NAVEM evaluates \(J\) pointwise, requires no stabilization, and yields a smooth, oscillation-free solution comparable to a fine FEM reference [2507.05786].

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 \(\widetilde V_h^1(E)=\operatorname{span}\{\varphi_{j,E}^{NN}\}\), and Strang’s lemma for broken spaces yields
\[
\|u-u_h^{NN}\|_{H^1_h}\le C\bigl(h+h^{-2}\varepsilon\bigr),
\]
where
\[
\varepsilon=\max_{E,j}\|\varphi_{j,E}-\varphi_{j,E}^{NN}\|_{H^{1/2}(\partial E)}.
\]
Under decay of the boundary-fit error, one recovers the standard lowest-order rates
\[
\|u-u_h^{NN}\|_{H^1(\Omega)}=O(h),
\qquad
\|u-u_h^{NN}\|_{L^2(\Omega)}=O(h^2).
\]
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 [2409.15917].

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
\[
\mathcal B_{j,E}(v)(x)=\psi_E^0(x)\,v(x)+\psi_{j,E}(x),
\]
where \(\psi_E^0\) is a bubble vanishing on \(\partial E\) and \(\psi_{j,E}\) is a transfinite interpolant equal to the VEM basis on \(\partial E\). This construction enforces exact Dirichlet data on \(\partial E\) and exact continuity across adjacent elements [2601.09595].

The two continuous variants pursue different principles. B-NAVEM uses a PINN approximation of the local Laplace problem
\[
\Delta \varphi_{j,E}=0 \text{ in } E,\qquad \varphi_{j,E}=\psi_{j,E} \text{ on } \partial E,
\]
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,
\[
\sum_j \varphi_{j,E}(x)\equiv 1,\qquad
\sum_j (\mathbf v_j)_1\,\varphi_{j,E}(x)\approx x_1,\qquad
\sum_j (\mathbf v_j)_2\,\varphi_{j,E}(x)\approx x_2.
\]
The reported trade-off is explicit: B-NAVEM preserves approximate harmonicity inside but is \(\simeq 8\times\) 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 [2601.09595].

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 \(\varepsilon\) versus network size remain to be studied. Possible extensions include multi-patch adaptivity, unsteady PDEs, and coupling NAVEM with domain decomposition or multigrid [2409.15917]. 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.

Source: https://www.emergentmind.com/topics/neural-approximated-virtual-element-method