initialize

2-2_pendulum/PendulumSim.cpp

+#include "PendulumSim.h"
+
+typedef Eigen::Matrix<double, 6, 1> Vector6d;
+typedef Eigen::Matrix<double, 2, 6> Matrix26d;
+
+bool PendulumSim::advance() {
+    auto pos1 = p1->getPosition();
+    auto pos2 = p2->getPosition();
+    auto vel1 = p1->getLinearVelocity();
+    auto vel2 = p2->getLinearVelocity();
+
+    // TODO update positions and velocities of particle p1, p2
+    // c1 = 0.5 * (x1.dot(x1) - l1^2) = 0
+    // c2 = 0.5 * ((xi-x2).dot(xi-x2) - l2^2) = 0
+
+
+
+    // advance m_time
+    m_time += m_dt;
+    m_step++;
+
+    // log
+    if ((m_step % m_log_frequency) == 0) {
+        m_trajectories[0].push_back(p1->getPosition());
+        m_trajectories[1].push_back(p2->getPosition());
+    }
+
+    return false;
+}
\ No newline at end of file

2-2_pendulum/PendulumSim.h

+#include "Simulation.h"
+#include "igl/PI.h"
+
+using namespace std;
+
+struct DoublePendulum {
+public:
+    double m1, m2;
+    double l1, l2;
+};
+
+/*
+ */
+class PendulumSim : public Simulation {
+public:
+    PendulumSim() : Simulation() {
+        init();
+        m_trajectories.clear();
+        m_trajectoryColors.clear();
+    }
+
+    virtual void init() override {
+        std::string path = "sphere.off";
+        m_objects.clear();
+        m_objects.push_back(RigidObject(path));
+        m_objects.push_back(RigidObject(path));
+        p1 = &m_objects[0];
+        p2 = &m_objects[1];
+
+        m_dt = 1e-2;
+        m_log_frequency = 30;
+        m_gravity << 0, -9.81, 0;
+        // visualization color for sphere and trajectories
+        m_color1 << 1.0, 0.5, 0;
+        m_color2 << 0.0, 0.5, 1.0;
+        reset();
+    }
+
+    virtual void resetMembers() override {
+        p1->reset(); p2->reset();
+        p1->setScale(0.04); p2->setScale(0.04);
+        p1->setPosition(Eigen::Vector3d(m_pendulum.l1*cos(igl::PI/4), m_pendulum.l1*sin(igl::PI/4.0), 0));
+        p2->setPosition(Eigen::Vector3d((m_pendulum.l1+m_pendulum.l2)*cos(igl::PI/4), (m_pendulum.l1+m_pendulum.l2)*sin(igl::PI/4.0), 0));
+        p1->setColors(m_color1);
+        p2->setColors(m_color2);
+
+        clearTrajectories();
+    }
+
+    virtual void updateRenderGeometry() override {
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {    // negative id means newly created object, reverse memory for it
+                m_renderVs.emplace_back();
+                m_renderFs.emplace_back();
+            }
+
+            m_objects[i].getMesh(m_renderVs[i], m_renderFs[i]);
+        }
+    }
+
+	virtual bool advance() override;
+
+    virtual void renderRenderGeometry(igl::opengl::glfw::Viewer &viewer) override {
+        // render point
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {
+                int new_id = 0;
+                if (i > 0) {
+                    new_id = viewer.append_mesh();
+                    o.setID(new_id);
+                } else {
+                    o.setID(new_id);
+                }
+
+                size_t meshIndex = viewer.mesh_index(o.getID());
+                viewer.data_list[meshIndex].set_face_based(true);
+                viewer.data_list[meshIndex].point_size = 2.0f;
+                viewer.data_list[meshIndex].clear();
+            }
+            size_t meshIndex = viewer.mesh_index(o.getID());
+
+            viewer.data_list[meshIndex].set_mesh(m_renderVs[i], m_renderFs[i]);
+            viewer.data_list[meshIndex].compute_normals();
+
+            Eigen::MatrixXd color;
+            o.getColors(color);
+            viewer.data_list[meshIndex].set_colors(color);
+        }
+
+        for (size_t trajectory = 0; trajectory < m_trajectories.size(); trajectory++) {
+            if(m_trajectories[trajectory].size() <= 1) 
+                continue;
+            Eigen::MatrixXd p1s, p2s, cs;
+            p1s.resize(m_trajectories[trajectory].size()-1, 3);
+            p2s.resize(m_trajectories[trajectory].size()-1, 3);
+            cs.resize(m_trajectories[trajectory].size()-1, 3);
+            for (size_t point = 1; point < m_trajectories[trajectory].size(); point++) {
+                p1s.row(point-1) << m_trajectories[trajectory][point-1].transpose();
+                p2s.row(point-1) << m_trajectories[trajectory][point].transpose();
+                cs.row(point-1) << m_trajectoryColors[trajectory];
+                // viewer.data().add_edges(m_trajectories[trajectory][point-1].transpose(),
+                //         m_trajectories[trajectory][point].transpose(),
+                //         m_trajectoryColors[trajectory]);
+            }
+            viewer.data(0).add_edges(p1s, p2s, cs);
+        }
+
+		// draw the pendulum
+        viewer.data(0).add_edges(Eigen::RowVector3d::Zero(), p1->getPosition().transpose(), m_color);
+        viewer.data(0).add_edges(p1->getPosition().transpose(), p2->getPosition().transpose(), m_color);
+    }
+
+    void clearTrajectories() {
+        m_trajectories.clear();
+        m_trajectories.push_back(vector<Eigen::Vector3d>());    // for p1
+        m_trajectories.push_back(vector<Eigen::Vector3d>());    // for p2
+        m_trajectoryColors.clear();
+        m_trajectoryColors.push_back(m_color1);                  // p1 color
+        m_trajectoryColors.push_back(m_color2);                  // p1 color
+    }
+
+    void setPendulum(DoublePendulum &s) { m_pendulum = s; }
+    void setLogFrequency(int f) { m_log_frequency = f; }
+
+    void getTrajectories(int index, Eigen::MatrixX3d &mat) const {
+        int num_points = m_trajectories[index].size();
+        mat.resize(num_points, 3);
+        for (int i = 0; i < num_points; i++) {
+            mat.row(i) = m_trajectories[index][i];
+        }
+    }
+    int getNumTrajectories() const { return m_trajectories.size(); }
+
+private:
+
+    Eigen::Vector3d m_gravity;
+
+    DoublePendulum m_pendulum;
+    RigidObject *p1, *p2;
+
+    std::vector<Eigen::MatrixXd> m_renderVs;  // vertex positions for rendering
+    std::vector<Eigen::MatrixXi> m_renderFs;  // face indices for rendering
+
+    int m_log_frequency;  // how often should we log the COM in the GUI
+    vector<vector<Eigen::Vector3d>> m_trajectories;
+    Eigen::RowVector3d m_color, m_color1, m_color2;
+    vector<Eigen::RowVector3d> m_trajectoryColors;
+};
\ No newline at end of file

2-2_pendulum/main.cpp

