triangle mesh is now stored everywhere as a split list of vertices and triangles

10 files changed, 84 insertions, 72 deletions

diff --git a/detectors/lbne.py b/detectors/lbne.py
index 81f5dd7..3c28c58 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -1,20 +1,16 @@
-import numpy as np
-from copy import deepcopy
-from itertools import product
-
 import os
 import sys
+import numpy as np
 
 dir = os.path.split(os.path.realpath(__file__))[0]
-
 sys.path.append(dir + '/..')
 
 import models
-from mesh import *
-from transform import rotate
+from mesh import mesh_from_stl
 from solid import Solid
 from geometry import Geometry
-from materials import glass, h2o
+from transform import rotate, make_rotation_matrix
+from itertools import product
 
 endcap_spacing = .485
 
@@ -28,33 +24,35 @@ class LBNE(Geometry):
     def __init__(self):
         super(LBNE, self).__init__()
 
-        pmt_mesh = read_stl(models.dir + '/hamamatsu_12inch.stl')/1000.0
-        pmt_solid = Solid(pmt_mesh, glass, h2o)
+        pmt_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch.stl')
+        pmt_mesh.vertices /= 1000.0
+
+        pmtid = 0
 
         # construct the barrel
         for i in range(pmts_per_string):
             for j in range(nstrings):
-                pmt = deepcopy(pmt_solid)
-                pmt.mesh += (-radius,0,-height/2+i*height/(pmts_per_string-1))
-                pmt.mesh = rotate(pmt.mesh, j*2*np.pi/nstrings, (0,0,1))
-                self.add_solid(pmt)
+                rotation = make_rotation_matrix(j*2*np.pi/nstrings, (0,0,1))
+                displacement = rotate((-radius,0,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1))
+                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
+                pmtid += 1
 
         # construct the top endcap
         for x, y in np.array(tuple(product(\
                     np.arange(-radius+0.075, radius, endcap_spacing),
                     np.arange(-radius+0.075, radius, endcap_spacing)))):
             if np.sqrt(x**2 + y**2) <= radius:
-                pmt = deepcopy(pmt_solid)
-                pmt.mesh = rotate(pmt.mesh, -np.pi/2, (0,1,0))
-                pmt.mesh += (x,y,+height/2+height/(pmts_per_string-1)/2)
-                self.add_solid(pmt)
+                rotation = make_rotation_matrix(-np.pi/2, (0,1,0))
+                displacement = (x,y,+height/2+height/(pmts_per_string-1)/2)
+                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
+                pmtid += 1
 
         # construct the bottom endcap
         for x, y in np.array(tuple(product(\
                     np.arange(-radius+0.075, radius, endcap_spacing),
                     np.arange(-radius+0.075, radius, endcap_spacing)))):
             if np.sqrt(x**2 + y**2) <= radius:
-                pmt = deepcopy(pmt_solid)
-                pmt.mesh = rotate(pmt.mesh, +np.pi/2, (0,1,0))
-                pmt.mesh += (x,y,-height/2-height/(pmts_per_string-1)/2)
-                self.add_solid(pmt)
+                rotation = make_rotation_matrix(+np.pi/2, (0,1,0))
+                displacement = (x,y,-height/2-height/(pmts_per_string-1)/2)
+                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
+                pmtid += 1
diff --git a/geometry.py b/geometry.py
index d1bec2a..fdc361c 100644
--- a/geometry.py
+++ b/geometry.py
@@ -2,6 +2,7 @@ import numpy as np
 import numpy.ma as ma
 import pycuda.driver as cuda
 from pycuda import gpuarray
