An RBF-Based h-Adaptive Cartesian Grid Refinement Method for Arbitrary Single/Multi-Body Hull Modeling and Reconstruction

: Complex single/multi-body structures are generally found in ship and ocean engineering. They have the smooth, sharp, concave, and convex surface features in common. Precise modeling of the structures is the basis of numerical simulation. However, the most widely used explicit modeling method requires considerable manual operations. The result is also difﬁcult to reproduce. Therefore, this paper presents a Radial basis function (RBF) based hierarchical (h-) adaptive Cartesian grid method. The RBF is introduced for arbitrary implicit modeling over the Cartesian framework. Thanks to its natural properties, the RBF is easy to use, highly automated, and only needs a set of scatter points for modeling. The block-based h-adaptive mesh reﬁnement (AMR) combined with the RBF aims to enhance the local grid resolution. It accelerates the calculation of intersecting points compared with the uniform Cartesian grid. The accuracy, efﬁciency, and robustness of the proposed method are validated by the simulation of the 3D analytical ellipsoidal surface and the unclosed conic surface. To select suitable parameters, we thoroughly investigated the uncertainty factors including sample points, RBF functions, and h-AMR grids. The simulation results of the single/multi-body Wigley hull and KCS hull forms veriﬁed the proper selection of the factors and the feasibility of our method to solve practical problems.


Introduction
Complex geometries, such as various types of single/multi-body hull forms, are frequently found in ship and ocean engineering. To conduct the numerical study on such problems, accurately modeling the arbitrary geometries is the first step, and the key process is to perform the robust and efficient reconstruction of complex surfaces within the simulation. Additionally, the practical marine structures are three-dimensional (3D) and usually have the smooth, sharp, concave, and convex features in common. That makes the modeling and reconstruction more challenging.
Generally, the numerical simulation of the geometric problems is performed by the mesh-based CFD approaches, which can be divided into the body-fitted and the non-bodyfitted methods. The former, including the separate C-/O-/H-mesh and the hybrid mesh [1], requires the mesh to align with the solid surface. Thus, it can accurately represent the boundary, but is hard to generate the grid for the complex geometries. By contrast, the non-body-fitted method uses the regular background Cartesian grid. That gives the method the advantage of simple data structure, fast grid generation, and easy local refinement [2,3]. The cut cell method is one of the most popular non-body-fitted methods. Compared with the diffused interface methods, such as the immersed boundary method (IBM), it can represent the interface sharply. However, to exploit the benefits of the cut cell method, there are two challenges needed to be concerned. First, the simulation of the complex marine Figure 1 shows the main process of this work. It can be divided into three main steps: (1) initial grid generation and the geometric refinement of the h-AMR grids; (2) single/multi-body RBF geometric modeling through scatter points; and (3) the calculation of the intersecting points based on the combined RBF and h-AMR background grid. Specifically, the initial h-AMR Cartesian grid is generated, and its data structure is adopted throughout the work process. Meantime, the RBF implicit surface is automatically constructed by a given set of sample points. According to the characteristics of the marine structure, this work also includes the detailed comparisons of effects of uncertainty factors, which include different RBF functions and parameters, on the accuracy of results. Moreover, the RBF and h-AMR grids are combined to reconstruct the implicit isosurface and improve efficiency. We directly use the RBF equation to test whether the grid nodes of the candidate h-AMR grids are inside/outside/on the solid surface. Once the cut cells have been identified, they are subsequently refined by the h-AMR technique. Thereafter, Symmetry 2021, 13, 15 4 of 29 this work combines the Newton-Raphson iterative method. Then the intersections can be iteratively calculated based on the RBF value over the h-AMR cut cells.

Spatial Discretization of Block-Based H-Adaptive Cartesian Grids
The uniform grid maintains the same grid size throughout the domain, while the h-AMR method adjusts different grid sizes in different regions. Figure 2 represents the uniform Cartesian grids (left) and the geometric h-AMR grids (right) for a 2D rectangle (blue). It can be seen that the AMR grids can reduce the total amount of grids while reaching the same grid resolution at the surface boundary. The existence of multi-layered structure grids means the h-AMR grids cannot be directly indexed by the implicitly connected coordinates (i, j, k) in the uniform Cartesian grid storage. Therefore, the unstructured hierarchical tree structure needs to be adopted to store the connectivity and information on the grid of different levels.

