Vulkan_ Mass Spring Systems_ Mass spring model)

The common method of cloth simulation is to express cloth as three-dimensional mesh, and then simulate cloth by elastic force and external force. It has a wide range of applications in the game, such as the swing of sleeves, the flutter of flags; in addition, it can also be used to simulate fitting and so on.

(pictures of cloth under gravity and wind only)

1, Mass Spring Systems

The simplest way to simulate a deformed object is to represent it as a mass spring system. Because the method is simple, this framework is ideal for learning various techniques of physical simulation and time integration. A spring mass contains a series of mass connected by multiple springs, so the physical properties of the system are very direct, and the simulation program is easy to write.

1.1 Mass Spring Systems

The point spring model is a model that the mass of an object is dispersed to a particle, and the spring connection is used between the particles to represent the physical characteristics of an object, such as elasticity, damping force, etc. In order to improve the degree of simulation, researchers have made many improvements on the basis of this simple model. There are three kinds of springs that connect particles.

The first is structural spring, which is used to connect the transverse and longitudinal particles, and plays a role of fixing the model structure. The second is the torsion spring, also known as the shear spring, which connects the adjacent particles on the diagonal to prevent the distortion of the model. The third is the tensile spring, also known as the bending spring, which connects two particles separated by one particle horizontally and longitudinally to ensure the smooth edge of the model during deformation (such as cloth folding).

1.2 physical model (damping force and elasticity)

Suppose there are N particles with mass mi, position xi, velocity VI, 1 < I < n

These particles are connected by a set of springs S, and the spring parameters are (i,j, lo, ks, kd). i,j are connected spring particles, Lo is the length when the spring is completely relaxed, ks is the spring elastic coefficient, kd is the damping coefficient. According to huco'S law, the force exerted by the spring on the two vertices can be expressed as:

From the conservation of force, we know that fi + fj = 0. The size of FI and fj is proportional to the elongation of spring.
For the calculation of damping:

For a particle, the force equation is as follows:

1.3 simulation algorithm

In the simulation algorithm, Newton's second law is the key, f = ma. In the case of known mass and external force, the acceleration can be obtained by a=f/m. The second order ordinary differential equation is written as two first order equations:

The analytical solution can be obtained:

The initial state is v(to) = vo, x(to)=xo

The integration adds all the changes in time t, and the simulation process starts from to to to continuously calculate x(t) and v(t), and then update the position of the particle.

The simplest numerical solution is to use the finite difference approximate differentiation:

Where delta t is the time step, and t is the frame number. If the previous formula is carried in, we can get:

This is the explicit Euler solution. The next state is completely determined by the current state.

The pseudo code of the whole process is as follows:

 // Initialization data
 forall particles i
         initialize xi , vi and mi
 endfor
 // Cyclic simulation
 loop
        forall particles i
               fi ← fg + fcoll +f(xi , vi , x j , v j )
        endfor
        forall particles i
             vi ← vi + ∆t fi /mi
             xi ← xi + ∆t vi
         endfor
         display the system every nth time
 endloop

fg is gravity, fcoll is the external force of collision.

Explicit integration is the simplest integration method, but there is a serious problem - it can only be stable when the step size is very small, which means that multiple simulation steps need to be carried out between two visible frames. The main reason for the instability is that the Euler method only depends on the current state. It assumes that the external force is a constant in the whole step size.

1.4 cloth simulation

In the mass spring model, the whole cloth is represented by mesh points, and the line between points is called spring.

The position of a cloth point is controlled by 12 points around it. Eight of them are adjacent points, and the other four are adjacent points of the interval.

For spring, there are three kinds of spring, one is structural spring, one is shear spring, the other is bending spring. They correspond to the above 12 springs respectively.

Among them, the structural spring simulates the transverse force, in order to prevent the cloth from generating large tension and compression deformation in the transverse direction.

Similarly, the shear spring is to prevent the cloth from producing large tension and compression deformation in the oblique direction.

The bending spring connects two separated points. Its main function is to prevent the excessive bending of the cloth when it is bent and deformed, resulting in the collapse of the cloth angle. Sometimes this coefficient can be ignored.

In the specific code implementation, we initialize the position, speed (initialized to 0), quality (set to a unified value), and whether it is fixed (fixed will not be subject to external forces, that is, gravity).

2, vulkan implementation

In this part, we use GPU to calculate the shader and integrate the spring force to realize the cloth simulation model. At the same time, we also realize the basic collision between a cloth and a fixed scene object.

2.1 scene data structure