+from mesh import Mesh
 
 def interleave(arr, bits):
     """
@@ -56,16 +57,23 @@ class Geometry(object):
         self.solids.append(solid)
 
     def build(self, bits=8):
-        self.mesh = np.concatenate([solid.build() for solid in self.solids])
+        vertices = []
+        triangles = []
+        for solid in self.solids:
+            triangles.extend(solid.mesh.triangles + len(vertices))
+            vertices.extend(np.inner(solid.mesh.vertices, solid.rotation) + \
+                                solid.displacement)
 
-        zvalues_mesh = morton_order(self.mesh, bits)
+        self.mesh = Mesh(vertices, triangles)
+
+        zvalues_mesh = morton_order(self.mesh[:], bits)
         reorder = np.argsort(zvalues_mesh)
         zvalues_mesh = zvalues_mesh[reorder]
 
         if (np.diff(zvalues_mesh) < 0).any():
             raise Exception('zvalues_mesh out of order.')
 
-        self.mesh = self.mesh[reorder]
+        self.mesh.triangles = self.mesh.triangles[reorder]
 
         material1 = np.concatenate([solid.material1 for solid in self.solids])
         material2 = np.concatenate([solid.material2 for solid in self.solids])
@@ -96,9 +104,8 @@ class Geometry(object):
         for i, surface in enumerate(surface2[reorder]):
             self.surface2_index[i] = self.surfaces.index(surface)
 
-        self.solid_index = \
-            np.concatenate([np.tile(i, len(solid.mesh)) \
-                                for i in range(len(self.solids))])[reorder]
+        self.solid_id = np.concatenate([np.tile(solid.id, len(solid.mesh)) \
+                                            for solid in self.solids])[reorder]
 
         unique_zvalues = np.unique(zvalues_mesh)
         zvalues = np.empty(unique_zvalues.size, dtype=np.uint64)
@@ -112,13 +119,13 @@ class Geometry(object):
             i1 = np.searchsorted(zvalues_mesh, z)
             i2 = np.searchsorted(zvalues_mesh, z, side='right')
 
-            self.lower_bounds[i] = [np.min(self.mesh[i1:i2,:,0]),
-                                    np.min(self.mesh[i1:i2,:,1]),
-                                    np.min(self.mesh[i1:i2,:,2])]
+            self.lower_bounds[i] = [np.min(self.mesh[i1:i2][:,:,0]),
+                                    np.min(self.mesh[i1:i2][:,:,1]),
+                                    np.min(self.mesh[i1:i2][:,:,2])]
             
-            self.upper_bounds[i] = [np.max(self.mesh[i1:i2,:,0]),
-                                    np.max(self.mesh[i1:i2,:,1]),
-                                    np.max(self.mesh[i1:i2,:,2])]
+            self.upper_bounds[i] = [np.max(self.mesh[i1:i2][:,:,0]),
+                                    np.max(self.mesh[i1:i2][:,:,1]),
+                                    np.max(self.mesh[i1:i2][:,:,2])]
             
             self.node_map[i] = i1
             self.node_length[i] = i2-i1
@@ -178,10 +185,16 @@ class Geometry(object):
         bind it to the appropriate textures, and return a list
         of the texture references.
         """
-        mesh = np.empty(self.mesh.shape[0]*3, dtype=gpuarray.vec.float4)
-        mesh['x'] = self.mesh[:,:,0].flatten()
-        mesh['y'] = self.mesh[:,:,1].flatten()
-        mesh['z'] = self.mesh[:,:,2].flatten()
+        vertices = np.empty(len(self.mesh.vertices), dtype=gpuarray.vec.float4)
+        vertices['x'] = self.mesh.vertices[:,0]
+        vertices['y'] = self.mesh.vertices[:,1]
+        vertices['z'] = self.mesh.vertices[:,2]
+
+        triangles = \
+            np.empty(len(self.mesh.triangles), dtype=gpuarray.vec.uint4)
+        triangles['x'] = self.mesh.triangles[:,0]
+        triangles['y'] = self.mesh.triangles[:,1]
+        triangles['z'] = self.mesh.triangles[:,2]
 
         lower_bounds = np.empty(self.lower_bounds.shape[0], dtype=gpuarray.vec.float4)
         lower_bounds['x'] = self.lower_bounds[:,0]
