Graph Neural Networks for Mesh Generation and Adaptation in Structural and Fluid Mechanics

1. Introduction

The generation of meshes constitutes a crucial stage in numerically simulating diverse problems within the domain of computational science [1]. Unstructured meshes find prevalent application, notably in fields like computational fluid dynamics (CFD) [2,3] and computational solid mechanics (CSM) [4,5], wherein the finite element (FE) method is employed to approximate solutions for general partial differential equations (PDE), particularly in domains featuring challenging geometries [6,7]. The accuracy of the FE solution is intricately linked to the quality of the associated mesh.

Mesh generation depends on a number of parameters, such as the local metrics at given points, but also, in the general case, the maximum and minimum sizes of cell elements, the Hausdorff coefficient, the acceptable gradation of the metric, and so on [8]. A set of parameters whose choice is far from trivial and whose poor choice can lead to meshes of poor quality, a generation process that is too long (because it is over-constrained), or even a generation failure.

A mesh can be adapted and refined to suit the physical problem being solved. The adaptation of unstructured meshes has proven very effective for stationary calculations in Fluid Mechanics and Solid Mechanics, with the aim of capturing the behaviour of physical phenomena while obtaining the desired accuracy for the numerical solution. In addition, this method considerably reduces the computational cost of the numerical solution by reducing the number of degrees of freedom. The a priori and a posteriori error estimators generally give an upper bound in the approximation error of the solved PDE. These estimates depend on the approximation space, the mesh size, the equation, the approximate solution in the case of a posteriori estimators and the exact solution in the case of a priori estimators.

Theoretical frameworks for anisotropic error estimation have been extensively developed, resulting in a degree of standardisation in the adaptation process. Various works [9,10,11], have contributed to the estimation of approximation and interpolation errors. Recent analyses of interpolation error [12,13,14,15,16] have refocused attention on metric-based mesh adaptation, particularly utilising a metric derived from a recovered Hessian. Notably, Hessian-based metric mesh adaptation offers several advantages, including a general computational framework, a relatively straightforward implementation, and, most importantly, robustness.

The use of Machine Learning (ML) or Deep Learning (DL) in the field of numerical simulation has until now focused mainly on the prediction of stationary or non-stationary fields as solutions to a given differential equation. This work has been motivated by the desire to reduce the calculation time of an approximate solution by relying on a neural network (NN) for the prediction task, rather than on a FE solver [17,18]. The field has benefited from the emergence of neural networks, and more specifically convolution neural networks (CNNs), which are generally used to extract local and global information from an image, represented as an array of pixels [19]. To reproduce this framework, the simulation domains are represented by Cartesian grids, onto which the geometry and certain initial conditions are projected [20,21,22,23,24,25]. This operation results in the loss of the information carried by the initial unstructured mesh, and produces regular meshes that can contain many more points, resulting in higher computation costs.

To address this constraint, a model can be formulated to perform convolutions on graph structures. The extension of CNNs to graphs holds significant importance, especially given that numerous solvers rely on FE methods utilising unstructured meshes for discretisation in both CFD and CSM.

Message passing neural networks (MPNNs) is a subclass of GNNs [26] that have gained prominence in the domain of mesh-based simulations, offering a framework for handling complex spatial data. The use of MPNNs in the context of mesh-based simulations represents a departure from traditional convolutional approaches, providing a more flexible methodology.

One notable instance of MPNN application is found in the work of Gilmer et al. (2017) [27], where a message-passing neural network is proposed for quantum chemistry simulations. This approach divides the convolution process into two stages: the aggregation of information from neighbouring nodes and edges into a hidden node state, followed by the use of this hidden state to update node features.

Battaglia et al. (2018) [28] introduce a slightly different approach known as graph network (GN). Messages are first transmitted from nodes to edges through an edge convolution kernel, leading to the update of edge features. Subsequently, updated edge features are aggregated to nodes using permutation-invariant operations, forming the edge messages. Finally, a node convolution kernel processes the original node features along with the edge messages to produce updated node features. This method proves to be particularly adept at capturing intricate relationships within graph-structured data.

Pfaff et al. (2021) [29] extended the application of graph networks to mesh-based simulations, focusing on scenarios such as incompressible flow around cylinders and compressible flow around airfoils. Their proposed neural network functions as an accurate incremental simulator with the added capability of adapting to the mesh structure. An intriguing feature of this approach is its ability to generalise well to mesh shapes and sizes not encountered during the training phase, emphasising its robustness in handling diverse simulation conditions.