Characteristics of Different h-AMR Grids
For the specific use, the proper h-AMR method needs to be selected. Depending on the grid forms stored in data tree nodes, the h-AMR can be further divided into three typical methods, the patch-based [4,26,28], cell-based [2,25,29] and block-based [30][31][32] method. Figure 3 represents the geometric refinement process of the 2D rectangle by the three typical h-AMR grid methods. Each group of grids surrounded by yellow lines forms a data tree node. In Figure 3a, the patch-based adaptive grid consists of arbitrary rectangular patches. It needs the clustering algorithm to organize the marked cells into grid patches, and each grid patch forms a tree node. Compared with other h-AMR methods, it generates the fewest Symmetry 2021, 13, 15 5 of 29 number of tree nodes. However, it may not be easy to apply the appropriate clustering algorithm [32].
The cell-based methods in Figure 3b, on the other hand, directly bisect the marked cells into 2 n dim (n dim is the space dimension) finer children cells. That makes it more flexible in refining arbitrary tagged grids. It does not need the use of clustering methods. However, due to the one-to-one correspondence between grid and tree nodes, the hierarchical data tree requires more storage capacity than a patch-based approach.
The block-based method, as illustrated in Figure 3c, has the same characteristics of the other two h-AMR methods. It includes both the grid block (heavy lines) and grid cell structures (light lines). The grid blocks are used for tagging, refining and forming tree nodes. Meanwhile, the grid cells are uniformly arranged in each grid block. The blockbased h-AMR avoids the clustering procedure compared with the patch-based method. In the meanwhile, it not only retains the flexibility of the cell-based method but also partially preserves the structured grid characteristics. Therefore, this work applies the block-based h-AMR method [30]. That means the grid generation, geometric modeling, and reconstruction techniques must be fully compatible with the characteristics of the block-based h-AMR. Figure 4 shows a simple 2D block-based h-AMR grid and its data tree. In Figure 4a, the heavy lines and the light lines present the grid blocks and the grid cells, separately. After the refining procedure, the grid blocks need to be reordered, as shown in the red lines in Figure 4a. Each block is stored in the matching tree node and has a unique indexing identifier, as shown in Figure 4b. The domain (root block) forms the first tree level. The blocks without any children named leaf blocks perform the refinement and calculation. For example, in Figure 4, block 1 is the root block as well as the first data tree node located at Level 1. Block 2 (yellow) is one child of block 1, and also the parent of block 3 (purple). Block 10 (green) in Level 4 and block 3 (green) in Level 3 are both leaf blocks. The grid cells are independent of the tree structure, and only participate in the practical calculations. Figure 5 represents the index of grid cells in the grid block. Each grid block contains n cx × n cy = 4 × 4 interior structured Cartesian grid cells. They can be indexed by the implicitly connected Cartesian coordinate (i, j), where i = 1 . . . n cx and j = 1 . . . n cy refer to the x-and y-directions, respectively. When the number of grid cells in each block is one, the block-based method is similar to the cell-based AMR method. Before the simulation, the root block needs to be discretized into smaller uniform grid blocks. The initial uniform block-based h-AMR Cartesian grids are mainly controlled by three parameters, (1) the size of the computational domain (root block), (2) the coarsest refinement level l min of the root block, and (3) the number of the grid cells in each grid block. The initial domain contains 2 n dim (l min −1) sub-grid blocks. In this work, the computational domain size should be determined by the given scatter points. The geometric h-AMR grid layout will be further controlled by the user-defined AMR level l AMR .

Theory of RBF Isosurface Equation
The RBF technique uses a set of sample points to implicitly define the surface Γ based on basis functions. This makes it feasible to model both single and multi-body geometries. The bodies in the multi-body system can be separately reconstructed by their own sample points in the same domain, as shown in Figure 6. Hence, the main operation is to construct the continuous and differentiable RBF isosurface equation for each object. The surface Γ is represented by the RBF isosurface equation. Given a set of 3D distinct scatter points {p i } ⊂ R 3 , i = 1, 2, . . . , N (also known as RBF centers), the RBF equation can be expressed as follows: where N is the number of sample points. φ : [0, ∞) → R 3 is the real-valued basis function, and it has various types. p − p i indicates the Euclidean distance from the space point p(x, y, z) to the given RBF center p i (x i , y i , z i ), as shown in Equation (2). λ i is the weight coefficient to be calculated, which needs to satisfy the orthogonality in Equation (3).

of 29
Once the F(p) = 0 (p ∈ Γ ⊂ R 3 ), Equation (1) becomes the zero isosurface equation which represents the actual geometric surface. Equation (4) is the matrix form of Equation (1). Because of the zero terms on the right-hand side of the equation, it might lead to a trivial solution, which means that the F(p) might be zero everywhere. To avoid this, it is often necessary to add some non-zero off-surface points to the sample points. However, these additional points also increase the calculations.

Non-Zero Distance d
This work utilizes the constant d on the right side of the RBF equation instead of using 0. Therefore, the original number of sample points can be maintained. As shown in Figure 7, the non-zero d represents a fixed distance from the numerical surface to the actual surface. Thereby, we can directly calculate the non-zero scalar field. Then Equation (1) can be transformed as, and the non-zero RBF equation for numerical calculation is represented as, Equation (7) is the matrix notation form of Equation (6). Then the coefficient λ i can be solved by matrix solving method such as the Gaussian elimination.
Afterward, Equation (6) can be reconstructed by the coefficient λ i and the sample points. The actual zero-valued RBF isosurface in Equation (5) can be obtained by subtracting the distance d. In our previous 2D work [17], d = 1 is suggested. This work further validates the correctness and robustness of the distance d in the simulation.

