ex3

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

CMakeLists.txt

@@ -33,3 +33,4 @@ endif()
 # add_subdirectory(1-2_spring)
 add_subdirectory(2-1_bead)
 add_subdirectory(2-2_pendulum)
+add_subdirectory(3_cloth)
\ No newline at end of file