+#include <igl/writeOFF.h>
+#include <thread>
+#include "Gui.h"
+#include "Simulator.h"
+#include "PendulumSim.h"
+
+/*
+ */
+class PendulumGui : public Gui {
+   public:
+    // simulation parameters
+    float m_dt = 1e-2;
+    int m_log_frequency = 30;
+
+    DoublePendulum m_pendulum;  // stores properties of the pendulum
+    PendulumSim *p_pendulumSim = NULL;
+
+    // const vector<char const *> m_integrators = {"Analytic", "Explicit Euler", "Symplectic Euler", "Midpoint"};
+    // int m_selected_integrator = 0;
+
+    PendulumGui() {
+        // initialize the pendulum to be used
+        m_pendulum.l1 = 5.0;
+        m_pendulum.l2 = 6.0;
+        m_pendulum.m1 = 0.1;
+        m_pendulum.m2 = 0.3;
+
+        p_pendulumSim = new PendulumSim();
+        p_pendulumSim->setPendulum(m_pendulum);
+        setSimulation(p_pendulumSim);
+
+        // show vertex velocity instead of normal
+        callback_clicked_vertex =
+            [&](int clickedVertexIndex, int clickedObjectIndex, Eigen::Vector3d &pos, Eigen::Vector3d &dir) {
+                RigidObject &o = p_pendulumSim->getObjects()[clickedObjectIndex];
+                pos = o.getVertexPosition(clickedVertexIndex);
+                dir = o.getVelocity(pos);
+            };
+        start();
+    }
+
+    virtual void updateSimulationParameters() override {
+        // change all parameters of the simulation to the values that are set in the GUI
+        p_pendulumSim->setPendulum(m_pendulum);
+        p_pendulumSim->setTimestep(m_dt);
+        p_pendulumSim->setLogFrequency(m_log_frequency);
+    }
+
+    virtual void clearSimulation() override {
+        p_pendulumSim->clearTrajectories();
+    }
+
+    /*
+     * Writes each trajectory to an individual off-file.
+     */
+    void exportTrajectories() {
+        Eigen::MatrixX3d mat;
+        for (int i = 0; i < p_pendulumSim->getNumTrajectories(); i++) {
+            string filename = "trajectory" + to_string(i) + ".off";
+            p_pendulumSim->getTrajectories(i, mat);
+            if (mat.rows() <= 1) {
+                continue;
+            }
+            if (igl::writeOFF(filename, mat, Eigen::MatrixXi())) {
+                cout << "Wrote trajectory to " << filename << endl;
+            } else {
+                cout << "Failed to write trajectory to " << filename << endl;
+            }
+        }
+    }
+
+    virtual bool childKeyCallback(igl::opengl::glfw::Viewer &viewer,
+                                  unsigned int key, int modifiers) override {
+        switch (key) {
+            case 'e':
+            case 'E':
+                exportTrajectories();
+                return true;
+        }
+        return false;
+    }
+
+    virtual void drawSimulationParameterMenu() override {
+        if (ImGui::Button("Export Trajectories", ImVec2(-1, 0))) {
+            exportTrajectories();
+        }
+        ImGui::InputDouble("mass 1", &m_pendulum.m1, 0, 0);
+        ImGui::InputDouble("mass 2", &m_pendulum.m2, 0, 0);
+        ImGui::InputDouble("L1", &m_pendulum.l1, 0, 0);
+        ImGui::InputDouble("L2", &m_pendulum.l2, 0, 0);
+        ImGui::InputFloat("dt", &m_dt, 0, 0);
+        // ImGui::Combo("Integrator", &m_selected_integrator, m_integrators.data(), m_integrators.size());
+        ImGui::InputInt("Log Frequency", &m_log_frequency, 0, 0);
+    }
+};
+
+int main(int argc, char *argv[]) {
+    // create a new instance of the GUI for the pendulum simulation
+    new PendulumGui();
+
+    return 0;
+}
\ No newline at end of file

3_cloth/CMakeLists.txt

+
+cmake_minimum_required(VERSION 3.1)
+project(3_cloth)
+
+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

3_cloth/ClothSim.cpp

+#include "ClothSim.h"
+#include "igl/Timer.h"
+
+using Eigen::Vector2i;
+typedef Eigen::Triplet<double> Triplet;
+typedef Eigen::SparseMatrix<double> SparseMatrix;
+typedef Eigen::Matrix<double, 3, 3> Matrix33;
+
+// add matrix into triplet array with index [i*3, i*3+1, i*3+2]*[j*3, j*3+1, j*3+2]
+static void fillInTriplet(int i, int j, std::vector<Triplet>& triplets, const Matrix33& matrix) {
+    for(int a = 0; a < 3; ++a)
+        for(int b = 0; b < 3; ++b)
+            triplets.push_back(Triplet(i*3+a, j*3+b, matrix(a, b)));
+}
+
+// add v into vector b in index [i*3, i*3+1, i*3+2]
+static void addInVector(int i, Eigen::VectorXd& b, const Eigen::Vector3d& v) {
+    for(int a = 0; a < 3; ++a)
+        b[i*3+a] += v[a];
+}
+
+// To test if a sparse matrix A is symmetric
+static bool isSymmetrical(Eigen::SparseMatrix<double>& A) {
+	bool ret = true;
+	Eigen::SparseMatrix<double> res = Eigen::SparseMatrix<double>(A.transpose()) - A;
+	for (unsigned int k = 0; k < res.outerSize(); k++) {
+		for (Eigen::SparseMatrix<double>::InnerIterator it(res, k); it; ++it)
+		{
+			if (std::fabs(it.value()) > 1e-6) {
+				std::cout<< "row:{} col:{} value: {}" << it.row() << it.col() << A.coeffRef(it.row(), it.col());
+				return false;
+			}
+		}
+	}
+	return ret;
+}
+
+// useful functions:
+// 1. particleCoordinateToIdx()
+// 2. isValidCoor()
+bool ClothSim::advance() {
+    igl::Timer t;
+    t.start();
+    if(m_method == 0) { // euler
+
+        std::vector<Eigen::Vector3d> f;
+        f.resize(particles.size());
+
+        // TODO compute force
+        for(int idx = 0; idx < particles.size(); ++idx) {
+            
+        }
+
+        // TODO constrain fixed particles
+        for(int i = 0; i < fixedParticleIdx.size(); ++i) {
+
+        }
+
+        // TODO update velocity and position of particles
+        for(int idx = 0; idx < particles.size(); ++idx) {
+        }
+
+    } else if (m_method == 1) { // implicit euler
+        
+        std::vector<Triplet> triplets;
+        Eigen::VectorXd b(particles.size()*3); b.setZero();
+        SparseMatrix A(particles.size()*3, particles.size()*3);
+
+        // TODO compute triplets arrays and right-hand-size b
+        // NOTE: remember to process fixed particles, use isFixedParticle to test if a particle is fixed
+        for(int idx = 0; idx < particles.size(); ++idx) {
+
+        }
+        
+        A.setFromTriplets(triplets.begin(), triplets.end());
+        // NOTE: this is for debug purpose, make sure your A matrix is symmetric
+        // if(!isSymmetrical(A)) {std::cout<<"wrong"<<std::endl; return true;}
+        cgSolver.compute(A);
+        Eigen::VectorXd x = cgSolver.solve(b);
+        // std::cout<< "error:" << cgSolver.error() << ", iter:" << cgSolver.iterations() << std::endl;
+
+        // TODO with x computed, update position and velocity of particles
+        for(int i = 0; i < particles.size(); ++i) {
+
+        }
+    } else {    // combination of a half step of explicit Euler and a half step of implicit Euler
+        // TODO
+    }
+    
+    // advance m_time
+    m_time += m_dt;
+    m_step++;
+
+    // log
+    if ((m_step % m_log_frequency) == 0) { 
+        std::cout<<t.getElapsedTimeInMicroSec()<<std::endl;
+    }
+    
+    return false;
+}
\ No newline at end of file

