<figure><img src="image_1.png" /></figure>

HAL
open science

# The Neural Approximated Virtual Element Method for Elasticity Problems

Stefano Berrone, Moreno Pintore, Gioana Teora

► To cite this version:

Stefano Berrone, Moreno Pintore, Gioana Teora. The Neural Approximated Virtual Element Method for Elasticity Problems. 2025. hal-05150590

HAL Id: hal-05150590

https://hal.sorbonne-universite.fr/hal-05150590v1

Preprint submitted on 8 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THE NEURAL APPROXIMATED VIRTUAL ELEMENT METHOD FOR ELASTICITY PROBLEMS

## A PREPRINT

<figure><img src="image_3.png" /></figure>

Stefano Berrone
Dipartimento di Scienze Matematiche
“G. L. Lagrange”
Politecnico di Torino, TO, 10129
stefano.berrone@polito.it

<figure><img src="image_4.png" /></figure>

Moreno Pintore
Laboratoire Jacques-Louis Lions,
Sorbonne Université,
MEGAVOLT Team, Inria,
4 place Jussieu, 75005 Paris, France
moreno.pintore@sorbonne-universite.fr

<figure><img src="image_5.png" /></figure>

Gioana Teora
Dipartimento di Scienze Matematiche
“G. L. Lagrange”
Politecnico di Torino, TO, 10129
gioana.teora@polito.it

July 8, 2025

## ABSTRACT

We present the Neural Approximated Virtual Element Method to numerically solve elasticity problems. This hybrid technique combines 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. We present the discrete formulation of the problem and provide numerical tests on both linear and non-linear elasticity problems, demonstrating the advantages of having a simple discretization, particularly in handling non-linearities.

**Keywords** NAVEM · Elasticity · Virtual Element Method · Neural Network · Basis Functions · Polygonal Meshes

## 1 Introduction

Discretization methods relying on meshes comprising general polygons are receiving a great deal of attention from the scientific community. Indeed, polytopal methods can be very useful for a large variety of reasons, including the automatic handle of hanging nodes, the use of moving meshes, the ease of meshing complex geometries, and adaptivity [1, 2, 3]. Moreover, these methods have been successfully applied to several fields of computational mechanics, revealing several advantages over methods employing classical triangular/tetrahedral and quadrilateral elements. For instance, polygonal methods have shown great advantages in handling crack propagation and branching [4], contact problems [5], and topology optimization [6].

In this framework, an important polygonal method is the Virtual Element Method (VEM), introduced in [7, 8] for elliptic problems and then extended to elastic problems in [9, 10]. VEM can be seen as a generalization of the standard Finite Element Method (FEM) [11, 12] on more general meshes. One of the main features of the VEM, is that the local basis functions are not known in closed form. Indeed, even if they can be defined as the solution of local differential problems, they are never explicitly computed or approximated in the virtual element framework. It is therefore necessary to introduce suitable projection and stabilization operators to derive computable discrete bilinear forms associated with the Partial Differential Equation (PDE) of interest. In particular, the usage of a stabilization term is an undesirable aspect of the VEM method, as it lacks a theoretical foundation, remains problem-dependent, and an improper choice
can significantly affect the method accuracy [13]. Moreover, when dealing with non-linearities, the presence of the stabilization term makes the construction of the method more involved and increases the computational effort needed to solve the problem. Finally, the need for a projector to access the point-wise evaluation of virtual functions makes the computation of post-process quantities more difficult [14].

Because of these reasons, in the last few years, several numerical methods have emerged to address these issues while preserving the main advantages of the method. A non-exhaustive list of such methods is the following. In [15, 16], a stabilization-free formulation of the method, derived by enlarging the local projection polynomial order and enhancing the function spaces according to the polygon geometry, is presented. In [17, 18], two approaches are proposed to locally approximate the local virtual basis functions and to get rid of any stabilization or projection operators: one based on the reduced basis method [19] and the other on the Laplace solver introduced in [20]. We observe that, in these last two methods, the basis functions are explicitly computed on each element, meaning that the “virtual” attribute only refers to the underlying discrete space.

Within the same context, the Neural Approximated Virtual Element Method (NAVEM) is introduced in [21, 22] for general elliptic problems. The main idea of this novel method is to use neural networks [23] to accurately predict on the fly the local virtual basis functions, which are then used like in a standard FEM framework to compute a discrete solution. For this reason, the NAVEM method falls within the general recent framework of Scientific Machine Learning (SciML) [24], by combining the advantages of the non-linear approximation capabilities of neural networks with the ones of more classical numerical methods. This framework is exponentially increasing in popularity because it represents the intersection between two active communities that rarely interacted before. We highlight that several libraries are available to easily implement and adapt complex neural networks. The most used are Tensorflow [25], PyTorch [26], and JAX [27].

In this paper, we present the NAVEM for linear and non-linear elastic problems. Moreover, we further highlight some new key aspects of the neural network that allow us to improve the NAVEM accuracy. It is very important to notice that the neural networks used in the NAVEM framework are trained using only the geometric information of the elements. Therefore, after a potentially expensive training phase, the same neural networks can be used multiple times to solve problems characterized by different material coefficients, constitutive laws, body forces, initial configurations but also, and most importantly, on various meshes. In this sense, it can be regarded as a method based on the so-called offline-online splitting.

The manuscript is organized as follows. In Section 2, we present the model problem and its VEM discretization, describing its main properties, which are essential to devise a proper neural network training strategy, and highlighting its main issues. In Section 3, we introduce the NAVEM for elasticity problems. In particular, in Section 3.1, we construct the function space from where the neural networks select the local basis functions in the NAVEM discretization, in Sections 3.2 and 3.3, we describe the neural networks architecture and training, respectively, and in Section 3.4 we discuss the general structure of the NAVEM as a hybrid method mixing ideas from the machine learning and scientific computing communities. Later, in Section 4, we prove through numerical experiments the effectiveness of the method for elasticity problems, and the advantages of its simple formulation. The conclusion of the present work and the future perspectives are discussed in Section 5.

## 2 The Model Problem and The Virtual Element Method

Let us consider a homogeneous isotropic elastic body $\Omega \subset \mathbb{R}^2$ with boundary $\Gamma$. We assume that this body is clamped on $\Gamma_D \subseteq \Gamma$, $|\Gamma_D| \neq 0$, and subjected to a body load $f \in [L^2(\Omega)]^2$.

In the case of static problems, the governing equation of the solid in the initial configuration is

$$ \begin{cases} -\nabla \cdot \boldsymbol{\sigma}(\boldsymbol{x}, \nabla \boldsymbol{u}) = \boldsymbol{f}(\boldsymbol{x}) & \text{in } \Omega, \\ \boldsymbol{u} = \mathbf{0} & \text{on } \Gamma_D, \\ \boldsymbol{\sigma}(\boldsymbol{x}, \nabla \boldsymbol{u})\boldsymbol{n} = \mathbf{0} & \text{on } \Gamma_N := \partial\Omega \setminus \Gamma_D, \end{cases} \quad (1) $$

where $\boldsymbol{u} : \Omega \to \mathbb{R}^2$ represents the body displacement, whereas the second-order tensor $\boldsymbol{\sigma} = \boldsymbol{\sigma}(\boldsymbol{x}, \nabla \boldsymbol{u})$ represents a constitutive law for the material at each point $\boldsymbol{x} \in \Omega$.

Let us consider the space

$$ \mathbf{V} = [H_{0,\Gamma_D}^1(\Omega)]^2 := \{\boldsymbol{v} \in [H^1(\Omega)]^2 : \boldsymbol{v} = \mathbf{0} \text{ on } \Gamma_D\} $$

and the bilinear form

$$ a(\boldsymbol{u}, \boldsymbol{v}) = \int_{\Omega} \boldsymbol{\sigma}(\boldsymbol{x}, \nabla \boldsymbol{u}(\boldsymbol{x})) : \nabla \boldsymbol{v}, \quad (2) $$
where, given two second-order tensors $\sigma$ and $T$, we define $\sigma : T = \sum_{i,j=1}^{n} \sigma_{ij} T_{ij}$. The variational formulation of Problem (1) reads as: Find $u \in V$ such that

$$a(u, v) = \int_{\Omega} f \cdot v \quad \forall v \in V. \qquad (3)$$

We remark that here we are not considering internal constraints, such as incompressibility, which usually require suitable treatments.

## 2.1 The Virtual Element Space

Let us consider a tessellation $T_h$ of the domain $\Omega$ made up of polygonal elements $E$. We denote by $h_E$, $N_E^v$, $x_E$ and $\mathcal{E}_{h,E}$ the diameter, the number of vertices/edges, the centroid and the set of edges of $E$, respectively. We further set, as usual, $h = \max_{E \in \mathcal{T}_h} h_E$.

The virtual element space for the displacement is given by [9]:

$$\mathcal{V}_{h,1}(E) = \left\{ \begin{array}{l} v \in [H^1(E)]^2 : (i) \Delta v = 0 \\ (ii) \quad v|_e \in [\mathbb{P}_1(e)]^2, \forall e \in \mathcal{E}_{h,E} \text{ and } v \in [C^0(\partial E)]^2 \end{array} \right\},$$

