main.cpp

#include "Constants.h"
#include "MPMGrid.h"
#include "Particle.h"
#include <CircleSource.h>
#include <GlyphPolyDataMapper.h>
#include <ImageData.h>
#include <ImageMapper.h>
#include <PhongMaterial.h>
#include <PlaneSource.h>
#include <PolyData.h>
#include <PolyDataMapper.h>
#include <PolyDataPointCloud.h>
#include <Renderer.h>
#include <RenderWindow.h>
#include <SphereSource.h>
#include <TrackballCamera.h>
#include <TrackballCameraInteractor.h>

// Updates image with mass values from mpm grid
void updateGridImage(MPMGrid* mpmGrid, ImageMapper* mapper);
// Updates particle poly data with scalar values from particles
void updateParticlePoly(MPMGrid* mpmGrid, GlyphPolyDataMapper* mapper);

int main(int argc, char *argv[])
{
	// Create the window
	// This has to happen before any gl calls in other objects because 
	// glfw can only make the opengl context when creating the window.
	RenderWindow renWindow("MPM Simulation");

	// Create the camera for the renderer to use
	TrackballCamera cam(1.4f, 1.57f, 30.0f);

	// Create the renderer
	Renderer ren;
	ren.setCamera(&cam);
	ren.addMaterial(PhongMaterial(glm::vec3(0.2f, 0.4f, 0.2f), 0.5f));
	ren.addMaterial(PhongMaterial(glm::vec3(0.8f, 0.2f, 0.2f), 1.0f));
	renWindow.setRenderer(&ren);

	// Setup the camera interactor (maps user window input to camera)
	TrackballCameraInteractor iren;
	iren.setCamera(&cam);
	renWindow.setInteractor(&iren);

	// Setup a ground plane
	PlaneSource plane;
	plane.update();

	PolyDataMapper planeMapper;
	planeMapper.setInput(plane.getOutput());
	planeMapper.setMaterial(ren.getMaterial(0));
	planeMapper.setModelMatrix(MathHelp::translate(0.0f, -2.0f, 0.0f) * MathHelp::scale(10.0f));
	planeMapper.update();
	ren.addRenderItem(&planeMapper);

	// Generate point cloud within circle
	CircleSource particleShape;
	particleShape.setRadius(0.5f);
	particleShape.setDivisions(25);
	particleShape.update();
	const UINT particleCount = static_cast<UINT>(particleShape.getOutput()->getArea() / (PARTICLE_DIAMETER * PARTICLE_DIAMETER));
	PolyDataPointCloud polyDataPtCloud;
	polyDataPtCloud.setUse2d(true);
	polyDataPtCloud.setOptimizeByRadius(true);
	polyDataPtCloud.setRadius(0.05f);
	polyDataPtCloud.setNumberOfIterations(30);
	polyDataPtCloud.setNumberOfPoints(particleCount);
	polyDataPtCloud.setInput(particleShape.getOutput());
	polyDataPtCloud.update();
	PolyData* ptCloudPolyData = polyDataPtCloud.getOutput();

	// Create a uv sphere source for instancing
	SphereSource particleSphereSource;
	particleSphereSource.setRadius(0.02f);
	particleSphereSource.update();
	// Create the particle mapper
	GlyphPolyDataMapper particleMapper;
	particleMapper.setInput(particleSphereSource.getOutput());
	particleMapper.allocateOffsets(particleCount);
	particleMapper.allocateColorData(particleCount);
	glm::vec3* offsetData = reinterpret_cast<glm::vec3*>(particleMapper.getOffsetData());
	// Setup the polydata and particles
	glm::vec3* posData = reinterpret_cast<glm::vec3*>(ptCloudPolyData->getVertexData());
	Particle* particles = new Particle[particleCount];
	for (UINT i = 0; i < particleCount; i++)
	{
		offsetData[i] = posData[i];

		particles[i].pos = &offsetData[i];
		particles[i].mass = PARTICLE_MASS;
		particles[i].bulk = BULK_MODULUS;
		particles[i].shear = SHEAR_MODULUS;
	}

	ren.addRenderItem(&particleMapper);

	// Setup the MPMGrid for simulation
	MPMGrid mpmGrid;
	GLfloat padScale = 2.0f;
	GLfloat* bounds = polyDataPtCloud.getBounds();
	glm::vec2 padSize = glm::vec2(bounds[1] - bounds[0], bounds[3] - bounds[2]) * padScale;
	glm::vec2 boundsCenter = glm::vec2(bounds[1] + bounds[0], bounds[3] + bounds[2]) * 0.5f;
	glm::vec2 origin = boundsCenter - padSize * 0.5f;
	mpmGrid.initGrid(origin, padSize, GRID_DIM, GRID_DIM);
	mpmGrid.initParticles(particles, particleCount);

	updateParticlePoly(&mpmGrid, &particleMapper);

	// Setup a background image for visualizing the node values
	ImageMapper imageMapper;
	imageMapper.setModelMatrix(MathHelp::translate(0.0f, boundsCenter.y, -0.1f));
	ren.addRenderItem(&imageMapper);

	// Update loop
	UINT frameCount = 0;
	float frameTime = 0.0f;
	while (renWindow.isActive())
	{
#ifdef STATS
		printf("Frame: %d\n", frameCount);
		auto start = std::chrono::steady_clock::now();
#endif
		// Do the actual simulation
		for (UINT i = 0; i < SUBSTEPS; i++)
		{
			mpmGrid.projectToGrid();
			// Split so user can apply their own forces to the velocities
			mpmGrid.update(TIMESTEP);
		}
		
#ifdef STATS
		frameTime += TIMESTEP * SUBSTEPS;
		printf("PostFrameTime: %f\n", frameTime);
		auto end = std::chrono::steady_clock::now();
		printf("Sim Time: %fs\n", std::chrono::duration<double, std::milli>(end - start).count() / 1000.0);
#endif

		updateGridImage(&mpmGrid, &imageMapper);
		updateParticlePoly(&mpmGrid, &particleMapper);
		renWindow.render();
#ifdef OUTPUTFRAMES
		// Get the frame
		GLint vp[4];
		glGetIntegerv(GL_VIEWPORT, vp);
		ImageData image;
		UINT dim[3] = { static_cast<UINT>(vp[2]), static_cast<UINT>(vp[3]), 1 };
		double spacing[3] = { 1.0, 1.0, 1.0 };
		double origin[3] = { 0.0, 0.0, 0.0 };
		image.allocate2DImage(dim, spacing, origin, 3, ScalarType::UCHAR_T);
		glReadPixels(0, 0, dim[0], dim[1], GL_RGB, GL_UNSIGNED_BYTE, image.getData());
		// Write the frame as png
		PNGWriter writer;
		if (frameCount < 10)
			writer.setFileName("output/frame_000" + std::to_string(frameCount) + ".png");
		else if (frameCount < 100)
			writer.setFileName("output/frame_00" + std::to_string(frameCount) + ".png");
		else if (frameCount < 1000)
			writer.setFileName("output/frame_0" + std::to_string(frameCount) + ".png");
		else if (frameCount < 10000)
			writer.setFileName("output/frame_" + std::to_string(frameCount) + ".png");
		writer.setInput(&image);
		writer.update();
#endif
		frameCount++;
	}

	delete[] particles;

	return 1;
}

