Fixes #119 (#120)

* avoiding NaN in doublearea

* added zero area check in RFTS

src/gpytoolbox/doublearea_intrinsic.py

@@ -39,6 +39,9 @@ def doublearea_intrinsic(l_sq,F):
     # Using Kahan's formula
     # https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
     a,b,c = l[:,0], l[:,1], l[:,2]
-    dblA = 0.5 * np.sqrt((a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c)))
-
+    # previously (gave NaNs for very small triangles)
+    # dblA = 0.5 * np.sqrt((a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c)))
+    arg = (a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c))
+    dblA = 0.5 * np.sqrt(np.maximum(arg, 0.))
+
     return dblA

src/gpytoolbox/reach_for_the_spheres.py

@@ -663,7 +663,12 @@ def reach_for_the_spheres_iteration(state,
         state.V = sp.sparse.linalg.spsolve(Q,b)
 
     # catching flow singularities so we fail gracefully
-    if np.any((np.isnan(state.V))):
+    if np.any((np.isnan(state.V))) or np.any(doublearea(state.V, state.F)==0): 
+        print(A)
+        print(R)
+        print(b)
+        print(state.V)
+
         if verbose:
             print("we found a flow singularity. Returning the last converged solution.")
         state.V = state.V_last_converged.copy()

test/test_doublearea_intrinsic.py

@@ -23,6 +23,18 @@ def test_many_equilateral_triangles(self):
         A = gpy.doublearea_intrinsic(l_sq,f)
         self.assertTrue(np.isclose(A, 2.*np.full(f.shape[0],np.sqrt(3.)/4. * c)).all())
 
+    def test_very_tiny_triangle(self):
+        c = np.random.default_rng().random() + 0.1
+
+        l_sq = c * np.array([[0.0000001,0.0000001,1.]])
+        f = np.array([[0,1,2]])
+
+
+        # this returns NaN as of v0.2.0
+        A = gpy.doublearea_intrinsic(l_sq,f)
+        # it shouldn't anymore
+        self.assertTrue(np.isnan(A[0]) == False)
+
 
 
 

test/test_reach_for_the_spheres.py

@@ -72,5 +72,20 @@ def test_segfault(self):
         # this should not segfault
         gpy.write_mesh("test_last_converged.obj", Vr, Fr)
 
+    def test_singularity(self):
+        V,F = gpy.read_mesh("test/unit_tests_data/horse.obj")
+        # is mesh normalized? print corners
+        # print(np.min(V, axis=0))
+        # print(np.max(V, axis=0))
+        V = gpy.normalize_points(V)
+        j = 32
+        sdf = lambda x: gpy.signed_distance(x, V, F)[0]
+        gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
+        U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T
+        V0, F0 = gpy.icosphere(2)
+        # this should not crash, we should catch the singularity and output a the last converged mesh
+        Vr,Fr = gpy.reach_for_the_spheres(U, sdf, V0, F0, min_h = .01, verbose = False)
+
+
 if __name__ == '__main__':
     unittest.main()