Adapt existing channel example to new SDF-based rendering and voxelization

5 files changed, 224 insertions, 211 deletions

diff --git a/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py b/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py
index 166c36d..8806029 100644
--- a/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py
+++ b/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py
@@ -22,9 +22,9 @@ from utility.projection import Projection, Rotation
 from utility.opengl     import MomentsTexture
 from utility.mouse      import MouseDragMonitor, MouseScrollMonitor
 
-lattice_x = 200
-lattice_y = 90
-lattice_z = 90
+lattice_x = 180
+lattice_y = 100
+lattice_z = 100
 
 updates_per_frame = 5
 
@@ -52,47 +52,52 @@ boundary = Template("""
 )
 
 grid_fin = """
-v = rotate_z(translate(v, v3(center.x/2.5, center.y, center.z)), -0.5);
+v = rotate_z(translate(v, v3(center.x/2, center.y, center.z)), -0.8);
 float width = 1;
-float rotat = 0.64;
-
-return rounded(unify(
-	sub(
-		box(v, v3(5, 30, 40)),
-		box(v, v3(6, 28, 38))
-	),
-	intersect(
-		box(v, v3(5, 30, 40)),
-		unify(
-			unify(
-				box(rotate_x(v, rotat), v3(10, width, 100)),
-				box(rotate_x(v, -rotat), v3(10, width, 100))
-			),
-			unify(
-				unify(
-					unify(
-						box(rotate_x(translate(v, v3(0,0,25)), rotat), v3(10, width, 100)),
-						box(rotate_x(translate(v, v3(0,0,25)), -rotat), v3(10, width, 100))
-					),
-					unify(
-						box(rotate_x(translate(v, v3(0,0,-25)), rotat), v3(10, width, 100)),
-						box(rotate_x(translate(v, v3(0,0,-25)), -rotat), v3(10, width, 100))
-					)
-				),
-				unify(
-					unify(
-						box(rotate_x(translate(v, v3(0,0,50)), rotat), v3(10, width, 100)),
-						box(rotate_x(translate(v, v3(0,0,50)), -rotat), v3(10, width, 100))
-					),
-					unify(
-						box(rotate_x(translate(v, v3(0,0,-50)), rotat), v3(10, width, 100)),
-						box(rotate_x(translate(v, v3(0,0,-50)), -rotat), v3(10, width, 100))
-					)
-				)
-			)
-		)
-	)
-), 0.4);
+float angle = 0.64;
+
+return add(
+    sadd(
+        sub(
+            rounded(box(v, v3(5, 28, 38)), 1),
+            rounded(box(v, v3(6, 26, 36)), 1)
+        ),
+        cylinder(translate(v, v3(0,0,-45)), 5, 12),
+        1
+    ),
+    sintersect(
+        box(v, v3(5, 28, 38)),
+        add(
+            add(
+                box(rotate_x(v, angle), v3(10, width, 100)),
+                box(rotate_x(v, -angle), v3(10, width, 100))
+            ),
+            add(
+                add(
+                    add(
+                        box(rotate_x(translate(v, v3(0,0,25)), angle), v3(10, width, 100)),
+                        box(rotate_x(translate(v, v3(0,0,25)), -angle), v3(10, width, 100))
+                    ),
+                    add(
+                        box(rotate_x(translate(v, v3(0,0,-25)), angle), v3(10, width, 100)),
+                        box(rotate_x(translate(v, v3(0,0,-25)), -angle), v3(10, width, 100))
+                    )
+                ),
+                add(
+                    add(
+                        box(rotate_x(translate(v, v3(0,0,50)), angle), v3(10, width, 100)),
+                        box(rotate_x(translate(v, v3(0,0,50)), -angle), v3(10, width, 100))
+                    ),
+                    add(
+                        box(rotate_x(translate(v, v3(0,0,-50)), angle), v3(10, width, 100)),
+                        box(rotate_x(translate(v, v3(0,0,-50)), -angle), v3(10, width, 100))
+                    )
+                )
+            )
+        ),
+        2
+    )
+);
 """
 
 def glut_window(fullscreen = False):
@@ -179,6 +184,10 @@ vec3 fabs(vec3 x) {
     return abs(x);
 }
 
+float fabs(float x) {
+    return abs(x);
+}
+
 <%include file="template/sdf.lib.glsl.mako"/>
 
 float sdf(vec3 v) {
@@ -209,29 +218,28 @@ float maxRayLength(vec3 origin, vec3 ray) {
     return max(1.0, ${max_ray_length} - distanceToLattice(origin + ${max_ray_length}*ray));
 }
 
-vec4 trace_obstacle(vec3 pos, vec3 ray) {
-    bool hit = false;
-    vec3 color  = vec3(0.05);
-    vec3 sample_pos = pos;
+vec4 trace_obstacle(vec3 origin, vec3 ray, float delta) {
+    vec3 color      = vec3(0);
+    vec3 sample_pos = origin;
+    float ray_dist = 0.0;
 
     for (int i = 0; i < OBSTACLE_STEPS; ++i) {
-        float dist = sdf(sample_pos);
-        if (abs(dist) < EPSILON) {
-            vec3 n = normalize(sdf_normal(sample_pos));
-            float brightness = abs(dot(n, ray));
-            color += brightness;
-            hit = true;
-            break;
+        const float sdf_dist = sdf(sample_pos);
+        ray_dist += sdf_dist;
+
+        if (ray_dist > delta) {
+            return vec4(0.0);
+        }
+
+        if (abs(sdf_dist) < EPSILON) {
+            const vec3 n = normalize(sdf_normal(sample_pos));
+            return vec4(color + abs(dot(n, ray)), 1.0);
         } else {
-            sample_pos += dist*ray;
+            sample_pos = origin + ray_dist*ray;
         }
     }
 
-    if (hit) {
-        return vec4(color, 1.0);
-    } else {
-        return vec4(color, 0.0);
-    }
+    return vec4(0.0);
 }
 
 vec3 trace(vec3 pos, vec3 ray) {
@@ -240,17 +248,16 @@ vec3 trace(vec3 pos, vec3 ray) {
     vec3 color = vec3(0.0);
 
     for (int i=0; i < RAYMARCH_STEPS; ++i) {
-        const vec3 pos = pos + i*delta*ray;
-        const vec4 data = texture(moments, unit(pos));
-        if (sdf(pos) > delta) {
-            color += 1./RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
+        const vec3 sample_pos = pos + i*delta*ray;
+        const vec4 data = texture(moments, unit(sample_pos));
+        if (sdf(sample_pos) > delta) {
+            color += 0.5/RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
         } else {
-            vec4 obstacle_color = trace_obstacle(pos, ray);
+            const vec4 obstacle_color = trace_obstacle(sample_pos, ray, delta);
             if (obstacle_color.w == 1.0) {
-                color += obstacle_color.xyz;
-                break;
+                return color + obstacle_color.xyz;
             } else {
-                color += 1./RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
+                color += 0.5/RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
             }
         }
     }
diff --git a/channel_3d_volumetric_rendering_gl_interop.py b/channel_3d_volumetric_rendering_gl_interop.py
index 7c78d9f..1b149da 100644
--- a/channel_3d_volumetric_rendering_gl_interop.py
+++ b/channel_3d_volumetric_rendering_gl_interop.py
@@ -1,5 +1,9 @@
 import numpy
-from string import Template
+
+from mako.template import Template
+from mako.lookup import TemplateLookup
+
+from pathlib import Path
 
 from simulation         import Lattice, Geometry
 from utility.particles  import Particles
@@ -12,78 +16,76 @@ from OpenGL.GLUT import *
 
 from OpenGL.GL import shaders
 
-from geometry.sphere   import Sphere
-from geometry.box      import Box
-from geometry.cylinder import Cylinder
+from geometry.box import Box
 
 from utility.projection import Projection, Rotation
 from utility.opengl     import MomentsTexture
 from utility.mouse      import MouseDragMonitor, MouseScrollMonitor
 
-lattice_x = 256
+lattice_x = 200
 lattice_y = 64
 lattice_z = 64
 
-updates_per_frame = 5
+updates_per_frame = 8
 
 inflow = 0.01
 relaxation_time = 0.51
 
 lbm = LBM(D3Q27)
 
-def get_cavity_material_map(g):
-    return [
-        (lambda x, y, z: x > 0 and x < g.size_x-1 and
-                         y > 0 and y < g.size_y-1 and
-                         z > 0 and z < g.size_z-1,     1), # bulk fluid
-        (lambda x, y, z: x == 1 or x == g.size_x-2 or
-                         y == 1 or y == g.size_y-2 or
-                         z == 1 or z == g.size_z-2,    2), # walls
-        (lambda x, y, z: x == 1,                       3), # inflow
-        (lambda x, y, z: x == g.size_x-2,              4), # outflow
-
-        # wall with hole
-        (Box(2*g.size_x//20, 2.5*g.size_x//20, 0, g.size_y, 0, g.size_z), 5),
-        (Sphere(2.5*g.size_x//20, g.size_y//2, g.size_z//2, 10),          1),
-
-        (Box(6.0*g.size_x//20, 6.5*g.size_x//20, 0, g.size_y, 0, g.size_z), 5),
-        (Sphere(6.5*g.size_x//20, 2*g.size_y//3, 2*g.size_z//3, 10),         1),
-        (Sphere(6.5*g.size_x//20, 1*g.size_y//3, 1*g.size_z//3, 10),         1),
-
-        (Box(10.0*g.size_x//20, 10.5*g.size_x//20, 0, g.size_y, 0, g.size_z), 5),
-        (Sphere(10.5*g.size_x//20, 1*g.size_y//3, 2*g.size_z//3, 10),         1),
-        (Sphere(10.5*g.size_x//20, 2*g.size_y//3, 1*g.size_z//3, 10),         1),
-
-        (Box(14.0*g.size_x//20, 14.5*g.size_x//20, 0, g.size_y, 0, g.size_z), 5),
-        (Sphere(14.5*g.size_x//20, 2*g.size_y//3, 2*g.size_z//3, 10),         1),
-        (Sphere(14.5*g.size_x//20, 1*g.size_y//3, 1*g.size_z//3, 10),         1),
-
-        (lambda x, y, z: x == 0 or x == g.size_x-1 or
-                         y == 0 or y == g.size_y-1 or
-                         z == 0 or z == g.size_z-1,   0)  # ghost cells
-    ]
-
 boundary = Template("""
-    if ( m == 2 || m == 5 ) {
+    if ( m == 2 ) {
         u_0 = 0.0;
         u_1 = 0.0;
         u_2 = 0.0;
     }
     if ( m == 3 ) {
-        u_0 = min(time/5000.0 * $inflow, $inflow);
+        u_0 = min(time/5000.0 * ${inflow}, ${inflow});
         u_1 = 0.0;
         u_2 = 0.0;
     }
     if ( m == 4 ) {
         rho = 1.0;
     }
-""").substitute({
-    "inflow": inflow
-})
+""").render(
+    inflow = inflow
+)
+
+channel = """
+return add(
+    add(
+        ssub(
+            box(translate(v, v3(25,center.y,center.z)), v3(5,32,32)),
+            add(
+                sphere(translate(v, v3(20,0.5*center.y,1.5*center.z)), 10),
+                sphere(translate(v, v3(30,1.5*center.y,0.5*center.z)), 10)
+            ),
+            2
+        ),
+        ssub(
+            box(translate(v, v3(85,center.y,center.z)), v3(5,32,32)),
+            add(
+                sphere(translate(v, v3(90,1.5*center.y,1.5*center.z)), 10),
+                sphere(translate(v, v3(80,0.5*center.y,0.5*center.z)), 10)
+            ),
+            2
+        )
+    ),
+    ssub(
+        box(translate(v, v3(145,center.y,center.z)), v3(5,32,32)),
+        add(
+            cylinder(rotate_y(translate(v, v3(145,1.5*center.y,0.5*center.z)), 1), 10, 10),
+            cylinder(rotate_y(translate(v, v3(145,0.5*center.y,1.5*center.z)), -1), 10, 10)
+        ),
+        2
+    )
+);
+"""
 
 def glut_window(fullscreen = False):
     glutInit(sys.argv)
-    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH | GLUT_MULTISAMPLE)
+    glutSetOption(GLUT_MULTISAMPLE, 8)
+    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_MULTISAMPLE)
 
     if fullscreen:
         window = glutEnterGameMode()