First, customize the calculation pipeline, graphic pipeline and cloth related structures for the scene:

	//scene
	uint32_t sceneSetup = 1; //0 cloth collision scene, 1 pure cloth simulation
	bool specializedComputeQueue = false; //Follow up on its use
	//Resources for the drawing route section
	struct {
		VkDescriptorSetLayout descriptorSetLayout;
		VkDescriptorSet descriptorSet;
		VkPipelineLayout pipelineLayout;
		struct Pipelines {
			VkPipeline cloth;
			VkPipeline sphere;
		} pipelines;
		vks::Buffer indices;
		vks::Buffer uniformBuffer;
		struct graphicsUBO {
			glm::mat4 projection;
			glm::mat4 view;
			glm::vec4 lightPos = glm::vec4(-1.0f, 3.0f, -1.0f, 1.0f);
		} ubo;
	} graphics;

	//Calculate resources for pipeline section
	struct {
		struct StorageBuffers {
			vks::Buffer input;
			vks::Buffer output;
		} storageBuffers;
		struct Semaphores {
			VkSemaphore ready{ 0L };
			VkSemaphore complete{ 0L };
		} semaphores;
		vks::Buffer uniformBuffer;
		VkQueue queue;
		VkCommandPool commandPool;
		std::array<VkCommandBuffer,2> commandBuffers;
		VkDescriptorSetLayout descriptorSetLayout;
		std::array<VkDescriptorSet,2> descriptorSets;
		VkPipelineLayout pipelineLayout;
		VkPipeline pipeline;
		struct computeUBO {
			float deltaT = 0.0f;                                     //time step 
			float particleMass = 0.1f;                               //Particle size
			float springStiffness = 8000.0f;                         //Spring stiffness
			float damping = 0.25f;                                   //attenuation
			float restDistH;                                         //Fixed length of horizontal spring
			float restDistV;                                         //Fixed length of vertical spring
			float restDistD;				                         //Fixed length of diagonal spring
			float sphereRadius = 0.5f;                               //Sphere radius
			glm::vec4 spherePos = glm::vec4(0.0f, 0.0f, 0.0f, 0.0f); //Sphere position
			glm::vec4 gravity = glm::vec4(0.0f, 9.8f, 0.0f, 0.0f);   //gravity
			glm::ivec2 particleCount;                                //Number of particles
		} ubo;
	} compute;

	// SSBO cloth mesh particle declaration
	struct Particle {
		glm::vec4 pos;
		glm::vec4 vel;
		glm::vec4 uv;
		glm::vec4 normal;
		float pinned;
		glm::vec3 _pad0;
	};
	
	//Clothing grid properties
	struct Cloth {
		glm::uvec2 gridsize = glm::uvec2(50, 50);
		glm::vec2 size = glm::vec2(2.0f, 2.0f);
	} cloth;

2.2 storage buffer data creation

