UGM2N: An Unsupervised and Generalizable Mesh Movement Network via M-Uniform Loss

We sincerely thank the reviewer YRmJ for their valuable comments and suggestions, which help significantly improve the quality of this manuscript. Below, we provide a point-by-point response to your feedback, with annotations indicating whether each item pertains to a Question (Q) , Weakness (W) , or Limitation (L) in our original submission.

Response to W1\L1

Writing and Typo Issues.

We sincerely appreciate your feedback. In the expression $σ =∫_K ρ(x) dx$, the function $\rho(x)$ should actually be $m(x)$. We apologize for this notational error in our manuscript. For the clarity of our writing，we recognize that certain sections may have been overly technical or insufficiently structured, which could hinder readability. To address this concern, we will take the following steps in revision:

1. Improved Organization – We will refine the logical flow of ideas, ensuring a clearer progression from motivation to methodology and results.
2. Enhanced Explanations – Key concepts and assumptions will be expanded with additional intuitive examples or illustrative diagrams where appropriate.
3. Language Edits – Complex sentences will be simplified, and terminology will be consistently defined to improve accessibility without sacrificing rigor.

We are committed to making the paper more accessible to a broader audience while maintaining its technical depth. Please let us know if specific sections require further attention, and we would be glad to provide additional clarifications.

Response to W2

Justification of the assumption.

