# Adaptive finite element interpolated neural networks

Santiago Badia1,2,∗<EMAIL_ADDRESS>,  Wei Li1<EMAIL_ADDRESS>and
Alberto F. Martín3<EMAIL_ADDRESS>1 School of Mathematics, Monash
University, Clayton, Victoria 3800, Australia.  2 Centre Internacional de
Mètodes Numèrics a l’Enginyeria, Campus Nord, UPC, 08034, Barcelona, Spain.  3
School of Computing, The Australian National University, Canberra ACT 2600,
Australia.  ∗ Corresponding author.

###### Abstract.

The use of neural networks to approximate partial differential equations
(PDEs) has gained significant attention in recent years. However, the
approximation of PDEs with localised phenomena, e.g., sharp gradients and
singularities, remains a challenge, due to ill-defined cost functions in terms
of pointwise residual sampling or poor numerical integration. In this work, we
introduce $h$-adaptive finite element interpolated neural networks. The method
relies on the interpolation of a neural network onto a finite element space
that is gradually adapted to the solution during the training process to
equidistribute a posteriori error indicator. The use of adaptive interpolation
is essential in preserving the non-linear approximation capabilities of the
neural networks to effectively tackle problems with localised features. The
training relies on a gradient-based optimisation of a loss function based on
the (dual) norm of the finite element residual of the interpolated neural
network. Automatic mesh adaptation (i.e., refinement and coarsening) is
performed based on a posteriori error indicators till a certain level of
accuracy is reached. The proposed methodology can be applied to indefinite and
nonsymmetric problems. We carry out a detailed numerical analysis of the
scheme and prove several a priori error estimates, depending on the
expressiveness of the neural network compared to the interpolation mesh. Our
numerical experiments confirm the effectiveness of the method in capturing
sharp gradients and singularities for forward PDE problems, both in 2D and 3D
scenarios. We also show that the proposed preconditioning strategy (i.e.,
using a dual residual norm of the residual as a cost function) enhances
training robustness and accelerates convergence.

###### Key words and phrases:

neural networks, PINNs, finite elements, PDE approximation, $h$-adaptivity

##  1\. Introduction

pde are an effective mathematical tool for modelling a variety of physical
phenomena in science and engineering, such as heat conduction, fluid dynamics
and electromagnetism. Since analytical solutions of partial differential
equations are only available in few special cases, numerical methods, such as
the finite element method (FEM), become the primary approach for approximating
PDE solutions using computers. The FEM has solid mathematical foundations [1]
and it can readily handle a broad range of PDE problems. Over the years,
advanced discretisation techniques have been proposed [2], and highly
efficient (non)linear solvers that can exploit large-scale supercomputers have
been designed and polished [3, 4, 5].

The FEM is a linear method that provides the best approximation in some
specific measure on a finite-dimensional vector space. However, this space is
not adapted to local features (e.g., sharp gradients or singularities) and the
convergence is slow for problems that exhibit multiple scales. Adaptive FEs,
e.g., the $h$-adaptive FEM, adapt the space to the solution by gradually
increasing resolution in regions with larger _a posteriori_ error estimations.
It adds an external loop, known as the solve, estimate, mark and refine loop,
to the standard FEM workflow. Robust _a posteriori_ error estimators have been
designed for diverse problems (see e.g., [6, 7]), making the adaptive FEM a
widely employed technique in scientific disciplines.

Recently, deep learning, especially NNs, has demonstrated outstanding
proficiency in many challenging tasks, including biomedical image segmentation
[8], natural language processing [9] and speech recognition [10], just to
mention a few. The striking success achieved by deep learning has encouraged
researchers to explore the application of NNs in solving PDEs. In contrast to
standard discretisation methods, NNs are parametrised manifolds that enjoy
non-linear approximation power, i.e., NNs have built-in adaptation
capabilities essential in multi-scale problems. Notably, PINNs have been
proposed in [11] as a novel approach to tackle PDE-constrained problems. The
PINN method integrates the PDE into the loss function used to train the NN.
This loss evaluates the strong PDE residual on collocation points, and the
derivatives within this loss can be computed using automatic differentiation.
pinn are easy to implement and have been examined in different types of PDEs,
such as fractional [12] and stochastic PDEs [13].

Many recent works have tried to address the limitations of PINNs, aiming to
improve their robustness, efficiency and accuracy. The authors in [14]
introduce VPINNs, which make use of a polynomial test space and take the norm
of the corresponding residual vector as the loss function. To overcome the
integration challenge in VPINNs, the authors in [15] propose the interpolated
(IVPINN) method. This method interpolates a NN onto a FE space and uses this
polynomial interpolation in the variational residual. In PINNs, the imposition
of Dirichlet boundary conditions is a challenging task; as opposed to FE
spaces, NNs are not interpolatory objects. The addition of a penalty term to
weakly impose the Dirichlet boundary conditions has a very negative effect on
the training. In [15], the authors compose the NN with a polynomial function
to design a NN that vanishes on the boundary. This approach is applicable to
simple geometries but its extension to general geometries is unclear [16].
Besides, while better than the penalty method, such composition negatively
impacts the generalisation of NNs themselves [17].

Our previous work [17] introduces the  interpolated  (FEINN) method. This
method interpolates the NN onto a FE space that strongly satisfies the
boundary condition, avoiding the issues discussed above. Besides, we propose
to use the dual norm of the PDE residual as the loss function, which can be
re-stated as a preconditioned loss function that is well-posed at the
continuous level; see also [18] for the application of dual norms in VPINNs.
The preconditioning requires the application of the inverse of the inner
product matrix (e.g., a Laplacian in $H^{1}(\Omega)$-conforming cases), which
we efficiently approximate using geometric multigrid techniques. The FEINN
method has been successfully applied to, e.g., forward PDEs with smooth or
non-smooth solutions and inverse heat conduction design problems. We show that
the trained NNs can outperform the FEM solution by several orders of magnitude
when the exact solution is smooth.

Despite these efforts, the training of PINNs and related methods often fail
when dealing with target solutions containing high-frequency modes [19], e.g.,
solutions with sharp gradients or discontinuities. Some works aim to improve
the expressivity of the NN to capture these solutions, adding more
neurons/layers [20], or adapting the non-linear activation functions [21].
However, none of these works solve the integration issue; a quadrature on a
uniform mesh is used in [20] and random collocation points in [21]. To solve
this issue, several approaches propose an adaptive distribution of points in
collocation-based PINNs. The authors in [22] design a residual-based adaptive
refinement process to add more collocation points in the locations with the
largest PDE residual, while a deep generative model that guides the sampling
of collocation points is proposed in [23]. The adaptive collocation point
movement method in [24] resamples the collocation points and uses an adaptive
loss weighting strategy. In [25], in addition to the strong residual of the
PDE, the gradient information of the NN is also used to guide the sampling
procedure. All these adaptive methods have demonstrated enhanced performance
and robustness compared to standard PINNs when solving PDEs with steep
gradients. Still, they suffer the applicability and robustness issues of
collocated PINNs mentioned above. The provably accurate adaptive quadratures
in [26] address the integration issue of NNs by automatically creating an
accurate quadrature rule for loss functions that involve the integration of
PDE residuals. However, the computation of the quadrature is computationally
expensive beyond 2D.

A posteriori error estimations are crucial for adaptive methods. Since the NN
in our method is interpolated onto a FE space, we can readily leverage the
well-studied FE error estimators [6, 27]. In recent years, NN-based a
posteriori error estimators have been developed. In [28], the authors propose
a NN error estimator based on the dual weighted residual method for PINNs.
Another related work is [29], where an error estimator for VPINNs consisting
of a residual-type term, a loss-function term, and data oscillation terms is
introduced. Functional error control techniques have been derived in [30] to
estimate the error in deep Galerkin simulations.

Remarkably, the state-of-the-art contains few investigations into the
application of NNs for problems with singular solutions. In [15, 17], the
authors use the FE interpolation of the trained NNs as the final solution of a
2D PDE with singularity. The deep Fourier residual (DFR) method [31] tackles
problems with singular solutions using NNs via adaptive domain decomposition.
However, the method is still prototypical in the sense that automatic
refinement is currently limited to 1D scenarios and relies on tensor products
of 1D orthonormal basis functions and an overkill integration procedure is
used.

In our previous study [17], we have demonstrated FEINNs outstanding
performance in solving various forward and inverse problems. In this work, we
carry out a detailed numerical analysis of the method for the FE residual
minimisation in a dual norm. Since the training of FEINNs on a fixed
underlying mesh limits their ability to leverage the full potential of NNs
non-linear capabilities, we propose the $h$-adaptive FEINN method to overcome
this limitation by introducing an _automatic_ adaptation of the interpolation
space (mesh) during training, driven by _a posteriori_ error estimation. We
propose a NN error estimator based on the strong PDE residual (which is
straightforward to implement and can be accurately computed via Gaussian
quadratures) and a standard FE Kelly indicator. These advances allow for a
more comprehensive exploitation of the non-linear nature of NNs, thus
enhancing the capability of the method to tackle problems with localised
features. The minimisation of the loss function relies on a train, estimate,
mark, adapt procedure, which is repeated till a certain level of accuracy is
reached. We leverage hierarchically adapted non-conforming octree-based meshes
and adaptive FE spaces with hanging node constraints to keep the conformity of
the interpolation.

We conduct a comprehensive set of numerical experiments, focusing on forward
problems featuring steep gradients in 2D and singularities in both 2D and 3D.
We verify that expressive enough NNs can replicate the rate of convergence of
the $h$-adaptive FEM solutions. In addressing the sharp-gradient problem, the
generalisation capabilities of the trained NNs yield remarkable results. The
non-interpolated NN solutions can be orders of magnitude more accurate than
their interpolation counterparts on the same mesh. This observation is
consistent with the findings in our previous work [17] for smooth solutions.
Our findings reveal that the performance of the NN error indicator is on par
with, if not surpassing, the traditional Kelly FE error estimator in various
situations. Moreover, we confirm that a well-posed residual norm in the loss
function has a tremendously positive impact on the convergence of the
$h$-adaptive training process.

The outline of the article is the following. Sec. 2 gives details on the
method that we propose, the model problem that we tackle, the FE
discretisation, the NN architecture, and the loss functions we aim to minimise
in the FEINN discretisation. In Sec. 3 we perform relevant numerical analysis
of the interpolated NNs on a fixed FE space and discuss the generalisation
error of the NNs. Sec. 4 describes the implementation of the $h$-adaptive
FEINN method and Sec. 5 presents numerical experiments on several forward
problems with sharp gradients and singularities. Finally, Sec. 6 draws
conclusions and lists potential directions for further research.

##  2\. Methodology

###  2.1. Continuous problem

Let us consider a Hilbert space $\tilde{U}$ on a Lipschitz polyhedral domain
$\Omega\subset\mathbb{R}^{d}$. We can consider a partition of the boundary
$\partial\Omega$ into a Dirichlet $\Gamma_{D}$ and Neumann $\Gamma_{N}$
component. Let $U$ be the subspace of $\tilde{U}$ with zero trace on
$\Gamma_{D}$. We denote by $(\cdot,\cdot)_{U}:U\times U\rightarrow\mathbb{R}$
an inner product in $U$ and by $\|\cdot\|_{U}$ the norm induced by the inner
product. We can define an abstract weak formulation of a PDE problem: find

$u\in U\ :\ a(u,v)=\ell(v)-a(\bar{u},v),\quad\forall v\in U,$ (1)

where $a:U\times U\rightarrow\mathbb{R}$ is a bilinear form and $\ell\in
U^{\prime}$ is a functional in the dual space of $U$. Non-homogeneous
Dirichlet boundary conditions are handled by a lifting $\bar{u}\in\tilde{U}$
while Neumann boundary conditions and body forces are included in $\ell$. The
final solution of the problem is $u+\bar{u}\in\tilde{U}$. In order for this
problem to be well-posed, the following inf-sup condition must be satisfied:
there exists a constant $\beta_{0}>0$ such that

$\inf_{w\in U}\sup_{v\in U}\frac{a(w,v)}{\|w\|_{U}\|v\|_{U}}\geq\beta_{0}.$
(2)

Besides, the bilinear form must be continuous, i.e., there exists a constant
$\gamma>0$ such that

$a(u,v)\leq\gamma\|u\|_{U}\|v\|_{U},\quad\forall u,v\in U.$ (3)

The solution of (1) vanishes the PDE residual

$\mathcal{R}({u})\doteq\ell(\cdot)-a({u}+\bar{u},\cdot)\in{U}^{\prime},$ (4)

In this work, we propose an $h$-adaptive FEINN formulation to discretise (1),
which can be applied to any PDE that fits in the previous framework. It covers
not only coercive problems but also indefinite and nonsymmetric problems. We
consider linear PDEs in the analysis, but the method can readily be applied to
non-linear problems. We stress, however, that we focus on _low-dimensional_
PDEs, e.g., $d\in\\{2,3\\}$. The numerical approximation of high-dimensional
PDEs face other problems, due to the curse of dimensionality, and the
requirements are different.

###  2.2. Finite element method

We consider a shape-regular partition $\mathcal{T}_{h}$ of $\Omega$. On
$\mathcal{T}_{h}$, we can define a trial FE space
$\tilde{U}_{h}\subset\tilde{U}$ and the subspace $U_{h}\subset U$ with zero
trace. We define a FE interpolant $\pi_{h}:\mathcal{C}^{0}\rightarrow U_{h}$
obtained by evaluation of the degrees of freedom of $U_{h}$. Similarly, we can
define the interpolant $\tilde{\pi}_{h}$ onto $\tilde{U}_{h}$. We can readily
define a FE lifting $\bar{u}_{h}\in\tilde{U}_{h}$ by interpolation on
$\Gamma_{D}$.