3_cloth/ClothSim.h

+#include "Simulation.h"
+#include "igl/list_to_matrix.h"
+
+using namespace std;
+
+struct Spring {
+public:
+    double damping;
+    double k_struct, k_bend, k_shear;
+    double L_struct, L_bend, L_shear;
+};
+
+struct Particle {
+public:
+    Particle(const Eigen::Vector3d& p = Eigen::Vector3d::Zero(), const Eigen::Vector3d& v=Eigen::Vector3d::Zero())
+        : position(p), velocity(v) {}
+    Eigen::Vector3d position;
+    Eigen::Vector3d velocity;
+};
+
+/*
+ */
+class ClothSim : public Simulation {
+public:
+    ClothSim() : Simulation(), n(10), m(10), domainX(5) {
+        init();
+        // m_trajectories.clear();
+        // m_trajectoryColors.clear();
+    }
+
+    virtual void init() override {
+        std::string path = "sphere.off";
+        m_objects.clear();
+        m_objects.push_back(RigidObject(path));
+        m_objects.push_back(RigidObject(path));
+        p1 = &m_objects[0];
+        p2 = &m_objects[1];
+
+        m_dt = 1e-2;
+        m_gravity << 0, -9.81, 0;
+        n = 10, m = 10;
+        domainX = 5;
+        dx = domainX / (n-1);
+
+        cgSolver.setTolerance(1e-7);
+	    cgSolver.setMaxIterations(1000);
+        reset();
+    }
+
+    virtual void resetMembers() override {
+        double domainY = dx * (m-1);
+        // initialize fixed point visualizations
+        p1->reset(); p2->reset();
+        p1->setScale(0.01); p2->setScale(0.01);
+        p1->setPosition(Eigen::Vector3d(-domainX/2.0, domainY, 0));
+        p2->setPosition(Eigen::Vector3d(domainX/2.0, domainY, 0));
+
+        // initialize particles
+        particles.resize(n * m);
+        for(int i = 0; i < n; ++i) {
+            for(int j = 0; j < m; ++j) {
+                particles[particleCoordinateToIdx(i, j)] = 
+                    Particle(Eigen::Vector3d(dx * i - domainX/2, domainY, dx * j), Eigen::Vector3d::Zero());
+            }
+        }
+
+        // initialize fixed particle
+        fixedParticleIdx.resize(2);
+        fixedParticleIdx[0] = particleCoordinateToIdx(0, 0);
+        fixedParticleIdx[1] = particleCoordinateToIdx(n-1, 0);
+        isFixedParticle.resize(n*m, false);
+        isFixedParticle[particleCoordinateToIdx(0, 0)] = true;
+        isFixedParticle[particleCoordinateToIdx(n-1, 0)] = true;
+
+        // initialize cloth mesh triangles
+        std::vector<std::vector<int>> faces;
+        faces.reserve((n-1)*(m-1)*2);
+        for(int i = 0; i < n-1; ++i) {
+            for(int j = 0; j < m-1; ++j) {
+                faces.push_back({particleCoordinateToIdx(i, j), particleCoordinateToIdx(i, j+1), particleCoordinateToIdx(i+1, j+1)});
+                faces.push_back({particleCoordinateToIdx(i, j), particleCoordinateToIdx(i+1, j+1), particleCoordinateToIdx(i+1, j)});
+            }
+        }
+        igl::list_to_matrix(faces, m_cloth_renderF);
+    }
+
+    virtual void updateRenderGeometry() override {
+        // create cloth geometry
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {    // negative id means newly created object, reverse memory for it
+                m_renderVs.emplace_back();
+                m_renderFs.emplace_back();
+            }
+
+            m_objects[i].getMesh(m_renderVs[i], m_renderFs[i]);
+        }
+
+        m_cloth_renderV.resize(n*m, 3);
+        for(int i = 0; i < n*m; ++i) {
+            m_cloth_renderV.row(i) << particles[i].position.transpose();
+        }
+    }
+
+	virtual bool advance() override;
+
+    virtual void renderRenderGeometry(igl::opengl::glfw::Viewer &viewer) override {
+        // render point
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {
+                int new_id = 0;
+                if (i > 0) {
+                    new_id = viewer.append_mesh();
+                    o.setID(new_id);
+                } else {
+                    o.setID(new_id);
+                }
+
+                size_t meshIndex = viewer.mesh_index(o.getID());
+                viewer.data_list[meshIndex].set_face_based(true);
+                viewer.data_list[meshIndex].point_size = 2.0f;
+                viewer.data_list[meshIndex].clear();
+            }
+            size_t meshIndex = viewer.mesh_index(o.getID());
+
+            viewer.data_list[meshIndex].set_mesh(m_renderVs[i], m_renderFs[i]);
+            viewer.data_list[meshIndex].compute_normals();
+
+            Eigen::MatrixXd color;
+            o.getColors(color);
+            viewer.data_list[meshIndex].set_colors(color);
+        }
+
+        if(cloth_render_id < 0) {
+            cloth_render_id = viewer.append_mesh();
+            viewer.data_list[cloth_render_id].set_face_based(true);
+            viewer.data_list[cloth_render_id].point_size = 2.0f;
+            viewer.data_list[cloth_render_id].clear();
+        }
+        viewer.data_list[cloth_render_id].set_mesh(m_cloth_renderV, m_cloth_renderF);
+        viewer.data_list[cloth_render_id].compute_normals();
+    }
+
+    void setMethod(int m) { m_method = m; }
+    void setSpring(Spring &s) { m_spring = s; }
+    void setLogFrequency(int f) { m_log_frequency = f; }
+    int n, m;
+    double m_mass;
+    double domainX, dx;
+
+private:
+
+    int particleCoordinateToIdx(int i, int j) { return i*m + j; }
+    int particleCoordinateToIdx(const Eigen::Vector2i& coor) { return particleCoordinateToIdx(coor.x(), coor.y()); }
+    Eigen::Vector2i particleIdxToCoordinate(int idx) { return Eigen::Vector2i(int(idx/m), idx%m); }
+    bool isValidIdx(int idx) { return idx >= 0 && idx < particles.size(); }
+    bool isValidCoor(const Eigen::Vector2i& coor) {
+        return coor.x() >= 0 && coor.x() < n && coor.y() >= 0 && coor.y() < m;
+    }
+
+    int m_method;   // (0: euler, 1: implicit euler)
+    Eigen::Vector3d m_gravity;
+    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper, Eigen::DiagonalPreconditioner<double>> cgSolver;
+
+    Spring m_spring;
+    RigidObject *p1, *p2;
+
+    std::vector<int> fixedParticleIdx;
+    std::vector<bool> isFixedParticle;
+    std::vector<Particle> particles;
+    std::vector<Eigen::MatrixXd> m_renderVs;  // vertex positions for rendering
+    std::vector<Eigen::MatrixXi> m_renderFs;  // face indices for rendering
+    int cloth_render_id = -1;
+    Eigen::MatrixXd m_cloth_renderV;  // vertex positions for rendering
+    Eigen::MatrixXi m_cloth_renderF;  // face indices for rendering
+
+    int m_log_frequency;  // how often should we log the COM in the GUI
+    // vector<vector<Eigen::Vector3d>> m_trajectories;
+    // Eigen::RowVector3d m_color, m_color1, m_color2;
+    // vector<Eigen::RowVector3d> m_trajectoryColors;
+};
\ No newline at end of file

