README.md

How-to use

1. get bgfx and follow https://bkaradzic.github.io/bgfx/build.html#quick-start for your platform

• you need the two repos bimg and bx next to bgfx
• inside bgfxrun something like ../bx/tools/bin/darwin/genie --with-combined-examples xcode10(for your platform, xcode is for mac)
• in the end you should run something like make osx-release64 (or similar for your platform)

2. copy the code in one of the examples (I used bgfx/examples/05-instancing)
3. cd bgfx/examples/platform and run from there the application in the build folder

• for mac you can use ../../.build/osx64_clang/bin/examples.app/Contents/MacOS/examplesRelease

Or use another visualization software of your choosing.. Qt is quite nice, but I had issues with Metal support on my Mac.

sph.cpp

/*
 * Copyright 2011-2020 Branimir Karadzic. All rights reserved.
 * License: https://github.com/bkaradzic/bgfx#license-bsd-2-clause
 */

#include "common.h"
#include "bgfx_utils.h"
#include "imgui/imgui.h"
#include <math.h>
#include <iostream>
#include <map>
#include <vector>

namespace
{


const uint32_t numInstances   = 800;
const float dt = 5E-3;
const float mass = 1;
const float invMass = 1/mass;
const float g = -9.81;
const float wallStiffness = 1000;
const float damping = 0.001;
const float cs = 5;
const float gamma = 1.0;
const float p0 = 0.05;
// const float gasConstant = 10;

const float initialSpacing = 0.8;
// sph constants
// kernel influence radius
const float h = 0.45 * initialSpacing;
// constant factor for kernel
const float sigma = 1/(sqrt(pow(M_PI,3))*pow(h,3));
// spacing should usually be consistent with the kernel size

struct PosColorVertex
{
	float m_x;
	float m_y;
	float m_z;
	uint32_t m_abgr;

	static void init()
	{
		ms_layout
			.begin()
			.add(bgfx::Attrib::Position, 3, bgfx::AttribType::Float)
			.add(bgfx::Attrib::Color0,   4, bgfx::AttribType::Uint8, true)
			.end();
	};

	static bgfx::VertexLayout ms_layout;
};

bgfx::VertexLayout PosColorVertex::ms_layout;

static PosColorVertex s_cubeVertices[8] =
{
	{-1.0f,  1.0f,  1.0f, 0xff000000 },
	{ 1.0f,  1.0f,  1.0f, 0xff0000ff },
	{-1.0f, -1.0f,  1.0f, 0xff00ff00 },
	{ 1.0f, -1.0f,  1.0f, 0xff00ffff },
	{-1.0f,  1.0f, -1.0f, 0xffff0000 },
	{ 1.0f,  1.0f, -1.0f, 0xffff00ff },
	{-1.0f, -1.0f, -1.0f, 0xffffff00 },
	{ 1.0f, -1.0f, -1.0f, 0xffffffff },
};

static const uint16_t s_cubeIndices[36] =
{
	0, 1, 2, // 0
	1, 3, 2,
	4, 6, 5, // 2
	5, 6, 7,
	0, 2, 4, // 4
	4, 2, 6,
	1, 5, 3, // 6
	5, 7, 3,
	0, 4, 1, // 8
	4, 5, 1,
	2, 3, 6, // 10
	6, 3, 7,
};


class ExampleSPH : public entry::AppI
{
public:
	ExampleSPH(const char* _name, const char* _description, const char* _url)
		: entry::AppI(_name, _description, _url)
	{
	}

