G-Adaptivity: optimised graph-based mesh relocation for finite element methods

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 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.

Scalability of the GNN Model

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.

GNNs have been widely observed to scale well. For instance:

• 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.
• 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.

To 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). 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.

One 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.

Extending to 3D Simulations

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 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.

Scalability of the FEM solver

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.

The 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.

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.

Claims and Evidence

We 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).

Theoretical claims

The 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.

Supplementary Material

We 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.

Training 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.

Relation to Broader Scientific Literature

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.

Questions for Authors

1. 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, and the following table (note out-of-the-box UM2N cannot be applied to 3D problems):

2. 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.

Table 1: The effect of $d\tau$ and diff. timesteps on the Error reduction (%)

Table 2: The effect of $d\tau$ and diff. timesteps on inference time (ms)

Table 3: The effect of equi-dist loss regularisation on Error reduction (%)

3. 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.

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.

Methods and evaluation criteria

1. 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 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.

Table 1: Error reduction (%) for various methods and non-convex domains (n.b. negative indicates error increase)

Mesh deformation times and aspect ratios can be found in Figure D.

2. 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.

Other strengths and weaknesses

Through 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.

Other comments or suggestions

We appreciate the point raised concerning the motivation and advantages of our approach.

1. 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.
2. 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).
3. 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).

Questions for Authors

1. 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.

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.

Methods and evaluation criteria

1. We have expanded our experiments to include more complex, non-convex domains, and high resolution meshes (see Figures C & D). 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 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. |Model|Error Red. (%)|Time (ms)|Aspect| |-|-|-|-| |MA|12.71 ± 0.00|41049|2.97 ± 0.00| |G-Adapt|28.08 ± 0.36|494|6.91 ± 0.20|
2. 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.

Theoretical claims

The 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.

Questions

1. 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.

2. 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: |Model|Error Reduction (%)|Time (ms)|Aspect| |-|-|-|-| |DirectOpt|27.40 ± 0.00|126,028| 33.99 ± 0.00| |MA|12.69 ± 0.00|3,780|2.11 ± 0.00| |UM2N| 6.83 ± 1.10|70|1.99 ± 0.03| |UM2N-G|16.40 ± 2.65|30|2.61 ± 0.17| |G-Adapt|21.01 ± 0.33|88|2.92 ± 0.03|

3. 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).