Let us consider a Petrov-Galerkin discretisation of (1) with a test space
$V_{h}\subset U$. The FE problem reads: find $\tilde{u}_{h}=\bar{u}_{h}+u_{h}$
where

$u_{h}\in U_{h}\ :\
a(u_{h},v_{h})=\ell(v_{h})-a(\bar{u}_{h},v_{h})\doteq\bar{\ell}(v_{h}),\quad\forall
v_{h}\in V_{h}.$ (5)

Problem (5) is well-posed under the assumption that the trial and test FE
spaces $U_{h}$, $V_{h}$ have the same dimension and satisfy the following inf-
sup condition [32]: there exists a constant $\beta>0$ independent of the mesh
size $h$ such that

$\inf_{w_{h}\in U_{h}}\sup_{v_{h}\in
V_{h}}\frac{a(w_{h},v_{h})}{\|w_{h}\|_{U}\|v_{h}\|_{U}}\geq\beta,\qquad\beta>0.$
(6)

We represent with $\mathcal{R}_{h}$ the restriction
$\mathcal{R}|_{{U}_{h}\times V_{h}}$.

###  2.3. Neural networks

In this work, we use fully-connected feed-forward NNs. These kind of networks
are constructed through a series of affine transformations and non-linear
activation functions. However, the only requirement in our analysis is that
the output can be written as a linear combination of the neurons at the
previous layer, allowing for the application of other architectures as well.
We represent the NN architecture by a tuple $(n_{0},\ldots
n_{L})\in\mathbb{N}^{(L+1)}$, where $L$ denotes the total number of layers and
$n_{k}$ means the number of neurons in layer $k$ with $0\leq k\leq L$. We take
$n_{0}=d$ and, for the scalar-valued PDEs we consider, we have $n_{L}=1$.

We denote by
$\boldsymbol{\Theta}_{k}:\mathbb{R}^{n_{k-1}}\to\mathbb{R}^{n_{k}}$ the affine
map from layer $k-1$ to $k$ ($1\leq k\leq L$), defined as
$\boldsymbol{\Theta}_{k}\boldsymbol{x}=\boldsymbol{W}_{k}\boldsymbol{x}+\boldsymbol{b}_{k}$,
where $\boldsymbol{W}_{k}\in\mathbb{R}^{n_{k}\times n_{k-1}}$ represents the
weight matrix and $\boldsymbol{b}_{k}\in\mathbb{R}^{n_{k}}$ the bias vector.
The activation function $\rho:\mathbb{R}\to\mathbb{R}$ is applied element-wise
after each affine map, excluding the final one. With these definitions in
place, the NN can be conceptualised as a parametrisable function
$\mathcal{N}(\boldsymbol{\theta}):\mathbb{R}^{d}\to\mathbb{R}$, expressed as:

$\mathcal{N}(\boldsymbol{\theta})=\boldsymbol{\Theta}_{L}\circ\rho\circ\boldsymbol{\Theta}_{L-1}\circ\ldots\circ\rho\circ\boldsymbol{\Theta}_{1},$
(7)

where $\boldsymbol{\theta}$ denotes the collection of all trainable parameters
$\boldsymbol{W}_{k}$ and $\boldsymbol{b}_{k}$. Besides, we apply the same
activation function across all layers. In this work, we represent the NN
architecture as $\mathcal{N}$ and a specific instance of the NN with
parameters $\boldsymbol{\theta}$ as $\mathcal{N}(\boldsymbol{\theta})$. We
extend the notation to also use $\mathcal{N}$ to represent the manifold of all
possible realisations of the NN, i.e.,
$\left\\{\mathcal{N}(\boldsymbol{\theta})\ :\
\boldsymbol{\theta}\in\mathbb{R}^{n_{\mathcal{N}}}\right\\}$,
$n_{\mathcal{N}}$ being the total number of parameters.

The NN structure in (7) can be split into two parts: the last-layer operation,
$\boldsymbol{\Theta}_{L}$, defined by parameters
$\boldsymbol{\theta}_{ll}\in\mathbb{R}^{n_{L}}$, and the hidden-layer
operations,
$\boldsymbol{\Theta}_{L-1}\circ\ldots\circ\rho\circ\boldsymbol{\Theta}_{1}$,
defined by parameters $\boldsymbol{\theta}_{hl}\in\mathbb{R}^{n_{hl}}$, with
$n_{hl}=n_{\mathcal{N}}-n_{L}$. We refer to the input of the last-layer
operation as the last-layer functions, consisting of $n_{L-1}$ output
functions from the hidden-layer operations. The output of a NN realisation can
be expressed as a linear combination of the last-layer functions, with
coefficients $\boldsymbol{W}_{L}$, plus a constant, the bias
$\boldsymbol{b}_{L}$. Thus, the hidden-layer parameters
$\boldsymbol{\theta}_{hl}$ define a linear space, i.e., the one spanned by the
last-layer functions. We denote this linear space as
$S(\boldsymbol{\theta}_{hl})$, and refer to it as to the last-layer space.
With this notation, the manifold of NN realisations can alternatively be
expressed as the union of all possible last-layer vector spaces:

$\mathcal{N}=\bigcup_{\boldsymbol{\theta}_{hl}\in\mathbb{R}^{n_{hl}}}S(\boldsymbol{\theta}_{hl}).$
(8)

In the FEINN method, we apply the $U_{h}$ interpolation operator
$\pi_{h}(\cdot)$ to the NN realisations. Since the interpolation operator is
linear, we can define the last-layer interpolated space
$S_{h}(\boldsymbol{\theta}_{hl})\subset U_{h}$ as the span of the
interpolation of the last-layer functions. The interpolated NN manifold
$\mathcal{N}_{h}\doteq\left\\{\pi_{h}(v)\ :\ v\in\mathcal{N}\right\\}\subset
U_{h}$ can be represented as the union of all last-layer interpolated spaces.
These definitions are used in Sec. 3.

###  2.4. Finite element interpolated neural networks

The FEINN method combines the NN in (7) and the FE problem in (5). The goal is
to find

$u_{\mathcal{N}}\in\underset{w_{\mathcal{N}}\in\mathcal{N}}{\mathrm{arg\,min}}\,\mathscr{L}(w_{\mathcal{N}}),\qquad\mathscr{L}(w_{\mathcal{N}})\doteq\left\lVert\mathcal{R}_{h}({\pi}_{h}(w_{\mathcal{N}}))\right\rVert_{X^{\prime}},$
(9)

where $u_{\mathcal{N}}$ is the NN approximation of the solution, and
$\left\lVert\cdot\right\rVert_{X^{\prime}}$ is some suitable dual norm.
Different norms will be considered in the numerical experiments in Sec. 5.
However, we only consider $X=U$ in the numerical analysis.

In general, it is unclear whether $\mathcal{N}$ and $\mathcal{N}_{h}$ are
closed spaces. If $\mathcal{N}_{h}$ is not closed, we cannot attain the
minimum in the optimisation problem (9), since we work on $\mathcal{N}_{h}$,
not its closure. To be mathematically precise, the FEINN problem must be
stated as: find

$u_{\mathcal{N}}\in\mathcal{N}\ :\
\pi_{h}(u_{\mathcal{N}})\in\underset{w_{\mathcal{N},h}\in\overline{\mathcal{N}_{h}}}{\mathrm{arg\,min}}\,\left\lVert\mathcal{R}_{h}(w_{\mathcal{N},h})\right\rVert_{X^{\prime}}.$
(10)

However, for any tolerance $\epsilon>0$, we can find

$u_{\mathcal{N},\epsilon}\in\mathcal{N}\ :\
\left\lVert\mathcal{R}_{h}(\pi_{h}(u_{\mathcal{N},\epsilon}))\right\rVert_{X^{\prime}}\leq\underset{w_{\mathcal{N},h}\in\overline{\mathcal{N}_{h}}}{\mathrm{arg\,min}}\,\left\lVert\mathcal{R}_{h}(w_{\mathcal{N},h})\right\rVert_{X^{\prime}}+\epsilon.$
(11)

######  Remark 2.1.

We refer to the solution of (11) as an $\epsilon$-_quasiminimiser_ of the
FEINN problem. The concept of quasiminimisers in the frame of PINNs has
already been previously discussed [33, 34]. This is of relative practical
relevance since the problem is solved iteratively up to a certain level of
accuracy and one can produce a sequence of $\epsilon$-_quasiminimisers_ that
converges to one of the solutions of the FEINN problem (10).

######  Remark 2.2.

It has been proven in [35] that the space of functions that can be generated
by a NN is not closed in $L^{p}$-norms for $0<p<\infty$ for most activation
functions. For the $L^{\infty}$-norm, and (parameterised) ReLU activation
functions, the NN space is proven to be closed in some cases. ReLU activation
functions are not smooth enough for PINNs or Deep Galerkin methods [26].
However, in FEINNs, $\mathcal{N}_{h}$ only depends on the pointwise values of
the neural network on the nodes of $U_{h}$ and ReLU activation functions are
suitable, since the method computes spatial derivatives of FE functions. Thus,
using the results in [35, Sec. 3.2], one can check that $\mathcal{N}_{h}$ is
closed for ReLU activation functions and the minimum in (10) is attained. This
is not the case for other activation functions (see Ex. 3.11).

As discussed in [17], we can define a discrete Riesz projector
$\mathcal{B}_{h}^{-1}:V_{h}^{\prime}\to V_{h}$ such that

$\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(w_{h})\in V_{h}\
:\left(\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(w_{h}),v_{h}\right)_{U}=\mathcal{R}_{h}(w_{h})(v_{h}),\quad\forall
v_{h}\in V_{h},$

and $\mathcal{B}_{h}^{-1}$ is (an approximation of) the inverse of the
corresponding Gramm matrix in $V_{h}$. Then, we can rewrite the loss function
in primal norm as follows:

$\mathscr{L}(u_{\mathcal{N}})=\left\lVert\mathcal{B}_{h}^{-1}\mathcal{R}_{h}({\pi}_{h}(u_{\mathcal{N}}))\right\rVert_{X}.$
(12)

Unlike loss functions that rely on discrete norms of the residual vector, this
loss function, _preconditioned_ by $\mathcal{B}_{h}$, is well-defined in the
limit $h\downarrow 0$. In practice, $\mathcal{B}_{h}^{-1}$ can be replaced by
any spectrally equivalent approximation to reduce computational costs; e.g.
one geometric multigrid (GMG) cycle has been proposed in [17].

######  Remark 2.3.

We note that the FE residual $\mathcal{R}_{h}$ in (9) is isomorphic to the
vector
$[\mathbf{r}_{h}(w_{h})]_{i}=\left<\mathcal{R}(w_{h}),\varphi^{i}\right>\doteq\mathcal{R}(w_{h})(\varphi^{i})$,
where $\\{\varphi^{i}\\}_{i=1}^{N}$ are the FE shape functions that span the
test space $V_{h}$. As a result, we can use the following loss function:

$\mathscr{L}(u_{\mathcal{N}})=\left\lVert\mathbf{r}_{h}(\pi_{h}(u_{\mathcal{N}}))\right\rVert_{\chi},$
(13)

where $\left\lVert\cdot\right\rVert_{\chi}$ is an algebraic norm of a vector.
The algebraic norm is ill-posed at the continuous level and, as a consequence,
the convergence of the minimisation process is expected to deteriorate as
$h\downarrow 0$ [36].

###  2.5. $h$-Adaptive finite element interpolated neural networks

The $h$-adaptive FEINN method is essentially the FEINN method equipped with a
dynamically/automatically adapting mesh. The method is sketched in Alg. 1. It
consists of a loop of the form:

$\text{TRAIN}\rightarrow\text{ESTIMATE}\rightarrow\text{MARK}\rightarrow\text{ADAPT}.$

At the beginning of the $i$-th loop iteration, the NN realisation
$u_{\mathcal{N}(\boldsymbol{\theta}_{i})}$ has initial parameters
$\boldsymbol{\theta}_{i}$. After training of
$u_{\mathcal{N}(\boldsymbol{\theta}_{i})}$ using the weak PDE loss built upon
the partition $\mathcal{T}_{i}$, the parameters of the NN become
$\boldsymbol{\theta}_{i+1}$. Then, we estimate the local error for each
element $K\in\mathcal{T}_{i}$ using an error indicator $\zeta_{K}$. Next, we
mark a set of elements $\mathcal{M}_{i}^{\rm r}$ for refinement and another
set $\mathcal{M}_{i}^{\rm c}$ for coarsening. We denote the refinement ratio
as $\delta^{\rm r}\in(0,1)$, and the coarsening ratio as $\delta^{\rm
c}\in(0,1)$. To determine $\mathcal{M}_{i}^{\rm r}$ and $\mathcal{M}_{i}^{\rm
c}$, we use Dörfler marking. Essentially, we select the top $\delta^{\rm r}$
elements with the largest estimator for refinement and the bottom $\delta^{\rm
c}$ with the smallest estimator for coarsening. We then adapt the marked
elements to construct a new partition $\mathcal{T}_{i+1}$. Finally, we start a
new iteration of the loop by interpolating
$u_{\mathcal{N}(\boldsymbol{\theta}_{i+1})}$ onto $\mathcal{T}_{i+1}$. This
iterative procedure continues until a stopping criterion is met.

Requirements: $k_{U}$, stopping criterion, $\delta^{\rm r}\in(0,1)$,
$\delta^{\rm c}\in(0,1)$

Input: $u_{\mathcal{N}(\boldsymbol{\theta}_{0})}$, $\mathcal{T}_{0}$

Output: $u_{\mathcal{N}(\boldsymbol{\theta}_{I+1})}$, $\mathcal{T}_{I}$