@@ -193,28 +206,32 @@ class Geometry(object):
         upper_bounds['y'] = self.upper_bounds[:,1]
         upper_bounds['z'] = self.upper_bounds[:,2]
 
-        self.mesh_gpu = cuda.to_device(mesh)
+        self.vertices_gpu = cuda.to_device(vertices)
+        self.triangles_gpu = cuda.to_device(triangles)
         self.lower_bounds_gpu = cuda.to_device(lower_bounds)
         self.upper_bounds_gpu = cuda.to_device(upper_bounds)
         self.node_map_gpu = cuda.to_device(self.node_map)
         self.node_length_gpu = cuda.to_device(self.node_length)
 
-        mesh_tex = module.get_texref('mesh')
+        vertices_tex = module.get_texref('vertices')
+        triangles_tex = module.get_texref('triangles')
         lower_bounds_tex = module.get_texref('lower_bounds')
         upper_bounds_tex = module.get_texref('upper_bounds')
         node_map_tex = module.get_texref('node_map')
         node_length_tex = module.get_texref('node_length')
 
-        mesh_tex.set_address(self.mesh_gpu, mesh.nbytes)
+        vertices_tex.set_address(self.vertices_gpu, vertices.nbytes)
+        triangles_tex.set_address(self.triangles_gpu, triangles.nbytes)
         lower_bounds_tex.set_address(self.lower_bounds_gpu, lower_bounds.nbytes)
         upper_bounds_tex.set_address(self.upper_bounds_gpu, upper_bounds.nbytes)
         node_map_tex.set_address(self.node_map_gpu, self.node_map.nbytes)
         node_length_tex.set_address(self.node_length_gpu, self.node_length.nbytes)
 
-        mesh_tex.set_format(cuda.array_format.FLOAT, 4)
+        vertices_tex.set_format(cuda.array_format.FLOAT, 4)
+        triangles_tex.set_format(cuda.array_format.UNSIGNED_INT32, 4)
         lower_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
         upper_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
         node_map_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
         node_length_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
 
-        return [mesh_tex, lower_bounds_tex, upper_bounds_tex, node_map_tex, node_length_tex]
+        return [vertices_tex, triangles_tex, lower_bounds_tex, upper_bounds_tex, node_map_tex, node_length_tex]
diff --git a/gputhread.py b/gputhread.py
index 12ae145..f573916 100644
--- a/gputhread.py
+++ b/gputhread.py
@@ -54,7 +54,7 @@ class GPUThread(threading.Thread):
             triangles = dest[(dest != -1)]
 
             bincount = np.zeros(len(self.geometry.solids))
-            gpu_bincount = np.bincount(self.geometry.solid_index[triangles])
+            gpu_bincount = np.bincount(self.geometry.solid_id[triangles])
             bincount[:gpu_bincount.size] = gpu_bincount
 
             self.output.put(bincount)
diff --git a/materials.py b/materials.py
index 9160437..f71a583 100644
--- a/materials.py
+++ b/materials.py
@@ -22,3 +22,6 @@ air = Material('air')
 h2o = Material('h2o')
 glass = Material('glass')
 vacuum = Material('vacuum')
+
+photocathode = Surface('photocathode')
+aluminum = Surface('aluminum')
diff --git a/mesh.py b/mesh.py
index a4b9048..a3efae5 100644
--- a/mesh.py
+++ b/mesh.py
@@ -26,11 +26,8 @@ class Mesh(object):
     def build(self):
         return self.vertices[self.triangles]
 
-    def __add__(self, other):
-        vertices = np.concatenate((self.vertices, other.vertices))
-        triangles = np.concatenate((self.triangles, other.triangles+len(self.vertices)))
-
-        return Mesh(vertices, triangles)
+    def __getitem__(self, key):
+        return self.vertices[self.triangles[key]]
 
     def __len__(self):
         return len(self.triangles)
