## 2 Materials and Methods
### Parallel Ray Generation  
The starting point distribution of rays is generated according to the calculation task, which is determined by the mesh size and directly affects the quality of mesh generation. Because of the large number of rays, it is more efficient to generate starting points directly in the device memory than the host memory and then to transfer them to the device memory. In this subsection, the mesh size is determined by optimizing algorithm, and then a method of generating ray coordinates in a single GPU thread is proposed. Supposing that the computational domain is a rectangular region. The computational domain can be fixed by two points:\(p_{1}(x_{1},y_{1},z_{1})\) and \(p_{2}(x_{2},y_{2},z_{2})\). There are some geometries in the domain, and all the geometries in the domain will be meshed in order. The first step is to calculate the bounding box of each geometry. The bounding box of \(G_{i}\) can be defined by a minimum point \(p_{min}\) and a maximum point \(p_{max}\). The mesh size of the uniform Cartesian mesh equals to a constant value \(\mathbf{s_{opt}}\) in three dimensions of the whole computational domain, respectively. Taking the \(X\) dimension as an example, the \(X\) dimension is divided into several intervals by bounding boxes of geometries, and the interval boundaries are defined by \(X_{i}\) and \(n_{x}\) is the number of endpoints of intervals. Two restrictions must be satisfied [30]: (a) endpoints of intervals have to be part of the mesh lines; (b) the cell size of each geometry \(G_{i}\) does not exceed a maximum size \(s_{max}\). According to the geometric characteristics and calculation requirements, \(s_{max}\) can be selected artificially. Firstly, the maximum mesh cell size \(s_{opt}\) which satisfies the above restrictions must be calculated. To obtain the optimized mesh cell size, the cell size function in \(X\), \(Y\) and \(Z\) dimension can be expressed as Equations (1)-(3), and the optimized cell size can be calculated through minimizing these three functions, respectively. \[f_{x}(s_{x,opt})=\sum_{i=1}^{n_{x}-1}(\frac{\Delta X_{i}}{s_{x,opt}}) \tag{1}\]  
\[f_{y}(s_{y,opt})=\sum_{j=1}^{n_{y}-1}(\frac{\Delta Y_{j}}{s_{y,opt}}) \tag{2}\]  
\[f_{z}(s_{z,opt})=\sum_{k=1}^{n_{z}-1}(\frac{\Delta Z_{k}}{s_{z,opt}}) \tag{3}\]  
with  
\[\Delta X_{i}=X_{i+1}-X_{i}\quad i=1,...,n_{x}-1 \tag{4}\]  
\[\Delta Y_{j}=Y_{j+1}-Y_{j}\quad j=1,...,n_{y}-1 \tag{5}\]  
\[\Delta Z_{k}=Z_{k+1}-Z_{k}\quad k=1,...,n_{z}-1 \tag{6}\]  
Then, through the minimum point of bounding box and the optimized mesh cell size \(s_{opt}\), the coordinates of mesh lines can be calculated with the Equations (7)-(9). Once the mesh lines are obtained, the start point of rays can be calculated. The process of ray and mesh line generation can be illustrated in Figure 3.  
\[x_{i}=x_{min}+i\times s_{x,opt}\quad i=0,...,n_{x}-1 \tag{7}\]  
\[y_{j}=y_{min}+j\times s_{y,opt}\quad j=0,...,n_{y}-1 \tag{8}\]  
\[z_{k}=z_{min}+k\times s_{z,opt}\quad k=0,...,n_{z}-1 \tag{9}\]  
When a parallel program runs on the GPU, the initial data should be transferred from the host memory to the device memory. Thus, the less data transmission occurs in the program, the more efficient the program is. To reduce the data transmission, ray generation, intersection calculation and  
\[f_{y}(s_{y,opt})=\sum_{j=1}^{n_{y}-1}(\frac{\Delta Y_{j}}{s_{y,opt}}) \tag{2}\]0 Cartesian mesh generation baseline, (**a**) computational domain; (**b**) a ray passing through the computational domain; (**c**) result of mesh generation. property mapping are all executed on the GPU. Basically, only two data transmissions, i.e., the initial computational domain information transferred from host to device and the mesh information transferred from device to host, are necessary. The initial computational domain information contains: the coordinates of two points representing the bounding box, the coordinate array and normal vector array of triangular facets, and the optimized mesh cell size \(\mathbf{s_{opt}}\).