Using MLS MPM to combine shaders and reduce sources of errors

Signed-off-by: Tobias Frisch <tfrisch@uni-koblenz.de>

projects/wobble_bobble/shaders/particle.inc

@@ -14,6 +14,7 @@ struct ParticleMinimal {
 struct Particle {
     ParticleMinimal minimal;
     mat4 deformation;
+	mat4 mls;
 };
 
 const float PI = 3.1415926535897932384626433832795;

projects/wobble_bobble/shaders/transform_particles_to_grid.comp

@@ -27,6 +27,12 @@ void main()	{
 
         if (localOffset < particles.length()) {
             shared_particles[gl_LocalInvocationIndex] = particles[localOffset].minimal;
+
+            vec3 distance = (position - shared_particles[gl_LocalInvocationIndex].position);
+
+            shared_particles[gl_LocalInvocationIndex].velocity += (
+                mat3(particles[localOffset].mls) * distance
+            );
         } else {
             shared_particles[gl_LocalInvocationIndex].position = vec3(0.0f);
             shared_particles[gl_LocalInvocationIndex].size = 0.0f;

projects/wobble_bobble/shaders/update_grid_forces.comp

@@ -14,6 +14,8 @@ layout(set=0, binding=2, std430) readonly buffer particleBuffer {
 layout( push_constant ) uniform constants {
     float lame1;
     float lame2;
+    float alpha;
+    float beta;
     float t;
     float dt;
 };
@@ -41,63 +43,11 @@ void main()	{
     vec3 velocity = gridSample.xyz;
     float mass = gridSample.w;
 
-    vec3 force = vec3(0.0f);
-    uint offset = 0;
-
     barrier();
     memoryBarrierBuffer();
 
-    for (offset = 0; offset < particles.length(); offset += SHARED_PARTICLES_BATCH_SIZE) {
-        uint localOffset = offset + gl_LocalInvocationIndex;
-
-        if (localOffset < particles.length()) {
-            shared_particles[gl_LocalInvocationIndex] = particles[localOffset];
-        } else {
-            shared_particles[gl_LocalInvocationIndex].minimal.position = vec3(0.0f);
-            shared_particles[gl_LocalInvocationIndex].minimal.size = 0.0f;
-            shared_particles[gl_LocalInvocationIndex].minimal.velocity = vec3(0.0f);
-            shared_particles[gl_LocalInvocationIndex].minimal.mass = 0.0f;
-            shared_particles[gl_LocalInvocationIndex].deformation = mat4(0.0f);
-        }
-
-        barrier();
-        memoryBarrierShared();
-
-        for (uint i = 0; i < SHARED_PARTICLES_BATCH_SIZE; i++) {
-            float volume = sphere_volume(shared_particles[i].minimal.size);
-            mat3 F = mat3(shared_particles[i].deformation);
-
-            float J = determinant(F);
-
-            if ((J > 0.0f) && (volume > 0.0f)) {
-                mat3 F_invT = inverse(transpose(F));
-
-                vec3 v = shared_particles[i].minimal.velocity;
-                float m = shared_particles[i].minimal.mass;
-
-                mat3 delta = lame2 * (F - F_invT) + lame1 * log(J) * F_invT;
-
-                vec3 weight_gradient = voxel_particle_grad_weight(
-                    position,
-                    shared_particles[i].minimal
-                );
-
-                force -= (
-                    volume *
-                    delta * transpose(F) *
-                    weight_gradient
-                );
-            }
-        }
-
-        barrier();
-        memoryBarrierShared();
-    }
-
-    if (dt > 0.0f) {
-        if (mass > 0.0f) {
-            velocity += force * dt / mass;
-        }
+    if (mass > 0.0f) {
+        velocity += vec3(0.0f, -9.81f * dt * 0.0f, 0.0f);
     }
 
     bvec3 lowerID = lessThanEqual(gl_GlobalInvocationID, ivec3(0));

projects/wobble_bobble/shaders/update_particle_deformation.comp

-#version 450
-#extension GL_GOOGLE_include_directive : enable
-
-layout(local_size_x = 64, local_size_y = 1, local_size_z = 1) in;
-
-#include "particle.inc"
-
-layout(set=0, binding=0, std430) restrict buffer particleBuffer {
-    Particle particles [];
-};
-
-layout(set=0, binding=1) uniform texture3D gridImage;
-layout(set=0, binding=2) uniform sampler gridSampler;
-
-layout( push_constant ) uniform constants {
-    float K;
-    float E;
-    float t;
-    float dt;
-};
-
-void main()	{
-    if (gl_GlobalInvocationID.x < particles.length()) {
-        ParticleMinimal minimal = particles[gl_GlobalInvocationID.x].minimal;
-        mat3 deformation = mat3(particles[gl_GlobalInvocationID.x].deformation);
-
-        ivec3 gridResolution = textureSize(sampler3D(gridImage, gridSampler), 0);
-        ivec3 gridWindow = ivec3(minimal.size * 2.0f * gridResolution);
-
-        mat3 velocity_gradient = mat3(0.0f);
-
-        int i, j, k;
-
-        for (i = -gridWindow.x; i <= gridWindow.x; i++) {
-            for (j = -gridWindow.y; j <= gridWindow.y; j++) {
-                for (k = -gridWindow.z; k <= gridWindow.z; k++) {
-                    vec3 offset = vec3(i, j, k) / gridResolution;
-                    vec3 voxel = minimal.position + offset;
-
-                    vec4 gridSample = texture(sampler3D(gridImage, gridSampler), voxel);
-
-                    velocity_gradient += outerProduct(
-                        gridSample.xyz,
-                        voxel_particle_grad_weight(voxel, minimal)
-                    );
-                }
-            }
-        }
-
-        barrier();
-        memoryBarrierBuffer();
-
-        mat3 deformation_elastic = deformation;
-
-        if (dt > 0.0f) {
-            deformation_elastic = (mat3(1.0f) + dt * velocity_gradient) * deformation;
-        }
-
-        barrier();
-        memoryBarrierBuffer();
-
-        particles[gl_GlobalInvocationID.x].deformation = mat4(deformation_elastic);
-    }
-}
\ No newline at end of file

projects/wobble_bobble/shaders/update_particle_positions.comp

-#version 450
-#extension GL_GOOGLE_include_directive : enable
-
-layout(local_size_x = 64, local_size_y = 1, local_size_z = 1) in;
-
-#include "particle.inc"
-
-layout(set=0, binding=0, std430) restrict buffer particleBuffer {
-    Particle particles [];
-};
-
-layout( push_constant ) uniform constants {
-    float K;
-    float E;
-    float t;
-    float dt;
-};
-
-void main()	{
-    memoryBarrierBuffer();
-
-    if (gl_GlobalInvocationID.x < particles.length()) {
-        vec3 position = particles[gl_GlobalInvocationID.x].minimal.position;
-        float size = particles[gl_GlobalInvocationID.x].minimal.size;
-        vec3 velocity = particles[gl_GlobalInvocationID.x].minimal.velocity;
-
-        if (dt > 0.0f) {
-            position = position + velocity * dt;
-        }
-
-        for (uint i = 0; i < 3; i++) {
-            if (position[i] - size < 0.0f) {
-                position[i] = -position[i] + 2.0f * size;
-                velocity[i] *= -1.0f;
-            } else
-            if (position[i] + size > 1.0f) {
-                position[i] = 2.0f * (1.0f - size) - position[i];
-                velocity[i] *= -1.0f;
-            }
-        }
-
-        barrier();
-        memoryBarrierBuffer();
-
-        particles[gl_GlobalInvocationID.x].minimal.position = position;
-        particles[gl_GlobalInvocationID.x].minimal.velocity = velocity;
-    }
-
-    barrier();
-    memoryBarrierBuffer();
-}
\ No newline at end of file

projects/wobble_bobble/shaders/update_particle_velocities.comp

@@ -15,8 +15,12 @@ layout(set=0, binding=2) uniform texture3D gridOldImage;
 layout(set=0, binding=3) uniform sampler gridSampler;
 
 layout( push_constant ) uniform constants {
+    float lame1;
+    float lame2;
     float alpha;
     float beta;
+    float t;
+    float dt;
 };
 
 void main()	{
@@ -24,16 +28,21 @@ void main()	{
     memoryBarrierImage();
 
     if (gl_GlobalInvocationID.x < particles.length()) {
-        ParticleMinimal minimal = particles[gl_GlobalInvocationID.x].minimal;
+        Particle particle = particles[gl_GlobalInvocationID.x];
+
+        vec3 position = particle.minimal.position;
+        float size = particle.minimal.size;
+        float volume = sphere_volume(size);
+        float mass = particle.minimal.mass;
 
         ivec3 gridResolution = textureSize(sampler3D(gridImage, gridSampler), 0);
-        ivec3 gridWindow = ivec3(minimal.size * 2.0f * gridResolution);
+        ivec3 gridWindow = ivec3(size * 2.0f * gridResolution);
 
         mat3 affine_D = mat3(0.0f);
         mat3 affine_B = mat3(0.0f);
 
         vec3 velocity_pic = vec3(0.0f);
-        vec3 velocity_flip = vec3(minimal.velocity);
+        vec3 velocity_flip = vec3(particle.minimal.velocity);
 
         int i, j, k;
 
@@ -41,12 +50,12 @@ void main()	{
             for (j = -gridWindow.y; j <= gridWindow.y; j++) {
                 for (k = -gridWindow.z; k <= gridWindow.z; k++) {
                     vec3 offset = vec3(i, j, k) / gridResolution;
-                    vec3 voxel = minimal.position + offset;
+                    vec3 voxel = position + offset;
 
                     vec4 gridSample = texture(sampler3D(gridImage, gridSampler), voxel);
                     vec4 gridOldSample = texture(sampler3D(gridOldImage, gridSampler), voxel);
 
-                    float weight = voxel_particle_weight(voxel, minimal);
+                    float weight = voxel_particle_weight(voxel, particle.minimal);
                     vec3 velocity = gridSample.xyz * weight;
 
                     affine_D += outerProduct(weight * offset, offset);
@@ -61,30 +70,52 @@ void main()	{
         barrier();
         memoryBarrierBuffer();
 
+        mat3 mls_Q = mat3(0.0f);
         mat3 affine_C = mat3(0.0f);
-        vec3 velocity_apic = vec3(0.0f);
+
+        mat3 F = mat3(particle.deformation);
 
         if (abs(determinant(affine_D)) > 0.0f) {
-            affine_C = affine_B * inverse(affine_D);
+            mat3 D_inv = inverse(affine_D);
+            float J = determinant(F);
+
+            if ((J > 0.0f) && (mass > 0.0f)) {
+                mat3 F_T = transpose(F);
+                mat3 F_invT = inverse(F_T);
+
+                mat3 delta = lame2 * (F - F_invT) + lame1 * log(J) * F_invT;
+
+                //mls_Q += beta * dt * volume * particle.minimal.weight_sum * delta * F_T * D_inv / mass;
+            }
+
+            affine_C = affine_B * D_inv;
+            mls_Q += beta * affine_C;
         }
 
-        for (i = -gridWindow.x; i <= gridWindow.x; i++) {
-            for (j = -gridWindow.y; j <= gridWindow.y; j++) {
-                for (k = -gridWindow.z; k <= gridWindow.z; k++) {
-                    vec3 offset = vec3(i, j, k) / gridResolution;
+        F = (mat3(1.0f) + dt * affine_C) * F;
 
-                    velocity_apic += affine_C * offset;
-                }
+        vec3 velocity = mix(velocity_pic, velocity_flip, alpha);
+
+        position = position + velocity * dt;
+
+        for (uint i = 0; i < 3; i++) {
+            if (position[i] - size < 0.0f) {
+                position[i] = -position[i] + 2.0f * size;
+                velocity[i] *= -1.0f;
+            } else
+            if (position[i] + size > 1.0f) {
+                position[i] = 2.0f * (1.0f - size) - position[i];
+                velocity[i] *= -1.0f;
             }
         }
 
         barrier();
         memoryBarrierBuffer();
 
-        vec3 velocity_alpha = mix(velocity_pic, velocity_flip, alpha);
-        vec3 velocity_beta = mix(velocity_alpha, velocity_apic, beta);
-
-        particles[gl_GlobalInvocationID.x].minimal.velocity = velocity_beta;
+        particles[gl_GlobalInvocationID.x].minimal.position = position;
+        particles[gl_GlobalInvocationID.x].minimal.velocity = velocity;
+        particles[gl_GlobalInvocationID.x].deformation = mat4(F);
+        particles[gl_GlobalInvocationID.x].mls = mat4(mls_Q);
     }
 
     barrier();

projects/wobble_bobble/src/main.cpp

@@ -17,18 +17,16 @@ struct Particle {
 	float weight_sum;
 	
 	glm::mat4 deformation;
+	glm::mat4 mls;
 };
 
 struct Physics {
 	float lame1;
     float lame2;
-	float t;
-	float dt;
-};
-
-struct Tweaks {
 	float alpha;
 	float beta;
+	float t;
+	float dt;
 };
 
 float sphere_volume(float radius) {
@@ -80,6 +78,8 @@ void distributeParticles(Particle *particles, size_t count, const glm::vec3& cen
 		
 		particles[i].pad = glm::vec3(0.0f);
 		particles[i].weight_sum = 1.0f;
+		
+		particles[i].mls = glm::mat4(0.0f);
 	}
 }
 
@@ -268,21 +268,6 @@ int main(int argc, const char **argv) {
 		core.writeDescriptorSet(updateGridForcesSets[0], writes);
 	}
 	
-	std::vector<vkcv::DescriptorSetHandle> updateParticleDeformationSets;
-	vkcv::ComputePipelineHandle updateParticleDeformationPipeline = createComputePipeline(
-			core, compiler,
-			"shaders/update_particle_deformation.comp",
-			updateParticleDeformationSets
-	);
-	
-	{
-		vkcv::DescriptorWrites writes;
-		writes.storageBufferWrites.push_back(vkcv::BufferDescriptorWrite(0, particles.getHandle()));
-		writes.sampledImageWrites.push_back(vkcv::SampledImageDescriptorWrite(1, grid.getHandle()));
-		writes.samplerWrites.push_back(vkcv::SamplerDescriptorWrite(2, gridSampler));
-		core.writeDescriptorSet(updateParticleDeformationSets[0], writes);
-	}
-	
 	std::vector<vkcv::DescriptorSetHandle> updateParticleVelocitiesSets;
 	vkcv::ComputePipelineHandle updateParticleVelocitiesPipeline = createComputePipeline(
 			core, compiler,
@@ -299,19 +284,6 @@ int main(int argc, const char **argv) {
 		core.writeDescriptorSet(updateParticleVelocitiesSets[0], writes);
 	}
 	
-	std::vector<vkcv::DescriptorSetHandle> updateParticlePositionsSets;
-	vkcv::ComputePipelineHandle updateParticlePositionsPipeline = createComputePipeline(
-			core, compiler,
-			"shaders/update_particle_positions.comp",
-			updateParticlePositionsSets
-	);
-	
-	{
-		vkcv::DescriptorWrites writes;
-		writes.storageBufferWrites.push_back(vkcv::BufferDescriptorWrite(0, particles.getHandle()));
-		core.writeDescriptorSet(updateParticlePositionsSets[0], writes);
-	}
-	
 	vkcv::ShaderProgram gfxProgramGrid;
 	
 	compiler.compile(
@@ -600,19 +572,14 @@ int main(int argc, const char **argv) {
 		Physics physics;
 		physics.lame1 = material.m_lame1;
 		physics.lame2 = material.m_lame2;
+		physics.alpha = alpha;
+		physics.beta = beta;
 		physics.t = static_cast<float>(0.000001 * static_cast<double>(time.count()));
 		physics.dt = static_cast<float>(0.000001 * static_cast<double>(deltatime.count()));
 		
 		vkcv::PushConstants physicsPushConstants (sizeof(physics));
 		physicsPushConstants.appendDrawcall(physics);
 		
-		Tweaks tweaks;
-		tweaks.alpha = alpha;
-		tweaks.beta = beta;
-		
-		vkcv::PushConstants tweakPushConstants (sizeof(tweaks));
-		tweakPushConstants.appendDrawcall(tweaks);
-		
 		cameraManager.update(physics.dt);
 
 		glm::mat4 mvp = cameraManager.getActiveCamera().getMVP();
@@ -676,23 +643,6 @@ int main(int argc, const char **argv) {
 		core.recordImageMemoryBarrier(cmdStream, grid.getHandle());
 		core.recordEndDebugLabel(cmdStream);
 		
-		core.recordBeginDebugLabel(cmdStream, "UPDATE PARTICLE DEFORMATION", { 0.78f, 0.89f, 0.94f, 1.0f });
-		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
-		core.prepareImageForSampling(cmdStream, grid.getHandle());
-		
-		core.recordComputeDispatchToCmdStream(
-				cmdStream,
-				updateParticleDeformationPipeline,
-				dispatchSizeParticles,
-				{ vkcv::DescriptorSetUsage(
-						0, core.getDescriptorSet(updateParticleDeformationSets[0]).vulkanHandle
-				) },
-				physicsPushConstants
-		);
-		
-		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
-		core.recordEndDebugLabel(cmdStream);
-		
 		core.recordBeginDebugLabel(cmdStream, "UPDATE PARTICLE VELOCITIES", { 0.78f, 0.89f, 0.94f, 1.0f });
 		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
 		core.prepareImageForSampling(cmdStream, grid.getHandle());
@@ -705,22 +655,6 @@ int main(int argc, const char **argv) {
 				{ vkcv::DescriptorSetUsage(
 						0, core.getDescriptorSet(updateParticleVelocitiesSets[0]).vulkanHandle
 				) },
-				tweakPushConstants
-		);
-		
-		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
-		core.recordEndDebugLabel(cmdStream);
-		
-		core.recordBeginDebugLabel(cmdStream, "UPDATE PARTICLE POSITIONS", { 0.78f, 0.89f, 0.94f, 1.0f });
-		core.recordBufferMemoryBarrier(cmdStream, particles.getHandle());
-		
-		core.recordComputeDispatchToCmdStream(
-				cmdStream,
-				updateParticlePositionsPipeline,
-				dispatchSizeParticles,
-				{ vkcv::DescriptorSetUsage(
-						0, core.getDescriptorSet(updateParticlePositionsSets[0]).vulkanHandle
-				) },
 				physicsPushConstants
 		);