After defining the scene data, we need to create a prepareStorageBuffers function to fill the calculation shader storage buffer:

	//Set and fill the compute shader store buffer with particles
	void prepareStorageBuffers()
	{
		std::vector<Particle> particleBuffer(cloth.gridsize.x *  cloth.gridsize.y);

		float dx =  cloth.size.x / (cloth.gridsize.x - 1);
		float dy =  cloth.size.y / (cloth.gridsize.y - 1);
		float du = 1.0f / (cloth.gridsize.x - 1);
		float dv = 1.0f / (cloth.gridsize.y - 1);
		switch (sceneSetup) {
			case 0 :
			{
				//  Horizontal adjustment. Cloth falling on the ball
				glm::mat4 transM = glm::translate(glm::mat4(1.0f), glm::vec3(- cloth.size.x / 2.0f, -2.0f, - cloth.size.y / 2.0f));
				for (uint32_t i = 0; i <  cloth.gridsize.y; i++) {
					for (uint32_t j = 0; j <  cloth.gridsize.x; j++) {
						particleBuffer[i + j * cloth.gridsize.y].pos = transM * glm::vec4(dx * j, 0.0f, dy * i, 1.0f);
						particleBuffer[i + j * cloth.gridsize.y].vel = glm::vec4(0.0f);
						particleBuffer[i + j * cloth.gridsize.y].uv = glm::vec4(1.0f - du * i, dv * j, 0.0f, 0.0f);
					}
				}
				break;
			}
			case 1:
			{
				// Vertical. Fixed cloth
				glm::mat4 transM = glm::translate(glm::mat4(1.0f), glm::vec3(- cloth.size.x / 2.0f, - cloth.size.y / 2.0f, 0.0f));
				for (uint32_t i = 0; i <  cloth.gridsize.y; i++) {
					for (uint32_t j = 0; j <  cloth.gridsize.x; j++) {
						particleBuffer[i + j * cloth.gridsize.y].pos = transM * glm::vec4(dx * j, dy * i, 0.0f, 1.0f);
						particleBuffer[i + j * cloth.gridsize.y].vel = glm::vec4(0.0f);
						particleBuffer[i + j * cloth.gridsize.y].uv = glm::vec4(du * j, dv * i, 0.0f, 0.0f);
						// Fixed particle
						//particleBuffer[i + j *  cloth.gridsize .y].pinned = (i == 0) && ((j == 0) || (j ==   cloth.gridsize .x / 4) || (j ==   cloth.gridsize .x -   cloth.gridsize .x / 4) || (j ==   cloth.gridsize . X - 1)); / / four points
						//particleBuffer[i + j *  cloth.gridsize .y].pinned = (i == 0) && ((j == 0) || (j ==   cloth.gridsize .x / 2) || (j ==   cloth.gridsize . X - 1)); / / three points
						particleBuffer[i + j * cloth.gridsize.y].pinned = (i == 0) && ((j == 0) || (j ==  cloth.gridsize.x - 1));  //Two points
						// Moving the sphere to eliminate the impact of sphere collision
						compute.ubo.spherePos.z = -100.0f;
					}
				}
				break;
			}
		}

		VkDeviceSize storageBufferSize = particleBuffer.size() * sizeof(Particle);
		
		// Store data
		// After uploading, the SSBO on the host will not be changed, so it will be copied to the local memory of the device
		vks::Buffer stagingBuffer;

		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_TRANSFER_SRC_BIT,
			VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT,
			&stagingBuffer,
			storageBufferSize,
			particleBuffer.data());

		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_VERTEX_BUFFER_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT | VK_BUFFER_USAGE_TRANSFER_DST_BIT,
			VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT,
			&compute.storageBuffers.input,
			storageBufferSize);

		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_VERTEX_BUFFER_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT | VK_BUFFER_USAGE_TRANSFER_DST_BIT,
			VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT,
			&compute.storageBuffers.output,
			storageBufferSize);

		// Copy from staging buffer
		VkCommandBuffer copyCmd = vulkanDevice->createCommandBuffer(VK_COMMAND_BUFFER_LEVEL_PRIMARY, true);
		VkBufferCopy copyRegion = {};
		copyRegion.size = storageBufferSize;
		vkCmdCopyBuffer(copyCmd, stagingBuffer.buffer, compute.storageBuffers.input.buffer, 1, &copyRegion);
		vkCmdCopyBuffer(copyCmd, stagingBuffer.buffer, compute.storageBuffers.output.buffer, 1, &copyRegion);

		//Add an initial release barrier to the graph queue so that when the compute command buffer executes for the first time, it does not lack a release corresponding to its get
		addGraphicsToComputeBarriers(copyCmd);
		vulkanDevice->flushCommandBuffer(copyCmd, queue, true);

		stagingBuffer.destroy();

		// Indexes index data
		std::vector<uint32_t> indices;
		for (uint32_t y = 0; y <  cloth.gridsize.y - 1; y++) {
			for (uint32_t x = 0; x <  cloth.gridsize.x; x++) {
				indices.push_back((y + 1) *  cloth.gridsize.x + x);
				indices.push_back((y)*  cloth.gridsize.x + x);
			}
			// Primitive restart (filled with special value 0xFFFFFFFF)
			indices.push_back(0xFFFFFFFF);
		}
		uint32_t indexBufferSize = static_cast<uint32_t>(indices.size()) * sizeof(uint32_t);
		indexCount = static_cast<uint32_t>(indices.size());

		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_TRANSFER_SRC_BIT,
			VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT,
			&stagingBuffer,
			indexBufferSize,
			indices.data());

		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_INDEX_BUFFER_BIT | VK_BUFFER_USAGE_TRANSFER_DST_BIT,
			VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT,
			&graphics.indices,
			indexBufferSize);

		// Copy from staging buffer
		copyCmd = vulkanDevice->createCommandBuffer(VK_COMMAND_BUFFER_LEVEL_PRIMARY, true);
		copyRegion = {};
		copyRegion.size = indexBufferSize;
		vkCmdCopyBuffer(copyCmd, stagingBuffer.buffer, graphics.indices.buffer, 1, &copyRegion);
		vulkanDevice->flushCommandBuffer(copyCmd, queue, true);

		stagingBuffer.destroy();
	}

In this function, we can see that we call an addGraphicsToComputeBarriers method to add an initial release barrier for the graph queue. The specific implementation is as follows:

	//Add a barrier to the graphics queue
	void addGraphicsToComputeBarriers(VkCommandBuffer commandBuffer) 
	{
		if (specializedComputeQueue) {
			VkBufferMemoryBarrier bufferBarrier = vks::initializers::bufferMemoryBarrier();
			bufferBarrier.srcAccessMask = VK_ACCESS_VERTEX_ATTRIBUTE_READ_BIT;
			bufferBarrier.dstAccessMask = VK_ACCESS_SHADER_WRITE_BIT;
			bufferBarrier.srcQueueFamilyIndex = vulkanDevice->queueFamilyIndices.graphics;
			bufferBarrier.dstQueueFamilyIndex = vulkanDevice->queueFamilyIndices.compute;
			bufferBarrier.size = VK_WHOLE_SIZE;

			std::vector<VkBufferMemoryBarrier> bufferBarriers;
			bufferBarrier.buffer = compute.storageBuffers.input.buffer;
			bufferBarriers.push_back(bufferBarrier);
			bufferBarrier.buffer = compute.storageBuffers.output.buffer;
			bufferBarriers.push_back(bufferBarrier);
			vkCmdPipelineBarrier(commandBuffer,
				VK_PIPELINE_STAGE_VERTEX_INPUT_BIT,
				VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
				VK_FLAGS_NONE,
				0, nullptr,
				static_cast<uint32_t>(bufferBarriers.size()), bufferBarriers.data(),
				0, nullptr);
		}
	} 