3_cloth/main.cpp

+#include <igl/writeOFF.h>
+#include <thread>
+#include "Gui.h"
+#include "Simulator.h"
+#include "ClothSim.h"
+
+/*
+ */
+class ClothGui : public Gui {
+   public:
+    // simulation parameters
+    float m_dt = 1e-2;
+    int n = 10, m = 10;
+    double mass = 1;
+    double domainX = 5;
+    int m_log_frequency = 30;
+
+    Spring m_spring;  // stores properties of the spring
+    ClothSim *p_clothSim = NULL;
+    const vector<char const *> m_integrators = {"Explicit Euler", "Implicit Euler", "Combination"};
+    int m_selected_integrator = 0;
+
+
+    ClothGui() {
+        m_spring.damping = 10;
+        m_spring.k_struct = 3000;
+        m_spring.k_shear = 1000;
+        m_spring.k_bend = 500;
+
+        // initialize the pendulum to be used
+        p_clothSim = new ClothSim();
+        p_clothSim->setSpring(m_spring);
+        setSimulation(p_clothSim);
+
+        start();
+    }
+
+    virtual void updateSimulationParameters() override {
+        // change all parameters of the simulation to the values that are set in the GUI
+        double dx = domainX / (n-1);
+        m_spring.L_struct = dx;
+        m_spring.L_shear = dx * sqrt(2);
+        m_spring.L_bend = 2 * dx;
+        
+        p_clothSim->setSpring(m_spring);
+        p_clothSim->setTimestep(m_dt);
+        p_clothSim->setLogFrequency(m_log_frequency);
+        p_clothSim->setMethod(m_selected_integrator);
+        p_clothSim->m_mass = mass;
+        p_clothSim->n = n;
+        p_clothSim->m = m;
+        p_clothSim->domainX = domainX;
+        p_clothSim->dx = dx;
+    }
+
+    virtual void clearSimulation() override {
+        // 
+    }
+
+
+    virtual void drawSimulationParameterMenu() override {
+        ImGui::InputInt("n", &n, 0, 0);
+        ImGui::InputInt("m", &m, 0, 0);
+        ImGui::InputDouble("mass", &mass, 0, 0);
+        ImGui::InputDouble("Domain X", &domainX, 0, 0);
+        ImGui::InputDouble("struct stiffness", &m_spring.k_struct, 0, 0);
+        ImGui::InputDouble("shear stiffness", &m_spring.k_shear, 0, 0);
+        ImGui::InputDouble("bend stiffness", &m_spring.k_bend, 0, 0);
+        ImGui::InputDouble("damping", &m_spring.damping, 0, 0);
+        ImGui::InputFloat("dt", &m_dt, 0, 0);
+        ImGui::Combo("Integrator", &m_selected_integrator, m_integrators.data(), m_integrators.size());
+        ImGui::InputInt("Log Frequency", &m_log_frequency, 0, 0);
+    }
+};
+
+int main(int argc, char *argv[]) {
+    // create a new instance of the GUI for the pendulum simulation
+    new ClothGui();
+
+    return 0;
+}
\ No newline at end of file

4-1_spinning/CMakeLists.txt

+
+cmake_minimum_required(VERSION 3.1)
+project(4-1_spinning)
+
+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

4-1_spinning/SpinSim.cpp

+#include "SpinSim.h"
+
+bool SpinSim::advance() {
+	Eigen::Vector3d w = p_body->getAngularVelocity();
+	
+	// TODO
+	// update orientation
+	switch (m_method) {
+	case 0: {
+		// matrix-based angular velocity
+		break;
+	}
+
+	case 1: {
+		// quaternion-based angular velocity
+		break;
+	}
+	default:
+            std::cerr << m_method << " is not a valid rotation method."
+                        << std::endl;
+	}
+
+	// advance time
+	m_time += m_dt;
+	m_step++;
+
+	return false;
+}
\ No newline at end of file

4-1_spinning/SpinSim.h

+#include "Simulation.h"
+
+using namespace std;
+
+/*
+ */
+class SpinSim : public Simulation {
+   public:
+    SpinSim() : Simulation() { init(); }
+
+    virtual void init() override {
+        std::string file = "cube.off";
+
+        m_objects.clear();
+        m_objects.push_back(RigidObject(file));
+        p_body = &m_objects.back();
+
+        m_dt = 1e-3;
+
+        reset();
+    }
+
+    virtual void resetMembers() override {
+        p_body->reset();
+        p_body->setMass(m_mass);
+
+        p_body->applyForce(m_force, Eigen::Vector3d(1, 1, 1));
+        p_body->applyForce(-m_force, Eigen::Vector3d(-1, -1, -1));
+
+        // compute linear and angular momentums
+        p_body->setLinearMomentum(p_body->getLinearMomentum() +
+                                  m_dt * p_body->getForce());
+        p_body->setAngularMomentum(p_body->getAngularMomentum() +
+                                   m_dt * p_body->getTorque());
+        p_body->resetForce();
+        p_body->resetTorque();
+    }
+
+    virtual void updateRenderGeometry() override {
+        p_body->getMesh(m_renderV, m_renderF);
+    }
+
+	virtual bool advance() override;
+
+    virtual void renderRenderGeometry(
+        igl::opengl::glfw::Viewer &viewer) override {
+        viewer.data().set_mesh(m_renderV, m_renderF);
+        viewer.data().compute_normals();
+    }
+
+#pragma region SettersAndGetters
+    void setForce(const Eigen::Vector3d &f) { m_force = f; }
+
+    void setMass(double m) { m_mass = m; }
+
+    void setMethod(int m) { m_method = m; }
+
+    Eigen::Vector3d getRotationalEnergy() {
+        return 0.5 * p_body->getInertia().diagonal().cwiseProduct(
+                         p_body->getAngularVelocity());
+    }
+#pragma endregion SettersAndGetters
+
+   private:
+    double m_mass;
+    Eigen::Vector3d m_force;
+    int m_method = 0;
+
+    RigidObject *p_body;
+
+    Eigen::MatrixXd m_renderV;  // vertex positions for rendering
+    Eigen::MatrixXi m_renderF;  // face indices for rendering
+};
\ No newline at end of file

4-1_spinning/main.cpp