	void init(int32_t _argc, const char* const* _argv, uint32_t _width, uint32_t _height) override
	{
		Args args(_argc, _argv);

		m_width  = _width;
		m_height = _height;
		m_debug  = BGFX_DEBUG_TEXT;
		m_reset  = BGFX_RESET_VSYNC;

		bgfx::Init init;
		init.type     = args.m_type;
		init.vendorId = args.m_pciId;
		init.resolution.width  = m_width;
		init.resolution.height = m_height;
		init.resolution.reset  = m_reset;
		bgfx::init(init);

		// Enable debug text.
		bgfx::setDebug(m_debug);

		// Set view 0 clear state.
		bgfx::setViewClear(0
			, BGFX_CLEAR_COLOR|BGFX_CLEAR_DEPTH
			, 0x303030ff
			, 1.0f
			, 0
			);

		// Create vertex stream declaration.
		PosColorVertex::init();

		// Create static vertex buffer.
		m_vbh = bgfx::createVertexBuffer(
					  bgfx::makeRef(s_cubeVertices, sizeof(s_cubeVertices) )
					, PosColorVertex::ms_layout
					);

		// Create static index buffer.
		m_ibh = bgfx::createIndexBuffer(
					bgfx::makeRef(s_cubeIndices, sizeof(s_cubeIndices) )
					);

		// Create program from shaders.
		m_program = loadProgram("vs_instancing", "fs_instancing");

		m_timeOffset = bx::getHPCounter();

		imguiCreate();

		positions = new float[numInstances*3];
		velocities = new float[numInstances*3];
		forces = new float[numInstances*3];
		densities = new float[numInstances];
		for (uint32_t i = 0; i<numInstances; i++) {
			positions[i*3] = (i % 10)*initialSpacing + 0.1;
			positions[i*3+2] = (((i / 10) % 10) * initialSpacing) + 0.1;
			positions[i*3+1] = (i / 10 / 10) * initialSpacing + 0.1;
			velocities[i*3] = 0.0;
			velocities[i*3+1] = 0.0;
			velocities[i*3+2] = 0.0;
			neighbors.push_back(std::vector<uint32_t>());
		}
	}

	int shutdown() override
	{
		imguiDestroy();

		// Cleanup.
		bgfx::destroy(m_ibh);
		bgfx::destroy(m_vbh);
		bgfx::destroy(m_program);

		// Shutdown bgfx.
		bgfx::shutdown();

		return 0;
	}

	bx::Vec3 diff(bx::Vec3 v1, bx::Vec3 v2) {
		float dx = v1.x - v2.x;
		float dy = v1.y - v2.y;
		float dz = v1.z - v2.z;
		return bx::Vec3{dx, dy, dz};
	}

	bx::Vec3 add(bx::Vec3 v1, bx::Vec3 v2) {
		float dx = v1.x + v2.x;
		float dy = v1.y + v2.y;
		float dz = v1.z + v2.z;
		return bx::Vec3{dx, dy, dz};
	}

	bx::Vec3 scale(bx::Vec3 v1, float s) {
		float dx = v1.x * s;
		float dy = v1.y * s;
		float dz = v1.z * s;
		return bx::Vec3{dx, dy, dz};
	}

	float kernel1(bx::Vec3 r) {
		double d = sqrt(r.x*r.x + r.y*r.y + r.z*r.z);
		return sigma*pow(M_E,-(pow(d,2)/pow(h,2)));
	}

	bx::Vec3 kernel1Grad(bx::Vec3 r) {
		double d = sqrt(r.x*r.x + r.y*r.y + r.z*r.z);
		return bx::Vec3{
			float(sigma*(-2*r.x/pow(h,2))*pow(M_E,-(pow(d,2)/pow(h,2)))),
			float(sigma*(-2*r.y/pow(h,2))*pow(M_E,-(pow(d,2)/pow(h,2)))),
			float(sigma*(-2*r.z/pow(h,2))*pow(M_E,-(pow(d,2)/pow(h,2))))
		};
	}

	void calculateNeighbors(uint32_t i) {
		neighbors[i].clear();
		for (uint32_t j = 0; j < numInstances; j++) {
			if (j == i) {
				continue;
			}
			bx::Vec3 p1 = bx::Vec3{positions[i*3],positions[i*3+1],positions[i*3+2]};
			bx::Vec3 p2 = bx::Vec3{positions[j*3],positions[j*3+1],positions[j*3+2]};
			bx::Vec3 r = diff(p1, p2);
			float distSq = r.x*r.x + r.y*r.y + r.z*r.z;
			if (distSq < 10 * h) {
				neighbors[i].push_back(j);
			}
		}
	}

