G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Authors: James Rowbottom
Georg Maierhofer
Teo Deveney
Eike Hermann Müller
Alberto Paganini
Katharina Schratz
Pietro Lio
Carola-Bibiane Schönlieb
Chris Budd
Venue: ICML 2025 spotlightposter (2025)
Keywords: Mesh Adaptation
r-Adaptivity
Monge–Ampere
Deep Learning
Finite Element Methods
Partial Differential Equations
Neural Network
Graph Neural Networks
Primary area: deep_learning->graph_neural_networks
TL;DR: A graph neural network approach to mesh relocation in finite element methods that directly minimises solution error, achieving higher accuracy compared to classical and prior ML-based r-adaptivity techniques.

Abstract: We present a novel, and effective, approach to achieve optimal mesh relocation in finite element methods (FEMs). The cost and accuracy of FEMs is critically dependent on the choice of mesh points. Mesh relocation (r-adaptivity) seeks to optimise the mesh geometry to obtain the best solution accuracy at given computational budget. Classical r-adaptivity relies on the solution of a separate nonlinear ``meshing'' PDE to determine mesh point locations. This incurs significant cost at remeshing, and relies on estimates that relate interpolation- and FEM-error. Recent machine learning approaches have focused on the construction of fast surrogates for such classical methods. Instead, our new approach trains a graph neural network (GNN) to determine mesh point locations by directly minimising the FE solution error from the PDE system Firedrake to achieve higher solution accuracy. Our GNN architecture closely aligns the mesh solution space to that of classical meshing methodologies, thus replacing classical estimates for optimality with a learnable strategy. This allows for rapid and robust training and results in an extremely efficient and effective GNN approach to online r-adaptivity. Our method outperforms both classical, and prior ML, approaches to r-adaptive meshing. In particular, it achieves lower FE solution error, whilst retaining the significant speed-up over classical methods observed in prior ML work.

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
James Rowbottom *1 Georg Maierhofer
Teo Deveney? 23 Eike Mueller Alberto Paganini*
Katharina Schratz Pietro Liô
Carola-Bibiane Schönlieb Chris Budd
Abstract
1. Introduction
We present novel, and effective, approach to achieve optimal mesh relocation in finite element methods (FEMs). The cost and accuracy of FEMs is critically dependent on the choice of mesh points. Mesh relocation (r-adaptivity) seeks to optimise the mesh geometry obtain the best solution accuracy given computational budget. Classical r-adaptivity relies on the solution of a separate nonlinear "meshing" PDE to determine mesh point locations. This incurs significant cost at remeshing, and relies on estimates that relate interpolation- and FEM-error. Recent machine learning approaches have focused the construc- tion of fast surrogates for such classical methods. Instead, our new approach trains graph neural network (GNN) to determine mesh point locations by directly minimising the FE solution error from the PDE system Firedrake to achieve higher so- lution accuracy. Our GNN architecture closely aligns the mesh solution space to that of classical meshing methodologies, thus replacing classical estimates for optimality with learnable strategy. This allows for rapid and robust training and re- sults in an extremely efficient and effective GNN approach to online r-adaptivity. Our method out- performs both classical, and prior MI., approaches to r-adaptive meshing. In particular, it achieves lower FE solution error, whilst retaining the sig- nificant speed-up over classical methods observed in prior ML. work.
Finite element methods (FEM) are currently the most widely-used tool for the large scale solution of partial differ- ential equations (PDEs) (Ainsworth & Oden, 1997; Cotter, 2023). Central advantages are robustness, reliable error estimates, and thoroughly developed code bases (such as deal.II (Africa et al., 2024), DUNE (Bastian et al., 2008), Fenics (Logg et al., 2012), and Firedrake (Ham et al., 2023)), which are highly parallelisable and efficient. However, even with such optimised software the simulation of large scale problems (c.g. weather forecasting, structural simulations in engineering systems) is computationally costly. An im- portant ingredient that determines the cost is the number of degrees of freedom (DOFs) required by a FEM to satisfy chosen error tolerance. Since this cost depends on the number N₂ of mesh points (x of the underlying compu- tational mesh, it is desirable to keep N₂ moderate. Mesh adaptivity based on mesh refinement and/or relocation to capture important solution features at the right scale can bal- ance computational cost and aceuracy. However, classical mesh-adaptive methods can be difficult to implement and require significant computational resources, In contrast, in the present work, we introduce a cheap, stable and highly efficient graph neural network (GNN) architecture to imple- ment learnable mesh relocation (r-adaptivity). Our method keeps N₂ fixed and adapts the mesh point locations to reduce the overall FE error. Many classical mesh relocation meth- ods have focused finding and minimising mathematical substitutes of the FE error (usually simplified upper bounds) and solving additional (differential) equations to relocate the mesh points. For example, one can solve the Monge- Ampère (MA) equation (Budd et al., 2009) to find mesh point locations that minimise the interpolation error, which is an upper bound (up to some parameters and constants) of the FEM error arising from Céa's lemma (Huang & Russell, 2011). Recent Machine Learning (ML) approaches rely on similar mathematical simplifications and learn surro- gate for the mesh equations (Song et al., 2022; Zhang et al., 2024) leading to significant speed-up with comparable error reduction. In the present work we take an entirely different approach. We present G-adaptivity, an approach to mesh adaptivity that trains a GNN to generate meshes that directly minimise the error of the corresponding FEM solution. We
Equal contribution 'Department of Applied Mathematics and Theoretical Physics, University of Cambridge, UK Department of Computer Science, University of Bath, UK 'Department of Mathematical Sciences, University of Bath, UK *School of Computing and Mathematical Sciences, University of Leices- ter, UK 'Laboratoire Jacques-Louis Lions, Sorbonne Université, France 'Department of Computer Science and Technology, Uni- versity of Cambridge, UK. Correspondence to: James Rowbottom <jr908@ ca0o@@cam..c. auk..>, Georg Maierhofer <gam37 cam.ac.aac.mmmmmmmmmmmmmmm ac.uk>..
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

G-Adaptivity; optimised graph-based mesh relocation for finite element methods
couple backpropagation through novel diffusion-based GNN-deformer, with mesh point gradients obtained through an application of Firedrake adjoint (Mitusch et al., 2019; Ham et al., 2023), to minimise the FEM approximation error directly (as opposed to the upper bound considered in (Huang & Russell, 201 The result is a model capa- ble outperforming the current state-of-the-art r-adaptive methods, whilst retaining the significant acceleration ML based approaches (ef. Figure 1).
Thorough numerical comparison with classical and recent approaches in terms of accuracy, mesh quality and computational time. Our experiments include both stationary and time-dependent test cases.
2. Related work
The effeetive approximation of PDE solutions is one of the central problems in computational mathematics. Over the recent decade, extensive work has been devoted to using ML. for the numerical approximation of PDEs. This includes physics informed neural networks (PINNS) (Raissi et al., 2019; Raissi, 2018), Fourier Neural Operators (FNOs) (Li et., 2020b; 2023), graph neural operators (Li et al., 2020a), DeepONets (Lu et al., 2021), Message Passing Neural PDE Solvers (Brandstetter al., 2022) and the deep Ritz method (E & Yu, 2018). The majority of such approaches try to directly approximate the PDE, or the associated solution operator, with machine learning surrogate. Such methods offer certain advantages (for example in high dimensional settings (Han et al., 2018)), but are typically outperformed by traditional numerical methods in accuracy in most set- tings (Grossmann al., 2023). Our approach is different. We use MI central ingredient of finite element dis- cretisation to construct improved computational mesh, which is then coupled to a classical PDE solver. The crucial advantage is that we retain convergence guarantees and ro- bustmess of FEMs, something that is often lacking in direct ML-based PDE approximations. At the same time our ap- proach achieves a significant speed up in the calculation of the improved mesh compared to classical approaches.
Solution values
MA: 74% ER in 5183ms
DirectOpt: 42.96% ER in 9911ms
G-Adlpt: 39.99% ER in 120ms
Adaptive mesh methods are widely used tool for improv- ing the performance of a classical FEM by varying the local density of the mesh points. This necessary the PDE solution has small length scales or singularities. Adaptivity allows achieving high accuracy without resorting to uniform mesh refinement. The most popular form h-adaptivity (Ainsworth & Oden, 1997), in which mesh cells are subdi- vided when an a-posteriori estimate of the solution error is large. Such methods have complex data structures (see e.g. (Burstedde et al., 201 and, possibly, poor mesh regularity. Alternatively, the relocation based r-adaptive methods con- sidered in this paper move fixed number of mesh points to achieve high density of points where a monitor mz() of the solution error large. Done correctly this lead significant error reduction but at some extra cost (Huang & Russell, 2011).
Figure 1. Optimised meshes from our new approach (G-Adapt) on the example of Poissons' equation: the eтror reduction (ER) achieved by classical Monge-Ampère (MA) can significantly improved with direct optimisation (DirectOpt) of the FEM loss with respect to the mesh points, but prohibitive additional cost. Our new approach achieves near optimal meshes in fraction of the inference time.
Contributions Our work improves earlier MI. based ap- proaches to mesh relocation in the following ways:
A novel training mechanism capitalising recent ad- vances in FEM systems, which leads to fast meshing algorithm that reduces the FE error even over state- of-the-art classical meshing methods. This was not possible in any prior surrogate ML approach;
GNNs are the dominant approach to applying machine learn- ing to irregularly structured data (Bronstein et al., 2017: Battaglia et al., 2018). There has been proliferation of architectures inspired by spectral graph theory (Defferrard et al., 2016), convolutional (GCN) (Kipf & Welling, 2022), message passing (MPNN) (Gilmer al., 2017) and atten-
An improved GNN architecture based on a diffusion deformer, which allows for improved mesh relocation quality and provable non-tangling of generated meshes;
A novel equidistribution loss regularizer, which en- forces mesh regularity in unsupervised GNN training:
2

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
tional (GAT) (Veličković et al., 2018) approaches. More recently range of differential equation inspired architec- tures (Chamberlain et al., 20;;b;a; Giovanni et al., 2023) apply analytical tools to solve known problems with GNNs including stability, over smoothing and bottleneck phenom- ena. This algorithmic aligmment along with the powerful message passing paradigm provide new solutions to some of the most pressing problems in science, including protein folding (Jumper al., 2021), weather prediction (Lam et al., 2023), dynamics learning (Pfaff et al., 2023) and new nu- merical PDE solvers (Brandstetter et al., 2022; Lieen & Günnemann, 2022; Alet et al., 2019).
meshing modifies the location of the node coordinates. Of- ten, and in this work, r-adaptivity is particularly concerned with the reduction of the squared L2-error
E(Z,U2) == 112 ul (2
(3)
For transient problems, we tacitly assume that (3) is evalu- ated at the final time
T.
3.2. Adaptive Meshing
Relocation based r-adaptivity turns a mesh with a certain topology into another mesh with the same topology. For this, the mesh points z are moved, but their connectivity (and hence the associated data structures) is unaltered. Such methods typically map a fixed mesh in a computational do- main (i.e. a representation of the mesh graph) to a deformed mesh in the physical domain where the PDE is posed. Note that, for transient problems, Q should be replaced by ft and d should be replaced by d - 1 in the following explanation (cf. Section 3.1). We denote the mesh points in the physical domain by Z = (z( Ji=i ) which form a triangulation of Ω with Ny mesh elements, i.c.
Fast ML based methods reduce the significant computa- tional cost of classical methods for adaptive meshing. This includes work on h-adaptive mesh refinement (Foucart et al., 2023; Freymuth et al., 2023), and many contributions to adaptivity based on surrogate ML solvers of classical mesh movement PDEs (Yang et al., 2023; Hu et al., 2024) and supervised learning for mesh adaptivity using Graph Neu- ral Networks (GNNs) (Song et al., 2022). A notable re- cent development is the universal mesh movement network (UM2N) (Zhang et al., 2024), which achieves error reduc- tion on par with MA, but at significant speed up, and can also be applied to multiply connected domains.
T ( c 2: 22
IT NT.
We define the following domains and coordinates: the "computational" domain Qe is mapped to the "physical" domain 22p C Idd and the "computational" coordinates
3. Preliminaries and background
& E EBC are mapped to the "physical" coordinates # E Qpp. To construct an adaptive mesh we consider differentiable, possibly time-dependent, deformation map F Qg + S2p. so that F(E.t) t) and F(ACc) = 22pp. If g are the fixed mesh points in the computational domain then 2(0) F(g(1) t). Assuming the mesh in the computational domain is regular, then determining the (properties of the) mesh in the physical domain, reduces to finding, (and analysing), F.
3.1. Problem specification
We consider finite element solutions to nonlinear second- order PDEs on general domains Q of dimension d. In ab- stract form, we can write these PDEs as follows:
Flu,, u, Vu, V2u)= f in 2, a0 = g on 02
(1)
where a, 3 E R. For transient problems, l2 = ă x (0, T),
Location based methods find F by solving a PDE, or a linked variational principle. Monge-Ampére (MA) methods assume that F is a Legendre transform with 'mesh po- tential' Φ{E t) for which F = Ve The linearisation of F is given by J = OF/06- = (()) where H is the Hessian of Φ. Relocation methods usually equidistribute a monitor function m(2) so that e satisfies the MA equation
where 9 is a (d - 1)-dimensional spatial domain and (0, T)
is the time-interval of interest. In this case, we employ the method of lines and combine the FEM with suitable time- stepping schemes (Hairer& & Wanner, 1996). To compute finite element solutions, we introduce mesh T of the spa- tial domain with N₂ nodes, which we collect in the node set Z. The mesh T is used to construct trial and test func- tions with local support to discretize (1). For example, for a Poisson problem with homogeneous Dirichlet boundary conditions we consider the space of piecewise linear func- tions (vanishing on 02) Sz on T and solve: Find Uz E Sz such that
m(z)]H(() (V (VV))(())))))))))))())))))))))() = =0, for constant
0.
(4)
For example, in (Huang & Russell, 2011) m(z) is an a- priori monitor of the interpolation error. The PDE (4) has a unique, convex, solution (Budd et al., 2013) which avoids mesh tangling. However, (4) is expensive solve and, in its pure form, only applicable to simply connected domains. Solution procedures include relaxation methods (Budd et al., 2009), quasi- measi-Newton methods (McRae et al., 2018), surro- gates (Song et al., 2022; Zhang et al., 2024), and PINNs (Yang et al., 2023).
VO) = (J, (J))) Ve E S2;
(2)
where (c))) denotes the inner-product in (?))..
To minimise the error E(Z, U2) between the exact solution of (1) and its finite element approximation Uz, rdadppivve
3

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Velocity based methody find an ODE describing the mesh point evolution in pseudo-time that
Integral the G-adaptivity framework is the construction of the feature matrix such that the GNN can act as mesh de- former. Similar to (Zhang et al., 2024) we construct the fea- ture matrix by concatenating coordinates of regular mesh & E RRN xd with a learnable feature encoding h (2. H) which, motivated by (4), is dependent on the Frobenius norm of the Hessian H (U22)) - 1|0,0,u u | $2 + of the FEM solution Uza on the undeformed mesh 2 When higher order finite element functions are used in the approx- imation space of U2 the Hessian can be obtained by simple differentiation, but even in the case of linear elements this in- formation is recoverable using widely-used techniques such as the one described in Appendix A.2. The final input fea- ture matrix is then X = (2 || xX) E IR x(d+ Xy. = he (z), H),, where À denotes the index set for the node fea- tures, which is passed into a GNN mesh deformer that then outputs the relocated mesh points. Previous works (Song et al., 2022; Zhang et al., 2024) used a graph attention net- work (GAT) (Veličković et al., 2018) as the GNN mesh deformer
0z 10r = v0z=1(2(1), (t)
(5)
The choice of velocity function is critical to the success of such methods, and is often motivated by natural Lagrangian structures of the underlying PDE. These methods provide the basis of our diffusion-based deformer (diffformer) de- scribed in section 4.1 and while they often lead to mesh tangling where mesh lines cross (cf. Ch. in (Huang & Russell, 1)), our architecture is specifically designed to enforce non-tangling of the mesh (cf. Theorem 4.2).
4. The G-adaptivity framework
The G-adaptive mesh relocation method described below essentially a velocity based method with learnable coef- ficients that are trained by calculating the rate of change of the FE solution error E with respect to the mesh point location. As we explain below, the G-adaptivity framework combines feature selection, structural regularization and direct optimisation to learn optimal mesh relocation in an unsupervised manner whilst avoiding mesh tangling.
Zk+1 x++1
A (x = a (Ax
(6)
where X = (2 Xx) and WA is a learnable linear transfor- mation matrix. To prevent mesh crossing the non-linearity and channel mixing are excluded from the positional chan- nels in (6). In the above Ao (X+) is row-stochastic meaning that the top row of (6) corresponds to graph-bases aver- aging over graph neighbours. Motivated by (Chamberlain et al., 2021b) and velocity- based methods for meshpoint re- location introduced in section 3.2, in our G-Adaptive frame- work this average is replaced by a diffusion based deformer (henceforth referred to as Diffformer)
4.1. Graph-based adaptive mesh refinement
For simplicity of exposition we focus our discussion on the 2D case, but note that this approach generalises in a straight- forward manner to 3D cases as shown in Section 5.5. A mesh T (i.c. a triangulation of the domain 22 with mesh- points Z gives rise to a natural graph, with the nodeset V Z and the edgeset & ((z',z) E V x V; E TT s.t. z',2/ E A), i.e. two nodes share an edge if there is triangle in the mesh T which has both nodes as ver- tices. The graph (V,E) can be enriched with node fea- tures {x, E Rd : V} represented in matrix notation as X E X.d. For example we could associate to each mesh point (0)) the value of the solution field u() as feature. Likewise, we can introduce latent features that propagate through repeated application of a map the graph, this is used in our architecture (cf. Figure 2). Mesh connectivity is stored in the adjacency matrix A (where ag = if (i, E E and zero otherwise). A graph neu- ral network (GNN) Ma : N. E >> N.dx is a map from features features constructed with layers Lo, N.dx N acting node wise as
E(T = (A+(X*) 1)2(7), 2(0) 2²
(7)
which is solved to finite end time Tend and leads to the meshpoint update zh+1 = B(Tend), i.e. an overall deformer of the form
zk+1 (x = Oz (ax(xx(X)))
Z(Tend)
(8)
As before, the learnable attention Ag row-stochastic, meaning (7) is essentially a diffusion equation on the graph V. Further details on our Diffformer are provided in Ap- pendix A. We can stack multiple layers of (8), each time varying the number of hidden feature dimensions which are updated using the second row of (8). We denote the overall GNN by the map Ma and a schematic overview of the components of Ma is provided in Figure 2.
x k+1 x+1 LLo = doo
Σ spoe (xx) JEN;
x
4.2. Structural regularization
where Po. is the learnable edge-wise operation, do, is a learnable node-wise aggregation and N,, = {j E V; (i.j) E
The architectural changes between (6) and (8) lead to several regularity properties that we refer to as structural regulariza- tion. The Diffformer based architecture has a key advantage
E} is the set of nodes adjacent to the meshpoint i.
4

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Curvature information: H(
Input mesh coordinates )
Zm (5m 1=1
Feature extractor X hoy Zin; H)
Computation of attention - Ax
Figure 3. The action of the graph diffusion pulls nodes into the convex hull of their graph neighbours.
Diffformer:
Feature deformer:
8 = -
x Jx((*)XX
classical relocation methods.
Next layer: k k
Deformed mesh coordinates Bon o0s/im1 1
4.3. Firedrake adjoint optimal gradient computation and direct FEM loss
Training the GNN Ma requires computing the derivative of F(2) E Uz) with respect to the node coordinates Z. Since evaluating E'z, Uz) requires solving the PDE
Figure 2. Schematic overview of our new graph diffusion-based architecture.
(1) first, a naive application of automated differentiation would result in the solution of additional (Ng x d) PDEs (N; x (d- 1) in the transient case). To avoid the additional computational cost, we employ the well-established method of adjoints. Specifically, adjoint capabilities implemented in pyadjoint (Mitusch et al., 2019). With pyadjoint, derivatives with respect to mesh coordinates can be computed in an automated fashion as shape derivatives of E Uz) ) in directions discretized with vector-valued linear Lagrangian FEs (Ham et al., 2019a).
over other velocity based methods in the generation of reg- ular meshes. A requirement of FEM meshes is that they are not 'tangled', i.e. that they form a well-posed triangu- lation the domain 52 (i.e. no triangles overlap). This follows if each mesh point in the interior of the convex hull of its neighbours on the graph and can equivalently be characterised using the Jacobian of the deformation map M.
employ Firedrake's automated
Definition 4.1. Let (+) be the Jacobian of the deformation map M at simplex (6), A mesh is said to be tangled if there exists simplex where the determinant of the Jacobian (J)) < 0 (Huang & Russell, 2011).
Remark 4.3. Often, the exact solution to (1) is not known and must itself be approximated with the FEM (e.g.. by interpolating onto Uz a FE solution computed on much finer mesh). In this case, essential correct the direc- tional derivatives computed with Firedrake by adding the corrections terms stemming from evaluating the formula
Velocity based methods often lead to tangled meshes due to the local way in which the mesh point movement is defined. However, this does not arise in our method.
1 f(u Uz) dr
Theorem 4.2 Eiscrete-Time Non-Tangling). If the diffu- sion equation (7) is solved with the forward Euler method, then for sufficiently small pseudo-timestep dr < 1/2, the discrete mesh evolution under the deformation map M pre- serves element orientations, ensuring that no mesh tangling occurs.
along each finite element direction V E St (V S s(d-1) in the transient case). This is necessary because the shape derivative of a FE function a direction discretized with FEs is zero (Ham et al., 2019a).
full proof is given in Appendix F but in essence the dif- fusion process ensures that meshpoints are simultaneously moved along directions that point into the convex hull of neighbouring meshpoints, thus ensuring that tangling cannot occur (cf. Figure 3).
4.4. Regularized gradients
In line with the concept of equidistribution discussed in Appendix E.2 we introduce a regularizing term train- ing which further enforces mesh regularity in an unsuper- vised manner and leads to improved training of the mesh deformer using only information about predefined mon- itor function m( (we follow (Zhang et al., 2024) and pro- vide a global signal that moves mesh points into regions of the domain where the solution varies and likely requires
The proof relies on the softmax of the attention mechanism normalising the adjacency to be row stochastic and for the time step of the residual connection to be controllable. This a benefit over (6) and allows (8) to learn an anisotropic diffusion which is akin a learnable monitor function from
use m(Uz) = 1 maxp
[VU2] ( The motivation
5

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
more meshpoints to resolve. For this we add the following regularizing term to our loss:
For each experiment described below we build a randomised dataset with held-out test data. This means in each case we specify set of solution values through varying source terms or boundary conditions (adjusted to the PDE at hand). We then generate training and test sets of coarse, undeformed, meshes T of varying resolution with associated FEM solution values and, for each mesh, we also compute a reference solution (on finer reference mesh) which serves as comparison for the error computation. For most test cases the solution values are generated by randomly sampling Gaussians in the domain Ω. For the example of the Navier- Stokes flow around a cylinder we used snapshots of a time- series simulation of vortex shedding. Full details on the specific configurations for each experiment are provided in Appendix C.
Datasets
Lequi (2) = Σ
IS
JAn m(x)da m
2
AGIET
where m- = 171- Σ SaU m(x)dr. Given the area ADTET
of simplex in the mesh a() E T the terms in the above loss are approximated by SAU) m(z)dx x a( (a) E EZEaU) m(z), where C2 = 1/3, cz 1/4. This leads to the following full regularized loss which we use in the training of our Diffformer:
Co- E(Mo) + (Mo))(Mo)
During training the weighted graph Laplacian (Ao - in (7)) will adaptively adjust to minimize both terms, meaning the mesh evolution will not purely follow the degree weighted graph Laplacian dynamics but will now be biased towards error reduction and equidistribution.
Baselines We compare our algorithm against two adap- tive mesh algorithms: classical MA as described and im- plemented by (Wallwork et al., 2024) and the ML based surrogate GNN method UM2N (Zhang et al., 2024). These two state-of-the-art approaches serve as baseline for the FEM error reduction and deformation time. As a third base- line we train UM2N on our regularized PDE loss (9), which we denote by UM2N-G in the tables below, in order to high- light the performance improvements gained from both our new architecture (Figure 2) and our new training (9).
5. Experimental results
We evaluate G-adaptivity on three classical meshing prob- lems in two-dimensions: an elliptic PDE (Poisson's 's equa- tion) in a variety of convex domains, a nonlinear time- evolution PDE (Burger's equation), and the time dependent Stwier-Stokes equations in a multiply connected domain. In the following we present the performance improvements obtained in terms of the FEM T.2-erron reduction (ef. (3)) and compute time, using our novel approach for adaptive meshing on each of these problems. Additional experiments and sensitivity analysis is provide in the Appendix D. Full code to build the datasets and reproduce our results can be found at https://github..o http://githbb.com//owwbotoot wwtp///wwwbotooniggggaaaaaaoeeee
Experiment details Here we refer to our framework G- adaptivity and our model Diffformer synonymously, which we train using the regularized PDE-loss (9) (which is the regularized L FEM upproximation emror + equidistribution regularizer). Calculation of the 1.2-error is obtained by cal- culating an FE solution on the moved mesh and comparing this to the projection of the reference solution of a fine regu- lar reference mesh onto sufficiently higher order elements. We train our new model (G-Adapt) and UM2N-G for 300 epochs using an Adam optimiser and learning rate of 0.001, Our model has 4 diffformer blocks as described in Appendix C. Each blocked is rolled out using explicit Euler integration for 32 steps with a step size of 0. For the baseline UM2N we trained using 1000 epochs in order to achieve good per- formance but we believe further tuning of the training may be required in order to achieve a similar performance to the one reported in (Zhang et al., 2024). UM2N remains important baseline and we expect that with appropriate training the method would be able to achieve similar error reduction (ER) as the Monge-Ampee (MA) solver, but we would like to highlight that even in the best reported results of the original paper, UM2N never achieved a larger error reduction than MA.
5.1. Experimental details
Method Our experimental pipeline consists of three parts: (i) we build datasets containing information about the PDE, FEM solution on a regular (i.e. not relocated) grid and the corresponding approximation of the Hessian of the solu- tion; (ii) then we train either our model or the baseline to predict a relocated mesh on which we perform another FE solve to obtain the improved solution approximation; (iii) finally, we compare this FEM solution to a reference solu- tion (calculated on a fine reference mesh) and determine the change in -error over the original undeformed mesh, i.c. in the above notation we look at relative error reduction of E(UM,) over E(U These steps are repeated for the three PDE datasets as described below. Note that the Fire- drake adjoint solve is only required during training of our G-Adaptive network and not during inference. The times reported in the below numerical examples thus contain a true reflection of the fast online mesh adaption times.
Evaluation We report three metrics to evaluate the per- formance of the mesh relocation methods at hand: (i) the relative 1,2 error reduction (ER) of the FE solution on the
6

---

G-Adaptivity; optimised graph-based mesh relocation for finite element methods
relocated mesh versus the FE solution on the initial coarse mesh (larger error reduction means improved performance); (ii) the time taken to relocate the mesh (shorter times means faster relocation); (iii) the aspect ratio of the deformed mesh as a measure of mesh quality as described Appendix F.4 (a single digit aspect ratio is generally acceptable and a smaller aspect ratio indicates a more regular mesh). Each experiment performed in full five times (training and eval- uation) with different random seeds to provide the error bars.
data which is available during inference, thus it does constitute practical adaptive meshing strategy. In contrast, once trained, our G-Adaptive approach yields fast online mesh movement without needing reference solution values.
5.2. Benchmarking on Poisson's equation
Our first benchmark is the classical Poisson problem -Vau fz) with Dirichlet boundary conditions. Full details of the FEM formulation are given in Appendix C.1. We benchmark the results against two datasets on a square (cf. Figure 1) and polygonal domain respectively (ef. Figure 4). Further evaluations on five additional (non- convex) geometries are provided in Appendix D.4. For both aforementioned examples sample source terms and boundary conditions corresponding underlying Gaussian fields. On the square domain Ω [0, 1j2 we initialise the mesh-deformation with a regular grid and to showcase G-adaptivty' s ability to work on irregular domains with un- structured meshes we apply a similar methodology to the convex polygonal dataset (a sample is shown Figure 4). The results for both datasets are presented in Table 9, The central observation in these results is that our methodol- ogy provides the very first MI. approach to mesh relocation which is able to outperform MA in terms of eтror reduction, while retaining the fast mesh relocation times given by the state-of-the-an GNN UM2N (Zhang et al., 2024).
Figure 4. The solution fields and corresponding G-Adaptive mesh with 23.10 % error reduction in 59ms polygonal domain.
5.3. l'mee-eepndent Burgers' equation
In second example we highlight that our approach can equally well be applied to deneneennd problems, in particular the viscous Burgers' equation:
du at
+ (u V)u vVu
Further details on the specific FEM implementation (and implicit time-stepper) used given in Appendix C.2. We randomly sample Gaussians the square domain Ω [0, 1|2 as initial conditions for the evolution in Burgers' equation and perform the following experiments.
Table 1. Benchmarking results on Poisson Square and Poisson Convex Polygon datasets.
Burgers' square rollout: We train the models on a set Gaussian initial conditions for a timestep = 0.02 with steps and evaluate by following
trajectories of randomly sampled Gaussians in the Burgers equation for 20 timesteps, remeshing after every 2 steps (ef. Figure 5 and Appendix D.3). The results in the top part of Table 2 show the average error reduction over achieved over every
POISSON SQUARE
MODEL
ERROR RED. (%)
TIME (MS)
ASPECT
DIRECTOPT
27.40 + 0.00
126028
33.99 +
MA
12.69 + 0.00
3780
2.11 + 0.00
block of two timesteps. While the MA performs well on this task, we note that the UM2N and UM2N-G baselines appear to lead a negative error reduction (i.e. an increase), which is likely due the fact that the Burgers' equation changes the solution shape and thus trajectories will lead iittof-distribution cases for methods that are trained only on initial conditions. Due to the structural regularity of our new approach (cf. Section 4.2) our approach is able to deal with cot-of-distribution data very well, and most importantly is able to outperform MA in terms of error reduction while retaining a fast mesh relocation time.
UM2N
6.83 + 1.10
70
1.99 + 0.03
UM2N-G
16.40 +2.65
30
2.61
G-ADAPT
21.01 0.33
2.92 + 0.03
POISSON CONVEX POLYGON
DIRECTOPT
20.51 0.00
56280
1.000
MA
10.97 + 0.00
4446
1.95
UM2N
3.12 + 0.38
36
1.66+0.03
UM2N-G
15,00 + 0.13
16
1.88 0.03
G-ADAPT
16.84 1 0.10
55
1.86 1 0.02
I The direct optimization method is included here purely for exposition, showing that MA-meshes are not necessarily optimal. DirectOpt computes the optimal mesh for given PDE with known solution but is extremely slow and relies on
Burgers' square 10 steps: The interpolation error in remeshing is significant and provides a central limitation to current mesh relocation techniques (cf. (Budd al.,

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
2009)). It is thus desirable relocate meshes only after several timesteps. It turns that our approach lends it- self to targeted training not just of GNN that reduces the FEM error in a stationary sense, but GNN that seeks to find an optimal mesh given a specified remeshing frequency. The classical method MA has no means of inferring this information or adjusting the meshes accordingly. On this example we trained the GNN on collection of random Gaussian initial conditions with the loss attained by solving the corresponding FEM problem for 10 timesteps of size 81 = 0.02. The results in Table 2 highlight that in this way we can achieve even more significant ER over MA thus lead- ing to efficient meshes that require less frequent changes in time-evolving systems.
Table Benchmarking results Burgers' Square datasets.
BURGERS' SQUARE ROLLOUT
MODEL
ERROR RED. (%)
TIME (MS)
ASPECT
MA
25.78 + 0.00
18884
1.99
0,00
UM2N
-11.24 2.52
314
2.52
0.39
UM2N-G
-2.33 + 3.73
316
2.83
0.11
G-ADAPT
27.17 + 0.34
717
2.85 + 0.12
BURGERS' SQUARE STEPS
MA
12.28 0.00
11566
1.99 0.00
UM2N
3.52 + 0.60 0.60
30
2.33 0.01
UM2N-G
16,38 + 2.85
41
1.83 = 0.10
G-ADAPT
21.66 + 3.13
93
2.82 = 0.06
5.4. Navier-Stokes equation and flow past cylinder
Our final example is the canonical flow past cylinder prob- lem we simulate data using an FE solution for 400 time steps of size öt 0.01 of the time series evolution ex- pressed in Gaussian basis function expansions (ef. Figure 6 and Appendix A.3). The training and test data are 25 and 50 respectively random snapshots from the range E [1, 4] with remeshing after every 5 timesteps. Full details of the PDE and FEM formulation are provided in Appendix C.3. Again we observe good error reduction and fast mesh relocation times in our new methodology.
Velocity field t=0.2
Velocity field
Solution velocity field (x-component)
mesh at t-0.2
MA mesh t- 0.6
Undeformed mesh
G-Adaptivity deformed mesh (23.57% ER in 32ms)
G-Adapt. mesh at t=0.
G-Adapt. mesh at t=0.6
Figure 6. The G--Adaptivity-defomaa mesh the Navier-Stokes equation (23.57% error reduction in 32ms). The adapted mesh correctly recognises
Figure Snapshots the velocity field (x-component) together with the corresponding deformed meshes provided by Monge- Ampère (MA) with 46.52% average error reduction over the full solution path compared to the deformed meshes provided by approach (G-Adapt.) with 19.15% error reduction.
of large solution curvature and resolves them more finely (on the upstream side of the cylinder resolving the stagnation point singularity and along the path of shed vorticity).
8

---

G-Adaptivity; optimised graph-based mesh relocation for finite element methods
Our model consistently achieved significant mesh adapta- tion, accuracy improvement, and computational acceleration compared to Monge-Ampiss matching the performance ob- served on smaller-scale experiments.
Table 3. Benchmarking results Navier-Stokes datasets,
NAVIER-STOKES
MODEL
ERROR RED. (%)
TIME (MS)
ASPECT
Solution
MA
MA
NA 1.34 + 0.57
UM2N
44
1.65 = 0.06
00x00 serres
UM2N-G
25.55 + 0.81
30
2.32 2 0.06
G-ADAPT
26.36*1.37
49
3.51 # 0.81
Standard Monge-Ampèrr solvers do not converge on mul- tiply connected domains. i Since MA data available use the best UM2N model from Section 5.2.
150x150 points
5.5. 3D adaptive meshing
The G- G-Aaapiivity framework and diffusion deformer model also easily adapted the 3D setting. To demonstrate this 10x10x unit cube for the 3D Poisson problem with Dirichlet boundary conditions and Gaussian solutions, analogous to Section 5.2. An example of the corresponding results can be seen in Figure 7. In the interest of brevity, the full numerical results are presented in Appendix 5.6 and show that the method outperforms MA significantly (out-of-the-box UM2N does not apply in 3D) and that it leads to effective mesh point concentration in regions of interest.
150x150 soints
we perform an experiment
Figure 8. G-Adaptivity large-scale fine meshes
6. Conclusions and future work
We have presented a novel, and effective, approach the classical problem of r-adaptive meshing for FEM solutions of PDEs. In particular, with differentiable FEM solver (Firedrake), and loss function given by the regularized solution error, can be ef- fectively used to optimise the location mesh points to minimise the FEM error. Hence we can take an entirely different route from prior work (both classical and ML ap- proaches) which determine good choices of mesh points by mualysis-inspires heuristics using location based approach. We demonstrate the advantages of our method on challeng- ing test problems in two dimensions, including in multiply connected domain, and find that, on those examples, we are able to outperform both classical and ML methods in terms of error reduction while retaining similar computational prior MI. work. We note that the direct FEM error optimisation approach extends naturally more complex domains, and PDEs, where classical methods may struggle providing a basis for future extensions of this work.
demonstrate, that GNNs together
Solution field
Rogular mesh
mesh
G-Adapt mesh
Figure 7. Examples of 3D solution fields and adapted meshes.
5.6. Scalability of the G-Adaptivity framework
The G-Adaptivity framework is able to scale to very large meshes. In particular the inductive learning property of GNNs ensures the ability of GNNs to transfer to unseen graphs this case meaning we can perform super-resolution. In Table we report experiments where the model trained on 15x15 mesh and inference is performed larger 60x60 (3,600 nodes) and 150x (22,500 nodes) meshes for the Poisson problem with 128 sampled Gaussians (see Figure 8). In order to scale the transformer encoder, which in naive form scales with O edges we use sliding win- dow SWIN (Liu et al., 2021) style transformer to capture the monitor function embedding the mid-length scales.
Finally, we note that any machine learning- -based approach is inherently statistical in nature, meaning that GNN- based meshing tools are likely to perform worse on out-of- distribution data. We observed this in our experiments with G-Adaptive
both pre- mre-trinee UM2N models and approach when applied to PDEs whose solutions exhibited markedly different scales and features from those seen dur- ing training. Enhancing the scale-gened generalisation capabilities of ML-based adaptive meshing therefore remains an impor- tant open problem for future investigation.

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Acknowledgements
P., Buschelman, K., Dalcin, L.,, Dener, A., Bijkhout, V., Gropp, W., et al. Petsc users manual. 2019.
The authors would like to thank Patrick Farrell and David Ham for helpful advice related to capabilities of Firedrake, as well Joseph G. Wallwork and Mingrui Zhang for shar- ing code and advice the use of the UM2N baseline model and the Python Movement package. JR. GM, TD, CBS & CB gratefully acknowledge support from the EPSRC programme grant in 'The Mathematics of Deep Learning', under the project EP/VO26259/1. GM gratefully acknowl- edges funding from the Mathematical Institute, University of Oxford. KS gratefully acknowledges funding from the Euro- pean Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 850941).
Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Kornubee, R., Ohlberger, M., and Sander, O. A generic grid interface for parallel and adaptive scientific computing. part ii: implementation and tests in dune. Computing, 82:121-138, 2008.
Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-
Gonzalez, A., Zambaldi, V. F., Malinowski, M., Tacchetti,
A., Raposo, D., Santoro, A., Faulkner, R., Gülçehre, Ç., Song, H. F., Ballard, A. J., Gilmer, J., Dahl, G. E., Vaswani, A., Allen, K. R., Nash, C., Langston, V., Dyer, C., Heess, N., Wierstra, D., Kohli, P., Botvinick, M. M., Vinyals, O., Li, Y., and Pascanu, R. Relational induc- tive biases, deep learning, and graph networks, CoRR, abs/1 bbss1806.00161,, 2018.
Impact Statement
This paper presents work that aims accelerate and improve the performance of mesh relocation methods for FEM using Machine Learning. FEMs are omnipresent in scientific computing and the generation and adaption of an effective mesh is paramount to any large-scale application of FEMs. The potential applications and benefits of advancements in this field are thus significant as complex nonlinear PDEs appear everywhere in nature: from weather and climate forecasting over oceanography up to general relativity.
Bouziani, N. and Ham, D. A. Physics-driven machine learn- ing models coupling Py PyTorch and Firedrake. 2023.
Bouziani, N., Ham, D. A., and Farsi, A. Differentiable pro-
gramming across the PDE and Machine Learning barrier. arXiv preprint arXiv:2409. 060850 2024.
Brandstetter, J., Worrall, D., and Welling, M. Message Pass- ing Neural PDE Solvers. arXiv:2202.03376 [cs, math], February 2022.
References
Bronstein, M. M., Bruna, J., LeCun, Y.. Szlam, A., and Van- dergheynst, P. Geometric deep learning: Going beyond euclidean data. IEEE Signal Processing Magazine, 34(4): 18-42, 2017. doi: 10. 10.1109//PP201172999911
Africa, P. C., Arndt, D., Bangerth, W., Blais, B., Fehling, M., Gassmöller, R., Heister, T., Heltai, L.,, Kinnewig, S., Kronbichler, M., Maier, M., Munch, P., Schreter- Fleischhacker, M., Thiele, J. P.. Turcksin, B., Wells, D., and Yushutin, V. The deal.ii library, version 9.6. Journal of Numerical Mathematics, 92(4):369-380 2024. doi: 10. 11155/jjma2222-01133
Budd, C., Cullen, M., and Walsh, E. Ampér- based moving mesh methods for numerical weather prediction, with applications to the Eady problem. Journal of Com- putational Physics, 236:247-270, March 2013. ISSN 00219991. doi: j0.1016/...p..0122.11.
Ainsworth, M. and Oden, J. T. A posteriori error estimation in finite element analysis. Computer Methods in Applied Mechanics and Engineering, 142(1):1-88, March 1997. ISSN 0045-7825. doi: 24.001660045--885599)))000010
Budd, C. J., Huang, W., and Russell, R. D. Adaptivity with moving grids. Acta Numerica, 18:111-241, May 2009. ISSN 1474-0508, 0962-4929. S0962492906400015.
Alet, F., Jeewajee, A. K., Villalonga, M. B., Rodriguez, A., Lozano-Perez, T., and Kaelbling, I.. Graph Element Networks: Adaptive, structured computation and memory. In Proceedings of the 36th International Conference on Machine Learning, pp. 212-222. PMLR, May 2019.
doi: 10.1017/
Burstedde, C., Wilcox, L. C., and Ghattas, O. p4est: Scal- able algorithms for parallel adaptive mesh refinement on forests of octrees. SIAM Journal Scientific Computing, 2011.
Alnæs, M. S., Logg, A.. Ølgaard, K. B., Rognes, M. E., and Wells, G. N. Unified form language: A domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematicul Software (TOMS), 40(2):1-37, 2014.
33(3):
Chamberlain, B., Rowbottom, J., Eynard, D., Di Giovanni, F., Dong, X., and Bronstein, M. Beltrami Flow and Neural Diffusion on Graphs. In Advances in Neural Information Processing Systems, volume 34, pp. 1594-1609 Curran Associates, Inc., 2021a.
Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune,
10

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Chamberlain, B., Rowbottom, J., Gorinova, M. I., Bronstein, M., Webb, S., and Rossi, E. GRAND: Graph Neural Diffusion. In Proceedings of the 38th International Con- ference on Machine Learning, pp. 4418- PMLR, July 2021b.
beat the Finite Element Method? 2023.
arXiv:2302.04107,
Hairer, E. and Wanner, G. Solving Ordinary Differen- tial Equations II, volume 14 of Springer Series in Com- putational Mathematicc. Springer. Berlin, Heidelberg, 1996. ISBN 978-9-642-05220-D 978-3-642-05221-77 doi: 10. 10.1007/978-3-422052222
Chorin, A. J. The numerical solution of the navier-stokes equations for an incompressible fluid. Bulletin of the American Mathematical Society, (3(6):928-931) 1967.
Ham, D. A., Mitchell, L., Paganini, A., and Wechsung, F. Automated shape differentiationn in the unified form language. Structural and Multidisciplinary Op- timization, 60(5):1813-1820, 2019a. doi: 10.1007/
Chorin, J. Numerical solution of the navier-stokes equa- tions. Mathematics of Compulation, 22(104):755-76, 1968.
Cotter, C. Compatible finite element methods for geophysical fluid dynamics. Acta Numerica, 32:291- 393, May 2023. ISSN 0962-4929, 1474-0508. doi: 10.1017/5092249223300000
Ham, D. A., Mitchell, L., Paganini, A., and Wechsung, F. Automated shape differentiation in the Unified Form Language. Structural and multidisciplinary optimization, 60:1813 10:1813-1820 2019b.
Defferrard, M., Bresson, X., and Vandergheynst, P. Convo- lutional Neural Networks on Graphs with Fast Localized Spectral Filtering. In Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
Ham, D. A., Kelly, H. J., Mitchell, L.., Cotter, C. J., Kirby, R. C., Sagiyama, K., Bouziani, N., Vorderwuelbecke, S., Gregory, J., Betteridge, J.. Shapero, D. R., Nixon-Hill, R. W., Ward, C. J., Farrell, P. E., Brubeck, P. D., Mars- den, I., Gibson, T. H. Homolya, M., Sun, T., McRae, A. T. T., Luporini, F., Gregory, A., Lange, M., Funke, S. W., Rathgeber. F., Bercea, G.-T., and Markall, G. R. Firedrake User Manual. Imperial College London and University of Oxford and Baylor University and Univer- sity of Washington, first edition edition, 5 2023.
E, W. and Yu, B. The Deep Ritz Method: Deep Learning-Baae Based Numerical Algorithm for Solving Vari- ational Problems. Communications in Mathematies and Statistics, 1-1) March 2018. ISSN X...-671X. doi: 111007//40304--18-0122
Farrell, E., Ham, D. A., Funke, S. W., and Rognes, M. E. Automated derivation of the adjoint of high-level transient finite element programs. SIAM Journal on Scientific Computing, 35(4):C369-C35 C393, 2013.
Han, J., Jentzen, A., and E, W. Solving high- -ish-dimensional partial differential equations using deep learning. Pro- ceedings of the National Academy of Sciences, 115(34); 8510,, 2018
Foucart, C., Charous, A., and Lermusiaux, F. Deep rein- forcement learning for adaptive mesh refinement. Journal of Computational Physics, 491:1 12381, October 2023. ISSN 0021-9991. doi: 10. 1e..00///.cp.222..1222
Hu, P., Wang. Y., and Ma, Z.-M. Better Neural PDE Ssolvers Through Data-Free Mesh Movers. The Twelfth Interna- tional Conference on Learning Representations, 2024.
Freymuth, N., Dahlinger, P., Würth, T., Reisch, S., Kärger, L., and Neumann, G. Swarm reinforcement learning for adaptive mesh refinement. Advances in Neural Informa- tion Processing Systems, 2023,
Huang, W. and Russell, R. Adaptive Moving Mesh Methods, volume 174. January 2011. ISBN 978-1-4419-7915.5 doi: 10.10077978-11444997166
Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K. Bates, R., Zidek, A., Potapenko, A., Bridgland, A., Meyer, C., Kohl, S. A. A., Ballard, A. J., Cowie, A., Romera-Paredes, B., Nikolov, S., Jain, R., Adler, J., Back. T.. Petersen, S., Reiman, D., Clancy, E., Zielinski, M., Steinegger, M., Pacholska, M., Berghammer, T., Bodenstein, S., Silver, D., Vinyals, O., Senior, A. W., Kavukcuoglu, K., Kohli, P., and Hassabis, D. Highly accurate pro- tein structure prediction with AlphaFold. Nature, 596 (7873):583-58) 5899 August 2021. ISSN 1476-4687. doi: 10. 111/8//41586-022 021-199---
Gilmer, J., Schoenbolz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. Neural Message Passing for Quantum Chem- istry. In Proceedings of the 34th International Conference Machine Learning, pp. 1263-1272. PMLR, July 2017.
Giovanni, D., Rowbottom, J., Chamberlain, B. P., Markovich, T., and Bronstein, M. M. Understanding convolution on graphs via energies. Transactions on Ma- chine Learning Research, June 2023. ISSN 2835-8856.
Grossmann, T. G., Komorowska, U. J., Latz, J., and Schönlieb, C.B.. Can Physics-Informe Neural Networks

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Kipf, T. N. and Welling, M. Semi-Supervised Classification with Graph Convolutional Networks. In International Conference on Learning Representations, July 2022.
Sphere Using Finite Elements. SIAM Journal on Sci- entific Computing, 40(2)A1121-A111 2018.
Mitusch, S., Funke, S., and Dokken, dolfin-adjoint 2018. 1: automated adjoints for fenies and firedrake. Journal of Open Source Software, 4(38):1292, 2019.
Lam, R., Sanchez-Gonzaler A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., Ravuri, S., Ewalds, T., Eaton- Rosen, Z., Hu, W., Merose, A., Hoyer. S., Holland, G., Vinyals, O., Stott, J., Pritzel, A., Mohamed, S., and Battaglia, P. Learning skillful medium-range global weather forecasting. Science, (92(6677):141161144 De- cember 2023. doi: 10. 00.1126/sciemcc.adi23336
Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., and Battaglia, P. Learning Mesh-Based Simulation with Graph Networks. In International Conference on Learn- ing Representations, February 2023.
Li, Z.,, Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., and Anandkumar, A. Multipole Graph Neural Operator for Parametric Partial Differential Equations. In Advances in Neural Information Processing Systems, volume 33, pp. 6755-6766. Curran Associates, Inc., 2020a.
Picasso, M., Alauzet, F,, Borouchaki, H., and George, P.-L. A numerical study of some hessian recovery techniques on isotropic and anisotropic meshes. SIAM Journal on Scientific Computing,
2011.
Piccolo, C. and Cullen, M. A new implementation of the adaptive mesh transform in the Met Office 3D-Var system. Quarterly Journal of the Royal Meteorological Society,
Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Fourier Neural Operator for Parametric Partial Differential Equa- tions. In International Conference on Learning Represen- tations, September 2020b.
138(667) 2012.
Raissi, M. Forward-Backward Stochastic Neural Networks: Deep Learning of High-dimensional Partial Differential Equations, April 2018.
Li, Z., Kovachki, N. B., Choy, C., Li, B., Kossaifi, J., Otta, S. P., Nabian, M. A., Stadler, M., Hundt, C., Azizzade- nesheli, K., and Anandkumar, A. Geometry-Informed Neural Operator for Large-Scale 3D PDEs. In Thirty- Seventh Conference on Neural Information Processing Systems, November 2023.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686 -777 February 2019. ISSN 0021-999 doi: 10. 1.p.006/jjjpp22118 1.0.04..
Lienen, M. and Günnemann, S. Learning the Dynamics of Physical Systems from Sparse Observations with Fi- nite Element Networks. In International Conference Learning Representationss March 2022.
Rathgeber, F, Markall, G. R., Mitchell, L., Loriant, N., Ham, D. A., Bertolli, C., and Kelly, P. H. Pyop2: A high-level framework for performance-portie porttacl- simula- tions on unstructured meshes. In 2012 SC Companion: High Performance Computing, Networking Storage and Analysis, pp. 1116-1123, IEEE, 2012.
Liu, Zz,, Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang. Z., Lin, S., and Guo, B. Swin transformer: Hierarchical vision trans- former using shifted windows. In 2021 IEEE/CVF In- ternational Conference Computer Vision, ICCV 2021, Montreal, QC, Canada, October 10-17, 2021, pp. 9992- 10002. IEEE, 2021.
Song, W., Zhang, M., Wallwork, J. G., Gao, J., Tian, 2., Sun, F., Piggott, M. D., Chen, J., Shi, Z., Chen, X., and Wang, J. M2N: Mesh Movement Networks for PDE Solvers. In Advances in Neural Information Processing Systems, May 2022.
Logg, A., Mardal, K.-A,, and Wells, G. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.
Uribe, D., Durand, C., Baudouin, C., and Bigot, R. Enhanc- ing data representation in forging processes: Investigating discretization and R-adaptivity strategies with Proper Or- thogonal Decomposition reduction. Finite Elements in Analysis and Design, 242: 2024.
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3): 229, March 2021. ISSN 2522-5839. doi: 10. 00.10088/62256--011000000
Vallet, M.-G., Manole, C.-M., Dompierre. J., Dufour, S., and Guibault, F. Numerical comparison of some hessian recovery techniques. International Journal for Numerical Methods in Engineering, 72(8):987-1077, 2007.
McRae, A. T. T., Cotter, C. J., and Budd, C. J. Optimal- Transport-Ban Based Mesh Adaptivity on the Plane and
12

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y. Graph Attention Networks. In International Conference on Learning Representations, February 2018
Wallwork, J. G., Kramer, S. C., Zhang, M., and Dundovic, D. Movement, May 2024.
Yang, Y., Yang, Q., Deng, Y., and He, Q. MMPDE-Net and Moving Sampling Physics-informed Neural Networks Based On Moving Mesh Method, November 2023.
Zhang, M., Wang, C., Kramer, S., Wallwork, J. G., Li, S.,
Liu, J., Chen, X., and Piggott, M. D. Towards universal mesh movement networks. 2024.
13

