## 3 3D SPH Surface-to-Surface Contact Algorithm Using Surface Reconstruction
### Detection of Contacts  
Through the surface reconstruction procedure, reconstructed surfaces can be obtained individually for each object. To construct a general surface-to-surface contact algorithm, all the triangular faces on the reconstructed surfaces are collected for contact detection. Each triangular face is called a master segment. Meanwhile, all the particles on the reconstructed surfaces along with isolated particles (particles not connected to any remaining tetrahedrons), which possibly appear in extremely large deformation problems, are collected. Each of these particles is called a slave particle. With such a treatment, it is easy to deal with complex contact problems without troublesome identification and definition of possible contacts between different objects. However, this comes with a relatively expensive computational cost. The task of contact detection is to identify all slave particle-master segment pairs in contact. To achieve this, the master segments are processed one by one, and, for each master segment, a strategy of a box test followed by a crossover test is used to search for slave particles in contact with the master segment. The search strategy used here is similar to that used for the finite element method [25]. For a given master segment, the box test is to find the slave particles possibly in contact with it from all the slave particles. Figure 2 shows a schematic diagram of the box test. A  
Figure 1: Surface reconstruction for an 3D SPH object: (**a**) Particles representing the object; (**b**) Tetrahedrons obtained by Delaunay triangulation of particles; (**c**) Remaining tetrahedrons (including internal faces); (**d**) Reconstructed surface of the object (no internal faces). capture box is constructed by extending the minimum bounding box of the master segment outward for a distance \(d_{\rm e}\). The extending distance \(d_{\rm e}\) is calculated as  
\[d_{\rm e}=\alpha_{\rm e}V_{\rm r}\Delta t, \tag{11}\]  
where \(\alpha_{\rm e}\) is the coefficient of extending distance, for which 1.2 is sufficient for most problems, \(\Delta t\) is the time step, and \(V_{\rm r}\) is the maximum relative velocity of all particles which is defined as  
\[V_{\rm r}=\sqrt{\left(v_{\rm max}^{x}-v_{\rm min}^{x}\right)^{2}+\left(v_{\rm max }^{y}-v_{\rm min}^{y}\right)^{2}+\left(v_{\rm max}^{z}-v_{\rm min}^{z}\right)^{ 2}}, \tag{12}\]  
where \(v_{\rm max}^{x}\) and \(v_{\rm min}^{x}\) are, respectively, the maximum and minimum velocities of all particles in \(x\)-direction, and \(v_{\rm max}^{y}\), \(v_{\rm min}^{y}\), \(v_{\rm max}^{z}\), \(v_{\rm min}^{z}\) have the same meanings for \(y\)- and \(z\)-directions. To accelerate the box test process, the efficient point-in-box algorithm [26] is employed to identify the slave particles in the capture box. Those salve particles identified by the box test are contact candidates for the master segment. The crossover test is to analyze whether the salve particle candidates found by the box test are actually in contact with the master segment. Figure 3 shows a schematic diagram of the crossover test. For a salve particle candidate, its penetration distance into the master segment and the position of its projection point to the master segment need to be determined. The penetration distance \(\delta_{\rm P}\) is calculated as  
\[\delta_{\rm P}=-\big{[}n^{1}(x_{s}-x_{1})+n^{2}(y_{s}-y_{1})+n^{3}(z_{s}-z_{1}) \big{]}, \tag{13}\]  
where \(\mathbf{n}=\big{(}n^{1},\ n^{2},\ n^{3}\big{)}=\mathbf{n}^{\prime}\big{/}\|\mathbf{n}^{ \prime}\|\) is the unit outward normal vector of the master segment, and \(\mathbf{n}^{\prime}\) is the outward normal vector which can be calculated as  
\[\mathbf{n}^{\prime}=\mathbf{x}_{21}\times\mathbf{x}_{32}=\big{(}A^{\prime},B^{\prime},C^{ \prime}\big{)}, \tag{14}\]  
where \(\mathbf{x}_{21}=\mathbf{x}_{2}-\mathbf{x}_{1}\), \(\mathbf{x}_{32}=\mathbf{x}_{3}-\mathbf{x}_{2}\), \(A^{\prime}=y_{21}z_{32}-y_{32}z_{21}\), \(B^{\prime}=x_{32}z_{21}-x_{21}z_{32}\), and \(C^{\prime}=x_{21}y_{32}-x_{32}y_{21}\).