for  _$i=0,1,2,...$_ do

TRAIN $u_{\mathcal{N}(\boldsymbol{\theta}_{i})}$ by interpolating it into a FE
space of order $k_{U}$ built out of $\mathcal{T}_{i}$ using the loss defined
in (13) (or (12) for the preconditioned case); update NN parameters from
$\boldsymbol{\theta}_{i}$ to $\boldsymbol{\theta}_{i+1}$;

if  _stopping criterion is met_ then

break

end if

ESTIMATE the local error for each element $K\in\mathcal{T}_{i}$ by
$\zeta_{K}=\zeta_{K}(u_{\mathcal{N}(\boldsymbol{\theta}_{i+1})})$;

MARK a set $\mathcal{M}_{i}^{\rm r}\subset\mathcal{T}_{i}$ with largest
$\zeta_{K}$ such that $\\#\mathcal{M}_{i}^{\rm r}\approx\delta^{\rm
r}\\#\mathcal{T}_{i}$ and another set $\mathcal{M}_{i}^{\rm
c}\subset\mathcal{T}_{i}$ with smallest $\zeta_{K}$ such that
$\\#\mathcal{M}_{i}^{\rm c}\approx\delta^{\rm c}\\#\mathcal{T}_{i}$;

ADAPT the mesh: refine $K\in\mathcal{M}_{i}^{\rm r}$, and coarsen
$K\in\mathcal{M}_{i}^{\rm c}$ to construct a new partition of the domain
$\mathcal{T}_{i+1}$;

end for

$I\leftarrow i$;

return _$u_{\mathcal{N}(\boldsymbol{\theta}_{I+1})}$ , $\mathcal{T}_{I}$_

Algorithm 1 The $h$-adaptive FEINN procedure

Alg. 1 can be in principle combined with any kind of adaptive mesh approach.
Mainly because of performance and scalability reasons, in this work we use
hierarchically-adapted (i.e., nested) non-conforming octree-based meshes. We
provide more details on our particular implementation of Alg. 1 in Sec. 4.2.

###  2.6. Gradient-conforming discretisation

As a model problem in this work, we consider the Poisson equation. The problem
reads as (1) with $\tilde{U}=H^{1}(\Omega)$, $U=H^{1}_{0}(\Omega)$ and the
bilinear form and functional defined as follows:

$a(u,v)=\int_{\Omega}\boldsymbol{\nabla}(u)\cdot\boldsymbol{\nabla}(v),\quad\ell(v)=\int_{\Omega}fv,$

We consider an _exotic_ Petrov-Galerkin FE method. The well-posedness of the
discretisation is guaranteed by the discrete inf-sup condition (6) and the
same dimension for trial and test spaces (see [32]). We use a grad-conforming
nodal Lagrangian FE space $U_{h}$ of order $k_{U}$ built upon a (possibly non-
conforming) partition $\mathcal{T}_{h}$ of $\Omega$. The FE space $U_{h}$ is
built upon the partition $\mathcal{T}_{h}$. We use a _linearised_ test FE
space $V_{h}$ ($k_{V}=1$), built upon the partition $\mathcal{T}_{h}$ obtained
after $k_{U}$ levels of uniform refinement applied to $\mathcal{T}_{h}$. It is
easy to check that the dimension of ${U}_{h}$ and $V_{h}$ are identical. The
reason behind this choice of the test space is the dramatic improvement in the
convergence of the optimisation process compared to a Galerkin discretisation,
as observed, e.g., in [15, 17].

###  2.7. Error indicators

We introduce three types of error indicators to estimate the local errors of
each element. In this section, we restrict the exposition to the Laplacian
problem for the sake of simplicity. However, the error indicators can be
extended to other PDEs.

The first error indicator is the classic Kelly error estimator [27, 6], a
well-established indicator in $h$-adaptive FEM, particularly for the Poisson
equation. This error indicator estimates the error per element by integrating
the gradient jump of the NN interpolation along the faces of each element.
Specifically, for the element $K$, the Kelly error indicator is defined as

$\zeta_{K}^{k}=\sqrt{\sum_{F\in\partial
K}c_{F}\int_{F}\left\llbracket\frac{\partial({\pi}_{h}(u_{\mathcal{N}})+\bar{u}_{h})}{\partial
n}\right\rrbracket^{2}},$ (14)

where $c_{F}$ is a scale factor, $\partial K$ is the boundary of $K$, and
$\llbracket\cdot\rrbracket$ denotes the jump across the face $F$. The FE
offset $\bar{u}_{h}$ term in $\zeta_{K}^{k}$ can be excluded if the Dirichlet
boundary is exactly approximated by the employed finite elements, i.e.
$\bar{u}_{h}=g$. In the notation $\zeta_{K}^{k}$, the superscript $k$ stands
for “Kelly”, and the subscript $K$ denotes an element of the mesh. We use
$c_{F}=h_{K}/24$ as proposed in [6], where $h_{K}$ is the characteristic size
of element $K$.

Then we introduce another error indicator that leverages the NN and the strong
form of the PDE. The error estimator at element $K$ is defined as

$\zeta_{K}^{n}=\sqrt{\int_{K}|\boldsymbol{\Delta}u_{\mathcal{N}}+f|^{2}}.$
(15)

The superscript $n$ in $\zeta_{K}^{n}$ stands for “network”. We note that
$\zeta_{K}^{n}$ is similar to the mean PDE residual $\varepsilon_{r}$ defined
in [22]. The main difference is that $\zeta_{K}^{n}$ is evaluated with a high-
degree Gauss quadrature, while $\varepsilon_{r}$ is computed by Monte Carlo
integration. At each mesh adaptation step, the strategy in [22] adds more
collocation points only on the region with the highest $\varepsilon_{r}$.
Instead, our approach can refine multiple elements with the largest errors and
coarsen those with the smallest errors simultaneously.

Finally, we include the real error of the interpolated NN as the third
indicator, which is defined for element $K$ as:

$\zeta_{K}^{r}=\sqrt{\int_{K}|{\pi}_{h}(u_{\mathcal{N}})+\bar{u}_{h}-u|^{2}},$

where $u$ is the exact solution. The superscript $r$ in $\zeta_{K}^{r}$ stands
for “real”. In practice, the real error is unavailable. In this study, we use
it as a reference to assess and compare the performance of the other two
indicators.

##  3\. Numerical analysis

In this section, we extend the analysis in [17] in multiple ways. In
particular, in Sec. 3.1, we bound the error of the interpolated NN under the
general scenario in which the NN cannot exactly emulate the FE space $U_{h}$,
i.e., $\mathcal{N}_{h}\subsetneq U_{h}$ (or, equivalently, the interpolation
operator $\pi_{h}:\mathcal{N}\mapsto U_{h}$ is not surjective). Remarkably, we
prove that the error of the interpolated NN is bounded by the interpolation
error among all possible vector spaces that result from the interpolation of
the last-layer functions in the NN architecture. In Sec. 3.2, we perform error
analysis under a weakened assumption of exact emulation which we refer to as
quasi-emulation. This analysis is relevant as exact emulation might not be
feasible for some activation functions, regardless of the expressivity of the
NN. Finally, in Sec. 3.3, we also discuss the well-posedness of the FEINN
problem in a quotient space and discuss the accuracy of the non-interpolated
NN by considering the maximum distance between elements of an equivalence
class.

###  3.1. Error analysis

The error analysis in this section relies on the discrete inf-sup condition
(6), and can readily be applied to indefinite systems. In the following
theorem, we state the well-posedness of the residual minimisation problem in a
FE setting.

######  Theorem 3.1.

Let us consider a pair of FE spaces $U_{h}$ and $V_{h}$ that satisfy the inf-
sup condition (6). The discrete quadratic minimisation problem

$u_{h}=\underset{w_{h}\in
U_{h}}{\mathrm{arg}\,\mathrm{min}}\|r_{h}(w_{h})\|_{U},\qquad
r_{h}(w_{h})\doteq\mathcal{B}_{h}^{-1}\mathcal{R}_{h}(w_{h}),$ (16)

is well-posed. The problem can be re-stated as: find $u_{h}\in U_{h}$ and
$r_{h}=r_{h}(u_{h})\in V_{h}$ such that

$\displaystyle(r_{h},v_{h})_{U}+a(u_{h},v_{h})$
$\displaystyle=\bar{\ell}(v_{h}),\quad$ $\displaystyle\forall v_{h}\in V_{h},$
(17) $\displaystyle a(w_{h},r_{h})$ $\displaystyle=0,\quad$
$\displaystyle\forall w_{h}\in U_{h}.$ (18)

The solution satisfies the following a priori estimates:

$\|r_{h}\|_{U}\leq\underset{w_{h}\in
U_{h}}{\mathrm{inf}}4\gamma\|u-w_{h}\|_{U},\qquad\|u-u_{h}\|_{U}\leq\left(1+\frac{4\gamma}{\beta}\right)\underset{w_{h}\in
U_{h}}{\mathrm{inf}}\|u-w_{h}\|_{U}.$ (19)

###### Proof.

First, we note that the loss function is differentiable (see [17, Prop. 3.1]).
To obtain (17)-(18), we compute the variational derivative of (16) to
variations of $u_{h}$ and recall the definition of the discrete Riesz
projector. Under the assumption of the discrete inf-sup condition (6), we can
prove the well-posedness of the problem using the Babuska-Brezzi theory.
Namely, using (6), we can readily find $k_{h}\in V_{h}$ such that
$\|k_{h}\|_{U}=\|u_{h}\|_{U}$ and $a(u_{h},k_{h})\geq\beta\|u_{h}\|_{U}^{2}$.
We can take $v_{h}=r_{h}+\alpha k_{h}$ for an arbitrary $\alpha>0$ and
$w_{h}=-u_{h}$, to get

$\displaystyle\|r_{h}\|^{2}_{U}+\alpha a(u_{h},k_{h})+\alpha(r_{h},k_{h})_{U}$
$\displaystyle\geq\|r_{h}\|^{2}_{U}+\alpha\beta\|u_{h}\|^{2}_{U}-\frac{\alpha\beta}{2}\|u_{h}\|_{U}^{2}-\frac{\alpha}{2\beta}\|r_{h}\|^{2}_{U}$
(20)
$\displaystyle=\left(1-\frac{\alpha}{2\beta}\right)\|r_{h}\|^{2}_{U}+\frac{\alpha\beta}{2}\|u_{h}\|^{2}_{U}.$
(21)

Taking $\alpha=\beta$, we get, using (17):

$\displaystyle\frac{1}{4}(\|r_{h}\|_{U}+\beta\|u_{h}\|_{U})^{2}$
$\displaystyle\leq\frac{1}{2}\|r_{h}\|^{2}_{U}+\frac{\beta^{2}}{2}\|u_{h}\|^{2}_{U}\leq\|r_{h}\|^{2}_{U}+\beta
a(u_{h},k_{h})+\beta(r_{h},k_{h})_{U}$ (22)
$\displaystyle\leq\|\bar{\ell}\|_{U^{\prime}}(\|r_{h}\|_{U}+\beta\|u_{h}\|_{U}).$
(23)

Using the fact that $\bar{\ell}(v_{h})=a(u,v_{h})$, and subtracting
$a(w_{h},v_{h})$ on both sides of (17) for an arbitrary $w_{h}\in U_{h}$ and
using the continuity of the bilinear form $a$, we get the following error
estimate from (22):

$\|r_{h}\|_{U}+\beta\|u_{h}-w_{h}\|_{U}\leq 4\gamma\|u-w_{h}\|_{U}.$ (24)

Using the triangle inequality, we readily get the second estimate in (19). It
proves the theorem. ∎

######  Remark 3.2.

We note that $dim(U_{h})\neq dim(V_{h})$ in general, and the trial space can
be of lower dimension than the test space. In the case in which
$dim(U_{h})=dim(V_{h})$, the solution $u_{h}$ of (17)-(18) is the solution of
(1) and $r_{h}=0$.

The previous theorem readily applies to any subspace $W_{h}\subseteq U_{h}$,
since the inf-sup condition readily holds for $W_{h}\times V_{h}$. We use the
notation $u(W_{h})$ and $r(W_{h})$ to denote the unique solution of (17)-(18)
on $W_{h}\times V_{h}$. We can prove the following result.

######  Proposition 3.3.

Let us consider $u_{h,1},\,u_{h,2}\in U_{h}$. The following error estimate
holds:

$\frac{1}{\beta}\|u_{h,1}-u_{h,2}\|_{U}\leq\|r(u_{h,1})-r(u_{h,2})\|_{U}\leq\gamma\|u_{h,1}-u_{h,2}\|_{U}.$
(25)

###### Proof.

Let us denote for conciseness $u_{h}\doteq u_{h,1}$, $u_{h}^{\prime}\doteq
u_{h,2}$, $r_{h}\doteq r(u_{h,1})$, $r_{h}^{\prime}\doteq r(u_{h,2})$. Using
the definition of $r_{h}(\cdot)$ in (16), we have:

$\displaystyle(r_{h},v_{h})_{U}$
$\displaystyle=-a(u_{h},v_{h})+\bar{\ell}(v_{h}),$
$\displaystyle(r_{h}^{\prime},v_{h})_{U}$
$\displaystyle=-a(u_{h}^{\prime},v_{h})+\bar{\ell}(v_{h}),\qquad\forall
v_{h}\in V_{h}.$

Subtracting the previous equations, we get

$(r_{h}-r_{h}^{\prime},v_{h})_{U}=-a(u_{h}-u_{h}^{\prime},v_{h}),\qquad\forall
v_{h}\in V_{h}.$ (26)

The upper bound can readily be obtained using the continuity of $a$ in (3) ,
while the lower bound is obtained using the discrete inf-sup condition (6) and
the Cauchy-Schwarz inequality. ∎