Yet, this recent literature has not been applied to predict the complete process of AMR, from the initial geometry to the adapted mesh, in three dimensions. To this end, we propose Adaptnet, which is a framework for learning mesh parameters and tetrahedral mesh-based simulations using GNNs (Figure 1). This framework consists of two networks trained to generate unstructured meshes adapted to the physics under study:

• Meshnet is trained to predict mesh parameters of a given geometry (CAD file) later used to generate an initial mesh.
• Graphnet is trained to predict a metric field, either directly or indirectly by predicting the velocity field, later used to generate an anisotropic adapted mesh.

Both models are message passing graph neural networks that rely on the architecture of MeshGraphNets [29,30] that was further adapted to steady-state predictions. We will explain in detail material, methods and results for the case of the steady-state Stokes problem. The linear elasticity problem will be presented at the end to demonstrate the model’s ability to switch from one physics to another. We made use of open-source libraries to direct our work towards open-source solutions and to promote reproducibility.

The remainder of this paper is structured as follows:

• Section 2: Materials and Methods - This section presents the problem statement, describes the dataset used, elaborates on the model architecture, and details the training configuration.
• Section 3: Results - Here, we provide the results of our experiments, focusing on the performance of the proposed models Meshnet, Graphnet, and Adaptnet.
• Section 4: Discussion - This section discusses the implications of our findings, in particular the generalization of our approach to unseen data and configurations.
• Section 5: Conclusion - Finally, we summarize the key contributions of our work and suggest directions for future research.

2. Materials and Methods

2.1. Problem Statement

Both the Stokes flow and linear elasticity problems will be examined using the same geometric configuration. Specifically, each problem will be analyzed within a domain featuring a rectangular box with circular obstacles at its center, allowing for a consistent comparison of the effects of fluid flow and structural deformation on a shared geometry.

2.1.1. Stokes Problem

We consider the steady-state Stokes problem of a fluid flow around a cylinder. The domain $\Omega$ is a box of size $(L\times h\times 1)$ with circular obstacles of radius r in the centre. The boundary $\partial \Omega$ is divided into four parts: $\partial {\Omega }_{in}$, $\partial {\Omega }_{out}$, $\partial {\Omega }_{wall}$ and $\partial {\Omega }_{obs}$, corresponding to the boundaries of inlet, outlet, wall, and obstacles, respectively. The problem is defined as follows:

$$\left\{\begin{array}{cc}-\nu \Delta \mathbf{u}+\nabla p=\mathbf{f}\hfill & \mathrm{in}\Omega \hfill \\ \nabla ·\mathbf{u}=0\hfill & \mathrm{in}\Omega \hfill \\ \mathbf{u}={\mathbf{u}}_{in}\hfill & \mathrm{on}\partial {\Omega }_{in}\hfill \\ \mathbf{u}=\mathbf{0}\hfill & \mathrm{on}\partial {\Omega }_{obs}\hfill \end{array}\right.$$

(1)

where $\mathbf{u}$ is the velocity field, p is the pressure field, $\mathbf{f}$ is the body force, $\nu$ is the viscosity and ${\mathbf{u}}_{in}$ is the inlet velocity. The inlet velocity is set to ${\mathbf{u}}_{in}=(1,0,0)$ and the viscosity is set to $\nu =0.1$. The body force is set to $\mathbf{f}=(0,0,0)$.

2.1.2. Linear Elasticity Problem

We consider the structural mechanics problem of an elastic plate, fixed on both sides, deformed by the action of two cylindrical actuators. The domain $\Omega$ is a plate of size $(L\times h\times 1)$. The boundary $\partial \Omega$ is divided into three parts: $\partial {\Omega }_{load}$, $\partial {\Omega }_{walls}$ and $\partial {\Omega }_{border}$, corresponding to the loaded cylinders, the fixed walls and the rest of the borders of the plate . The problem is defined as follows:

$$\left\{\begin{array}{cc}-\nabla ·\sigma \left(\mathbf{u}\right)=\mathbf{f}\hfill & \mathrm{in}\Omega \hfill \\ \sigma \left(\mathbf{u}\right)=\lambda \mathrm{tr}\left(ϵ\right(\mathbf{u}\left)\right)\mathbf{I}+2\mu ϵ\left(\mathbf{u}\right)\hfill & \mathrm{in}\Omega \hfill \\ ϵ\left(\mathbf{u}\right)={\displaystyle \frac{1}{2}}\left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right)\hfill & \mathrm{in}\Omega \hfill \\ \mathbf{u}=0\hfill & \mathrm{on}\partial {\Omega }_{walls}\hfill \\ \sigma ·\mathbf{n}=\mathbf{T}\hfill & \mathrm{on}\partial {\Omega }_{load}\hfill \end{array}\right.$$

(2)

The load is set to $\mathbf{T}=(0,-10,0)$, the Young modulus to $21.5\times {10}^4\mathrm{Pa}$ and the Poisson ratio to $0.29$. Body force is set to $\mathbf{f}=(0,0,0)$.

2.2. Dataset

The data set is composed of 500 geometries of size $(L\times h\times 1)$ with 2 circular obstacles of radius $(r_1,r_2)$ inside (Figure 2). The geometries are generated by randomly setting the box dimensions and the obstacles radius around control values.

The data sets are divided into 375 training samples, 75 validation samples, and 50 test samples. The training and validation samples are used to train the model, whereas the test samples are used to evaluate the model’s performance. Each sample is numbered from 0 to 499. Most illustrative figures will use sample 106 from the test set.

2.3. Mesh Parameters

CAD files are automatically generated under .geo_unrolled extension using a python script. In this work, the mesh parameters are limited to the mesh size defined at each point of the CAD file. Mesh size is set randomly around $0.1$ at the circular obstacles while the one at the box extremities is set accordingly to their proximity to an obstacle, in order to test the network under real-case mesh parametrization, rather than a uniform one. Figure 3 highlights in red the 12 points representing each CAD file and labelled with a mesh size used to generate the initial mesh. Mesh generation is performed using Gmsh [31] software.

2.4. Mesh-Based Simulations

Mesh-based simulations are performed using the open-source FE library FreeFem++ [32]. For the Stokes problem, it outputs the velocity field $\mathbf{u}$ and the pressure field p on the mesh (Figure a). This library uses an Hessian-based metric to perform anisotropic mesh adaptation. Given a triangulation ${\Omega }_h$, one can derive an upper bound of the approximation error using an interpolation error analysis. This upper bound is expressed thanks to the recovered Hessian of the approximated solution $u_h$[9]. In fact, using $\mathcal{P}1$ linear elements, we usually cannot directly compute the Hessian of the solution. Instead, we compute an approximation called the recovered Hessian matrix $H_R$($u_h\left(x\right)$) [10].

The recovered Hessian matrix is not a metric because it is not positive definite. Therefore, we define the following metric tensor:

$$\mathcal{M}=\mathcal{R}\Lambda {\mathcal{R}}^T$$

(3)

where $\mathcal{R}$ is the orthogonal matrix built with the eigenvectors $(e_1,e_2,e_3)$ of $H_R$($u_h\left(x\right)$) and $\Lambda =\mathit{diag}\left(\right|{\lambda }_1|,|{\lambda }_2|,|{\lambda }_3\left|\right)$ is the diagonal matrix of absolute value of the eigenvalues of $H_R$($u_h\left(x\right)$). By construction in FreeFem++ and Mmg[33], the prescribed size is the inverse of the square root of the metric eigenvalues, so ${\displaystyle \frac{1}{{\lambda }_1^2}}$ is the size imposed in the first direction, ${\displaystyle \frac{1}{{\lambda }_2^2}}$, in the second direction, and ${\displaystyle \frac{1}{{\lambda }_3^2}}$ in the third direction. Ultimately, $\mathcal{M}$ is symmetric positive definite and is characterised by only six components: $(m_{11},m_{21},m_{22},m_{31},m_{32},m_{33})$.

A similar set of figures can be found in Appendix A.1 for the linear elasticity problem.

2.5. Graph Representation

CAD files and unstructured meshes are essentially a list of points and edges. We can take advantage of this structure to encode these input data into graphs. Each node and edge of the graph will carry an initial information (based on coordinates, boundary conditions, and initial conditions) and the target to predict (either mesh size or physical field). For each data, we can access:

• Geometry: 3D coordinates
• Topology: node connections
• Boundary Conditions (BC): inlet, outlet, walls, obstacles or fluid
• Initial Conditions (IC): velocity and pressure field (mesh-based simulations only)

A graph is denoted by $\mathcal{G}=(\mathcal{V},\mathcal{E})$ where $\mathcal{V}$ and $\mathcal{E}$ are vertex and edge sets, respectively. Each node vertex i is represented by a graph vertex $i\in \mathcal{V}$. Two nodes i and j connected are represented by two graph edges, with edge $(i,j)\in \mathcal{E}$ pointing to node j and edge $(j,i)$ pointing to node i. In other words, the mesh nodes become graph vertices $\mathcal{V}$, and mesh edges become bidirectional mesh-edges $\mathcal{E}$ in the graph. Properties associated with node i are termed as node feature, ${\mathbf{x}}_i^n:\mathcal{V}⟶{\mathbb{R}}^{n_N}$. The properties corresponding to edge $(i,j)$ are referred to as edge feature, $\mathbf{f}(i,j):\mathcal{E}⟶{\mathbb{R}}^{n_E}$.

For both networks, node features ${\mathbf{x}}_i^n$ consists of the belonging to a physical surface (BC) ${\mathbf{n}}_i$, a one-hot vector of dimension 5, and, for Graphnet, the velocity field, a vector of dimension 3, such that the dimension of the node features is $n_N^{Meshnet}=5$ and $n_N^{Graphnet}=8$. The edge feature $\mathbf{f}(i,j)$ of edge $(i,j)$ is constructed to enrich the graph connectivity information with the (signed) distance between nodes, ${\mathbf{u}}_j-{\mathbf{u}}_i\in {\mathbb{R}}^3$, and its absolute value, such that $n_E=4$. Figure 5 and Figure 6 illustrate this process for Meshnet and Graphnet, respectively.

2.6. Model

In both Meshnet and Graphnet, the behaviour of the targeted variables are captured by learning directly on the graph using latent node and edge features derived from the physical node and edge features reviewed in Section 2.5. We use a graph neural network model with an Encoder-Processer-Decoder architecture [29,34]. See Appendix A.2 for more details.

The architecture proposed by Pfaff et al. [29] is essentially designed to predict the dynamical states of a system over time. Since our test cases are steady-state predictions, we adapt this architecture following the work of Harsch et al. [35] to better capture local and global features.

2.6.1. Encoder

Using the graph representation of Section 2.5, we compute the initial latent feature vectors, ${\mathbf{v}}_i^0\in {\mathbb{R}}^{n_H}$ and ${\mathbf{e}}_{(i,j)}^0\in {\mathbb{R}}^{n_H}$ from the physical feature vectors, ${\mathbf{x}}_n^i$ and $\mathbf{f}(i,j)$. The hyperparameter $n_H$ denotes the size of the latent vectors. The computation of the latent vectors is done using the node and edge multilayer perceptrons (MLPs) (Table A1), denoted respectively by ${\mathrm{mlp}}_v^0$ and ${\mathrm{mlp}}_e^0$, as follows:

$$\left\{\begin{array}{c}{\mathbf{v}}_i^0={\mathrm{mlp}}_v^0\left({\mathbf{x}}_n^i\right)\hfill \\ {\mathbf{e}}_{(i,j)}^0={\mathrm{mlp}}_e^0\left(\mathbf{f}(i,j)\right)\hfill \end{array}\right.$$

(4)

2.6.2. Processor

The processor consists of m message-passing steps computed in sequence. At step $l\in \{1,\cdots ,m\}$ in the sequence, each graph edge feature ${\mathbf{e}}_{(i,j)}^l$ is updated using its value at the previous message-passing step $l-1$ and the values of the adjacent node features at step $l-1$, as follows

$${\mathbf{e}}_{(i,j)}^l={\mathrm{mlp}}_e^l\left(\left[{\mathbf{e}}_{(i,j)}^{l-1},{\mathbf{v}}_i^{l-1},{\mathbf{v}}_j^{l-1}\right]\right),l\in \{1,\cdots ,m\}$$

(5)

to obtain the updated value. In Equation (5), the operator $\left[·\right]$ concatenates the given arguments on the feature dimension. Then, each graph node ${\mathbf{v}}_i^l$ is updated using its value at the previous message-passing step, $l-1$, and the aggregation of its incident edge features at step $l-1$:

$${\mathbf{v}}_i^l={\mathrm{mlp}}_v^l\left(\left[{\mathbf{v}}_i^{l-1},\sum _{j\in \mathrm{adj}\left(i\right)}{\mathbf{e}}_{(i,j)}^l\right]\right),l\in \{1,\cdots ,m\}$$

(6)

where $\mathrm{adj}\left(i\right)$ is the set of nodes connected to node i. Using Equations (5) and (6), the processor computes an updated set of node features.

The update of the edge-based messages, ${\mathbf{e}}_{(i,j)}^l$, is the key to the accuracy of the flow predictions as it propagates information between neighbouring graph nodes (i.e., between neighbouring mesh cell elements). This design choice differentiates MeshGraphNets [29] from other classical GNN frameworks relying only on node features, such as GCN and GraphSAGE [36].

2.6.3. Decoder

After performing m steps of the previously described process, the decoder maps the updated graph nodes latent features ${\mathbf{v}}_i^l$ to the node-based properties in physical state using MLP as follows:

$${\widehat{y}}_i={\mathrm{mlp}}^{m+1}\left({\mathbf{v}}_i^l\right)$$

(7)

where m is the number of steps performed by the processor.

2.6.4. Local and Global Features

To increase the ability to extract local and global properties, we adopt the structure proposed by Wang et al. [37] for Dynamic Graph CNN and adapted by Harsch et al. [35]. This adaptation aggregates the information given by each message passing block to calculate a global feature vector. In that sense, the number of layers is reduced to 5, instead of 15 previously. To avoid potential overfitting, a pooling operation is added, and the result is concatenated with the previous local features. This adapted version is called Graphnet (Pool) (Figure 8).

2.7. Loss Function and Optimiser

We train both networks by minimising the error between the true node label and the predicted node label. We use the per-node root mean square error (RMSE) loss to quantify the data error for each simulation. The loss function reads:

$$\mathcal{L}=\sqrt{{\displaystyle \frac{1}{n_B}}{\sum }_{i=1}^{n_B}{\parallel y_i-{\widehat{y}}_i\parallel }_2^2}$$

(8)

where $n_B$ denotes the number of nodes in a batch of training meshes, $y_i$ denotes the true output in the data set, ${\widehat{y}}_i$ is the output predicted by the network, as formalised in Equation (7).

ADAM optimiser is used and the learning rate is constant, equal to $0.001$. With the data set and a fixed set of hyper-parameters at hand, each epoch takes $30s$ on two NVIDIA A100 GPUs. The value of specific training hyperparameters is given in Table A3.

3. Results

3.1. Meshnet

3.1.1. Fine Tuning Input Data for Message-Passing GNN

Meshnet was trained to learn to predict the mesh size at each point of a given geometry. These mesh sizes are then used to generate a first mesh using Gmsh [31]. The loss curve highlights the message-passing nature of this GNN. However, when we first trained the network, the case did not converge, giving the blue cyan loss curve on Figure 9a. During the graph-encoding step (Figure 5), we pass to the network the distance of each edge of the geometry in the edge feature matrix. However, in our geometrical representation, the points of the box are not linked to the obstacles (Figure 9b). Thus, the network is blind and cannot learn this specific mapping. As a workaround, we manually added controlled edges so that the new edge feature will keep the distance between the points of the box and the obstacles (Figure 9b).

3.1.2. Evaluating Meshes Similarity and Quality

To evaluate the network, we are also interested in the mesh similarities between the mesh generated with the ground truth parameters, and the one generated by the predictions. The comparison of the number of points and tetrahedra between the ground truth and the prediction can only give us a small overview (Figure 10a).

We can take, for instance, the cases of simulations 106 and 130. The former shows a good similarity in terms of points and tetrahedra, while the latter has a prediction with twice as many points and tetrahedra than the ground truth. To better assess the quality of each mesh, we use "Glassmaier" parameter [38] defined as

$$Q_G=1+{\displaystyle \frac{\mathrm{True}\mathrm{Surf}.}{\mathrm{Ideal}\mathrm{Surf}.}}+{\displaystyle \frac{\mathrm{True}\mathrm{Vol}.}{\mathrm{Ideal}\mathrm{Vol}.}}$$

(9)

and takes on values between 1 and 3. It tends to describe the dimensionality of the figure, as listed in Table 1. The ideal volume and surface are calculated for a regular tetrahedron with a side length equal to the average of the 6 distances between the 4 points.

The distribution of this quality estimator among the ground truth and predicted meshes are very similar for both cases (Figure 10b), and the visual observation confirm this conclusion (Figure 11). In fact, the number of points and tetrahedra is greatly sensitive to the the mesh size defined on the obstacles, since it is the smallest one, thus giving the greatest number of points of tetrahedra. Thus, in the case of simulation 130, the ground truth and the prediction have a similar pattern of mesh, but one is finer at the obstacles.

3.2. Graphnet

3.2.1. Direct Predictions: Velocity and Anisotropic Metric Fields

Graphnet was trained to learn mesh-based simulations. First, we trained one network to predict the three components of the velocity field, in order to use it to compute an anisotropic metric from the Hessian of the velocity field (Figure 12a). Then, we decided to predict directly the metric field by training one network for each of the six components of the metric tensor(Figure 12b). Finally, we use the adapted architecture Graphnet (Section 2.6.4) to improve the training (Figure 12c).

The case of simulation 106 highlights that the velocity field, in addition to be a 3-component field, is a more complex learning task than the metric field. As expected, error is concentrated around the obstacles, where the gradients are higher and the boundary conditions changing (Figure 13). The metric field prediction still shows some error, but catches precisely the evolution of this variable across the simulation domain (Figure 14).

3.2.2. Mesh Adaptation

Ultimately, we are looking for the resulting adapted mesh generated through this metric field. We use Mmg software [33] to realise such adaptation in 3D. Again, generated meshes show great visual similarity (Figure 15), confirmed by quality analysis (Figure 16b). As expected, the number of points and tetrahedra has been divided by two through the AMR process, both for ground truth and prediction in the case of simulation 106.

This result can be extended to the remaining simulations of the test sets. The relative errors in number of points and tetrahedra (Figure 16a) show great improvement in comparison to Meshnet (Figure 10a), meaning that the gap in points an tetrahedra created by the prediction of the initial mesh has very little influence on the ultimate adapted mesh. Mesh quality is more important, as it is almost conserved through the AMR process (Figure 16b).

3.3. Adaptnet

Adaptnet is the framework combining Meshnet and Graphnet (Figure 1. Given an initial geometry, it will predict initial mesh sizes (Meshnet) to produce an initial mesh and predict on this mesh an metric field to adapt it to the Stokes problem. This process has already been illustrated through the example of the 3-hole geometry. First, we showed how Meshnet produces an initial mesh (Figure 17), and then, we use this mesh to predict the anisotropic metric tensor at each point of the simulation domain and adapt the mesh (Figure 18). By predicting directly the metric tensor, we are able to skip the calculation of the velocity field, unlike a FE solver, speeding up calculation time accordingly (Table 3). For Graphnet, we compare the prediction time of the six models (one for each component of the metric tensor) with the computation time of FreeFem++, on the 3-hole test case, using a single CPU. For Meshnet, the prediction time should be compared with the engineering time needed to choose the mesh parameters.

4. Discussion

4.1. Meshnet Generalisation

The architectural choice of using relative encoding on graphs has been shown to be very conducive to generalisation [34]. Our model generalises well, as can be seen when we want it to predict the mesh size at each point on a 3-hole geometry (Figure 17).

4.2. Graphnet Generalisation

The same observation can be made when we want to predict the adaptation metric in a 3-hole geometry (Figure 18).

Finally, we tested our model with a different geometry from those seen previously. This time, we chose a square (Figure 20). The point here is not to test our model with every possible shape but simply to highlight its ability to interpolate in a wide spectrum of graph structures.

4.3. Linear Elasticity Problem

The attempt to train the graph neural network (GNN) on the linear elasticity problem proved to be less successful than expected (Figure 22), highlighting the challenges inherent in this complex task. Despite meticulous parameter tuning and a comprehensive dataset, the network struggled to effectively capture and generalise the underlying patterns within the graph structures. The intricate interplay of nodes and edges posed a formidable challenge, and the model exhibited difficulties in discerning meaningful relationships and dependencies.

These results underscore the need for further investigation into refining GNN architectures and training methodologies, emphasising the intricate nature of graph-based data and the nuances associated with their representation in neural networks.

At first glance, the quality of the mesh produced may seem even better than the reference (Figure 23). But quantitative analysis highlights the poorer quality of prediction in this case (Figure 24). Still, it can be noticed that the plate is very thin along z - axis. As a consequence, the tetrahedra generated will necessarily be stretched in this direction, resulting in lower quality. But the overall quality distribution is more satisfactory, as the first visual impression suggested.

5. Conclusions

We propose a Graph Neural Network mesh-based method to model fluid dynamics and solid mechanics problems for an accurate and efficient prediction of anisotropic adaptation metric fields. We extended the approach of Pfaff et al. [29] to tetrahedral meshes in three dimensions. Following the work of Harsch et al. [35], we adapted the GNN architecture to better suit the prediction of steady states. Our method may allow for more efficient simulations than traditional simulators, and because it is differentiable, it could be used to retrieve the Hessian matrix directly, when traditional solvers have to compute a recovered Hessian matrix.

The experiments demonstrate the model’s strong understanding of the geometric structure. The method is independent of the structure of the simulation domain, being able to capture highly irregular meshes. We show that the model does not require any a priori domain information, e.g. inflow velocity or material parameters. Thus, the model can be used for any other systems represented as field data.

Training phase underlined the model’s strong understanding of the geometric structure of the unstructured meshes. It was shown that it can achieve effective prediction with the sole knowledge of connectivity and the belonging to a physical surface. It doesn’t rely on any prior information such as the inflow velocity or the material parameters. This, it can be easily adapted and scaled to other physical problems governed by PDEs. However, our test show that some specific tuning might be necessary to achieve correct predictions when adapting from a fluid problem to a solid one for instance.

The strength of this work lies in the tetrahedral mesh generation and adaptation pipeline. The emphasis was placed on producing work based on open source tools [31,32,33] to promote reproducibility and enable others to build on it. Future work in this area should address several topics. It could explore more complex geometries to prove its applicability to industrial problems. Alternative architectures could also be explored. Adding concatenation and pooling has been shown to improve performance, but other options could show even better improvements. Finally, this framework could be further improved by adding a specific physical constraints to the loss function, similar to Physics-Informed Neural Networks (PINNs) [39], which could enable learning on sparse datasets, reducing the need for generating hundreds of simulations.

Appendix A.1. Linear elasticity problem

Figure A1. Finite Element solving of Elasticity problem and anisotropic metric computation for sample 106.

Figure A1. Finite Element solving of Elasticity problem and anisotropic metric computation for sample 106.

Appendix A.2. Model

The MLPs of the Encoder, Processor, and Decoder are ReLUactivated two-hidden-layer MLPs with layer and output size of 128, except for the Decoder MLP whose output size matches the prediction. All MLPs outputs except the Decoder one are normalized by a LayerNorm. All input and target features are normalised to zero mean unit variance, using dataset statistics.

Table A1. Details of node and edge encoders.

Node MLP (${\mathbf{mlp}}_{\mathit{v}})$	Edge MLP (${\mathbf{mlp}}_{\mathit{e}})$
Input: x	Input: x
$x=\mathrm{Linear}(5/8,128)\left(x\right)$	$x=\mathrm{Linear}(3,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{LayerNorm}\left(128\right)\left(x\right)$	$x=\mathrm{LayerNorm}\left(128\right)\left(x\right)$

Table A1. Details of node and edge encoders.

Node MLP (${\mathbf{mlp}}_{\mathit{v}})$	Edge MLP (${\mathbf{mlp}}_{\mathit{e}})$
Input: x	Input: x
$x=\mathrm{Linear}(5/8,128)\left(x\right)$	$x=\mathrm{Linear}(3,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{RELU}\left(x\right)$	$x=\mathrm{RELU}\left(x\right)$
$x=\mathrm{Linear}(128,128)\left(x\right)$	$x=\mathrm{Linear}(128,128)\left(x\right)$
$x=\mathrm{LayerNorm}\left(128\right)\left(x\right)$	$x=\mathrm{LayerNorm}\left(128\right)\left(x\right)$

Table A2. Hyperparameters used for processor.

Parameter name	Value
Number of GNN layer for the processor	15
Latent size for the processor	128
Activation	ReLu
Type of normalization	Layer normalization
Input feature size of Node MLP	256
Input feature size of Edge MLP	384

Table A2. Hyperparameters used for processor.

Parameter name	Value
Number of GNN layer for the processor	15
Latent size for the processor	128
Activation	ReLu
Type of normalization	Layer normalization
Input feature size of Node MLP	256
Input feature size of Edge MLP	384