where $\Delta v = \begin{bmatrix} \Delta v_1 \\ \Delta v_2 \end{bmatrix}$ is the component-wise Laplacian operator. The degrees of freedom for a vector-valued function $v \in \mathcal{V}_{h,1}(E)$ are the values of $v$ at the vertices of the element $E$. Thus, we have

$$N_E^{\text{dof}} = 2N_E^v.$$

Let us define the set of Lagrange basis functions $\{\varphi_{i,E}\}_{i=1}^{N_E^{\text{dof}}}$ related to the aforementioned degrees of freedom. It is easy to check that

$$\varphi_{i,E} = \begin{cases} \begin{bmatrix} \varphi_{i,E} \\ 0 \end{bmatrix} & \text{if } i = 1, \dots, N_E^v, \\ \begin{bmatrix} 0 \\ \varphi_{i-N_E^v,E} \end{bmatrix} & \text{if } i = N_E^v + 1, \dots, N_E^{\text{dof}}, \end{cases} \qquad (4)$$

where $\{\varphi_{j,E}\}_{j=1}^{N_E^v}$ corresponds to the set of Lagrange basis functions related to degrees of freedom for the scalar virtual element space [7]:

$$\mathcal{V}_{h,1}(E) = \left\{ \begin{array}{l} v \in H^1(E) : (i) \Delta v = 0 \\ (ii) \quad v|_e \in \mathbb{P}_1(e), \forall e \in \mathcal{E}_{h,E} \text{ and } v \in C^0(\partial E) \end{array} \right\}. \qquad (5)$$

Indeed, it holds $\mathcal{V}_{h,1}(E) = [\mathcal{V}_{h,1}(E)]^2$.

## 2.2 The Virtual Element Discretization

For each polygon $E$, we introduce the *computable* local polynomial projector $\Pi_1^{\nabla,E} : [H^1(E)]^2 \to [\mathbb{P}_1(E)]^2$, which is defined as

