## 3 3D SPH Surface-to-Surface Contact Algorithm Using Surface Reconstruction
### 3D Surface Reconstruction of SPH Objects
For some extremely distorted tetrahedrons, the quantity under the root square sign in the expression of \(S_{I}\) may be calculated as a very small minus value or zero due to the round-off error of numerical calculation. For this case, special treatment should be adopted to avoid numerical calculation difficulty and obtain reasonable results. In our implementation, the value \(S_{I}\) is set to a very large positive value, then a large \(R_{I}\) will be calculated from Equation (8) and the corresponding tetrahedron will be dropped. Similarly, the value of \(V_{I}\) calculated with Equation (9) may be minus or zero. If it is zero, then \(R_{I}\) is not calculated from Equation (8), but rather directly set to a very large positive value, so the corresponding tetrahedron will be dropped. If it is minus, then it is set to the absolute value and \(R_{I}\) is still calculated from Equation (8). In order to search the triangular faces on the boundaries of the remaining tetrahedrons and obtain the reconstructed surfaces of the object, an algorithm consisting of two stages (see Algorithm 1) is employed. The first stage is to determine the remaining tetrahedrons connected to each particle, and the second stage is to search the boundary faces of the remaining tetrahedrons using the information obtained in the first stage. The detailed search algorithm is given as follows:  
```
Stage 1 Determine remaining tetrahedrons connected to each particle 1.1. Initialize total number of remaining tetrahedrons connected each particle: nTetPar = 0.1.2. While \(i\leq n_{\rm T}\)do
1.3. While \(j\leq 4\)do
1.4. Get number of particle at the \(j\)th vertice of tetrahedron \(i\): \(k=\) Tet(\(i\), \(j\))
1.5. Update total number of tetrahedrons connected to particle \(k\): nTetPar(\(k\)) = nTetPar(\(k\)) + 1.6. Record number of tetrahedron connected to particle \(k\): TetPar(\(k\), nTetPar(\(k\))) = \(i\)
1.7. \(j=j+1\)
1.8. \(i=i+1\)

Stage 2 Search boundary faces of remaining tetrahedrons
2.1. Initialize total number of boundary faces and internal face flag: nBndFac = 0, IntFacFlag = False
2.2. While \(i\leq n_{\rm T}\)do
2.3. While \(j\leq 4\)do
2.4. If IntFacFlag(\(i\), \(j\)) is True, then the \(j\)th face of tetrahedron \(i\) is an internal face, go to step 2.17
2.5. Get numbers of the three particles on the \(j\)th face of tetrahedron \(i\): \(p_{1}\), \(p_{2}\), \(p_{3}\)
2.6. Get the total number of tetrahedrons connected to particle \(p_{1}\): \(num\) = nTetPar(\(p_{1}\))
2.7. While \(k\leq num\)do
2.8. Get the number of the \(k\)th tetrahedron connected to particle \(p_{1}\): \(l=\) TetPar(\(p_{1}\), \(k\))
2.9. If \(i=l\), then skip tetrahedron \(l\), and go to step 2.15
2.10. While \(m\leq 4\)do
2.11. Get numbers of the three particles on the \(m\)th face of tetrahedron \(l\): \(q_{1}\), \(q_{2}\), \(q_{3}\)
2.12. Compare \(p_{1}\), \(p_{2}\), \(p_{3}\) with \(q_{1}\), \(q_{2}\), \(q_{3}\). If they coincide, then go to Step 2.13, else go to 2.14
2.13. Set internal face flags: IntFacFlag(\(i\), \(j\)) = True, IntFacFlag(\(l\), \(m\)) = True; and go to step 2.17
2.14. \(m=m+1\)
2.15. \(k=k+1\)
2.16. Identify the \(j\)th face of tetrahedron \(i\) as a boundary face and record the face
2.17. \(j=j+1\)
2.18. \(i=i+1\)

\({}^{*}\)\(n_{\rm T}\) is the total number of remaining tetrahedrons; Loop variables are initialized to be 1 just before loops begin. ```  
**Algorithm 1**Searching of triangular faces on boundaries of remaining tetrahedrons  
The present surface reconstruction procedure is available for arbitrarily shaped 3D objects. It can produce 3D explicit surfaces which properly describe the shape of an object. It eliminates the need for predefining surfaces with connectivity between boundary particles in the initial phase of the SPH simulation. It enables us to develop surface-to-surfacecontact algorithms available for problems with arbitrary deformation, including material failure problems, where new surfaces may be generated during the deformation process.