Where specializedComputeQueue is a global function.

		//Ensure that the code works for different drawing and calculation queue families and the same queue family
#ifdef DEBUG_FORCE_SHARED_GRAPHICS_COMPUTE_QUEUE
		vulkanDevice->queueFamilyIndices.compute = vulkanDevice->queueFamilyIndices.graphics;
#endif
		//Check whether the calculation queue family is different from the drawing queue family
		specializedComputeQueue = vulkanDevice->queueFamilyIndices.graphics != vulkanDevice->queueFamilyIndices.compute;

At the same time, we need to set up the other two pipeline barriers to deal with resource problems when different pipelines are used. In the future, we will use them one by one and detail:

	//Compute shader assignment do business compute command add barrier
	void addComputeToComputeBarriers(VkCommandBuffer commandBuffer) 
	{
		VkBufferMemoryBarrier bufferBarrier = vks::initializers::bufferMemoryBarrier();
		bufferBarrier.srcAccessMask = VK_ACCESS_SHADER_WRITE_BIT;
		bufferBarrier.dstAccessMask = VK_ACCESS_SHADER_READ_BIT;
		bufferBarrier.srcQueueFamilyIndex = vulkanDevice->queueFamilyIndices.compute;
		bufferBarrier.dstQueueFamilyIndex = vulkanDevice->queueFamilyIndices.compute;
		bufferBarrier.size = VK_WHOLE_SIZE;
		std::vector<VkBufferMemoryBarrier> bufferBarriers;
		bufferBarrier.buffer = compute.storageBuffers.input.buffer;
		bufferBarriers.push_back(bufferBarrier);
		bufferBarrier.buffer = compute.storageBuffers.output.buffer;
		bufferBarriers.push_back(bufferBarrier);
		vkCmdPipelineBarrier(
			commandBuffer,
			VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
			VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
			VK_FLAGS_NONE,
			0, nullptr,
			static_cast<uint32_t>(bufferBarriers.size()), bufferBarriers.data(),
			0, nullptr);
	} 
	//Use when the storage buffer is released back to the drawing queue
	void addComputeToGraphicsBarriers(VkCommandBuffer commandBuffer)
	{
		if (specializedComputeQueue) {
			VkBufferMemoryBarrier bufferBarrier = vks::initializers::bufferMemoryBarrier();
			bufferBarrier.srcAccessMask = VK_ACCESS_SHADER_WRITE_BIT;
			bufferBarrier.dstAccessMask = VK_ACCESS_VERTEX_ATTRIBUTE_READ_BIT;
			bufferBarrier.srcQueueFamilyIndex = vulkanDevice->queueFamilyIndices.compute;
			bufferBarrier.dstQueueFamilyIndex = vulkanDevice->queueFamilyIndices.graphics;
			bufferBarrier.size = VK_WHOLE_SIZE;
			std::vector<VkBufferMemoryBarrier> bufferBarriers;
			bufferBarrier.buffer = compute.storageBuffers.input.buffer;
			bufferBarriers.push_back(bufferBarrier);
			bufferBarrier.buffer = compute.storageBuffers.output.buffer;
			bufferBarriers.push_back(bufferBarrier);
			vkCmdPipelineBarrier(
				commandBuffer,
				VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
				VK_PIPELINE_STAGE_VERTEX_INPUT_BIT,
				VK_FLAGS_NONE,
				0, nullptr,
				static_cast<uint32_t>(bufferBarriers.size()), bufferBarriers.data(),
				0, nullptr);
		}
	} 

2.3 creating uniformbuffers