$$\left\{ \begin{array}{l} \left( \nabla v - \nabla \Pi_1^{\nabla,E} v, \nabla p \right)_E = 0 \quad \forall p \in [\mathbb{P}_1(E)]^2, \\ \int_{\partial E} \Pi_1^{\nabla,E} v = \int_{\partial E} v \quad \forall v \in [H^1(E)]^2. \end{array} \right.$$

Moreover, we denote by $\Pi_0^{0,E}$ both the vector-valued and the tensor-valued $L^2$-projectors onto constants. In particular, $\forall v \in [H^1(E)]^2$, it holds [28]

$$\Pi_0^{0,E} \nabla v = \nabla \Pi_1^{\nabla,E} v,$$

and

$$\Pi_0^{0,E} \mathrm{div} v = \Pi_0^{0,E} \mathrm{tr}(\nabla v) = \mathrm{tr}(\Pi_0^{0,E} \nabla v).$$
Specifically, by computable we mean that the projection of virtual functions can be computed only through the information provided by the degrees of freedom of such a function. Now, let us split the continuous bilinear form (2) according to $\mathcal{T}_h$, i.e.

$$a(\mathbf{u}, \mathbf{v}) = \sum_{E \in \mathcal{T}_h} a^E(\mathbf{u}, \mathbf{v}), \quad \text{with} \quad a^E(\mathbf{u}, \mathbf{v}) = \int_E \sigma(\mathbf{x}, \nabla \mathbf{u}) : \nabla \mathbf{v}. \qquad (6)$$

For all $E \in \mathcal{T}_h$, we define the discrete virtual bilinear form as [10]:

$$a_h^E(\mathbf{u}_h, \mathbf{v}_h; \mathbf{w}_h) := \int_E \sigma(\mathbf{x}, \Pi_0^{0,E} \nabla \mathbf{u}_h) : \Pi_0^{0,E} \nabla \mathbf{v}_h \\ + \alpha_E(\mathbf{w}_h) S^E((I - \Pi_1^{\nabla,E})\mathbf{u}_h, (I - \Pi_1^{\nabla,E})\mathbf{v}_h). \qquad (7)$$

where the stability term $S^E(\cdot, \cdot)$ is chosen in such a way it scales like $a^E(\cdot, \cdot)$ for a unitary material constants on the kernel of $\Pi_1^{\nabla,E}$. On the other hand, the stabilizing parameter $\alpha_E(\mathbf{w}_h)$ takes into account different material constants as well as non-linear materials.

Finally, to solve Problem (7), we must address the issue of non-linearity. Among the commonly used numerical methods for solving systems of non-linear equations, the Newton-Raphson method is the fastest, offering a quadratic convergence provided that the initial estimate is sufficiently close to the solution. Therefore, selecting an initial estimate close to the solution is crucial for accelerating the convergence. In solid mechanics, the initial estimate is usually set to the undeformed shape of the structure, i.e. the initial displacement is set to the all-zeros vector. Since the initial estimate usually starts from the zero displacement, the Newton-Raphson method converges quickly when the applied load is small. However, convergence becomes challenging when large applied loads induce large displacements. The incremental force method addresses this issue by applying the load in increments. Within each increment, the Newton-Raphson procedure remains unchanged. A new load increment is applied after the solution corresponding to the previous load increment has successfully converged [29].

Thus, as in [10], for each $n = 1, \dots, N$, we solve the following problem: Given $\mathbf{w}_h \in \mathbf{V}_{h,1}$, find $\mathbf{u}_h^n \in \mathbf{V}_{h,1}$ such that

$$\sum_{E \in \mathcal{T}_h} a_h^E(\mathbf{u}_h^n, \mathbf{v}_h; \mathbf{w}_h) = \sum_{E \in \mathcal{T}_h} \left( \mathbf{f}^n, \Pi_0^{0,E} \mathbf{v}_h \right)_E \quad \forall \mathbf{v}_h \in \mathbf{V}_{h,1}, \qquad (8)$$

where $\mathbf{f}^n = \frac{n}{N}\mathbf{f}$ and $\mathbf{u}_h^0 = \mathbf{0}$ on $\Omega$. For each $n$, Problem (8) is solved with the Newton-Raphson method. Setting

$$\mathcal{G}_h^E(\mathbf{u}_h^n, \mathbf{v}_h; \mathbf{w}_h) = a_h^E(\mathbf{u}_h^n, \mathbf{v}_h; \mathbf{w}_h) - \left( \mathbf{f}^n, \Pi_0^{0,E} \mathbf{v}_h \right)_E,$$

we introduce

$$D\mathcal{G}_h(\mathbf{u}_h^{n,k}, \mathbf{v}_h; \mathbf{w}_h) \cdot [\delta_h^{n,k}] = \int_E \left( \mathcal{A}(\mathbf{x}, \Pi_0^{0,E} \nabla \mathbf{u}_h^{n,k}) : \nabla \Pi_0^{0,E} \delta_h^{n,k} \right) : \Pi_0^{0,E} \nabla \mathbf{v}_h \\ + \alpha_E(\mathbf{w}_h) S^E((I - \Pi_1^{\nabla,E})\delta_h^{n,k}, (I - \Pi_1^{\nabla,E})\mathbf{v}_h), \qquad (9)$$

where the elastic modulus $\mathcal{A}(\mathbf{x}, \nabla \mathbf{u})$ is defined as the fourth-order tensor

$$\mathcal{A}(\mathbf{x}, \nabla \mathbf{u}) = \frac{\partial \sigma(\mathbf{x}, \nabla \mathbf{u})}{\partial \nabla \mathbf{u}}. \qquad (10)$$

Finally, as in [30], we generate a sequence $\mathbf{u}_h^{n,k+1} = \mathbf{u}_h^{n,k} + \delta_h^{n,k}$ for each $k \ge 0$ by solving: Given $\mathbf{w}_h \in \mathbf{V}_{h,1}$, find $\delta_h^{n,k} \in \mathbf{V}_{h,1}$ such that

$$\sum_{E \in \mathcal{T}_h} \left[ \mathcal{G}_h^E(\mathbf{u}_h^{n,k}, \mathbf{v}_h; \mathbf{w}_h) + D\mathcal{G}_h^E(\mathbf{u}_h^{n,k}, \mathbf{v}_h; \mathbf{w}_h) \cdot [\delta_h^{n,k}] \right] = 0 \quad \forall \mathbf{v}_h \in \mathbf{V}_{h,1}.$$

Usually, $\mathbf{w}_h$ is set equal to the solution found at the last loading increment step, i.e. $\mathbf{u}_h^{n-1}$. This choice is taken to simplify the Newton-Raphson application by avoiding the computation of derivatives of $\alpha_E$ in the related tangent matrix (9). Indeed, although $\mathbf{w}_h = \mathbf{u}_h^n$ is a more intuitive choice, it is also more computationally demanding [10, 14]. Other possible choices include setting $\mathbf{w}_h = \mathbf{0}$, i.e. evaluating $\alpha_E$ in the undeformed configuration. However, this might lead to unsatisfactory results, also in the small deformations regime [14]. Moreover, different choices for $\alpha_E$ have been introduced:
• norm-based stabilization [10]:

$$ \alpha_E(\boldsymbol{w}_h) = \left\| \mathbf{A}(\boldsymbol{x}, \Pi_0^{0,E} \nabla \boldsymbol{w}_h) \right\|, $$

which is chosen to produce a strictly positive stabilization term. The operator $|||\cdot|||$ could be chosen as any fourth-order tensor norm.

• trace-based stabilization [14]:

$$ \alpha_E(\boldsymbol{w}_h) = \frac{1}{4} \mathrm{tr} \left( \mathbf{A}(\boldsymbol{x}, \Pi_0^{0,E} \nabla \boldsymbol{w}_h) \right). $$

We observe that, with this choice, the stabilization may assume negative values.

• stiffness-based stabilization [31]:

$$ \alpha_E(\boldsymbol{w}_h) = \frac{1}{4N^{\mathrm{dof}}} \left( \sum_{i=1}^{N^{\mathrm{dof}}} ((\mathbf{K}_T^C)_{ii})^2 \right)^{\frac{1}{2}} \quad \text{or} \quad \alpha_E(\boldsymbol{w}_h) = \frac{1}{4N^{\mathrm{dof}}} \sum_{i=1}^{N^{\mathrm{dof}}} (\mathbf{K}_T^C)_{ii}, $$

where $\mathbf{K}_T^C$ is the global tangent system matrix related to the consistency term.

# 3 The Neural Approximated Virtual Element Method

In the NAVEM method, we propose to approximate $\varphi_{j,E} \in \mathcal{V}_{h,1}(E)$ by means of a neural network-based approach and to assemble the final linear system as in the more classical Finite Element Method [22].

To do so, we approximate $\varphi_{j,E}$ by a function $\varphi_{j,E}^{NN}$ in a suitable functional space $\mathcal{H}_{j,E}^{NN}$. A similar procedure is repeated to approximate the gradient of $\varphi_{j,E}$ (or, more generally, all the required derivatives). Finally, since the approximation $\varphi_{j,E}^{NN}$ of $\varphi_{j,E}$ is known in closed form as well as its derivatives, it can be evaluated point-wise inside each element $E \in \mathcal{T}_h$ to numerically compute the integrals in (3), without introducing any projection or stability operator.

## 3.1 The Local Approximation Space

In this section, we describe the local function space $\mathcal{H}_{j,E}^{NN}$ in which the approximation $\varphi_{j,E}^{NN}$ of $\varphi_{j,E}$ is sought. Such space coincides with the one adopted in [22]. Here, for the sake of clarity, we briefly recall the main concepts and refer to [22] for a more detailed description.

All the functions in $\mathcal{H}_{j,E}^{NN}$ are defined on a reference squared region $S^{NN} = [-R^{NN}, R^{NN}]^2 \subset \mathbb{R}^2$ independent of $E$. Therefore, we assume that a suitable affine map is available to map each element $E \in \mathcal{T}_h$ to a reference element $\hat{E} \subset S^{NN}$. A convenient map is the one described in [32], but other choices are admissible.

Let $\mathbb{H}_{\ell^{NN}}(S^{NN})$ be the space of harmonic polynomials up to order $\ell^{NN} \ge 0$. A basis of such a space is the set of scaled polynomials

$$ \{1, \Re\left(\left(\frac{z}{R^{NN}}\right)^{\ell}\right), \Im\left(\left(\frac{z}{R^{NN}}\right)^{\ell}\right), \ell = 1, \dots, \ell^{NN}\}. \quad (11) $$

Here $z$ is the complex number $z = x_1 + ix_2$ associated with the point $\boldsymbol{x} = [x_1, x_2]^T \in \mathbb{R}^2$, and $\Re$ and $\Im$ denote the real and imaginary part of a complex quantity. Note that all the polynomials in (11) can be recursively retrieved as in [33].

As discussed in [22], harmonic polynomials are not enough to accurately approximate the VEM basis functions, we thus introduce additional functions mimicking the target virtual functions. Let us introduce the Laplace problem

$$ \left\{ \begin{array}{ll} \Delta \tilde{\Phi} = 0 & \text{in } \Omega_{\Phi} = (-1, 1)^2, \\ \tilde{\Phi} = 1 + x_2 & \text{on } \Gamma_{\Phi,1} = \{x_1 = 1 \text{ and } -1 \le x_2 \le 0\}, \\ \tilde{\Phi} = 1 - x_2 & \text{on } \Gamma_{\Phi,2} = \{x_1 = 1 \text{ and } 0 \le x_2 \le 1\}, \\ \tilde{\Phi} = 0 & \text{on } \partial \Omega_{\Phi} \setminus \{\Gamma_{\Phi,1} \cup \Gamma_{\Phi,2}\}, \end{array} \right. $$

and approximate its solution with a function $\Phi$ of the form

$$ \Phi(z) = \sum_{\alpha=1}^{N^1} c_{\alpha}^{1} \Re \left( \frac{d_{\alpha}}{z - z_{\alpha}} \right) + \sum_{\beta=0}^{N^2} c_{\beta}^{2} \Im \left( \left( \frac{z}{2} \right)^{\beta} \right). \quad (12) $$
In this formula, $\Re((\frac{z}{2})^{\beta})$, $\beta = 0, \dots, N^2$, are known harmonic polynomials and $\Re(\frac{d_{\alpha}}{z-z_{\alpha}})$ are known harmonic rational functions associated with known poles $z_{\alpha} = 1 + 2 \exp(-4(\sqrt{N_1} - \sqrt{\alpha}))$, for $\alpha = 1, \dots, N^1$. The sought coefficients $c_{\alpha}^1$, $\alpha = 1, \dots, N^1$, and $c_{\beta}^1$, $\beta = 1, \dots, N^2$, are computed by solving a linear least square problem to minimize the difference between $\Phi$ and $\tilde{\Phi}$ on a given set of points on $\partial\Omega_{\Phi}$. The function $\Phi$ is then used to compute three auxiliary functions $\Phi_{j,E}^{j-1}$, $\Phi_{j,E}^j$, and $\Phi_{j,E}^{j+1}$ through three different affine maps. Such maps map the point $(1, 0)$ to the $(j-1)$, $j$, or $(j-1)$-th vertex of $E$, and ensure that $E$ is included in the image of $\Omega_{\Phi}$. We note that, since $\Phi$ does not depend on $E$ or $\hat{E}$, it is computed only once.

Finally, for any element $E \in T_h$ and each vertex $v_j$ of $E$, with $= 1, \dots, N_E^v$, we can define the space $\mathcal{H}_{j,E}^{NN}$, as:

$$ \mathcal{H}_{j,E}^{NN} = \operatorname{span} \left\{ \{\tilde{p}_{\beta}\}_{\beta=1}^{2\ell^{NN}+1}, \Phi_{j,E}^{j-1}, \Phi_{j,E}^{j}, \Phi_{j,E}^{j+1} \right\}, \quad (13) $$

where $\{\tilde{p}_{\beta}\}_{\beta=1}^{2\ell^{NN}+1}$ are the harmonic polynomials in (11). Therefore, the dimension of $\mathcal{H}_{j,E}^{NN}$ is

$$ \dim \mathcal{H}_{j,E}^{NN} = 2\ell^{NN} + 4. \quad (14) $$

We observe that this dimension does not depend on the number of functions $N^1 + N^2$ used to approximate $\Phi$.

## 3.2 The Neural Network Architecture

As previously stated, the main ingredient of the NAVEM method is the approximation $\varphi_{j,E}^{NN}$ of the local VEM basis function $\varphi_{j,E}$. This approximation is a suitable function in $\mathcal{H}_{j,E}^{NN}$ which is built through a neural network trained by mimicking the function $\varphi_{j,E}$ at the boundary of $E$, where such function is well-known. For the ease of notation, we define

$$ \mathcal{H}_{j,E}^{NN} = \operatorname{span}\{h_k\}_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}}, \quad (15) $$

where $\{h_k\}_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}}$ coincides with the set of basis functions in (13). Therefore, the function $\varphi_{j,E}^{NN}$ can be expressed as

$$ \varphi_{j,E}^{NN} = \sum_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}} c_k^{\varphi}(v_j, E)h_k \quad (16) $$

for suitable coefficients $c_k^{\varphi}(v_j, E)$ depending on the pair $(v_j, E)$. In this context, the involved neural network is a parametric function mapping the pair $(v_j, E)$ to the vector $[c_k^{\varphi}(v_j, E)]_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}}$. We remark that, since a neural network accepts only vectors as inputs, the pair $(v_j, E)$ has to be suitably encoded. Multiple input encoding strategies are admissible. In this work, we adopt the approach described in [22], where an input encoding strategy is detailed, along with variability and input reduction strategies for generic polygons.

