fast voxel centers and interpolating stress data to grid data

utils/planner.py

@@ -5,6 +5,7 @@ from pyvista import examples
 from math import pi
 import trimesh as tm
 import pymeshfix as pf
+from scipy.spatial import KDTree
 
 
 class Planner:
@@ -43,14 +44,24 @@ class Planner:
         print(c-s)
         print("Getting Voxel centers")
         points_to_add = []
-        for i in range(self.grid.n_points):
-            # TODO This logic is not necessarily accurate - voxel could still be empty
-            # but should be fine for our case
-            if np.any(
-                    np.all(self.grid.points[i] + [voxel_dim, voxel_dim, layer_height] == self.grid.points,
-                        axis=1)):
-                points_to_add.append(i)
-
-        self.centers = self.grid.points[points_to_add] + [voxel_dim/2, voxel_dim/2, layer_height/2]
-        print("Done")
-        print(time.time()-c)
+        temptree = KDTree(self.grid.points, leafsize=5) #kd tree for fast lookup of adjacent grid points
+        # TODO This logic is not necessarily accurate - voxel could still be empty
+        # but should be fine for our case
+        d1, ii  = temptree.query(self.grid.points + [voxel_dim, voxel_dim, layer_height])
+        d2, ii  = temptree.query(self.grid.points + [voxel_dim, 0, layer_height])
+        d3, ii  = temptree.query(self.grid.points + [0, voxel_dim, 0])
+        self.centers = self.grid.points[d1 + d2 + d3 < 1e-4] + [voxel_dim/2, voxel_dim/2, layer_height/2]
+
+        gs = time.time()
+        print(gs-c)
+        print("Interpolating pointcloud to grid stresses")
+        # get mean stress at each voxel
+
+        self.centtree = KDTree(self.centers, leafsize=5) #kd tree for mapping grid points to fea mesh
+        print(time.time()- gs)
+        self.gridstresses = np.zeros(self.centers.shape[0])
+        dd, ii = self.centtree.query(self.centers, k=4)
+        self.gridstress = np.sum(1/dd * self.mesh['StressValues'][ii], axis=1)/np.sum(1/dd, axis=1)
+#        print(time.time()- gs)
+
+