	template<typename F>
	void forNeighbors(uint32_t i, F lambda) {
		for (uint32_t n_idx = 0; n_idx < neighbors[i].size(); n_idx++) {
			uint32_t j = neighbors[i][n_idx];
			bx::Vec3 p1 = bx::Vec3{positions[i*3],positions[i*3+1],positions[i*3+2]};
			bx::Vec3 p2 = bx::Vec3{positions[j*3],positions[j*3+1],positions[j*3+2]};
			bx::Vec3 r = diff(p1, p2);
			float distSq = r.x*r.x + r.y*r.y + r.z*r.z;
			if (distSq < 10 * h) {
				lambda(j, r);
			}
		}
	}

	float pressure(float density) {
		return ((p0*cs*cs)/gamma)*(pow(density/p0,gamma)-1);
	}

	void sph() {
		for (uint32_t i = 0; i<numInstances; i++) {
			densities[i] = 0;
			calculateNeighbors(i);
		}
		for (uint32_t i = 0; i<numInstances; i++) {
			forNeighbors(i, [i, this] (uint32_t, bx::Vec3 r) {
				this->densities[i] += mass * kernel1(r);
			});
		}
		for (uint32_t i = 0; i<numInstances; i++) {
			float pressureI = pressure(densities[i]);
			forNeighbors(i, [i, pressureI,this] (uint32_t j, bx::Vec3 r) {
				float pressureJ = pressure(densities[j]);

				float pressureFactor = ((pressureJ / (densities[j] * densities[j])) + (pressureI / (densities[i] * densities[i])));

				bx::Vec3 kernelGrad = kernel1Grad(r);

				bx::Vec3 pressureForce = scale(kernelGrad, pressureFactor * mass * mass);
				this->forces[i * 3] -= pressureForce.x;
				this->forces[i * 3 + 1] -= pressureForce.y;
				this->forces[i * 3 + 2] -= pressureForce.z;
			});
		}
	}

	void simulate() {
		// gather general forces
		for (uint32_t i = 0; i<numInstances*3; i+=3) {
			forces[i] = 0.0;
			forces[i+1] = mass * g;
			forces[i+2] = 0.0;
			if (positions[i] < 0) {
				forces[i] -= wallStiffness * positions[i];
			}
			if (positions[i+1] < 0) {
				forces[i+1] -= wallStiffness * positions[i+1];
			}
			if (positions[i+2] < 0) {
				forces[i+2] -= wallStiffness * positions[i+2];
			}
			if (positions[i] > 8) {
				forces[i] -= wallStiffness * (positions[i] - 8);
			}
			if (positions[i+2] > 16) {
				forces[i+2] -= wallStiffness * (positions[i+2] - 16);
			}
		}

		sph();

		for (uint32_t i = 0; i<numInstances*3; i++) {
			velocities[i] += forces[i]*invMass * dt;
			velocities[i] *= 1-damping;
			positions[i] += velocities[i]*dt;
		}
	}