Let $x_0$ be the encoding of a pair $(v_j, E)$. In the following, we drop the subscripts $j, E$ to lighten the notation, while also admitting that the neural network can be evaluated on multiple pairs at once. The operator mapping $x_0$ to the coefficients $\mathbf{c}^{\varphi}(x_0) := [c_k^{\varphi}]_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}}$ is expressed as:

$$ \begin{align} \mathbf{x}_0^{\varphi} &:= x_0, \nonumber \\ \mathbf{x}_{\ell}^{\varphi} &:= \rho(A_{\ell}^{\varphi} \mathbf{x}_{\ell-1}^{\varphi} + b_{\ell}^{\varphi}), \tag{17} \\ \mathbf{c}^{\varphi}(x_0) &:= A_L^{\varphi} \mathbf{x}_{L-1}^{\varphi} + b_L^{\varphi}. \nonumber \end{align} $$

Equation (17) describes the so-called architecture of the neural network [34]. Here, $L$ denotes the number of layers, $\rho$ is a non-linear activation function, and $A_{\ell}^{\varphi} \in \mathbb{R}^{N_{\ell} \times N_{\ell-1}}$ and $b_{\ell}^{\varphi} \in \mathbb{R}^{N_{\ell}}$ are matrices and vectors containing the trainable coefficients of the neural network, which are optimized by minimizing a suitable loss function. Specifically, to optimize the coefficients of this neural network, we propose to minimize the distance between the traces of the functions $\varphi_{j,E}^{NN}$ and $\varphi_{j,E}$ on $\partial E$, i.e. we want to minimize the following quantity

$$ \epsilon_{j,E}^{\varphi} = \|\varphi_{j,E}^{NN} - \varphi_{j,E}\|_{H^{1/2}(\partial E)}, \quad (18) $$

for each pair $(j, E)$ in a given training dataset. We remark that this quantity is computable because it involves only the evaluations of $\varphi_{j,E}$ on the boundary of $E$, where this function is known in closed form.
Let $q_{j,E} := \nabla \varphi_{j,E}$ be the gradient of $\varphi_{j,E}$. We observe that his gradient can be approximated by the gradient of the function $\varphi_{j,E}^{NN}$ defined in (16) [21], which can be exactly computed as:

$$ \nabla \varphi_{j,E}^{NN} = \sum_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}} c_k^{\varphi}(v_j, E) \nabla h_k. \quad (19) $$

However, such a strategy is observed to be suboptimal. For this reason, it is convenient to introduce a second neural network to approximate $q_{j,E}$ with vector-functions contained in $\nabla \mathcal{H}_{j,E}^{NN}$. This neural network, similarly to (17), can be expressed as:

$$ \begin{aligned} x_0^q &:= x_0, \\ x_\ell^q &:= \rho(A_\ell^q x_{\ell-1}^q + b_\ell^q), & \ell = 1, \dots, L-1, \\ c^q(x_0) &:= A_L^q x_{L-1}^q + b_L^q. \end{aligned} \quad (20) $$

The output $c^q(x_0) := [c_k^q]_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}}$ of (20) represents the set of coefficients that allows us to approximate $q_{j,E}$ through the following vector

$$ q_{j,E}^{NN} = \sum_{k=1}^{\dim \mathcal{H}_{j,E}^{NN}} c_k^q(j, E) \nabla h_k. \quad (21) $$

The trainable coefficients $A_\ell^q$ and $b_\ell^q$ are optimized by minimizing a *computable* loss function. This cost term is defined through the following quantity

$$ \epsilon_{j,E}^q = \|(q_{j,E}^{NN} - q_{j,E}) \cdot t\|_{L^2(\partial E)}, \quad (22) $$

where $t$ denotes the tangential vector to the boundary $\partial E$. In [22], the authors proved that small values of $\epsilon_{j,E}^\varphi$ and $\epsilon_{j,E}^q$ are enough to ensure a good approximation of $\varphi_{j,E}$ and of $q_{j,E}$ on the entire element $E$, even if they are minimized only on $\partial E$. We observe that such an approach introduces a consistency error because (21) does not exactly coincide with the derivative of (16), but it significantly reduces the oscillations in the gradients of the NAVEM basis functions, leading to a more stable method.

## 3.3 The Training Phase

We highlight that employing the input encoding strategy presented in [22] requires defining a different network for each class of polygons characterized by a different number of vertices $N^v \ge 4$. For $N^v = 3$, the virtual element basis functions, as well as the NAVEM basis functions, coincide with the ones of the standard FEM. Thus, each training set has to contain only the class of polygons corresponding to the network that is trained.

For each class of polygons, let $\mathcal{T}^{\text{train}}$ be a training dataset of size $\#\mathcal{T}^{\text{train}}$ containing polygons $E$ with $N^v$ vertices. Then, to train the neural networks (17) and (20), i.e. to approximate the basis functions and their gradients, we minimize the following loss functions