diff --git a/photon.py b/photon.py
index f8ece6a..b494384 100644
--- a/photon.py
+++ b/photon.py
@@ -9,7 +9,7 @@ def uniform_sphere(size=None, dtype=np.double):
             Number of points to generate. If no size is specified, a single
             point is returned.
 
-    source: Weisstein, Eric W. "Sphere Point Picking." Mathworld. 
+    Source: Weisstein, Eric W. "Sphere Point Picking." Mathworld. 
     """
 
     theta, u = np.random.uniform(0.0, 2*np.pi, size), \
diff --git a/solid.py b/solid.py
index fe22463..6bf3790 100644
--- a/solid.py
+++ b/solid.py
@@ -54,5 +54,5 @@ class Solid(object):
         else:
             self.color = np.tile(color, len(self.mesh))
 
-    def build(self):
-        return np.inner(self.mesh.build(), self.rotation) + self.displacement
+    def __len__(self):
+        return len(self.mesh)
diff --git a/src/kernel.cu b/src/kernel.cu
index 5794057..a020180 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -10,7 +10,8 @@
 #define STACK_SIZE 500
 
 /* flattened triangle mesh */
-texture<float4, 1, cudaReadModeElementType> mesh;
+texture<float4, 1, cudaReadModeElementType> vertices;
+texture<uint4, 1, cudaReadModeElementType> triangles;
 
 /* lower/upper bounds for the bounding box associated with each node/leaf */
 texture<float4, 1, cudaReadModeElementType> upper_bounds;
@@ -79,11 +80,11 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 		{
 			for (i=0; i < length; i++)
 			{
-				int mesh_idx = 3*(index + i);
+				uint4 triangle_data = tex1Dfetch(triangles, index+i);
 
-				float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
-				float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
-				float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
+				float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
+				float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y));
+				float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z));
 
 				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
 				{
@@ -183,11 +184,11 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int
 	}
 	else
 	{
-		int mesh_idx = 3*intersection_idx;
+		uint4 triangle_data = tex1Dfetch(triangles, intersection_idx);
 
-		float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
-		float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
-		float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
+		float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
+		float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y));
+		float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z));
 
 		*(pixels+idx) = get_color(direction, v0, v1, v2);
 	}
diff --git a/transform.py b/transform.py
index ef8c96c..4793c0f 100644
--- a/transform.py
+++ b/transform.py
@@ -17,8 +17,4 @@ def rotate(x, phi, n):
     Rotate an array of points `x` through an angle phi counter-clockwise
     around the axis `n` (when looking towards +infinity).
     """
-    x = np.asarray(x)
-
-    r = make_rotation_matrix(phi, n)
-
-    return np.inner(x,r)
+    return np.inner(np.asarray(x),make_rotation_matrix(phi, n))
diff --git a/view.py b/view.py
index 1859195..2007973 100755
--- a/view.py
+++ b/view.py
@@ -26,13 +26,13 @@ def view(geometry, name=''):
         - move: shift+click and drag the mouse
     """
 
-    lower_bound = np.array([np.min(geometry.mesh[:,:,0]),
-                            np.min(geometry.mesh[:,:,1]),
-                            np.min(geometry.mesh[:,:,2])])
+    lower_bound = np.array([np.min(geometry.mesh[:][:,:,0]),
+                            np.min(geometry.mesh[:][:,:,1]),
+                            np.min(geometry.mesh[:][:,:,2])])
 
-    upper_bound = np.array([np.max(geometry.mesh[:,:,0]),
-                            np.max(geometry.mesh[:,:,1]),
-                            np.max(geometry.mesh[:,:,2])])
+    upper_bound = np.array([np.max(geometry.mesh[:][:,:,0]),
+                            np.max(geometry.mesh[:][:,:,1]),
+                            np.max(geometry.mesh[:][:,:,2])])
 
     scale = np.linalg.norm(upper_bound-lower_bound)