#include "ParticleSimApplication.h" #include "ParticleSimGLCanvas.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace gti320; namespace { static const float deltaT = 0.01667f; /** * Crée un système masse-ressort qui simule un tissu suspendu */ static inline void createHangingCloth(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 16; const int x_start = 240; const int y_start = 80; const int dx = 32; const int dy = 32; int index = 0; for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { const int x = x_start + j * dx; const int y = y_start + i * dy; Particle particle(Vector2f(x, y), Vector2f(0, 0), Vector2f(0, 0), 1.0); if (j == 0 && i == (N - 1)) particle.fixed = true; if (j == (N - 1) && i == (N - 1)) particle.fixed = true; particleSystem.addParticle(particle); if (i > 0) { Spring s(index - N, index, k, (float) dy); particleSystem.addSpring(s); } if (j > 0) { Spring s(index - 1, index, k, (float) dx); particleSystem.addSpring(s); } if (i > 0 && j > 0) { Spring s(index - N - 1, index, k, std::sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s); } ++index; } } } /** * Crée un système masse-ressort qui simule un grand tissu suspendu */ static inline void createLargeHangingCloth(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 32; const int x_start = 240; const int y_start = 80; const int dx = 16; const int dy = 16; int index = 0; for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { const int x = x_start + j * dx; const int y = y_start + i * dy; Particle particle(Vector2f(x, y), Vector2f(0, 0), Vector2f(0, 0), 1.0); if (j == 0 && i == (N - 1)) particle.fixed = true; if (j == (N - 1) && i == (N - 1)) particle.fixed = true; particleSystem.addParticle(particle); if (i > 0) { Spring s(index - N, index, k, (float) dy); particleSystem.addSpring(s); } if (j > 0) { Spring s(index - 1, index, k, (float) dx); particleSystem.addSpring(s); } if (i > 0 && j > 0) { Spring s(index - N - 1, index, k, std::sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s); } ++index; } } } /** * Crée un système masse-ressort qui simule une corde suspendu par ses * extrémités. */ static inline void createHangingRope(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 20; const int x_start = 200; const int dx = 32; int index = 0; for (int j = 0; j < N; ++j) { const int x = x_start + j * dx; const int y = 480; Particle particle(Vector2f(x, y), Vector2f(0, 0), Vector2f(0, 0), 1.0); particle.fixed = (index == 0) || (index == N - 1); particleSystem.addParticle(particle); if (j > 0) { Spring s(index - 1, index, k, (float) dx); particleSystem.addSpring(s); } ++index; } } /** * Crée un système masse-ressort qui simule une poutre flexible */ static inline void createBeam(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 20; const int x_start = 200; const int y_start = 400; const int dx = 32; const int dy = 32; int index = 0; for (int j = 0; j < N; ++j) { const int x = x_start + j * dx; // Bottom particle { Particle particle(Vector2f(x, y_start), Vector2f(0, 0), Vector2f(0, 0), 1.0); particle.fixed = (j == 0); particleSystem.addParticle(particle); if (j > 0) { Spring s(index - 1, index, k, (float) sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s); Spring s2(index - 2, index, k, (float) dx); particleSystem.addSpring(s2); } ++index; } // Top particle { Particle particle(Vector2f(x, y_start + dy), Vector2f(0, 0), Vector2f(0, 0), 1.0); particle.fixed = (j == 0); particleSystem.addParticle(particle); Spring s(index - 1, index, k, (float) dy); particleSystem.addSpring(s); if (j > 0) { Spring s2(index - 2, index, k, (float) dx); particleSystem.addSpring(s2); Spring s3(index - 3, index, k, (float) sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s3); } ++index; } } } /** * Créez votre propre exemple */ static inline void createVotreExemple(ParticleSystem &particleSystem, float k) { particleSystem.clear(); // Amusez-vous. Rendu ici, vous le méritez. // Pont const int x_start = 200; const int y_start = 200; const int pillar_n = 6; const int pillar_dy = 64; const int bridge_n = 10; const int bridge_dx = 32; const int bridge_dy = 32; // Piliers de gauche et droite for (int pillar = 0; pillar < 2; pillar++) { const auto pillar_x = (float) (x_start + pillar * bridge_dx * (bridge_n + 1)); const int spring_offset = pillar * pillar_n; for (int i = 0; i < pillar_n; i++) { Vector2f position(pillar_x, (float) (y_start + i * pillar_dy)); Particle p(position, Vector2f(0, 0), Vector2f(0, 0), 1.0f); p.fixed = i == 0 || i == pillar_n - 1; particleSystem.addParticle(p); if (i > 0) { // Le ressort a déjà de la tension pour pas que tout s'écroule immédiatement Spring s(spring_offset + i - 1, spring_offset + i, 0, pillar_dy / 1.2f); particleSystem.addSpring(s); } } } const int pillar_connection_index = pillar_n / 2; const int bridge_y = y_start + pillar_connection_index * pillar_dy; const auto l0_diag = (float) std::sqrt(std::pow(bridge_dx, 2) + std::pow(bridge_dy, 2)); const int bridge_index_start = pillar_n * 2; Vector bridgePositions(bridge_n * 2); // Tablier for (int i = 0; i < bridge_n; i++) { const auto current_x = (float) (x_start + (i + 1) * bridge_dx); for (int j = 0; j < 2; j++) { const auto current_y = (float) (bridge_y - j * bridge_dy); Vector2f position(current_x, current_y); Particle p(position, Vector2f(0, 0), Vector2f(0, 0), 1.0f); particleSystem.addParticle(p); bridgePositions(i * 2 + j) = position; const int index = bridge_index_start + i * 2 + j; if (j == 0) { // Ressorts internes au tablier if (i > 0) { Spring s1(index - 1, index, 0, l0_diag); Spring s2(index - 2, index, 0, bridge_dx); particleSystem.addSpring(s1); particleSystem.addSpring(s2); } // Ressorts entre les extrémités du tablier et les piliers if (i == 0) { Spring s(pillar_connection_index, index, 0, bridge_dx); particleSystem.addSpring(s); } if (i == bridge_n - 1) { Spring s(pillar_n + pillar_connection_index, index, 0, bridge_dx); particleSystem.addSpring(s); } } else { Spring s1(index - 1, index, 0, bridge_dy); particleSystem.addSpring(s1); // Ressorts internes au tablier if (i > 0) { Spring s2(index - 2, index, 0, bridge_dx); particleSystem.addSpring(s2); } // Ressorts entre les extrémités du tablier et les piliers if (i == 0) { Spring s3(pillar_connection_index, index, 0, bridge_dx); particleSystem.addSpring(s3); Spring s4(pillar_connection_index - 1, index, 0, l0_diag); particleSystem.addSpring(s4); } if (i == bridge_n - 1) { Spring s3(pillar_n + pillar_connection_index, index, 0, bridge_dx); particleSystem.addSpring(s3); Spring s4(pillar_n + pillar_connection_index - 1, index, 0, l0_diag); particleSystem.addSpring(s4); } } } } // Toit du tablier en demi-cercle float r = bridge_n - 1; const float a = r / 2; for (int i = 0; i < bridge_n; i++) { // Calcul de y en utilisant l'équation du haut d'un cercle float y = (float) std::sqrt(std::pow(r, 2) - std::pow(i - a, 2)) + a; Vector position = Vector2f(x_start + bridge_dx + i * bridge_dx, bridge_dx + y * bridge_dx); Particle p(position, Vector2f(0, 0), Vector2f(0, 0), 1.0f); particleSystem.addParticle(p); // Ressorts entre les particules du toit int index = pillar_n * 2 + bridge_n * 2 + i; if (i > 0) { Spring s(index - 1, index, 0, bridge_dx); particleSystem.addSpring(s); } // Ressorts jusqu'au tablier int bridge_index = bridge_index_start + i * 2; float d = (position - bridgePositions(i * 2)).norm(); Spring s1(bridge_index, index, 0, d); particleSystem.addSpring(s1); // Ressorts liés aux piliers if (i == 1 || i == bridge_n - 2) { int p_index_m = i == 1 ? 1 : 2; int p_index = pillar_n * p_index_m - 1; Vector2f pp(i == 1 ? x_start : x_start + bridge_dx * (bridge_n + 1), y_start + pillar_n * pillar_dy); Spring p1(p_index, index, 0, (position - pp).norm() / 1.5f); particleSystem.addSpring(p1); } // Ressorts aux deux extrémités if (i == 0 || i == bridge_n - 1) { int p_index_m = i == 0 ? 0 : 1; int p_index = p_index_m * pillar_n + (pillar_n / 2); Spring p2(p_index, index, 0, (float) std::sqrt(std::pow(bridge_dx, 2) + std::pow(bridge_dy, 2))); particleSystem.addSpring(p2); } } } } ParticleSimApplication::ParticleSimApplication() : nanogui::Screen(nanogui::Vector2i(1280, 720), "GTI320 Labo Physique lineaire", true, false, true, true, false, 4, 1), m_particleSystem(), m_stepping(false), m_fpsCounter(0), m_fpsTime(0.0), m_maxIter(10), m_solverType(kGaussSeidel) { initGui(); createBeam(m_particleSystem, m_stiffness); // le modèle "poutre" est sélectionné à l'initialisation m_particleSystem.pack(m_p0, m_v0, m_f0); perform_layout(); reset(); } void ParticleSimApplication::initGui() { // Initialisation de la fenêtre m_window = new nanogui::Window(this, "Particle sim"); m_window->set_position(nanogui::Vector2i(8, 8)); m_window->set_layout(new nanogui::GroupLayout()); // initialisation du canvas où est affiché le système de particules m_canvas = new ParticleSimGLCanvas(this); m_canvas->set_background_color({255, 255, 255, 255}); m_canvas->set_size({1000, 600}); m_canvas->set_draw_border(false); // Initialisation de la fenêtre de contrôle nanogui::Window *controls = new nanogui::Window(this, "Controls"); controls->set_position(nanogui::Vector2i(960, 10)); controls->set_layout(new nanogui::GroupLayout()); Widget *tools = new Widget(controls); tools->set_layout(new nanogui::BoxLayout(nanogui::Orientation::Vertical, nanogui::Alignment::Middle, 0, 20)); // Intervalles des curseur const auto stiffnessMinMax = std::make_pair(0.0f, logf(5000.f)); const auto iterMinMax = std::make_pair(1.f, 100.f); // Affichage du FPS m_panelFPS = new Widget(tools); m_panelFPS->set_layout(new nanogui::BoxLayout(nanogui::Orientation::Horizontal, nanogui::Alignment::Middle, 0, 5)); m_labelFPS = new nanogui::Label(m_panelFPS, "FPS :"); m_textboxFPS = new nanogui::TextBox(m_panelFPS); m_textboxFPS->set_fixed_width(60); m_textboxFPS->set_value("0"); // Affichage du numéro de frame m_panelFrames = new Widget(tools); m_panelFrames->set_layout( new nanogui::BoxLayout(nanogui::Orientation::Horizontal, nanogui::Alignment::Middle, 0, 5)); m_labelFrames = new nanogui::Label(m_panelFrames, "Frame :"); m_textboxFrames = new nanogui::TextBox(m_panelFrames); m_textboxFrames->set_fixed_width(60); m_textboxFrames->set_value("0"); // Boutons pour le choix du solveur m_panelSolver = new Widget(tools); m_panelSolver->set_layout(new nanogui::BoxLayout(nanogui::Orientation::Vertical, nanogui::Alignment::Middle, 0, 5)); new nanogui::Label(m_panelSolver, "Solver : "); nanogui::Button *b = new nanogui::Button(m_panelSolver, "Gauss-Seidel"); b->set_flags(nanogui::Button::RadioButton); b->set_pushed(true); b->set_callback([this] { m_solverType = kGaussSeidel; }); b = new nanogui::Button(m_panelSolver, "Gauss-Seidel (coloration)"); b->set_callback([this] { m_solverType = kColorGaussSeidel; }); b->set_flags(nanogui::Button::RadioButton); b = new nanogui::Button(m_panelSolver, "Cholesky"); b->set_callback([this] { m_solverType = kCholesky; }); b->set_flags(nanogui::Button::RadioButton); b = new nanogui::Button(m_panelSolver, "None"); b->set_callback([this] { m_solverType = kNone; }); b->set_flags(nanogui::Button::RadioButton); // Curseur de rigidité Widget *panelSimControl = new Widget(tools); panelSimControl->set_layout( new nanogui::BoxLayout(nanogui::Orientation::Vertical, nanogui::Alignment::Middle, 0, 5)); m_panelStiffness = new Widget(panelSimControl); m_panelStiffness->set_layout( new nanogui::BoxLayout(nanogui::Orientation::Horizontal, nanogui::Alignment::Middle, 0, 5)); m_labelStiffness = new nanogui::Label(m_panelStiffness, "Stiffness : "); m_sliderStiffness = new nanogui::Slider(m_panelStiffness); m_sliderStiffness->set_range(stiffnessMinMax); m_textboxStiffness = new nanogui::TextBox(m_panelStiffness); m_sliderStiffness->set_callback([this](float value) { m_stiffness = std::exp(value); onStiffnessSliderChanged(); }); m_sliderStiffness->set_value(std::log(300.f)); // Curseur du nombre maximum d'itération pour Jacobi et Gauss-Seidel Widget *panelMaxIter = new Widget(panelSimControl); panelMaxIter->set_layout( new nanogui::BoxLayout(nanogui::Orientation::Horizontal, nanogui::Alignment::Middle, 0, 5)); new nanogui::Label(panelMaxIter, "Max iterations : "); nanogui::Slider *sliderMaxIter = new nanogui::Slider(panelMaxIter); sliderMaxIter->set_range(iterMinMax); nanogui::TextBox *textboxMaxIter = new nanogui::TextBox(panelMaxIter); textboxMaxIter->set_value(std::to_string(m_maxIter)); sliderMaxIter->set_value(m_maxIter); sliderMaxIter->set_callback([this, textboxMaxIter](float value) { m_maxIter = (int) value; textboxMaxIter->set_value(std::to_string(m_maxIter)); }); // Bouton «Simulate» nanogui::Button *startStopButton = new nanogui::Button(panelSimControl, "Simulate"); startStopButton->set_flags(nanogui::Button::ToggleButton); startStopButton->set_change_callback([this](bool val) { m_stepping = val; if (val) { m_prevTime = glfwGetTime(); draw_all(); } }); // Bouton «Step» nanogui::Button *stepButton = new nanogui::Button(panelSimControl, "Step"); stepButton->set_callback([this] { if (!m_stepping) step(deltaT); }); // Bouton «Reset» nanogui::Button *resetButton = new nanogui::Button(panelSimControl, "Reset"); resetButton->set_callback([this] { reset(); }); // Boutons pour le choix du modèle Widget *panelExamples = new Widget(tools); panelExamples->set_layout(new nanogui::BoxLayout(nanogui::Orientation::Vertical, nanogui::Alignment::Middle, 0, 5)); new nanogui::Label(panelExamples, "Examples : "); nanogui::Button *loadClothButton = new nanogui::Button(panelExamples, "Cloth"); loadClothButton->set_callback([this] { createHangingCloth(m_particleSystem, m_stiffness); m_particleSystem.pack(m_p0, m_v0, m_f0); reset(); }); nanogui::Button *loadLargeClothButton = new nanogui::Button(panelExamples, "Large cloth"); loadLargeClothButton->set_callback([this] { createLargeHangingCloth(m_particleSystem, m_sliderStiffness->value()); m_particleSystem.pack(m_p0, m_v0, m_f0); reset(); }); nanogui::Button *loadBeamButton = new nanogui::Button(panelExamples, "Beam"); loadBeamButton->set_callback([this] { createBeam(m_particleSystem, m_stiffness); m_particleSystem.pack(m_p0, m_v0, m_f0); reset(); }); nanogui::Button *loadRopeButton = new nanogui::Button(panelExamples, "Rope"); loadRopeButton->set_callback([this] { createHangingRope(m_particleSystem, m_stiffness); m_particleSystem.pack(m_p0, m_v0, m_f0); reset(); }); nanogui::Button *loadVotreExemple = new nanogui::Button(panelExamples, "Le vôtre"); loadVotreExemple->set_callback([this] { createVotreExemple(m_particleSystem, m_stiffness); m_particleSystem.pack(m_p0, m_v0, m_f0); reset(); }); } /** * Réaction aux événements déclenchés par le clavier */ bool ParticleSimApplication::keyboard_event(int key, int scancode, int action, int modifiers) { if (Screen::keyboard_event(key, scancode, action, modifiers)) return true; if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS) { set_visible(false); return true; } return false; } /** * Boucle principale * * Cette fonction est appelée périodiquement lorsque le programme est actif. * C'est ici que tout se passe. Si la simulation est en cours, la fonction * `step` est appelée pour faire avancer le système d'un intervalle de temps * DELTA_T. Ensuite, l'affichage est mis à jour. */ void ParticleSimApplication::draw_contents() { Screen::draw_contents(); if (m_stepping) { auto now = glfwGetTime(); float dt = now - m_prevTime; step(deltaT); // Update frames per second // m_fpsTime += dt; ++m_fpsCounter; if (m_fpsCounter > 30) { const float fps = (float) m_fpsCounter / m_fpsTime; char buf[64]; snprintf(buf, sizeof(buf), "%3.1f", fps); m_fpsCounter = 0; m_fpsTime = 0.0; m_textboxFPS->set_value(buf); } m_prevTime = now; updateFrameCounter(); } redraw(); } /** * Appelée lorsque le curseur de rigidité est modifié. La nouvelle rigidité est * affectée à tous les ressorts */ void ParticleSimApplication::onStiffnessSliderChanged() { // Update all springs with the slider value for (Spring &s: getParticleSystem().getSprings()) { s.k = m_stiffness; } char buf[16]; snprintf(buf, sizeof(buf), "%4.0f", m_stiffness); m_textboxStiffness->set_value(buf); } /** * Effectue un pas de simulation de taille dt. */ void ParticleSimApplication::step(float dt) { // Construction des matrices de masse et de rigidité // m_particleSystem.buildMassMatrix(m_M); m_particleSystem.buildDfDx(m_dfdx); // Calcul des forces actuelles sur chacune de sparticules // m_particleSystem.computeForces(); m_canvas->applyMouseSpring(); // Assemblage des vecteurs d'états. // m_particleSystem.pack(m_x, m_v, m_f); ////////////////////////////////////////////////////////////////////////////////// // Construire le système d'équation linéaire sous la forme `A*v_plus = b`. // la construction de A et b est donnée dans les diapos du Cours 8. // // Note : lors du calcul de b, NE PAS calculer `Mg + Kx` ce // calcul est inutilement coûteux. Pour être plus efficace, on utilise // directement le vecteur d'état m_f. // ////////////////////////////////////////////////////////////////////////////////// // Remarque : A et b sont déclarés `const` et ce n'est pas une erreur. C'est // pour vous forcer à optimiser votre code. // // Considérez les exemples suivants : // // # Version 1 // Matrix A; // constructeur par défaut // A = B + C; // operator+ on construit une matrice et elle est suite copiée avec operator= // // # Version 2 // Matrix A = B + C; // la matrice construite dans operator+ est la matrice A. // // Bilan : Version 1 utilise 2 constructeurs et 1 opérateur de copie // Version 2 utilise un seul constructeur et aucune copie ////////////////////////////////////////////////////////////////////////////////// // const Matrix A = m_M + -1.0f * std::pow(deltaT, 2.0f) * m_dfdx; const Vector b = deltaT * m_f + m_v; // Résolution du système d'équations `A*v_plus = b`. // Vector v_plus; Vector acc; // vecteur d'accélérations switch (m_solverType) { case kGaussSeidel: gaussSeidel(A, b, v_plus, m_maxIter); break; case kColorGaussSeidel: gaussSeidelColor(A, b, v_plus, m_graphColor.getPartitions(), m_maxIter); break; case kCholesky: cholesky(A, b, v_plus); break; default: case kNone: // N'utilise pas de solveur, il s'agit de l'implémentation naive de // l'intégration d'Euler. acc.resize(m_M.rows()); // vecteur d'accélérations for (int i = 0; i < m_M.rows(); ++i) acc(i) = (1.0 / m_M(i, i)) * m_f(i); v_plus = m_v + dt * acc; break; } // Mise à jour du vecteur d'état de position via l'intégration d'Euler // implicite. Les nouvelles position sont calculées à partir des position // actuelles m_x et des nouvelles vitesses v_plus. Les nouvelles positions // sont stockées directement dans le vecteur m_x. m_x = m_x + dt * v_plus; // Affecte les valeurs calculées dans le vecteurs d'états aux particules du // système m_particleSystem.unpack(m_x, v_plus); } /** * Réinitialisation du système de particules */ void ParticleSimApplication::reset() { m_frameCounter = 0; m_particleSystem.unpack(m_p0, m_v0); m_graphColor.color(m_particleSystem); onStiffnessSliderChanged(); } /** * Mise à jour du compteur de frames */ void ParticleSimApplication::updateFrameCounter() { ++m_frameCounter; char buf[16]; snprintf(buf, sizeof(buf), "%d", m_frameCounter); m_textboxFrames->set_value(buf); }