From 847964fd1e37ac345e13d71e364e4ff17bc1915a Mon Sep 17 00:00:00 2001 From: william Date: Sat, 9 Mar 2024 16:19:14 -0500 Subject: [PATCH] =?UTF-8?q?Solveurs,=20probl=C3=A8mes=20de=20perf?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- labo01/DenseStorage.h | 2 +- labo01/Operators.h | 2 +- labo_physique/ParticleSimApplication.cpp | 323 ++++++++++------------- labo_physique/ParticleSystem.cpp | 19 +- labo_physique/Solvers.h | 96 +++++-- labo_physique/Vector2d.h | 8 + 6 files changed, 252 insertions(+), 198 deletions(-) diff --git a/labo01/DenseStorage.h b/labo01/DenseStorage.h index c0394b6..22f6e8f 100644 --- a/labo01/DenseStorage.h +++ b/labo01/DenseStorage.h @@ -137,7 +137,7 @@ namespace gti320 { * Constructeur de copie */ DenseStorage(const DenseStorage &other) - : m_data(new _Scalar[m_size]), m_size(other.m_size) { + : m_data(new _Scalar[other.m_size]), m_size(other.m_size) { memcpy(m_data, other.m_data, m_size * sizeof(_Scalar)); } diff --git a/labo01/Operators.h b/labo01/Operators.h index 0075e2c..eff76bd 100644 --- a/labo01/Operators.h +++ b/labo01/Operators.h @@ -27,7 +27,7 @@ namespace gti320 { operator*(const Matrix<_Scalar, RowsA, ColsA, StorageA> &A, const Matrix<_Scalar, RowsB, ColsB, StorageB> &B) { assert(A.cols() == B.rows()); - auto result = Matrix<_Scalar, Dynamic, Dynamic>(A.rows(), B.cols()); + auto result = Matrix<_Scalar, RowsA, ColsB>(A.rows(), B.cols()); for (int col = 0; col < B.cols(); col++) { for (int row = 0; row < A.rows(); row++) { diff --git a/labo_physique/ParticleSimApplication.cpp b/labo_physique/ParticleSimApplication.cpp index cdacbbe..6ffbe00 100644 --- a/labo_physique/ParticleSimApplication.cpp +++ b/labo_physique/ParticleSimApplication.cpp @@ -27,16 +27,14 @@ using namespace gti320; -namespace -{ +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) - { + static inline void createHangingCloth(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 16; @@ -46,10 +44,8 @@ namespace const int dy = 32; int index = 0; - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { + 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; @@ -58,20 +54,17 @@ namespace if (j == (N - 1) && i == (N - 1)) particle.fixed = true; particleSystem.addParticle(particle); - if (i > 0) - { - Spring s(index - N, index, k, (float)dy); + 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); + 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)); + if (i > 0 && j > 0) { + Spring s(index - N - 1, index, k, std::sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s); } ++index; @@ -84,8 +77,7 @@ namespace /** * Crée un système masse-ressort qui simule un grand tissu suspendu */ - static inline void createLargeHangingCloth(ParticleSystem& particleSystem, float k) - { + static inline void createLargeHangingCloth(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 32; @@ -95,10 +87,8 @@ namespace const int dy = 16; int index = 0; - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { + 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; @@ -106,20 +96,17 @@ namespace 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); + 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); + 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)); + if (i > 0 && j > 0) { + Spring s(index - N - 1, index, k, std::sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s); } ++index; @@ -132,8 +119,7 @@ namespace * 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) - { + static inline void createHangingRope(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 20; @@ -141,17 +127,15 @@ namespace const int dx = 32; int index = 0; - for (int j = 0; j < N; ++j) - { + 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); + if (j > 0) { + Spring s(index - 1, index, k, (float) dx); particleSystem.addSpring(s); } ++index; @@ -162,8 +146,7 @@ namespace /** * Crée un système masse-ressort qui simule une poutre flexible */ - static inline void createBeam(ParticleSystem& particleSystem, float k) - { + static inline void createBeam(ParticleSystem &particleSystem, float k) { particleSystem.clear(); const int N = 20; @@ -173,8 +156,7 @@ namespace const int dy = 32; int index = 0; - for (int j = 0; j < N; ++j) - { + for (int j = 0; j < N; ++j) { const int x = x_start + j * dx; // Bottom particle @@ -182,11 +164,10 @@ namespace 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)); + 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); + Spring s2(index - 2, index, k, (float) dx); particleSystem.addSpring(s2); } @@ -199,13 +180,12 @@ namespace 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); + Spring s(index - 1, index, k, (float) dy); particleSystem.addSpring(s); - if (j > 0) - { - Spring s2(index - 2, index, k, (float)dx); + 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)); + Spring s3(index - 3, index, k, (float) sqrt((float) dx * dx + (float) dy * dy)); particleSystem.addSpring(s3); } ++index; @@ -218,8 +198,7 @@ namespace /** * TODO Créez votre propre exemple */ - static inline void createVotreExemple(ParticleSystem& particleSystem, float k) - { + static inline void createVotreExemple(ParticleSystem &particleSystem, float k) { particleSystem.clear(); // TODO Amusez-vous. Rendu ici, vous le méritez. @@ -229,9 +208,11 @@ namespace } -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) -{ +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 @@ -241,8 +222,7 @@ ParticleSimApplication::ParticleSimApplication() : nanogui::Screen(nanogui::Vect reset(); } -void ParticleSimApplication::initGui() -{ +void ParticleSimApplication::initGui() { // Initialisation de la fenêtre m_window = new nanogui::Window(this, "Particle sim"); m_window->set_position(nanogui::Vector2i(8, 8)); @@ -250,16 +230,16 @@ void ParticleSimApplication::initGui() // 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_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"); + 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); + Widget *tools = new Widget(controls); tools->set_layout(new nanogui::BoxLayout(nanogui::Orientation::Vertical, nanogui::Alignment::Middle, 0, 20)); // Intervalles des curseur @@ -276,7 +256,8 @@ void ParticleSimApplication::initGui() // 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_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); @@ -286,7 +267,7 @@ void ParticleSimApplication::initGui() 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"); + 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; }); @@ -301,115 +282,106 @@ void ParticleSimApplication::initGui() 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)); + 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_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_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)); + 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); + nanogui::Slider *sliderMaxIter = new nanogui::Slider(panelMaxIter); sliderMaxIter->set_range(iterMinMax); - nanogui::TextBox* textboxMaxIter = new nanogui::TextBox(panelMaxIter); + 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)); - }); + 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"); + 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(); - } - }); + 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); - }); + 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(); - }); + 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); + 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 *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 *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 *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(); - }); + 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) -{ +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) { @@ -427,12 +399,10 @@ bool ParticleSimApplication::keyboard_event(int key, int scancode, int action, i * `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() -{ +void ParticleSimApplication::draw_contents() { Screen::draw_contents(); - if (m_stepping) - { + if (m_stepping) { auto now = glfwGetTime(); float dt = now - m_prevTime; @@ -442,9 +412,8 @@ void ParticleSimApplication::draw_contents() // m_fpsTime += dt; ++m_fpsCounter; - if (m_fpsCounter > 30) - { - const float fps = (float)m_fpsCounter / m_fpsTime; + 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; @@ -463,11 +432,9 @@ void ParticleSimApplication::draw_contents() * Appelée lorsque le curseur de rigidité est modifié. La nouvelle rigidité est * affectée à tous les ressorts */ -void ParticleSimApplication::onStiffnessSliderChanged() -{ +void ParticleSimApplication::onStiffnessSliderChanged() { // Update all springs with the slider value - for (Spring& s : getParticleSystem().getSprings()) - { + for (Spring &s: getParticleSystem().getSprings()) { s.k = m_stiffness; } @@ -479,8 +446,7 @@ void ParticleSimApplication::onStiffnessSliderChanged() /** * Effectue un pas de simulation de taille dt. */ -void ParticleSimApplication::step(float dt) -{ +void ParticleSimApplication::step(float dt) { // Construction des matrices de masse et de rigidité // m_particleSystem.buildMassMatrix(m_M); @@ -520,33 +486,32 @@ void ParticleSimApplication::step(float dt) // Version 2 utilise un seul constructeur et aucune copie ////////////////////////////////////////////////////////////////////////////////// // - const Matrix A; - const Vector b; + 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; + 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; } // TODO Mise à jour du vecteur d'état de position via l'intégration d'Euler @@ -563,8 +528,7 @@ void ParticleSimApplication::step(float dt) /** * Réinitialisation du système de particules */ -void ParticleSimApplication::reset() -{ +void ParticleSimApplication::reset() { m_frameCounter = 0; m_particleSystem.unpack(m_p0, m_v0); m_graphColor.color(m_particleSystem); @@ -575,8 +539,7 @@ void ParticleSimApplication::reset() /** * Mise à jour du compteur de frames */ -void ParticleSimApplication::updateFrameCounter() -{ +void ParticleSimApplication::updateFrameCounter() { ++m_frameCounter; char buf[16]; snprintf(buf, sizeof(buf), "%d", m_frameCounter); diff --git a/labo_physique/ParticleSystem.cpp b/labo_physique/ParticleSystem.cpp index dd9d44b..a2171f7 100644 --- a/labo_physique/ParticleSystem.cpp +++ b/labo_physique/ParticleSystem.cpp @@ -95,7 +95,7 @@ void ParticleSystem::unpack(const Vector &_pos, /** - * Construction de la matirce de masses. + * Construction de la matrice de masses. */ void ParticleSystem::buildMassMatrix(Matrix &M) { const int numParticles = m_particles.size(); @@ -131,6 +131,9 @@ void ParticleSystem::buildDfDx(Matrix &dfdx) { dfdx.resize(dim, dim); dfdx.setZero(); + auto identity = Matrix(); + identity.setIdentity(); + // Pour chaque ressort... for (const Spring &spring: m_springs) { // TODO @@ -141,6 +144,20 @@ void ParticleSystem::buildDfDx(Matrix &dfdx) { // Astuce: créer une matrice de taille fixe 2 par 2 puis utiliser la classe SubMatrix pour accumuler // les modifications sur la diagonale (2 endroits) et pour mettre à jour les blocs non diagonale (2 endroits). + auto p0 = m_particles[spring.index0]; + auto p1 = m_particles[spring.index1]; + Vector distance = p1.x - p0.x; + Matrix distance_m = distance.as_matrix(); + Matrix distance_t = distance_m.transpose(); + float l = distance.norm(); + Matrix l2_m = std::pow(l, 2.0f) * identity; + float l3 = std::pow(l, 3.0f); + + Matrix dd = -1.0f * (distance_m * distance_t); + Matrix term1 = spring.k * identity; + Matrix term2 = -1.0f * (1 / l3) * spring.k * spring.l0 * (l2_m + dd); + + dfdx.block(spring.index0, spring.index1, 2, 2) = term1 + term2; } } diff --git a/labo_physique/Solvers.h b/labo_physique/Solvers.h index 4e2d84d..3815509 100644 --- a/labo_physique/Solvers.h +++ b/labo_physique/Solvers.h @@ -14,10 +14,11 @@ #include "Math3D.h" -namespace gti320 -{ +namespace gti320 { // Identification des solveurs - enum eSolverType { kNone, kGaussSeidel, kColorGaussSeidel, kCholesky }; + enum eSolverType { + kNone, kGaussSeidel, kColorGaussSeidel, kCholesky + }; // Paramètres de convergences pour les algorithmes itératifs static const float eps = 1e-4f; @@ -26,21 +27,54 @@ namespace gti320 /** * Résout Ax = b avec la méthode Gauss-Seidel */ - static void gaussSeidel(const Matrix& A, - const Vector& b, - Vector& x, int k_max) - { + static void gaussSeidel(const Matrix &A, + const Vector &b, + Vector &x, int k_max) { // TODO // // Implémenter la méthode de Gauss-Seidel + int n = b.size(); + + x.resize(n); + x.setZero(); + + bool converged = false; + int k = 0; + + do { + Vector nx = x; + + for (int i = 0; i < n; i++) { + nx(i) = b(i); + + for (int j = 0; j < i; j++) { + nx(i) = nx(i) - A(i, j) * nx(j); + } + + for (int j = i + 1; j < n; j++) { + nx(i) = nx(i) - A(i, j) * x(j); + } + + nx(i) = nx(i) / A(i, i); + } + + k++; + Vector r = A * x - b; + + converged = k >= k_max || + (nx - x).norm() / nx.norm() < tau || + r.norm() / b.norm() < eps; + + x = nx; + } while (!converged); } /** * Résout Ax = b avec la méthode Gauss-Seidel (coloration de graphe) */ - static void gaussSeidelColor(const Matrix& A, const Vector& b, Vector& x, const Partitions& P, const int maxIter) - { + static void gaussSeidelColor(const Matrix &A, const Vector &b, + Vector &x, const Partitions &P, const int maxIter) { // TODO // // Implémenter la méthode de Gauss-Seidel avec coloration de graphe. @@ -51,31 +85,63 @@ namespace gti320 /** * Résout Ax = b avec la méthode de Cholesky */ - static void cholesky(const Matrix& A, - const Vector& b, - Vector& x) - { + static void cholesky(const Matrix &A, + const Vector &b, + Vector &x) { + int n = A.rows(); + x.resize(n); + x.setZero(); + // TODO // // Calculer la matrice L de la factorisation de Cholesky + auto L = Matrix(n, n); + for (int j = 0; j < n; j++) { + for (int i = j; i < n; i++) { + float s = 0; + for (int k = 0; k < j; k++) { + s += L(i, k) * L(j, k); + } + + if (i == j) { + L(i, i) = std::sqrt(A(i, i) - s); + } else { + L(i, j) = (A(i, j) - s) / L(j, j); + } + } + } // TODO // // Résoudre Ly = b + auto y = Vector(n); + for (int i = 0; i < n; i++) { + y(i) = b(i); + for (int j = 0; j < i; j++) { + y(i) -= L(i, j) * y(j); + } + y(i) /= L(i, i); + } // TODO // // Résoudre L^t x = y // - // Remarque : ne pas caculer la transposer de L, c'est inutilement + // Remarque : ne pas calculer la transposée de L, c'est inutilement // coûteux. + for (int i = n - 1; i >= 0; i--) { + x(i) = y(i); + for (int j = i + 1; j < n; j++) { + x(i) -= L(j, i) * x(j); + } - + x(i) /= L(i, i); + } } } diff --git a/labo_physique/Vector2d.h b/labo_physique/Vector2d.h index 94745d9..5f16e1e 100644 --- a/labo_physique/Vector2d.h +++ b/labo_physique/Vector2d.h @@ -11,6 +11,7 @@ * */ +#include "Matrix.h" #include "Math3D.h" namespace gti320 @@ -102,6 +103,13 @@ namespace gti320 { return sqrt(dot(*this)); } + + Matrix<_Scalar, 2, 1> as_matrix() const { + Matrix<_Scalar, 2, 1> mat; + mat(0, 0) = (*this)(0); + mat(1, 0) = (*this)(1); + return mat; + } }; typedef Vector Vector2f;