Simulation complète et fonctionnelle

labo_physique/GraphColoring.cpp

@@ -4,6 +4,8 @@
 using namespace gti320;
 
 void GraphColoring::color(ParticleSystem &particleSystem) {
+    m_partitions.clear();
+
     // La palette de couleurs
     ColorList C;
 
@@ -28,9 +30,10 @@ void GraphColoring::color(ParticleSystem &particleSystem) {
         p.color = color;
 
         if (m_partitions.size() <= color) {
-            m_partitions.push_back(std::vector<int>());
+            m_partitions.emplace_back();
         }
-        m_partitions[color].push_back(i);
+        m_partitions[color].push_back(i * 2);
+        m_partitions[color].push_back(i * 2 + 1);
         i++;
     }
 }
@@ -38,8 +41,9 @@ void GraphColoring::color(ParticleSystem &particleSystem) {
 int
 GraphColoring::findColor(const Particle &p, const std::vector<Particle> &particles, const std::vector<Spring> &springs,
                          ColorList &C) const {
-    size_t n = C.size();
+    int n = (int) C.size();
-    auto count = new int[n];
+    int count[n];
+    memset(count, 0, sizeof(int) * n);
 
     // TODO Trouver la premi<6D>re couleur de la palette C qui n'est pas attribu<62>e <20> une particule voisine.
     //      Si une couleur est introuvable, ajouter une nouvelle couleur <20> la palette et retournez la couleur.
@@ -56,6 +60,7 @@ GraphColoring::findColor(const Particle &p, const std::vector<Particle> &particl
         }
     }
 
+    C.push_back(n);
     return n;
 }
 

labo_physique/ParticleSimApplication.cpp

@@ -202,7 +202,50 @@ namespace {
         particleSystem.clear();
 
         // TODO Amusez-vous. Rendu ici, vous le méritez.
+        const int N = 20;
+        const int x_start = 200;
+        const int y_start = 400;
+        const int dx = 32;
+        const int dy = 32;
 
+        Particle p1(Vector2f(x_start + dx * 2, y_start + dy * 4), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p1);
+
+        Particle p2(Vector2f(x_start, y_start), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p2);
+
+        Particle p3(Vector2f(x_start + dx * 4, y_start), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p3);
+
+        Particle p4(Vector2f(x_start + dx * 2, y_start + dy * 2), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p4);
+
+        Particle p5(Vector2f(x_start + dx, y_start + dy), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p5);
+
+        Particle p6(Vector2f(x_start + dx * 3, y_start + dy), Vector2f(0, 0), Vector2f(0, 0), 1.0);
+        particleSystem.addParticle(p6);
+
+        Spring s1(0, 1, 0, dx);
+        particleSystem.addSpring(s1);
+        Spring s2(0, 2, 0, dx);
+        particleSystem.addSpring(s2);
+        Spring s3(1, 2, 0, dx);
+        particleSystem.addSpring(s3);
+
+        Spring s4(3, 4, 0, dx);
+        particleSystem.addSpring(s4);
+        Spring s5(3, 5, 0, dx);
+        particleSystem.addSpring(s5);
+        Spring s6(4, 5, 0, dx);
+        particleSystem.addSpring(s6);
+
+        Spring s7(1, 4, 0, dx);
+        particleSystem.addSpring(s7);
+        Spring s8(0, 3, 0, dx);
+        particleSystem.addSpring(s8);
+        Spring s9(2, 5, 0, dx);
+        particleSystem.addSpring(s9);
     }
 
 

labo_physique/ParticleSystem.cpp

@@ -147,17 +147,19 @@ void ParticleSystem::buildDfDx(Matrix<float, Dynamic, Dynamic> &dfdx) {
         auto p0 = m_particles[spring.index0];
         auto p1 = m_particles[spring.index1];
 
-        Vector<float, 2> distance = p1.x - p0.x;
+        float l = (p1.x - p0.x).norm();
-        Matrix<float, 2, 1> distance_m = distance.as_matrix();
+        float a = spring.k * (1 - spring.l0 / l);
-        Matrix<float, 1, 2> distance_t = distance_m.transpose<float, 1, 2, ColumnStorage>();
-        float l = distance.norm();
-        Matrix<float, 2, 2> l2_m = std::pow(l, 2.0f) * identity;
-        float l3 = std::pow(l, 3.0f);
 
-        Matrix<float, 2, 2> dd = -1.0f * (distance_m * distance_t);
+        Matrix<float, 2, 2> diag;
-        Matrix<float, 2, 2> term1 = spring.k * identity;
+        diag.setZero();
-        Matrix<float, 2, 2> term2 = -1.0f * (1 / l3) * spring.k * spring.l0 * (l2_m + dd);
+        diag(0, 0) = -a;
+        diag(1, 1) = -a;
 
-        dfdx.block(spring.index0, spring.index1, 2, 2) = term1 + term2;
+        Matrix<float, 2, 2> ndiag = -1.0f * diag;
+
+        dfdx.block(spring.index0, spring.index1, 2, 2) = diag;
+        dfdx.block(spring.index0 + 2, spring.index1 + 2, 2, 2) = diag;
+        dfdx.block(spring.index0 + 2, spring.index1, 2, 2) = ndiag;
+        dfdx.block(spring.index0, spring.index1 + 2, 2, 2) = ndiag;
     }
 }

labo_physique/Solvers.h

@@ -86,12 +86,13 @@ namespace gti320 {
 
         bool converged = false;
         int k = 0;
+        Vector<float, Dynamic> nx;
 
         do {
-            Vector<float, Dynamic> nx = x;
+            nx = x;
             for (const auto &indices: P) {
 #pragma omp parallel for
-                for (int i = 0; i < indices.size(); i++) {
+                for (const int &i: indices) {
                     nx(i) = b(i);
 
                     for (int j = 0; j < i; j++) {
@@ -99,7 +100,7 @@ namespace gti320 {
                     }
 
                     for (int j = i + 1; j < n; j++) {
-                        nx(i) = nx(i) - A(i, j) * x(j);
+                        nx(i) = nx(i) - A(i, j) * nx(j);
                     }
 
                     nx(i) = nx(i) / A(i, i);