adding surface extraction pass

.clang-format

@@ -128,7 +128,7 @@ BreakFunctionDefinitionParameters: false
 BreakInheritanceList: BeforeColon
 BreakStringLiterals: true
 BreakTemplateDeclarations: MultiLine
-ColumnLimit:     140
+ColumnLimit:     160
 CommentPragmas:  '^ IWYU pragma:'
 CompactNamespaces: false
 ConstructorInitializerIndentWidth: 4

res/shaders/Extract.slang

@@ -0,0 +1,47 @@
+import FluidGridData;
+
+struct Parameter
+{
+    float iso;
+    FluidGridData<float> density;
+    RWStructuredBuffer<float3> vertices;
+    RWStructuredBuffer<uint> indices;
+    RWStructuredBuffer<uint> surfaceCount;
+    RWStructuredBuffer<uint> vertexCount;
+};
+ParameterBlock<Parameter> params;
+
+float getDensity(uint x, uint y, uint z)
+{
+    return params.density[x + gridSize.x * y + gridSize.x * gridSize.y * z];
+}
+
+[shader("compute")]
+[numthreads(32, 8, 1)]
+void main(uint3 dispatchThreadID : SV_DispatchThreadID)
+{
+    uint x = dispatchThreadID.x+1;
+    uint y = dispatchThreadID.y+1;
+    uint z = dispatchThreadID.z+1;
+    FluidGridData<float> density = params.density;
+    if(x >= gridSize.x-3 || y >= gridSize.y-3 || z >= gridSize.z-3) return;
+
+    float phi[8] = {
+        getDensity(x, y, z),
+        getDensity(x + 1, y, z),
+        getDensity(x, y + 1, z),
+        getDensity(x + 1, y + 1, z),
+        getDensity(x, y, z + 1),
+        getDensity(x + 1, y, z + 1),
+        getDensity(x, y + 1, z + 1),
+        getDensity(x + 1, y + 1, z + 1)
+    };
+
+    uint cubeIndex = 0;
+    for (uint i = 0; i < 8; ++i) {
+        if (phi[i] < params.iso) {
+            cubeIndex |= (1 << i);
+        }
+    }
+    if(cubeIndex == 0 || cubeIndex == 255) return;
+}

res/shaders/Render.slang

@@ -27,6 +27,6 @@ float4 fragmentMain(VertexOut input) : SV_Target
     FluidGridData<float> density = params.density;
     int x = int(input.uv.x * (gridSize.x - 2) + 1);
     int y = int(input.uv.y * (gridSize.y - 2) + 1);
-    float d = density[x, y, 0];
+    float d = density[x, y, 1];
     return float4(d, d, d, 1);
 }

src/CMakeLists.txt

@@ -1,9 +1,11 @@
 target_sources(FluidSimulation
   PRIVATE
-    SimulationComputePass.h
-    SimulationComputePass.cpp
     FluidRenderPass.h
     FluidRenderPass.cpp
     FluidSimulationView.h
     FluidSimulationView.cpp
-    main.cpp)
+    main.cpp
+    SimulationComputePass.h
+    SimulationComputePass.cpp
+    SurfaceExtractPass.h
+    SurfaceExtractPass.cpp)

src/Common.h

@@ -0,0 +1,8 @@
+#pragma once
+#include "MinimalEngine.h"
+
+constexpr static uint32 gridSize = 64;
+constexpr static float viscosity = 0.001f;
+constexpr static float diffusion = 0.001f;
+constexpr static float iso = 0.3f;
+constexpr static float dt = 0.1f;

src/SimulationComputePass.cpp

@@ -1,23 +1,27 @@
 #include "SimulationComputePass.h"
+#include "Graphics/Buffer.h"
 #include "Graphics/Command.h"
 #include "Graphics/Descriptor.h"
 #include "Graphics/Enums.h"
 #include "Graphics/Graphics.h"
 #include "Graphics/Initializer.h"
 #include "Graphics/Shader.h"
+#include "Common.h"
 
 using namespace Seele;
 