RBF Functions and Shape Parameters
There are various types of functions and parameters that can be selected to construct the RBF isosurface equation. Based on our research field, this work mainly analyzes two kinds of RBFs: multiquadric RBF (MQ-RBF) and Wendland's CSRBF.
The MQ-RBF [9] is a simple GSRBF, as represented in Equation (8), where c is the shape parameter. Frank [10] compared the scattered data interpolation algorithms and pointed out that the MQ-RBF has a better fitting result. It can be used to solve partial Symmetry 2021, 13, 15 8 of 29 differential equations [33] and has been widely used in geologic modeling [6]. It has also been adopted in our previous works [17][18][19].
The capabilities and effects of MQ-RBFs are mainly controlled by the shape parameter c. The c is used as the smoothing factor to make the reconstructed surface smoother. However, if c is too large, it may cause the ill-conditioned matrix and increased numerical errors [34]. Our previous works [17][18][19] use c = 0 in hull modeling. This paper further tests two types of shape parameters c: the constant c, and the variable c of MQ-RBF. Table 1 lists some shape parameters c used in MQ-RBF. Table 1. Some shape parameters c in multiquadric (MQ)-RBF [10,34].

Name
Definition Type In Table 1, the constant c assigns the fixed value to all the columns of the RBF coefficient matrix. The D of Franke's empirical formula [10] is recommended to be the smallest diameter of the bounding sphere containing all data points. For simplicity, this work uses the diagonal of the bounding box L diag instead of D.
On the other hand, the variable c adopts different shape parameters c j for each matrix column according to the formula. In Table 1, the variable c in [35,36] has been applied to solve 2D function interpolations and linear boundary value problems. Here we test and compare their effects on 3D arbitrary modeling. c min and c max can be determined by Equation (9). The function rand() in Random c is employed to produce random numbers. We should notice that the calculated pseudo-random numbers practically rely on programming languages. Table 2 lists Wendland's CSRBFs. The essential influencing factor of the CSRBF is the parameter continuities and the support size r cs .

Order of Continuity Definition
In Table 2, (1 − r) + means the control function of the compactly supported domain, which satisfies Equation (10). Othake [22] suggests using σ s = 0.75 times the diagonal of the bounding box of the sample points as the r cs , as shown in Equation (11). In practice, the selection of σ s should depend on specific problems. Therefore, we analyze the uncertainty of σ s in the simulations.

RBF-Based Geometric h-AMR and Calculation of Intersection
The RBF-based geometric modeling method is actually independent of the block-based h-adaptive Cartesian background grids. In this work, we combined them to efficiently calculate the intersections.

Surface Identification and RBF-Based Geometric AMR
First, the RBF isosurface equation is directly used as the geometric h-AMR criterion. Based on the value F(p) of the RBF isosurface equation, the positional relationship between an arbitrary space point p(x, y, z) and the geometric surface can be determined, as shown in Equation (12) and Figure 8. Since the refinement procedure is presented on the grid blocks, the next step is to find cut blocks. To help subtract the off-surface grid blocks, this work introduces the bounding box. Figure 9 represents the bounding boxes for multi-body hulls and the candidate cut blocks. As illustrated in Figure 9a, the axis-aligned rectangular region Ω is established by the minimum and maximum coordinates of the given sample points. Additionally, the scale factor k s is used to properly expand the range of the object, thus avoiding the incomplete reconstruction of the geometric surface. Figure 9b represents the expanded bounding box Ω . Here we select k s = 1.05 [7] in the simulations. Within the bounding box, the RBF value of each block node (coordinates of grid vertices; 2D: four; 3D: eight) needs to be calculated. As shown in Figure 10, once the RBF value of one node is 0, or any two nodes have opposite signs, we stop calculating the block and mark it to perform h-AMR.