+#include <deque>
+#include "Gui.h"
+#include "Simulator.h"
+#include "SpinSim.h"
+
+class SpinGui : public Gui {
+   public:
+    // simulation parameters
+    Eigen::Vector3f m_force;
+    float m_mass = 1.0;
+    float m_dt = 1e-3;
+    int m_maxHistory = 200;
+    std::vector<float> m_energy_history;
+
+    const vector<char const *> m_methods = {"Matrix", "Quaternion"};
+    int m_selected_method = 0;
+
+    SpinSim *p_spinSim = NULL;
+
+    SpinGui() {
+        m_force << 10000, 0, 0;
+        m_energy_history.clear();
+
+        p_spinSim = new SpinSim();
+        p_spinSim->setForce(m_force.cast<double>());
+        p_spinSim->setMass(m_mass);
+        setSimulation(p_spinSim);
+
+        // show vertex velocity instead of normal
+        callback_clicked_vertex =
+            [&](int clickedVertexIndex, int clickedObjectIndex,
+                Eigen::Vector3d &pos, Eigen::Vector3d &dir) {
+                RigidObject &o = p_spinSim->getObjects()[clickedObjectIndex];
+                pos = o.getVertexPosition(clickedVertexIndex);
+                dir = o.getVelocity(pos);
+            };
+        start();
+    }
+
+    virtual void updateSimulationParameters() override {
+        // change all parameters of the simulation to the values that are set in
+        // the GUI
+        p_spinSim->setForce(m_force.cast<double>());
+        p_spinSim->setMass(m_mass);
+        p_spinSim->setTimestep(m_dt);
+        p_spinSim->setMethod(m_selected_method);
+    }
+
+    virtual void drawSimulationParameterMenu() override {
+        ImGui::Combo("Method", &m_selected_method, m_methods.data(),
+                     m_methods.size());
+        ImGui::InputFloat3("Force", m_force.data(), 0);
+        ImGui::InputFloat("Mass", &m_mass, 0, 0, 6);
+        ImGui::InputFloat("dt", &m_dt, 0, 0);
+    }
+
+    virtual void drawSimulationStats() override {
+        Eigen::Vector3d E = p_spinSim->getRotationalEnergy();
+        m_energy_history.push_back(E.cast<float>().cwiseAbs().sum());
+        if (m_energy_history.size() > m_maxHistory)
+            m_energy_history.erase(m_energy_history.begin(),
+                                   m_energy_history.begin() + 1);
+        ImGui::Text("E_x: %.3f", E(0));
+        ImGui::Text("E_y: %.3f", E(1));
+        ImGui::Text("E_z: %.3f", E(2));
+        ImGui::PlotLines("Total Energy", &m_energy_history[0],
+                         m_energy_history.size(), 0, NULL, 0, 1000,
+                         ImVec2(0, 200));
+    }
+};
+
+int main(int argc, char *argv[]) {
+    // create a new instance of the GUI for the spinning simulation
+    new SpinGui();
+
+    return 0;
+}
\ No newline at end of file

4-2_collision/CMakeLists.txt

+
+cmake_minimum_required(VERSION 3.1)
+project(4-2_collision)
+
+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

4-2_collision/CollisionDetection.cpp

+#include "CollisionDetection.h"
+
+void CollisionDetection::computeBroadPhase(int broadPhaseMethod) {
+    // compute possible collisions
+    m_overlappingBodys.clear();
+    
+	switch (broadPhaseMethod) {
+	case 0: { // none
+		for (size_t i = 0; i < m_objects.size(); i++) {
+			for (size_t j = i + 1; j < m_objects.size(); j++) {
+				m_overlappingBodys.push_back(std::make_pair(i, j));
+			}
+		}
+		break;
+	}
+     
+	case 1: {  // AABB
+        // compute bounding boxes
+        std::vector<AABB> aabbs(m_objects.size());
+        for (size_t i = 0; i < aabbs.size(); i++) {
+            aabbs[i].computeAABB(m_objects[i]);
+        }
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            for (size_t j = i + 1; j < m_objects.size(); j++) {
+                // add pair of objects to possible collision if their
+                // bounding boxes overlap
+                if (aabbs[i].testCollision(aabbs[j])) {
+                    m_overlappingBodys.push_back(std::make_pair(i, j));
+                }
+            }
+        }
+        break;
+    }
+	
+	case 2: {
+		// TODO: implement other broad phase algorithm
+		break;
+	}
+	}
+}
+
+void CollisionDetection::computeNarrowPhase(int narrowPhaseMethod) {
+    switch (narrowPhaseMethod) {
+    case 0: {
+        // exhaustive
+        // iterate through all pairs of possible collisions
+        for (auto overlap : m_overlappingBodys) {
+            std::vector<Contact> temp_contacts[2];
+            // compute intersection of a with b first and intersectino
+            // of b with a and store results in temp_contacts
+            for (int switcher = 0; switcher < 2; switcher++) {
+                RigidObject* a =
+                    &m_objects[(!switcher) ? overlap.first
+                                            : overlap.second];
+                RigidObject* b =
+                    &m_objects[(!switcher) ? overlap.second
+                                            : overlap.first];
+
+                Eigen::MatrixXd Va, Vb;
+                Eigen::MatrixXi Fa, Fb;
+                a->getMesh(Va, Fa);
+                b->getMesh(Vb, Fb);
+
+                // iterate through all faces of first object
+                for (int face = 0; face < Fa.rows(); face++) {
+                    // iterate through all edges of given face
+                    for (size_t j = 0; j < 3; j++) {
+                        int start = Fa(face, j);
+                        int end = Fa(face, (j + 1) % 3);
+
+                        // check if there is a collision
+                        ContactType ct = isColliding(
+                            Va.row(start), Va.row(end), Vb, Fb);
+
+                        // find collision and check for duplicates
+                        switch (ct) {
+                            case ContactType::VERTEXFACE: {
+                                auto ret = m_penetratingVertices.insert(
+                                    Fa(face, j));
+                                // if not already in set
+                                if (ret.second) {
+                                    Contact temp_col =
+                                        findVertexFaceCollision(
+                                            Va.row(Fa(face, j)), Vb,
+                                            Fb);
+                                    temp_col.a = a;
+                                    temp_col.b = b;
+                                    temp_col.type =
+                                        ContactType::VERTEXFACE;
+                                    temp_contacts[switcher].push_back(
+                                        temp_col);
+                                }
+                                break;
+                            }
+                            case ContactType::EDGEEDGE: {
+                                int orderedStart = std::min(start, end);
+                                int orderedEnd = std::max(start, end);
+                                auto ret = m_penetratingEdges.insert(
+                                    std::make_pair(orderedStart,
+                                                    orderedEnd));
+                                // if not already in set
+                                if (ret.second) {
+                                    Contact temp_col =
+                                        findEdgeEdgeCollision(
+                                            Va.row(orderedStart),
+                                            Va.row(orderedEnd), Vb, Fb);
+                                    temp_col.a = a;
+                                    temp_col.b = b;
+                                    temp_col.type =
+                                        ContactType::EDGEEDGE;
+                                    temp_contacts[switcher].push_back(
+                                        temp_col);
+                                }
+                                break;
+                            }
+                            case ContactType::NONE: {
+                                break;
+                            }
+                        }
+                    }
+                }
+            }
+
+            // look for vertexFace
+            bool found = false;
+            for (int i = 0; i < 2; i++) {
+                for (auto cont : temp_contacts[i]) {
+                    if (cont.type == ContactType::VERTEXFACE) {
+                        m_contacts.push_back(cont);
+                        found = true;
+                        break;
+                    }
+                }
+                if (found) {
+                    break;
+                }
+            }
+            if (found) {
+                continue;
+            }
+
+            // take single edgeedge if possible
+            if (temp_contacts[0].size() > 0 &&
+                temp_contacts[0].size() < temp_contacts[1].size()) {
+                for (int i = 0; i < temp_contacts[0].size(); i++) {
+                    m_contacts.push_back(temp_contacts[0][i]);
+                }
+            } else if (temp_contacts[1].size() > 0 &&
+                        temp_contacts[0].size() >
+                            temp_contacts[1].size()) {
+                for (int i = 0; i < temp_contacts[1].size(); i++) {
+                    m_contacts.push_back(temp_contacts[1][i]);
+                }
+            } else if (temp_contacts[0].size() > 0) {
+                for (int i = 0; i < temp_contacts[0].size(); i++) {
+                    m_contacts.push_back(temp_contacts[0][i]);
+                }
+            } else if (temp_contacts[1].size() > 0) {
+                for (int i = 0; i < temp_contacts[1].size(); i++) {
+                    m_contacts.push_back(temp_contacts[1][i]);
+                }
+            }
+        }
+        break;
+    }
+
+	case 1: {
+		// TODO: implement other narrow phase algorithm
+		break;
+	}
+    }
+}
+
+void CollisionDetection::applyImpulse(double eps) {
+    // compute impulse for all contacts
+    for (auto contact : m_contacts) {
+        Eigen::Vector3d vrel_vec = contact.a->getVelocity(contact.p) -
+                                    contact.b->getVelocity(contact.p);
+        double vrel = contact.n.dot(vrel_vec);
+        if (vrel > 0) {
+            // bodies are moving apart
+            continue;
+        }
+
+		// TODO: compute impulse and update the following momentums
+		//contact.a->setLinearMomentum();
+		//contact.b->setLinearMomentum();
+		//contact.a->setAngularMomentum();
+		//contact.b->setAngularMomentum();
+    }
+}
\ No newline at end of file