After creating the storage buffer data, we need to prepare and initialize UniformBuffer for shader provisioning:

	// Preparing and initializing a Uniform buffer containing shaders
	void prepareUniformBuffers()
	{
		// Calculate shading uniform buffer block
		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT,
			VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT,
			&compute.uniformBuffer,
			sizeof(compute.ubo));
		VK_CHECK_RESULT(compute.uniformBuffer.map());
	
		// Initial values
		float dx = cloth.size.x / (cloth.gridsize.x - 1);
		float dy = cloth.size.y / (cloth.gridsize.y - 1);
		//Spring system length definition
		compute.ubo.restDistH = dx;
		compute.ubo.restDistV = dy;
		compute.ubo.restDistD = sqrtf(dx * dx + dy * dy);
		compute.ubo.particleCount = cloth.gridsize;

		updateComputeUBO();

		// Vertex shading uniform buffer block
		vulkanDevice->createBuffer(
			VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT,
			VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT,
			&graphics.uniformBuffer,
			sizeof(graphics.ubo));
		VK_CHECK_RESULT(graphics.uniformBuffer.map());

		updateGraphicsUBO();
	}
	//Update Uniform data in compute shaders
	void updateComputeUBO()
	{
		if (!paused) {
			compute.ubo.deltaT = 0.000005f;
			//compute.ubo.deltaT = frameTimer * 0.0075f;
			//Simulate adding wind
			if (simulateWind) {
				std::default_random_engine rndEngine(benchmark.active ? 0 : (unsigned)time(nullptr));
				std::uniform_real_distribution<float> rd(1.0f, 6.0f);
				compute.ubo.gravity.x = cos(glm::radians(-timer * 360.0f)) * (rd(rndEngine) - rd(rndEngine));
				compute.ubo.gravity.z = sin(glm::radians(timer * 360.0f)) * (rd(rndEngine) - rd(rndEngine));
			}
			else {
				compute.ubo.gravity.x = 1.0f;
				compute.ubo.gravity.z = 1.0f;
			}
		}
		else {
			compute.ubo.deltaT = 0.0f;
		}
		memcpy(compute.uniformBuffer.mapped, &compute.ubo, sizeof(compute.ubo));
	}
	//Update Uniform data in graphics shaders
	void updateGraphicsUBO()
	{
		graphics.ubo.projection = camera.matrices.perspective;
		graphics.ubo.view = camera.matrices.view;
		memcpy(graphics.uniformBuffer.mapped, &graphics.ubo, sizeof(graphics.ubo));
	}

2.4 graphic pipeline and shader

The next step is to create descriptors, descriptor pools, descriptor sets, graphic pipelines and other conventional operations. Among the cloth and ball pipelines, there is only a simple von's lighting model, which is not detailed

2.5 mass spring system treatment

Now, we need to create a separate calculation pipeline to deal with the mass spring. For example, at the initial time, the mass state of a cloth is as follows:

After calculating the force of adjacent particles processed by the pipeline, all particles will deform under the force as shown in the following figure:

We first need to create a separate calculation shader to handle the calculation of the spring system:

	void prepareCompute()
	{
		// Create a device queue with computing power
		vkGetDeviceQueue(device, vulkanDevice->queueFamilyIndices.compute, 0, &compute.queue);

		// Create calculation pipeline
		std::vector<VkDescriptorSetLayoutBinding> setLayoutBindings = {
			vks::initializers::descriptorSetLayoutBinding(VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0),
			vks::initializers::descriptorSetLayoutBinding(VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 1),
			vks::initializers::descriptorSetLayoutBinding(VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 2),
		};

		VkDescriptorSetLayoutCreateInfo descriptorLayout =
			vks::initializers::descriptorSetLayoutCreateInfo(setLayoutBindings);

		VK_CHECK_RESULT(vkCreateDescriptorSetLayout(device,	&descriptorLayout, nullptr,	&compute.descriptorSetLayout));

		VkPipelineLayoutCreateInfo pipelineLayoutCreateInfo =
			vks::initializers::pipelineLayoutCreateInfo(&compute.descriptorSetLayout, 1);

		// A push constant used to pass some parameters
		VkPushConstantRange pushConstantRange =
			vks::initializers::pushConstantRange(VK_SHADER_STAGE_COMPUTE_BIT, sizeof(uint32_t), 0);
		pipelineLayoutCreateInfo.pushConstantRangeCount = 1;
		pipelineLayoutCreateInfo.pPushConstantRanges = &pushConstantRange;

		VK_CHECK_RESULT(vkCreatePipelineLayout(device, &pipelineLayoutCreateInfo, nullptr,	&compute.pipelineLayout));

		VkDescriptorSetAllocateInfo allocInfo =
			vks::initializers::descriptorSetAllocateInfo(descriptorPool, &compute.descriptorSetLayout, 1);

		// Create two descriptor sets for switching input and output buffers
		VK_CHECK_RESULT(vkAllocateDescriptorSets(device, &allocInfo, &compute.descriptorSets[0]));
		VK_CHECK_RESULT(vkAllocateDescriptorSets(device, &allocInfo, &compute.descriptorSets[1]));

		std::vector<VkWriteDescriptorSet> computeWriteDescriptorSets = {
			vks::initializers::writeDescriptorSet(compute.descriptorSets[0], VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 0, &compute.storageBuffers.input.descriptor),
			vks::initializers::writeDescriptorSet(compute.descriptorSets[0], VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 1, &compute.storageBuffers.output.descriptor),
			vks::initializers::writeDescriptorSet(compute.descriptorSets[0], VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER, 2, &compute.uniformBuffer.descriptor),

			vks::initializers::writeDescriptorSet(compute.descriptorSets[1], VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 0, &compute.storageBuffers.output.descriptor),
			vks::initializers::writeDescriptorSet(compute.descriptorSets[1], VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 1, &compute.storageBuffers.input.descriptor),
			vks::initializers::writeDescriptorSet(compute.descriptorSets[1], VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER, 2, &compute.uniformBuffer.descriptor)
		};

		vkUpdateDescriptorSets(device, static_cast<uint32_t>(computeWriteDescriptorSets.size()), computeWriteDescriptorSets.data(), 0, NULL);

		// Create pipe
		VkComputePipelineCreateInfo computePipelineCreateInfo = vks::initializers::computePipelineCreateInfo(compute.pipelineLayout, 0);
		computePipelineCreateInfo.stage = loadShader(getAssetPath() + "shaders/computecloth/cloth.comp.spv", VK_SHADER_STAGE_COMPUTE_BIT);
		VK_CHECK_RESULT(vkCreateComputePipelines(device, pipelineCache, 1, &computePipelineCreateInfo, nullptr, &compute.pipeline));

		// A separate command pool as a family of computing queues may not be the same as a drawing
		VkCommandPoolCreateInfo cmdPoolInfo = {};
		cmdPoolInfo.sType = VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO;
		cmdPoolInfo.queueFamilyIndex = vulkanDevice->queueFamilyIndices.compute;
		cmdPoolInfo.flags = VK_COMMAND_POOL_CREATE_RESET_COMMAND_BUFFER_BIT;
		VK_CHECK_RESULT(vkCreateCommandPool(device, &cmdPoolInfo, nullptr, &compute.commandPool));

		// Create a command buffer for the calculation operation
		VkCommandBufferAllocateInfo cmdBufAllocateInfo =
			vks::initializers::commandBufferAllocateInfo(compute.commandPool, VK_COMMAND_BUFFER_LEVEL_PRIMARY, 2);	

		VK_CHECK_RESULT(vkAllocateCommandBuffers(device, &cmdBufAllocateInfo, &compute.commandBuffers[0]));

		// Graph / compute synchronized semaphores
		VkSemaphoreCreateInfo semaphoreCreateInfo = vks::initializers::semaphoreCreateInfo();
		VK_CHECK_RESULT(vkCreateSemaphore(device, &semaphoreCreateInfo, nullptr, &compute.semaphores.ready));
		VK_CHECK_RESULT(vkCreateSemaphore(device, &semaphoreCreateInfo, nullptr, &compute.semaphores.complete));

		// Building a single command buffer that contains compute dispatch commands
		buildComputeCommandBuffer();
	}