Calculation of Intersections
The next step is to calculate the intersections for geometric reconstruction. The calculation of the intersection is performed on the grid cells inside the h-AMR cut block. The method used to detect cut cells is similar to the method of cutting blocks. The difference lies in their traversal methods. The candidate cut cells are searched by the implicitly Cartesian connection. As shown in Figure 5, if the cell node is (i, j, k), its neighboring cell nodes are (i + 1, j, k), (i, j + 1, k), and (i, j, k + 1), respectively.
Due to the RBF isosurface being a continuous function, there should exist a zero-valued intersection point between the two opposite-valued nodes. Besides, with the sufficiently fine Cartesian background grid, there will be at most one intersection on each edge of the cut cell. Figure 11a shows a 3D cut cell. p 0 is a cell node. p 1 , p 2 , and p 3 are the neighbor cell nodes of p 0 in x-, yand zdirections, respectively. We take the intersections on the XY plan as an example, as shown in Figure 11b. Based on the RBF value, it can be seen that F(p 0 ) > 0, F(p 1 ) = 0, and F(p 2 ) = 0. It indicates that p 0 is the fluid point (p f luid ) outside the surface, p 1 is the solid point (p solid ) inside the surface, and p 2 is an intersection point on the surface. In this situation, the intersection p I between p 0 and p 1 needs to be calculated. Based on the Newton-Raphson iterative method, the weighting coefficient s (s ∈ [0, 1]) is applied to associate p f luid with p solid . When s = 1 (s f luid )/0 (s solid ), the intersection point p I becomes p f luid /p solid , respectively. Then the coordinate of p I can be defined as Equation (13), in which the subscript n represents the n th iterative step.
Subsequently, we take the initial guessing value s1 = 0.5 into Equation (13). This assumes the intersection p I 1 is at the midpoint p mid of the cell edge. Here we use a very small value ε = 10 −6 as the convergence condition. If F(p I 1 ) − ε > 0, Equation (14) is used to solve the sn+1 in the next iterative step. where F (p I n ) can be calculated from Equation (15): Equation (15) contains two conditions. If F(p I n ) > 0, it means p I n is a fluid point p f luid , and the actual intersection p I is closer to the p solid . Therefore, the range of sn+1 reduces to s n+1 ∈ (0, s n ), and the p f luid is replaced by the p I n in the next iterative step. The procedure is repeatedly executed by Equations (13)-(15) until F(p I n ) converges to ε.

Simulation and Results
This section simulates the 3D analytic geometries and practical hull forms. It includes the ellipsoidal surface, unclosed conic surface, and single/multi-body Wigley and KCS hull forms. The goal is to verify the accuracy, efficiency, and robustness of the RBF-based block-based h-AMR method. All simulations in this work use dimensionless parameters and CPU time in seconds. The numerical tests are carried out on an Intel Core i7 CPU 2.6 GHz processor.
In the simulation, we compared and analyzed the influencing factors of uncertainty in detail. In this work, the factors include: (1)  Based on the methodology, the uncertainty factor can be analyzed from two aspects: RBF implicit modeling and h-AMR-based reconstruction.

Case 1: Ellipsoidal Surface
Case 1 simulates an ellipsoidal surface. It mainly analyzes the uncertainty of influencing factors of the RBF method. Equations (16) and (17) represents the parametric equation and standard equation of the ellipsoidal surface, respectively.
(x 0 , y 0 , z 0 ) is the center of the ellipsoid. In Case 1, the ellipsoid is located at the center of the computational domain. a, b, c are the lengths of the principal axes. Table 3 lists the parameters of the ellipsoid and the h-AMR grids. To analyze the uncertainty of sample points, four sets of uniformly distributed sample points are selected to reconstruct the ellipsoidal surface. The sample points are obtained by the interpolation in Equation (16), and their numbers are approximately increased twofold. Table 4 compare the reconstructed ellipsoidal surface and the intersection points by the MQ-RBF of c = 0 and different numbers of sample points (indicated by P S ), respectively. In Table 4, P I is the number of intersections. The errors conclude the root mean square error (RMSE), mean absolute error (meanE), and maximum absolute error (maxE). T M and T R denote the time by RBF implicit modeling and the intersection calculation, respectively. T AVG indicates the reconstructed time T R consumed per intersection point. In Case 1 and Case 2, the Equation (7) is calculated by the Gaussian elimination.   In Figure 12, we can see that all sets of sample points can completely reconstruct the surface. Among the four reconstructed ellipsoidal surfaces, the deformation of (a) Ps = 26 is most obvious. There is no visible difference between the results by (c) Ps = 114 and (d) Ps = 182, indicating that the results may have converged. In Table 4, even a small number of sample points can give an acceptable result: the RMSE of Ps = 26 is smaller than 5 × 10 −2 . The errors gradually decreased with the increasing number of sample points. While the number of reconstructed intersection points P I slightly increased. Moreover, by taking Figure 13 into account, this indicates that with more sample points, the reconstructed ellipsoidal surface is fuller and more convex, which leads to an increase in the number of intersections based on the same cut cell size.

Figures 12 and 13 and
On the other hand, in Table 4, T M , T R , and T AVG all increase with the number of sample points. T M is mainly cost by the construction of the coefficient matrix and the solution of the weight coefficients of the RBF equation. T R is mainly related to the number of intersections and the RBF equation. T AVG is only affected by the RBF equation. Since the increasing Ps increases the size of the RBF coefficient matrix and RBF equation, all the time increases. Based on the results, we select the P S = 62 for further discussions.

Validation of d
To verify the robustness of d = 1 in the RBF equation, here we analyze and compare the effect of different distances d, as shown in Table 5. The RMSE% compares the percentage of the RMSE solved by different d with d = 1. It can be seen that the errors by different d are almost the same, in which the variation of the RMSE% is smaller than 1%. This indicates that the value distance d is generally independent of the accuracy of the RBF method. Additionally, there is no considerable difference between the different T AVG and T R . When d gradually reduces, the P I slightly increases. Otherwise, the P I tends to be stable. Meantime, the reconstructed surface by d = 1 requires the lowest T M . Therefore, it is feasible to use d = 1 as a general given condition for the 3D geometric modeling and reconstruction by the MQ-RBF.

