# The Gradient and the Hessian of the Distance between Point and Triangle in 3D
## 4 Algorithm and Testing  
The calculations are implemented in double-precision floating-point arithmetic in C, and the tests are performed with squared distance \(s\) and its first and second derivatives. The first step is to determine the partition to which \(\mathbf{p}_{c}\) belongs (Figure 2). This step coincides with the original point-triangle algorithm [1], but subsequent calculations are extended to evaluate \(\frac{\partial s}{\partial x_{i}}\) and \(\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}\). If \(\mathbf{p}_{c}\) belongs to partitions 1, 3, 5, then the point-line algorithm is invoked. For partitions 2, 4, 6, point-point calculations are performed. For partition 0, the point-plane algorithm is used.  
For testing, 10 million cases were generated, about \(2/3\) of which are random coordinates that come from the uniform distribution in the range \([-1,1]\). The remaining cases come from the finite element simulation of fracture, where the point \(\mathbf{p}_{c}\) is often close to one of the triangle's vertices. Most of the cases that come from the simulation have a low ratio of the distance to the shortest edge of the triangle, in the range between \(10^{-8}\) and \(10^{-5}\). Such test cases result in a lower accuracy of the final answer than the random arrangements of points.  
Figure 3: The closest point \(\mathbf{p}_{c}\) may coincide with (**a**) one of triangle’s vertices or (**b**) belong to one of the edges.  
Figure 2: Partitioning of the plane by the triangle domain. Different domains are distinguished by the values of barycentric coordinates of the projection of \(\mathbf{p}_{0}\) onto the plane of the triangle.  
To evaluate the accuracy of the proposed algorithm, calculations are first performed using arbitrary precision arithmetic with at least 17 digits calculated precisely. These results are denoted \(s_{P}\frac{\partial s_{P}}{\partial x_{i}}\), \(\frac{\partial^{2}s_{P}}{\partial x_{i}\partial x_{j}}\) and are compared to the results obtained in floating-point arithmetic \(s_{c}\), \(\frac{\partial s_{c}}{\partial x_{i}}\), \(\frac{\partial^{2}s_{c}}{\partial x_{i}\partial x_{j}}\). Relative errors are computed separately for the squared distance, its first derivatives and its second derivatives:  
\[E_{0}=\left|(s_{c}-s_{p})/s_{p}\right|, \tag{21}\]  
\[E_{1}=\max_{i}\left|\left(\frac{\partial s_{c}}{\partial x_{i}}-\frac{ \partial s_{p}}{\partial x_{i}}\right)/\frac{\partial s_{p}}{\partial x_{i}} \right|, \tag{22}\]  
\[E_{2}=\max_{i,j}\left|\left(\frac{\partial^{2}s_{c}}{\partial x_{i}\partial x _{j}}-\frac{\partial^{2}s_{p}}{\partial x_{i}\partial x_{j}}\right)/\frac{ \partial^{2}s_{p}}{\partial x_{i}\partial x_{j}}\right|. \tag{23}\]  
\(E_{0}\) is the relative error of the squared distance and is used as a baseline to compare with the relative errors of the first and second derivatives \(E_{1}\) and \(E_{2}\). Ideally, the values of \(E_{0}\), \(E_{1}\) and \(E_{2}\) would have the same order of magnitude. However, due to the large number of algebraic operations, the accuracy of the calculation of the second derivative is lower. The values \(E_{0}\), \(E_{1}\) and \(E_{2}\) cover a wide range of values. In about half of the cases, the precision for calculating second derivatives is better than \(10^{-13}\). However, the practical interest lies in investigating the worst cases because one incorrect calculation could affect a whole scientific study.