ex2

.gitignore

@@ -5,4 +5,6 @@ build/
 *.aux
 *.log
 *.gz
-*.out
\ No newline at end of file
+*.out
+
+(old)*
\ No newline at end of file

2-1_bead/BeadSim.cpp

+#include "BeadSim.h"
+
+
+bool BeadSim::advance() {
+
+
+    auto v = p_bead->getLinearVelocity();
+    auto p = p_bead->getPosition();
+
+    // TODO update position and velcity of p_bead
+    // constraint C(x) = 0.5*(p.dot(p) - m_radius^2) = 0;
+
+
+    m_time += m_dt;
+    m_step++;
+
+    // log
+    if ((m_step % m_log_frequency) == 0) {
+        m_trajectories.back().push_back(p_bead->getPosition());
+    }
+    return false;
+
+}
\ No newline at end of file

2-1_bead/BeadSim.h

+#include "Simulation.h"
+#include "igl/PI.h"
+
+using namespace std;
+
+
+class BeadSim : public Simulation {
+public:
+    BeadSim() : 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));
+		p_bead = &m_objects[0];
+
+		m_dt = 5e-2;
+        m_mass = 1.0;
+        m_log_frequency = 30;
+		m_radius = 10;
+        m_gravity << 0, -9.81, 0;
+
+		reset();
+	}
+
+    virtual void resetMembers() override {
+        p_bead->reset();
+        p_bead->setScale(0.1);
+        // initial position, should be valid for constraints
+        p_bead->setPosition(Eigen::Vector3d(m_radius*cos(igl::PI/3), m_radius*sin(igl::PI/3.0), 0));
+
+        if (m_trajectories.size() == 0 || m_trajectories.back().size() > 1) {
+            m_trajectories.push_back(vector<Eigen::Vector3d>());
+            m_trajectoryColors.push_back(m_color);
+        } else {
+            m_trajectoryColors.back() = m_color;
+        }
+    }
+
+    virtual void updateRenderGeometry() override {
+        p_bead->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);
+       
+        for (size_t trajectory = 0; trajectory < m_trajectories.size(); trajectory++) {
+            for (size_t point = 1; point < m_trajectories[trajectory].size(); point++) {
+                viewer.data().add_edges(m_trajectories[trajectory][point-1].transpose(),
+                        m_trajectories[trajectory][point].transpose(),
+                        m_trajectoryColors[trajectory]);
+                // viewer.data().add_points(
+                //     m_trajectories[trajectory][point].transpose(),
+                //     m_trajectoryColors[trajectory]);
+            }
+        }
+    }
+
+    void clearTrajectories() {
+        m_trajectories.clear();
+        m_trajectories.push_back(vector<Eigen::Vector3d>());
+        m_trajectoryColors.clear();
+        m_trajectoryColors.push_back(m_color);
+    }
+    void setRadius(double r) { m_radius = r; }
+    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:
+    RigidObject *p_bead;
+    double m_radius;
+    double m_mass;
+
+    Eigen::Vector3d m_gravity;
+
+    Eigen::MatrixXd m_renderV;  // vertex positions for rendering
+    Eigen::MatrixXi m_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;
+    vector<Eigen::RowVector3d> m_trajectoryColors;
+};
\ No newline at end of file

2-1_bead/CMakeLists.txt

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

2-1_bead/main.cpp

+#include "BeadSim.h"
+#include "Gui.h"
+
+/*
+ * GUI for the bead simulation.
+ */
+class BeadGui : public Gui {
+   public:
+    // simulation parameters
+    float m_dt = 1e-2;
+    float m_radius = 10.0;
+    int m_log_frequency = 30;
+
+    BeadSim *p_BeadSim = NULL;
+
+    BeadGui() {
+        p_BeadSim = new BeadSim();
+        setSimulation(p_BeadSim);
+
+        start();
+    }
+
+    virtual void updateSimulationParameters() override {
+        // change all parameters of the simulation to the values that are set in the GUI
+        p_BeadSim->setTimestep(m_dt);
+        p_BeadSim->setRadius(m_radius);
+        p_BeadSim->setLogFrequency(m_log_frequency);
+    }
+
+
+    virtual void drawSimulationParameterMenu() override {
+        ImGui::InputFloat("Radius", &m_radius, 0, 0);
+        ImGui::InputFloat("dt", &m_dt, 0, 0);
+    }
+
+    virtual void clearSimulation() override {
+        p_BeadSim->clearTrajectories();
+    }
+};
+
+int main(int argc, char *argv[]) {
+    new BeadGui();
+
+    return 0;
+}
\ No newline at end of file

2-2_pendulum/CMakeLists.txt

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

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

CMakeLists.txt

@@ -28,6 +28,8 @@ if (NOT LIBIGL_FOUND)
 endif()
 
 # add_subdirectory(0_dummy)
-add_subdirectory(1-0_earth)
-add_subdirectory(1-1_cannonball)
-add_subdirectory(1-2_spring)
\ No newline at end of file
+# add_subdirectory(1-0_earth)
+# add_subdirectory(1-1_cannonball)
+# add_subdirectory(1-2_spring)
+add_subdirectory(2-1_bead)
+add_subdirectory(2-2_pendulum)