Uncertainty of RBF and Parameters
Tables 6 and 7 compare the intersections reconstructed by the MQ-RBF with different shape parameters c, and the CSRBF with different continuities and support coefficients σ s , respectively. In Table 7, the minimum coefficient σ s = 0.75 is recommended by [22]. The RMSE% compares the percentage of the RMSE of different c with the c = 0.
From the results, we can see that both the MQ-RBF and CSRBF can successfully solve the coefficient matrix equation without singularity solutions. Since the ellipsoid is a typical smooth surface. The smaller error can indicate that the applied RBF has smooth characteristics. In Table 6, the errors by the MQ-RBF with non-zero shape parameters are smaller than those of c = 0. It is consistent with the idea mentioned in [37] that c works to increase the smoothness and decrease the fidelity. Specifically, the fixed Franke's empirical c and the varied Random c lead to high reconstruction accuracy, reaching 10.87% and 1.35% of the RMSE of MQ-RBF of c = 0, respectively. As to the consumption of time, the T M , and T R both slightly increase with the non-zero c. That's because the non-zero c requires extra calculations to construct the RBF coefficient matrix. Since the Ps in Case 1 is very small, the influence of the RBF coefficient matrix on the consumed time is limited. The different T AVG are almost the same. It can be concluded that the use of non-zero c can relatively reduce the errors of the reconstructed ellipsoidal surface while increasing the T M .
In Table 7, the results by the C 0 − CSRBF are very close to the MQ-RBF of c = 0. The errors of the CSRBF decrease rapidly with higher continuities (from C 2 to C 6 ), and they are even smaller than the errors of Ps = 114 and 182 in Table 4. The RMSE% of CSRBFs with continuities of C 2~C6 and σ s ≥ 1.0 are lower than 10%. It is consistent that the increase in continuities can increase the smoothness of the surface. Additionally, the errors also decrease with the increase of the support size coefficient σ s , and then gradually converge. The magnitude of the difference between the RMSE by different σ s is 10 −4 . However, in the case of using the same matrix solving method, the CSRBFs cost more T M, T AVG , and T R than the MQ-RBF of c = 0. Moreover, the T AVG slightly increases with the higher order of continuities and the larger σ s when σ s ≥ 0.75. That is because the CSRBF equation becomes complex with the higher continuities, thus increasing the matrix constructing time.
The results of Case 1 can verify the accuracy and effectiveness of our method for 3D surface reconstruction. In Case 1, for the MQ-RBF, Frank's empirical c and the Random c provide better results. Besides, the CSRBF with continuities greater than C 0 can greatly smoothen the reconstructed surface but also increase the time. The increasing support size coefficient σ s also contributes to the smoothly reconstructed surface. Taking into account the accuracy and time consumed, we only analyze the C 2 -CSRBF with σ s ≥ 1.0 in further simulations.

Case 2: Unclosed Conic Surface
Normally, only the offsets of hull lines are given before hull modeling and simulation. The deck can be closed by a plane according to the draft T and depth D. Therefore, the conic surface is selected to verify the accuracy of the proposed method on the unclosed objects. It has features of a smooth surface and sharp vertex. The parametric equations and standard equations for conic surfaces are shown in Equations (18) and (19), respectively.
(x 0 , y 0 , z 0 ) is the vertex coordinates of the cone at the center of the computational domain. a, b, c are the length of the principal axes. Table 8 presents the parameters of the cone and the h-AMR grids. In practical simulations we usually use a higher density of sample points in local areas, such as the fore and aft of the hull, to improve the global accuracy and fidelity of the surface. Thus here we use the unevenly distributed sample points of Ps = 401. In this case, the results mainly discuss the effect of RBFs and the influencing factors of the block-based h-AMR grids.

Uncertainty of RBF and Parameters
First, we analyze the effects of RBF functions and parameters on the accuracy and efficiency of the 3D unclosed surface. Based on the results of Case 1, Table 9 compares the accuracy and efficiency of the reconstructed surface by the MQ-RBF with different shape parameters c and C 2 -CSRBF with different support coefficients σ s . In Table 9, the surface reconstructed by MQ-RBF with Franke's c obtained a relatively higher accuracy. However, the MQ-RBF with Random c leads to significant errors, which is contradicting the results of Case 1. That indicates the MQ-RBF with Random c is not robust for the reconstruction of the 3D unclosed surface. Additionally, the smoothness of C 2 -CSRBF is higher than that of MQ-RBF under different shape parameters. In the case of σ s = 1.0, the RMSE% already reaches 57.3% of the MQ-RBF of c = 0. However, since the cone has a sharp top, the smooth shape parameters of RBFs have a smaller effect on the reconstructed surface than Case 1.
On the other hand, MQ-RBF of c = 0 consumes the least time. The increment of the T M , T AVG, and T R by the MQ-RBF with non-zero c is smaller than the C 2 -CSRBF. The MQ-RBF of Random c is abnormally small, indicating that it has failed to reconstruct the surface correctly and leads to an unacceptably large error. Additionally, the time consumed by the C 2 -CSRBFs increases by the augment of the σ s . Therefore, by taking consideration of the consumed time, as well as the smaller σ s being beneficial for dispersing the coefficient matrix, we utilize the CSRBF of σ s = 1.0 in further discussions.

