fix: ⚡ changing aabb to sparse matrix

previous numpy array used a lot of ram and was slow. Significant memory and speed improvements

LoopStructural/interpolators/supports/_3d_unstructured_tetra.py

@@ -4,6 +4,7 @@
 import logging
 
 import numpy as np
+from scipy.sparse import csr_matrix
 
 from ._3d_base_structured import BaseStructuredSupport
 from LoopStructural.utils import getLogger
@@ -62,10 +63,9 @@ def __init__(self, nodes, elements, neighbours, aabb_nsteps=None):
         # make a big table to store which tetra are in which element.
         # if this takes up too much memory it could be simplified by using sparse matrices or dict but
         # at the expense of speed
-        self.aabb_table = np.zeros(
+        self.aabb_table = csr_matrix(
             (self.aabb_grid.n_elements, len(self.elements)), dtype=bool
         )
-        self.aabb_table[:] = False
         self.shared_element_relationships = np.zeros((self.elements.shape[0]*3,2),dtype=int)
         self.shared_elements = np.zeros((self.elements.shape[0]*3,3),dtype=int)
         self._init_face_table()
@@ -168,7 +168,7 @@ def _initialise_aabb(self):
         logic = np.logical_and(x_logic, y_logic)
         logic = np.logical_and(logic, z_logic)
 
-        self.aabb_table = logic
+        self.aabb_table = csr_matrix(logic)
 
     @property
     def ntetra(self):
@@ -331,6 +331,10 @@ def evaluate_gradient(self, pos, property_array):
         return values
 
     def inside(self, pos):
+        if pos.shape[1] > 3:
+            logger.warning(f'Converting {pos.shape[1]} to 3d using first 3 columns')
+            pos = pos[:, :3]
+
         inside = np.ones(pos.shape[0]).astype(bool)
         for i in range(3):
             inside *= pos[:, i] > self.origin[None, i]
@@ -358,18 +362,24 @@ def get_element_for_location(self, points):
         -------
 
         """
+
         cell_index = np.array(self.aabb_grid.position_to_cell_index(points)).swapaxes(
             0, 1
         )
+        inside = self.aabb_grid.inside(points)
         global_index = (
             cell_index[:, 0]
             + self.aabb_grid.nsteps_cells[None, 0] * cell_index[:, 1]
             + self.aabb_grid.nsteps_cells[None, 0]
             * self.aabb_grid.nsteps_cells[None, 1]
             * cell_index[:, 2]
         )
-        tetra_indices = self.aabb_table[global_index, :]
-        row, col = np.where(tetra_indices)
+
+        tetra_indices = self.aabb_table[global_index[inside], :].tocoo()
+        # tetra_indices[:] = -1
+        row = tetra_indices.row
+        col = tetra_indices.col
+        # using returned indexes calculate barycentric coords to determine which tetra the points are in
         vertices = self.nodes[self.elements[col, :4]]
         pos = points[row, :]
         vap = pos[:, :] - vertices[:, 0, :]