$$ \begin{aligned} \mathcal{L}_{\varphi} &= \sqrt{\frac{1}{\#\mathcal{T}^{\text{train}} N^v} \sum_{E \in \mathcal{T}^{\text{train}}} \sum_{j=1}^{N^v} \left( \epsilon_{j,E}^{\varphi} \right)^2}, \\ \mathcal{L}_{q} &= \sqrt{\frac{1}{\#\mathcal{T}^{\text{train}} N^v} \sum_{E \in \mathcal{T}^{\text{train}}} \sum_{j=1}^{N^v} \left( \epsilon_{j,E}^{q} \right)^2}, \end{aligned} \quad (23) $$

respectively. Finally, a standard regularization term penalizing the $L^2$-norm of the trainable weights, with regularization coefficients $10^{-8}$, is added to both losses. We observe that these formulas differ from the ones defined in [22]. Indeed, in [22], the used loss functions coincide with the squared losses $\mathcal{L}_{\varphi}^2$ and $\mathcal{L}_{q}^2$. Note that the global minima of these two pairs of losses are the same, but we empirically observed better performance minimizing (23). Moreover, we here decide to employ a new optimizer: after performing 5000 epochs with the first-order ADAM optimizer [35] as in [22], we carry out 5000 updates with the self-scaled BFGS optimizer described in [36, 37]. This last advanced optimization technique enhances convergence by adaptively scaling the weights updates provided by the BFGS optimizer.

To appreciate the effectiveness of these new techniques, we compare the final values of the training losses by training the neural networks as described here and in the previous manuscript [22]. More precisely, we here build new neural networks characterized by the same architecture used in [22], i.e. they all have 5 layers, with 50 neurons in each layer, and the *tanh* activation function in each hidden layer. For the network associated with the class of polygons with $N^v$ vertices, the input dimension is $2(N^v - 1)$, while the output is always $2\ell^{NN} + 4$, as discussed in Section 3.1. We fix
$ℓ^{NN} = 20$, independently of the associated class of polygons, and set $R^{NN} = 3$. A fine-tuning strategy of all these values should be advisable, but it is beyond the scope of the current paper. Moreover, in order to train such networks, we use the same datasets RDQM, for $N^v = 4$, and VM, for $N^v \ge 5$, presented in [22]: the former contains randomly generated quadrilaterals, and the latter contains polygons sampled from Voronoi meshes generated through [38]. The entire code is written in Python, and the implementation of the definition, optimization, and evaluation of the neural networks heavily relies on the TensorFlow library [25].

The differences between the accuracy of the two approaches can be appreciated in Table 1. In particular, we report the value of $\mathcal{L}_{\varphi}$ and $\mathcal{L}_{q}$ for the networks used in [22] and for the networks used in the current manuscript. It is evident that the change of the loss and of the optimizer allows one to more accurately approximate both the basis functions and their gradients. We remark that these changes do not influence the overall computational cost of the training phase because the cost of the additional operations is negligible.

<table><caption>Table 1: Final training loss functions computed using the neural networks in [22] and the ones in the current manuscript.</caption><thead><tr><th>N<sup>v</sup></th><th>LC<sub>φ</sub> [22]</th><th>LC<sub>φ</sub></th><th>LC<sub>q</sub> [22]</th><th>LC<sub>q</sub></th></tr></thead><tbody><tr><td>4</td><td>4.62e-03</td><td>5.13e-04</td><td>1.00e-02</td><td>2.94e-03</td></tr><tr><td>5</td><td>1.97e-03</td><td>2.81e-04</td><td>4.97e-03</td><td>1.84e-03</td></tr><tr><td>6</td><td>1.02e-03</td><td>1.12e-04</td><td>2.96e-03</td><td>1.07e-03</td></tr><tr><td>7</td><td>2.15e-03</td><td>3.40e-04</td><td>5.22e-03</td><td>2.07e-03</td></tr><tr><td>8</td><td>1.21e-03</td><td>3.26e-04</td><td>4.88e-03</td><td>2.47e-03</td></tr></tbody></table>

We would like to remark that the comparison is meaningful because we are aiming to approximate the same scalar virtual element space $V_{h,1}(E)$.

## 3.4 The Online Phase and the NAVEM Discretization

In this section, assuming that suitable trained neural networks are available to approximate the basis functions related to $T_h$, we present the NAVEM discretization of the elastic problem (1). The main steps of the entire procedure can be subdivided in two main parts. The first one (steps S.1–S.3) is specific to the NAVEM approach and it is summarized in Figure 1. The second one (steps S.4–S.5) depends only on the problem and can be handled as in the standard Finite Element Method.

<figure><img src="image_6.png" /><figcaption>Figure 1: A sketch of the NAVEM online phase (steps S.1–S.3).</figcaption></figure>

Let us describe the first part. For each element $E$ in the tessellation $T_h$ of the domain $\Omega$, and for each vertex $v_j$ of $E$, with $j = 1, \dots, N_E^v$, we

S.1 classify the pair $(v_j, E)$, as discussed in Section 3.3, to select the neural networks that have to be used to approximate the corresponding basis function and its gradient;

S.2 encode the information $(v_j, E)$ to define a vector $x_0^{j,E}$ that represents the input vector for the neural networks;

S.3 predict the coefficients $c^{\varphi}(v_j, E)$ and $c^q(v_j, E)$ to compute $\varphi_{j,E}^{NN}$ and $q_{j,E}^{NN}$ according to (16) and (21), respectively.

Now, we define the local scalar NAVEM space $V_{h,1}^{NN}(E)$ and the related gradient space $\nabla V_{h,1}^{NN}(E)$ as:

$$V_{h,1}^{NN}(E) = \operatorname{span}\{\varphi_{j,E}^{NN}, j = 1, \dots, N_E^v\},$$

$$\nabla V_{h,1}^{NN}(E) = \operatorname{span}\{q_{j,E}^{NN}, j = 1, \dots, N_E^v\},$$
which represent the approximation spaces for the scalar VEM space $V_{h,1}(E)$, defined in 5, and the related gradient space $\nabla V_{h,1}(E)$, respectively.

The local NAVEM space for the displacement and its gradient can be defined following a standard strategy, i.e. it can be defined as

$$V_{h,1}^{NN}(E) := V_{h,1}^{NN}(E) \times V_{h,1}^{NN}(E), \quad \nabla V_{h,1}^{NN}(E) = \nabla V_{h,1}^{NN}(E) \times \nabla V_{h,1}^{NN}(E). \quad (24)$$

for each element $E \in \mathcal{T}_h$. By gluing together these local spaces, we obtain the global NAVEM discrete spaces for the displacement. Finally, we observe that, even if the NAVEM functions are not exactly continuous across elements, we do not double the degrees of freedom. Thus,

$$N^{\text{dof}} = \dim(\mathbf{V}_{h,1}^{NN}) = \dim(\mathbf{V}_{h,1}), \quad (25)$$

i.e. the dimension of $V_{h,1}^{NN}$ is twice the number of vertices of $\mathcal{T}_h$ that do not belong to the Dirichlet boundary.

Specifically, each NAVEM function $v_h^{NN} \in V_{h,1}^{NN}$ and its gradient can be written as:

$$v_h^{NN} = \sum_{i=1}^{N^{\text{dof}}/2} v_i \begin{bmatrix} \varphi_i^{NN} \\ 0 \end{bmatrix} + \sum_{i=N^{\text{dof}}/2}^{N^{\text{dof}}} v_i \begin{bmatrix} 0 \\ \varphi_{i-N^{\text{dof}}/2}^{NN} \end{bmatrix},$$

$$\nabla v_h^{NN} = \sum_{i=1}^{N^{\text{dof}}/2} v_i \begin{bmatrix} q_{1,i}^{NN} & q_{2,i}^{NN} \\ 0 & 0 \end{bmatrix} + \sum_{i=N^{\text{dof}}/2}^{N^{\text{dof}}} v_i \begin{bmatrix} q_{1,i-N^{\text{dof}}/2}^{NN} & 0 \\ 0 & q_{2,i-N^{\text{dof}}/2}^{NN} \end{bmatrix},$$

where $\{v_i\}_{i=1}^{N^{\text{dof}}}^N$ represents the set of degrees of freedom for the vector-function $v_h^{NN}$.

Given the output of the online phase S.1–S.3, the NAVEM implementation follows the exact steps performed to solved the elastic problem (3) with a standard Finite Element Method. In particular, the NAVEM discretization for the elastic problem (3) reads as: Find $u_h^{NN} \in V_{h,1}^{NN}$ such that

$$\sum_{E \in \mathcal{T}_h} (\sigma(\mathbf{x}, \nabla u_h^{NN}), \nabla v_h^{NN})_E = \sum_{E \in \mathcal{T}_h} (f, v_h^{NN})_E \quad \forall v_h^{NN} \in V_{h,1}^{NN}.$$

This problem is a non-linear problem that can be solved using a Newton-Raphson solver, eventually paired with an incremental force strategy when large displacements are taken into account. Thus, given an initial guess $u^{NN,0} \in V_{h,1}^{NN}$, for each non-linear step $k \ge 0$, we must solve the following problem: Find $\delta_h^{NN,k} \in V_{h,1}$ such that

$$\sum_{E \in \mathcal{T}_h} \left[ G_{NN}^{E}(u_h^{NN,k}, v_h^{NN}) + D G_{NN}^{E}(u_h^{NN,k}, v_h^{NN}) \cdot [\delta_h^{NN,k}] \right] = 0 \quad \forall v_h \in V_{h,1},$$

where

$$G_{NN}^{E}(u_h^{NN,k}, v_h^{NN}) = \int_{E} \sigma(\mathbf{x}, \nabla u_h^{NN,k}) : v_h^{NN} - \int_{E} f \cdot v_h^{NN}, \quad (26)$$

$$D G_{NN}^{E}(u_h^{NN,k}, v_h^{NN}) \cdot [\delta_h^{NN,k}] = \int_{E} (\mathcal{A}(\mathbf{x}, \nabla u_h^{NN,k}) : \nabla \delta_h^{NN,k}) : \nabla v_h^{NN}, \quad (27)$$

and $\mathcal{A}$ is the elastic modulus defined in (10). In the following, we will denote with the same bold symbol both the vector-function $v_h^{NN} \in V_{h,1}^{NN}$ and the vector of the related degrees of freedom, assuming the meaning is clear from the context. In conclusion, the second part of the NAVEM discretization can be summarized as:

**S.4 assemble** the tangent stiffness matrix $K^{NN,k}$ corresponding to (27) and the right-hand side $\mathbf{r}^{NN,k}$ related to (26).

**S.5 solve** the system

$$K^{NN,k} \delta^{NN,k} = -\mathbf{r}^{NN,k}, \quad (28)$$

and set $u^{NN,k+1} = \delta^{NN,k} + u^{NN,k}$.

The process is repeated until convergence. The exit vector $u^{NN,k+1}$ collects the coefficients $\{u_i\}_{i=1}^{N^{\text{dof}}}^N$ that describe the NAVEM solution $u_h^{NN}$.
# 4 Numerical Results

In this section, we test the discussed method on linear and non-linear elasticity problems and on different meshes.

We highlight that in all the experiments, the NAVEM basis functions are predicted by exploiting the same neural networks built to perform the experiment in Section 3.3.

Denoting by $u$ the exact displacement, we test the performance of the NAVEM by looking at the behaviour of the following errors

$$ \text{err}_{0}^{NN} = \sqrt{\sum_{E \in \mathcal{T}_h} \|u - u_h^{NN}\|_{L^2(E)}^2}, \qquad (29) $$

$$ \text{err}_{1}^{NN} = \sqrt{\sum_{E \in \mathcal{T}_h} \|\nabla u - \nabla u_h^{NN}\|_{L^2(E)}^2}, $$

as the mesh parameter $h$ decreases. For comparison purposes, we solve these problems also with the standard virtual element method. In this case, since we can not access the point-wise evaluation of virtual functions, we define the VEM errors as usual, that is [8]

$$ \text{err}_{0}^{VEM} = \sqrt{\sum_{E \in \mathcal{T}_h} \|u - \Pi_{1}^{E,0} u_{h}^{VEM}\|_{L^2(E)}^2}, \qquad (30) $$

$$ \text{err}_{1}^{VEM} = \sqrt{\sum_{E \in \mathcal{T}_h} \|\nabla u - \Pi_{0}^{E,0} \nabla u_{h}^{VEM}\|_{L^2(E)}^2}, $$

where we denote by $u_h^{VEM}$ the virtual element solution. We remark that to evaluate the high-order projection $\Pi_{1}^{E,0} u_h^{VEM}$ in virtual element, an enhanced version of the local space should be taken into account [39].

## 4.1 Test 1: The linear case - convergence test

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 2: Test 1: Left: Last refinement of the Distorted-Square family of meshes. Right: Behaviour of errors (29)–(30) as the mesh parameter $h$ decreases.</figcaption></figure>

We consider the square domain $\Omega = (0, 1)^2$ and the linear version of the problem (3), where

$$ a(u, v) = 2\mu \int_{\Omega} \epsilon(u) : \epsilon(v) + \lambda \int_{\Omega} \operatorname{div} u \operatorname{div} v, \qquad (31) $$

and $\mu, \lambda$ are the positive Lamé coefficients, whereas $\epsilon(u) = \frac{1}{2}(\nabla u + (\nabla u)^T)$. We highlight that the algorithm detailed in the previous sections is valid for any constitutive law and, therefore, it offers a flexible setting to deal with elastic problems [10].
In this test, we take the material constant values $\mu = 1.5$, $\lambda = 3$ and we set the boundary conditions and the body load $f$ in such a way the exact displacement is:

$$ \mathbf{u}(\mathbf{x}) = \begin{bmatrix} u_1(\mathbf{x}) \\ 5u_1(\mathbf{x}) \end{bmatrix}, \quad u_1(\mathbf{x}) = 16x_1(1-x_1)x_2(1-x_2) + 1.1. \qquad (32) $$

In particular, we choose to clamp the displacement at $\Gamma_D = \{\mathbf{x} \in \Gamma : x_2 = 0\} \cup \{\mathbf{x} \in \Gamma : x_1 = 0\}$, whereas the remaining boundary is free.

We build a family of four distorted quadrilateral meshes obtained by randomly distorting square grids of $4 \times 4$, $8 \times 8$, $16 \times 16$, and $32 \times 32$ elements. The last refinement is shown in Figure 2. In the Virtual Element Method, the stabilizing parameter is set to be $2\mu$.

The behavior of errors (29)-(30) are shown in Figure 2. We observe that the expected order of convergence for both methods is achieved, i.e. the errors decrease as $O(h^2)$ for the $L^2$-error and as $O(h)$ for the error in the $H^1$-seminorm. Moreover, coherently with the results obtained for the scalar version of the NAVEM method in [22], we observe lower error constants with respect to VEM for both errors. We also note that the elements in this family of meshes are well represented by the training set. Thus, the network accuracy does not influence these results.

Finally, we want to remark that the evaluation of a neural network has an intrinsic cost that tends to be negligible for fine meshes, where the cardinality of the input of the neural networks grows, as shown in [22]. In a vectorial setting, as the one offered by the elasticity problems, this intrinsic cost is even lower because the neural networks are only used to construct the basis functions of $V_{h,1}^{NN}(E)$. When a function in $V_{h,1}^{NN}(E)$ has to be evaluated, the functions in $V_{h,1}^{NN}(E)$ associated with each dimension are reused multiple times. This is possible because the vector space is just the Cartesian product of a scalar space, approximated with the neural network, with itself.

## 4.2 Test 2: A benchmark non-linear elastic problem with analytical solution

In this section, we propose a test taken from [10] to show some advantages of using an approach that does not require any stabilization term, such as the NAVEM. Specifically, we consider $\Omega = (0, 1)^2$ and set $\Gamma_D = \Gamma$. We select the following constitutive law

$$ \sigma(\mathbf{x}, \nabla \mathbf{u}(\mathbf{x})) = \mu(\epsilon(\mathbf{u}))\epsilon(\mathbf{u}) \qquad (33) $$

where the scalar function $\mu$ is given by

$$ \mu(\epsilon(\mathbf{u})) = 3(1 + \|\epsilon(\mathbf{u})\|^2) = 3(1 + \sum_{n,m=1}^{2} [\epsilon(\mathbf{u})]_{nm}^2)^2. \qquad (34) $$

We further consider two external body forces, compatible with the following two analytical solutions:

$$ \textbf{Case 1: } \mathbf{u}(\mathbf{x}) = x_1(1-x_1)x_2(1-x_2) \begin{bmatrix} 1 \\ 1 \end{bmatrix}; $$

$$ \textbf{Case 2: } \mathbf{u}(\mathbf{x}) = 80x_1(1-x_1)x_2(1-x_2) \begin{bmatrix} 1 \\ 1 \end{bmatrix}. $$

In particular, we notice that in **Case 1** the solution is characterized by deformations of moderate magnitude, whereas much larger deformations occur in **Case 2**. For the current experiment, we consider two families of meshes. The former comprises four square meshes of $4 \times 4$, $8 \times 8$, $16 \times 16$, and $32 \times 32$ square elements each. The second family is made up of 4 regular Voronoi meshes subjected to sine distortions realized with the MATLAB library mVEM [38]. The last refinement for each family of meshes is shown in Figure 3. Preliminary, we remark that the neural networks used to predict the NAVEM basis functions are the ones used to make the experiment reported in Section 3.3.

Due to non-linearity and possible large deformations, we here decide to test the NAVEM and VEM performance with and without the incremental force method described in Section 2.2, along with the Newton-Raphson scheme. Specifically, we set $N = 20$ for the incremental procedure and, for each $n = 1, \dots, N-1$, we stop the Newton-Raphson scheme after 20 non-linear iteration steps. At the incremental step $n = N$, we allow the Newton-Raphson scheme to perform all the non-linear iteration steps needed to reach the desired accuracy, i.e. we require that the relative norm of the residual is below $10^{-10}$.

We observe that, when the incremental method is applied, the VEM stabilization is chosen as the norm-based stabilization described in Section 2.2, where $w_h = u_h^{n-1}$. On the other hand, when we apply the plain Newton-Raphson method,
<figure><img src="image_9.png" /><figcaption>Figure 3: Test 2: Last refinements related to the two families of meshes.</figcaption></figure>

without the incremental force approach, we set $w_h = u_h^0 = 0$. In particular, we observe that this last choice coincides with the “fixed scaling” described in [10].

Figure 4 reports the behaviour of the $H^1$-errors in (29) and (30), as the mesh parameter $h$ decreases, for the test **Case 1**. Concerning this test case, we observe that the expected order of convergence is attained for both the methods with and without employing the incremental strategy. Specifically, we note that the choice of stabilization does not seem to influence the accuracy of the VEM method for both the families of meshes. Moreover, we observe that the NAVEM error curves are downward shifted with respect to the VEM error curves, as usual. Moreover, we highlight that the $L^2$-error curves follow a very similar behaviour. Thus we decide to not report this verbose information.

<figure><img src="image_10.png" /><figcaption>Figure 4: Test 2, **Case 1**: Behaviour of the $H^1$-errors (29)–(30) as the mesh parameter $h$ decreases. “(U)” means that the incremental force method is adopted. Left: Cartesian grids family. Right: Distorted-Voronoi grids family.</figcaption></figure>

Figure 4 reports the behaviour of $H^1$-error (29)–(30) as the mesh parameter $h$ decreases for the test **Case 2**, where larger displacements are recorded. Unlike the previous test case, a good choice of the stabilization term is essential to obtain good quality results when using the VEM method. Coherently with the results shown in [10], we observe that adopting the incremental strategy is mandatory for VEM to avoid to level-off the stabilization term with respect to the consistency term, losing accuracy, when Voronoi meshes are employed. On the other hand, we observe that this behaviour is less stressed in the presence of quadrilateral elements, where the stabilization term is generally less pronounced. Concerning the NAVEM method, we observe that it does not need at all the incremental strategy, even in the presence of the larger displacement related to **Case 2**. Indeed, it behaves in the right way independently of the family of meshes and of the magnitude of the displacement for this test, while again showing a lower error constants with respect to VEM.

Moreover, in Figure 6, we observe that NAVEM does not suffer of the same unstable spurious oscillating behaviour, remarked also in [10], that characterizes the VEM solution when a fixed scaling is employed for the **Case 2**. In this
<figure><img src="image_11.png" /><figcaption>Figure 5: Test 2, **Case 2**: Behaviour of the $H^1$-errors (29)–(30) as the mesh parameter $h$ decreases. “(U)” means that the incremental force method is adopted. Left: Cartesian grids family. Right: Distorted-Voronoi grids family.</figcaption></figure>

figure, we draw a slice of the first component of the VEM solution, of the NAVEM solution, and of the exact solution. It is evident that the VEM solution is characterized by large oscillations that ruin its accuracy, whereas the other two curves are almost exactly overlapped and are indistinguishable in the plot.

<figure><img src="image_12.png" /><figcaption>Figure 6: Test 2: Slice of the component $u_1$ related to **Case 2** and the last refinement of the Distorted-Voronoi family. Red: VEM. Blue: NAVEM. Black: exact displacement. We observe that the blue curve overlaps the black one.</figcaption></figure>

Finally, we observe that the robustness that characterizes NAVEM translates to a higher computational efficiency with respect to VEM. Indeed, when the incremental strategy is not employed, the Newton-Raphson method requires 125 non-linear steps to reach the desired accuracy for the last refinement of Distorted-Voronoi family and **Case 2**. Whereas, since VEM requires the incremental strategy to perform a correct choice for the stabilization, the whole non-linear incremental strategy needs about 500 steps to converge.

## 4.3 Test 3: Neo-Hookean Hyperelasticity

In this last test case, we consider Problem (1) with $\sigma$ corresponding to the non-symmetric first Piola-Kirchhoff stress tensor $\mathbf{P}$. More precisely, we focus on a hyperelastic material of Neo-Hookean type by setting

$$ \mathbf{P}(\mathbf{F}) = \mu(\mathbf{F} - \mathbf{F}^{-T}) + \lambda J \Theta(J) \Theta'(J) \mathbf{F}^{-T}, \quad (35) $$
where $\mathbf{F}$ is the deformation gradient $\mathbf{F} = \mathbf{I} + \nabla \mathbf{u}$, $J \in \mathbb{R}_{>0}$ is the determinant of $\mathbf{F}$, and $\Theta: \mathbb{R}_{>0} \to \mathbb{R}$ is a smooth function such that $\Theta(J) = 0 \Leftrightarrow J = 1$ and $\Theta'(1) \neq 0$ [40]. In equations (35), $\lambda$ and $\mu$ are given positive constants.

We observe that, to carry out the discretization related to a Neo-Hookean material, a discrete version of the determinant $J$ must be provided. In the Virtual Element framework, two possible discretization are considered. The simplest choice coincides with taking on each element $E \in \mathcal{T}_h$, for each $k \ge 0$ and $n \ge 1$:

$$J = \det(\mathbf{I} + \nabla \mathbf{u}_h^{\text{VEM},n,k}) \approx \det(\mathbf{I} + \Pi_0^{0,E} \nabla \mathbf{u}_h^{\text{VEM},n,k}).$$

Again, since the virtual functions are not known in a closed form in the interior of the element, a projector operator should be introduced to access the point-wise evaluations of such functions. With this kind of approximation, the determinant $J$ is constant on each element $E \in \mathcal{T}_h$. A more stable version is proposed in [14], where $J$ is approximated with its mean value over $E \in \mathcal{T}_h$. More precisely, the authors in [14] approximate $J$ as follows:

$$J \approx \frac{1}{|E|} \int_E J = \int_{\tilde{E}} d\tilde{\mathbf{x}} = \frac{|\tilde{E}|}{|E|}, \qquad (36)$$

where $|E|$ denotes the measure of the element $E \in \mathcal{T}_h$, whereas $\tilde{\cdot}$ defines a quantity $\cdot$ in the current deformed configuration. In the NAVEM framework, we evaluate the determinant simply as

$$J = \det(\mathbf{I} + \nabla \mathbf{u}_h^{\mathcal{NN},n,k}), \qquad (37)$$

since we are able to evaluate NAVEM functions and their gradients everywhere.

We test the method considering

$$\Theta(J) = \log J, \quad \Theta'(J) = \frac{1}{J}, \quad \lambda = 5.1, \quad \mu = 1.0,$$

whose values of Lamé coefficients correspond to a Poisson ratio $\approx 0.45$. We observe that in this case, the elastic modulus $\mathbb{A}$ in (10) in this specific case depends only on the deformation gradient and it is given by [40]:

$$\mathbb{A}(\mathbf{F}) = \frac{\partial \mathbf{P}(\mathbf{F})}{\partial \mathbf{F}} = \mu(\mathbf{I} \otimes \mathbf{I} + \mathbf{F}^{-T} \otimes \mathbf{F}^{-1}) - \lambda J \Theta(J) \Theta'(J) \mathbf{F}^{-T} \otimes \mathbf{F}^{-1} \\ + \lambda J [\Theta(J) \Theta'(J) + J(\Theta'(J))^2 + J \Theta(J) \Theta''(J)] \mathbf{F}^{-T} \otimes \mathbf{F}^{-T},$$

where $\otimes$, $\otimes$ and $\overline{\otimes}$ are defined such that $\{\mathbf{A} \otimes \mathbf{B}\}_{ijkl} = \mathbf{A}_{ij}\mathbf{B}_{kl}$, $\{\mathbf{A} \otimes \mathbf{B}\}_{ijkl} = \mathbf{A}_{il}\mathbf{B}_{jk}$ and $\{\mathbf{A} \otimes \mathbf{B}\}_{ijkl} = \mathbf{A}_{ik}\mathbf{B}_{jl}$, for all $i, j, k, l = 1, \dots, 2$ and for each pair of second-order tensors $\mathbf{A}$ and $\mathbf{B}$.

The initial configuration $\Omega$ is set as the unit square $[0, 1]^2$ and the displacement is clamped at $\Gamma_D = \{0\} \times [0, 1]$, while the remaining part of the boundary is free. The body load is given by

$$f(x_1, x_2) = [100x_2^3, 0] \quad \forall (x_1, x_2) \in \Omega. \qquad (38)$$

We simulate this test on the mesh shown in Figure 7a, generated as in [41], performing $N = 50$ iterations of the incremental force method. At each step $n = 1, \dots, N$, we iterate the Newton-Raphson method until the residual norm is below $1.0e-09$. Figures 7b and 7c show the $x_2$-component of the discrete displacement obtained with NAVEM and VEM, respectively, whereas Figure 8b shows the result of the Finite Element discretization of the problem on the finer triangular mesh 8a. We first observe that, even though we employ the more stable approximation of determinant, defined in (36), to produce virtual element results, we do not appreciate a significant difference between the two approximation procedures for the determinant in this test case. Moreover, we highlight that all the aforementioned methods take about 5 Newton steps for each incremental step $n$ to reach the desired accuracy. Despite this, we can observe some spurious oscillations in the $x_2$-component of the virtual element displacement that are not present in the NAVEM or FEM solution.

To conclude, we observe that the NAVEM has a simpler formulation than the VEM one and, in the provided numerical tests, it is more stable in the presence of large deformations. On the other hand, despite having a formulation similar to the FEM one, it ensures greater flexibility by allowing the use of more general meshes.

# 5 Conclusions

In this paper, we present the Neural Approximated Virtual Element Method to solve linear and non-linear elasticity problems. The NAVEM is a polygonal method that approximates via neural networks the function spaces used in the
<figure><img src="image_13.png" /><figcaption>Figure 7: Test 3: Left: Mesh used in the simulation. Center: Displacement in the $x_2$-direction obtained with the NAVEM method. Right: Displacement in the $x_2$-direction obtained with the VEM method.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 8: Test 3: Left: Mesh used to compute the FEM reference solution. Right: Displacement in the $x_2$-direction obtained with the FEM method.</figcaption></figure>

Virtual Element Method. This neural approximation provides closed-form local basis functions, eliminating the need for additional stabilization and projection operators.

The absence of a stabilization term is particularly advantageous in nonlinear problems, where designing a suitable and computable stabilization operator that preserves accuracy can be challenging. Elasticity problems represent a particular domain in which numerous different non-linearities may be present, depending on the simulated materials. For this class of problems, we discuss the NAVEM formulation, emphasizing the simpler form with respect to the standard VEM formulation. Numerical results are in agreement with what has been observed for elliptic problems in [21, 22] and highlight the benefits of bypassing stabilization operators.

Future perspectives include extending NAVEM to more physically relevant problems, such as handling internal constraints like incompressibility. Additionally, crack propagation and adaptivity present promising applications for this method. Indeed, as observed in [22], the method can be easily optimized for mesh comprising triangles with hanging nodes. This feature opens the possibility of integrating NAVEM with standard finite element methods by refining only elements where non-linearities are most intense: introducing hanging nodes instead of globally remeshing neighboring elements.

## Acknowledgements

The author S.B. kindly acknowledges partial financial support provided by PRIN project "Advanced polyhedral discretisations of heterogeneous PDEs for multiphysics problems" (No. 20204LN5N5_003), by PNRR M4C2 project of CN00000013 National Centre for HPC, Big Data and Quantum Computing (HPC) (CUP: E13C22000990001) and the funding by the European Union through project Next Generation EU, M4C2, PRIN 2022 PNRR project
P2022BH5CB_001 "Polyhedral Galerkin methods for engineering applications to improve disaster risk forecast and management: stabilization-free operator-preserving methods and optimal stabilization methods.". The author M.P. kindly acknowledges financial support provided by PEPR/IA (https://www.pepr-ia.fr/). The author G.T. kindly acknowledges the financial support provided by INdAM-GNCS Project "Metodi numerici efficienti per problemi accoppiati in sistemi complessi" (CUP: E53C24001950001).

## References

[1] C. Canuto, D. Fassino, Higher-order adaptive virtual element methods with contraction properties, Mathematics in Engineering 5 (6) (2023) 1-33. doi:10.3934/mine.2023101.
URL https://www.aimspress.com/article/doi/10.3934/mine.2023101

[2] S. Berrone, S. Ferraris, D. Grappein, G. Teora, F. Vicini, A 3d-1d virtual element method for modeling root water uptake (2024). arXiv:2412.12884.
URL https://arxiv.org/abs/2412.12884

[3] S. Berrone, F. Vicini, Effective polygonal mesh generation and refinement for vem, Mathematics and Computers in Simulation 231 (2025) 239-258. doi:https://doi.org/10.1016/j.matcom.2024.12.007.
URL https://www.sciencedirect.com/science/article/pii/S0378475424004786

[4] S. E. Leon, D. W. Spring, G. H. Paulino, Reduction in mesh bias for dynamic fracture using adaptive splitting of polygonal finite elements, International Journal for Numerical Methods in Engineering 100 (8) (2014) 555-576.
arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.4744, doi:https://doi.org/10.1002/nme.4744.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.4744

[5] S. O. R. Biabanaki, A. R. Khoei, P. Wriggers, Polygonal finite element methods for contact-impact problems on non-conformal meshes, Computer Methods in Applied Mechanics and Engineering 269 (2014) 198-221.
URL https://apisemanticscholar.org/CorpusID:119486202

[6] A. L. Gain, G. H. Paulino, L. S. Duarte, I. F. Menezes, Topology optimization using polytopes, Computer Methods in Applied Mechanics and Engineering 293 (2015) 411-430. doi:https://doi.org/10.1016/j.cma.2015.05.007.
URL https://www.sciencedirect.com/science/article/pii/S0045782515001735

[7] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, A. Russo, Basic principles of Virtual Element Methods, Mathematical Models and Methods in Applied Sciences 23 (01) (2013) 199-214.

[8] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, Virtual element method for general second order elliptic problems on polygonal meshes, Mathematical Models and Methods in Applied Sciences 26 (04) (2016) 729-750.
arXiv:https://doi.org/10.1142/S0218202516500160, doi:10.1142/S0218202516500160.
URL https://doi.org/10.1142/S0218202516500160

[9] L. Beirão da Veiga, F. Brezzi, L. D. Marini, Virtual elements for linear elasticity problems, SIAM Journal on Numerical Analysis 51 (2) (2013) 794-812. arXiv:https://doi.org/10.1137/120874746, doi:10.1137/120874746.
URL https://doi.org/10.1137/120874746

[10] L. Beirão da Veiga, C. Lovadina, D. Mora, A virtual element method for elastic and inelastic problems on polytope meshes, Computer Methods in Applied Mechanics and Engineering 295 (2015) 327-346. doi:https://doi.org/10.1016/j.cma.2015.07.013.
URL https://www.sciencedirect.com/science/article/pii/S004578251500225X

[11] S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Vol. 15 of Texts in Applied Mathematics, Springer, 2008. doi:10.1007/978-0-387-75934-0.
URL http://dx.doi.org/10.1007/978-0-387-75934-0

[12] Ph. Ciarlet, The Finite Element Method for Elliptic Problems, Vol. 40 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2002. arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9780898719208, doi:10.1137/1.9780898719208.
URL https://epubs.siam.org/doi/abs/10.1137/1.9780898719208

[13] S. Berrone, S. Scialò, G. Teora, The mixed virtual element discretization for highly-anisotropic problems: the role of the boundary degrees of freedom, Mathematics in Engineering 5 (6) (2023) 1-32. doi:10.3934/mine.2023099.
URL https://www.aimspress.com/article/doi/10.3934/mine.2023099
[14] H. Chi, L. B. da Veiga, G. Paulino, Some basic formulations of the virtual element method (vem) for finite deformations, Computer Methods in Applied Mechanics and Engineering 318 (2017) 148–192. doi:https://doi.org/10.1016/j.cma.2016.12.020.
URL https://www.sciencedirect.com/science/article/pii/S0045782516309094

[15] A. Chen, N. Sukumar, Stabilization-free serendipity virtual element method for plane elasticity, Computer Methods in Applied Mechanics and Engineering 404 (2023) 115784. doi:https://doi.org/10.1016/j.cma.2022.115784.
URL https://www.sciencedirect.com/science/article/pii/S004578252200740X

[16] A. Chen, N. Sukumar, Stabilization-free virtual element method for plane elasticity, Computers & Mathematics with Applications 138 (2023) 88–105. doi:https://doi.org/10.1016/j.camwa.2023.03.002.
URL https://www.sciencedirect.com/science/article/pii/S0898122123000913

[17] F. Credali, S. Bertoluzza, D. Prada, Reduced basis stabilization and post-processing for the virtual element method, Computer Methods in Applied Mechanics and Engineering 420 (2024) 116693. doi:https://doi.org/10.1016/j.cma.2023.116693.
URL https://www.sciencedirect.com/science/article/pii/S0045782523008162

[18] M. Trezzi, U. Zerbinati, When rational functions meet virtual elements: the lightning virtual element method, Calcolo 61 (3) (2024) 35. doi:10.1007/s10092-024-00585-1.
URL https://doi.org/10.1007/s10092-024-00585-1

[19] J. S. Hesthaven, G. Rozza, B. Stamm, et al., Certified reduced basis methods for parametrized partial differential equations, Vol. 590, Springer, 2016.

[20] A. Gopal, L. N. Trefethen, Solving laplace problems with corner singularities via rational functions, SIAM Journal on Numerical Analysis 57 (5) (2019) 2074–2094. arXiv:https://doi.org/10.1137/19M125947X,
doi:10.1137/19M125947X.
URL https://doi.org/10.1137/19M125947X

[21] S. Berrone, D. Oberto, M. Pintore, G. Teora, The lowest-order neural approximated virtual element method, in: Numerical Mathematics and Advanced Applications ENUMATH 2023, Vol. 1, Springer, 2025, pp. 129–138.
URL https://doi.org/10.1007/978-3-031-86173-4_13

[22] S. Berrone, M. Pintore, G. Teora, The lowest-order neural approximated virtual element method on polygonal elements, Computers & Structures 314 (2025) 107753.
URL https://doi.org/10.1016/j.compstruc.2025.107753

[23] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (7553) (2015) 436–444.

[24] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific machine learning through physics-informed neural networks: Where we are and what’s next, Journal of Scientific Computing 92 (3) (2022) 88. doi:10.1007/s10915-022-01939-z.
URL https://doi.org/10.1007/s10915-022-01939-z

[25] M. Abadi, et al., TensorFlow: Large-scale machine learning on heterogeneous systems, software available from tensorflow.org (2015).
URL http://tensorflow.org/

[26] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, PyTorch: An imperative style, high-performance deep learning library, in: Proceedings of the 33rd International Conference on Neural Information Processing Systems, Curran Associates Inc., 2019, pp. 8024–8035.
URL https://dl.acm.org/doi/10.5555/3454287.3455008

[27] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, Q. Zhang, JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax (2018).

[28] L. Beirão da Veiga, F. Brezzi, L. D. Marini, A. Russo, The virtual element method, Acta Numerica 32 (2023) 123–202. doi:10.1017/S0962492922000095.

[29] N. Kim, Introduction to Nonlinear Finite Element Analysis, Springer US, 2014.
URL https://books.google.it/books?id=PDeSBQAAQBAJ

[30] P. Wriggers, Nonlinear Finite Element Methods, Springer Berlin Heidelberg, 2008.
URL https://books.google.it/books?id=IaV0wgG2jacC
[31] P. Wriggers, F. Aldakheel, B. Hudobivnik, Virtual Element Methods in Engineering Sciences, Springer International Publishing, 2023.
URL https://books.google.it/books?id=KdnfEAAAQBAJ

[32] S. Berrone, G. Teora, F. Vicini, Improving high-order vem stability on badly-shaped elements, Mathematics and Computers in Simulation 216 (2024) 367-385. doi:https://doi.org/10.1016/j.matcom.2023.10.003.
URL https://www.sciencedirect.com/science/article/pii/S0378475423004287

[33] J. B. Perot, C. Chartrand, A mimetic method for polygons, Journal of Computational Physics 424 (2021) 109853.
doi:https://doi.org/10.1016/j.jcp.2020.109853.
URL https://www.sciencedirect.com/science/article/pii/S0021999120306276

[34] I. Goodfellow, Y. Bengio, A. Courville, Y. Bengio, Deep learning, Vol. 1, MIT press Cambridge, 2016.

[35] D. Kingma, J. Ba, Adam: a method for stochastic optimization, International Conference on Learning Representations (2014). doi:https://doi.org/10.48550/arXiv.1412.6980.

[36] S. Wright, J. Nocedal, et al., Numerical Optimization, Vol. 35, Springer, 1999.

[37] E. Kiyani, K. Shukla, J. F. Urbán, J. Darbon, G. E. Karniadakis, Which optimizer works best for physics-informed neural networks and kolmogorov-arnold networks?, arXiv preprint arXiv:2501.16371 (2025).

[38] Y. Yu, mvem: A matlab software package for the virtual element methods (2022). arXiv:2204.01339.
URL https://arxiv.org/abs/2204.01339

[39] B. Ahmad, A. Alsaedi, F. Brezzi, L. Marini, A. Russo, Equivalent Projectors for Virtual Element Methods, Comput Math Appl 66 (3) (2013) 376-391.

[40] M. Abbas, A. Ern, N. Pignet, Hybrid high-order methods for finite deformations of hyperelastic materials, Computational Mechanics 62 (4) (2018) 909-928. doi:10.1007/s00466-018-1538-0.
URL https://doi.org/10.1007/s00466-018-1538-0

[41] A. Lamperti, M. Cremonesi, U. Perego, A. Russo, C. Lovadina, A hu-washizu variational approach to self-stabilized virtual elements: 2d linear elastostatics, Computational Mechanics 71 (2023) 935-955. doi:10.1007/s00466-023-02282-2.
URL https://doi.org/10.1007/s00466-023-02282-2