Effect of Grid Size
This section tests the effect of grid cell size. Figure 14 compares the geometric blockbased h-AMR grids and intersections (cyan) with sample points (red) of the cone reconstructed by different cut cell sizes. Table 10 lists the reconstructed errors and time by C 2 -CSRBF of σ s = 1.0.  The results show that the density of grid size has a small effect on the reconstructed errors. The difference of the RMSE by different grid sizes is smaller than 5 × 10 −4 . Therefore, this proves the independence between the 3D reconstructed accuracy and the background grid size. Meantime, the generated intersections P I and the reconstructed T R increase with the increasing number of grids. The more P I can help to capture the more detailed features of the reconstructed surface.

Effectiveness of Geometric h-AMR Grids
To verify the effect of the block-based h-AMR on the accuracy and efficiency of 3D geometric reconstruction, we first analyze the effect of initial refinement level l min and AMR level l AMR on the reconstruction efficiency under the same cut cell size. The total refinement level is 6. Table 11 and Figure 15 compare the accuracy and efficiency of the reconstructed surface by C 2 -CSRBF of σ s = 1.0 and different l AMR . In this simulation, each grid block contains the same number of grid cells (n cx × n cy × n cz = 4 × 4 × 4).  From the results, we can see that when the cut cell size is the same and each grid block contains the same number of grid cells, the errors by different AMR levels are almost the same. This can verify the robustness of the proposed method. Additionally, the number of grid blocks, the reconstruction time T R , the initial grid generation time T IG, and the total consumed time T total are generally reduced with the higher AMR levels. When l AMR = 4, the time slightly increases. It proves the effectiveness of the h-AMR procedure in improving the efficiency of 3D reconstruction problems. The layouts of grids in Figure 15d and (e) are identical. It is due to the limitation of the block-based h-AMR grid that the gap of the level of adjacent blocks is not bigger than 1. Thus, the different AMR levels may result in an identical final h-AMR background grid.
Furthermore, we compare the accuracy and the efficiency of the cone by different numbers of grid cells contained in each grid block, as illustrated in Table 12 and Figure 16. Except for the uniform grid cell in initial level 1, the other h-AMR grid cells use the initial level 3. The cut cell size remains the same.  It can be seen from Table 12 that the errors are almost the same. The total refinement levels and block numbers correspondingly increase by the reducing number of grid cells. Additionally, T R decreases by the reduced number of grid cells. When the number of grid cells n cx × n cy × n cz decreases from 4 × 4 × 4 to 2 × 2 × 2, the time slightly increases. This proves the block-based h-AMR process can improve the reconstructed efficiency compared with the uniform Cartesian grid. Meantime, it also indicates that the increasing number of grid cells in each grid block can generally increase the flexibility and the reconstructed efficiency of the h-AMR procedure but requires more computer storage resources. Therefore, the following simulation uses the n cx × n cy × n cz = 4 × 4 × 4.

Case 3: Wigley Hull Form
In case 3, we simulate the Wigley hull to verify the effectiveness and efficiency of the proposed method for the practical ship reconstruction. As a general-purpose mathematical ship type, the Wigley hull has a simple hull form, and the body can be represented by a cluster of parabolas. It has a sharp bottom, aft, and fore, in which the fore and aft profiles are symmetrical. Equation (20) shows the exact mathematical expression of the Wigley hull.
where L W L is the length of the waterline, B is the beam, and T is the draught. They refer to the x-, y-, and zdirections, respectively. In Case 3, B/L W L = 0.1 and B/T = 1.6. The main dimensions of the Wigley hull and the parameters of h-AMR background grids are represented in Table 13. Table 13. Main dimensions of Wigley hull and parameters of h-AMR background grids.

Name Value
Length of waterline (L W L ) 60 Maximum beam of waterline (B) 6 Draft (T) 3.75 Domain length (L d ) 160 Number of grid cells (n cx = n cy = n cz ) 4 Initial grid level (l min ) 5 AMR level (l AMR ) 4 The given hull offsets are usually unclosed. In this work, we directly reconstruct the hull body by the unclosed sample points and the Boolean operation rather than adding additional deck points.

Geometric Modeling
According to the workflow in Figure 1, the sample points need to be determined first. Figure 17 shows the RBF isosurface of the Wigley hull by the MQ-RBF of c = 0 and two sets of sample points. Since the initial offset points of ships are usually not enough for accurate implicit modeling, the two sets of sample points are both obtained by interpolation from the original offsets.
In Figure 17, we can see that both sets of sample points can result in the correct reconstruction of the Wigley hull. However, due to the relatively sparse given points, the surface reconstructed by (a) Ps = 2808 sample points is comparatively not smooth at the midship section. Its waterline profile shows pronounced depression. By using the (b) Ps = 5422 sample points, the distortion of the reconstructed surface can be clearly reduced.

