diff options
Diffstat (limited to 'geometry.py')
| -rw-r--r-- | geometry.py | 75 |
1 files changed, 25 insertions, 50 deletions
diff --git a/geometry.py b/geometry.py index a98b5cd..0657a77 100644 --- a/geometry.py +++ b/geometry.py @@ -1,8 +1,5 @@ import sys import numpy as np -import pycuda.driver as cuda -from pycuda import gpuarray -from itertools import * from itertoolset import * from tools import timeit @@ -159,7 +156,7 @@ def interleave(arr, bits): if len(arr.shape) != 2 or arr.shape[1] != 3: raise Exception('shape mismatch') - z = np.zeros(arr.shape[0], dtype=np.uint32) + z = np.zeros(arr.shape[0], dtype=np.uint64) for i in range(bits): z |= (arr[:,2] & 1 << i) << (2*i) | \ (arr[:,1] & 1 << i) << (2*i+1) | \ @@ -174,13 +171,13 @@ def morton_order(mesh, bits): """ lower_bound, upper_bound = mesh.get_bounds() - if bits <= 0 or bits > 10: - raise Exception('number of bits must be in the range (0,10].') + if bits <= 0 or bits > 21: + raise Exception('number of bits must be in the range (0,21].') max_value = 2**bits - 1 def quantize(x): - return np.uint32((x-lower_bound)*max_value/(upper_bound-lower_bound)) + return np.uint64((x-lower_bound)*max_value/(upper_bound-lower_bound)) mean_positions = quantize(np.mean(mesh.assemble(), axis=1)) @@ -212,7 +209,6 @@ class Geometry(object): return len(self.solids)-1 @timeit - #@profile def build(self, bits=8, shift=3): offsets = [ (0,0) ] for solid in self.solids: @@ -275,61 +271,40 @@ class Geometry(object): self.upper_bounds = np.empty((unique_zvalues.size,3), dtype=np.float32) assembled_mesh = self.mesh.assemble(group=False) - zi1_array = np.searchsorted(zvalues_mesh, unique_zvalues) - zi2_array = np.searchsorted(zvalues_mesh, unique_zvalues, side='right') + self.node_map = np.searchsorted(zvalues_mesh, unique_zvalues) + self.node_map_end = np.searchsorted(zvalues_mesh, unique_zvalues, side='right') - for i, (zi1, zi2) in enumerate(zip(zi1_array,zi2_array)): - self.lower_bounds[i] = assembled_mesh[zi1*3:zi2*3].min(0) - self.upper_bounds[i] = assembled_mesh[zi1*3:zi2*3].max(0) + for i, (zi1, zi2) in enumerate(izip(self.node_map, self.node_map_end)): + self.lower_bounds[i] = assembled_mesh[zi1*3:zi2*3].min(axis=0) + self.upper_bounds[i] = assembled_mesh[zi1*3:zi2*3].max(axis=0) - zvalues = np.array(unique_zvalues, copy=True) - self.node_map = np.array(zi1_array, copy=True) - self.node_length = zi2_array - zi1_array self.layers = np.zeros(unique_zvalues.size, dtype=np.uint32) self.first_node = unique_zvalues.size begin_last_layer = 0 - layer = 1 + for layer in count(1): + bit_shifted_zvalues = unique_zvalues >> shift + unique_zvalues = np.unique(bit_shifted_zvalues) - while True: - bit_shifted_zvalues = zvalues >> shift - unique_bit_shifted_zvalues = np.unique(bit_shifted_zvalues) - zvalues = np.empty(unique_bit_shifted_zvalues.size, dtype=np.uint64) + i0 = begin_last_layer + bit_shifted_zvalues.size - self.lower_bounds.resize(\ - (self.lower_bounds.shape[0]+unique_bit_shifted_zvalues.size,3)) - self.upper_bounds.resize(\ - (self.upper_bounds.shape[0]+unique_bit_shifted_zvalues.size,3)) + self.node_map.resize(self.node_map.size+unique_zvalues.size) + self.node_map[i0:] = np.searchsorted(bit_shifted_zvalues, unique_zvalues) + begin_last_layer + self.node_map_end.resize(self.node_map_end.size+unique_zvalues.size) + self.node_map_end[i0:] = np.searchsorted(bit_shifted_zvalues, unique_zvalues, side='right') + begin_last_layer - self.node_map.resize(\ - self.node_map.size+unique_bit_shifted_zvalues.size) - self.node_length.resize(\ - self.node_length.size+unique_bit_shifted_zvalues.size) + self.layers.resize(self.layers.size+unique_zvalues.size) + self.layers[i0:] = layer - self.layers.resize(self.layers.size+unique_bit_shifted_zvalues.size) + self.lower_bounds.resize((self.lower_bounds.shape[0]+unique_zvalues.size,3)) + self.upper_bounds.resize((self.upper_bounds.shape[0]+unique_zvalues.size,3)) - for i, z in enumerate(unique_bit_shifted_zvalues): - i1 = np.searchsorted(bit_shifted_zvalues, z) + \ - begin_last_layer - i2 = np.searchsorted(bit_shifted_zvalues, z, side='right') + \ - begin_last_layer - - zvalues[i] = z - - i += begin_last_layer + bit_shifted_zvalues.size - - self.lower_bounds[i] = self.lower_bounds[i1:i2].min(axis=0) - self.upper_bounds[i] = self.upper_bounds[i1:i2].max(axis=0) - - self.node_map[i] = i1 - self.node_length[i] = i2-i1 - self.layers[i] = layer + for i, zi1, zi2 in izip(count(i0), self.node_map[i0:], self.node_map_end[i0:]): + self.lower_bounds[i] = self.lower_bounds[zi1:zi2].min(axis=0) + self.upper_bounds[i] = self.upper_bounds[zi1:zi2].max(axis=0) begin_last_layer += bit_shifted_zvalues.size - layer += 1 - - if unique_bit_shifted_zvalues.size == 1: + if unique_zvalues.size == 1: break - |