*2.4. Node Generation Algorithm*

Benito et al. (2017) pointed out that when the discretization (cloud of nodes) becomes more regular, the obtained results become more logical [21]. Furthermore, to avoid the diffraction caused by the difference between the node and the interface positions, it is necessary to adjust the node position at the inner part of the interface and the boundary, so that the nodes are located at the interface and the boundary. In this study, we proposed an algorithm combining the Poisson disk node generation algorithm and the centroid Voronoi adjustment algorithm, to obtain a reasonable node distribution according to the velocity model.

### (1) Randomly Distributed Node Generation

The Poisson disk node generation algorithm proposed by Fornberg et al. (2015) was used [28]. In the two-dimensional case, the calculation steps for the random node generation algorithm (Figure 2) are as follows:


**Figure 2.** Generation process for the Poisson disk algorithm [28]. The red solid point is the DDP, and the black hollow points are PDPs.

## (2) Node Adjustment

To adjust the nodes, overlay the random nodes obtained in the previous step with the interface and boundary nodes, and discard the nodes within a half radius of any boundary node, as shown in Figure 3b. After the redundant nodes are discarded, the nodes at the internal interface and the boundary will be unevenly distributed. Therefore, we propose a node adjustment algorithm based on the centroid Voronoi structure. This algorithm divides Voronoi polygons into multiple polygonal regions formed by the boundary and interfaces, and uses the centroid point of the Voronoi polygon as the new node location. Since the distribution of the nodes before adjustment is relatively uniform, only 10–15 iterations are needed to evenly distribute the nodes near the interface. This method has good stability and can obtain a more uniform node distribution, however, since the conventional Voronoi polygon division is borderless, the Voronoi polygon division requires computing the intersection of the Voronoi diagram and the irregular boundary, which leads to a high computation cost. To solve this problem, GPU was used to accelerate the Voronoi polygon division and the adjustment of the centroid points in the irregular boundary area.

**Figure 3.** Generation steps for the Voronoi element under the restriction of an irregular boundary (step (a) and step (b)) [30]. X*<sup>i</sup>* is the point to be adjusted, and X1-X5 are the five points surrounding X*i*.

The algorithm takes the internal node and boundary node coordinates as its inputs, and its main steps are as follows:


Figure 4 shows the node generation steps for the two-layer simple model. It can be seen from the figure that the algorithm first generates variable density nodes through Poisson disks in the boundary region (Figure 4a). Then, the nodes that are too close to the interface nodes are deleted, and nodes with an uneven distribution at the boundary are obtained (Figure 4c). Finally, variable density internal nodes, and a uniform distribution of boundary nodes, are obtained through the centroid Voronoi adjustment algorithm (Figure 4d). The figure shows not only that the final node distribution is accurately arranged at the interface location, but also that the overall nodes layout is evenly distributed across the different regions.

**Figure 4.** Node generation of two-layer model: (**a**) Variable density nodes; (**b**) Variable density nodes overlap with boundary nodes and interface nodes; (**c**) Delete the nodes in (**b**) that are too close to the interface or boundary; (**d**) Nodes after adjustment based on Voronoi centroid algorithm.