-#define SWAP(a, b)                                                                                                                         \
-  {                                                                                                                                        \
-    auto temp = std::move(a);                                                                                                              \
-    a = std::move(b);                                                                                                                      \
-    b = std::move(temp);                                                                                                                   \
+#define SWAP(a, b)                                                             \
+  {                                                                            \
+    auto temp = std::move(a);                                                  \
+    a = std::move(b);                                                          \
+    b = std::move(temp);                                                       \
   }
 
-SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : RenderPass(graphics) {
+SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics)
+    : RenderPass(graphics) {
   {
-    setBounds.pipelineLayout = graphics->createPipelineLayout("SetBoundsPipelineLayout");
+    setBounds.pipelineLayout =
+        graphics->createPipelineLayout("SetBoundsPipelineLayout");
     setBounds.descriptorLayout = graphics->createDescriptorLayout("params");
     setBounds.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
         .name = "b",
@@ -32,28 +36,34 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
     graphics->beginShaderCompilation(ShaderCompilationInfo{
         .name = "SetBound",
         .modules = {"SetBound"},
-        .entryPoints = {{"setBound", "SetBound"}, {"setBoundEdges", "SetBound"}, {"setBoundCorners", "SetBound"}},
+        .entryPoints = {{"setBound", "SetBound"},
+                        {"setBoundEdges", "SetBound"},
+                        {"setBoundCorners", "SetBound"}},
         .rootSignature = setBounds.pipelineLayout,
     });
     setBounds.pipelineLayout->create();
     setBounds.shader = graphics->createComputeShader({0});
     setBounds.shaderEdges = graphics->createComputeShader({1});
     setBounds.shaderCorners = graphics->createComputeShader({2});
-    setBounds.pipeline = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = setBounds.shader,
-        .pipelineLayout = setBounds.pipelineLayout,
-    });
-    setBounds.pipelineEdges = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = setBounds.shaderEdges,
-        .pipelineLayout = setBounds.pipelineLayout,
-    });
-    setBounds.pipelineCorners = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = setBounds.shaderCorners,
-        .pipelineLayout = setBounds.pipelineLayout,
-    });
+    setBounds.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = setBounds.shader,
+            .pipelineLayout = setBounds.pipelineLayout,
+        });
+    setBounds.pipelineEdges =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = setBounds.shaderEdges,
+            .pipelineLayout = setBounds.pipelineLayout,
+        });
+    setBounds.pipelineCorners =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = setBounds.shaderCorners,
+            .pipelineLayout = setBounds.pipelineLayout,
+        });
   }
   {
-    linearSolve.pipelineLayout = graphics->createPipelineLayout("LinearSolvePipelineLayout");
+    linearSolve.pipelineLayout =
+        graphics->createPipelineLayout("LinearSolvePipelineLayout");
     linearSolve.descriptorLayout = graphics->createDescriptorLayout("params");
     linearSolve.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
         .name = "a",
@@ -76,7 +86,8 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
         .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
     });
     linearSolve.descriptorLayout->create();
-    linearSolve.pipelineLayout->addDescriptorLayout(linearSolve.descriptorLayout);
+    linearSolve.pipelineLayout->addDescriptorLayout(
+        linearSolve.descriptorLayout);
     graphics->beginShaderCompilation(ShaderCompilationInfo{
         .name = "LinearSolve",
         .modules = {"LinearSolve"},
@@ -85,36 +96,45 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
     });
     linearSolve.pipelineLayout->create();
     linearSolve.shader = graphics->createComputeShader({0});
-    linearSolve.pipeline = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = linearSolve.shader,
-        .pipelineLayout = linearSolve.pipelineLayout,
-    });
+    linearSolve.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = linearSolve.shader,
+            .pipelineLayout = linearSolve.pipelineLayout,
+        });
   }
   {
-    computeDivergence.pipelineLayout = graphics->createPipelineLayout("ComputeDivergencePipelineLayout");
-    computeDivergence.descriptorLayout = graphics->createDescriptorLayout("params");
-    computeDivergence.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-        .name = "velocityX",
-        .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-    });
-    computeDivergence.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-        .name = "velocityY",
-        .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-    });
-    computeDivergence.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-        .name = "velocityZ",
-        .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-    });
-    computeDivergence.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-        .name = "divergence",
-        .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-    });
-    computeDivergence.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
-        .name = "pressure",
-        .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
-    });
+    computeDivergence.pipelineLayout =
+        graphics->createPipelineLayout("ComputeDivergencePipelineLayout");
+    computeDivergence.descriptorLayout =
+        graphics->createDescriptorLayout("params");
+    computeDivergence.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "velocityX",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    computeDivergence.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "velocityY",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    computeDivergence.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "velocityZ",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    computeDivergence.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "divergence",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
+    computeDivergence.descriptorLayout->addDescriptorBinding(
+        Gfx::DescriptorBinding{
+            .name = "pressure",
+            .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+        });
     computeDivergence.descriptorLayout->create();