######  Proposition 3.4.

Let us consider $y_{h}\in U_{h}$ and $W_{h}\subseteq U_{h}$ such that
$\|r(y_{h})\|_{U}\leq\|r_{h}(W_{h})\|_{U}+\epsilon$, $\epsilon\geq 0$. The
following a priori error estimate holds:

$\|u-y_{h}\|_{U}\leq\rho\inf_{w_{h}\in
W_{h}}\|u-w_{h}\|_{U}+\beta\epsilon,\qquad\rho\doteq
1+\frac{4\gamma}{\beta}+{8\gamma}{\beta}.$ (27)

###### Proof.

Using Prop. 3.3, the triangle inequality, the statement of the proposition and
the error estimate in (19), we get:

$\frac{1}{\beta}\|y_{h}-u(W_{h})\|_{U}\leq\|r(y_{h})-r(W_{h})\|_{U}\leq
2\|r(W_{h})\|_{U}+\epsilon\leq 8\gamma\inf_{w_{h}\in
W_{h}}\|u-w_{h}\|_{U}+\epsilon.$

Then, using the triangle inequality and (19), we get:

$\|u-y_{h}\|_{U}\leq\|_{U}u-u(W_{h})\|_{U}+\|y_{h}-u(W_{h})\|_{U}\leq\rho\inf_{w_{h}\in
W_{h}}\|u-w_{h}\|_{U}+\epsilon.$

∎

Next, we prove that the infimum (10) (or in (11)) is the same if we restrict
to the space of solutions of (17)-(18) for all last-layer interpolated spaces,
i.e.,
$\mathcal{M}_{h}\doteq\left\\{u_{h}(S_{h}(\boldsymbol{\theta}_{hl})):\boldsymbol{\theta}_{hl}\in\mathbb{R}^{n_{hl}}\right\\}$.

######  Proposition 3.5.

The following equalities hold:

$\underset{w_{\mathcal{N}}\in\mathcal{N}}{\mathrm{inf}}\|r_{h}(\pi_{h}(w_{\mathcal{N}}))\|_{U}=\underset{w_{h}\in\mathcal{N}_{h}}{\mathrm{inf}}\|r_{h}(w_{h})\|_{U}=\underset{u_{h}\in\mathcal{M}_{h}}{\mathrm{inf}}\|r_{h}(u_{h})\|_{U}.$
(28)

###### Proof.

The first equality is an obvious result of the definition for
$\mathcal{N}_{h}$. The second equality uses the fact that $\mathcal{N}_{h}$
can be represented as the union of all last-layer spaces. Due to Th. 3.1,
among all functions in $S_{h}$, there exists a unique solution $u_{h}(S_{h})$
that minimises the functional. It is the solution of (17)-(18) for
$S_{h}\times V_{h}$ and its residual is denoted by $r_{h}(S_{h})$. It proves
the result. ∎

######  Remark 3.6.

To consider the minimisation on $\mathcal{M}_{h}$ is practically possible. It
involves to solve a linear system for the last-layer weights of the NN, e.g.,
using an iterative Krylov method. It motivates the usage of hybrid solvers
that combine non-convex optimisation (e.g., Adam) and linear solvers. If not
at every epoch, a final linear solve can be performed to obtain the optimal
solution in the last-layer space.

Let us consider the closure $\overline{\mathcal{N}_{h}}$ of $\mathcal{N}_{h}$
in $U$. $U_{h}$ is closed as it is a finite-dimensional vector space on a
complete normed field ($\mathbb{R}$ or $\mathbb{C}$). Thus,
$\overline{\mathcal{N}_{h}}\subseteq U_{h}$. We can prove the following
result.

######  Theorem 3.7.

For any $\epsilon>0$ and $\epsilon$-quasiminimiser $u_{\mathcal{N},\epsilon}$
of (11), the following a priori error estimate holds:

$\|u-\pi_{h}(u_{\mathcal{N},\epsilon})\|_{U}\leq\rho\inf_{w_{h}\in\mathcal{N}_{h}}\|u-w_{h}\|_{U}+\beta\epsilon.$
(29)

###### Proof.

Let us consider the minimisation problem on $\overline{\mathcal{N}_{h}}$. The
infimum is attained since the space is closed.

Thanks to Prop. 3.5, the minimum in (10) is attained for some
$u_{h}\in\overline{\mathcal{M}_{h}}$. Thus, for any $\epsilon$-quasiminimiser
$u_{h,\epsilon}$ of (11), it holds
$\|r_{h}(u_{h,\epsilon})\|_{U}\leq\|r_{h}(u_{h})\|_{U}+\epsilon\leq\|r_{h}(v_{h})\|_{U}+\epsilon$
for any $v_{h}\in\mathcal{M}_{h}$. Invoking Prop. 3.4 and (8), we prove the
result. ∎

###  3.2. Quasi-emulation

In the previous section, we have considered the case in which the NN space
cannot emulate the FE space, i.e., $\mathcal{N}_{h}\subsetneq U_{h}$. If the
NN is expressive enough, we can prove that $\mathcal{N}_{h}=U_{h}$ following
the construction in [17, Prop. 3.2].

######  Proposition 3.8.

Let us consider an adaptive Lagrangian FE space $U_{h}$ in $\mathbb{R}^{d}$
with $N$ nodes, which can be defined on a non-conforming mesh and include
hanging nodes. Let $\mathcal{N}$ be a neural network architecture with ReLU
activation function, 3 hidden layers and $(3dN,dN,N)$ neurons per layer. The
interpolation operator $\pi_{h}:\mathcal{N}\mapsto U_{h}$ is surjective.

###### Proof.

The proof follows the same lines as in [17, Prop. 3.2]. The neural network
above can represent in the last layer all the shape functions in the adaptive
FE space $U_{h}$, both the free ones as the ones that are constrained and
associated to hanging nodes. The linear combination of these last layer shape
functions can represent any function in $U_{h}$. As a result, the
interpolation operator $\pi_{h}$ is surjective. ∎

######  Remark 3.9.

The previous result is overly pessimistic and can be improved when using
meshes that are obtained after adaptive mesh refinement of a regular mesh. In
this case, the FE space $U_{h}$ can be represented by a NN with a much smaller
number of neurons, since one can exploit a tensor product structure in the
proof (see [17]).

Exact emulation is not attainable for some activation functions, like $\tanh$,
regardless of the expressivity of the NN. In this case, we introduce the
concept of quasi-emulation.

######  Definition 3.10.

A NN is a $\delta$-emulator of a vector space $U_{h}$ if for all $w_{h}\in
U_{h}$, there exists a $w\in\mathcal{N}$ such that
$\|\pi_{h}(w)-w_{h}\|_{U}\leq\delta$.

######  Example 3.11.

Let us consider a 1D problem for simplicity and a fully-connected NN as in
Sec. 2.3 with $\rho=\tanh$ activation, $m\in\mathbb{N}$ and the weights and
biases:

$\boldsymbol{W}_{1}=\begin{bmatrix}m,\ -m\end{bmatrix}^{T},\
\boldsymbol{b}_{1}=\begin{bmatrix}\frac{m}{4},\
\frac{m}{4}\end{bmatrix}^{T},\quad\boldsymbol{W}_{2}=\begin{bmatrix}m,\
m\end{bmatrix},\ \boldsymbol{b}_{2}=\begin{bmatrix}0\end{bmatrix}.$

One can check that the corresponding realisation
$f_{m}({x})=\rho\circ(\boldsymbol{W}_{2}(\rho\circ(\boldsymbol{W}_{1}{x}+\boldsymbol{b}_{1})+\boldsymbol{b}_{2}))$
satisfies $f_{m}({x}=0)=1-\delta$, $f_{m}({x}=-1)=f_{m}({x}=1)=c\delta$ for
some constant $c$, and monotonically decreases as $x$ goes to $\pm\infty$.
Thus, $\lim_{m\to\infty}f_{m}(x)=h(x)$, where $h(x)$ is the hat function with
support in $[-1,1]$ and value 1 at $x=0$. We can readily obtain a hat function
in any interval by applying translation and/or scaling. Thus, we can construct
a NN with a set of last-layer functions that emulates a FE space basis up to a
desired tolerance $\delta>0$. Similar constructions are possible in higher
dimensions.

Due to the well-posedness of the FE problem in Th. 3.1, there is a unique
minimum $u_{h}\in U_{h}=\overline{\mathcal{N}_{h}}$ of (17)-(18). For exact
emulation, $\mathcal{N}_{h}=\overline{\mathcal{N}_{h}}$. Thus, we can attain
minimisers

$u_{\mathcal{N}}\in\mathcal{N}\ :\
\pi_{h}(u_{\mathcal{N}})\in\underset{w_{\mathcal{N},h}\in{\mathcal{N}}_{h}}{\mathrm{arg\,min}}\,\left\lVert\mathcal{R}_{h}(w_{\mathcal{N},h})\right\rVert_{X^{\prime}}.$

In fact, the minimisers are all the $u_{\mathcal{N}}\in\mathcal{N}$ such that
$\pi_{h}(u_{\mathcal{N}})=u_{h}$ and the a priori error estimate (19) holds.
For quasi-emulation, the following result is a direct consequence of Th. 3.1
and the previous definition.

######  Corollary 3.12.

Let us consider a NN $\mathcal{N}$ that is a $\delta$-emulator of the FE space
$U_{h}$. For any $\epsilon>0$, an $\epsilon$-quasiminimiser $u_{\mathcal{N}}$
of (11) satisfies the following a priori error estimate:

$\|u-\pi_{h}(u_{\mathcal{N}})\|_{U}\leq\rho\inf_{w_{h}\in
U_{h}}\|u-w_{h}\|_{U}+\beta\epsilon.$ (30)

###  3.3. Generalisation error and equivalence classes

In the discussions above, we have only bounded the error of the interpolated
NN. Analogously, we have discussed the well-posedness in the space of
interpolations of the NN. The total error of the NN is harder to bound, since
it is _transparent_ to the loss.

In order to better bound the total error of the NN, we can make use of
quotient spaces (see a related approach in [18] for VPINNs). Let us define the
following equivalence relation $\sim$: two NN realisations
$\mathcal{N}(\boldsymbol{\theta})$ and
$\mathcal{N}(\boldsymbol{\theta}^{\prime})$ are equivalent if their
interpolations on the adaptive FE space $U_{h}$ are identical, i.e.,
$\pi_{h}(\mathcal{N}(\boldsymbol{\theta}))=\pi_{h}(\mathcal{N}(\boldsymbol{\theta}^{\prime}))$.
The equivalence class of $\mathcal{N}(\boldsymbol{\theta})$ is denoted by
$[\mathcal{N}(\boldsymbol{\theta})]$. We can define the quotient space
$\mathcal{N}/\sim$ as the set of all equivalence classes
$[\mathcal{N}(\boldsymbol{\theta})]$.

Let us consider the case of exact emulation for simplicity. In this situation,
the problem (10) admits a unique solution $u_{h}$ in $U_{h}$. Then, the
problem is also well-posed in the quotient space $\mathcal{N}/\sim$; the
unique solution is the equivalence class $[u_{h}]\in\mathcal{N}/\sim$.
Analogously, in the case of quasi-emulation, quasi-minimisers in the quotient
space will be arbitrarily close to $[u_{h,\epsilon}]$.

Assuming that the training process converges to the global minimiser, we can
prove the following result.

######  Proposition 3.13.

Let $u_{h}\in U_{h}$ be the solution of (5). Let us assume that the training
process has attained an $\epsilon$-quasiminimiser
$u_{h,\epsilon}\in\mathcal{N}_{h}$, or analogously, a
$[u_{h,\epsilon}]\in\mathcal{N}/\sim$. Any possible realisation
$u_{\mathcal{N}}\in[u_{h,\epsilon}]\subset\mathcal{N}$ obtained after training
satisfies the following a priori error result:

$\left\lVert u-u_{\mathcal{N}}\right\rVert_{U}\leq\left\lVert
u-u_{h}\right\rVert_{U}+\sup_{v\in[\pi_{h}(u_{\mathcal{N}})]}\left\lVert\pi_{h}^{\perp}(v)\right\rVert_{U}+\epsilon,$
(31)

where $\pi_{h}^{\perp}(u)=u-\pi_{h}(u)$.

###### Proof.

Using the triangle inequality and taking the supremum of the error between the
NN and the FE solution in the NN architecture, we prove the result. ∎

The first term in the error estimate is the FE error, the second term is the
generalisation error, and the third term is a training error. The
generalisation error describes the maximum distance between the FE solution
and the NN realisations in the corresponding equivalence class. Since the loss
function does not _feel_ $\|\pi_{h}^{\perp}(u_{\mathcal{N}})\|$, there is no
explicit control on it. This error is related to the expressiveness of the NN
with respect to the FE mesh. The more expressive the NN, the larger the
generalisation error can be, as it is well-known in the machine learning
literature (over-fitting). Besides, if the NN is too coarse, the emulation
assumption will not hold and the FEINN will not be able to capture the FE
solution. 111 This discussion motivates the usage of a regularisation term
that controls the distance between the NN and the FE solution. For example,
one could consider a term of the order
$h^{-2}\|\pi^{\perp}(u_{\mathcal{N}})\|_{L^{2}(\Omega)}$ or
$\|\pi^{\perp}(u_{\mathcal{N}})\|_{H^{1}(\Omega)}$, integrated with a high
order quadrature. Such a term is weakly consistent, i.e., it does not affect
the convergence of FEINNs compared to FEM. This stabilisation resembles the
so-called orthogonal subscale stabilisation common in the FE literature to
circumvent the inf-sup condition (see, e.g., [37, 38]). We have not considered
stabilisation in this work.