@@ -96,19 +98,17 @@ def glut_window(fullscreen = False):
 
 window = glut_window(fullscreen = False)
 
-vertex_shader = shaders.compileShader(Template("""
+vertex_shader = shaders.compileShader("""
 #version 430
 
 layout (location=0) in vec4 vertex;
-                   out vec3 color;
 
 uniform mat4 projection;
 uniform mat4 rotation;
 
 void main() {
     gl_Position = projection * rotation * vertex;
-    color = vec3(1.0,1.0,1.0);
-}""").substitute({}), GL_VERTEX_SHADER)
+}""", GL_VERTEX_SHADER)
 
 fragment_shader = shaders.compileShader("""
 #version 430
@@ -116,7 +116,7 @@ fragment_shader = shaders.compileShader("""
 in vec3 color;
 
 void main(){
-    gl_FragColor = vec4(color.xyz, 0.0);
+    gl_FragColor = vec4(vec3(0.5), 1.0);
 }""", GL_FRAGMENT_SHADER)
 
 raycast_vertex_shader = shaders.compileShader("""
@@ -145,7 +145,43 @@ uniform sampler3D moments;
 out vec4 result;
 
 vec3 unit(vec3 v) {
-    return vec3(v[0] / $size_x, v[1] / $size_y, v[2] / $size_z);
+    return vec3(v[0] / ${size_x}, v[1] / ${size_y}, v[2] / ${size_z});
+}
+
+const vec3 center = vec3(${size_x/2.5}, ${size_y/2}, ${size_z/2});
+
+#define EPSILON 1e-1
+#define RAYMARCH_STEPS 64
+#define OBSTACLE_STEPS 16
+
+vec3 v3(float x, float y, float z) {
+	return vec3(x,y,z);
+}
+
+vec2 v2(float x, float y) {
+	return vec2(x,y);
+}
+
+vec3 fabs(vec3 x) {
+    return abs(x);
+}
+
+float fabs(float x) {
+    return abs(x);
+}
+
+<%include file="template/sdf.lib.glsl.mako"/>
+
+float sdf(vec3 v) {
+    ${sdf_source}
+}
+
+vec3 sdf_normal(vec3 v) {
+    return normalize(vec3(
+        sdf(vec3(v.x + EPSILON, v.y, v.z)) - sdf(vec3(v.x - EPSILON, v.y, v.z)),
+        sdf(vec3(v.x, v.y + EPSILON, v.z)) - sdf(vec3(v.x, v.y - EPSILON, v.z)),
+        sdf(vec3(v.x, v.y, v.z + EPSILON)) - sdf(vec3(v.x, v.y, v.z - EPSILON))
+    ));
 }
 
 vec3 palette(float x) {
@@ -156,65 +192,55 @@ vec3 palette(float x) {
     );
 }
 
-vec3 blueBlackRedPalette(float x) {
-	if ( x < 0.5 ) {
-		return mix(
-			vec3(0.0, 0.0, 1.0),
-			vec3(0.0, 0.0, 0.0),
-			2*x
-		);
-	} else {
-		return mix(
-			vec3(0.0, 0.0, 0.0),
-			vec3(1.0, 0.0, 0.0),
-			2*(x - 0.5)
-		);
-	}
+float distanceToLattice(vec3 v) {
+  return box(v, vec3(${size_x},${size_y},${size_z}));
 }
 
-vec3 v(float x, float y, float z) {
-    return texture(moments, unit(vec3(x,y,z))).yzw;
+float maxRayLength(vec3 origin, vec3 ray) {
+    return max(1.0, ${max_ray_length} - distanceToLattice(origin + ${max_ray_length}*ray));
 }
 
-vec3 curl(vec3 p) {
-    const float h = 1./$size_x;
+vec4 trace_obstacle(vec3 origin, vec3 ray, float delta) {
+    vec3 color      = vec3(0);
+    vec3 sample_pos = origin;
+    float ray_dist = 0.0;
 
-    const float dyvz = (v(p.x,p.y+1,p.z).z - v(p.x,p.y-1,p.z).z) / (2.0*h);
-    const float dzvy = (v(p.x,p.y,p.z+1).y - v(p.x,p.y,p.z-1).y) / (2.0*h);
+    for (int i = 0; i < OBSTACLE_STEPS; ++i) {
+        const float sdf_dist = sdf(sample_pos);
+        ray_dist += sdf_dist;
 
-    const float dzvx = (v(p.x,p.y,p.z+1).x - v(p.x,p.y,p.z-1).x) / (2.0*h);
-    const float dxvz = (v(p.x+1,p.y,p.z).z - v(p.x-1,p.y,p.z).z) / (2.0*h);
+        if (ray_dist > delta) {
+            return vec4(0.0);
+        }
 
-    const float dxvy = (v(p.x+1,p.y,p.z).y - v(p.x-1,p.y,p.z).y) / (2.0*h);
-    const float dyvx = (v(p.x,p.y+1,p.z).x - v(p.x,p.y-1,p.z).x) / (2.0*h);
+        if (abs(sdf_dist) < EPSILON) {
+            const vec3 n = normalize(sdf_normal(sample_pos));
+            return vec4(color + abs(dot(n, ray)), 1.0);
+        } else {
+            sample_pos = origin + ray_dist*ray;
+        }
+    }
 
-    return vec3(
-        dyvz - dzvy,
-        dzvx - dxvz,
-        dxvy - dyvx
-    );
+    return vec4(0.0);
 }
 
 vec3 trace(vec3 pos, vec3 ray) {
-    const float ray_length = $max_ray_length;
-    const float delta      = 1./ray_length;
+    const float delta = maxRayLength(pos, ray) / RAYMARCH_STEPS;
 
-    vec3 color  = vec3(0.0);
+    vec3 color = vec3(0.0);
 
-    for (float t = 0.0; t < 1.0; t += delta) {
-        const vec3 sample_pos = unit(pos + t*ray_length*ray);
-        if (length(sample_pos) < sqrt(3.1)) {
-            const vec4 data = texture(moments, sample_pos);
-            if (data[3] != 1.0) {
-                color += delta * palette(length(data.yzw) / $inflow);
+    for (int i=0; i < RAYMARCH_STEPS; ++i) {
+        const vec3 sample_pos = pos + i*delta*ray;
+        const vec4 data = texture(moments, unit(sample_pos));
+        if (sdf(sample_pos) > delta) {
+            color += 0.5/RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
+        } else {
+            const vec4 obstacle_color = trace_obstacle(sample_pos, ray, delta);
+            if (obstacle_color.w == 1.0) {
+                return color + obstacle_color.xyz;
             } else {
-                const vec3 n = data.xyz; // recover surface normal
-                const float brightness = clamp(dot(n, ray), 0, 1);
-                color += vec3(max(0.3,brightness));
-                break;
+                color += 0.5/RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow});
             }
-        } else {
-            break;
         }
     }
 
@@ -226,13 +252,16 @@ void main(){
 
     result = vec4(trace(frag_pos, ray), 1.0);
 }
-""").substitute({
-    "size_x": lattice_x,
-    "size_y": lattice_y,
-    "size_z": lattice_z,
-    "inflow": inflow,
-    "max_ray_length": lattice_x
-}), GL_FRAGMENT_SHADER)
+""", lookup = TemplateLookup(directories = [
+    Path(__file__).parent
+])).render(
+    size_x = lattice_x,
+    size_y = lattice_y,
+    size_z = lattice_z,
+    inflow = inflow,
+    max_ray_length = max(lattice_x,lattice_y,lattice_z)**3,
+    sdf_source = channel
+), GL_FRAGMENT_SHADER)
 
 domain_program   = shaders.compileProgram(vertex_shader, fragment_shader)
 domain_projection_id = shaders.glGetUniformLocation(domain_program, 'projection')
@@ -247,17 +276,14 @@ lattice = Lattice(
     descriptor   = D3Q27,
     geometry     = Geometry(lattice_x, lattice_y, lattice_z),
     moments      = lbm.moments(optimize = True),
-    collide      = lbm.bgk(f_eq = lbm.equilibrium(), tau = relaxation_time),
+    collide      = lbm.bgk(f_eq = lbm.equilibrium(), tau = relaxation_time, optimize = True),
     boundary_src = boundary,
     opengl       = True
 )
 
-moments_texture = MomentsTexture(lattice)
+lattice.setup_channel_with_sdf_obstacle(channel)
 
-material_map = get_cavity_material_map(lattice.geometry)
-primitives   = list(map(lambda material: material[0], filter(lambda material: not callable(material[0]), material_map)))
-lattice.apply_material_map(material_map)
-lattice.sync_material()
+moments_texture = MomentsTexture(lattice)
 
 projection = Projection(distance = 2*lattice_x)
 rotation = Rotation([-0.5*lattice_x, -0.5*lattice_y, -0.5*lattice_z])
@@ -280,7 +306,7 @@ def on_display():
     shaders.glUseProgram(domain_program)
     glUniformMatrix4fv(domain_projection_id, 1, False, numpy.ascontiguousarray(projection.get()))
     glUniformMatrix4fv(domain_rotation_id,   1, False, numpy.ascontiguousarray(rotation.get()))
-    glLineWidth(2)
+    glLineWidth(3)
     glBegin(GL_LINES)
     for i, j in cube_edges:
         glVertex(cube_vertices[i])
@@ -310,6 +336,6 @@ glutDisplayFunc(on_display)
 glutReshapeFunc(lambda w, h: projection.update_ratio(w, h))
 glutMouseFunc(lambda *args: list(map(lambda m: m.on_mouse(*args), mouse_monitors)))
 glutMotionFunc(lambda *args: list(map(lambda m: m.on_mouse_move(*args), mouse_monitors)))
-glutTimerFunc(10, on_timer, 10)
+glutTimerFunc(30, on_timer, 30)
 
 glutMainLoop()
diff --git a/template/opengl.mako b/template/opengl.mako
index 181101c..47264e8 100644
--- a/template/opengl.mako
+++ b/template/opengl.mako
@@ -65,7 +65,6 @@ def neighbor_offset(c_i):
 %>
 
 __kernel void collect_gl_moments_to_texture(__global ${float_type}* f,
-                                            __global int* material,
 % if descriptor.d == 2:
                                             __write_only image2d_t moments)
 % elif descriptor.d == 3:
@@ -86,38 +85,17 @@ __kernel void collect_gl_moments_to_texture(__global ${float_type}* f,
 
     float4 data;
 
-    if (material[gid] != 5) {
 % if descriptor.d == 2:
-      data.x = ${ccode(moments_assignment[0].rhs)};
-      data.y = ${ccode(moments_assignment[1].rhs)};
-      data.z = ${ccode(moments_assignment[2].rhs)};
-      data.w = sqrt(data.y*data.y + data.z*data.z);
+    data.x = ${ccode(moments_assignment[0].rhs)};
+    data.y = ${ccode(moments_assignment[1].rhs)};
+    data.z = ${ccode(moments_assignment[2].rhs)};
+    data.w = sqrt(data.y*data.y + data.z*data.z);
 % elif descriptor.d == 3:
-      data.x = ${ccode(moments_assignment[0].rhs)};
-      data.y = ${ccode(moments_assignment[1].rhs)};
-      data.z = ${ccode(moments_assignment[2].rhs)};
-      data.w = ${ccode(moments_assignment[3].rhs)};
+    data.x = ${ccode(moments_assignment[0].rhs)};
+    data.y = ${ccode(moments_assignment[1].rhs)};
+    data.z = ${ccode(moments_assignment[2].rhs)};
+    data.w = ${ccode(moments_assignment[3].rhs)};
 % endif
-    } else {
-      const int material_west  = material[gid + ${neighbor_offset((-1,0,0))}];
-      const int material_east  = material[gid + ${neighbor_offset((1,0,0))}];
-      const int material_north = material[gid + ${neighbor_offset((0,1,0))}];
-      const int material_south = material[gid + ${neighbor_offset((0,-1,0))}];
-      const int material_up    = material[gid + ${neighbor_offset((0,0, 1))}];
-      const int material_down  = material[gid + ${neighbor_offset((0,0,-1))}];
-
-      // recover surface normal approximation using surrounding materials
-      float3 n;
-      if (material_west  != 5) { n.x =  1; }
-      if (material_east  != 5) { n.x = -1; }
-      if (material_north != 5) { n.y = -1; }
-      if (material_south != 5) { n.y =  1; }
-      if (material_up    != 5) { n.z = -1; }
-      if (material_down  != 5) { n.z =  1; }
-
-      data.xyz = normalize(n);
-      data.w   = 1.0; // signal impermeable material to raytracer
-    }
 
     write_imagef(moments, ${moments_cell()}, data);
 }
diff --git a/template/sdf.lib.glsl.mako b/template/sdf.lib.glsl.mako
index 515f959..c58cc7b 100644
--- a/template/sdf.lib.glsl.mako
+++ b/template/sdf.lib.glsl.mako
@@ -12,6 +12,10 @@ float box(vec3 v, vec3 b) {
   return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
 }
 
+float cylinder(vec3 v, float r, float h) {
+	return max(length(v.xy) - r, fabs(v.z) - 0.5*h);
+}
+
 vec3 flip_xy(vec3 v) {
 	return v3(v.y, v.x, v.z);
 }
@@ -43,7 +47,7 @@ float rounded(float a, float r) {
 	return a - r;
 }
 
-float sunify(float a, float b, float k) {
+float sadd(float a, float b, float k) {
 	float h = clamp(0.5 + 0.5 * (b - a) / k, 0.0, 1.0);
 	return mix(b, a, h) - k * h * (1 - h);
 }
@@ -62,7 +66,7 @@ float sub(float a, float b) {
 	return max(-b, a);
 }
 
-float unify(float a, float b) {
+float add(float a, float b) {
 	return min(a, b);
 }
 
diff --git a/utility/opengl.py b/utility/opengl.py
index ff88b17..ec8ab57 100644
--- a/utility/opengl.py
+++ b/utility/opengl.py
@@ -123,7 +123,6 @@ class MomentsTexture:
                 self.lattice.grid.size(),
                 self.lattice.layout,
                 self.lattice.memory.cl_pop_b,
-                self.lattice.memory.cl_material,
                 self.cl_gl_moments)
         else:
             self.program.collect_gl_moments_to_texture(
@@ -131,6 +130,5 @@ class MomentsTexture:
                 self.lattice.grid.size(),
                 self.lattice.layout,
                 self.lattice.memory.cl_pop_a,
-                self.lattice.memory.cl_material,
                 self.cl_gl_moments)