-    computeDivergence.pipelineLayout->addDescriptorLayout(computeDivergence.descriptorLayout);
+    computeDivergence.pipelineLayout->addDescriptorLayout(
+        computeDivergence.descriptorLayout);
     graphics->beginShaderCompilation(ShaderCompilationInfo{
         .name = "Divergence",
         .modules = {"Divergence"},
@@ -124,13 +144,15 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
     computeDivergence.pipelineLayout->create();
 
     computeDivergence.shader = graphics->createComputeShader({0});
-    computeDivergence.pipeline = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = computeDivergence.shader,
-        .pipelineLayout = computeDivergence.pipelineLayout,
-    });
+    computeDivergence.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = computeDivergence.shader,
+            .pipelineLayout = computeDivergence.pipelineLayout,
+        });
   }
   {
-    project.pipelineLayout = graphics->createPipelineLayout("ProjectPipelineLayout");
+    project.pipelineLayout =
+        graphics->createPipelineLayout("ProjectPipelineLayout");
     project.descriptorLayout = graphics->createDescriptorLayout("params");
     project.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
         .name = "velocityX",
@@ -158,13 +180,15 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
     });
     project.pipelineLayout->create();
     project.shader = graphics->createComputeShader({0});
-    project.pipeline = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = project.shader,
-        .pipelineLayout = project.pipelineLayout,
-    });
+    project.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = project.shader,
+            .pipelineLayout = project.pipelineLayout,
+        });
   }
   {
-    advect.pipelineLayout = graphics->createPipelineLayout("AdvectPipelineLayout");
+    advect.pipelineLayout =
+        graphics->createPipelineLayout("AdvectPipelineLayout");
     advect.descriptorLayout = graphics->createDescriptorLayout("params");
     advect.descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
         .name = "density",
@@ -201,24 +225,44 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
     advect.pipelineLayout->create();
 
     advect.shader = graphics->createComputeShader({0});
-    advect.pipeline = graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
-        .computeShader = advect.shader,
-        .pipelineLayout = advect.pipelineLayout,
-    });
+    advect.pipeline =
+        graphics->createComputePipeline(Gfx::ComputePipelineCreateInfo{
+            .computeShader = advect.shader,
+            .pipelineLayout = advect.pipelineLayout,
+        });
   }
 
   uint64_t bufferSize = gridSize * gridSize * gridSize * sizeof(float);
-  Array<float> initialData(bufferSize / sizeof(float));
-  for(int i = 0; i < 20; ++i) {
-    int x = rand() % gridSize;
-    int y = rand() % gridSize;
-    initialData[x + y * gridSize] = 1.0f;
+  Array<float> initialDensity(bufferSize / sizeof(float));
+  Array<float> initialVelocityX(bufferSize / sizeof(float));
+  Array<float> initialVelocityY(bufferSize / sizeof(float));
+  Array<float> initialVelocityZ(bufferSize / sizeof(float));
+  Array<float> initialPhi(bufferSize / sizeof(float));
+  for (uint32 x = 1; x < gridSize - 2; ++x) {
+    for (uint32 y = 1; y < gridSize - 2; ++y) {
+        initialDensity[x + gridSize * y + gridSize * gridSize] = 1.0f;
+        initialPhi[x + gridSize * y + gridSize * gridSize] = iso - 1.0f;
+    }
   }
+
+  phiBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = bufferSize,
+          },
+      .name = "PhiBuffer",
+  });
+  phi0Buffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = bufferSize,
+          },
+      .name = "Phi0Buffer",
+  });
   densityBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
       .sourceData =
           {
               .size = bufferSize,
-              .data = (uint8*)initialData.data(),
           },
       .name = "DensityBuffer",
   });
@@ -227,6 +271,14 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
           {
               .size = bufferSize,
           },
+      .name = "ScratchBuffer",
+  });
+  density0Buffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = bufferSize,
+              .data = (uint8 *)initialDensity.data(),
+          },
       .name = "Density0Buffer",
   });
   velocityXNextBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
@@ -275,6 +327,7 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
       .sourceData =
           {
               .size = bufferSize,
+              .data = (uint8 *)initialVelocityX.data(),
           },
       .name = "VelocityX0Buffer",
   });
@@ -282,6 +335,7 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
       .sourceData =
           {
               .size = bufferSize,
+              .data = (uint8 *)initialVelocityY.data(),
           },
       .name = "VelocityY0Buffer",
   });
@@ -289,6 +343,7 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
       .sourceData =
           {
               .size = bufferSize,
+              .data = (uint8 *)initialVelocityZ.data(),
           },
       .name = "VelocityZ0Buffer",
   });
@@ -296,7 +351,9 @@ SimulationComputePass::SimulationComputePass(Seele::Gfx::PGraphics graphics) : R
 
 SimulationComputePass::~SimulationComputePass() {}
 