The previous analysis is FE-centric, in the sense that we consider the NN
error by comparing it with the FE solution and use FE theory to quantify the
error. Even though the FEINN can capture the FE solution with a much smaller
number of parameters than the number of DoFs in the FE space [17], it can be
argued that if the only aim is to attain the FE error, it is more effective to
solve the FE problem. However, the FEINN approach can naturally deal with
inverse problems by simply adding data misfit terms to the loss function and
the NN parametrisation regularises the inverse problem and can readily be used
in a one-shot fashion, avoiding the need of adjoint problems [17]. It has also
been observed experimentally in [17] that the FE-centric analysis of FEINNs is
overly pessimistic in general. The trained NN using the FEINN approach can be
up to two orders of magnitude more accurate than the FE solution for forward
problems. In Sec. 5, we analyse this fact in the case of singular solution and
sharp gradients combined with $h$-adaptivity and a posteriori error
estimation. A mathematical justification of the superior accuracy of the FEINN
trained NN compared to the FEM solution is still open.

The two error estimates in Sec. 2.7 are related to these two different
perspectives. The a posteriori error estimate (14) aims to reduce the error of
the interpolated NN, following the FE-centric perspective. Instead, the
estimate (15) aims to directly reduce the error between the NN and the exact
solution, considering the interpolation just as a way to enforce the Dirichlet
condition, define a trial space to compute the weak residual and perform an
approximate numerical integration.

##  4\. Implementation

###  4.1. Automatic mesh adaptation using forest-of-octrees

Automatic mesh adaptation is a key component of $h$-adaptive FEINNs. In this
work, we leverage a particular kind of hierarchically-adapted meshes known as
forest-of-octrees (see, e.g., [39, 40]). Forest-of-octrees meshes can be seen
as a two-level decomposition of $\Omega$, referred to as macro and micro
level, respectively. The macro level is a suitable conforming partition
$\mathcal{C}$ of $\Omega$ into quadrilateral ($d=2$) or hexahedral cells
($d=3$). This mesh, which may be generated, e.g., using an unstructured mesh
generator, is referred to as the coarse mesh. At the micro level, each of the
cells of $\mathcal{C}$ becomes the root of an adaptive octree with cells that
can be recursively and dynamically refined or coarsened using the so-called
$1:2^{d}$ uniform partition rule. If a cell is marked for refinement, then it
is split into $2^{d}$ children cells by subdividing all parent cell edges. If
all children cells of a parent cell are marked for coarsening, then they are
collapsed into the parent cell. The union of all leaf cells in this hierarchy
forms the decomposition of the domain at the micro level, i.e., $\mathcal{T}$.

In this work, we utilise the GridapP4est.jl [41] Julia package in order to
handle such kind of meshes. This package, built upon the p4est meshing engine
[39], is endowed with the so-called Morton space-filling curves, and it
provides high-performance and low-memory footprint algorithms to handle
forest-of-octrees.

###  4.2. Non-conforming interpolation of neural networks

Unlike our previous work [17], which relies on (uniformly-refined) conforming
meshes, in this work we interpolate the NN onto (grad-conforming) FE spaces
built out of non-conforming forest-of-octrees meshes. With this kind of
meshes, generating a cell-conforming global enumeration of DoFs is not
sufficient to impose the required continuity at the interface of cells with
different refinement levels. In order to remedy this, it is standard practice
to equip the FE space with additional linear multi-point constraints at the
DoFs lying on hanging vertices, edges and faces. The structure and set up of
such constraints is well-established knowledge, and thus not covered here;
see, e.g., [40, 42] for further details. In any case, the presence of these
linear multi-point constraints in the FE space implies that the interpolation
of the NN at hanging DoFs cannot by given by its evaluation at the
corresponding nodes, but rather by a linear combination of the NN evaluations
at the nodes of the adjacent coarser cell. To this end, GridapP4est.jl
provides the required functionality to build and apply these linear multi-
point constraints, and to resolve these constraints during FE assembly.

###  4.3. Finite element interpolated neural networks

In this section we briefly overview the main building blocks of our Julia
implementation of $h$-adaptive FEINNs. For a more complete coverage of these
(including the inverse problem case), we refer to [17]. At a purely algebraic
level, we can rewrite the preconditioned loss (12) as:

$\mathscr{L}(\boldsymbol{\theta})=\left\lVert
r_{h}(\boldsymbol{\theta})\right\rVert_{X},$

where $r_{h}\in V_{h}$ is the discrete Riesz projection of the FE residual.
The vector of DoFs of $r_{h}$ is
$\mathbf{B}^{-1}(\mathbf{A}\mathbf{u}_{h}-\mathbf{f})$, where $\mathbf{B}$ is
the discretisation in the test space $V_{h}$ of the Gram matrix associated to
the inner product in $U$, e.g., the Laplacian matrix for the Poisson problem,
$\mathbf{A}$ is the coefficient matrix resulting from discretisation,
$\boldsymbol{\theta}$ is the parameters for the NN $u_{\mathcal{N}}$,
$\mathbf{u}_{h}$ is the vector of DoFs of $U_{h}$, and $\mathbf{f}$ is the
right-hand side (RHS) vector. We use the abstract data structures in Gridap.jl
[43, 44] as implemented by GridapP4est.jl to compute $\mathbf{A}$,
$\mathbf{f}$, and the FE nodal coordinates $\mathbf{X}$ of those DoFs which
are not subject to linear multi-point constraints. We utilise Flux.jl [45, 46]
to build $u_{\mathcal{N}}$ and evaluate $u_{\mathcal{N}}$ at $\mathbf{X}$ to
compute $\mathbf{u}_{h}$. Without preconditioning, i.e.,
$\mathbf{B}=\mathbf{I}$, we can minimise the loss
$\left\lVert\mathbf{A}\mathbf{u}_{h}(\boldsymbol{\theta})-\mathbf{f}\right\rVert_{\chi}$
directly. During training, we rely on Zygote.jl [47] to compute the gradient
of the loss with respect to $\boldsymbol{\theta}$ (see [17] for details). We
note that in the context of $h$-adaptive FEINNs, we need to recompute
$\mathbf{A}$, $\mathbf{f}$, and $\mathbf{u}_{h}$ at each adaptation step.
Furthermore, depending on the specific problem, these arrays can be stored
cell-wise to improve efficiency.

If a preconditioner is employed, i.e., $\mathbf{B}\neq\mathbf{I}$, we need to
perform these two additional steps at each iteration:

  1. (1)

Compute $r_{h}$ by solving
$\mathbf{B}^{-1}(\mathbf{A}\mathbf{u}_{h}(\boldsymbol{\theta})-\mathbf{f})$ on
$V_{h}$.

  2. (2)

Compute the preconditioned loss $\mathscr{L}(\boldsymbol{\theta})=\left\lVert
r_{h}(\boldsymbol{\theta})\right\rVert_{X}$.

The forward and backward passes of the above steps are both handled by
Gridap.jl. It is noteworthy that in step (1), we do not explicitly compute the
inverse of $\mathbf{B}$, but rather solve the linear system
$\mathbf{B}\mathbf{r}_{h}=\mathbf{A}\mathbf{u}_{h}-\mathbf{f}$ to compute the
vector of DoFs $\mathbf{r}_{h}$ that uniquely determine $r_{h}\in V_{h}$.
Considering the problem linear, we can reuse the factorisation of
$\mathbf{B}$, leading to linear complexity computational cost in the
application of $\mathbf{B}^{-1}$ per iteration. Alternatively, as explored in
[17], we can use a cheaper preconditioner $\mathbf{B}$, such as the GMG
preconditioner, which has linear complexity in the number of DoFs of the FE
space. Besides, given that $r_{h}\in V_{h}$ and $V_{h}$ is a linearised FE
space, and the gradient of $r_{h}$ is piecewise constant, the computation of
$\left\lVert r_{h}\right\rVert_{X}$ is very cheap.

##  5\. Numerical experiments

Tab. 1 lists the NN architectures and training iterations for the different
problems tackled in this section. For 2D experiments, we adopt the NN
architecture in [17], namely $L=5$ layers and $n=50$ neurons for each hidden
layer. For 3D experiments, since we work with a more complex problem with more
DoFs, we increase the number of layers to $L=10$. In addition, we use
$\rho=\tanh$ as the activation function, the Glorot uniform method [48] as the
NN parameter initialisation strategy. We employ the BFGS optimiser in Optim.jl
[49] to train the NNs, as recommended in [17]. Moreover, in our preliminary
experiments, we have observed that the results are generally robust to
variations in NN initialisations, provided that the initial mesh is
sufficiently fine and the NN is adequately trained. As a result, at each
adaptation step, we heuristically set the number of the training iterations as
a function of the number of DoFs in the trial FE space. In particular, the
higher the number of DoFs, the higher the number of training iterations. We
acknowledge that advanced/automatic optimiser termination conditions could
result in better accuracy/performance trade-offs than the ones shown in this
section. However, developing such conditions is beyond the scope of this
paper.

Section number and problem | # layers | # neurons | # BFGS iterations |  DoF milestones
---|---|---|---|---
Sec. 5.1 2D arc wavefront | 5 | 50 | [500, 1000, 1500] | [5000, 10000]
Sec. 5.2 2D L-shaped singularity | 5 | 50 | [3000, 4000, 5000] | [10000, 20000]
Sec. 5.3 3D singularities | 10 | 50 | [3000, 4000, 5000] | [50000, 100000]
Table 1. Example of selected hyperparameter values for the three types of
problems tested in this study. The numbers shown correspond to $k_{U}=4$
assuming no preconditioner. For different values of $k_{U}$ and/or
preconditioner strategy, please see the corresponding subsection. The number
of training iterations at each adaptation step are determined by DoF
milestones (last column of the table). For instance, in the 2D arc wavefront
problem, we use a fixed number of 500 iterations up to 5,000 DoFs, 1,000
iterations within 5,000-10,000 DoFs, and 1,500 iterations over 10,000 DoFs.

Regarding the refinement and coarsening ratios, for simplicity, we just use
fixed ratios for all adaptation steps, even though these ratios might be
different for each step in practice. We use refinement ratio $\delta^{\rm
r}=0.15$ for 2D cases and $\delta^{\rm r}=0.1$ for 3D cases. For both space
dimensions, the coarsening ratio is chosen to be $\delta^{\rm c}=0.01$.
Besides, for the stopping criterion, we simply set a maximum number of
adaptation steps: 7 for 2D problems and 4 for 3D problems.

We use the $L^{2}(\Omega)$ and $H^{1}(\Omega)$ error norms to evaluate the
accuracy of the identified solution $u^{id}$:

$e_{L^{2}(\Omega)}(u^{id})=\left\lVert
u-u^{id}\right\rVert_{L^{2}(\Omega)},\qquad
e_{H^{1}(\Omega)}(u^{id})=\left\lVert u-u^{id}\right\rVert_{H^{1}(\Omega)},$

where $u$ is the true solution,
$\left\lVert\cdot\right\rVert_{L^{2}(\Omega)}=\sqrt{\int_{\Omega}|\cdot|^{2}}$,
and
$\left\lVert\cdot\right\rVert_{H^{1}(\Omega)}=\sqrt{\int_{\Omega}|\cdot|^{2}+|\boldsymbol{\nabla}(\cdot)|^{2}}$.
The integrals in these terms are evaluated with sufficient Gauss quadrature
points to guarantee accuracy. It should be noted that $u^{id}$ can either be a
NN (i.e., $u_{\mathcal{N}}$) or its interpolation (i.e.,
$\tilde{\pi}_{h}(u_{\mathcal{N}})+\bar{u}_{h}$) onto a suitable FE space. In
the different plots in this section, we label the NN interpolation as “FEINN”
and refer to the NN itself with an additional tag “NN only”. From now on, we
drop $\Omega$ in the notation for the norms for conciseness.

The unpreconditioned loss (13) has proven effective in various cases, as shown
in [17]. Therefore, we use it as the default loss function in the experiments.
The preconditioned loss function (12) is only used in Sec. 5.1.2 to examine
the effectiveness of the proposed preconditioning strategy. With regards to
the norm choice in the loss (13), we use the $\ell^{1}$ norm. We note that
other norms can also be explored. For example, our previous work [17] shows
that both $\ell^{2}$ and $\ell^{1}$ norms in (13) produce similar results on
(quasi-)uniform meshes. Besides, we explore different norm options in the
preconditioned loss (12), as detailed in Sec. 5.1.2.

The main motivation of this work is to showcase the capability of the proposed
method in handling PDEs with sharp gradients or singularities. For this
purpose, we have considered standard tests on simple geometries. However, as
demonstrated in our previous work [17], the FEINN method can also be applied
to problems posed on more complex domains, and $h$-adaptive FEINNs inherit
this capability seamlessly. To attack this kind of multi-scale problems, one
can generate a forest-of-octrees built out of a coarse mesh $\mathcal{C}$
generated by an unstructured mesh generator such as, e.g., Gmsh [50]; see Sec.
4.1.

###  5.1. The 2D arc wavefront problem with sharp gradients

((a)) $u$

((b)) Initial mesh

((c)) $h$-Adaptive FEM final mesh

Figure 1. The true solution, initial mesh, and $h$-adaptive FEM final mesh for
the 2D arc wavefront problem. The final mesh is obtained through $h$-adaptive
FEM with 7 iterative adaptation steps using the real FEM error as the error
indicator.

Let us first consider a smooth problem with sharp gradients. Specifically, the
problem is defined on a square domain $\Omega=[0,1]^{2}$. We pick suitable $f$
and $g$ such that the exact solution is:

$u(x,y)=\arctan(100(\sqrt{(x+0.05)^{2}+(y+0.05)^{2}}-0.7)).$