In addition, we need to create a separate command buffer to calculate the mass spring:

	void buildComputeCommandBuffer()
	{
		VkCommandBufferBeginInfo cmdBufInfo = vks::initializers::commandBufferBeginInfo();
		cmdBufInfo.flags = VK_COMMAND_BUFFER_USAGE_SIMULTANEOUS_USE_BIT;

		for (uint32_t i = 0; i < 2; i++) {

			VK_CHECK_RESULT(vkBeginCommandBuffer(compute.commandBuffers[i], &cmdBufInfo));

			// Get storage buffer from graph queue
			addGraphicsToComputeBarriers(compute.commandBuffers[i]);

			vkCmdBindPipeline(compute.commandBuffers[i], VK_PIPELINE_BIND_POINT_COMPUTE, compute.pipeline);

			uint32_t calculateNormals = 0;
			vkCmdPushConstants(compute.commandBuffers[i], compute.pipelineLayout, VK_SHADER_STAGE_COMPUTE_BIT, 0, sizeof(uint32_t), &calculateNormals);

			// Assign calculation jobs
			const uint32_t iterations = 32;
			for (uint32_t j = 0; j < iterations; j++) {
				readSet = 1 - readSet;
				vkCmdBindDescriptorSets(compute.commandBuffers[i], VK_PIPELINE_BIND_POINT_COMPUTE, compute.pipelineLayout, 0, 1, &compute.descriptorSets[readSet], 0, 0);

				if (j == iterations - 1) {
					calculateNormals = 1;
					vkCmdPushConstants(compute.commandBuffers[i], compute.pipelineLayout, VK_SHADER_STAGE_COMPUTE_BIT, 0, sizeof(uint32_t), &calculateNormals);
				}

				vkCmdDispatch(compute.commandBuffers[i], cloth.gridsize.x / 10, cloth.gridsize.y / 10, 1);

				//Do not add obstacles to the last iteration of the loop, because we will explicitly release the graph queue
				if (j != iterations - 1) {
					addComputeToComputeBarriers(compute.commandBuffers[i]);
				}

			}

			//Release the storage buffer back to the drawing queue
			addComputeToGraphicsBarriers(compute.commandBuffers[i]);
			vkEndCommandBuffer(compute.commandBuffers[i]);
		}
	}

2.6 calculation shader processing spring system

Next let's focus on the compute shaders used in 2.5 cloth.comp :

#version 450
struct Particle {
	vec4 pos;
	vec4 vel;
	vec4 uv;
	vec4 normal;
	float pinned;
};

layout(std430, binding = 0) buffer ParticleIn {
	Particle particleIn[ ];
};

layout(std430, binding = 1) buffer ParticleOut {
	Particle particleOut[ ];
};

layout (local_size_x = 10, local_size_y = 10) in;