4-2_collision/CollisionDetection.h

+#ifndef COLLISIONDETECTION_H
+#define COLLISIONDETECTION_H
+
+#include <igl/ray_mesh_intersect.h>
+#include <set>
+#include <utility>
+#include <vector>
+#include "AABB.h"
+#include "RigidObject.h"
+
+enum class ContactType { VERTEXFACE, EDGEEDGE, NONE };
+
+struct Contact {
+    RigidObject* a;      // body containing vertex
+    RigidObject* b;      // body containing face
+    Eigen::Vector3d n;   // outwards pointing normal of face
+    Eigen::Vector3d p;   // world-space vertex location
+    Eigen::Vector3d ea;  // edge direction for A
+    Eigen::Vector3d eb;  // edge direction for B
+
+    ContactType type;  // type of contact
+};
+
+class CollisionDetection {
+   public:
+    CollisionDetection(std::vector<RigidObject>& world) : m_objects(world) {}
+
+    // pass objects in scene to collision detection
+    void setObjects(std::vector<RigidObject>& world) { m_objects = world; }
+
+	void computeBroadPhase(int broadPhaseMethod);
+
+    // test if ray from start towards end intersects object with vertices V and
+    // faces F
+    ContactType isColliding(const Eigen::Vector3d& start,
+                            const Eigen::Vector3d& end,
+                            const Eigen::MatrixXd& V,
+                            const Eigen::MatrixXi& F) {
+        std::vector<igl::Hit> hits;
+        igl::ray_mesh_intersect(start, end - start, V, F, hits);
+
+        // count number of intersections with object that happen between start
+        // and end (so along the edge)
+        int cntr = 0;
+        for (auto hit : hits) {
+            if (hit.t > 0 && hit.t <= 1) {
+                cntr++;
+            }
+        }
+
+        // if hits is odd then the ray enters through one face and
+        // does not leave again through another, hence the starting point is
+        // inside the object, if the number of intersections between start and
+        // end of edge are even, then the edge entering through one and leaving
+        // through another face, hence the edge is intersection the object
+        ContactType ret = ContactType::NONE;
+        if (hits.size() % 2 == 1) {
+            ret = ContactType::VERTEXFACE;
+        }
+        if (cntr % 2 == 0 && cntr > 0) {
+            ret = ContactType::EDGEEDGE;
+        }
+        return ret;
+    }
+
+    // compute normal for all faces and use it to find closesest face to
+    // given vertex
+    Contact findVertexFaceCollision(const Eigen::Vector3d& vertex,
+                                    const Eigen::MatrixXd& V,
+                                    const Eigen::MatrixXi& F) {
+        double minDist = std::numeric_limits<double>::infinity();
+        Eigen::Vector3d minNormal;
+        for (int i = 0; i < F.rows(); i++) {
+            Eigen::Vector3d a = V.row(F(i, 0));
+            Eigen::Vector3d b = V.row(F(i, 1));
+            Eigen::Vector3d c = V.row(F(i, 2));
+
+            Eigen::Vector3d n = -(b - a).cross(c - a).normalized();
+
+            Eigen::Vector3d v = vertex - a;
+
+            double distance = v.dot(n);
+            // if vertex inside
+            if (distance >= 0) {
+                if (distance < minDist) {
+                    minDist = distance;
+                    minNormal = -n;
+                }
+            }
+        }
+        Contact ret;
+        ret.n = minNormal;
+        ret.p = vertex + minDist * minNormal;
+        return ret;
+    }
+
+    // loop over edges of faces and compute closest to given edge
+    Contact findEdgeEdgeCollision(const Eigen::Vector3d& start,
+                                  const Eigen::Vector3d& end,
+                                  const Eigen::MatrixXd& V,
+                                  const Eigen::MatrixXi& F) {
+        double minDist = std::numeric_limits<double>::infinity();
+        Contact ret;
+        for (int i = 0; i < F.rows(); i++) {
+            Eigen::Vector3d a = V.row(F(i, 0));
+            Eigen::Vector3d b = V.row(F(i, 1));
+            Eigen::Vector3d c = V.row(F(i, 2));
+
+            Eigen::Vector3d n_face = -(b - a).cross(c - a).normalized();
+
+            for (int j = 0; j < 3; j++) {
+                Eigen::Vector3d s = V.row(F(i, j));
+                Eigen::Vector3d e = V.row(F(i, (j + 1) % 3));
+
+                Eigen::Vector3d ea = end - start;
+                Eigen::Vector3d eb = e - s;
+                Eigen::Vector3d n =
+                    (ea).cross(eb);  // direction of shortest distance
+                double distance = n.dot(start - s) / n.norm();
+
+                Eigen::Vector3d plane_normal = n.cross(eb).normalized();
+                double t =
+                    (s - start).dot(plane_normal) / (ea.dot(plane_normal));
+                if (n_face.dot(n) < 0 && distance < 0 && -distance < minDist &&
+                    t >= 0 && t <= 1) {
+                    ret.ea = ea;
+                    ret.eb = eb;
+                    ret.n = n;
+                    ret.p = start + t * ea;
+                    minDist = -distance;
+                }
+            }
+        }
+        return ret;
+    }
+
+	void computeNarrowPhase(int narrowPhaseMethod);
+
+	void applyImpulse(double eps = 1.0);
+
+    void clearDataStructures() {
+        m_penetratingEdges.clear();
+        m_penetratingVertices.clear();
+        m_overlappingBodys.clear();
+        m_contacts.clear();
+    }
+
+    void computeCollisionDetection(int broadPhaseMethod = 0,
+                                   int narrowPhaseMethod = 0,
+                                   double eps = 1.0) {
+        clearDataStructures();
+
+        computeBroadPhase(broadPhaseMethod);
+
+        computeNarrowPhase(narrowPhaseMethod);
+
+        applyImpulse(eps);
+    }
+
+    std::vector<Contact> getContacts() { return m_contacts; }
+
+    std::vector<RigidObject>& m_objects;  // all objects in scene
+    // result of broadphase, pairs of objects with possible collisions
+    std::vector<std::pair<size_t, size_t>> m_overlappingBodys;
+
+    // set of vertex indices that penetrate a face, used to avoid duplicates
+    std::set<int> m_penetratingVertices;
+    // set of pairs of vertex indices that represent a penetrating edge, used to
+    // avoid duplicates
+    std::set<std::pair<int, int>> m_penetratingEdges;
+
+    // computed contact points
+    std::vector<Contact> m_contacts;
+};
+
+#endif

4-2_collision/CollisionSim.h

