1 files changed, 48 insertions, 0 deletions
diff --git a/tangle/LLBM/kernel/collect_streamlines.h b/tangle/LLBM/kernel/collect_streamlines.h
new file mode 100644
index 0000000..c787713
--- /dev/null
+++ b/tangle/LLBM/kernel/collect_streamlines.h
@@ -0,0 +1,48 @@
+#pragma once
+
+__device__ void writeDot(cudaSurfaceObject_t texture, int x, int y, uchar4 pixel) {
+  surf2Dwrite(pixel, texture, (x-1)*sizeof(uchar4), (y  ));
+  surf2Dwrite(pixel, texture, (x-1)*sizeof(uchar4), (y+1));
+  surf2Dwrite(pixel, texture, (x  )*sizeof(uchar4), (y  ));
+  surf2Dwrite(pixel, texture, (x  )*sizeof(uchar4), (y+1));
+} 
+
+__device__ float2 readVelocity(cudaSurfaceObject_t u_x, cudaSurfaceObject_t u_y, float2 p) {
+  return make_float2(tex2D<float>(u_x, p.x, p.y), tex2D<float>(u_y, p.x, p.y));
+}
+
+template <typename DESCRIPTOR, typename T>
+__global__ void renderStreamlinesToTexture(
+    descriptor::Cuboid<DESCRIPTOR> cuboid
+  , cudaSurfaceObject_t u_x
+  , cudaSurfaceObject_t u_y
+  , T* origins
+  , unsigned nSteps
+  , float charU
+  , cudaSurfaceObject_t texture
+)  {
+  const int iOrigin = threadIdx.x + blockIdx.x * blockDim.x;
+
+  float h = 1;
+  float2 p = make_float2(origins[2*iOrigin+0], origins[2*iOrigin+1]);
+
+  for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
+    float2 u = make_float2(tex2D<float>(u_x, p.x, p.y), tex2D<float>(u_y, p.x, p.y));
+
+    float2 k1 = readVelocity(u_x, u_y, p) / charU;
+    float2 k2 = readVelocity(u_x, u_y, p + 0.5f*h * k1) / charU;
+    float2 k3 = readVelocity(u_x, u_y, p + 0.5f*h * k2) / charU;
+    float2 k4 = readVelocity(u_x, u_y, p +      h * k3) / charU;
+
+    p += h * (1.f/6.f*k1 + 1.f/3.f*k2 + 1.f/3.f*k3 + 1.f/6.f*k4);
+
+    const int screenX = std::nearbyint(p.x);
+    const int screenY = cuboid.nY-1 - std::nearbyint(p.y);
+
+    if (screenX < 1 || screenY < 1 || screenX > cuboid.nX-2 || screenY > cuboid.nY-2) {
+      break;
+    }
+
+    writeDot(texture, screenX, screenY, {255,255,255,255});
+  }
+}