layout (binding = 2) uniform UBO 
{
	float deltaT;                //time step 
	float particleMass;          //Particle size 
	float springStiffness;		 //Spring stiffness
	float damping;				 //attenuation
	float restDistH;			 //Fixed length of horizontal spring
	float restDistV;			 //Fixed length of vertical spring
	float restDistD;			 //Fixed length of diagonal spring
	float sphereRadius;			 //Sphere radius
	vec4 spherePos;				 //Sphere position
	vec4 gravity;				 //gravity
	ivec2 particleCount;		 //Number of particles
} params;

layout (push_constant) uniform PushConsts {
	uint calculateNormals;
} pushConsts;

//Find out the force of p1 point on p0 point
vec3 springForce(vec3 p0, vec3 p1, float restDist) 
{
	vec3 dist = p0 - p1;
	return normalize(dist) * params.springStiffness * (length(dist) - restDist);
}

void main() 
{
	uvec3 id = gl_GlobalInvocationID; 

	uint index = id.y * params.particleCount.x + id.x;
	if (index > params.particleCount.x * params.particleCount.y) 
		return;

	// fixed
	if (particleIn[index].pinned == 1.0) {
		particleOut[index].pos = particleOut[index].pos;
		particleOut[index].vel = vec4(0.0);
		return;
	}

	// Initialize gravity
	vec3 force = params.gravity.xyz * params.particleMass;

	vec3 pos = particleIn[index].pos.xyz;
	vec3 vel = particleIn[index].vel.xyz;

	// Spring force from neighboring particles
	// Left
	if (id.x > 0) {
		force += springForce(particleIn[index-1].pos.xyz, pos, params.restDistH);
	} 
	// right
	if (id.x < params.particleCount.x - 1) {
		force += springForce(particleIn[index + 1].pos.xyz, pos, params.restDistH);
	}
	// upper
	if (id.y < params.particleCount.y - 1) {
		force += springForce(particleIn[index + params.particleCount.x].pos.xyz, pos, params.restDistV);
	} 
	// lower
	if (id.y > 0) {
		force += springForce(particleIn[index - params.particleCount.x].pos.xyz, pos, params.restDistV);
	} 
	// Left
	if ((id.x > 0) && (id.y < params.particleCount.y - 1)) {
		force += springForce(particleIn[index + params.particleCount.x - 1].pos.xyz, pos, params.restDistD);
	}
	// Bottom left
	if ((id.x > 0) && (id.y > 0)) {
		force += springForce(particleIn[index - params.particleCount.x - 1].pos.xyz, pos, params.restDistD);
	}
	// Upper right
	if ((id.x < params.particleCount.x - 1) && (id.y < params.particleCount.y - 1)) {
		force += springForce(particleIn[index + params.particleCount.x + 1].pos.xyz, pos, params.restDistD);
	}
	// lower right
	if ((id.x < params.particleCount.x - 1) && (id.y > 0)) {
		force += springForce(particleIn[index - params.particleCount.x + 1].pos.xyz, pos, params.restDistD);
	}

	force += (-params.damping * vel);

	// integration
	vec3 f = force * (1.0 / params.particleMass);
	particleOut[index].pos = vec4(pos + vel * params.deltaT + 0.5 * f * params.deltaT * params.deltaT, 1.0);
	particleOut[index].vel = vec4(vel + f * params.deltaT, 0.0);

	// Sphere collision (applicable to scene 2)
	vec3 sphereDist = particleOut[index].pos.xyz - params.spherePos.xyz;
	if (length(sphereDist) < params.sphereRadius + 0.01) {
		//If the particle is inside the sphere, push it to the outside radius
		particleOut[index].pos.xyz = params.spherePos.xyz + normalize(sphereDist) * (params.sphereRadius + 0.01);		
		// Cancellation speed
		particleOut[index].vel = vec4(0.0);
	}

	// Normal processing
	if (pushConsts.calculateNormals == 1) {
		vec3 normal = vec3(0.0);
		vec3 a, b, c;
		if (id.y > 0) {
			if (id.x > 0) {
				a = particleIn[index - 1].pos.xyz - pos;
				b = particleIn[index - params.particleCount.x - 1].pos.xyz - pos;
				c = particleIn[index - params.particleCount.x].pos.xyz - pos;
				normal += cross(a,b) + cross(b,c);
			}
			if (id.x < params.particleCount.x - 1) {
				a = particleIn[index - params.particleCount.x].pos.xyz - pos;
				b = particleIn[index - params.particleCount.x + 1].pos.xyz - pos;
				c = particleIn[index + 1].pos.xyz - pos;
				normal += cross(a,b) + cross(b,c);
			}
		}
		if (id.y < params.particleCount.y - 1) {
			if (id.x > 0) {
				a = particleIn[index + params.particleCount.x].pos.xyz - pos;
				b = particleIn[index + params.particleCount.x - 1].pos.xyz - pos;
				c = particleIn[index - 1].pos.xyz - pos;
				normal += cross(a,b) + cross(b,c);
			}
			if (id.x < params.particleCount.x - 1) {
				a = particleIn[index + 1].pos.xyz - pos;
				b = particleIn[index + params.particleCount.x + 1].pos.xyz - pos;
				c = particleIn[index + params.particleCount.x].pos.xyz - pos;
				normal += cross(a,b) + cross(b,c);
			}
		}
		particleOut[index].normal = vec4(normalize(normal), 0.0f);
	}
}

