From 342d554e2da1f80aa93d5d28d2d87573db87bbb0 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Fri, 27 Dec 2019 23:45:51 +0100 Subject: Add SDF-based grid fin example --- ...sdf_grid_fin_volumetric_rendering_gl_interop.py | 352 +++++++++++++++++++++ 1 file changed, 352 insertions(+) create mode 100644 channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py (limited to 'channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py') diff --git a/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py b/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py new file mode 100644 index 0000000..166c36d --- /dev/null +++ b/channel_3d_sdf_grid_fin_volumetric_rendering_gl_interop.py @@ -0,0 +1,352 @@ +import numpy + +from mako.template import Template +from mako.lookup import TemplateLookup + +from pathlib import Path + +from simulation import Lattice, Geometry +from utility.particles import Particles +from symbolic.generator import LBM + +import symbolic.D3Q19 as D3Q19 + +from OpenGL.GL import * +from OpenGL.GLUT import * + +from OpenGL.GL import shaders + +from geometry.box import Box + +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 + +updates_per_frame = 5 + +inflow = 0.075 +relaxation_time = 0.51 + +lbm = LBM(D3Q19) + +boundary = Template(""" + 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_1 = 0.0; + u_2 = 0.0; + } + if ( m == 4 ) { + rho = 1.0; + } +""").render( + inflow = inflow +) + +grid_fin = """ +v = rotate_z(translate(v, v3(center.x/2.5, center.y, center.z)), -0.5); +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); +""" + +def glut_window(fullscreen = False): + glutInit(sys.argv) + glutSetOption(GLUT_MULTISAMPLE, 8) + glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_MULTISAMPLE) + + if fullscreen: + window = glutEnterGameMode() + else: + glutInitWindowSize(800, 500) + glutInitWindowPosition(0, 0) + window = glutCreateWindow("LBM") + + return window + +window = glut_window(fullscreen = False) + +vertex_shader = shaders.compileShader(""" +#version 430 + +layout (location=0) in vec4 vertex; + +uniform mat4 projection; +uniform mat4 rotation; + +void main() { + gl_Position = projection * rotation * vertex; +}""", GL_VERTEX_SHADER) + +fragment_shader = shaders.compileShader(""" +#version 430 + +in vec3 color; + +void main(){ + gl_FragColor = vec4(vec3(0.5), 1.0); +}""", GL_FRAGMENT_SHADER) + +raycast_vertex_shader = shaders.compileShader(""" +#version 430 + +layout (location=0) in vec4 vertex; + out vec3 frag_pos; + +uniform mat4 projection; +uniform mat4 rotation; + +void main() { + gl_Position = projection * rotation * vertex; + frag_pos = vertex.xyz; +}""", GL_VERTEX_SHADER) + +raycast_fragment_shader = shaders.compileShader(Template(""" +#version 430 + +in vec3 frag_pos; + +uniform vec4 camera_pos; + +uniform sampler3D moments; + +out vec4 result; + +vec3 unit(vec3 v) { + 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); +} + +<%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) { + return mix( + vec3(0.0, 0.0, 0.0), + vec3(1.0, 0.0, 0.0), + x + ); +} + +float distanceToLattice(vec3 v) { + return box(v, vec3(${size_x},${size_y},${size_z})); +} + +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; + + 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; + } else { + sample_pos += dist*ray; + } + } + + if (hit) { + return vec4(color, 1.0); + } else { + return vec4(color, 0.0); + } +} + +vec3 trace(vec3 pos, vec3 ray) { + const float delta = maxRayLength(pos, ray) / RAYMARCH_STEPS; + + 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}); + } else { + vec4 obstacle_color = trace_obstacle(pos, ray); + if (obstacle_color.w == 1.0) { + color += obstacle_color.xyz; + break; + } else { + color += 1./RAYMARCH_STEPS * palette(length(data.yzw) / ${inflow}); + } + } + } + + return color; +} + +void main(){ + const vec3 ray = normalize(frag_pos - camera_pos.xyz); + + result = vec4(trace(frag_pos, ray), 1.0); +} +""", 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 = grid_fin +), GL_FRAGMENT_SHADER) + +domain_program = shaders.compileProgram(vertex_shader, fragment_shader) +domain_projection_id = shaders.glGetUniformLocation(domain_program, 'projection') +domain_rotation_id = shaders.glGetUniformLocation(domain_program, 'rotation') + +raycast_program = shaders.compileProgram(raycast_vertex_shader, raycast_fragment_shader) +raycast_projection_id = shaders.glGetUniformLocation(raycast_program, 'projection') +raycast_rotation_id = shaders.glGetUniformLocation(raycast_program, 'rotation') +raycast_camera_pos_id = shaders.glGetUniformLocation(raycast_program, 'camera_pos') + +lattice = Lattice( + descriptor = D3Q19, + geometry = Geometry(lattice_x, lattice_y, lattice_z), + moments = lbm.moments(optimize = True), + collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = relaxation_time, optimize = True), + boundary_src = boundary, + opengl = True +) + +lattice.setup_channel_with_sdf_obstacle(grid_fin) + +moments_texture = MomentsTexture(lattice) + +projection = Projection(distance = 2*lattice_x) +rotation = Rotation([-0.5*lattice_x, -0.5*lattice_y, -0.5*lattice_z]) + +cube_vertices, cube_edges = lattice.geometry.wireframe() + +def on_display(): + for i in range(0,updates_per_frame): + lattice.evolve() + + moments_texture.collect() + + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + + glEnable(GL_DEPTH_TEST) + glDepthFunc(GL_LESS) + + camera_pos = numpy.matmul([0,-projection.distance,0,1], rotation.get_inverse()) + + 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(3) + glBegin(GL_LINES) + for i, j in cube_edges: + glVertex(cube_vertices[i]) + glVertex(cube_vertices[j]) + glEnd() + + shaders.glUseProgram(raycast_program) + glUniformMatrix4fv(raycast_projection_id, 1, False, numpy.ascontiguousarray(projection.get())) + glUniformMatrix4fv(raycast_rotation_id, 1, False, numpy.ascontiguousarray(rotation.get())) + glUniform4fv(raycast_camera_pos_id, 1, camera_pos) + moments_texture.bind() + Box(0,lattice.geometry.size_x,0,lattice.geometry.size_y,0,lattice.geometry.size_z).draw() + + glutSwapBuffers() + +mouse_monitors = [ + MouseDragMonitor(GLUT_LEFT_BUTTON, lambda dx, dy: rotation.update(0.005*dy, 0.005*dx)), + MouseDragMonitor(GLUT_RIGHT_BUTTON, lambda dx, dy: rotation.shift(0.25*dx, 0.25*dy)), + MouseScrollMonitor(lambda zoom: projection.update_distance(5*zoom)) +] + +def on_timer(t): + glutTimerFunc(t, on_timer, t) + glutPostRedisplay() + +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(30, on_timer, 30) + +glutMainLoop() -- cgit v1.2.3