Thank you for raising this insightful question. In our M-Uniform loss, we introduced $Var(L_K) ≈ Var(L'_K)$ as a mild assumption. Since our goal is to handle triangular meshes, we also assume that our meshes possess a certain level of quality, which is achievable with modern mesh generation software—even for highly irregular or complex geometries.

In fact, our requirement that each mesh element appears in a similar number of local patches essentially means that each mesh node should have a similar degree (i.e., the number of connections to neighboring nodes). Generally, boundary nodes are relatively few compared to interior nodes, and most interior nodes have a degree of 6 or 5 (with a small fraction being 4, which already indicates a more distorted element). In a 2D test case with 60K mesh elements generated by mature mesh generation software, we observed that 92.96% of nodes had a degree of 6, while 2.96% had a degree of 5—meaning around 95% of nodes fell into these two categories. Therefore, we consider this a reasonable assumption.

To further demonstrate our model's performance when this assumption does not hold, we have tested it on real-world cases with arbitrary geometry boundaries, specifically in the Tohoku tsunami simulation (ocean modeling), where only about 69% of nodes had degrees of 6 or 5. The results showed that our model could still accurately capture wave propagation, with mesh element sizes dynamically adapting to wave propagation. This also highlights our model's generalization capability for highly irregular or complex geometries. Additionally, we have evaluated our model's performance on highly irregular meshes where mesh elements do not appear in a similar number of local patches. The flow field solution for each sample was defined as $\cos(2 \pi (x - \mu_x))$, where $\mu_x$ is a random value in [0, 1]. The error reduction results achieved by our model are presented below:

The experimental results demonstrate that even when mesh elements do not appear in a similar number of local patches, our model can still achieve accuracy improvement in most cases(achieving success in 8 out of 10 cases).

Response to W3\L2

Real-world scenarios. Thank you for your suggestions. Based on your feedback and Reviewer AXEC’s comments, we have added a real-world case with arbitrary geometry boundaries, sourced from Thetis project in ocean modeling. The experimental results demonstrate that our model dynamically adapts mesh elements to accurately simulate wave propagation while maintaining high-quality, distortion-free meshes without tangling, even under complex boundary conditions.

Due to technical constraints (e.g., one case requiring additional data access and another involving website maintenance), we were unable to include the other two real-world scenarios (Columbia River plume simulation and North Sea tidal simulation) at this stage. However, we will make every effort to incorporate all three cases in the revised manuscript to further validate the robustness of our method.

Validation in 3D large-scale scenarios

Thank you for your feedback. Theoretically, our method imposes no restrictions on mesh scale, as our GNN-based model can accept graphs of arbitrary size as input. However, in practice, this often leads to exceeding GPU memory limits.

On a single RTX TITAN GPU with 24GB memory, our model can handle approximately 40,000 mesh elements, with model inference taking about 345ms and interpolation requiring around 61ms per epoch. In comparison, the MA method takes approximately 25510ms, demonstrating that our approach offers significant efficiency advantages over the MA method. We have also verified the proportion of model inference recommendation, interpolation time, and I/O time under different numbers of mesh elements, and the results are as follows. It can be seen that I/O time is the main performance bottleneck, which inspired our future optimization direction.

On the other hand, our method is highly suitable for parallelization because there are no data dependencies between individual node patches. This allows for easy distribution across N GPUs to enable large-scale mesh processing. Additionally, we can trade time for space by processing patches of mesh points incrementally, thereby increasing the manageable data scale. Adopting this approach (without considering additional overheads like patch grouping and aggregation), our model's processing time for meshes on the order of 1 million elements is roughly 10s, which still outperforms the MA method (which already takes 25s for 40,000 elements).

Response to Q

Geometric rotation invariance property.

We sincerely appreciate the reviewer for raising this intriguing question. First, we attribute these minor discrepancies to inherent asymmetries in the underlying mesh triangulation—rotating the test mesh by 90° fails to achieve perfect overlap with the original mesh. More importantly, this observation prompts a fundamental investigation into whether intelligent mesh adaptation methods possess rotation invariance: does rotating the mesh (while preserving nodal flow variables) affect error reduction performance?

In principle, arbitrary angular rotations should not influence error reduction when flow field variables remain unchanged at mesh nodes. However, the black-box nature of neural networks makes this difficult to guarantee theoretically. To empirically validate this property, we have conducted systematic experiments rotating both the Helmholtz case mesh and its flow field through 360° (in 10° increments), evaluating error reduction across all angles for different ground truth solutions. The comprehensive results are presented below:

The results clearly demonstrate that, compared to the MA method, intelligent mesh adaptation approaches exhibit poorer rotation invariance, as evidenced by their larger performance variance ranges. This insight suggests a valuable direction for future research to enhance geometric invariance in learning-based adaptation methods.

While our method shows inferior rotation invariance compared to UM2N, we emphasize that our node-patch processing framework inherently possesses two critical advantages:

1. Scale invariance: The adaptation results remain unchanged under uniform scaling of the mesh.
2. Translation invariance: The output are consistent under arbitrary translational transformations of the grid.

These properties, which M2N and UM2N fundamentally lack, provide significant practical benefits in real-world applications where mesh scales and positions may vary. The limitation in rotation invariance remains an open challenge that requires further investigation in the future work.

Thank you very much for your constructive and insightful feedback. We have carefully considered each point and provided detailed responses. We welcome any further comments and look forward to the next stage of the review process.

Dear Reviewer YRmJ,

I hope this message finds you well.

We would like to sincerely thank you for your initial review of our submission to NeurIPS 2025 (Paper ID: 12196, Title: "UGM2N: An Unsupervised and Generalizable Mesh Movement Network via M-Uniform Loss"). We truly appreciate the time and effort you have already invested in evaluating our work.

Following the rebuttal phase, several of the reviewers engaged in further discussion and updated their assessments considering our responses. We are grateful for the constructive feedback and the opportunity to clarify key points.

We noticed that you have not yet participated in the post-rebuttal discussion. If you are available and willing, we would be deeply grateful if you could join the conversation and share any additional thoughts or questions you may have in light of the authors’ response. Your perspective is highly valued and would be instrumental in helping the AC and the committee reach a well-informed decision.

Thank you again for your time and consideration. We fully respect your schedule and judgment, and we appreciate any further input you might be able to provide.

Warm regards,

The Authors.

Thank you for your thorough rebuttal. I would like to keep the original score unchanged.

Dear Reviewer YRmJ,

First and foremost, we would like to express our sincere gratitude for your valuable review comments. The issues you raised in your review report are highly insightful and have significantly contributed to improving our work. We have provided detailed responses to your concerns in our rebuttal.

In the rebuttal, we have made additional clarifications regarding the points you highlighted—such as the clarity of our formulations, the justification for the assumptions in the M-Uniform Loss, and experiments under real scenarios. Furthermore, we have included three new experiments to: (1) validate the reasonableness of the assumptions in M-Uniform Loss (Response to W2), (2) address the efficiency of our method (Response to W3\L2), and (3) discuss the geometric invariance of our intelligent mesh adaptation approach (Response to Q). We would greatly appreciate it if you could let us know whether our responses have adequately addressed your concerns.

We deeply value your professional judgment, and with the utmost respect, we are writing this follow-up message in the hope of gaining a deeper understanding of any remaining concerns—particularly aspects that may not have been sufficiently clarified in our rebuttal—and to explore whether there might be further opportunities for improvement to help our paper meet the acceptance criteria.

Thank you once again for your time and thoughtful feedback.

Best regards,

The Authors

Response to Concerns Regarding Novelty and Originality

Thank you for your detailed evaluation and thoughtful comments on our manuscript. We appreciate your concerns regarding the novelty of our contributions and would like to clarify the originality of our approach through a point-by-point response.

First, we acknowledge the relevance of the recent work by Rowbottom et al. and appreciate the opportunity to clarify the citation timeline. Our research began in September 2024, during which the latest available version of their paper was the July 2024 preprint, which we cited. At that time, the preprint did not include elements overlapping with our locally-developed M-Uniform loss function. The ICML 2025 version was announced as accepted in May 2025 (the final published version will be available even later), and we only became aware of its final content after submitting our manuscript. Since then, we have carefully reviewed the latest version of Rowbottom et al.'s work. Rowbottom et al.'s work represents a significant and timely contribution to the field, demonstrating impressive theoretical depth and empirical rigor. Their integration of GNNs with Firedrake to minimize finite element solution errors and the Graph DiffFormer architecture represent a notable contribution, particularly in addressing mesh entanglement. We will ensure proper citation of this work in our revised manuscript.

Second, our two key innovations - the M-Uniform Loss and node-patch based mesh adaptation approach - remain fundamentally distinct from Rowbottom et al.'s work. As you noted, our method focuses entirely on local node optimization, while they continue using M2N/UM2N's whole-mesh input paradigm. Although their the ICML version's Section 5.6 mentions patch-based processing, two crucial differences exist: (1) Their approach resembles mesh partitioning rather than our fine-grained node-level processing, which works better for complex geometries where partition may vary significantly; (2) Their patch-based method was actually introduced during the review process (see their rebuttal on OpenReview, which is visible to everyone only after June 18, 2025 ) and was absent from their February 2025 version. This timeline strongly confirms the originality of our approach.

Finally, we would like to highlight additional contributions of our work, including: (1) a fully label-free training approach that relies solely on monitor function values, in contrast to methods requiring FEM gradient computations coupled with Firedrake; (2) a mesh-smoothing-inspired iterative optimization scheme that demonstrates consistent improvement across a variety of problem settings; and (3) zero-shot generalization capability, with our model trained on a limited dataset achieving effective performance across multiple PDE cases, compared to approaches that often require case-specific training.

We hope these explanations will help reassess our work's novelty and originality. Below we address your other constructive comments point-by-point (marked as Q for questions).

Response to Q1

We appreciate the clarification request. As described in Section 3.2 (lines 200–205), our training data consists of pre-generated local node patches extracted from high-quality reference meshes (see Figure 10). During training, these patches are treated as static samples—there is no dynamic mesh update or policy rollout. Therefore, the training distribution remains stable, and there is no distribution shift between epochs.

Response to Q2

In fact, in Figure 4 we have presented the time efficiency of different mesh adaptation methods across various mesh sizes, where our method's timing already includes interpolation operations. The results in Figure 4 demonstrate that our interpolation process does not introduce significant computational overhead, enabling our method to maintain highly efficient single-iteration performance while also achieving substantial reductions in multi-iteration adaptation time.

To provide further evidence, we have conducted timing experiments across different mesh sizes. The results, presented below, show that the main bottleneck is I/O between CPU and GPU, while both inference and interpolation times remain low. These findings confirm the scalability of our method, and we plan to further optimize I/O bottlenecks in future work.

Response to Q3

Thank you for your valuable feedback. Although we have not theoretically analyzed how to guarantee error reduction without mesh tangling, we have demonstrated the excellent performance of our model in error reduction and avoiding mesh tangling through extensive and diverse experiments. Specifically, we have verified the model's performance under varying flow field solutions, mesh resolutions, and mesh shapes, all of which confirm that our method achieves improved solution accuracy without any observed mesh tangling. Overall, our proposed method exhibits stronger robustness and generalization compared to existing traditional methods and ML-based mesh adaptation approaches.

To further illustrate the generalization capability of our model across different scenarios, we have added the following new experiments:

1. Random Geometry Test: On 1,000 randomly generated meshes (with 3–6 edge polygons and complex flow fields), our method improved solution accuracy in 807 cases — the highest among all baselines.
2. Real-world Case: On a challenging tsunami simulation case with complex coastal boundaries, our model successfully adapted the mesh without producing any invalid elements.
3. Non-convex Test: On 100 non-convex domains prone to mesh entanglement, our method achieved the best error reduction performance and produced no tangled meshes.

These results indicate that our model offers not only strong generalization but also practical robustness in avoiding mesh tangling under a wide range of conditions.

We are truly grateful for your detailed and thoughtful comments. They have given us a clearer perspective on both the strengths and limitations of our work. We look forward to your follow-up and are happy to provide further clarification if needed.

The authors gratefully acknowledge the reviewers' valuable feedback, which helped us refine and strengthen the paper. Below we provide our responses to the reviewers' comments, where questions are labeled as Q and weaknesses are denoted as W.

Response to W1\Q1

We sincerely appreciate the reviewer’s valuable suggestions regarding convergence and approximation guarantees. We fully acknowledge the concerns and recognize that Appendix A.2 currently provides only a preliminary theoretical analysis. As noted in our paper, our method draws inspiration from conventional mesh smoothing techniques, which also optimize mesh quality by minimizing a scalar function defined over local elements. However, our approach differs fundamentally in that we optimize the monitor function integral approximation $\mathcal{L}\_{K} =m_{K}|K|$ (Eq. 6), whereas traditional smoothing focuses on geometric quality metrics (e.g., aspect ratio, minimum angle).

While this similarity might suggest that convergence or approximation guarantees could be readily derived, such analysis proves challenging in practice. Under a serial implementation—where nodes are updated sequentially—each iteration provably improves $\mathcal{L}\_{\mathrm{var}} ( \mathcal{P}\_{i} )$ within a local neighborhood, theoretically ensuring convergence. However, the final approximation accuracy is inherently limited by the model’s training error, as the learned mapping cannot guarantee exact minimization (i.e., $\mathcal{L}\_{\mathrm{var}} ( \mathcal{P}\_{i} )=0$) for arbitrary node patches. Moreover, this serial approach is computationally inefficient due to its non-vectorizable nature.

In contrast, our practical implementation employs Jacobi-style parallel updates, which significantly improve efficiency but introduce inherent analytical difficulties. Although the model effectively optimizes local node distributions, the convergence behavior and approximation error remain challenging to characterize theoretically. Unlike Laplacian smoothing, which permits linear analysis, our transformation $T$ operates as a black-box nonlinear function, making it challenging to verify whether it possesses certain properties (e.g., whether it constitutes a contraction mapping).

Further complicating the analysis are the practical challenges of parallel updates, including data races (concurrent modifications of shared regions), update conflicts (local improvements causing global degradation), and information delays (reliance on stale neighbor positions). Additionally, while local node neighborhoods are often geometrically convex, the Hessian of our objective function may not be, precluding standard convex optimization arguments.

These challenges reflect broader theoretical difficulties inherent to adaptive mesh optimization, and are not unique to our model. Nevertheless, as demonstrated empirically in Section 5, our method achieves robust convergence and superior adaptation quality in practice. We have expanded Appendix A.2 to further clarify these theoretical considerations in response to the reviewer’s insightful feedback.

Despite the aforementioned theoretical challenges, we have empirically evaluated the convergence behavior of our model by monitoring the evolution of the nodal displacement norm $\delta_k=||\mathbf{X}^{k+1}-\mathbf{X}^{k}||$ across iterations, where $\mathbf{X}^{k}$ denotes the mesh node positions at optimization step $k$. As demonstrated in the results for the Helmholtz equations (shown below), the consistent decay of $||\mathbf{X}^{k+1}-\mathbf{X}^{k}||$ toward zero with increasing $k$ provides empirical evidence of algorithmic convergence. However, when $k$ increases (greater than 10), we also observed the phenomenon that $\delta_k$ fails to converge or even increases in test cases. In such scenarios, our proposed dynamic termination method can effectively address this issue, thereby enhancing the model's convergence.

Response to W2

We sincerely appreciate your valuable feedback. When dealing with complex boundaries, ensuring the geometric consistency of mesh boundaries is a challenge. The fundamental difficulty stems from the inevitable prediction errors in ML models, where even minor deviations of predicted boundary nodes can significantly compromise accuracy by displacing boundary layers - the critical thin regions near surfaces where viscous, thermal or electromagnetic effects dominate. This justifies our design choice to avoid direct boundary node prediction.

Our new experimental validation reveals that while simple clip operations (projecting nodes onto boundaries) show minimal impact (list below) for internally-dominated flows like Helmholtz and Burgers cases in our paper, they provide substantial improvements for cases involving rapidly varying boundary conditions (detailed in Response Q3 to Reviewer mhZx). These findings indicate that, while our current treatment provides reasonable robustness, there remains ample opportunity to further improve boundary prediction through learned or hybrid approaches.

Response to W3\Q2

Thank you for your feedback. Following your suggestions, we have added two new experiments to evaluate the performance of our method on both highly irregular and anisotropic meshes. The detailed experimental setup and results are presented below.

Highly irregular mesh

We uniformly sampled mesh points within a unit rectangular domain and generated meshes via Delaunay triangulation, followed by mesh smoothing to optimize mesh quality. The meshes generated this way contain numerous mesh points with different degrees, making them fail to satisfy mesh regularity.

We generated 10 samples in total, with the flow field solution on each sample given by $\cos(2 \pi (x - \mu_x))$, where $\mu_x$ is a random value in [0, 1]. The error reduction results of our model are shown below:

The results demonstrate that, although our model was trained on regular meshes, it can generalize well to highly irregular meshes, achieving effective error reduction in 8 out of 10 cases.

Anisotropic meshes

We have added a new test case for anisotropic meshes, where the flow field solution is given by $u(x,y)=(1-e^{\frac{-(x-0.5)^2}{0.01}})((x-0.5)^2-1), [x,y] \in [0,1] \times [0,1]$. First, we generate preliminary stretched mesh elements through anisotropic adaptation to create anisotropic meshes. These generated anisotropic meshes are then fed into our model to further optimize the node distribution to match the characteristics of the solution. In this solution, the Hessian values reach their maxima at approximately x = 0.3, 0.5, and 0.7. Our UGM2N model can capture these three locations more accurately compared to the initial anisotropic meshes. After solving with the optimized anisotropic meshes, the error reduction (absolute value) for each model is shown in the table below, demonstrating that our model has the potential to be applied to anisotropic meshes.

Response to W4

We sincerely appreciate your valuable feedback. In fact, our ablation studies in Section 4.4 have already evaluated the impact of different loss functions on model performance under identical architectures. The results demonstrate that our proposed loss function outperforms those used in previous methods when maintaining the same network structure. We would be happy to conduct further comparisons of loss functions and architectures (including M2N and UM2N) to validate our approach, if of interest.

Furthermore, compared to prior approaches, our model architecture is significantly simpler. This design choice enables highly efficient inference per adaptation epoch (as shown in Figure 4), while maintaining substantially better computational efficiency than MA methods even for multiple adaptation rounds.

We greatly appreciate the time and effort the reviewers have devoted to evaluating our work. Your constructive feedback has deepened our understanding of the key issues and has been extremely helpful. We remain open to further suggestions and look forward to your continued guidance.

Thank you for your constructive comments on our paper. We sincerely appreciate the effort you have dedicated to reviewing our work. Below, we address each of your questions and concerns in detail, with clarifications, new experiments, and implementation insights, where Q represents the question and W represents the weakness. We hope that our clarifications and additional experiments may encourage you to reconsider your score.

Response to Q1

We apologize for the unclear description of the reconstruction process. According to Fig. 1, our method appears to reposition all mesh nodes (including boundary nodes) of each patch back onto the mesh; however, as stated in Line 165 and Algorithm 1, we actually only reposition the central nodes, thereby avoiding issues arising from patch boundary conflicts.

Response to Q2

We greatly appreciate your constructive feedback. First, we would like to address the issue of mesh tangling. Our model utilizes an iterative mesh adaptation process; however, we observed that excessive iteration counts lead to progressive instability and the generation of tangled elements. To resolve this issue, we have implemented a monitoring mechanism for $\mathcal{L}_{\mathrm{var}}\left(\mathcal{M}^{\prime}\right)$ on the mesh and adopted dynamic termination to prematurely halt the adaptation process.

Second, we have conducted the tangle test you recommended in Section 5.3 of UM2N for further validation. We generated 100 meshes, each containing 8 nodes, with four located at the corners of a unit square and the remaining four randomly sampled within the unit square. The flow field was constructed using a mixture of Gaussian distributions. The results are presented below, demonstrating that both our model and UM2N consistently produced valid meshes. In contrast, the MA method struggles with non-convex cases and only generates valid meshes in limited scenarios.

Response to Q3

This is an extremely constructive suggestion. We have added the Tohoku tsunami simulation case adopted in UM2N. The results show that:

1. Our model accurately captures wave propagation, with mesh element sizes dynamically adapting to the wave movement;
2. In this complex boundary simulation, our model exhibits no mesh tangling, ensuring mesh validity throughout;
3. The adapted meshes generated show no significant element distortion, further guaranteeing mesh quality.

We also intended to include Columbia River and North Sea cases, but their bathymetry datasets were temporarily inaccessible. In the revised version, we will make every effort to include these cases and present their results in the supplementary materials.

Response to Q4

Thank you for your comments and feedback.

First, we analyzed the reason for the high variance in the results on the Poisson equation. We have found that this discrepancy primarily stems from the fact that we did not update the mesh boundaries, while the distribution of those solutions leads to large Hessian values at the mesh boundaries. To illustrate this point, we allowed mesh nodes to move tangentially along the boundaries and retested the model's performance. The new experimental results show that when further permitting bound node movement, the model's performance improved significantly.

Secondly, we would like to point out that the r-adaptive mesh method struggles to maintain a consistent error reduction ratio across different problems due to significant variations in physical properties, error distributions, and length scales among different flow fields. Its effectiveness heavily depends on the local characteristics of the flow field and the accuracy of the error indicator. Moreover, constrained by the requirement to preserve mesh topology, it cannot adapt to more complex error structures by adjusting the number of mesh elements. While the r-adaptation is effective, it faces challenges in ensuring consistent error reduction across varying flow fields.

Beyond highlighting this issue, we believe a more critical question is whether the method can reliably reduce errors (e.g., some approaches may degrade accuracy in certain scenarios), as this reflects the algorithm's robustness. To evaluate this, we tested the adaptation performance on 1,000 randomized cases (with random boundaries and flow fields). The results show that our method achieved effective accuracy improvement in 807 cases, demonstrating its robustness across diverse scenarios.

Response to Q5

The adapted mesh is computed from a uniform mesh at each step (the same as M2N and UM2N). At every time step, we project or interpolate the flow field from the previous time step onto the uniform mesh, then compute the current time step's mesh and use it as the input for the next step.

An autoregressive approach is promising, as it could avoid an interpolation step. However, since the input mesh would become more diverse, the model would need to handle a wider variety of mesh distributions—making it challenging to construct a dataset that covers sufficient cases. This issue is certainly worth further investigation in future research.

Response to Q6

Thank you for your suggestion. We have evaluated the performance of the benchmark methods on these two cases, with the following results:

1. MA: For both cases, the MA method exhibited divergence behavior during the MA equation solving and failed to generate valid adapted meshes.
2. M2N: In the moving cylinder case, the M2N did not demonstrate consistent behavior. While it could generate well-adapted meshes at certain time steps, it produced highly distorted mesh elements at others. For the supersonic case, the M2N could capture the shock position, but with lower mesh node density at the shock compared to UGM2N.
3. UM2N: In the supersonic case, although the UM2N could identify the shock position, it introduced element distortion in irrelevant regions. For the moving cylinder case, while it could fit the flow field distribution, excessive adjustments resulted in distorted mesh elements.

These comparisons further demonstrate the robustness of UGM2N under diverse flow conditions and complex boundaries. Additional results will be included in the revised supplementary materials.

Response to Q7

Efficiency Analysis

We appreciate your reminder regarding the insufficient discussion on time efficiency. In Fig. 4, we present the mesh adaptation time (multiple epochs) of UGM2N under various mesh sizes, where the dashed line indicates the time for a single epoch. It can be observed that the single-epoch adaptation efficiency of our model is higher than that of the M2N and UM2N, allowing us to maintain good efficiency even under multi-epoch adaptation.

To further analyze the time allocation of different operations within a single round, we have constructed meshes with varying numbers of mesh elements to test model inference time (as "infer. time"), mesh node interpolation time (as "inter. time"), and I/O time, and the results are shown below:

It can be observed that, in a single adaptation epoch, the majority of the time is consumed by I/O operations (data transfer between CPU and GPU), while model inference time and interpolation time remain relatively low. The I/O time consists of two operations:

1. Transferring mesh data (typically stored on the CPU) to the GPU before adaptive adaptation begins.
2. Transferring the newly adapted mesh nodes from the GPU back to the CPU for interpolation.

If a GPU-based Delaunay triangulation interpolation function were available, the second data transfer could be avoided, significantly improving our model's inference efficiency. However, we have not found a mature GPU implementation of such a function, and developing one falls outside the scope of this work. In the future, we will focus on developing this functionality to further enhance the efficiency of our method.

Convergence

This is an insightful and important observation. We suggest that this phenomenon could potentially relate to data shifts occurring during adaptation. Our model was trained on uniform meshes, and the distribution of adapted mesh elements progressively deviates with increasing iterations. This distributional shift introduces untrained mesh node distributions, causing randomized model outputs that may subsequently lead to mesh distortion and eventual tangling.

The most straightforward solution is our dynamic termination strategy. This strategy automatically halts the adaptation process when the $\mathcal{L}_{\mathrm{var}}\left(\mathcal{M}^{\prime}\right)$ ceases to decrease, thereby preventing further quality degradation. This solution effectively addresses the convergence challenges while maintaining computational efficiency.

Thank you again for your careful review and valuable feedback. We hope our responses have addressed your concerns, and we welcome any additional questions or suggestions you may have.

We sincerely appreciate the reviewers' insightful comments and constructive suggestions, which have greatly helped us improve the quality of our paper. Below, we will address each of the weaknesses (W) and questions (Q) raised in your review comments point by point.

Response to W1

Terminology

We fully acknowledge your concern regarding insufficient explanations of specialized terms. Since our work involves AI for physics, the field encompasses a wide range of specialized terminology due to its interdisciplinary nature, combining concepts from physics, applied mathematics, and computer science. We recognize that additional context would benefit readers with diverse backgrounds. In the revised version, we commit to expand explanations in the main text or appendices and include targeted references for further reading.

Suitability

Machine learning for physics-based simulation is a thriving subfield within AI for Science, as evidenced by numerous seminal works published at NeurIPS (e.g., SchNet, Hamiltonian Neural Networks, and Neural ODEs). The prior AI-based r-adaptation works, M2N and UM2N, also appeared in NeurIPS. We believe this work aligns well with NeurIPS’s scope.

Response to W2

Robustness

We acknowledge that our model does not achieve optimal performance in every single case. However, as evaluated from multiple perspectives in the paper, our approach consistently outperforms existing methods across various scenarios, including:

1. Flow Field Variations:
   • Achieves 9 best results out of 14 cases (64.3%).
   • Outperforms existing ML-based methods (M2N and UM2N) in 12 out of 14 cases (85.7%).
2. Mesh Scale Adaptability (Fig. 3):
   • Maintains error reduction capability across different flow fields.
   • Delivers 14 best results out of 20 test cases (70%).
   • Surpasses ML-based baselines in 15 cases (75%).
3. Variable Geometry Tests: Exceeds both traditional MA methods and ML-based approaches.

These results confirm that our method offers the most consistent performance overall, despite minor limitations in isolated cases. We also want to emphasize that achieving perfect performance in every scenario is inherently challenging—even traditional methods fail in certain cases (see Tab. 1).

Additional Robustness Validation To further validate robustness, we have conducted a new large-scale test on 1,000 mesh samples, where each sample contains random polygonal shapes (3–6 edges) and Gaussian mixture-based flow fields.

We solved the Helmholtz equation using these flow fields as solutions, with the comparative results shown in the table below. The results show that both MA and UM2N suffer from performance variance across different flow fields and geometries. In contrast, our method demonstrates significantly greater robustness in mesh adaptation.

Notably, our approach achieves the highest Positive ER Ratio (807 successful cases out of 1000), substantially outperforming both he traditional method and UM2N. The M2N results have been excluded from the table due to its inherent limitation in generalizing across arbitrary geometries and flow fields without case-specific training.

Response to W3\Q1

Thank you for pointing out the limitations in our experimental design. We have conducted an ablation study to systematically evaluate the influence of parameter $\lambda$ on model performance. We define $\lambda$ = 10**scale, with scale values spanning the range from -1 to 3 in increments of 0.5. To ensure reliability, each $\lambda$ value was tested with 5 repeated experiments. The converged test losses are shown in the table below, which indicate a linear relationship between $\lambda$ and the final converged loss, which aligns with the Eq. (8) of our loss function. This suggests that the value of $\lambda$ has minimal impact on the model's convergence.

The results on the the Poisson and Helmholtz equations are summarized in the below table. Overall, the ER remain relatively stable across different $\lambda$ values, suggesting that our loss function is not highly sensitive to $\lambda$. However, certain cases (e.g., Helmholtz1, Poisson4 and Poisson7) exhibit greater sensitivity to $\lambda$, which we attribute to case-specific effects rather than model limitations, as they prove challenging to optimize across all models.

Response to Q2

This is indeed an important point that warrants clarification. In solving PDEs, we typically know the governing PDE and boundary/initial conditions, while the solution u remains unknown. As you correctly noted, for simple problems we may obtain explicit expressions for u, but for complex flow problems such explicit solutions are generally unavailable.

However, our method does not require explicit expressions for u to compute the Hessian matrix, since we can numerically recover solution derivatives through established numerical differentiation techniques, which is detailed in book Adaptive Moving Mesh Methods P266. This numerical method only requires flow field variables at mesh nodes and computes the Hessian (and gradient) values via the least squares fitting method or the Galerkin formulation method, enabling our approach to generalize to arbitrary flow fields.

Response to Q3

This is an excellent question that provides valuable insights for further model optimization. Our analysis reveals that M2N's superior performance on the first Burgers equation stems from its tendency to generate more deformed meshes (compared to our smoother UGM2N outputs), which better capture abrupt flow variations in this specific case. However, as Fig. 2 demonstrates, this characteristic often leads to mesh distortion in most other cases, ultimately compromising accuracy. This suggests M2N's advantage is case-specific rather than methodological.

Regarding the suboptimal results on Poisson4-7, we identified that our boundary treatment - specifically, not updating boundary nodes - becomes problematic when large Hessian values concentrate near boundaries while interior values remain small. New experimental validation (see table below) shows significant improvement when allowing boundary-constrained node movement. Notably, even without this modification, our model outperforms both M2N and UM2N (which natively permit boundary movement), confirming our approach's fundamental advantages.

Further tests reveal boundary adaptation's differential impact:

• Helmholtz cases: Minimal improvement as flow variations predominantly occur in interior regions
• Burgers cases: Substantial ER boost to 35.27% in Case1 but limited gains in Case2

These findings highlight two critical factors: (1) solution characteristics and (2) boundary treatment quality, particularly for curved boundaries where simple assignment proves inadequate and projection becomes necessary. Developing enhanced boundary handling methods for complex geometries emerges as a key direction for future research.

Response to Q4

We appreciate your thoughtful comments. First, while our model was trained on only 4 flow fields, the use of node patches as fundamental data units resulted in a substantial dataset of 10,440 samples, which proved sufficient for effective model training. The model's demonstrated generalization capability across multiple test cases, enabled by its few-shot learning ability, eliminates the need to prepare large-scale training datasets.

To further investigate the impact of training data volume on model performance, we have conducted additional experiments with expanded flow field configurations (4, 8, and 16 fields), corresponding to training nodes of [10440, 20880, 41760]. The experimental results demonstrate consistent improvement in our model's generalization capability with increasing training data across Poisson, Helmholtz and Burgers equations, with particularly significant gains observed for Poisson and Helmholtz cases. This inspired us to expand richer datasets in subsequent work to further enhance the model's generalization ability.

We sincerely thank the reviewers for their thoughtful comments and insightful questions. Your feedback has been invaluable in helping us reflect on and refine various aspects of our work. We look forward to any further discussion or clarification you may require.

We are very pleased that we have addressed your concerns and achieved an improved score. Your comments have greatly enhanced the quality of our paper. We will include the additional content in the revised version of the manuscript. Finally, we would like to thank you once again for your efforts during the review process. (＾▽＾)