+#include "CollisionDetection.h"
+#include "Simulation.h"
+
+#include <deque>
+
+using namespace std;
+
+/*
+ */
+class CollisionSim : public Simulation {
+   public:
+    CollisionSim() : Simulation(), m_collisionDetection(m_objects) { init(); }
+
+    const int NUM_CUBES = 5;  // complexity of scene
+
+    virtual void init() override {
+        std::string path = "cube.off";
+        for (int i = 0; i < NUM_CUBES; i++) {
+            m_objects.push_back(RigidObject(path));
+        }
+		path = "cube_shallow.off";
+		for (int i = 0; i < 5; ++i) {
+			m_objects.push_back(RigidObject(path));
+		}
+		
+        m_collisionDetection.setObjects(m_objects);
+
+        m_dt = 1e-3 * 3;
+        m_gravity << 0, -9.81, 0;
+        m_mass = 1.0;
+        m_showContacts = false;
+        m_broadPhaseMethod = 0;
+        m_narrowPhaseMethod = 0;
+        m_eps = 1.0;
+
+        reset();
+    }
+
+    virtual void resetMembers() override {
+        for (auto &o : m_objects) {
+            o.reset();
+        }
+
+		double x1 = 5;
+		double y1 = 4;
+		m_objects[0].setPosition(Eigen::Vector3d(0, y1, 0));
+		m_objects[0].setRotation(
+			Eigen::Quaterniond(0, -0.3444844, -0.3444844, -0.8733046)); 
+		
+		Eigen::RowVector3d c0(204.0 / 255.0, 0, 0);
+		Eigen::RowVector3d c1(0, 128.0 / 255.0, 204.0 / 255.0);
+		for (int i = 1; i < NUM_CUBES; ++i) {
+			m_objects[i].setPosition(Eigen::Vector3d(x1, y1*(i+1), 0));
+			double a = double(i+1) / NUM_CUBES;
+			m_objects[i].setColors(c0*a + c1*(1 - a));
+		}
+
+        for (size_t i = 0; i < NUM_CUBES; i++) {
+            m_objects[i].setMass(m_mass);
+        }
+
+		for (int i = 0; i < 5; ++i) {
+			m_objects[i + NUM_CUBES].setScale(10);
+			m_objects[i + NUM_CUBES].setType(ObjType::STATIC);
+			m_objects[i + NUM_CUBES].setColors(Eigen::RowVector3d(0.5, 0.5, 0.5));
+			m_objects[i + NUM_CUBES].setMass(std::numeric_limits<double>::max());
+		}
+
+		m_objects[0 + NUM_CUBES].setPosition(
+			Eigen::Vector3d(0, -0.1, 0));
+
+		m_objects[1 + NUM_CUBES].setPosition(
+			Eigen::Vector3d(-10.1, 10, 0));
+		m_objects[1 + NUM_CUBES].setRotation(
+			Eigen::Quaterniond(0.7071, 0, 0, 0.7071));
+
+		m_objects[2 + NUM_CUBES].setPosition(
+			Eigen::Vector3d(10.1, 10, 0));
+		m_objects[2 + NUM_CUBES].setRotation(
+			Eigen::Quaterniond(0.7071, 0, 0, 0.7071));
+
+		m_objects[3 + NUM_CUBES].setPosition(
+			Eigen::Vector3d(0, 10, -10.1));
+		m_objects[3 + NUM_CUBES].setRotation(
+			Eigen::Quaterniond(0.5, 0.5, 0.5, 0.5));
+
+		m_objects[4 + NUM_CUBES].setPosition(
+			Eigen::Vector3d(0, 10, 10.1));
+		m_objects[4 + NUM_CUBES].setRotation(
+			Eigen::Quaterniond(0.5, 0.5, 0.5, 0.5));
+
+
+        updateVars();
+    }
+
+    virtual void updateRenderGeometry() override {
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {
+                m_renderVs.emplace_back();
+                m_renderFs.emplace_back();
+            }
+
+            m_objects[i].getMesh(m_renderVs[i], m_renderFs[i]);
+        }
+    }
+
+    virtual bool advance() override {
+        // compute the collision detection
+        m_collisionDetection.computeCollisionDetection(m_broadPhaseMethod, m_narrowPhaseMethod, m_eps);
+
+        // apply forces (only gravity in this case)
+        for (auto &o : m_objects) {
+            o.applyForceToCOM(m_gravity);
+        }
+
+        for (auto &o : m_objects) {
+            // integrate velocities
+            o.setLinearMomentum(o.getLinearMomentum() + m_dt * o.getForce());
+            o.setAngularMomentum(o.getAngularMomentum() + m_dt * o.getTorque());
+            o.resetForce();
+            o.resetTorque();
+
+            // angular velocity (matrix)
+            Eigen::Vector3d w = o.getAngularVelocity();
+            Eigen::Quaterniond wq;
+            wq.w() = 0;
+            wq.vec() = w;
+
+            // integrate position and rotation
+            o.setPosition(o.getPosition() + m_dt * o.getLinearVelocity());
+            Eigen::Quaterniond q = o.getRotation();
+            Eigen::Quaterniond dq = wq * q;
+            Eigen::Quaterniond new_q;
+            new_q.w() = q.w() + 0.5 * m_dt * dq.w();
+            new_q.vec() = q.vec() + 0.5 * m_dt * dq.vec();
+            o.setRotation(new_q.normalized());
+        }
+        // advance time
+        m_time += m_dt;
+        m_step++;
+
+        return false;
+    }
+
+    virtual void renderRenderGeometry(
+        igl::opengl::glfw::Viewer &viewer) override {
+        for (size_t i = 0; i < m_objects.size(); i++) {
+            RigidObject &o = m_objects[i];
+            if (o.getID() < 0) {
+                int new_id = 0;
+                if (i > 0) {
+                    new_id = viewer.append_mesh();
+                    o.setID(new_id);
+                } else {
+                    o.setID(new_id);
+                }
+
+                size_t meshIndex = viewer.mesh_index(o.getID());
+				if (i >= NUM_CUBES) {
+					viewer.data_list[meshIndex].show_lines = true;
+					viewer.data_list[meshIndex].show_faces = false;
+				}
+				else {
+					viewer.data_list[meshIndex].show_lines = false;
+				}
+                viewer.data_list[meshIndex].set_face_based(true);
+                viewer.data_list[meshIndex].point_size = 2.0f;
+                viewer.data_list[meshIndex].clear();
+            }
+            size_t meshIndex = viewer.mesh_index(o.getID());
+
+            viewer.data_list[meshIndex].set_mesh(m_renderVs[i], m_renderFs[i]);
+            viewer.data_list[meshIndex].compute_normals();
+
+            Eigen::MatrixXd color;
+            o.getColors(color);
+            viewer.data_list[meshIndex].set_colors(color);
+        }
+
+        if (m_showContacts) {
+            // number of timesteps to keep showing collision
+            int delay = 10;
+
+            // clear old points
+            viewer.data_list[1].points = Eigen::MatrixXd(0, 6);
+            viewer.data_list[1].point_size = 10.0f;
+
+            // remove expired points
+            while (m_contactMemory.size() > 0 &&
+                   m_contactMemory.front().second + delay < m_step) {
+                m_contactMemory.pop_front();
+            }
+
+            // get new points and add them to memory
+            auto contacts = m_collisionDetection.getContacts();
+            for (auto &contact : contacts) {
+                m_contactMemory.push_back(std::make_pair(contact, m_step));
+            }
+
+            // show points
+            for (auto &contact_int_p : m_contactMemory) {
+                viewer.data_list[1].add_points(
+                    contact_int_p.first.p.transpose(),
+                    (contact_int_p.first.type == ContactType::EDGEEDGE)
+                        ? Eigen::RowVector3d(0, 1, 0)
+                        : Eigen::RowVector3d(0, 0, 1));
+            }
+        }
+    }
+
+#pragma region SettersAndGetters
+    /*
+     * Compute magnitude and direction of momentum and apply it to o
+     */
+    void updateVars() {
+        Eigen::Vector3d momentum;
+        momentum << std::sin(m_angle), std::cos(m_angle), 0;
+        momentum *= m_force;
+        m_objects[0].setLinearMomentum(momentum);
+    }
+
+    void setAngle(double a) {
+        m_angle = a;
+        updateVars();
+    }
+
+    void setForce(double f) {
+        m_force = f;
+        updateVars();
+    }
+
+    void setMass(double m) { m_mass = m; }
+
+    void showContacts(bool s) {
+        if (!s) {
+            m_contactMemory.clear();
+        }
+        m_showContacts = s;
+    }
+
+    void setBroadPhaseMethod(int m) { m_broadPhaseMethod = m; }
+    void setNarrowPhaseMethod(int m) { m_narrowPhaseMethod = m; }
+
+    void setEps(double eps) { m_eps = eps; }
+
+    Eigen::Vector3d getKineticEnergy() {
+        Eigen::Vector3d res;
+        res.setZero();
+        for (auto o : m_objects) {
+            if (o.getType() == ObjType::STATIC) continue;
+            Eigen::Vector3d rotE = 0.5 * o.getInertia().diagonal().cwiseProduct(
+                                             o.getAngularVelocity());
+            Eigen::Vector3d kinE =
+                0.5 * o.getMass() * o.getLinearVelocity().array().square();
+            res += rotE + kinE;
+        }
+        return res;
+    }
+
+    Eigen::Vector3d getLinearMomentum() {
+        Eigen::Vector3d res;
+        res.setZero();
+        for (auto o : m_objects) {
+            if (o.getType() == ObjType::STATIC) continue;
+            res += o.getLinearMomentum();
+        }
+        return res;
+    }
+
+    Eigen::Vector3d getAngularMomentum() {
+        Eigen::Vector3d res;
+        res.setZero();
+        for (auto o : m_objects) {
+            if (o.getType() == ObjType::STATIC) continue;
+            res += o.getAngularMomentum();
+        }
+        return res;
+    }
+
+#pragma endregion SettersAndGetters
+
+   private:
+    [[maybe_unused]] int m_method;  // id of integrator to be used (0: analytical, 1: explicit
+                   // euler, 2: semi-implicit euler)
+    double m_angle;
+    double m_force;
+    double m_mass;
+
+    Eigen::Vector3d m_gravity;
+
+    CollisionDetection m_collisionDetection;
+    int m_broadPhaseMethod;
+    int m_narrowPhaseMethod;
+    double m_eps;
+
+    std::vector<Eigen::MatrixXd> m_renderVs;  // vertex positions for rendering
+    std::vector<Eigen::MatrixXi> m_renderFs;  // face indices for rendering
+
+    bool m_showContacts;
+    std::deque<std::pair<Contact, int>> m_contactMemory;
+};
\ No newline at end of file

