From d1453260590c1942e4bfc2dc4e2e7f58d88c3f60 Mon Sep 17 00:00:00 2001
From: Axel Henningsson <nilsaxelhenningsson@gmail.com>
Date: Sat, 1 Jun 2024 23:32:51 +0200
Subject: [PATCH] Fix floating precision in circumsphere_of_triangles

---
 tests/test_mesh.py     | 27 +++++++++++++++++++++++++++
 xrd_simulator/utils.py |  7 +++----
 2 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 72c3703..0b40cb8 100644
--- a/tests/test_mesh.py
+++ b/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
diff --git a/xrd_simulator/utils.py b/xrd_simulator/utils.py
index 44e5f02..e7b5b3b 100644
--- a/xrd_simulator/utils.py
+++ b/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'''