From 2e152217862e2156e5d73284dcb0fde6cf0b56fe Mon Sep 17 00:00:00 2001
From: Yuwei Xiao <xyw1105@126.com>
Date: Wed, 22 Apr 2020 10:05:44 +0800
Subject: [PATCH] ex5 fluid
---
5_fluid/CMakeLists.txt | 43 ++++++
5_fluid/FluidSim.cpp | 188 ++++++++++++++++++++++++++
5_fluid/FluidSim.h | 296 +++++++++++++++++++++++++++++++++++++++++
5_fluid/main.cpp | 78 +++++++++++
CMakeLists.txt | 5 +-
5 files changed, 608 insertions(+), 2 deletions(-)
create mode 100644 5_fluid/CMakeLists.txt
create mode 100644 5_fluid/FluidSim.cpp
create mode 100644 5_fluid/FluidSim.h
create mode 100644 5_fluid/main.cpp
diff --git a/5_fluid/CMakeLists.txt b/5_fluid/CMakeLists.txt
new file mode 100644
index 0000000..5a261ea
--- /dev/null
+++ b/5_fluid/CMakeLists.txt
@@ -0,0 +1,43 @@
+
+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
diff --git a/5_fluid/FluidSim.cpp b/5_fluid/FluidSim.cpp
new file mode 100644
index 0000000..c535541
--- /dev/null
+++ b/5_fluid/FluidSim.cpp
@@ -0,0 +1,188 @@
+#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);
+}
+
diff --git a/5_fluid/FluidSim.h b/5_fluid/FluidSim.h
new file mode 100644
index 0000000..24d2387
--- /dev/null
+++ b/5_fluid/FluidSim.h
@@ -0,0 +1,296 @@
+#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
diff --git a/5_fluid/main.cpp b/5_fluid/main.cpp
new file mode 100644
index 0000000..087e23d
--- /dev/null
+++ b/5_fluid/main.cpp
@@ -0,0 +1,78 @@
+#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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fa70829..1743953 100644
--- a/CMakeLists.txt
+++ b/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
--
GitLab