4-2_collision/main.cpp

+#include <igl/writeOFF.h>
+#include "CollisionSim.h"
+#include "Gui.h"
+
+class CollisionGui : public Gui {
+   public:
+    float m_angle = 1.047f;
+    float m_force = 10.0f;
+    float m_dt = 1e-1 * 3;
+    float m_mass = 1.0;
+    bool m_showContacts = false;
+    bool m_useAABB = false;
+    float m_eps = 0.6;
+
+    int m_maxHistory = 200;
+    std::vector<float> m_energy_history;
+
+    int m_selectedBroadPhase = 0;
+    const std::vector<char const *> m_broadphases = {"None", "AABB", "Own"};
+    int m_selectedNarrowPhase = 0;
+    const std::vector<char const *> m_narrowphases = {"Exhaustive", "Own"};
+
+    CollisionSim *p_CollisionSim = NULL;
+
+    CollisionGui() {
+        p_CollisionSim = new CollisionSim();
+        setSimulation(p_CollisionSim);
+
+        // show vertex velocity instead of normal
+        callback_clicked_vertex = [&](int clickedVertexIndex,
+                                      int clickedObjectIndex,
+                                      Eigen::Vector3d &pos,
+                                      Eigen::Vector3d &dir) {
+            RigidObject &o = p_CollisionSim->getObjects()[clickedObjectIndex];
+            pos = o.getVertexPosition(clickedVertexIndex);
+            dir = o.getVelocity(pos);
+        };
+        start();
+    }
+
+    virtual void updateSimulationParameters() override {
+        p_CollisionSim->setForce(m_force);
+        p_CollisionSim->setAngle(m_angle);
+        p_CollisionSim->setTimestep(m_dt);
+        p_CollisionSim->setMass(m_mass);
+        p_CollisionSim->setEps(m_eps);
+    }
+
+    virtual void clearSimulation() override {
+        p_CollisionSim->showContacts(false);
+        p_CollisionSim->showContacts(m_showContacts);
+    }
+
+    virtual void drawSimulationParameterMenu() override {
+        ImGui::SliderAngle("Angle", &m_angle, -180.0f, 180.0f);
+        ImGui::InputFloat("Force", &m_force, 0, 0);
+        ImGui::InputFloat("Mass", &m_mass, 0, 0);
+        ImGui::InputFloat("dt", &m_dt, 0, 0);
+        if (ImGui::Checkbox("Show contacts", &m_showContacts)) {
+            p_CollisionSim->showContacts(m_showContacts);
+        }
+        if (ImGui::Combo("Broadphase", &m_selectedBroadPhase,
+                         m_broadphases.data(), m_broadphases.size())) {
+            p_CollisionSim->setBroadPhaseMethod(m_selectedBroadPhase);
+        }
+        if (ImGui::Combo("Narrowphase", &m_selectedNarrowPhase,
+                         m_narrowphases.data(), m_narrowphases.size())) {
+            p_CollisionSim->setNarrowPhaseMethod(m_selectedNarrowPhase);
+        }
+        ImGui::InputFloat("eps", &m_eps, 0, 0);
+    }
+
+    virtual void drawSimulationStats() override {
+        Eigen::Vector3d E = p_CollisionSim->getKineticEnergy();
+        m_energy_history.push_back(E.cast<float>().cwiseAbs().sum());
+        if (m_energy_history.size() > m_maxHistory)
+            m_energy_history.erase(m_energy_history.begin(),
+                                   m_energy_history.begin() + 1);
+        ImGui::Text("E: %.3f, %.3f, %.3f", E(0), E(1), E(2));
+        ImGui::PlotLines("Total Energy", &m_energy_history[0],
+                         m_energy_history.size(), 0, NULL, 0, 1000,
+                         ImVec2(0, 200));
+        Eigen::Vector3d p = p_CollisionSim->getLinearMomentum();
+        ImGui::Text("M: %.3f, %.3f, %.3f", p(0), p(1), p(2));
+        Eigen::Vector3d l = p_CollisionSim->getAngularMomentum();
+        ImGui::Text("L: %.3f, %.3f, %.3f", l(0), l(1), l(2));
+    }
+};
+
+int main(int argc, char *argv[]) {
+    new CollisionGui();
+
+    return 0;
+}
\ No newline at end of file

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