As shown in Fig. 1a, the true solution steeply varies in the neighbourhood of
the arc-shaped front, resulting in sharp gradients along the arc. For the
domain discretisation, we use a forest-of-quadtrees built out of a coarse mesh
$\mathcal{C}$ with a single quadrilateral/adaptive quadtree.

We initiate the training of the NNs with $k_{U}=4$ on the initial mesh
depicted in Fig. 1b. This mesh results from the application of three levels of
uniform mesh refinement to $\mathcal{C}$. For illustration purposes, we solve
the problem using standard $h$-adaptive FEM with 7 mesh adaptation cycles, and
the resulting mesh is shown in Fig. 1c. The error indicator for the
$h$-adaptive FEM used to generate this mesh is the real FEM error.

####  5.1.1. The comparison of error indicators

In the first experiment, we investigate the impact of the choice of error
indicator on the results. We repeat the experiment 20 times, each time with a
different NN initialisation, for every error indicator. The errors versus
adaptation steps for both NNs and their interpolations are illustrated in Fig.
2. The solid line represents the median, and the band illustrates the range
from the 0th to the 90th percentiles across 20 runs. The red dashed line
corresponds to the $h$-adaptive FEM using the real FEM error for mesh
adaptation.

Figure 2. Convergence of $u^{id}$ in $L^{2}$ and $H^{1}$ errors for the 2D arc
wavefront problem, using different error indicators. The top row shows $L^{2}$
errors and the bottom row represents $H^{1}$ errors. The first column displays
errors of the interpolated NNs and the second column corresponds to errors of
the NNs themselves. The line denotes the median, and the band represents the
range from the minimum to the 90th percentile across 20 independent runs. The
red dashed line illustrates the error of the $h$-adaptive FEM solution, using
the real FEM error as the error indicator.

Let us first comment on the results for the interpolated NNs in the first
column of Fig. 2. When using the real error as the mesh adaptation criterion,
both $L^{2}$ and $H^{1}$ errors are consistently the smallest at each
adaptation step. Notably, for $H^{1}$ errors, the FEINN error line aligns with
the one of $h$-adaptive FEM. The other two error indicators are also very
effective for the interpolated NNs, as indicated by their stably decreasing
$L^{2}$ and $H^{1}$ error lines. Besides, the bands of the FEINN errors,
regardless of the error indicator used, are very narrow, underscoring the
robustness of FEINNs with respect to different error indicators and NN
initialisations.

Regarding the errors in NNs themselves, as shown in the second column of Fig.
2, we observe that the majority of the $L^{2}$ errors from the 20 experiments
are within an acceptable range. This is evident as all the lines (medians)
corresponding to different error indicators closely align with those of the
interpolated NNs. Nevertheless, the bands for $L^{2}$ and $H^{1}$ errors
become noticeably wider, suggesting larger variance in the results. The
network error indicator stands out for its robustness, as evidenced by the
narrowest bands in Fig. 2b and 2d. In our previous work [17], we have
highlighted that for a smooth true solution, NNs have the potential to
outperform the FEM solution on the same mesh. We observe a similar phenomenon
here concerning $H^{1}$ errors: with real error indicator, NNs without
interpolation can beat FEM solution (represented by the red dashed line) by
approximately 2 orders of magnitude at the last adaptation step. For the
remaining two error indicators, the majority of the $H^{1}$ errors fall below
the FEINN line, indicating improved performance of the NNs compared to their
interpolation.

In Fig. 3, we show the final meshes corresponding to the three different error
indicators used in this study. We represent, for instance, the mesh adapted
using the Kelly error indicator at adaptation step 7 as $\mathcal{T}_{7}^{k}$.
In this notation, the letter in the superscript indicates the error indicator
(“k” for Kelly, “n” for “network”, and “r” for “real”), while the number
denotes the adaptation step identifier. Notably, all error indicators
effectively identify the region with sharp gradients and appropriately refine
the mesh in that specific area. Nevertheless, differences also exist in these
meshes. Firstly, $\mathcal{T}_{7}^{k}$ (Fig. 3a) with 1,336 elements and
20,297 DoFs has the fewest elements, accounting for the largest errors
identified for Kelly in Fig. 2 compared to other error indicators. This
highlights that, in this problem, the network error indicator outperforms the
Kelly error indicator. Subsequently, $\mathcal{T}_{7}^{r}$ (Fig. 3c) with
2,146 elements and 32,893 DoFs closely resembles the $h$-adaptive FEM final
mesh with 2,119 elements and 32,457 DoFs in Fig. 3c, elucidating the
coincidence of the real indicator $H^{1}$ error line and the FEM line in Fig.
2c. Lastly, although $\mathcal{T}_{7}^{n}$ (Fig. 3b) with 2,293 elements and
35,089 DoFs is as dense as $\mathcal{T}_{7}^{r}$ (Fig. 3c), the patterns are
slightly different: $\mathcal{T}_{7}^{n}$ contains more elements in the middle
of the arc but fewer elements at two ends of the arc, while the opposite holds
true for $\mathcal{T}_{7}^{r}$. This discrepancy could explain the differences
observed in the convergence curves (Fig. 2) for real and network error
indicators.

((a)) $\mathcal{T}_{7}^{k}$

((b)) $\mathcal{T}_{7}^{n}$

((c)) $\mathcal{T}_{7}^{r}$

Figure 3. Comparison of final meshes obtained by training of the $h$-adaptive
FEINNs for the 2D arc wavefront problem. The meshes result from 7 mesh
adaptation steps using Kelly, network, and real error indicators,
respectively.

####  5.1.2. The effect of preconditioning

In the experiments so far, we have observed a notable variance with respect to
parameter initialisation in the errors corresponding to non-interpolated NNs,
particularly when using Kelly and real error indicators. In this experiment,
we investigate whether preconditioning can effectively enhance robustness
(i.e., reduce this variance) and accelerate the training process. We
concentrate on the impact of the preconditioner on the results using the Kelly
error indicator, as accessing to the true error is often impossible in
practice. Besides, we also note that the results obtained with the network
error indicator are very similar to those obtained with the Kelly error
indicator, thus we omit from this section the detailed results obtained with
the former indicator for brevity.

We consider the preconditioner $\mathbf{B}_{\rm inv\\_lin}$, which is defined
as the discrete Laplacian resulting from discretisation with the linearised
space $V_{h}$ used both as trial and test FE spaces. Thus, $\mathbf{B}_{\rm
inv\\_lin}$ is a sparse symmetric positive definite (SPD) matrix with reduced
computational cost compared to a regular sparse matrix. Note that we apply
$k_{U}$ levels of uniform refinement for the partition associated to $U_{h}$
to obtain the mesh for the _linearised_ test space $V_{h}$. Besides, in [17],
we have observed the effectiveness of a GMG preconditioner as an approximation
to the inverse of $\mathbf{B}_{\rm inv\\_lin}$, suggesting its potential
benefits for this problem as well. Nonetheless, for the purpose of this study,
we just exactly invert $\mathbf{B}_{\rm inv\\_lin}$ (by solving linear
systems), as it is sufficient to demonstrate the impact of preconditioning. In
these experiments, we cut the number of training iterations by reducing those
used in the unpreconditioned experiments by a factor of 4, except for the
final adaptation step, where these are divided by a factor of 2. This
adjustment accommodates for the increased demand for training on the finest
meshes.

In addition to evaluating to what extent preconditioners can expedite
convergence, we also investigate the impact of norm choice on the convergence.
We train the same NN (i.e., the same parameter initialisation) using the
preconditioned loss (12) equipped with 4 different norms: the $L^{1}$ and
$L^{2}$ norm of the discrete Riesz map of the residual and its gradient; we
denote the $L^{1}$ and $L^{2}$ norms on the gradient as $W^{1,1}$ and
$W^{1,2}$ norms, since $W^{k,p}$ represents the Sobolev norm with $k$
derivatives and $p$-th power. In Fig. 4, we report the $L^{2}$ and $H^{1}$
errors as a function of training iteration for the four norms. For reference,
we also include the results for the unpreconditioned case, labelled as “N/A”.
The first thing to note is that all the preconditioned curves are much steeper
than the unpreconditioned ones, indicating that the application of the
preconditioner significantly accelerates the convergence of both NNs and their
interpolations. Moreover, the $L^{2}$ convergence for both the NN and its
interpolation is highly unstable when a preconditioner is not applied, as
evidenced by the error curves’ fluctuating pattern. In contrast, the
preconditioned error curves are all steadily decreasing, suggesting that the
preconditioner effectively stabilises the training process. This also
reinforces the conclusion from our previous work [17] that using the $L^{2}$
norm in the preconditioned loss can accelerate the convergence of training in
terms of $L^{2}$ error. Furthermore, the $H^{1}$ convergence of the
interpolated NNs rapidly decreases in both unpreconditioned and preconditioned
cases, as indicated by the stair-like pattern of the convergence curves.
However, upon comparing Fig. 4c and Fig. 4d, it becomes apparent that the
convergence of NN and its interpolation are not at the same pace. Additional
iterations are necessary to ensure the convergence of the NN, even if its
interpolation has already reached a favourable position. To sum up, when using
the preconditioner $\mathbf{B}_{\rm inv\\_lin}$, all the four norms in (12)
are effective, with the $W^{1,1}$ norm being slightly better in terms of
convergence speed and stability.

Figure 4. Convergence history of $u^{id}$ in $L^{2}$ and $H^{1}$ errors for
the 2D arc wavefront problem, using different norms in the preconditioned loss
during training. Kelly error indicator is used for mesh adaptation. The top
row shows $L^{2}$ errors and the bottom row represents $H^{1}$ errors. The
first column displays errors of the interpolated NNs and the second column
corresponds to errors of the NNs themselves. Figure 5. Convergence of $u^{id}$
in $L^{2}$ and $H^{1}$ errors for the 2D arc wavefront problem, using
different preconditioners during training. Kelly error indicator is used for
mesh adaptation. Refer to the caption of Fig. 2 for details on the information
being displayed in this figure.

We now study how preconditioning impacts training robustness with respect to
NN initialisation. We use the same 20 NN initial parameters as in Sec. 5.1.1
and rerun the experiment using the preconditioned loss (12) equipped with the
preconditioner $\mathbf{B}_{\rm inv\\_lin}$ and the $W^{1,1}$ norm. The error
convergence using the preconditioned and unpreconditioned loss functions
versus mesh adaptation step is depicted in Fig. 5. The legend label “none”
means no preconditioning, while “inv_lin” represents that the preconditioner
at hand is applied. We can extract several key findings from the figure.
Firstly, the application of $\mathbf{B}_{\rm inv\\_lin}$ notably enhances the
performance of both the NNs and their FE interpolations. This is evident from
the consistently lower position of the preconditioned lines compared to the
unpreconditioned lines in Fig. 5. Although the preconditioned FEINN $H^{1}$
results may not surpass their unpreconditioned counterparts, they exhibit
close similarity, suggesting convergence to the same solution in both cases.
Secondly, the preconditioner $\mathbf{B}_{\rm inv\\_lin}$ significantly
improves the training robustness, evident in thinner bands for the
preconditioned results, especially notable for the NNs without interpolation
(second column in Fig. 5). Lastly, $\mathbf{B}_{\rm inv\\_lin}$ also helps NNs
surpass their interpolations in terms of $L^{2}$ results, suggested by the
lower position of the preconditioned line for NN $L^{2}$ errors in Fig. 5b
compared to the FEINN $L^{2}$ errors in Fig. 5a.

In summary, the proposed preconditioning strategy proves very effective for
$h$-adaptive FEINNs. It expedites the training process and enhances the
robustness of the results with respect to NN initialisations.

###  5.2. The 2D Fichera problem with singularity

In this subsection, we consider the Fichera problem defined on the L-shaped
domain $\Omega=[-1,1]^{2}\backslash[-1,0]^{2}$. For the discretisation, we use
a forest-of-quadtrees built out of a coarse mesh $\mathcal{C}$ with three
quadrilateral/adaptive quadtrees. We pick $f$ and $g$ such that the true
solution is:

$u(r,\theta)=r^{\frac{2}{3}}\sin(\frac{2}{3}(\theta+\frac{\pi}{2})),$

expressed in polar coordinates. The true solution and the domain are displayed
in Fig. 6a. Note that derivatives of $u$ are singular at the origin, in
particular, $u\in H^{5/3-\epsilon}(\Omega)$ for any $\epsilon>0$. We have
tackled the same problem on the unit square in [17], noting that while the
interpolated NNs can recover the FEM solution on a fixed mesh, there is still
room for their generalisation to be improved. Therefore, our goal in this
experiment is to assess whether the proposed adaptive training strategy
improves the NNs’ generalisation.

((a)) $u$

((b)) Initial mesh

((c)) $h$-Adaptive FEM final mesh

Figure 6. The true solution, initial mesh, and $h$-adaptive FEM final mesh for
the 2D problem with singularity defined on the L-shaped domain. The final mesh
is obtained through the $h$-adaptive FEM with 7 mesh adaptation steps using
the real FEM error as the error indicator for illustration purposes.

####  5.2.1. The comparison of error indicators