Surface Reconstruction
The following procedure is to generate the geometric block-based h-AMR grid and to calculate the intersection points. The computational domain (root grid block) covering the Wigley hull is generated based on the hull offsets, and then it is discretized into the initial uniform grid blocks. Based on the RBF isosurface generated by the two sets of sample points, we can obtain two layouts of geometric h-AMR grids, as illustrated in Figure 18. According to the workflow in Figure 1, the sample points need to be determined first. Figure 17 shows the RBF isosurface of the Wigley hull by the MQ-RBF of c = 0 and two sets of sample points. Since the initial offset points of ships are usually not enough for accurate implicit modeling, the two sets of sample points are both obtained by interpolation from the original offsets. In Figure 17, we can see that both sets of sample points can result in the correct reconstruction of the Wigley hull. However, due to the relatively sparse given points, the surface reconstructed by (a) Ps = 2808 sample points is comparatively not smooth at the midship section. Its waterline profile shows pronounced depression. By using the (b) Ps = 5422 sample points, the distortion of the reconstructed surface can be clearly reduced.   Based on the results of the previous cases, we select three RBFs: the MQ-RBF of c = 0, the Franke's empirical c, and the C 2 -CSRBF of σ s = 1.0 to compare the errors and time of the intersection points, as illustrated in Table 15. Figures 19 and 20 compares the RBF isosurfaces and intersections within the bounding box of the Wigley hull reconstructed by the three RBFs. From the results, we can see that the MQ-RBF of c = 0 completely reconstructs the Wigley hull with acceptable errors. It also costs the least T M and T R . Franke's empirical c produces remarkable errors compared with the MQ-RBF of c = 0. That is because the formula of Franke's empirical c is related to the distribution and the diameter of the enclosing sphere of the sample points. Since the hull has the unclosed sample points and a large aspect ratio, the effectiveness of Franke's empirical c has been affected.

Surface Reconstruction
On the other hand, the C 2 -CSRBF of σ s = 1.0 has the smallest RMSE but the reconstructed surface has some tiny leaks at the aft and bow, as shown in Figure 19c-2. This indicates that most of the 3D Wigley hull has the characteristics of convexity and smoothness. However, due to the sharp aft and bow, the reconstruction contains local inaccuracies by using the smoother RBF function. Meantime, the C 2 -CSRBF costs the most T M and T R . Overall, the MQ-RBF with c = 0 provides not only relatively highly accurate results but also stability and robustness.

Case 4: Multi-Body Hull Forms
In this section, we further validate the accuracy and efficiency of our method on 3D multi-body hull problems. The system includes the Wigley and KCS hulls. The KCS hull is the modern container ship with a bulbous bow and transom stern. The curvature of the bow and stern changes greatly compared with the hull body. It has been widely used in CFD verification. To commonly reconstruct the Wigley hull and the KCS hull, the computational domain should cover both sets of sample points. Table 16 lists the main dimensions of the KCS hull and the parameters of the h-AMR grids. First, we analyze the effect of the RBF function on the single KCS hulls, as shown in Figure 21. Based on the results of Case 3, Figure 21 compares the KCS standard explicit model and implicit RBF models reconstructed by MQ-RBF of c = 0 and C 2 -CSRBF of σ s = 1.0, respectively. It can be seen that there is no visible difference between the standard ship explicit model (a) and the implicit model by the MQ-RBF of c = 0 (b). However, our method only utilizes the scatter points for automatic modeling, which reduces the manual operation in explicit modeling. Besides, the implicit surface by the C 2 -CSRBF is smoother, except for the stern and bow which have local incorrect redundant fragments. That is because the C 2 -CSRBF possesses relatively high smoothness and convexity but the curvatures of the bow and stern of the KCS hull contain both concave and convex features. By comparison, the MQ-RBF of c = 0 provides higher fidelity. Therefore, we select the MQ-RBF of c = 0 to simulate the multi-hull system. Figure 22 represents the RBF isosurface and the sample points of the multi-body hull system. We can see that the multi-body system consists of multi-scale hull forms, and each hull is consistent with the corresponding single hull model. Due to the flexibility of the RBF modeling method, every single body can be separately reconstructed by the individual RBF function. Based on the previous results, we use MQ-RBF of c = 0 for both hulls in the multi-body system.