	bool update() override
	{
		if (!entry::processEvents(m_width, m_height, m_debug, m_reset, &m_mouseState) )
		{
			simulate();

			imguiBeginFrame(m_mouseState.m_mx
				,  m_mouseState.m_my
				, (m_mouseState.m_buttons[entry::MouseButton::Left  ] ? IMGUI_MBUT_LEFT   : 0)
				| (m_mouseState.m_buttons[entry::MouseButton::Right ] ? IMGUI_MBUT_RIGHT  : 0)
				| (m_mouseState.m_buttons[entry::MouseButton::Middle] ? IMGUI_MBUT_MIDDLE : 0)
				,  m_mouseState.m_mz
				, uint16_t(m_width)
				, uint16_t(m_height)
				);

			showExampleDialog(this);

			imguiEndFrame();

			// Set view 0 default viewport.
			bgfx::setViewRect(0, 0, 0, uint16_t(m_width), uint16_t(m_height) );

			// This dummy draw call is here to make sure that view 0 is cleared
			// if no other draw calls are submitted to view 0.
			bgfx::touch(0);

			float time = (float)( (bx::getHPCounter() - m_timeOffset)/double(bx::getHPFrequency() ) );

			// Get renderer capabilities info.
			const bgfx::Caps* caps = bgfx::getCaps();

			// Check if instancing is supported.
			if (0 == (BGFX_CAPS_INSTANCING & caps->supported) )
			{
				// When instancing is not supported by GPU, implement alternative
				// code path that doesn't use instancing.
				bool blink = uint32_t(time*3.0f)&1;
				bgfx::dbgTextPrintf(0, 0, blink ? 0x4f : 0x04, " Instancing is not supported by GPU. ");
			}
			else
			{
				const bx::Vec3 at  = { 4.0f, 2.0f,   4.0f };
				const bx::Vec3 eye = { 20.0f, 30.0f, -5.0f };

				// Set view and projection matrix for view 0.
				{
					float view[16];
					bx::mtxLookAt(view, eye, at);

					float proj[16];
					bx::mtxProj(proj, 60.0f, float(m_width)/float(m_height), 0.1f, 100.0f, bgfx::getCaps()->homogeneousDepth);
					bgfx::setViewTransform(0, view, proj);

					// Set view 0 default viewport.
					bgfx::setViewRect(0, 0, 0, uint16_t(m_width), uint16_t(m_height) );
				}

				// 80 bytes stride = 64 bytes for 4x4 matrix + 16 bytes for RGBA color.
				const uint16_t instanceStride = 80;

				if (numInstances == bgfx::getAvailInstanceDataBuffer(numInstances, instanceStride) )
				{
					bgfx::InstanceDataBuffer idb;
					bgfx::allocInstanceDataBuffer(&idb, numInstances, instanceStride);

					uint8_t* data = idb.data;

					// Write instance data for 11x11 cubes.
					for (uint32_t cube = 0; cube < numInstances; ++cube)
					{
						float* mtx = (float*)data;
						bx::mtxScale(mtx, 0.4);
						mtx[12] = positions[cube*3 + 0];
						mtx[13] = positions[cube*3 + 1];
						mtx[14] = positions[cube*3 + 2];

						float* color = (float*)&data[64];
						color[0] = densities[cube]*0.01;
						color[1] = 0.1f;
						color[2] = 0.7;
						color[3] = 1.0f;

						data += instanceStride;
					}

					// Set vertex and index buffer.
					bgfx::setVertexBuffer(0, m_vbh);
					bgfx::setIndexBuffer(m_ibh);

					// Set instance data buffer.
					bgfx::setInstanceDataBuffer(&idb);

					// Set render states.
					bgfx::setState(BGFX_STATE_DEFAULT);

					// Submit primitive for rendering to view 0.
					bgfx::submit(0, m_program);
				}
			}

			// Advance to next frame. Rendering thread will be kicked to
			// process submitted rendering primitives.
			bgfx::frame();

			return true;
		}

		return false;
	}

	entry::MouseState m_mouseState;

	uint32_t m_width;
	uint32_t m_height;
	uint32_t m_debug;
	uint32_t m_reset;
	bgfx::VertexBufferHandle m_vbh;
	bgfx::IndexBufferHandle  m_ibh;
	bgfx::ProgramHandle m_program;

	float* positions;
	float* velocities;
	float* forces;
	float* densities;

	std::vector<std::vector<uint32_t>> neighbors;

	int64_t m_timeOffset;
};

} // namespace

ENTRY_IMPLEMENT_MAIN(
	  ExampleSPH
	, "05-instancing"
	, "SPH."
	, ""
	);