-void SimulationComputePass::beginFrame(const Seele::Component::Camera &camera, const Seele::Component::Transform &transform) {
+void SimulationComputePass::beginFrame(
+    const Seele::Component::Camera &camera,
+    const Seele::Component::Transform &transform) {
   setBounds.descriptorLayout->reset();
   linearSolve.descriptorLayout->reset();
   computeDivergence.descriptorLayout->reset();
@@ -305,12 +362,16 @@ void SimulationComputePass::beginFrame(const Seele::Component::Camera &camera, c
 }
 
 void SimulationComputePass::render() {
+  graphics->beginDebugRegion("Diffuse");
   for (uint iteration = 0; iteration < 4; ++iteration) {
     float a = dt * viscosity * (gridSize - 2) * (gridSize - 2) * (gridSize - 2);
     float c = 1 + 6 * a;
-    linearSolvePass(velocityXNextBuffer, velocityXBuffer, velocityX0Buffer, a, c);
-    linearSolvePass(velocityYNextBuffer, velocityYBuffer, velocityY0Buffer, a, c);
-    linearSolvePass(velocityZNextBuffer, velocityZBuffer, velocityZ0Buffer, a, c);
+    linearSolvePass(velocityXNextBuffer, velocityXBuffer, velocityX0Buffer, a,
+                    c);
+    linearSolvePass(velocityYNextBuffer, velocityYBuffer, velocityY0Buffer, a,
+                    c);
+    linearSolvePass(velocityZNextBuffer, velocityZBuffer, velocityZ0Buffer, a,
+                    c);
 
     setBoundsPass(velocityXNextBuffer, 1);
     setBoundsPass(velocityYNextBuffer, 2);
@@ -324,11 +385,14 @@ void SimulationComputePass::render() {
   SWAP(velocityX0Buffer, velocityXBuffer);
   SWAP(velocityY0Buffer, velocityYBuffer);
   SWAP(velocityZ0Buffer, velocityZBuffer);
+  graphics->endDebugRegion();
 
+  graphics->beginDebugRegion("Project");
   Gfx::PShaderBuffer divergence = velocityXBuffer;
   Gfx::PShaderBuffer pressureCurrent = velocityYBuffer;
   Gfx::PShaderBuffer pressureNext = velocityZBuffer;
-  divergencePass(divergence, pressureCurrent, velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer);
+  divergencePass(divergence, pressureCurrent, velocityX0Buffer,
+                 velocityY0Buffer, velocityZ0Buffer);
 
   for (int iteration = 0; iteration < 4; ++iteration) {
     linearSolvePass(pressureNext, pressureCurrent, divergence, 1, 6);
@@ -339,19 +403,26 @@ void SimulationComputePass::render() {
     SWAP(pressureCurrent, pressureNext);
   }
 
-  projectPass(velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer, pressureCurrent);
+  projectPass(velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer,
+              pressureCurrent);
 
   setBoundsPass(velocityX0Buffer, 1);
   setBoundsPass(velocityY0Buffer, 2);
   setBoundsPass(velocityZ0Buffer, 3);
+  graphics->endDebugRegion();
 
-  advectPass(velocityXBuffer, velocityX0Buffer, velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer, dt);
-  advectPass(velocityYBuffer, velocityY0Buffer, velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer, dt);
-  advectPass(velocityZBuffer, velocityZ0Buffer, velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer, dt);
+  graphics->beginDebugRegion("Advect");
+  advectPass(velocityXBuffer, velocityX0Buffer, velocityX0Buffer,
+             velocityY0Buffer, velocityZ0Buffer, dt);
+  advectPass(velocityYBuffer, velocityY0Buffer, velocityX0Buffer,
+             velocityY0Buffer, velocityZ0Buffer, dt);
+  advectPass(velocityZBuffer, velocityZ0Buffer, velocityX0Buffer,
+             velocityY0Buffer, velocityZ0Buffer, dt);
 
   setBoundsPass(velocityXBuffer, 1);
   setBoundsPass(velocityYBuffer, 2);
   setBoundsPass(velocityZBuffer, 3);
+  graphics->endDebugRegion();
 
   // swap velocity buffers
   SWAP(velocityX0Buffer, velocityXBuffer);
@@ -359,34 +430,62 @@ void SimulationComputePass::render() {
   SWAP(velocityZ0Buffer, velocityZBuffer);
 
   // diffuse density
+  graphics->beginDebugRegion("DiffuseDensity");
+  Gfx::PShaderBuffer densityCurrent = densityBuffer;
+  Gfx::PShaderBuffer densityNext = scratchBuffer;
+  Gfx::PShaderBuffer density0 = density0Buffer;
+  graphics->copyBuffer(density0, densityCurrent);
+  densityCurrent->pipelineBarrier(Gfx::SE_ACCESS_TRANSFER_WRITE_BIT,
+                                  Gfx::SE_PIPELINE_STAGE_TRANSFER_BIT,
+                                  Gfx::SE_ACCESS_SHADER_READ_BIT,
+                                  Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
   for (int iteration = 0; iteration < 4; ++iteration) {
     float a = dt * diffusion * (gridSize - 2) * (gridSize - 2) * (gridSize - 2);
     float c = 1 + 6 * a;
-    linearSolvePass(scratchBuffer, densityBuffer, densityBuffer, a, c);
-    setBoundsPass(scratchBuffer, 0);
+    linearSolvePass(densityNext, densityCurrent, density0, a, c);
+    setBoundsPass(densityNext, 0);
+    SWAP(densityCurrent, densityNext);
   }
-  advectPass(densityBuffer, scratchBuffer, velocityX0Buffer, velocityY0Buffer, velocityZ0Buffer, dt);
-  setBoundsPass(densityBuffer, 0);
+  graphics->endDebugRegion();
+  graphics->beginDebugRegion("AdvectDensity");
+  advectPass(density0Buffer, densityCurrent, velocityX0Buffer, velocityY0Buffer,
+             velocityZ0Buffer, dt);
+  setBoundsPass(density0Buffer, 0);
+  graphics->endDebugRegion();
+
+  graphics->beginDebugRegion("AdjectLevelSet");
+  advectPass(phiBuffer, phi0Buffer, velocityX0Buffer, velocityY0Buffer,
+             velocityZ0Buffer, dt);
+  setBoundsPass(phiBuffer, 0);
+  SWAP(phiBuffer, phi0Buffer);
+  // TODO: reinitialize level set every few frames
+  graphics->endDebugRegion();
 }
 
 void SimulationComputePass::endFrame() {}
 
-void SimulationComputePass::publishOutputs() { resources->registerBufferOutput("DENSITY", densityBuffer); }
+void SimulationComputePass::publishOutputs() {
+  resources->registerBufferOutput("DENSITY", densityBuffer);
+}
 
 void SimulationComputePass::createRenderPass() {}
 
 void SimulationComputePass::setBoundsPass(Gfx::PShaderBuffer grid, int b) {
-  Gfx::ODescriptorSet bounds = setBounds.descriptorLayout->allocateDescriptorSet();
+  Gfx::ODescriptorSet bounds =
+      setBounds.descriptorLayout->allocateDescriptorSet();
   bounds->updateConstants("b", 0, &b);
   bounds->updateBuffer("grid", 0, grid);
   bounds->writeChanges();
   Gfx::OComputeCommand boundsCommand = graphics->createComputeCommand();
   boundsCommand->bindPipeline(setBounds.pipeline);
   boundsCommand->bindDescriptor(bounds);
-  boundsCommand->dispatch(gridSize / setBounds.threadGroupSize.x, gridSize / setBounds.threadGroupSize.y, 1);
+  boundsCommand->dispatch(gridSize / setBounds.threadGroupSize.x,
+                          gridSize / setBounds.threadGroupSize.y, 1);
   graphics->executeCommands(std::move(boundsCommand));
 
-  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_READ_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_READ_BIT,
+                        Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                        Gfx::SE_ACCESS_SHADER_WRITE_BIT,
                         Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 
   Gfx::OComputeCommand edgesCommand = graphics->createComputeCommand();
@@ -395,7 +494,9 @@ void SimulationComputePass::setBoundsPass(Gfx::PShaderBuffer grid, int b) {
   edgesCommand->dispatch(gridSize / setBounds.threadGroupSize.x, 1, 1);
   graphics->executeCommands(std::move(edgesCommand));
 
-  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_READ_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_READ_BIT,
+                        Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                        Gfx::SE_ACCESS_SHADER_WRITE_BIT,
                         Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 
   Gfx::OComputeCommand cornersCommand = graphics->createComputeCommand();
@@ -404,33 +505,43 @@ void SimulationComputePass::setBoundsPass(Gfx::PShaderBuffer grid, int b) {
   cornersCommand->dispatch(1, 1, 1);
   graphics->executeCommands(std::move(cornersCommand));
 
-  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  grid->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                        Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                        Gfx::SE_ACCESS_SHADER_READ_BIT,
                         Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }
 
-void SimulationComputePass::linearSolvePass(Gfx::PShaderBuffer current, Gfx::PShaderBuffer next, Gfx::PShaderBuffer grid0, float a,
+void SimulationComputePass::linearSolvePass(Gfx::PShaderBuffer next,
+                                            Gfx::PShaderBuffer current,
+                                            Gfx::PShaderBuffer grid0, float a,
                                             float c) {
-  Gfx::ODescriptorSet solve = linearSolve.descriptorLayout->allocateDescriptorSet();
-  solve->updateBuffer("grid0", 0, grid0);
-  solve->updateBuffer("current", 0, current);
+  Gfx::ODescriptorSet solve =
+      linearSolve.descriptorLayout->allocateDescriptorSet();
   solve->updateBuffer("next", 0, next);
+  solve->updateBuffer("current", 0, current);
+  solve->updateBuffer("grid0", 0, grid0);
   solve->updateConstants("a", 0, &a);
   solve->updateConstants("c", 0, &c);
   solve->writeChanges();
   Gfx::OComputeCommand solveCommand = graphics->createComputeCommand();
   solveCommand->bindPipeline(linearSolve.pipeline);
   solveCommand->bindDescriptor(solve);
-  solveCommand->dispatch(gridSize / linearSolve.threadGroupSize.x, gridSize / linearSolve.threadGroupSize.y,
+  solveCommand->dispatch(gridSize / linearSolve.threadGroupSize.x,
+                         gridSize / linearSolve.threadGroupSize.y,
                          gridSize / linearSolve.threadGroupSize.z);
   graphics->executeCommands(std::move(solveCommand));
 
-  next->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  next->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                        Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                        Gfx::SE_ACCESS_SHADER_READ_BIT,
                         Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }
-void SimulationComputePass::divergencePass(Seele::Gfx::PShaderBuffer divergence, Seele::Gfx::PShaderBuffer pressure,
-                                           Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY,
-                                           Seele::Gfx::PShaderBuffer velocityZ) {
-  Gfx::ODescriptorSet divergenceSet = computeDivergence.descriptorLayout->allocateDescriptorSet();
+void SimulationComputePass::divergencePass(
+    Seele::Gfx::PShaderBuffer divergence, Seele::Gfx::PShaderBuffer pressure,
+    Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY,
+    Seele::Gfx::PShaderBuffer velocityZ) {
+  Gfx::ODescriptorSet divergenceSet =
+      computeDivergence.descriptorLayout->allocateDescriptorSet();
   divergenceSet->updateBuffer("velocityX", 0, velocityX);
   divergenceSet->updateBuffer("velocityY", 0, velocityY);
   divergenceSet->updateBuffer("velocityZ", 0, velocityZ);
@@ -440,19 +551,27 @@ void SimulationComputePass::divergencePass(Seele::Gfx::PShaderBuffer divergence,
   Gfx::OComputeCommand divergenceCommand = graphics->createComputeCommand();
   divergenceCommand->bindPipeline(computeDivergence.pipeline);
   divergenceCommand->bindDescriptor(divergenceSet);
-  divergenceCommand->dispatch(gridSize / computeDivergence.threadGroupSize.x, gridSize / computeDivergence.threadGroupSize.y,
+  divergenceCommand->dispatch(gridSize / computeDivergence.threadGroupSize.x,
+                              gridSize / computeDivergence.threadGroupSize.y,
                               gridSize / computeDivergence.threadGroupSize.z);
   graphics->executeCommands(std::move(divergenceCommand));
   // divergence is now in velocityXBuffer, pressure is in velocityYBuffer
-  divergence->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  divergence->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                              Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                              Gfx::SE_ACCESS_SHADER_READ_BIT,
                               Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
-  pressure->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  pressure->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                            Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                            Gfx::SE_ACCESS_SHADER_READ_BIT,
                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }
 
-void SimulationComputePass::projectPass(Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY,
-                                        Seele::Gfx::PShaderBuffer velocityZ, Seele::Gfx::PShaderBuffer pressure) {
-  Gfx::ODescriptorSet projectSet = project.descriptorLayout->allocateDescriptorSet();
+void SimulationComputePass::projectPass(Seele::Gfx::PShaderBuffer velocityX,
+                                        Seele::Gfx::PShaderBuffer velocityY,
+                                        Seele::Gfx::PShaderBuffer velocityZ,
+                                        Seele::Gfx::PShaderBuffer pressure) {
+  Gfx::ODescriptorSet projectSet =
+      project.descriptorLayout->allocateDescriptorSet();
   projectSet->updateBuffer("velocityX", 0, velocityX);
   projectSet->updateBuffer("velocityY", 0, velocityY);
   projectSet->updateBuffer("velocityZ", 0, velocityZ);
@@ -461,21 +580,32 @@ void SimulationComputePass::projectPass(Seele::Gfx::PShaderBuffer velocityX, See
   Gfx::OComputeCommand projectCommand = graphics->createComputeCommand();
   projectCommand->bindPipeline(project.pipeline);
   projectCommand->bindDescriptor(projectSet);
-  projectCommand->dispatch(gridSize / project.threadGroupSize.x, gridSize / project.threadGroupSize.y,
+  projectCommand->dispatch(gridSize / project.threadGroupSize.x,
+                           gridSize / project.threadGroupSize.y,
                            gridSize / project.threadGroupSize.z);
   graphics->executeCommands(std::move(projectCommand));
 
-  velocityX->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  velocityX->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                             Gfx::SE_ACCESS_SHADER_READ_BIT,
                              Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
-  velocityY->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  velocityY->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                             Gfx::SE_ACCESS_SHADER_READ_BIT,
                              Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
-  velocityZ->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  velocityZ->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                             Gfx::SE_ACCESS_SHADER_READ_BIT,
                              Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }
 
-void SimulationComputePass::advectPass(Gfx::PShaderBuffer quantity, Gfx::PShaderBuffer quantity0, Gfx::PShaderBuffer velocityX,
-                                       Gfx::PShaderBuffer velocityY, Gfx::PShaderBuffer velocityZ, float dt) {
-  Gfx::ODescriptorSet advectSet = this->advect.descriptorLayout->allocateDescriptorSet();
+void SimulationComputePass::advectPass(Gfx::PShaderBuffer quantity,
+                                       Gfx::PShaderBuffer quantity0,
+                                       Gfx::PShaderBuffer velocityX,
+                                       Gfx::PShaderBuffer velocityY,
+                                       Gfx::PShaderBuffer velocityZ, float dt) {
+  Gfx::ODescriptorSet advectSet =
+      this->advect.descriptorLayout->allocateDescriptorSet();
   advectSet->updateBuffer("density", 0, quantity);
   advectSet->updateBuffer("density0", 0, quantity0);
   advectSet->updateBuffer("velocityX", 0, velocityX);
@@ -486,9 +616,13 @@ void SimulationComputePass::advectPass(Gfx::PShaderBuffer quantity, Gfx::PShader
   Gfx::OComputeCommand advectCommand = graphics->createComputeCommand();
   advectCommand->bindPipeline(advect.pipeline);
   advectCommand->bindDescriptor(advectSet);
-  advectCommand->dispatch(gridSize / advect.threadGroupSize.x, gridSize / advect.threadGroupSize.y, gridSize / advect.threadGroupSize.z);
+  advectCommand->dispatch(gridSize / advect.threadGroupSize.x,
+                          gridSize / advect.threadGroupSize.y,
+                          gridSize / advect.threadGroupSize.z);
   graphics->executeCommands(std::move(advectCommand));
 
-  quantity->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT, Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT, Gfx::SE_ACCESS_SHADER_READ_BIT,
+  quantity->pipelineBarrier(Gfx::SE_ACCESS_SHADER_WRITE_BIT,
+                            Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                            Gfx::SE_ACCESS_SHADER_READ_BIT,
                             Gfx::SE_PIPELINE_STAGE_COMPUTE_SHADER_BIT);
 }

src/SimulationComputePass.h

@@ -22,10 +22,6 @@ private:
     void divergencePass(Seele::Gfx::PShaderBuffer divergence, Seele::Gfx::PShaderBuffer pressure, Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ);
     void projectPass(Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ, Seele::Gfx::PShaderBuffer pressure);
     void advectPass(Seele::Gfx::PShaderBuffer quantity, Seele::Gfx::PShaderBuffer quantity0, Seele::Gfx::PShaderBuffer velocityX, Seele::Gfx::PShaderBuffer velocityY, Seele::Gfx::PShaderBuffer velocityZ, float dt);
-    float viscosity = 0.001f;
-    float diffusion = 0.0001f;
-    float dt = 0.1f;
-    constexpr static uint gridSize = 64;
     struct SetBounds
     {
         Seele::Gfx::OPipelineLayout pipelineLayout;
@@ -70,8 +66,11 @@ private:
         Seele::Gfx::PComputePipeline pipeline;
         constexpr static Seele::UVector threadGroupSize = {32, 8, 1};
     } advect;
+    Seele::Gfx::OShaderBuffer phiBuffer;
+    Seele::Gfx::OShaderBuffer phi0Buffer;
     Seele::Gfx::OShaderBuffer densityBuffer;
     Seele::Gfx::OShaderBuffer scratchBuffer;
+    Seele::Gfx::OShaderBuffer density0Buffer;
     Seele::Gfx::OShaderBuffer velocityXNextBuffer;
     Seele::Gfx::OShaderBuffer velocityYNextBuffer;
     Seele::Gfx::OShaderBuffer velocityZNextBuffer;

src/SurfaceExtractPass.cpp

@@ -0,0 +1,122 @@
+#include "SurfaceExtractPass.h"
+#include "Common.h"
+#include "Graphics/Descriptor.h"
+#include "Graphics/Enums.h"
+#include "Graphics/Graphics.h"
+#include "Graphics/Shader.h"
+
+using namespace Seele;
+
+SurfaceExtractPass::SurfaceExtractPass(Seele::Gfx::PGraphics graphics)
+    : RenderPass(graphics) {
+  pipelineLayout =
+      graphics->createPipelineLayout("SurfaceExtractPipelineLayout");
+  descriptorLayout = graphics->createDescriptorLayout("params");
+  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
+      .name = "density",
+      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+  });
+  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
+      .name = "vertexBuffer",
+      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+  });
+  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
+      .name = "indexBuffer",
+      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+  });
+  descriptorLayout->addDescriptorBinding(Gfx::DescriptorBinding{
+      .name = "surfaceCounts",
+      .descriptorType = Gfx::SE_DESCRIPTOR_TYPE_STORAGE_BUFFER,
+  });
+  descriptorLayout->create();
+  vertexBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = 1024 * 1024, // TODO: make this dynamic
+          },
+      .name = "SurfaceExtractVertexBuffer",
+  });
+  indexBuffer = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = 1024 * 1024, // TODO: make this dynamic
+          },
+      .name = "SurfaceExtractIndexBuffer",
+  });
+  surfaceCounts = graphics->createShaderBuffer(ShaderBufferCreateInfo{
+      .sourceData =
+          {
+              .size = sizeof(uint32) * 3, // TODO: make this dynamic
+          },
+  });
+}
+
+SurfaceExtractPass::~SurfaceExtractPass() {}
+
+void SurfaceExtractPass::beginFrame(
+    const Seele::Component::Camera &camera,
+    const Seele::Component::Transform &transform) {}
+
+void SurfaceExtractPass::render() {
+  Array<float> densityData(densityBuffer->getNumElements());
+  densityBuffer->readContents(0, densityData.size() * sizeof(float),
+                              densityData.data());
+  for (int x = 1; x < gridSize - 3; ++x) {
+    for (int y = 1; y < gridSize - 3; ++y) {
+      for (int z = 1; z < gridSize - 3; ++z) {
+        auto getDensity = [&](int x, int y, int z) {
+          return densityData[x + gridSize * (y + gridSize * z)];
+        };
+        float cube[8] = {
+            getDensity(x, y, z),
+            getDensity(x + 1, y, z),
+            getDensity(x, y + 1, z),
+            getDensity(x + 1, y + 1, z),
+            getDensity(x, y, z + 1),
+            getDensity(x + 1, y, z + 1),
+            getDensity(x, y + 1, z + 1),
+            getDensity(x + 1, y + 1, z + 1)
+        };
+
+      }
+    }
+  }
+}
+
+void SurfaceExtractPass::endFrame() {}
+
+void SurfaceExtractPass::publishOutputs() {
+  resources->registerBufferOutput("SURFACE_VERTICES", vertexBuffer);
+  resources->registerBufferOutput("SURFACE_INDICES", indexBuffer);
+  resources->registerBufferOutput("SURFACE_COUNTS", surfaceCounts);
+}
+
+void SurfaceExtractPass::createRenderPass() {
+  densityBuffer = resources->requestBuffer("DENSITY");
+}
+
+void SurfaceExtractPass::computeEdgeTable() {
+    for(uint i = 0; i < 256; ++i) {
+        uint32 edgeFlags = 0;
+        for(int j = 0; j < 8; ++j) {
+            if(i & (1 << j)) {
+                edgeFlags |= (1 << j);
+            }
+        }
+        edgeTable[i] = edgeFlags;
+    }
+}
+
+void SurfaceExtractPass::computeTriTable() {
+    for(uint i = 0; i < 256; ++i) {
+        uint32 edgeFlags = edgeTable[i];
+        int triIndex = 0;
+        for(int j = 0; j < 16; ++j) {
+            if(edgeFlags & (1 << j)) {
+                triTable[i][triIndex++] = j;
+            } else {
+                triTable[i][triIndex] = -1;
+            }
+        }
+    }
+}

src/SurfaceExtractPass.h

@@ -0,0 +1,36 @@
+#pragma once
+#include "Graphics/Graphics.h"
+#include "Graphics/RenderPass/RenderPass.h"
+#include "Graphics/Resources.h"
+#include <cstdint>
+
+class SurfaceExtractPass : public Seele::RenderPass {
+public:
+  SurfaceExtractPass(Seele::Gfx::PGraphics graphics);
+  virtual ~SurfaceExtractPass() override;
+  virtual void
+  beginFrame(const Seele::Component::Camera &camera,
+             const Seele::Component::Transform &transform) override;
+  virtual void render() override;
+  virtual void endFrame() override;
+  virtual void publishOutputs() override;
+  virtual void createRenderPass() override;
+
+private:
+  void computeEdgeTable();
+  void computeTriTable();
+  Seele::Gfx::OPipelineLayout pipelineLayout;
+  Seele::Gfx::ODescriptorLayout descriptorLayout;
+  Seele::Gfx::ODescriptorSet descriptorSet;
+  Seele::Gfx::OComputeShader computeShader;
+  Seele::Gfx::PComputePipeline pipeline;
+
+  Seele::Gfx::OShaderBuffer vertexBuffer;
+  Seele::Gfx::OShaderBuffer indexBuffer;
+  Seele::Gfx::OShaderBuffer surfaceCounts;
+  Seele::Gfx::PShaderBuffer densityBuffer;
+
+  static uint32_t edgeTable[256];
+  static uint32_t triTable[256][16];
+};
+DEFINE_REF(SurfaceExtractPass)