---

G-Adaptivity; optimised graph-based mesh relocation for finite element methods
A. Implementation details
A.1. Diffusion deformer (diffformer) details
We apply the diffformer in learned blocks
2h
D (1)
(X)
II + ( (8)
(x()) ( - 1)) x
(10)
12=0
where 0, T denotes the discrete time step index, N₂ is the number of blocks in the deformer, such that A (b) (X) the attentional adjacency matrix block b, dynamically learned as:
exp O X x2))
(b) ai
E CENN()) exp( O (X,, Xx))
(b)
Then the full G-adaptivity diffusion based deformation Map given by
No. II 6=0
M)(Xo,A) =
(b)
x (x)
(1I)
The process consists of: 1. Initializing the feature matrix (() as the feature positions. 2. Looping over N, deformer blocks, updating positions iteratively. Applying dd steps of discrete evolution refine the mesh over time.
The input feature matrix is Xo- = (E h) E RBNxxd+1y utilises the graph transformer encoder of (Zhang al., 2024) with the exact same hyperparameters. Similarly each attentional matrix A. is adapted from the same. We use N₂ = blocks and rollout using explicit Euler time integration for 32 timesteps with a step size of 0.1.
A.2. Hessian recovery
To identify parts of the domain Ω where the solution varies rapidly in space, we use an estimator for the local Hessian H H v) which is inspired by the approach in (Picasso al., 2011). For piecewise linear function u E V an approximation of components of is obtained by solving the weak problem
J 0,ud,o dx Hifu dr for all (Q2). vlan
(12)
for Hij subject to the strong Dirichlet boundary condition Hyylan = While there might be other Hessian recovery techniques (see e.g. (Vallet al., 2007)), observe empirically that our approach leads to good results if the Frobenius norm ||||| En. H2 is fed as an input to the GNN.
A.3. Gaussian basis function expansion for time-dynamic training
For technical reasons, in the Navier Stokes dataset it was necessary to provide the initial conditions used for training in analytical form an UFL (Alnæs et al., 2014) expression that can be fed to Firedrake. To achieve this, snapshots of the pressure and velocity fields are taken at specified times during the numerical solution of the time-dependent Navier Stokes equations. The fields w(x, y) are approximated as a sum of Gaussian basis functions in the form
Wear(r,)) = Σ aijo(æ - Fi.y - Uj) with Φ(x, v) exp
(13)
2
h2
where the n, x n,, = x nodal points (r,, Uj) are arranged in regular Cartesian grid over the domain with grid spacings h2 and hy. The expansion coefficients are chosen such that WoBF(2, ()), = w(2,, 81). The sum on the right hand side of (13) can be implemented an UFL expression.
14

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
B. Notes on the use of Firedrake in G-Adaptivity
Firedrake (Ham et al., 2023) is a Python framework for the automatic solution of finite element problems. The central design idea based on composable abstractions, which allow the expression of the partial differential equation in weak form at high level in Unified Form Language (UFL) (Alnæs et al., 2014). This abstraction is gradually lowered to generate C-kernels for matrix-assembly that can be executed in grid traversal with PyOP2 (Rathgeber et al., 2012). PETSc (Balay et al., 2019) provides a wide range of linear- and non-linear solvers for the resulting linear algebra problem. Firedrake supports broad collection of finite element discretisations and adjoint problem for a given forward equation. The recently added interface to PyTorch (Bouziani & Ham, 2023) is crucial for the work in this paper.
(Mitusch et al., 2019) allows the automatic construction of the
B.1. Additional details on implementation
Training the GNN requires computing the derivative of the loss function E(Z,Uz) U2) with respect to node coordinates Z. Since E(Z, U2) is a PDE-constrained functional, it is necessary to use adjoint models to compute these derivatives efficiently. The derivative and adjoint formulas depend on the loss function and its PDE constraints and automating their derivation crucial develop general r-adaptivity methodology that can be trained seamlessly on different test cases. Firedrake is the perfect tool for this because it can derive adjoint models (Farrell et al., 2013; Mitusch et al., 2019) and automatically compute derivatives of EZZ U2) with respect node coordinates (Ham et al., 2019b). Deriving these formulas by hand is nontrivial, tedious, and error prone. Firedrake is fully integrated with PyTTorch (Bouziani et al., 2024), and this is key to formulate hybrid FEM-torch architectures required to train the GNN. Implementing our approach in Firedrake required minimal adaptations: the GNN model must conform to the Firedrake external operator API, and a term must be added to the derivatives with respect to node coordinates when F Uz) comprises a finite element solution computed on a finer mesh.
As an example consider the shape derivative d. (Z, [2) of the functional J Uz) =I1uz (2) which is simplified version of E(Z, U2) in (3). The constraint on Uz is given by the simplest testcase: the weak Poisson equation in Appendix C.1 with f = 4. With Firedrake and Pdjjooin, we can compute dJ (Z, (2] as follows:
mesh = UnitSquareMesh (3, 3)
continue inniinneenntttiion/nnc
Q mesh mesh.coordd nates. function_space
Function (mesh. wowrdinate...ooo..oo............ funtttiosssuntiitsppppcpllllll mesh. wwwrdinates. (
(mesh. coordinates T)
FunctionSpace (mesh, \CG", 1)
= Function (V)
=
TestFunction (V)
solve (dot (grad (u), grad (v)) -4+v) =dx==0, u, bcs=Dirichlets( (v, 0, \on lon_boundary") )
assemble (d*+2*dx)
Jred ReducedFunctionall (J, Control )T))
Jred. derivative
Crucially, this only requires the implementation of the forward constraint equation in Appendix C. On the other hand, tedious manual derivation of (2, [Z] leads (2, Uz) SSa (Uz + VUz Vp - 4 T VUz(DT) + D da with p being the (weak) solution of the adjoint equation Ap- = 2Uz These formulae are problem dependent and will be significantly more complicated for other PDE constraints. For test cases such as the Navier Stokes equations in Appendix C.3 this approach quickly becomes intractable, as highlighted in (Ham al., 2019b, p. 1818).
In contrast, adapting the code above to the problems described in Appendix C.2 & C.3 requires only minor changes.
C. Mathematical description of the numerical experiments
C.1. Poisson's equation
Poisson's 's equation ~v = fz) is solved using the Finite Element method in the two-dimensional convex domain z Ω C R.2
We use the weak formulation (2) and seek piecewise linear functions u E Sz with Dirichlet Boundary
15

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
conditions wlan 0 such that
J
/Vu Vu
Vu dr = / fo dæ for all test functions v E Sz, vlon =
(14)
C.2. Burgers' equation
The non-linear viscous Burgers' equation describes the evolution of the vector- valued velocity field u as
Ōu - - 0 in f2
+ (u
(15)
Öt
where > is the kinematic viscosity and solve consider two-dimensional rectangular domain f R². The term (u V)u describes non-linear convection and v is the viscous diffusion.
We use piecewise linear Finite Element discretisation with E S2-1, A simple backward- Euler timestepping method with step-size At is employed to compute the velocity un+1 E 12-1 at the next timestep from the current velocity u, $2-1 The discretisdd weak form of (15) is given by: find U,+1 E 52-1 such that
J
un++1 (unt)
Un v (Un+1 Vun41) /V :V dx
(16)
A
for all test functions s2-1. The final two terms in (16) are the weak form of the nonlinear advection and viscous diffusion term respectively.
C.3. The Navier-Stokes Equations
We consider the incompressible Navier-Stokes equations in primitive form for time-dependent velocity field u and pressure p in the two-dimensional spatial domain S [0, 2.2] [0, 0.41]
du + (u V)u-vu- Vp f, in 8
(17)
V u = 0, in Ω2
(18)
Here is again the kinematic viscosity and is an external force term.
The Finite Element discretisation uses Taylor-Hood elements with piecewise linear pressure and vector-valued piecewise quadratic velocity functions (u,p) p) Q x Sz. The time- steepit procedure, which computes the velocity u, Q and pressure Pu+1 E $2 the next timestep from the current velocity u, E Q is variant of Chorin's projection method (Chorin, 1967; 1968). It consists of three steps, each of which requires the solution of weak problem.
Step 1: Compute tentative velocity u Find 11"
Q such that:
uu
Un At
vV) dx + / (pan v u(Vumid n) v) ds = / f v
v (u, Vumia) v + VVuiid
Jmn
da. dr. (19)
22
for all piecewise quadratic vector-valued test functions v E Q where Umid = (u, + u*). Homogeneous Dirichlet boundary conditions are applied at the top (y 0.41) and bottom (y 0) of the domain. The velocity field prescribed on the inflow boundary at the left side of the domain as u(z y) (4.0 1.5 0.41-y 0). The weak problem in (19) solved with GMRES iteration that preconditioned with successive overrelaxation (SOR).
Step 2: Solve for pressure correction To ensure that the velocity field at the next timestep is divergence-free, find Pa+1 E Sz which satisfies the following elliptic problem:
Vpn+i Vqdz= / Vpr. Vg 1 (V
(20)
At for all piecewise linear pressure test functions E Sz. To deal with the fact that the pressure is only determined up to an additive constant, homogencous Dirichlet boundary conditions are applied to Px+1, at the outflow boundary. The weak problem in (20) solved with a conjugate gradient iteration preconditioned with algebraic multigrid (AMG).
16

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Step 3: Update velocity
Find u 22- such that:
J Un+1 v dr u" v dx At / V(pn+1 - PPn) v dz.
(21)
for all piecewise quadratic vector-valued test functions v E Q2- The weak problem in (21) is solved with conjugate gradient iteration preconditioned with SOR.
D. Further details of numerical experiments
D.I. Model and data hyperparameters
Table 4 shows for each PDE and geometry the number of train and test samples will as the resolution or node count for the train, test dataset and evaluation mesh.
PDE
Poisson
Burgers
Navier-Stokes
Domain
Square
Polygonal
Square
Cylinder
Train/Test Samples
100/100
100/100
100/100
25/50
Train Resolution
[15x15, 20x20]
1 114 nodes
[15x15, 20x20]
201 nodes
Test Resolution
[12x|2...,23xx3]]
1 114 nodes
[12x12...,33x3]]
201 nodes
Eval Resolution
100x100
228 nodes
100x100
402 nodes
Table 4. Summary of PDE problem setups, including domains, sample sizes, and training/lesinn/eeaaltttinn resolutions.
17

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
D.2. Additional results in the Poisson square and polygon
In Figures 9 and 10 we present additional plots from the Poisson experiments on the square and polygonal domain detailed in the main paper exhibiting the types of meshes generated with novel G-Adaptive methodology.
FEM MA mesh
MA mesh
Mumodel mesh
FCM MA mech
Figure 9. Comparison of MA and ML model-generated meshes for Poisson problems on square and polygonal domains
18

---

G-Adaptivity; optimised graph-based mesh relocation for finite element methods
Mumodell mand
Figure 10. Comparison of MA and ML gener-geeerated meshes for Poisson problems square and polygonal domains
D.3. Additional examples of Burger's evolution
In Figure we include some additional mesh trajectories from the experiment performed Section 5.3.
Weccly held at 0-022
Weccty feld at t-0.4
Velcory Reld tw0.5
Webddt Ne di t-0.0
GuAdape t-0.2
G-00apt mech 1-04
G-Adapt 50
G-Aaatt mesh
Wilcity fiedd
Wecty Aeld 1-0
t-
Adaet. mesh at 1-002
G-Meept mesh 0.4
G-Adept mesh tt.
G-Acast mesh Mt 0.0
Figure Snapshots the velocity field (--comoneett together with corresponding deformed meshes provided by
19

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
D.4. Further experiments on non-convex domains
While section 5.4 already contains an example of non-convex domain we provide further evidence that our method extends to this case using domain data from UM2N (Zhang et al., 2024). In particular conducted experiments similar to the setup of 5.2 on five non-convex domains, cf. Figure 12. On each domain solve Poisson's equation for randomly sampled Gaussian solutions with 100 training datapoints and 100 unseen test datapoints.
The results (error reduction scores are listed in in Table 5 and the full results can be seen in Figure 12) confirm that our method performs robustly on non-convex geometries, achieving significantly greater error reduction than baselines and generating regular non-tangled meshes on all tested domains, succeeding even when some other approaches fail. Note that the UM2N results reported below were obtained using the pretrained model from the UM2N repository, since the MA meshes obtained using (Wallwork et al., 2024) were unsuitable for direct training in these cases.
Table 5. Error reduction scores (%) on Poisson's 's equation in non-convex domain.
DOMAIN
MA'
UM2N!
UM2N-G
G-ADAPT
GEOMETRY
0.23 + 0.00
-76.85 + 0.00
1.92 + 0.02
7.97 = 0.04
GEOMETRY 2
-1.00 0.00
-83.88 0.00
0.69 0.08
8.88 = 0.24
GEOMETRY
-7.82
-.09
6.62 0.09
H-GEOMETRY
108.31 + 0.00
-73.59 + 0.00
-0.92
7.51 0.00
L-GEOMETRY
89.40 + 0.00
-138.43 +: 0.00
13.94 1.18
16.25 + 0.25
Monge-Ampert solvers in general struggle with non-convex domains. t Since the MA data available is not suitable for accurate training we use the pretrained model from (Zhang et al., 2024) for our evaluation.
D.5. Model hyper-parameter sensitivity analysis
We have performed extensive sensitivity studies and found that our approach is robust to the particular choice of hyperpa- rameters for the diffformer blocks. We used N,, blocks and rollout using explicit Euler time integration for 32 timesteps with step size of 0.1. It should be noted that the hyperparameters were identical in all experiments performed in the paper and did not require finetuning to the specific problem. Tables 6 and 7 show the sensitivity of the model to the diffusion parameters terms of error reduction and inference time.
Table 6. The effect of pseudotime and diffusion timesteps the Error reduction.
DIFFUSION TIME AND TIME-STEP PERFORMANCE SENSITIVITY
dr' No. TIMESTEESSSSS
2
4
8
16
32
64
0.05
10.41
16.60
18.95
15.85
21.82
22.93
0.1
12.97
14.52
20.68
15.71
20.27
20.85
0.25
19.54
19.11
22.30
19.94
23.11
22.09
0.5
20.43
22.65
22.14
22.42
21.32
21.92
20.60
21.57
21.16
22.10
19.71
19.40
Table 7. The effect of pseudotime and diffusion timesteps on inference time (ms).
DIFFUSION TIME AND TIME-STEP INFERENCE TIME
dr' -IMESTMESSS
16
32
0.05
60
44
46
116
65
247
0.1
54
42
79
61
86
208
0.25
40
48
119
108
0.5
50
58
49
92
125
158
I
52
59
45
69
100
108
We investigated the sensitivity of G-Adaptivity to the weighting of the equidistribution regularizer in Section 4.4. We found
20

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Solution field
MA mesh
UM2N mesh
UM2N-G mesh
G-Adapt mesh
Geometry 1
Error red:
0.23%
Error red:
-6685%%
Eror red: 0.2%%
Error red: 7.97=0.04%
Time:
2433ms
98ms
Time:
Aspect: 2.43
Aspect:
3.41
Aspect:
1.48+0.00
Aspect:
1.79+0.01
Geometry 2
Error red: -1.00%
Error red:
-83.88%
Error red: 0.69
Error red: 8.88=0.24%
Time: 1648ms
Time:
84ms
Time:
58ms
Time: 154ms
Aspect: 69.55
Aspect:
3.49
Aspect:
1.37=0.00
Aspect:
1.75+0.01
N
Geometry 3
X
Error red:
Erior red:
-75.82%
Error -0.96
Error red: 6.62
Time:
Time:
94ms
Time:
84ms
Time:
171ms
H
Aspect:
Aspect:
4.16
Aspect:
1.39=0.00
Aspect:
1.59=0.00
H-Geometry
Error red: -108. 31%
Enor red:
-73.59%
Error -0.92=0.000
7.51=0.00%
Time: 550ms
Time:
83ms
Time:
59ms
Time:
158ms
Aspect:
Aspect:
7.05
Aspect:
1.36x0.00
Aspect:
1.64:0.00
L-Geometry
13.94=1.18%
16.25=0.25%
Time:
1306ms
Time:
94ms
Time:
88ms
Time:
160ms
Aspect:
90.75
Aspect:
59.43
Aspect:
3.68:0.41
Aspect:
4.26:0.75
Figure 12. Benchmarking results for Poisson's 's equation non-convex geometrics.
an optimal value in the Poisson square example was a loss weighting of 1, which was used throughout all experiments.
21

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Table 8. The effect of equidistribution loss regularisation Error reduction.
EQUIDISTRIBUTION LOSS REGULARIZER PERFORMANCE SENSITIVTY
REG. WEIGHT 0 (NO EQUI.-DIST. LOSS)
0.5
2
4
8
ERROR RED. (%)
22.42
22.95
23.99
23.21
22.14
20.96
D.6. Scalability of the G-Adaptivity framework
In Section 5.6 we claimed the G-Adaptivity framework is able to scale to very large meshes via super-resolution. In Table 9 we report experiments where the model is trained on 15x15 mesh and inference is performed on larger 60x60 (3,600 nodes) and 150x150 (22,500 nodes) meshes for the Poisson problem with 128 sampled Gaussians (see Figure 8). We also show the results for G-Adaptivity applied to 10x (1,000 node) cube (see Figure 7).
Table 9. Benchmarking results scaling G-Adaptivity on Poisson Square Cube datasets.
SCALING G-ADAPTIVITY
SCALE
MODEL
ERROR RED. (%)
TIME (MS)
60X60
MA
11.94 + 1,50
23084
G-ADAPT
27.47 # 0.89
452
150XI50
MA
17.96 =
II5395
G-ADAPT
25.70 1.51
2555
10X10X10
MA
12.71 = 0.00
41049
G-ADAPT
28.08 # 0.36
494
E. Further details on mesh relocation
E.1. Advantages and limitations of r-adaptivity
r-adaptivity is a newer technology than h-adaptivity and as such is not yet widely adapted in industry. However, it has certain significant advantages over h-adaptivityy. In particular works with a constant data structure, is easy to use on parallel architectures, it gives more regular mesh (often with guaranteed mesh regularity), it naturally inherits Lagrangian and scaling structures in PDE (which very useful for example in ocean modelling and studying PDEs with singularities), and can be easily linked to existing external software designed to solve a PDE discontinuous Galerkin solver). As a result, r-adaptive methods have recently been very successfully used, for example, in the operational data assimilation codes of national weather forecasting offices, which when coupled to the computational dynamical core, have led to a very significant increases in computational accuracy, particularly for resolving local weather features such as fog and ice (Piccolo & Cullen, 2012). r-adaptivity has also found natural applications in the steel industry where the Lagrangian nature of the approach is very well suited the evolving fine structures in the forging process (Uribe al., 2024). Possible disadvantages of r-adaptivity, such excessive mesh computation cost, and tendency mesh tangling, are exactly the issues address in this paper, proposing fast and accurate method which avoids tangling.
unstructured mesh (for example
E.2. The Equidistribution Principle
The equidistribution principle applied to mesh with cells C, used for an FE calculation of a function u(z), aims to minimise the total error over all the cells by equidistributing it over each cell. Typically the error over such cell can be measured (or estimated) by the integral of an appropriate monitor function over that cell, or more simply by the expression
mlz)[Cal
(22)
where z is representative point in the cell, and (in the two dimensional case) |C,| is the cell area. The equidistribution condition on the cells C then becomes
(z) = 0,
(23)
where f (to be determined) is a constant. The function m is usually a function of u. An important example is given by the problem of linearly interpolating uzz) as it follows from Céa's lemma that the resulting interpolation error is an (often tight)
22

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
upper bound for the FE solution error. In this case m will be function of the eurvature of (with the exact form dependent on the norm used to measure the error) (Huang & Russell, 2011).
In the context of r-adaptivity each such cell C, in the physical domain, will be the image, under the action of the deformation map F of reference cell (of fixed area) in the computational domain. The area [C], of C will then be proportional to dJ)(J) where J is the Jacobian of F. The equidistribution condition (23) then becomes:
m(z) (J) 0.
(24)
Note that the application of the monitor function in this way is equivalent to defining a measure on the physical space.
In one dimension the equidistribution condition (24) uniquely defines each cell length, and thus the cell shape, and hence the whole mesh. However in two dimensions it only gives the cell area but not the shape. To find the mesh uniquely additional conditions must be imposed. Noting the correspondence between the equidistribution condition and measure on the physical space, the deformation map can be viewed as mapping a uniform measure in the computational space to a new measure in the physical space. is natural to seek a map which minimises the cost of doing this, as this leads to meshes in the physical domain which are close to uniform and hence have minimal skewness and which avoid tangling. This gives an obvious link between mesh generation and optimal transport. In the continuous setting such map can be calculated by solving (either directly or by using surrogate solver) an associated Monge-Ampér equation, leading to the MA methods described in the main body of the text. Note that with modifications this procedure can also be used to generate meshes on non-planar manifolds (McRae et al., 2018).
F. Mesh Tangling Prevention
We provide formal proof that mesh evolution scheme based on the row-stochastic weighted graph Laplacian does not lead tangling, provided sufficiently small time step is chosen. The argument follows from the positivity of the determinant of the Jacobian of the deformation, which is preserved due to the eigen-structure of the graph Laplacian.
Remark F.1 (Iterative Application GNN Blocks). The below results extend to our GNN-based mesh deformer, which applies diffusion blocks iteratively. At each iteration, the network updates the node positions while resetting the adjacency weights and initial state Xo- Since each block follows the same form the results can be applied recursively. This ensures that stability and mesh preservation hold across multiple diffusion steps, allowing controlled adaptation of the mesh throughout the G-adaptivity pipeline.
Definition F.2 (Weighted Random Walk Normalized Graph Laplacian). Given a weighted graph g = (V,E,A)) with adjacency matrix A and a learnable weight matrix Ag, where (Ae)is represents the weighted edge between nodes i and j. the weighted degree matrix is defined as Di = Σx(Ag) (Ao)sy- The weighted random walk normalized graph Laplacian is given by:
As = I- Ag.
The operator Ag is symmetric positive semi-dcfinite, satisfying A z O. Its eigenvalues satisfy 0 x S... مح PL... with P&. I 2. The eigenvalues represent the frequecies, and the corresponding eigenvectors are denoted by In-1 ( Je=0
The weights (Aa)i satisfy a, > Of (i.g) e E and Σjen()) au =
In the degree-normalised graph (random walk) Laplacian (Ag A), sums are preserved due to the degree normalization, ensuring E, Aij 1, where À D
In the softmax--wighted case, weights are computed as
exp(f(X,, X,))
(Ao)is EKEN, eX,(
enforcing row stochasticity E, =.
23

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
Definition F.3 (Jacobian of mesh Deformation Map). Given the mesh deformation model Ma (X, + X, the Jacobian J of the transformation is given by:
= VM(X),,
where VM(() is the local derivative of the deformation map.
Definition E.4. (Mesh Tangling). We say that physical mesh is tangled at least one simplex the triangulation has negative determinant in its Jacobian matrix, i.c., det(J+) _ 0, for some where Jj is the Jacobian matrix of the affine transformation mapping the reference element to the physical element in the mesh. Equivalently, the mesh untangled all eigenvalues of the Hessian of the transformation function, its discrete counterpart given by the graph Laplacian, remain positive.
Proof of Mesh Tangling Prevention
We prove that Laplacian GNN-based mesh adaptation scheme prevents tangling, given sufficiently small time step. The argument follows from the positivity of the determinant of the mesh deformation Jacobian, which is preserved due to the eigen-structure of the graph Laplacian.
F.. Continuous-Time Evolution
The evolution of node positions follows the Laplacian- based update:
dX
A I| -AX.
where A is the weighted random-wall graph Laplacian. ordinary differential equation is:
As A is frozen over every diffusion block, the solution of this
((t) = = - ((),,
(25)
implying the determinant of the transformation Jacobian satisfies
[I ( [(0)== = e
((t) det(e-())J(())
(26)
Since tr(A) N₂ 0, we have J(t) > 0 for all 0, ensuring that no elements invert.
F.2 Time Step Constraints for Mesh Preservation
The discrete update for the mesh is x (I = propagating the determinant as J det([ dtA)J*, To prevent inversion, we require det(T dt > 0.
Theorem F.5 (Time Step Condition for Mesh Preservation). Given the discrete update x = (1 - d(A)x the mesh
remains untangled if dt
Proof. The determinant of the deformation Jacobian propagates det Jk+1 det dtt) det J To ensure det J > 0, we require det([- dt > 0. The eigenvalues of I dt A are μu; = 1 dt so the determinant condition reduces to
II dth,)) 0.
Noting that A has positive entries with row sum equal to 1, it follows by Gershgorin's 's theorem that the eigenvalues of A - are contained in the Gershgorin circle |A; - 1| < 1. Seeing as the coefficients of A - are real-valued the eigenvalues of - are either vell-aluedd or come in complex conjugate pairs. If X; is real-valued the contribution to the above determinant is dt > 2dt if > 1/2. If ImA, # then X, also an eigenvalue and the contribution to the determinant is (1 (1 dtA,)) dtReA; 2 + [mA: 2 > Hence we obtain det dt > if dt < 2 LI
24

---

G-Adaptivity: optimised graph-based mesh relocation for finite element methods
F.3. Monitor-Condditoned Time Step
To refine the time step bound, consider the propagation matrix M - dtA with eigenvalues μ4i 1 dt>j. The condition number of M is x(M) 1-dtAmaxx
1-dtAmin Similarly, by Gershgorin's 's theorem, Amax S 2.
Theorem F.6 Chonitor-Conditionel Time Step). Given the discrete update x (I dta)x* where the monitor function redistributes the mesh to improve spectral conditioning, the time step satisfies
(1 {(MM)
dt < min
2 2
Proof. The local mesh determinant propagates _k+1 det(T d/A)J.. Stability requires 1 df.Amax Since
Amax S 2, we obtain < NIM), completing the proof.
E.4. Mesh Quality Measures
Mesh quality measures are often used as indicators of whether a mesh will be effective when used to solve a PDE. In particular, the mesh-quality metrics are directly related to the conditioning FEM stiffness matrices, meaning poor mesh conditioning leads to numerical instabilities the FEM solvers. Our method is designed to minimise the FE solution error directly, but the inclusion of the equidistribution regularisation ensures that G-Adaptivity leads to meshes that maintain good conditioning while reducing the FEM error.In numerical experiments report aspect ratio strong indicator of this mesh conditioning, but in the relevant literature the following metrics are commonly used to often used to assess mesh scale, skewness, and regularity.
Two paradigms exist for evaluating mesh quality:
Known deformation map: Mesh quality is assessed directly using the eigenvalues λo a of the Jacobian.
Local geometric properties: Skewness can be measured the ratio of the circumeircle incircle radius, while regularity is inferred from element area variance.
Mesh quality can be quantified through:
Scale: Element size, measured λσd compared to natural length scale.
Skewness: The anisotropy elements, given by A1/No-
Regularity: Consistency of adjacent elements, e.g., variance in element areas.
Consistency: Stability of element shapes across the domain.
Aspect Ratio for our evaluation we use the aspect ratio of a triangular element, which defined as the ratio of the longest edge Iiax to the shortest altitude hmiu:
1max
AR
(27)
hmin
where Aain is the shortest perpendicular distance from the opposite vertex the longest edge. higher aspect ratio indicates more elongated elements.
25


Reviews:
[{"review_id": "OOBekUE3Lo", "reviewer": "Reviewer_MACS", "summary": "In this work, a GNN-based mesh relocation method is proposed by directly minimizing the finite element solution error. A diffusion-based GNN-deformer is applied which can reduce mesh tangling. Experiments show the proposed method achieves lower solution error, on Poisson's, Burgers', and Navier-Stokes equation problems.\n\n## update after rebuttal\nThe authors have carefully addressed the comments and suggestions from all reviewers and have provided abundant extra experimental results. I will remain my already positive score unchanged.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review of our manuscript and for your valuable comments and questions, particularly for encouraging us to perform **experiments on more complex geometries and larger scale problems,** and to provide **clarification on our theoretical results and the use of ML-based r-adaptivity in industry**. We appreciate your positive feedback, hope the following responses answer your questions and encourage you to consider raising your score for acceptance.\n\n### **Methods and evaluation criteria**\n1. We have expanded our experiments to include **more complex, non-convex domains, and high resolution meshes** (see [Figures C & D](https://imgur.com/a/rOdOAA0)). Additionally, our model **generalises naturally to 3D**, and we tested it on Poisson's equation on the unit cube with Dirichlet BCs and Gaussian solutions. The results presented below and in [Figure A](https://imgur.com/a/rOdOAA0) show that the method outperforms MA significantly (out-of-the-box UM2N does not apply in 3D) and leads to **effective mesh point concentration in regions of interest**.\n|Model|Error Red. (%)|Time (ms)|Aspect|\n|-|-|-|-|\n|MA|12.71 ± 0.00|41049|2.97 ± 0.00|\n|G-Adapt|28.08 ± 0.36|494|6.91 ± 0.20|\n2. We appreciate the reviewer’s note that aspect ratio is not essential for good meshes, which is a key aspect of our framework that aims to break the existing paradigm of classical meshing techniques: **minimal FEM-error meshes do not necessarily require a small aspect ratio**. While our primary loss term is the FEM error $E(\\mathcal{Z},U_\\mathcal{Z})$, **mesh quality also affects FEM-stiffness conditioning**, and high aspect ratios can introduce numerical instabilities in FEM solvers. We report the aspect ratio to show that G-Adaptive meshes achieves error reduction while maintaining reasonable conditioning. A clarification has been added in Appendix F.4.\n\n### **Theoretical claims**\nThe assumption that $d\\tau$ needs to be \"sufficiently small\" to avoid mesh tangling **only requires $d\\tau<0.5$**. This was shown in the original Appendix F.2 and we have now updated the statement of Theorem 4.2 to include this in the main body of the manuscript.\n\n### **Questions**\n1. r-adaptivity is a newer technology than h-adaptivity and as such is not yet widely adapted in industry. However, it has certain significant advantages over h-adaptivity. In particular it works with **a constant data structure, is easy to use on parallel architectures, it gives a more regular mesh** (often with guaranteed mesh regularity), it naturally **inherits Lagrangian and scaling structures in a PDE** (which is very useful for example in ocean modelling and studying PDEs with singularities), and can be **easily linked to existing external software** designed to solve a PDE on an unstructured mesh (for example a discontinuous Galerkin solver). As a result, **r-adaptive methods have recently been very successfully used**, for example, in the operational data assimilation codes of **national weather forecasting offices**, which when coupled to the computational dynamical core, have led to a very significant increases in computational accuracy, particularly for resolving local weather features such as fog and ice (Piccolo \\& Cullen, Q. J. R. Meteorol. Soc., 2012). r-adaptivity has also found natural applications in the **steel industry** where the Lagrangian nature of the approach is very well suited to the evolving fine structures in the forging process (Uribe et al., Finite Elem. Anal. Des., 2024). **Possible disadvantages of r-adaptivity, such as excessive mesh computation cost, and a tendency to mesh tangling, are exactly the issues we address in this paper, proposing a fast and accurate method which avoids tangling.**\n2. The **direct optimization method** is used in Fig. 1 **purely for exposition**, showing that MA-meshes are not necessarily optimal. DirectOpt computes the optimal mesh for a given PDE with known solution but is extremely slow and relies on data which is not available during inference. In contrast, once trained, **our G-Adaptive approach yields fast online mesh movement** without needing reference solution values. However, inspired by your comment we have added the DirectOpt results to Table 1:\n|Model|Error Reduction (%)|Time (ms)|Aspect|\n|-|-|-|-|\n|DirectOpt|27.40 ± 0.00|126,028| 33.99 ± 0.00|\n|MA|12.69 ± 0.00|3,780|2.11 ± 0.00|\n|UM2N| 6.83 ± 1.10|70|1.99 ± 0.03|\n|UM2N-G|16.40 ± 2.65|30|2.61 ± 0.17|\n|**G-Adapt**|21.01 ± 0.33|88|2.92 ± 0.03|\n\n3. We agree that our architecture is fundamentally a Neural ODE on a graph. However, the specific form of this differential equation is crucial to the success of our method: The governing equation, $\\dot{\\mathcal{Z}}(\\tau)=(\\mathbf{A}_{\\theta}(\\mathbf{X}^k)-\\mathbf{I}) \\mathcal{Z}(\\tau)$ resembles a discretized learnable diffusion equation, motivating our use of the term \"diffusion-based GNN-deformer.\". This terminology aligns with prior literature on similar architectures (Chamberlain et al., 2021a;b).", "claims_and_evidence": "The improvement in model structure and loss design makes sense, and is validated by experiments and ablation studies.", "methods_and_evaluation_criteria": "- The results on larger-scale, 3D, and more complicated geometry problems will further validate the performance of the proposed method.\n- In the experiments, aspect ratio is taken as a metric. However, this is not necessarily required for good meshes, especially for this work, the final target (lower FE error) is taken to supervise the training, so we may get anisotropic good meshes. Therefore, I am not sure if it is appropriate to use aspect ratio as an evaluation metric.", "theoretical_claims": "I think Theorem 4.2 is solid. Note that some strong assumptions are required, such as sufficiently small timesteps. Hence in practice mesh tangling can still happen in extreme cases.", "experimental_designs_or_analyses": "The experimental designs and analyses are sound.", "relation_to_broader_scientific_literature": "This work can be applied in downstream industrial physical simulations such as fluid, structural, heat simulations, etc.", "essential_references_not_discussed": "To my knowledge, I don't see any essential references not discussed.", "other_strengths_and_weaknesses": "None.", "other_comments_or_suggestions": "None.", "questions_for_authors": "- In industry, normally people don't care about guaranteeing the topology of the mesh, as long as efficient and accurate solutions can be provided. Can the authors provide some reasons or scenarios where r-adaptation is mandatory or better than h-adaptation?\n- Why is the DirectOpt method shown in Fig. 1 not compared in the experiments?\n- I am not sure if the proposed mesh deformer should be called \"diffusion-based\" or \"neural-ode-based\"?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3E48UVhjM7", "reviewer": "Reviewer_jV3d", "summary": "This paper focuses on using graph neural network to predict deformation of the computational domain in order to reduce error of solution obtained by finite element (FE) method. The philosophy follows [1]. The contribution of this paper is three-fold. 1. The authors proposed a new design on the model architecture based on flow/velocity-type method (eq. (7) and (8)). 2. In order to train the flow-type model, the training utilized a differentiable solver (Firedrake) and directly minimizing FE solution error. 3. A new regularization term is added to the loss function (eq. (9)).\n\nThe authors benchmarked on their own dataset. Particularly, their dataset is featured with convex domain, e.g., square and convex polygons.\n\n[1] Zhang, M., Wang, C., Kramer, S., Wallwork, J. G., Li, S.,\nLiu, J., Chen, X., and Piggott, M. D. Towards universal\nmesh movement networks. 2024.\n\n## update after rebuttal\n\nThe authors supplement experiments on non-convex domains and I changed my score to 3. However, as I commented in my review, why did not the author benchmark the dataset from [1] in the first place? I did not see any reason for not doing so. Then the author agreed the necessity of the comparison and conducted experiment in rebuttal. It seems to me that the design of experiment in the original version is lacking of thoughts to miss such basic and important benchmark. Without seeing a final version of the paper, the credibility of the supplemented experiment is weakened.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review of our manuscript and for your valuable comments and questions, following which we have performed several additional numerical experiments. **We appreciate your suggestion to benchmark more closely against UM2N, and we have now conducted extensive new experiments, including on non-convex domains.** We hope that we were able to address your concerns and that you can update your score accordingly.\n\n\n### **Methods and evaluation criteria**\n\n1. While we already benchmarked on examples from UM2N (Zhang et al., 2024), namely the flow past a cylinder, and on two datasets from M2N (Song et al., 2022), namely the Poisson and Burgers' equation on a square domain, **we have now obtained additional domain data from (Zhang et al., 2024).** This allowed us to conduct extensive additional experiments on **five non-convex domains** (four from the paper (Zhang et al., 2024) and an L-shaped domain). On each domain we solve Poisson's equation for randomly sampled Gaussian solutions with 100 training datapoints and 100 unseen test datapoints. The results (see [Figure D](https://imgur.com/a/rOdOAA0) and table below) confirm that our method performs robustly on non-convex geometries, achieving significantly greater error reduction than baselines and generating regular non-tangled meshes on all tested domains, succeeding even when some other approaches fail. Note that the UM2N results reported below were obtained using the pretrained model from the UM2N repository, and we continue to investigate and refine this baseline.\n\n#### Table 1: Error reduction (%) for various methods and non-convex domains (n.b. negative indicates error increase)\n|Domain|MA|UM2N|UM2N-G|G-Adapt|\n|-|-|-|-|-|\n|Geometry 1|0.23 ± 0.00|-76.85 ± 0.00|1.92 ± 0.02|**7.97 ± 0.04**|\n|Geometry 2|-1.00 ± 0.00|-83.88 ± 0.00|0.69 ± 0.08|**8.88 ± 0.24**|\n|Geometry 3|-| -75.82 ± 0.00| -0.96 ± 0.04|**6.62 ± 0.09**|\n|H-Geometry|-108.31 ± 0.00|-73.59 ± 0.00|-0.92 ± 0.00|**7.51 ± 0.00**|\n|L-Geometry|-89.40 ± 0.00| -138.43 ± 0.00|13.94 ± 1.18|**16.25 ± 0.25**|\n\nMesh deformation times and aspect ratios can be found in [Figure D](https://imgur.com/a/rOdOAA0).\n\n2. **Loss regularization:** We would like to clarify a key misunderstanding: the primary novelty in our loss is **not just the regularisation term in (9) but also the FEM error term** $E(\\mathcal{M}_{\\theta})$. Prior works relied on MSE loss to classically adapted meshes, whereas we directly minimise the FEM error, which is the **key driver of performance improvements**. The equidistribution loss is a secondary regulariser that **further enhances FEM error reduction** by ensuring that the monitor function (not the nodes) is equidistributed. This is in line with the reviewer's own intuition: \"the density of mesh nodes should increase in the area where error is high\", as the nodes will be distributed to areas where the curvature of the solution is high. Our experiments already highlight the advantages of our FEM loss component (UM2N vs UM2N-G). We further conducted an **ablation study** on the Poisson dataset showing that our choice of **weight 1 is optimal for this regularising loss term**.\n\n|Regularization weight|Error Reduction (%)|\n|-|-|\n| 0 (no equi-dist loss) | 22.42 |\n| 0.5 | 22.95 |\n| 1 | 23.99 |\n| 2 | 23.21 |\n| 4 | 22.14 |\n| 8 | 20.96 |\n\n### **Other strengths and weaknesses**\nThrough our extensive additional experiments we demonstrate that **our method is not restricted to convex domains**. Across all cases, **our approach significantly outperforms the two main baselines, UM2N (ML-based) and MA (classical method), in terms of FEM error reduction, while achieving comparable computational efficiency to UM2N.**\n\n### **Other comments or suggestions**\nWe appreciate the point raised concerning the motivation and advantages of our approach.\n\n1. A **central novelty and advantage** of our method is that it optimises the FEM solution error directly, which is in contrast to prior work, including UM2N, which had designed **surrogates** to classical meshing approaches (such as MA). These classical methods rely on heuristics and cannot directly minimise the FEM error.\n2. The flow-based approach, i.e. the diffusion deformer components in our GNN architecture are motivated by relaxation-based mesh movement and provide a **clear, quantifiable advantage over GAT-based deformers** used in UM2N as observed in our experiments (UM2N-G vs G-Adapt).\n3. In contrast to prior work, our approach **does not use MA-solutions as a basis for training** and as such can be seen as more data-efficient. It nevertheless requires the repeated solution of the training problems during the training part (no additional solve is required during inference).\n\n### **Questions for Authors**\n1. Thank you for raising this question. Our method is indeed **fully applicable to non-convex domains** as demonstrated by our additional experiments. We have updated our manuscript to explicitly state this and to showcase the above results.", "claims_and_evidence": "See below.", "methods_and_evaluation_criteria": "A more systematic comparison between your work and UM2N [1] should be made:\n\n1. Why not benchmarking the same dataset of UM2N [1]? I do not see how their dataset could not be your choice and why you have to choose your own dataset instead. Particularly, the dataset of UM2N includes non-convex domains. In contrast, your dataset contains only convex domains except a cylinder case for Navier-Stokes equations.\n\n2. From your experiment, all table 1,2,3 shows that  the regularization term in loss function (9) contributes most to your method, since UM2N-G (UM2N + the loss function (9)) also has great improvements from vanilla UM2N. Essentially, as one of the core novelties of this paper, this regularization term is a density control that encourages equidistribution of mesh nodes. Why this regularization is so effective on your dataset is under-explored. Is it as effective on the dataset of [1]? It is actually quite doubtful because the density of mesh nodes should increase in the area where error is high. Equidistribution heuristically contradicts with this strategy.\n\nTherefore, I strongly encourage the authors benchmark their method on the dataset of [1].\n\n[1] Zhang, M., Wang, C., Kramer, S., Wallwork, J. G., Li, S.,\nLiu, J., Chen, X., and Piggott, M. D. Towards universal\nmesh movement networks. 2024.", "theoretical_claims": "N.A.", "experimental_designs_or_analyses": "See above.", "relation_to_broader_scientific_literature": "N.A.", "essential_references_not_discussed": "N.A.", "other_strengths_and_weaknesses": "**Novelty** The methodology adopted in this work has certain novelty, i.e., flow-based model + training with differentiable solver.\n\nHowever, there is an apparent **weakness** in your experiment since your dataset is biased to convex domain, and therefore it is not so convincing that your method is superior to the baseline UM2N.", "other_comments_or_suggestions": "There could be a potential to fully justify any advantage of your flow-based model over existing machine learning method, e.g., UM2N. For example, is your method more data efficient? Since your method directly minimizes FE solution error with a differentiable solver, what benefit can such ``hybrid'' strategy can provide?", "questions_for_authors": "1. Is your method only applicable to convex domain? I see you have convex setup in your experiment. If yes, it'd be better to declare this requirement of convex domain clearly in your methodology as well.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SB4vWjW2eU", "reviewer": "Reviewer_Kufx", "summary": "This paper presents a novel approach to mesh relocation (r-adaptivity) in finite element methods (FEMs). Traditional r-adaptive methods optimize mesh geometry by solving additional meshing PDEs, which are computationally expensive. Recent machine learning (ML) methods focus on learning surrogates for these classical techniques. In contrast, this work introduces a graph neural network (GNN)-based approach that directly minimizes the FEM solution error rather than relying on surrogate error estimates. The proposed G-Adaptivity framework employs a diffusion-based GNN mesh deformer, which optimizes mesh point locations while enforcing non-tangling constraints. The paper claims that this outperforms both classical PDE-based methods and previous ML approaches in terms of FE solution accuracy while maintaining the computational efficiency of ML-based solutions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review of our manuscript and for your positive comments and questions. Below we provide further clarifications on the **role of Firedrake, hyperparameter choices** and **generalisability of our approach**. We hope the results of the additional experiments and our responses to your questions are sufficient to encourage you to consider raising your score for acceptance.\n\n### **Claims and Evidence**\nWe should clarify that a **fine grid reference solution is only required during training** of the neural network, not during inference. The Firedrake adjoint computations are only invoked during the training phase and **not** during mesh adaptivity. Thus **our method is indeed as efficient as reported** (online mesh movement in a few dozen milliseconds).\n\n### **Theoretical claims**\nThe statement that $d\\tau$ needs to be \"sufficiently small\" means that it **only needs to be smaller than 0.5**. This was shown in the original Appendix F.2 and we have now updated Theorem 4.2 to include this in the main body of the manuscript.\n\n\n### **Supplementary Material**\nWe have added a more detailed discussion of the role of Firedrake in the supplemenatry material of the revised manuscript, and include a shortened version of this here.\n\nTraining the GNN requires computing the derivative of the loss function $E(Z,U_Z)$ with respect to node coordinates $Z$, using adjoint models for efficiency. Automating their derivation is essential for a general $r$-adaptivity methodology that works across different test cases. **Firedrake is ideal for this, as it derives adjoint models and computes these derivatives automatically** (Ham et al., Struct. Multidiscip. Optim., 2019). Obtaining corresponding formulas by hand is difficult: e.g., for $J(Z,U_Z) = ||U_Z||^2_{L_2(\\Omega)}$, which is a simplified version of $E(Z,U_Z)$ in (3), the corresponding derivative takes the form $dJ(Z,U_Z)[T] = \\int_\\Omega (U_Z^2+\\nabla U_Z\\cdot \\nabla p - 4p) \\nabla\\cdot T - \\nabla U_Z (DT+DT^\\top)\\nabla p dx$ with $p$ being the (weak) solution of the adjoint equation $\\Delta p = 2U_Z$. This automatic approach can then be coupled with PyTorch using Firedrake’s ML integration (Bouziani et al., arXiv:2409.06085) to enable GNN training.\n\n\n### **Relation to Broader Scientific Literature**\n**Limitations of ML-based meshing:** The ML-based approach is inherently statistical, meaning that **GNN-based meshing tools are likely to perform worse on out-of-distribution test data.** We observed this in our experiments both with pre-trained UM2N models and our own G-Adaptive approach when applied to PDEs whose solutions featured vastly different scales and features than those on which the models were trained.\n\n### **Questions for Authors**\n1. Our method **extends naturally to 3D** adaptive meshing and we have performed an additional experiment on a cubic mesh with a selection of random Gaussians comparable to the experiment from Table 1 in our manuscript. The results are shown in [Figure A](https://imgur.com/a/rOdOAA0), and the following table (note out-of-the-box UM2N cannot be applied to 3D problems):\n\n|Model|Error Reduction (%)|Time (ms)|Aspect|\n|-|-|-|-|\n|MA|12.71 ± 0.00|41,049|2.97 ± 0.00|\n|G-Adapt|28.08 ± 0.36|494|6.91 ± 0.20|\n\n\n2. We have performed extensive ablation studies and found that our approach is **not very sensitive to the choice of hyperparameters**, as long as they remain in a reasonable range. In particular, we emphasize that **all experiments in the paper used identical hyperparameters without fine-tuning** to specific problems or PDEs.\n\n#### Table 1: The effect of $d\\tau$ and diff. timesteps on the Error reduction (%)\n|$d\\tau$\\No.-timesteps|2|4|8|16|32|64|\n|:-:|-|-|-|-|-|-|\n|0.05|10.41|16.60|18.95|15.85|21.82|22.93|\n|0.1|12.97|14.52|20.68|15.71|20.27|20.85|\n| 0.25 | 19.54 | 19.11 |22.30 |19.94 |23.11 |22.09 |\n| 0.5 | 20.43 | 22.65 |22.14 |22.42 |21.32 |21.92 |\n| 1 | 20.60 | 21.57 |21.16 |22.10 |19.71 |19.40 |\n\n#### Table 2: The effect of $d\\tau$ and diff. timesteps on inference time (ms)\n|$d\\tau$\\No.-timesteps|2|4|8|16|32|64|\n|:-:|-|-|-|-|-|-|\n| 0.05 | 60 | 44 | 46 | 116 | 65 | 247 |\n| 0.1 | 54 | 42 | 79 | 61 | 86 | 208 |\n| 0.25 | 41 | 40 | 48 | 56 | 119 | 108 |\n| 0.5 | 50 | 58 | 49 | 92 | 125 | 158 |\n| 1 | 52 | 59 | 45 | 69 | 100 | 108 |\n\n#### Table 3: The effect of equi-dist loss regularisation on Error reduction (%)\n|Reg. weight|Error Reduction (%)|\n|-|-|\n| 0 (no equi-dist loss) | 22.42 |\n| 0.5 | 22.95 |\n| 1 | 23.99 |\n| 2 | 23.21 |\n| 4 | 22.14 |\n| 8 | 20.96 |\n\n3. This is an outstanding suggestion and in fact part of ongoing work by the authors. It is natural to start by relocation to find an optimal meshpoint distribution followed by h-refinement in regions that require particularly close resolution, see (Dobrev et al., Eng. Comput., 2022) and (Piggot et al., Ocean Model., 2005). The flexibility of the current and ML approaches more generally makes them ideal candidates for such an hr-adaptive approach.", "claims_and_evidence": "The claims made in the submission are generally well-supported by empirical evidence. However, the Firedrake adjoint optimal gradient computation requires clarification. The authors remark that when the exact solution is unavailable, they approximate it using FEM interpolation on a high-resolution mesh. However, running a high-resolution FEM solver is computationally expensive. Is this additional computation included in the reported time comparisons? If not, the time efficiency claims should be adjusted accordingly.", "methods_and_evaluation_criteria": "The proposed method is well-motivated for r-adaptive meshing. It replaces classical PDE-based relocation with a GNN trained via backpropagation through a differentiable FEM solver (Firedrake). The evaluation criteria include: 1. Error reduction: FEM L2-error reduction relative to the baseline mesh. The results show that the proposed method achieves the best error reduction across all baselines. Notably, it succeeds in improving error reduction in the Burgers' Square Rollout, where both UM2N baselines fail. 2. Computational time:  The time required for mesh relocation. The proposed method outperforms the classical Monge-Ampère (MA) approach but is 2-3× slower than UM2N-G. 3. Mesh quality:  Evaluated using the aspect ratio of deformed meshes. The proposed method performs well overall but is outperformed by some baselines in specific cases (e.g., Poisson Convex Polygon problem).", "theoretical_claims": "The theoretical claims appear sound and well-supported. The proof of non-tangling (Theorem 4.2) is particularly valuable, as it ensures that mesh deformation does not lead to degenerate elements. No major concerns were identified.\n\nOne related question: while a sufficiently small pseudo-timestep dτ is theoretically guaranteed to prevent mesh tangling, what is the practical numerical threshold for dτ? Specifically, how is this value determined in experiments, and does it require tuning for different PDEs or mesh resolutions?", "experimental_designs_or_analyses": "The whole evaluation parts consists of three problems: the Poisson's equation, Burgers' equation and Navier-Stokes equation. The comparison with baselines considers error reduction, computational time, and aspect ratio, which are appropriate metrics. However, there are two notable gaps in the experimental design. first it lacks 3D evaluation: all of these problem are solved in the 2D space. Extending the method to 3D surfaces or volumetric meshes would better demonstrate its scalability and generalizability. Moreover, the paper does not analyze the effect of different loss function components (e.g., equidistribution loss) on performance. An ablation study would clarify the contributions of individual components.", "relation_to_broader_scientific_literature": "This paper is related to classic r-adaptivity, ML-based PDE solvers, and graph-based learning approaches.  It highlights the key limitation of prior ML-based r-adaptive meshing (UM2N), which relies on learning surrogates rather than directly optimizing the FEM loss. However, additional discussion on potential limitations of ML-driven meshing (e.g., generalization to higher dimensions) would be beneficial.", "essential_references_not_discussed": "Good to me.", "other_strengths_and_weaknesses": "This method proposes a novel GNN-based approach that directly minimizes the FEM error, avoiding heuristic error estimates. It also provides theoretical guarantees for mesh regularity and non-tangling. The method significantly improves computational efficiency over classical methods, supported by comprehensive experiments on various PDEs and mesh topologies.\n\nHowever, the proposed method is slower than UM2N baselines, which is not discussed in the paper. Besides, there is no explicit ablation study to quantify the impact of different model components, and limited discussion on generalization to 3D meshes or more complex PDEs.", "other_comments_or_suggestions": "I wonder like to see additional visualizations comparing mesh evolution across timesteps and discuss failure cases where G-adaptivity may not perform well.", "questions_for_authors": "1. How does G-Adaptivity generalize to 3D adaptive meshing? Have you tested it on volumetric meshes or 3D curved surface meshes on complex geometry?\n2. How sensitive is the model to changes in hyperparameters (e.g., number of GNN layers, training iterations)?\n3. Can the proposed approach be combined with h-adaptivity for further improvements?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "9RPuvpxqHH", "reviewer": "Reviewer_TRGq", "summary": "This paper proposes to apply graph neural network (GNN) architectures for mesh relocation in finite element methods (FEM). The technical contributions involve a novel training pipeline and improved network structures together with appropriate loss functions. Experiments demonstrate the effectiveness and potential of the proposed learning paradigm.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your positive feedback and recognition of the importance and potential impact of our work. Your comments encouraged us to **test our model at significantly larger scales, including 3D simulations**. We hope the [additional experiments](https://imgur.com/a/rOdOAA0) we performed adequately address your concerns, further strengthen the paper, and provide sufficient grounds to upgrade your recommendation to a clear acceptance. We agree that scalability is crucial for the relevance of any adaptive meshing approach. It turns out that our approach scales very well to larger problems, in three ways.\n\n### **Scalability of the GNN Model**\n**Firstly,** in forward mode, the diffusion deformer is able to scale to very large meshes by design. In particular the inductive learning property of GNNs ensures the ability of GNNs to transfer to unseen graphs in this case meaning we can perform super-resolution to scale to very large meshes. Our experiments demonstrate that our model can efficiently relocate **tens of thousands of nodes in just over two seconds on a standard laptop.**\n\nGNNs have been widely observed to scale well. For instance: \n\n- The \"OGB-LSC: A Large-Scale Challenge for Machine Learning on Graphs\" (Hu, 2021) demonstrates GNNs frequently being scaled to large graphs of millions/billions of nodes and edges.\n- Recent works like GraphCast (Lam, 2023), Aurora (Bodnar, 2024), MeshGraphNets (Pfaff, 2020) and Scalable Universal Physics Transformers (Alkin, 2025) all evidence the ability of transformer and GNN architectures to scale to very large scale simulations.\n\nTo substantiate our claims, we conducted **new experiments on a larger 150x150 mesh (22,500 nodes)** for the Poisson problem with 128 sampled Gaussians (see [Figures C](https://imgur.com/a/rOdOAA0)). Our model **consistently achieved significant mesh adaptation, accuracy improvement, and computational acceleration compared to Monge-Ampere (MA), matching the performance observed on smaller-scale experiments**.\n\n|Scale|Model|Error Reduction (%)|Time (ms)|\n|-|-|-|-|\n|60x60|MA|11.94 ± 1.50|23,084|\n||G-Adapt|**27.47 ± 0.89**|**452**|\n|150x150|MA|17.96 ± 1.27|115,395|\n||G-Adapt|**25.70 ± 1.51**|**2,555**|\n\n\nOne design challenge to note is that naive fully connected transformer encoder would create $22,500^2$ edges, making it computationally prohibitive. We overcame this by using a **sliding window (SWIN) style transformer** to capture the monitor function embedding at the mid-length scales. With this design choice the architecture can be scaled up much further thus providing a method that offers significant potential for real-life FEM applications. Moreover, for even larger-scale problems, **leveraging Firedrake and mesh hierarchy (multi-grid) techniques** could further enhance scalability, which we plan to explore in future work.\n\n### **Extending to 3D Simulations**\n**Secondly,** to further assess **scalability for real-world applications**, we expanded our method to 3D domains, prompted by the suggestions by reviewers Kufx and MACS. In [Figure A](https://imgur.com/a/rOdOAA0) and the below table we demonstrate that the G-Adaptivity framework and diffusion deformer model are easily adapted to the 3D setting performing an experiment on a 10x10x10 unit cube for the 3D Poisson problem. These results confirm that our approach leads to **highly competitive error reduction in the 3D setting while maintaining computational efficiency.**\n\n|Model|Error Reduction (%)|Time (ms)|Aspect|\n|-|-|-|-|\n|MA|12.71 ± 0.00|41,049|2.97 ± 0.00|\n|G-Adapt|28.08 ± 0.36|494|6.91 ± 0.20|\n\n### **Scalability of the FEM solver**\n**Finally,** the FEM solver which is used for the training of our GNN can also be efficiently scaled. **If a scalable solver is known for a particular PDE, then it can be easily adapted into the G-Adaptivity framework as Firedrake supports “sophisticated, programmable solvers through seamless coupling with PETSc”**. A good reference for this is (Kirby & Mitchell, SINUM, 2018). This allows us to leverage domain expertise and achieve scalability when this is needed. In our experiments, the default PETSc options of the Firedrake class `NonlinearVariationalSolver` suffices for competive FEM solution times.\n \nThe same applies to solving the adjoint equations generated with Pyadjoint. Firedrake allows passing the additional key \"adj_args\" to the state solver parameters to specify which solver should be used for the adjoint equation thus enabling **scaling in the same manner as the forward solver**.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "Yes.", "experimental_designs_or_analyses": "Yes.", "relation_to_broader_scientific_literature": "Well related to GNN and FEM literatures.", "essential_references_not_discussed": "No.", "other_strengths_and_weaknesses": "From my point of view, applying GNN for mesh deformation in FEMs is an interesting and highly promising approach. My major concern is the scalability of GNN architectures when the number of nodes increases significantly. This issue may restrict the applicability of the proposed method in other more complicated computation scenarios. The authors should provide more detailed discussions.", "other_comments_or_suggestions": "Overall, the technical scope of this paper is quite different from my expertise, making it difficult for me to provide detailed and in-depth evaluations of the specific technical implementations and experimental setups. However, I believe that the approach of applying GNNs for mesh relocation in this paper is rather valuable.", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: ## Strengths
- Learning-based fast r-adaptivity to optimize meshes for more accurate FE simulations with GNN and diffusion, which is important in practical applications
- Novelty in the technique directly optimizing meshes using GNN and diffusion trained with the FE loss using a differential solver, Firedrake, in contrast to former learning-based approaches, which learned the FE surrogates to optimize meshes
- Solid theoretical claim about discrete time non-tangling, ensuring that the mesh tangling cannot happen if the timestep is sufficiently small (<0.5)
- Versatility to non-convex and 3D meshes

## Weaknesses
- As the authors note, the generalization capability of learning-based methods is inherently bounded by the training data. More insights into the severity of the OOD inferencing in practical applications and how to overcome them is appreciated.
- Videos showing the optimization process and the simulated results of this method and competing methods will help the community to appreciate this method more
Decision: Accept (spotlight poster)