From e1f0c16a9b888a1a2f78806e7c9d96abf581c143 Mon Sep 17 00:00:00 2001
From: Yuwei Xiao <xyw1105@126.com>
Date: Tue, 24 Mar 2020 21:59:45 +0800
Subject: [PATCH] ex3
---
3_cloth/CMakeLists.txt | 43 ++++++++++
3_cloth/ClothSim.cpp | 100 ++++++++++++++++++++++
3_cloth/ClothSim.h | 182 +++++++++++++++++++++++++++++++++++++++++
3_cloth/main.cpp | 81 ++++++++++++++++++
CMakeLists.txt | 1 +
5 files changed, 407 insertions(+)
create mode 100644 3_cloth/CMakeLists.txt
create mode 100644 3_cloth/ClothSim.cpp
create mode 100644 3_cloth/ClothSim.h
create mode 100644 3_cloth/main.cpp
diff --git a/3_cloth/CMakeLists.txt b/3_cloth/CMakeLists.txt
new file mode 100644
index 0000000..d1a557c
--- /dev/null
+++ b/3_cloth/CMakeLists.txt
@@ -0,0 +1,43 @@
+
+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
diff --git a/3_cloth/ClothSim.cpp b/3_cloth/ClothSim.cpp
new file mode 100644
index 0000000..28bcf35
--- /dev/null
+++ b/3_cloth/ClothSim.cpp
@@ -0,0 +1,100 @@
+#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
diff --git a/3_cloth/ClothSim.h b/3_cloth/ClothSim.h
new file mode 100644
index 0000000..fd93b58
--- /dev/null
+++ b/3_cloth/ClothSim.h
@@ -0,0 +1,182 @@
+#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
diff --git a/3_cloth/main.cpp b/3_cloth/main.cpp
new file mode 100644
index 0000000..ee587a5
--- /dev/null
+++ b/3_cloth/main.cpp
@@ -0,0 +1,81 @@
+#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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c49e975..d67a197 100644
--- a/CMakeLists.txt
+++ b/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
--
GitLab