From 2cbdab7feb1cafe4491766f056b3fe58e32ac706 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 21 Feb 2022 14:34:46 +1100 Subject: [PATCH] fix: :zap: changing aabb to sparse matrix previous numpy array used a lot of ram and was slow. Significant memory and speed improvements --- .../supports/_3d_unstructured_tetra.py | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py index ac05f5084..368fb65fd 100644 --- a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py +++ b/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,9 +362,11 @@ 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] @@ -368,8 +374,12 @@ def get_element_for_location(self, points): * 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, :]