ex5 fluid

5_fluid/CMakeLists.txt

+
+cmake_minimum_required(VERSION 3.1)
+project(5_fluid)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/../cmake)
+set(CMAKE_CXX_FLAGS "-Wall")
+
+# libigl
+option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" OFF)
+option(LIBIGL_WITH_ANTTWEAKBAR       "Use AntTweakBar"    OFF)
+option(LIBIGL_WITH_CGAL              "Use CGAL"           OFF)
+option(LIBIGL_WITH_COMISO            "Use CoMiso"         OFF)
+option(LIBIGL_WITH_CORK              "Use Cork"           OFF)
+option(LIBIGL_WITH_EMBREE            "Use Embree"         OFF)
+option(LIBIGL_WITH_LIM               "Use LIM"            OFF)
+option(LIBIGL_WITH_MATLAB            "Use Matlab"         OFF)
+option(LIBIGL_WITH_MOSEK             "Use MOSEK"          OFF)
+option(LIBIGL_WITH_OPENGL            "Use OpenGL"         ON)
+option(LIBIGL_WITH_OPENGL_GLFW       "Use GLFW"           ON)
+option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui"          ON)
+option(LIBIGL_WITH_PNG               "Use PNG"            OFF)
+option(LIBIGL_WITH_PYTHON            "Use Python"         OFF)
+option(LIBIGL_WITH_TETGEN            "Use Tetgen"         OFF)
+option(LIBIGL_WITH_TRIANGLE          "Use Triangle"       OFF)
+option(LIBIGL_WITH_VIEWER            "Use OpenGL viewer"  ON)
+option(LIBIGL_WITH_XML               "Use XML"            OFF)
+
+if (NOT LIBIGL_FOUND)
+    find_package(LIBIGL REQUIRED QUIET)
+endif()
+
+# Add default project files
+file(GLOB LIBFILES ${PROJECT_SOURCE_DIR}/../include/*.*)
+source_group("Library Files" FILES ${LIBFILES})
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/../include)
+
+# Add your project files
+file(GLOB SRCFILES *.cpp)
+file(GLOB HFILES *.h)
+
+add_definitions(-DIGL_VIEWER_VIEWER_QUIET)
+add_executable(${PROJECT_NAME} ${SRCFILES} ${HFILES} ${LIBFILES} )
+target_link_libraries(${PROJECT_NAME} igl::core igl::opengl_glfw igl::opengl_glfw_imgui)
\ No newline at end of file

5_fluid/FluidSim.cpp

+#include "FluidSim.h"
+
+void FluidSim::solvePoisson() {
+	double dx2 = m_dx * m_dx;
+	double residual = m_acc + 1; // initial residual
+	double rho = 1;
+
+	Array2d& p = p_pressure->x();
+
+	for (int it = 0; residual > m_acc && it < m_iter; ++it) {
+		// Note that the boundaries are handles by the framework, so you iterations should be similar to:
+		for (int y = 1; y < m_res_y - 1; ++y) {
+			for (int x = 1; x < m_res_x - 1; ++x) {
+				double b = -p_divergence->x()(x, y) / m_dt * rho; // right-hand
+				// TODO: update the pressure values
+				p(x, y) = p(x, y);
+			}
+		}
+
+		// Compute the new residual, i.e. the sum of the squares of the individual residuals (squared L2-norm)
+		residual = 0;
+		for (int y = 1; y < m_res_y - 1; ++y) {
+			for (int x = 1; x < m_res_x - 1; ++x) {
+				double b = -p_divergence->x()(x, y) / m_dt * rho; // right-hand
+				// TODO: compute the cell residual
+				double cellResidual = 0;
+
+				residual += cellResidual * cellResidual;
+			}
+		}
+
+		// Get the L2-norm of the residual
+		residual = sqrt(residual);
+
+		// We assume the accuracy is meant for the average L2-norm per grid cell
+		residual /= (m_res_x - 2) * (m_res_y - 2);
+
+		// For your debugging, and ours, please add these prints after every iteration
+		if (it == m_iter - 1)
+			printf("Pressure solver: it=%d , res=%f \n", it, residual);
+		if (residual < m_acc)
+			printf("Pressure solver: it=%d , res=%f, converged \n", it, residual);
+	}
+}
+
+void FluidSim::correctVelocity() {
+	Array2d& p = p_pressure->x();
+	Array2d& u = p_velocity->x();	// x velocity
+	Array2d& v = p_velocity->y();	// y velocity
+
+	// Note: velocity u_{i+1/2}
+	for (int y = 1; y < m_res_y - 1; ++y)
+		for (int x = 1; x < m_res_x; ++x)
+			// TODO: update u
+			u(x, y) = u(x, y);
+
+	// Same for velocity v_{i+1/2}.
+	for (int y = 1; y < m_res_y; ++y)
+		for (int x = 1; x < m_res_x - 1; ++x)
+			// TODO: update v
+			v(x, y) = v(x, y);
+}
+
+void FluidSim::advectValues() {
+	// Densities live on the grid centers, the velocities on the MAC grid
+	// Separate their computation to avoid confusion
+
+	Array2d& d = p_density->x();
+	Array2d& u = p_velocity->x();
+	Array2d& v = p_velocity->y();
+
+	Array2d& d_ = p_density_tmp->x();
+	Array2d& u_ = p_velocity_tmp->x();
+	Array2d& v_ = p_velocity_tmp->y();
+
+	// Densities, grid centers
+	for (int y = 1; y < m_res_y - 1; ++y) {
+		for (int x = 1; x < m_res_x - 1; ++x) {
+			// TODO: Compute the velocity
+			double last_x_velocity = 0;
+			double last_y_velocity = 0;
+
+			// TODO: Find the last position of the particle (in grid coordinates)
+			double last_x = 0;
+			double last_y = 0;
+
+			// Make sure the coordinates are inside the boundaries
+			// Densities are known between 1 and res-2
+			if (last_x < 1) last_x = 1;
+			if (last_y < 1) last_y = 1;
+			if (last_x > m_res_x - 2) last_x = m_res_x - 2;
+			if (last_y > m_res_y - 2) last_y = m_res_y - 2;
+
+			// Determine corners for bilinear interpolation
+			int x_low = (int)last_x;
+			int y_low = (int)last_y;
+			int x_high = x_low + 1;
+			int y_high = y_low + 1;
+
+			// Compute the interpolation weights
+			double x_weight = last_x - x_low;
+			double y_weight = last_y - y_low;
+
+			// TODO: Bilinear interpolation
+			d_(x, y) = d(x, y);
+		}
+	}
+
+	// Velocities (u), MAC grid
+	for (int y = 1; y < m_res_y - 1; ++y) {
+		for (int x = 1; x < m_res_x; ++x) {
+			// TODO: Compute the velocity
+			double last_x_velocity = 0;
+			double last_y_velocity = 0;
+
+			// TODO: Find the last position of the particle (in grid coordinates)
+			double last_x = 0;
+			double last_y = 0;
+
+			// Make sure the coordinates are inside the boundaries
+			// Being conservative, one can say that the velocities are known between 1.5 and res-2.5
+			// (the MAC grid is inside the known densities, which are between 1 and res - 2)
+			if (last_x < 1.5) last_x = 1.5;
+			if (last_y < 1.5) last_y = 1.5;
+			if (last_x > m_res_x - 1.5) last_x = m_res_x - 1.5;
+			if (last_y > m_res_y - 2.5) last_y = m_res_y - 2.5;
+
+			// Determine corners for bilinear interpolation
+			int x_low = (int)last_x;
+			int y_low = (int)last_y;
+			int x_high = x_low + 1;
+			int y_high = y_low + 1;
+
+			// Compute the interpolation weights
+			double x_weight = last_x - x_low;
+			double y_weight = last_y - y_low;
+
+			// TODO: Bilinear interpolation
+			u_(x, y) = u(x, y);
+		}
+	}
+
+	// Velocities (v), MAC grid
+	for (int y = 1; y < m_res_y; ++y) {
+		for (int x = 1; x < m_res_x - 1; ++x) {
+			// TODO: Compute the velocity
+			double last_x_velocity = 0;
+			double last_y_velocity = 0;
+
+			// TODO: Find the last position of the particle (in grid coordinates)
+			double last_x = 0;
+			double last_y = 0;
+
+			// Make sure the coordinates are inside the boundaries
+			// Being conservative, one can say that the velocities are known between 1.5 and res-2.5
+			// (the MAC grid is inside the known densities, which are between 1 and res - 2)
+			if (last_x < 1.5) last_x = 1.5;
+			if (last_y < 1.5) last_y = 1.5;
+			if (last_x > m_res_x - 2.5) last_x = m_res_x - 2.5;
+			if (last_y > m_res_y - 1.5) last_y = m_res_y - 1.5;
+
+			// Determine corners for bilinear interpolation
+			double x_low = (int)last_x;
+			double y_low = (int)last_y;
+			double x_high = x_low + 1;
+			double y_high = y_low + 1;
+
+			// Compute the interpolation weights
+			double x_weight = last_x - x_low;
+			double y_weight = last_y - y_low;
+
+			// TODO: Bilinear interpolation
+			v_(x, y) = v(x, y);
+		}
+	}
+
+	// Copy the values in temp to the original buffers
+	for (int y = 1; y < m_res_y - 1; ++y)
+		for (int x = 1; x < m_res_x - 1; ++x)
+			d(x, y) = d_(x, y);
+	for (int y = 1; y < m_res_y - 1; ++y)
+		for (int x = 1; x < m_res_x; ++x)
+			u(x, y) = u_(x, y);
+	for (int y = 1; y < m_res_y; ++y)
+		for (int x = 1; x < m_res_x - 1; ++x)
+			v(x, y) = v_(x, y);
+}
+

5_fluid/FluidSim.h

+#include <igl/edges.h>
+#include "Simulation.h"
+#include "Grid2.h"
+#include "MACGrid2.h"
+
+using namespace std;
+
+/*
+ * Simulation of a simple smoke plume rising.
+ */
+class FluidSim : public Simulation {
+public:
+	FluidSim() : Simulation() { init(); }
+
+	virtual void init() override {
+		m_res_x = 128;
+		m_res_y = int(m_res_x*1.5); // 3:2 ratio
+		m_size_x = 1.0; //m_res_x; // or just 1.0
+		m_dx = m_size_x / m_res_x; // ! dx == dy
+		m_idx = m_res_x / m_size_x;
+		m_size_y = m_dx * m_res_y;
+		m_dt = 0.005 * sqrt((m_res_x + m_res_y) * 0.5);
+		m_acc = 1e-5;
+		m_iter = 2000;
+		m_field = 0;
+		m_velocityOn = false;
+		m_vScale = 20;
+		m_windOn = false;
+
+		p_density = new Grid2(m_res_x, m_res_y, m_dx);
+		p_density_tmp = new Grid2(m_res_x, m_res_y, m_dx);
+		p_pressure = new Grid2(m_res_x, m_res_y, m_dx);
+		p_divergence = new Grid2(m_res_x, m_res_y, m_dx);
+		p_vorticity = new Grid2(m_res_x, m_res_y, m_dx);
+		p_density->getMesh(m_renderV, m_renderF); // need to call once
+		
+		p_velocity = new MACGrid2(m_res_x, m_res_y, m_dx);
+		p_velocity_tmp = new MACGrid2(m_res_x, m_res_y, m_dx);
+		p_force = new MACGrid2(m_res_x, m_res_y, m_dx);
+		reset();
+	}
+
+	virtual void resetMembers() override {
+		p_density->reset();
+		p_density->applySource(0.45, 0.55, 0.1, 0.15);
+		p_pressure->reset();
+		p_divergence->reset();
+		p_velocity->reset();
+		p_force->reset();
+	}
+
+	virtual void updateRenderGeometry() override {
+		if (m_field == 0) {
+			p_density->getColors(m_renderC);
+		}
+		else if (m_field == 1) {
+			p_pressure->getColors(m_renderC, true);
+		}
+		else if (m_field == 2) {
+			p_divergence->getColors(m_renderC, true);
+		}
+		else if (m_field == 3) {
+			p_vorticity->getColors(m_renderC, true);
+		}
+		
+		if (m_velocityOn) {
+			p_velocity->updateEdges(m_vScale);
+		}
+	}
+
+	virtual bool advance() override {
+		// apply source in density field
+		p_density->applySource(0.45, 0.55, 0.1, 0.15);
+
+		// add in new forces
+		addBuoyancy();
+		if (m_windOn)
+			addWind();
+		addForce();
+
+		// remove divergence
+		solvePressure();
+
+		// advect everything
+		advectValues();
+
+		// reset forces
+		p_force->reset();
+
+		// advance m_time
+		m_time += m_dt;
+		m_step++;
+
+		return false;
+	}
+
+	virtual void renderRenderGeometry(
+		igl::opengl::glfw::Viewer& viewer) override {
+		viewer.data().set_mesh(m_renderV, m_renderF);
+		viewer.data().set_colors(m_renderC);
+
+		if (m_velocityOn) {
+			viewer.data().add_edges(p_velocity->s(), p_velocity->e(), Eigen::RowVector3d(0, 0, 0));
+			viewer.data().add_edges(p_velocity->vs(), p_velocity->ve(), p_velocity->vc());
+		}
+	}	
+#pragma region FluidSteps
+	void addBuoyancy() {
+		double scaling = 64.0 / m_res_x;
+
+		// add buoyancy
+		for (int i = 0; i < p_force->y().size(0); ++i) {
+			for (int j = 1; j < p_force->y().size(1) - 1; ++j) {
+				p_force->y()(i, j) += 0.1 * (p_density->x()(i, j - 1) + p_density->x()(i, j)) / 2.0 * scaling;
+			}
+		}
+	}
+
+	void addWind() {
+		double scaling = 64.0 / m_res_x;
+
+		static double r = 0.0;
+		r += 1;
+
+		const double fx = 2e-2 * cos(5e-2 * r) * cos(3e-2 * r) * scaling;
+
+		// add wind
+		for (int i = 0; i < p_force->x().size(0); ++i) {
+			for (int j = 0; j < p_force->x().size(1); ++j) {
+				p_force->x()(i, j) += fx;
+			}
+		}
+	}
+
+	void addForce() {
+		for (int i = 0; i < p_velocity->x().size(0); ++i) {
+			for (int j = 0; j < p_velocity->x().size(1); ++j) {
+				p_velocity->x()(i, j) += m_dt * p_force->x()(i, j);
+			}
+		}
+
+		for (int i = 0; i < p_velocity->y().size(0); ++i) {
+			for (int j = 0; j < p_velocity->y().size(1); ++j) {
+				p_velocity->y()(i, j) += m_dt * p_force->y()(i, j);
+			}
+		}
+	}
+	
+	void solvePressure() {
+		// copy out the boundaries 
+		setNeumann();
+		setZero();
+
+		computeDivergence();
+
+		// solve Poisson equation
+		copyBorder();
+		solvePoisson();
+
+		correctVelocity();
+
+		computeVorticity();
+
+		// for debugging
+		computeDivergence();
+	}
+
+	void setNeumann() {
+		// x-velocity
+		Array2d& u = p_velocity->x();
+		int sx = u.size(0);
+		int sy = u.size(1);
+		for (int y = 0; y < sy; ++y) {
+			u(0, y) = u(2, y);
+			u(sx - 1, y) = u(sx - 3, y);
+		}
+
+		// y-velocity
+		Array2d& v = p_velocity->y();
+		sx = v.size(0);
+		sy = v.size(1);
+		for (int x = 0; x < sx; ++x) {
+			v(x, 0) = v(x, 2);
+			v(x, sy - 1) = v(x, sy - 3);
+		}
+	}
+
+	void setZero() {
+		// x-velocity
+		Array2d& u = p_velocity->x();
+		int sx = u.size(0);
+		int sy = u.size(1);
+		for (int x = 0; x < sx; ++x) {
+			u(x, 0) = 0;
+			u(x, sy - 1) = 0;
+		}
+
+		// y-velocity
+		Array2d& v = p_velocity->y();
+		sx = v.size(0);
+		sy = v.size(1);
+		for (int y = 0; y < sy; ++y) {
+			v(0, y) = 0;
+			v(sx - 1, y) = 0;
+		}
+	}
+
+	void computeDivergence() {
+		// calculate divergence
+		for (int y = 1; y < m_res_y - 1; ++y) {
+			for (int x = 1; x < m_res_x - 1; ++x) {
+				double xComponent = (p_velocity->x()(x + 1, y) - p_velocity->x()(x, y)) * m_idx;
+				double yComponent = (p_velocity->y()(x, y + 1) - p_velocity->y()(x, y)) * m_idx;
+				p_divergence->x()(x, y) = xComponent + yComponent;
+			}
+		}
+	}
+
+	void computeVorticity() {
+		// calculate vorticity
+		for (int y = 1; y < m_res_y - 1; ++y) {
+			for (int x = 1; x < m_res_x - 1; ++x) {
+				double xComponent = (p_velocity->x()(x + 1, y) - p_velocity->x()(x, y)) * m_idx;
+				double yComponent = (p_velocity->y()(x, y + 1) - p_velocity->y()(x, y)) * m_idx;
+				p_vorticity->x()(x, y) = yComponent - xComponent;
+			}
+		}
+	}
+
+	void copyBorder() {
+		Array2d& p = p_pressure->x();
+		int sx = p.size(0);
+		int sy = p.size(1);
+		for (int y = 0; y < sy; ++y) {
+			p(0, y) = p(1, y);
+			p(sx - 1, y) = p(sx - 2, y);
+		}
+		for (int x = 0; x < sx; ++x) {
+			p(x, 0) = p(x, 1);
+			p(x, sy - 1) = p(x, sy - 2);
+		}
+	}
+
+#pragma endregion FluidSteps
+
+#pragma region Exercise
+	void solvePoisson();
+	void correctVelocity();
+	void advectValues();
+#pragma endregion Exercise
+
+#pragma region SettersAndGetters
+	void selectField(int field) { m_field = field; }
+	void selectVField(bool v) { m_velocityOn = v; }
+	void setVelocityScale(double s) { m_vScale = s; }
+	void setResX(int r) { m_res_x = r; }
+	void setResY(int r) { m_res_y = r; }
+	void setAccuracy(double acc) { m_acc = acc; }
+	void setIteration(int iter) { m_iter = iter; }
+	void setWind(bool w) { m_windOn = w; }
+
+	int getField() const { return m_field; }
+	bool getVField() const { return m_velocityOn; }
+	double getVelocityScale() const { return m_vScale; }
+	int getResX() const { return m_res_x; }
+	int getResY() const { return m_res_y; }
+	double getAccuracy() const { return m_acc; }
+	int getIteration() const { return m_iter; }
+	bool getWind() const { return m_windOn; }
+	double getTimestep() const { return m_dt; }
+#pragma endregion SettersAndGetters
+
+private:
+	int m_res_x, m_res_y;
+	double m_dx, m_idx;			// dx, inverse dx
+	double m_size_x, m_size_y;
+	double m_acc;				// solver accuracy
+	int m_iter;
+	int m_field;
+	bool m_velocityOn;
+	double m_vScale;
+	bool m_windOn;
+
+	Grid2* p_density;
+	Grid2* p_density_tmp;
+	Grid2* p_pressure;
+	Grid2* p_divergence;
+	Grid2* p_vorticity;
+	MACGrid2* p_velocity;
+	MACGrid2* p_velocity_tmp;
+	MACGrid2* p_force;
+
+	Eigen::MatrixXd m_renderV; // vertex positions, 
+	Eigen::MatrixXi m_renderF; // face indices 
+	Eigen::MatrixXd m_renderC; // face (or vertex) colors for rendering
+};
\ No newline at end of file

5_fluid/main.cpp

+#include <igl/writeOFF.h>
+#include <thread>
+#include "Gui.h"
+#include "Simulator.h"
+#include "FluidSim.h"
+
+/*
+ */
+class FluidGui : public Gui {
+public:
+	float m_dt;
+	float m_acc;
+	int m_iter;
+	float m_vScale;
+	bool m_windOn;
+
+	const vector<char const*> m_fields = {
+	   "Density", "Pressure", "Divergence", "Vorticity"
+	};
+	int m_selected_field;
+	bool m_velocityOn;
+
+	FluidSim* p_fluidSim = NULL;
+
+	FluidGui() {
+		turnOffLight(); // no light for field visualization
+
+		p_fluidSim = new FluidSim();
+		m_dt = p_fluidSim->getTimestep();
+		m_acc = p_fluidSim->getAccuracy();
+		m_iter = p_fluidSim->getIteration();
+		m_vScale = p_fluidSim->getVelocityScale();
+		m_windOn = p_fluidSim->getWind();
+		m_selected_field = p_fluidSim->getField();
+		m_velocityOn = p_fluidSim->getVField();
+		setSimulation(p_fluidSim);
+		start();
+	}
+
+	virtual void updateSimulationParameters() override {
+		// change all parameters of the simulation to the values that are set in
+		// the GUI
+		p_fluidSim->setVelocityScale(m_vScale);
+		p_fluidSim->setTimestep(m_dt);
+		p_fluidSim->setAccuracy(m_acc);
+		p_fluidSim->setIteration(m_iter);
+	}
+
+	virtual void clearSimulation() override {
+		p_fluidSim->setWind(m_windOn);
+	}
+
+	virtual void drawSimulationParameterMenu() override {
+		if (ImGui::Combo("Fields", &m_selected_field,
+			m_fields.data(), m_fields.size())) {
+			cout << m_selected_field << endl;
+			p_fluidSim->selectField(m_selected_field);
+			p_fluidSim->updateRenderGeometry();
+		}
+		if (ImGui::Checkbox("Velocity", &m_velocityOn)) {
+			p_fluidSim->selectVField(m_velocityOn);
+			p_fluidSim->updateRenderGeometry();
+		}
+		ImGui::InputFloat("v scale", &m_vScale, 0, 0);
+		ImGui::InputFloat("dt", &m_dt, 0, 0);
+		ImGui::InputFloat("accuracy", &m_acc, 0, 0, 5);
+		ImGui::InputInt("iter", &m_iter, 0, 0);
+		if (ImGui::Checkbox("Wind", &m_windOn))
+			p_fluidSim->setWind(m_windOn);
+	}
+};
+
+int main(int argc, char* argv[]) {
+	// create a new instance of the GUI for the spring simulation
+	new FluidGui();
+
+	return 0;
+}
\ No newline at end of file

CMakeLists.txt

@@ -34,5 +34,6 @@ endif()
 # add_subdirectory(2-1_bead)
 # add_subdirectory(2-2_pendulum)
 # add_subdirectory(3_cloth)
-add_subdirectory(4-1_spinning)
-add_subdirectory(4-2_collision)
\ No newline at end of file
+# add_subdirectory(4-1_spinning)
+# add_subdirectory(4-2_collision)
+add_subdirectory(5_fluid)
\ No newline at end of file