We can see that, according to the force of single mass spring in the first part, we have calculated and processed all the neighboring action points of a particle system, and finally reprocessed the normal for sampling, in which we also processed a ball collision (when the pure cloth model, we will distribute the ball principle, and the corresponding processing can be found in the above code).

2.7 rendering a scene

While calculating the shader to calculate the cloth simulation in real time, we need to perform the final rendering, in which we need to use the update vertex buffer to draw the cloth mesh data.

	void buildCommandBuffers()
	{
		VkCommandBufferBeginInfo cmdBufInfo = vks::initializers::commandBufferBeginInfo();

		VkClearValue clearValues[2];
		clearValues[0].color = { { 0.1f, 0.1f, 0.3f, 0.5f } };;
		clearValues[1].depthStencil = { 1.0f, 0 };

		VkRenderPassBeginInfo renderPassBeginInfo = vks::initializers::renderPassBeginInfo();
		renderPassBeginInfo.renderPass = renderPass;
		renderPassBeginInfo.renderArea.offset.x = 0;
		renderPassBeginInfo.renderArea.offset.y = 0;
		renderPassBeginInfo.renderArea.extent.width = width;
		renderPassBeginInfo.renderArea.extent.height = height;
		renderPassBeginInfo.clearValueCount = 2;
		renderPassBeginInfo.pClearValues = clearValues;

		for (int32_t i = 0; i < drawCmdBuffers.size(); ++i)
		{
			// Set target frame buffer
			renderPassBeginInfo.framebuffer = frameBuffers[i];

			VK_CHECK_RESULT(vkBeginCommandBuffer(drawCmdBuffers[i], &cmdBufInfo));

			// Get storage buffer from compute queue
			addComputeToGraphicsBarriers(drawCmdBuffers[i]);

			// Paint particle systems with update vertex buffers
			vkCmdBeginRenderPass(drawCmdBuffers[i], &renderPassBeginInfo, VK_SUBPASS_CONTENTS_INLINE);

			VkViewport viewport = vks::initializers::viewport((float)width, (float)height, 0.0f, 1.0f);
			vkCmdSetViewport(drawCmdBuffers[i], 0, 1, &viewport);

			VkRect2D scissor = vks::initializers::rect2D(width, height, 0, 0);
			vkCmdSetScissor(drawCmdBuffers[i], 0, 1, &scissor);

			VkDeviceSize offsets[1] = { 0 };

			// Render sphere
			if (sceneSetup == 0) {
				vkCmdBindPipeline(drawCmdBuffers[i], VK_PIPELINE_BIND_POINT_GRAPHICS, graphics.pipelines.sphere);
				vkCmdBindDescriptorSets(drawCmdBuffers[i], VK_PIPELINE_BIND_POINT_GRAPHICS, graphics.pipelineLayout, 0, 1, &graphics.descriptorSet, 0, NULL);
				vkCmdBindIndexBuffer(drawCmdBuffers[i], modelSphere.indices.buffer, 0, VK_INDEX_TYPE_UINT32);
				vkCmdBindVertexBuffers(drawCmdBuffers[i], 0, 1, &modelSphere.vertices.buffer, offsets);
				vkCmdDrawIndexed(drawCmdBuffers[i], modelSphere.indexCount, 1, 0, 0, 0);
			}

			// Rendering cloth
			vkCmdBindPipeline(drawCmdBuffers[i], VK_PIPELINE_BIND_POINT_GRAPHICS, graphics.pipelines.cloth);
			vkCmdBindDescriptorSets(drawCmdBuffers[i], VK_PIPELINE_BIND_POINT_GRAPHICS, graphics.pipelineLayout, 0, 1, &graphics.descriptorSet, 0, NULL);
			vkCmdBindIndexBuffer(drawCmdBuffers[i], graphics.indices.buffer, 0, VK_INDEX_TYPE_UINT32);
			vkCmdBindVertexBuffers(drawCmdBuffers[i], 0, 1, &compute.storageBuffers.output.buffer, offsets);
			vkCmdDrawIndexed(drawCmdBuffers[i], indexCount, 1, 0, 0, 0);

			drawUI(drawCmdBuffers[i]);

			vkCmdEndRenderPass(drawCmdBuffers[i]);

			// Release storage buffer to calculation queue
			addGraphicsToComputeBarriers(drawCmdBuffers[i]);

			VK_CHECK_RESULT(vkEndCommandBuffer(drawCmdBuffers[i]));
		}

	}

Run the program, as shown in the following figure (gravity only):

Or turn on wind simulation:

You can also switch to sphere collision (sceneSetup = 0), running as follows:

Tags: Spring

Posted on Sun, 14 Jun 2020 05:02:08 -0400 by Vorotaev