Let us first examine the effectiveness of all three error indicators. Starting
with an initial mesh shown in Fig. 6b (resulting from the application of 4
levels of uniform refinement to the initial coarse mesh $\mathcal{C}$), we
perform 20 independent runs for each indicator, and present the errors versus
the adaptation step in Fig. 7. We choose $k_{U}=4$ in these experiments. The
first observation from the figure is that the performance of all three error
indicators is comparable, as the lines representing different indicators
closely overlap. Besides, for interpolated NNs, the distinction among
different error indicators is visually indistinguishable. The $L^{2}$ error
curves (Fig. 7a) and $H^{1}$ error curves (Fig. 7c) for all indicators closely
follow the red dashed $h$-adaptive FEM line. The variance in errors at each
adaptation step is also very low, as indicated by very narrow bands.
Additionally, the NNs exhibit strong generalisation in terms of $L^{2}$ error,
with all lines in Fig. 7c closely aligning with the FEM line and featuring
very narrow bands. More importantly, the $h$-adaptive FEINN method enables NNs
to gradually learn the singularity and achieve the $h$-adaptive FEM rate of
convergence, as shown in Fig. 7d. Furthermore, the variance in the $H^{1}$
errors for NNs is minimal, as suggested by the narrow bands. In brief,
compared to the FEINN method defined on a fixed mesh in our previous work
[17], the $h$-adaptive FEINNs method improves the performance of the NNs
themselves and enable them to effectively capture the singularity,
irrespective of the error indicator employed.

Figure 7. Convergence of $u^{id}$ in $L^{2}$ and $H^{1}$ errors for the 2D
Fichera problem with singularity , using different error indicators. See
caption of Fig. 2 for details on the information being displayed in this
figure.

In Fig. 8, we show the final meshes obtained with different error indicators.
Once again, we notice the striking similarity between $\mathcal{T}_{7}^{r}$
(Fig. 8c) with 2,610 elements and 40,841 DoFs and the $h$-adaptive FEM final
mesh (Fig. 6c) with 2,637 elements and 41,313 DoFs. Although the mesh patterns
of $\mathcal{T}_{7}^{k}$ (Fig. 8a) with 2,637 elements and 40,705 DoFs and
$\mathcal{T}_{7}^{n}$ (Fig. 8b) with 2,703 elements and 41,781 DoFs are
different from $\mathcal{T}_{7}^{r}$, the elements around the singular point
are refined to the same level across all meshes. This is clearer when
comparing the zoomed view of the corner for the different meshes in these
figures. Besides, the differences in the number of elements and DoFs are
negligible, suggesting that training cost is similar for all error indicators.
In the context of this problem and experiment, with high $k_{U}$ and finer
initial mesh, the primary contribution to the $H^{1}$ errors arises from the
singularity. Consequently, despite the differences in mesh patterns, the same
element size around the singular point leads to very similar $H^{1}$ errors.
This observation is reinforced by the remarkably similar $H^{1}$ convergence
curves for different error indicators shown in Fig. 7d.

((a)) $\mathcal{T}_{7}^{k}$

((b)) $\mathcal{T}_{7}^{n}$

((c)) $\mathcal{T}_{7}^{r}$

Figure 8. Comparison of final meshes obtained by training of the $h$-adaptive
FEINNs for the 2D Fichera problem with singularity. The meshes result from 7
mesh adaptation steps using Kelly, network, and real error indicators,
respectively.

####  5.2.2. The effect of FE space interpolation order

So far, we have only conducted experiments using high interpolation order,
i.e., $k_{U}=4$. We now extend our exploration to lower polynomial orders in
solving the same problem, assessing the impact of $k_{U}$ on training results,
while comparing the error decay of $h$-adaptive FEINNs with FEM on uniformly
refined meshes. We experiment using the network error indicator and $k_{U}=2$.
To ensure a fair comparison, the DoFs of the initial trial FE spaces are kept
constant across different values of $k_{U}$. This results in a finer mesh for
$k_{U}=2$ compared to $k_{U}=4$. Besides, from our preliminary experiments, we
find that a lower $k_{U}$ requires fewer training iterations. Therefore, for
$k_{U}=2$, we determine the training iterations to be half of the iterations
for $k_{U}=4$ (specified in Tab. 1). It is important to point out that we
exclude $k_{U}=1$ experiments, as the approximation of a smooth NN with a
linear FE space is ineffective, thus making such experiments not very
meaningful.

Fig. 9 compares the errors of the NNs and their FEINNs versus number of DoFs
for various $k_{U}$. We use different colours to distinguish between orders in
the figure. It is evident that $h$-adaptive FEINNs can more effectively reduce
errors compared to uniform mesh refinement FEM, as all the solid lines
eventually fall below the dashed lines representing the convergences of
uniform FEM. Besides, we observe that the choice of $k_{U}$ has a limited
impact on the convergence order. Across different $k_{U}$, the FEM lines
exhibit almost identical slope, with only small distances between them.
Similarly, the lines representing different $k_{U}$ in both FEINN and NN
results are also very close to each other. In addition, for both interpolation
orders, the $L^{2}$ and $H^{1}$ convergences of both FEINN and NN are
satisfactory. It is worth noting that the results in Fig. 9 are reproducible
with regards to NN initialisations and error indicators. Numerous independent
runs with different NN initial parameters and error indicators have validated
the consistency of these results. We skip other similar results in the text
for conciseness.

Figure 9. Convergence of $u^{id}$ in $L^{2}$ and $H^{1}$ errors for the 2D
Fichera problem with singularity, using different interpolation orders.
Network error indicator is used for mesh adaptation. The dashed lines depict
the errors of the FEM solution with uniform refinement, obtained by using
different polynomial orders. See caption of Fig. 2 for other details in this
figure.

###  5.3. The 3D problem with singularities

((a)) $u$

((b)) Initial mesh

((c)) $h$-Adaptive FEM final mesh

Figure 10. The true solution, initial mesh, and $h$-adaptive FEM final mesh
for the 3D problem featuring singularities. The final mesh is obtained through
the $h$-adaptive FEM with 4 mesh adaptation steps using the real FEM error as
the error indicator.

In our final experiment, we showcase the capability of our method in
addressing 3D problems by considering a PDE defined on the unit cube
$\Omega=[0,1]^{3}$. We consider an initial mesh $\mathcal{C}$ of $3^{3}$
cells. We pick $f$ and $g$ such that the true solution (Fig. 10a) is:

$u(x,y,z)=\sqrt[3]{x^{2}+y^{2}+z^{2}}\sin(2\arcsin(xyz)+\pi/3).$

Note that the derivative of $u$ has two singularities: one at the origin,
stemming from the $\sqrt[3]{x^{2}+y^{2}+z^{2}}$ term, and the other at the
point $[1,1,1]$, coming from the $\arcsin(xyz)$ term. Thus, the final mesh is
expected to be denser around these two singular points, as shown in Fig. 10c.
We use $k_{U}=4$ and start with the initial mesh shown in Fig. 10b, performing
4 mesh adaptation steps. We note that the initial mesh is obtained from the
application of 1 level of uniform mesh refinement to $\mathcal{C}$. For each
error indicator, we conduct 10 independent runs with different initial NN
parameters.

Figure 11. Convergence of $u^{id}$ in $L^{2}$ and $H^{1}$ errors for the 3D
problem featuring singularities, using different error indicators. For each
error indicator, 10 independent runs are performed. See caption of Fig. 2 for
details on the information being displayed in this figure.

The error convergence of both FEINNs and NNs, as presented in Fig. 11,
illustrate an important finding: all three error indicators are equally
effective in solving this 3D problem. Notably, in the case of FEINNs, both
$L^{2}$ errors (Fig. 11a) and $H^{1}$ errors (Fig. 11c) are so closely aligned
with the red dashed FEM lines that they appear visually indistinguishable.
This highlights the exceptional accuracy of FEINNs in mirroring the
$h$-adaptive FEM solution, with negligible impact from NN initialisations.
Furthermore, the generalisation capability of NNs is convincing, with both
$L^{2}$ errors (Fig. 11b) and $H^{1}$ errors (Fig. 11d) closely tracking the
FEM lines. While there is observable variance in these errors, the associated
bands fall within an acceptable range and are not excessively wide. In short,
the results confirm the ability of $h$-adaptive FEINNs in effectively
attacking PDEs with singularities in 3D.

((a)) $\mathcal{T}_{4}^{k}$

((b)) $\mathcal{T}_{4}^{n}$

((c)) $\mathcal{T}_{4}^{r}$

Figure 12. Comparison of final meshes obtained by training of the $h$-adaptive
FEINNs for the 3D problem with singularities. The meshes result from 4 mesh
adaptation steps using Kelly, network, and real error indicators,
respectively.

Fig. 12 displays the final meshes obtained using different error indicators.
Notably, all three error indicators effectively identify the singular points
and refine the elements around them, despite slight variations in the mesh
patterns. In particular, $\mathcal{T}_{4}^{k}$ (Fig. 12a) with 1,819 elements
and 102,587 DoFs and $\mathcal{T}_{4}^{r}$ (Fig. 12c) with 1,819 elements and
101,835 DoFs look very similar to the $h$-adaptive FEM final mesh (Fig. 10c)
with 1,826 elements and 102,811 DoFs, while $\mathcal{T}_{4}^{n}$ (Fig. 12b)
with 2,050 elements and 109,771 DoFs tends to refine the elements along the
edges. Relatively, the number of elements and DoFs for these three meshes are
very close, indicating similar computational costs for training with different
error indicators. With a high interpolation order ($k_{U}=4$) and a fine
initial mesh, the specific patterns of the resulting meshes are less critical,
as long as the elements are refined around singular points. It is noteworthy
that, even in the existence of two singular points, the network error
indicator demonstrates comparable performance to the well-known Kelly error
indicator.

##  6\. Conclusions

In this work, we propose the $h$-adaptive FEINN method, an extension of the
FEINN method [17] that incorporates adaptive meshes. The method inherits the
desired properties of FEINNs, namely the applicability to general geometries,
simple imposition of Dirichlet boundary conditions and accurate integration of
the interpolated weak residual. Furthermore, since the interpolation of the NN
is performed on an adaptive FE space, the method preserves the non-linear
capabilities of the NN to capture sharp gradients and singularities. The dual
norm of the finite element residual in the loss function enhances training
robustness and accelerates convergence. We propose a train, estimate, mark and
adapt strategy to train the FEINN and adapt the mesh simultaneously. We have
carried out a detailed numerical analysis of the method, proving a priori
error estimates depending on the expressiveness of the neural network compared
to the interpolation mesh.

Our numerical findings confirm the capability of $h$-adaptive FEINNs in
capturing sharp gradients and singularities in 2D and 3D cases. Remarkably, in
the 2D arc wavefront experiments, NN exhibit excellent generalisation,
surpassing FEM solutions by two orders of magnitude in terms of $H^{1}$ error
at the final mesh adaptation step. Additionally, our results verify the
substantial improvements in training robustness and convergence with the
introduction of preconditioning techniques. The 2D L-shaped singularity
experiments showcase the superior error reduction capabilities of $h$-adaptive
FEINNs compared to uniform refinement FEM for various interpolation orders.
Furthermore, we demonstrate the method’s capability to capture multiple
singularities in the 3D numerical experiments. Alongside the development, we
introduce a NN error indicator and conduct a thorough comparison with the
standard Kelly error indicator in our experiments. The results show that the
NN error indicator is competitive with the Kelly error indicator in terms of
mesh adaptation efficiency and accuracy.

There are several promising directions for future exploration. One research
avenue is to examine the effectiveness of the adaptive training strategy in
solving inverse problems. In addition, we aim to investigate the applicability
of the proposed approach in attacking problems featuring singularities defined
in $H(\mathrm{curl})$ and $H(\mathrm{div})$ spaces [51]. Towards increasing
applicability, we would also like to explore the combination of $h$-adaptive
FEINNs with unfitted FE spaces [52], which can be used to tackle problems with
complex interfaces and geometries. Future work will also focus on extending
the framework to $p$ and $r$-refinement strategies, as well as to transient
PDEs. Finally, a mathematical analysis of generalisation error bounds for the
FEINN-trained NNs is a promising direction for future research.

##  7\. Acknowledgments

This research was partially funded by the Australian Government through the
Australian Research Council (project numbers DP210103092 and DP220103160).
This work was also supported by computational resources provided by the
Australian Government through NCI under the NCMAS and ANU Merit Allocation
Schemes. W. Li acknowledges the Monash Graduate Scholarship from Monash
University, and the support from the Laboratory for Turbulence Research in
Aerospace and Combustion (LTRAC) at Monash University through the use of their
HPC Clusters.

