diff options
| -rw-r--r-- | detectors/lbne.py | 2 | ||||
| -rw-r--r-- | detectors/sno.py | 16 | ||||
| -rw-r--r-- | itertoolset.py | 4 | ||||
| -rw-r--r-- | make.py | 144 | ||||
| -rw-r--r-- | solids/__init__.py | 24 | ||||
| -rw-r--r-- | solids/pmts.py | 18 | ||||
| -rw-r--r-- | src/intersect.h | 11 | ||||
| -rw-r--r-- | src/photon.h | 5 | ||||
| -rwxr-xr-x | view.py | 2 |
9 files changed, 135 insertions, 91 deletions
diff --git a/detectors/lbne.py b/detectors/lbne.py index d3c12da..8942cb2 100644 --- a/detectors/lbne.py +++ b/detectors/lbne.py @@ -21,7 +21,7 @@ def build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing, physic lbne = Geometry() # outer cylinder - cylinder_mesh = segmented_cylinder(radius, height+height/(pmts_per_string-1), theta=(2*np.pi/nstrings)/8, n=200) + cylinder_mesh = segmented_cylinder(radius, height+height/(pmts_per_string-1), nsteps=16*nstrings, nsegments=1200) cylinder_mesh.vertices = rotate(cylinder_mesh.vertices, np.pi/2, (-1,0,0)) lbne.add_solid(Solid(cylinder_mesh, water_wcsim, vacuum, black_surface, 0xff0000ff)) diff --git a/detectors/sno.py b/detectors/sno.py index 156dcbe..3e748b4 100644 --- a/detectors/sno.py +++ b/detectors/sno.py @@ -9,7 +9,7 @@ from solids import build_8inch_pmt_with_lc from transform import rotate, make_rotation_matrix import os -def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20): +def sno_vessel(sphere_radius, neck_radius, neck_top, nsteps=40): '''Compute the 2D coordinates of the profile of one side of a SNO-style acrylic vessel. The center of the sphere is at (0,0), with the neck extending along the positive y direction. @@ -17,7 +17,7 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20): sphere_radius: Radius of spherical part of profile neck_radius: Radius of neck part neck_top: y coordinate of top of neck - angle_step: angular step size (radius) between points on the sphere + nsteps: number of points around the circumference of the sphere Returns: Tuple of x and y coordinate numpy arrays. ''' @@ -31,7 +31,7 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20): raise ValueError('neck_top must be greater than the y-value where the sphere and cylinder intersect') # Start with point at bottom - angles = np.arange(-math.pi/2, max_angle, angle_step) + angles = np.linspace(-math.pi/2, max_angle, nsteps) x = list(np.cos(angles) * sphere_radius) y = list(np.sin(angles) * sphere_radius) x[0] = 0.0 # Round-off error might make cos(-pi/2) not exactly zero @@ -52,19 +52,19 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20): ##### SNO Parts -angle_step = math.pi/20 +nsteps = 40 av_outside_profile = sno_vessel(sphere_radius=6.0604, neck_radius=0.79375, - neck_top=10.50, angle_step=angle_step) + neck_top=10.50, nsteps=nsteps) # For simplicity, cap the top of the AV with acrylic av_inside_profile = sno_vessel(sphere_radius=6.0053, neck_radius=0.72898, - neck_top=10.00, angle_step=angle_step) + neck_top=10.00, nsteps=nsteps) av_outside_mesh = rotate_extrude(av_outside_profile[0], av_outside_profile[1], - angle_step) + nsteps) av_outside_mesh.vertices = rotate(av_outside_mesh.vertices, np.pi/2, (-1,0,0)) av_inside_mesh = rotate_extrude(av_inside_profile[0], av_inside_profile[1], - angle_step) + nsteps) av_inside_mesh.vertices = rotate(av_inside_mesh.vertices, np.pi/2, (-1,0,0)) dir = os.path.split(os.path.realpath(__file__))[0] diff --git a/itertoolset.py b/itertoolset.py index 9d5200c..a20d03b 100644 --- a/itertoolset.py +++ b/itertoolset.py @@ -47,3 +47,7 @@ def consume(iterator, n=None): else: # advance to the empty slice starting at position n next(islice(iterator, n, n), None) + +def flatten(listOfLists): + "Flatten one level of nesting" + return chain.from_iterable(listOfLists) @@ -1,67 +1,109 @@ import numpy as np from geometry import Mesh from transform import rotate +from itertoolset import * -def rotate_extrude(x, y, theta=np.pi/32, remove_duplicate_vertices=True): - x, y, = np.asarray(x), np.asarray(y) +def mesh_grid(grid): + return np.vstack(zip(grid[:-1].flatten(),grid[1:].flatten(),np.roll(grid[1:],-1,1).flatten()) + zip(grid[:-1].flatten(),np.roll(grid[1:],-1,1).flatten(),np.roll(grid[:-1],-1,1).flatten())) - if len(x.shape) != 1 or len(y.shape) != 1 or x.size != y.size: - raise ValueError('shape mismatch') +def linear_extrude(x1, y1, height, x2=None, y2=None): + """ + Return the solid mesh formed by linearly extruding the polygon formed by + the x and y points `x1` and `y1` by a distance `height`. If `x2` and `y2` + are given extrude by connecting the points `x1` and `y1` to `x2` and `y2`; + this allows the creation of tapered solids. - points = np.zeros((len(x),3), dtype=np.float32) - points[:,0] = x - points[:,1] = y + .. note:: + The path traced by the points `x` and `y` should go counter-clockwise, + otherwise the mesh will be inside out. - angles = np.arange(0, 2*np.pi, theta) + Example: + >>> # create a hexagon prism + >>> angles = np.linspace(0, 2*np.pi, 6, endpoint=False) + >>> m = linear_extrude(np.cos(angles), np.sin(angles), 5.0) + """ + if len(x1) != len(y1): + raise Exception('`x` and `y` arrays must have the same length.') - vertices = np.empty((len(points)*len(angles), 3), dtype=np.float32) - triangles = np.zeros((len(angles)*(len(points)-1)*2, 3), dtype=np.int32) + if x2 is None: + x2 = x1 - step = len(points) - 1 + if y2 is None: + y2 = y1 - for i, angle in enumerate(angles): - this_slice = i*len(points) - vertices[this_slice:this_slice+len(points)] = rotate(points, angle, (0,-1,0)) + if len(x2) != len(y2) or len(x2) != len(x1): + raise Exception('`x` and `y` arrays must have the same length.') - start = 2*i*step - next_slice = ((i+1) % len(angles))*len(points) + n = len(x1) - triangles[start:start+step,0] = np.arange(this_slice, this_slice+len(points)-1) - triangles[start:start+step,1] = np.arange(next_slice, next_slice+len(points)-1) - triangles[start:start+step,2] = np.arange(this_slice+1, this_slice+len(points)) + vertex_iterators = [izip(repeat(0,n),repeat(0,n),repeat(-height/2.0,n)),izip(x1,y1,repeat(-height/2.0,n)),izip(x2,y2,repeat(height/2.0,n)),izip(repeat(0,n),repeat(0,n),repeat(height/2.0,n))] - triangles[start+step:start+2*step,0] = np.arange(next_slice, next_slice+len(points)-1) - triangles[start+step:start+2*step,2] = np.arange(this_slice+1, this_slice+len(points)) - triangles[start+step:start+2*step,1] = np.arange(next_slice+1, next_slice+len(points)) + vertices = np.fromiter(flatten(roundrobin(*vertex_iterators)), float) + vertices = vertices.reshape((len(vertices)//3,3)) - return Mesh(vertices, triangles, remove_duplicate_vertices=remove_duplicate_vertices) + triangles = mesh_grid(np.arange(len(vertices)).reshape((len(x1),len(vertices)//len(x1))).transpose()[::-1]) + + return Mesh(vertices, triangles, remove_duplicate_vertices=True) + +def rotate_extrude(x, y, nsteps=64): + """ + Return the solid mesh formed by extruding the profile defined by the x and + y points `x` and `y` around the y axis. + + .. note:: + The path traced by the points `x` and `y` should go counter-clockwise, + otherwise the mesh will be inside out. + """ + if len(x) != len(y): + raise Exception('`x` and `y` arrays must have the same length.') + + points = np.array([x,y,np.zeros(len(x))]).transpose() + + steps = np.linspace(0, 2*np.pi, nsteps, endpoint=False) + vertices = np.vstack([rotate(points,angle,(0,-1,0)) for angle in steps]) + triangles = mesh_grid(np.arange(len(vertices)).reshape((len(steps),len(points))).transpose()[::-1]) + + return Mesh(vertices, triangles, remove_duplicate_vertices=True) def cube(size=1): - a = np.sqrt(size**2/2.0) - x = [0,a,a,0] - y = [-size/2.0,-size/2.0,size/2.0,size/2.0] - return rotate_extrude(x, y, np.pi/2) - -def cylinder(radius1=1, radius2=1, height=2, theta=np.pi/32): - x = [0,radius1,radius2,0] - y = [-height/2.0, -height/2.0, height/2.0, height/2.0] - return rotate_extrude(x, y, theta) - -def segmented_cylinder(radius, height=2, theta=np.pi/32, n=50): - x = np.concatenate((np.linspace(0, radius, n, endpoint=False), - [radius] * n, - np.linspace(radius, 0, n, endpoint=False), - [0.0])) - y = np.concatenate(([-height/2.0] * n, - np.linspace(-height/2.0, height/2.0, n, endpoint=False), - [height/2.0] * (n+1))) - return rotate_extrude(x, y, theta) - -def sphere(radius=1, theta=np.pi/32): - profile_angles = np.arange(-np.pi/2, np.pi/2+theta, theta) - return rotate_extrude(radius*np.cos(profile_angles), radius*np.sin(profile_angles), theta) - -def torus(radius=1, offset=3, phi=np.pi/32, theta=np.pi/32): - profile_angles = np.arange(0, 2*np.pi+phi, phi) - x, y = np.cos(profile_angles), np.sin(profile_angles) - return rotate_extrude(x + offset, y) + "Return a cube mesh whose sides have length `size`." + return linear_extrude([-size/2.0,size/2.0,size/2.0,-size/2.0],[-size/2.0,-size/2.0,size/2.0,size/2.0], height=size) + +def cylinder(radius=1, height=2, radius2=None, nsteps=64): + """ + Return a cylinder mesh with a radius of length `radius`, and a height of + length `height`. If `radius2` is specified, return a cone shaped cylinder + with bottom radius `radius`, and top radius `radius2`. + """ + if radius2 is None: + radius2 = radius + + return rotate_extrude([0,radius,radius2,0], [-height/2.0, -height/2.0, height/2.0, height/2.0], nsteps) + +def segmented_cylinder(radius, height=2, nsteps=64, nsegments=100): + """ + Return a cylinder mesh segmented into `nsegments` points along its profile. + """ + nsegments_radius = int((nsegments*radius/(2*radius+height))/2) + nsegments_height = int((nsegments*height/(2*radius+height))/2) + x = np.concatenate([np.linspace(0,radius,nsegments_radius,endpoint=False),[radius]*nsegments_height,np.linspace(radius,0,nsegments_radius,endpoint=False),[0]]) + y = np.concatenate([[-height/2.0]*nsegments_radius,np.linspace(-height/2.0,height/2.0,nsegments_height,endpoint=False),[height/2.0]*(nsegments_radius+1)]) + return rotate_extrude(x, y, nsteps) + +def sphere(radius=1, nsteps=64): + "Return a sphere mesh." + profile_angles = np.linspace(-np.pi/2, np.pi/2, nsteps) + return rotate_extrude(radius*np.cos(profile_angles), radius*np.sin(profile_angles), nsteps) + +def torus(radius=1, offset=3, nsteps=64, circle_steps=None): + """ + Return a torus mesh. `offset` is the distance from the center of the torus + to the center of the torus barrel. `radius` is the radius of the torus + barrel. `nsteps` is the number of steps in the rotational extrusion of the + circle. `circle_steps` if specified is the number of steps around the + circumference of the torus barrel, else it defaults to `nsteps`. + """ + if circle_steps is None: + circle_steps = nsteps + profile_angles = np.linspace(0, 2*np.pi, circle_steps) + return rotate_extrude(np.cos(profile_angles) + offset, np.sin(profile_angles), nsteps) diff --git a/solids/__init__.py b/solids/__init__.py index 46f3119..a24ce9c 100644 --- a/solids/__init__.py +++ b/solids/__init__.py @@ -12,11 +12,11 @@ from optics import * from view import buildable @buildable('12inch_pmt') -def build_12inch_pmt(outer_material=water, theta=np.pi/8): - return build_pmt(dir + '/hamamatsu_12inch.txt', 0.003, outer_material, theta) +def build_12inch_pmt(outer_material=water, nsteps=16): + return build_pmt(dir + '/hamamatsu_12inch.txt', 0.003, outer_material, nsteps) @buildable('12inch_pmt_shell') -def build_12inch_pmt_shell(outer_material=water, theta=np.pi/8): +def build_12inch_pmt_shell(outer_material=water, nsteps=16): return build_pmt_shell(dir + '/hamamatsu_12inch.txt') # from Jelena Maricic @@ -27,22 +27,22 @@ lc_12inch_rmin = 0.1524 lc_12inch_rmax = 0.235072 @buildable('12inch_pmt_with_lc') -def build_12inch_pmt_with_lc(outer_material=water, theta=np.pi/8): - pmt = build_12inch_pmt(outer_material, theta) +def build_12inch_pmt_with_lc(outer_material=water, nsteps=16): + pmt = build_12inch_pmt(outer_material, nsteps) return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax) @buildable('12inch_pmt_with_lc_hd') -def build_12inch_pmt_with_lc_hd(outer_material=water, theta=np.pi/64): - pmt = build_12inch_pmt(outer_material, theta) +def build_12inch_pmt_with_lc_hd(outer_material=water, nsteps=128): + pmt = build_12inch_pmt(outer_material, nsteps) return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax, npoints=100) @buildable('8inch_pmt') -def build_8inch_pmt(outer_material=water, theta=np.pi/12): - return build_pmt(dir + '/sno_pmt_reduced.txt', 0.003, outer_material, theta) +def build_8inch_pmt(outer_material=water, nsteps=24): + return build_pmt(dir + '/sno_pmt_reduced.txt', 0.003, outer_material, nsteps) @buildable('8inch_pmt_with_lc') -def build_8inch_pmt_with_lc(outer_material=water, theta=np.pi/12): - pmt = build_8inch_pmt(outer_material, theta) - lc = build_light_collector_from_file(dir + '/sno_cone.txt', outer_material, theta) +def build_8inch_pmt_with_lc(outer_material=water, nsteps=24): + pmt = build_8inch_pmt(outer_material, nsteps) + lc = build_light_collector_from_file(dir + '/sno_cone.txt', outer_material, nsteps) return pmt + lc diff --git a/solids/pmts.py b/solids/pmts.py index a30d007..fcbc8fc 100644 --- a/solids/pmts.py +++ b/solids/pmts.py @@ -25,11 +25,11 @@ def build_light_collector(pmt, a, b, d, rmin, rmax, npoints=10): lc_offset = np.interp(lc_radii[0], list(reversed(pmt_face_profile[:,0])), list(reversed(pmt_face_profile[:,1]))) - lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, pmt.theta) + lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, pmt.nsteps) return Solid(lc_mesh, pmt.outer_material, pmt.outer_material, surface=shiny_surface) -def build_pmt_shell(filename, outer_material=water, theta=np.pi/8): +def build_pmt_shell(filename, outer_material=water, nsteps=16): profile = read_csv(filename) # slice profile in half @@ -44,9 +44,9 @@ def build_pmt_shell(filename, outer_material=water, theta=np.pi/8): # convert mm -> m profile /= 1000.0 - return Solid(rotate_extrude(profile[:,0], profile[:,1], theta), glass, outer_material, color=0xeeffffff) + return Solid(rotate_extrude(profile[:,0], profile[:,1], nsteps), glass, outer_material, color=0xeeffffff) -def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8): +def build_pmt(filename, glass_thickness, outer_material=water, nsteps=16): profile = read_csv(filename) # slice profile in half @@ -63,8 +63,8 @@ def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8): offset_profile = offset(profile, -glass_thickness) - outer_envelope_mesh = rotate_extrude(profile[:,0], profile[:,1], theta) - inner_envelope_mesh = rotate_extrude(offset_profile[:,0], offset_profile[:,1], theta) + outer_envelope_mesh = rotate_extrude(profile[:,0], profile[:,1], nsteps) + inner_envelope_mesh = rotate_extrude(offset_profile[:,0], offset_profile[:,1], nsteps) outer_envelope = Solid(outer_envelope_mesh, glass, outer_material) @@ -78,16 +78,16 @@ def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8): # light collector pmt.profile = profile pmt.outer_material = outer_material - pmt.theta = theta + pmt.nsteps = nsteps return pmt -def build_light_collector_from_file(filename, outer_material, theta=np.pi/24): +def build_light_collector_from_file(filename, outer_material, nsteps=48): profile = read_csv(filename) # Convert mm to m profile /= 1000.0 - mesh = rotate_extrude(profile[:,0], profile[:,1], theta) + mesh = rotate_extrude(profile[:,0], profile[:,1], nsteps) solid = Solid(mesh, outer_material, outer_material, surface=shiny_surface) return solid diff --git a/src/intersect.h b/src/intersect.h index d2f54ce..26e1d7e 100644 --- a/src/intersect.h +++ b/src/intersect.h @@ -60,17 +60,20 @@ __device__ bool intersect_triangle(const float3 &origin, const float3 &direction angle between the ray and the plane normal to determine the brightness. `direction` must be normalized. */ -__device__ unsigned int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, const unsigned int base_color=0xFFFFFFFF) +__device__ unsigned int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, unsigned int base_color=0xFFFFFF) { float scale = dot(normalize(cross(v1-v0,v2-v1)),-direction); + if (scale < 0.0f) + { + base_color = 0xff0000; + scale = dot(-normalize(cross(v1-v0,v2-v1)),-direction); + } + unsigned int r = 0xFF & (base_color >> 16); unsigned int g = 0xFF & (base_color >> 8); unsigned int b = 0xFF & base_color; - if (scale < 0.0f) - scale = dot(-normalize(cross(v1-v0,v2-v1)),-direction); - r = floorf(r*scale); g = floorf(g*scale); b = floorf(b*scale); diff --git a/src/photon.h b/src/photon.h index fa3b85b..f76ca54 100644 --- a/src/photon.h +++ b/src/photon.h @@ -103,9 +103,6 @@ __device__ void fill_state(State &s, Photon &p) s.absorption_length = interp_property(p.wavelength, material1.absorption_length); s.scattering_length = interp_property(p.wavelength, material1.scattering_length); - printf("wavelength = %f", p.wavelength); - printf("scattering length = %f\n", s.scattering_length); - } // fill_state __device__ void rayleigh_scatter(Photon &p, curandState &rng) @@ -156,8 +153,6 @@ __device__ int propagate_to_boundary(Photon &p, State &s, curandState &rng) { if (scattering_distance <= s.distance_to_boundary) { - printf("scattering distance = %f\n", scattering_distance); - p.time += scattering_distance/(SPEED_OF_LIGHT/s.refractive_index1); p.position += scattering_distance*p.direction; @@ -206,7 +206,7 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False): cuda.Context.synchronize() elapsed = time.time() - t0 - print 'elapsed %f sec' % elapsed + #print 'elapsed %f sec' % elapsed pygame.surfarray.blit_array(screen, pixels_gpu.get().reshape(size)) pygame.display.flip() |