Multi-body Hulls Reconstruction
For the multi-body system, every single object generates its bounding box from a different set of sample points. Subsequently, the background grids within each bounding box can be locally refined through the separate RBF model and AMR levels. Figure 23 illustrates the geometric h-AMR grids of the multi-body hulls. It can be seen that each grid set follows the corresponding hull shape. Figure 24 represents the intersection points of the multi-body hulls. Specifically, we compare the errors and time between the multi-body hulls and the previous single hull, as shown in Table 17. From Table 17, it can be seen that the number of intersections P I of the comparative group is identical, and the reconstruction time T R and average reconstruction time T AVG are quite similar. The errors of the two Wigley hulls are very similar, in which the difference of the RMSE is at 10 −5 order magnitude. It also demonstrates the independence of the h-AMR background grids from the RBF model in surface reconstruction. Additionally, Figure 25 compares the typical cross-sections of the KCS hull, and the intersections are consistent with the sample points. Therefore, Case 4 verifies the accuracy and robustness of the proposed method to solve the multi-body reconstruction problem.

Discussion and Conclusions
This work designs an RBF-based h-adaptive refinement Cartesian grid method for modeling and reconstructing the complex single/multi-body geometries in ship and ocean engineering. The RBF technique is adopted within the non-body-fitted Cartesian background grid framework, in which the block-based h-AMR approach is utilized instead of the uniform Cartesian grid. Thus the accuracy of the geometric boundary can be improved, and the computational overhead can also be avoided. By combining the h-AMR grids, the RBF has achieved three functions, including the geometric modeling, geometric h-AMR grid indicator, and intersection calculation.
First, by giving the specific sets of scatter points, such as the unclosed ship offsets, the arbitrary single/multi-body RBF isosurface can be directly reconstructed without the need for additional manual operations. The usage of non-zero distance d = 1 keeps the number of samples in the RBF isosurface equation the same. Subsequently, within the automatically generated bounding box, the cut blocks and its cut cells can be searched and refined according to RBF equation value. Afterward, the Newton-Raphson iterative method is implemented to calculate the potential intersections between the h-AMR grids and the RBF isosurface. The effectiveness and robustness of the proposed method are validated by simulating the 3D analytical geometries and the practical single/multi-body hulls. The robustness of the non-zero distance d = 1 is also proved. Additionally, based on the characteristics of marine structures, the influencing factors of the proposed method are discussed in detail, including the number of sample points, RBF basis functions and parameters, grid size, and geometric h-AMR process.
The results are analyzed by the view of the RBF technique and the h-AMR background grids, respectively. For the RBF modeling method, the number of sample points P S is critical to the completeness and correctness of the reconstructed surface. Meantime, due to the characteristics of the RBF equation, the increasing number of sample points will also obviously increase the modeling time T M, and slightly cost more reconstruction time T R and average reconstruction time T AVG . That is due to the construction of the RBF equation and the matrix solving procedure. Additionally, under a certain number of sample points, using the appropriate parameters and RBFs can further improve the reconstruction accuracy. The main influencing factors include the shape parameters c for the MQ-RBF, the parameter continuities, and support coefficient σ s for the CSRBF. The non-zero shape parameters c provides the higher smoothness than c = 0. They will also cost more time due to the more complex coefficient matrix of the RBF equation, where the effect on the T R is smaller than the T M . Among non-zero c, Franke's empirical c is stable for reconstructing the smooth surface with uniformly distributed sample points. For the CSRBFs, the results of the C 0 -CSRBF are close to the MQ-RBF of c = 0. The CSRBFs with higher continuities and the σ s ≥ 1.0 can remarkably improve the accuracy of smooth surfaces. However, as for the objects containing sharp features, this smoothness and convexity will blur the sharp features and reduce the fidelity. Besides, within the same matrix solving methods, the CSRBFs will consume more T M , T R and T AVG compared with the MQ-RBFs. It can be concluded that the most concise, robust, and automatic way to model ship structures is to utilize the MQ-RBF with c = 0. It has high fidelity and is ideal for surface containing sharp features, or objects having both concave and convex curvatures. It also takes the least time.
On the other hand, the grid has a small effect on reconstruction accuracy. This proves the independence of the RBF technique from the background grid. Specifically, the finer grid size can describe more details of surface features. It also costs more T R due to the increase in intersection calculation. Besides, the effectiveness of the geometric block-based h-AMR grids is also validated. The h-AMR technique can maintain accuracy while reducing the time consumed compared with the uniform Cartesian grids. Meanwhile, within the same cut cell size, the T R can be reduced by using the higher x and the smaller number of grid cells in each grid block in a certain range.
Overall, the proposed RBF-based h-adaptive refinement Cartesian grid method is accurate, automatic, robust, and efficient for arbitrary single/multi-body geometric modeling and reconstruction. It only needs the given set of scatter points to simulate the surface, and the results are reproducible. In addition, we found the MQ-RBF of c = 0 is more robust and suitable to handle the complex single/multi-body hull structures containing sharp and smooth features. The method can be easily employed in other marine structures and other engineering fields. However, the traditional way to solve RBF is computationally expensive. In further research, we will focus on the local RBF technique. It can help reduce the geometric modeling time T M . Meantime, it is also possible to use the specific RBF in different areas, allowing for an improvement in local smoothness/convexity or high fidelity depending on the surface features.