Fix floating precision in circumsphere_of_triangles

tests/test_mesh.py

@@ -70,6 +70,33 @@ def test_move_mesh(self):
         for r in mesh.eradius:
             self.assertLessEqual(r, max_cell_circumradius*1.001)
 
+    def test_compute_mesh_spheres(self):
+
+        coord = np.array([
+                [-1.3856244e+02, 6.9698529e+02, -2.9481543e+02],
+                [5.1198740e+02, 5.7321143e+02, -1.3369661e+01],
+                [-2.4491163e-01, -2.4171574e-01, 1.4720881e-01],
+                [2.7638666e+02, 6.1436609e+02, -3.7048819e+02]
+            ])
+
+        enod = np.array([[0, 1, 2, 3]])
+        mesh = TetraMesh.generate_mesh_from_vertices(coord, enod)
+
+
+        theta = np.deg2rad(1e-4)
+        c, s = np.cos(theta), np.sin(theta)
+        R_z = np.array([
+            [c, -s, 0],
+            [s, c, 0],
+            [0, 0, 1]
+        ])
+
+        rotated_coord = R_z.dot(mesh.coord.T).T
+
+        eradius_rotated, _ = mesh._compute_mesh_spheres(rotated_coord, mesh.enod)
+
+        self.assertAlmostEqual(mesh.eradius[0], eradius_rotated[0])
+
     def test_update(self):
         rotation_axis=np.array([0, 0, 1.0])
         rotation_angle=np.pi/4.37

xrd_simulator/utils.py

@@ -341,7 +341,7 @@ def circumsphere_of_segments(segments):
     """Computes the minimum circumsphere of n segments given by a numpy array of vertices nx2x3"""
     centers = np.mean(segments,axis=1)
     radii = np.linalg.norm(centers-segments[:,0,:],axis=1)
-    return centers, radii
+    return centers, radii*1.0001 # because loss of floating point precision
 
 def circumsphere_of_triangles(triangles):
     """Computes the minimum circumsphere of n triangles given by a numpy array of vertices nx3x3. Prints a message if any tetrahedron has 0 volume."""
@@ -358,7 +358,7 @@ def circumsphere_of_triangles(triangles):
     centers = triangles[:,0,:]+ a_to_centre    
     radii = np.linalg.norm(a_to_centre,axis=1)
 
-    return centers, radii
+    return centers, radii*1.0001 # because loss of floating point precision
 
 def circumsphere_of_tetrahedrons(tetrahedra):
     """Computes the circumcenter of n tetrahedrons given by a numpy array of vertices nx4x3"""
@@ -375,8 +375,7 @@ def circumsphere_of_tetrahedrons(tetrahedra):
     centers = np.matmul(np.linalg.inv(A),B[:,:,np.newaxis])[:,:,0]
     radii = np.linalg.norm((tetrahedra.transpose(2,0,1)-centers.transpose(1,0)[:,:,np.newaxis])[:,:,0],axis=0)
 
-    return centers, radii
-
+    return centers, radii*1.0001 # because loss of floating point precision
 
 def sizeof_fmt(num, suffix='B'):
     ''' by Fred Cirera,  https://stackoverflow.com/a/1094933/1870254, modified'''