// Puts the mass values from the mpmgrid into an image
void updateGridImage(MPMGrid* mpmGrid, ImageMapper* mapper)
{
	// Extract the mass into an image
	ImageData* imageData = new ImageData();
	UINT dim[3] = { static_cast<UINT>(mpmGrid->gridWidth), static_cast<UINT>(mpmGrid->gridHeight), 1 };
	double spacing[3] = { mpmGrid->cellSize.x, mpmGrid->cellSize.y, 0.0 };
	double origin[3] = { 0.0, 0.0, 0.0 };
	imageData->allocate2DImage(dim, spacing, origin, 1, ScalarType::UCHAR_T);
	unsigned char* data = static_cast<unsigned char*>(imageData->getData());

	// Get the max mass
	GLfloat max = 0.0f;
	for (int i = 0; i < mpmGrid->gridHeight * mpmGrid->gridWidth; i++)
	{
		if (mpmGrid->nodes[i].mass > max)
			max = mpmGrid->nodes[i].mass;
	}
	// Set the image values
	GLfloat ratio = 255.0f / max;
	for (int i = 0; i < mpmGrid->gridHeight * mpmGrid->gridWidth; i++)
	{
		data[i] = static_cast<unsigned char>(mpmGrid->nodes[i].mass * ratio);
	}

	mapper->setInput(imageData);
	mapper->update();
}

// Puts the defGp values into color values for the poly
void updateParticlePoly(MPMGrid* mpmGrid, GlyphPolyDataMapper* mapper)
{
	glm::vec3* colors = reinterpret_cast<glm::vec3*>(mapper->getColorData());
	GLfloat max = 0.0f;
	for (UINT i = 0; i < mapper->getInstanceCount(); i++)
	{
		colors[i].y = glm::determinant(mpmGrid->particles[i].defGp);
		if (colors[i].y > max)
			max = colors[i].y;
		colors[i].x = 0.0f;
		colors[i].z = 0.0f;
	}
	for (UINT i = 0; i < mapper->getInstanceCount(); i++)
	{
		colors[i].y *= 1.0f / max;
	}
	mapper->update();
}