## References

  * [1] Alexandre Ern and Jean-Luc Guermond  “Finite Elements I”  Springer International Publishing, 2021  DOI: 10.1007/978-3-030-56341-7
  * [2] Douglas N. Arnold, Richard S. Falk and Ragnar Winther  “Finite element exterior calculus, homological techniques, and applications”  In _Acta Numerica_ 15 Cambridge University Press (CUP), 2006, pp. 1–155  DOI: 10.1017/s0962492906210018
  * [3] Peter R. Brune, Matthew G. Knepley, Barry F. Smith and Xuemin Tu  “Composing Scalable Nonlinear Algebraic Solvers”  In _SIAM Review_ 57.4 Society for Industrial & Applied Mathematics (SIAM), 2015, pp. 535–565  DOI: 10.1137/130936725
  * [4] Santiago Badia, Alberto F. Martín and Javier Principe  “Multilevel Balancing Domain Decomposition at Extreme Scales”  In _SIAM Journal on Scientific Computing_ 38.1 Society for Industrial & Applied Mathematics (SIAM), 2016, pp. C22–C52  DOI: 10.1137/15m1013511
  * [5] D. Drzisga, B. Gmeiner, U. Rüde, R. Scheichl and B. Wohlmuth  “Scheduling Massively Parallel Multigrid for Multilevel Monte Carlo Methods”  In _SIAM Journal on Scientific Computing_ 39.5 Society for Industrial & Applied Mathematics (SIAM), 2017, pp. S873–S897  DOI: 10.1137/16m1083591
  * [6] Mark Ainsworth and J.Tinsley Oden  “A posteriori error estimation in finite element analysis”  In _Computer Methods in Applied Mechanics and Engineering_ 142.1, 1997, pp. 1–88  DOI: https://doi.org/10.1016/S0045-7825(96)01107-3
  * [7] Joachim Schöberl  “A Posteriori Error Estimates for Maxwell Equations”  In _Mathematics of Computation_ 77.262 American Mathematical Society, 2008, pp. 633–649  URL: http://www.jstor.org/stable/40234527
  * [8] Olaf Ronneberger, Philipp Fischer and Thomas Brox  “U-Net: Convolutional Networks for Biomedical Image Segmentation”  In _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_ Cham: Springer International Publishing, 2015, pp. 234–241  DOI: https://doi.org/10.1007/978-3-319-24574-4_28
  * [9] Ronan Collobert et al.  “Natural language processing (almost) from scratch”  In _Journal of machine learning research_ 12.ARTICLE, 2011, pp. 2493–2537
  * [10] Alex Graves, Abdel-rahman Mohamed and Geoffrey Hinton  “Speech recognition with deep recurrent neural networks”  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2013, pp. 6645–6649  DOI: 10.1109/ICASSP.2013.6638947
  * [11] M. Raissi, P. Perdikaris and G.E. Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational Physics_ 378, 2019, pp. 686–707  DOI: https://doi.org/10.1016/j.jcp.2018.10.045
  * [12] Guofei Pang, Lu Lu and George Em Karniadakis  “fPINNs: Fractional Physics-Informed Neural Networks”  In _SIAM Journal on Scientific Computing_ 41.4, 2019, pp. A2603–A2626  DOI: https://doi.org/10.1137/18M1229845
  * [13] Dongkun Zhang, Lu Lu, Ling Guo and George Em Karniadakis  “Quantifying total uncertainty in physics-informed neural networks for solving forward and inverse stochastic problems”  In _Journal of Computational Physics_ 397, 2019, pp. 108850  DOI: https://doi.org/10.1016/j.jcp.2019.07.048
  * [14] Ehsan Kharazmi, Zhongqiang Zhang and George E.M. Karniadakis  “hp-VPINNs: Variational physics-informed neural networks with domain decomposition”  In _Computer Methods in Applied Mechanics and Engineering_ 374 Elsevier BV, 2021, pp. 113547  DOI: 10.1016/j.cma.2020.113547
  * [15] Stefano Berrone, Claudio Canuto and Moreno Pintore  “Variational Physics Informed Neural Networks: the Role of Quadratures and Test Functions”  In _Journal of Scientific Computing_ 92.3, 2022, pp. 100  DOI: 10.1007/s10915-022-01950-4
  * [16] N. Sukumar and Ankit Srivastava  “Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks”  In _Computer Methods in Applied Mechanics and Engineering_ 389 Elsevier BV, 2022, pp. 114333  DOI: 10.1016/j.cma.2021.114333
  * [17] Santiago Badia, Wei Li and Alberto F. Martín  “Finite element interpolated neural networks for solving forward and inverse problems”  In _Computer Methods in Applied Mechanics and Engineering_ 418, 2024, pp. 116505  DOI: https://doi.org/10.1016/j.cma.2023.116505
  * [18] Sergio Rojas, Paweł Maczuga, Judit Muñoz-Matute, David Pardo and Maciej Paszynski  “Robust Variational Physics-Informed Neural Networks”, 2023  arXiv:2308.16910 [math.NA]
  * [19] Sifan Wang, Xinling Yu and Paris Perdikaris  “When and why PINNs fail to train: A neural tangent kernel perspective”  In _Journal of Computational Physics_ 449, 2022, pp. 110768  DOI: https://doi.org/10.1016/j.jcp.2021.110768
  * [20] Zhiqiang Cai, Jingshuang Chen and Min Liu  “Self-adaptive deep neural network: Numerical approximation to functions and PDEs”  In _Journal of Computational Physics_ 455, 2022, pp. 111021  DOI: https://doi.org/10.1016/j.jcp.2022.111021
  * [21] Ameya D. Jagtap, Kenji Kawaguchi and George Em Karniadakis  “Adaptive activation functions accelerate convergence in deep and physics-informed neural networks”  In _Journal of Computational Physics_ 404, 2020, pp. 109136  DOI: https://doi.org/10.1016/j.jcp.2019.109136
  * [22] Lu Lu, Xuhui Meng, Zhiping Mao and George Em Karniadakis  “DeepXDE: A Deep Learning Library for Solving Differential Equations”  In _SIAM Review_ 63.1, 2021, pp. 208–228  DOI: 10.1137/19M1274067
  * [23] Kejun Tang, Xiaoliang Wan and Chao Yang  “DAS-PINNs: A deep adaptive sampling method for solving high-dimensional partial differential equations”  In _Journal of Computational Physics_ 476, 2023, pp. 111868  DOI: https://doi.org/10.1016/j.jcp.2022.111868
  * [24] Jie Hou, Ying Li and Shihui Ying  “Enhancing PINNs for solving PDEs via adaptive collocation point movement and adaptive loss weighting”  In _Nonlinear Dynamics_ 111.16, 2023, pp. 15233–15261  DOI: https://doi.org/10.1007/s11071-023-08654-w
  * [25] Zhiping Mao and Xuhui Meng  “Physics-informed neural networks with residual/gradient-based adaptive sampling methods for solving partial differential equations with sharp solutions”  In _Applied Mathematics and Mechanics_ 44.7, 2023, pp. 1069–1084  DOI: 10.1007/s10483-023-2994-7
  * [26] Alexandre Magueresse and Santiago Badia  “Adaptive quadratures for nonlinear approximation of low-dimensional PDEs using smooth neural networks”, 2024  DOI: https://doi.org/10.48550/arXiv.2303.11617
  * [27] D.. Kelly, J.. De S.., O.. Zienkiewicz and I. Babuska  “A posteriori error analysis and adaptive processes in the finite element method: Part I—error analysis”  In _International Journal for Numerical Methods in Engineering_ 19.11, 1983, pp. 1593–1619  DOI: https://doi.org/10.1002/nme.1620191103
  * [28] Piotr Minakowski and Thomas Richter  “A priori and a posteriori error estimates for the Deep Ritz method applied to the Laplace and Stokes problem”  In _Journal of Computational and Applied Mathematics_ 421 Elsevier BV, 2023, pp. 114845  DOI: 10.1016/j.cam.2022.114845
  * [29] Stefano Berrone, Claudio Canuto and Moreno Pintore  “Solving PDEs by variational physics-informed neural networks: an a posteriori error analysis”  In _ANNALI DELL’UNIVERSITA’ DI FERRARA_ 68.2, 2022, pp. 575–595  DOI: 10.1007/s11565-022-00441-6
  * [30] A.. Muzalevskiy and S.. Repin  “A Posteriori Error Control of Approximate Solutions to Boundary Value Problems Found by Neural Networks”  In _Journal of Mathematical Sciences_ 273.4 Springer ScienceBusiness Media LLC, 2023, pp. 492–510  DOI: 10.1007/s10958-023-06516-9
  * [31] Jamie M. Taylor, Manuela Bastidas, Victor M. Calo and David Pardo  “Adaptive Deep Fourier Residual method via overlapping domain decomposition”, 2024  arXiv:2401.04663 [math.NA]
  * [32] D.F. Griffiths and J. Lorenz  “An analysis of the petrov—galerkin finite element method”  In _Computer Methods in Applied Mechanics and Engineering_ 14.1, 1978, pp. 39–64  DOI: https://doi.org/10.1016/0045-7825(78)90012-9
  * [33] Yeonjong Shin, Zhongqiang Zhang and George Em Karniadakis  “Error estimates of residual minimization using neural networks for linear PDEs”  In _Journal of Machine Learning for Modeling and Computing_ 4.4 Begell House, 2023, pp. 73–101  DOI: 10.1615/jmachlearnmodelcomput.2023050411
  * [34] Ignacio Brevis, Ignacio Muga and Kristoffer G. Zee  “Neural control of discrete weak formulations: Galerkin, least squares and minimal-residual methods with quasi-optimal weights”  In _Computer Methods in Applied Mechanics and Engineering_ 402 Elsevier BV, 2022, pp. 115716  DOI: 10.1016/j.cma.2022.115716
  * [35] Philipp Petersen, Mones Raslan and Felix Voigtlaender  “Topological Properties of the Set of Functions Generated by Neural Networks of Fixed Size”  In _Foundations of Computational Mathematics_ 21.2 Springer ScienceBusiness Media LLC, 2020, pp. 375–444  DOI: 10.1007/s10208-020-09461-0
  * [36] Kent-Andre Mardal and Ragnar Winther  “Preconditioning discretizations of systems of partial differential equations”  In _Numerical Linear Algebra with Applications_ 18.1 Wiley, 2010, pp. 1–40  DOI: 10.1002/nla.716
  * [37] Ramon Codina, Santiago Badia, Joan Baiges and Javier Principe  “Variational Multiscale Methods in Computational Fluid Dynamics”  In _Encyclopedia of Computational Mechanics Second Edition_ Wiley, 2017, pp. 1–28  DOI: 10.1002/9781119176817.ecm2117
  * [38] Santiago Badia  “On stabilized finite element methods based on the Scott–Zhang projector. Circumventing the inf–sup condition for the Stokes problem”  In _Computer Methods in Applied Mechanics and Engineering_ 247–248 Elsevier BV, 2012, pp. 65–72  DOI: 10.1016/j.cma.2012.07.020
  * [39] Carsten Burstedde, Lucas C. Wilcox and Omar Ghattas  “p4est: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees”  In _SIAM Journal on Scientific Computing_ 33.3, 2011, pp. 1103–1133  DOI: 10.1137/100791634
  * [40] Santiago Badia, Alberto F. Martín, Eric Neiva and Francesc Verdugo  “A Generic Finite Element Framework on Parallel Tree-Based Adaptive Meshes”  In _SIAM Journal on Scientific Computing_ 42.6, 2020, pp. C436–C468  DOI: 10.1137/20M1328786
  * [41] Alberto F. Martín  “GridapP4est.jl”, 2024  URL: https://github.com/gridap/GridapP4est.jl
  * [42] Marc Olm, Santiago Badia and Alberto F. Martín  “On a general implementation of h- and p-adaptive curl-conforming finite elements”  In _Advances in Engineering Software_ 132 Elsevier BV, 2019, pp. 74–91  DOI: 10.1016/j.advengsoft.2019.03.006
  * [43] Santiago Badia and Francesc Verdugo  “Gridap: An extensible Finite Element toolbox in Julia”  In _Journal of Open Source Software_ 5.52 The Open Journal, 2020, pp. 2520  DOI: 10.21105/joss.02520
  * [44] Francesc Verdugo and Santiago Badia  “The software design of Gridap: A Finite Element package based on the Julia JIT compiler”  In _Computer Physics Communications_ 276 Elsevier BV, 2022, pp. 108341  DOI: 10.1016/j.cpc.2022.108341
  * [45] Michael Innes et al.  “Fashionable Modelling with Flux”  In _CoRR_ abs/1811.01457, 2018  arXiv: https://arxiv.org/abs/1811.01457
  * [46] Mike Innes  “Flux: Elegant Machine Learning with Julia”  In _Journal of Open Source Software_ , 2018  DOI: 10.21105/joss.00602
  * [47] Michael Innes  “Don’t Unroll Adjoint: Differentiating SSA-Form Programs”  In _CoRR_ abs/1810.07951, 2018  arXiv: http://arxiv.org/abs/1810.07951
  * [48] Xavier Glorot and Yoshua Bengio  “Understanding the difficulty of training deep feedforward neural networks”  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ 9, Proceedings of Machine Learning Research  Chia Laguna Resort, Sardinia, Italy: PMLR, 2010, pp. 249–256  URL: https://proceedings.mlr.press/v9/glorot10a.html
  * [49] Patrick Kofod Mogensen and Asbjørn Nilsen Riseth  “Optim: A mathematical optimization package for Julia”  In _Journal of Open Source Software_ 3.24, 2018, pp. 615  DOI: 10.21105/joss.00615
  * [50] Christophe Geuzaine and Jean-François Remacle  “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities”  In _International Journal for Numerical Methods in Engineering_ 79.11, 2009, pp. 1309–1331  DOI: https://doi.org/10.1002/nme.2579
  * [51] Marc Olm, Santiago Badia and Alberto F. Martín  “On a general implementation of h- and p-adaptive curl-conforming finite elements”  In _Advances in Engineering Software_ 132, 2019, pp. 74–91  DOI: https://doi.org/10.1016/j.advengsoft.2019.03.006
  * [52] Santiago Badia, Alberto F. Martín, Eric Neiva and Francesc Verdugo  “The Aggregated Unfitted Finite Element Method on Parallel Tree-Based Adaptive Meshes”  In _SIAM Journal on Scientific Computing_ 43.3 Society for Industrial & Applied Mathematics (SIAM), 2021, pp. C203–C234  DOI: 10.1137/20m1344512

  *[FEM]: finite element method
  *[PDE]: partial differential equation
  *[FEs]: finite element
  *[NNs]: neural network
  *[PDEs]: partial differential equation
  *[PINNs]: physics-informed   ()
  *[PINN]: physics-informed   ()
  *[NN]: neural network
  *[VPINNs]: variational   ()
  *[]:
  *[IVPINN]: interpolated   ()
  *[FE]: finite element
  *[FEINN]:   () interpolated
  *[DFR]: deep Fourier residual
  *[FEINNs]:   () interpolated
  *[GMG]: geometric multigrid
  *[DoFs]: degree of freedom
  *[RHS]: right-hand side
  *[DoF]: degree of freedom
  *[SPD]: symmetric positive definite