Discussion :

[guppyworld]

Bonjour,

Je suis actuellement entrain d'implémenter un système de particules dans un environement en 2D. Pour le moment le programme peut traiter environ 2 000 particules avant de commencer a ramer. Je commence a être a cours d'idée sur les optimisations que je peut faire pour augmenter le nombre de particules. J'ai donc plusieurs questions :

- Le conteneur choisis pour les particules vous semble t-il correct ?
- La répartition des particules dans une grille est-elle correct ?
- Est que une ligne de code vous parait bizarre ?

En vous remerciant d'avance pour vos aides

Le code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
void ParticleSystem::computeParticles(Map *pMap)
{
    int widthCell = 32;
    int heightCell = 32;
 
    int nbCellWidth = pMap->WIDTH * 32 / widthCell + 1;
    int nbCellHeight = pMap->HEIGHT * 32 / heightCell + 1;
 
     std::vector<int> grid[nbCellHeight][nbCellWidth];
 
    for(int i = 0; i<mParticles.size(); i++){
 
        int x = mParticles[i]->mPositionY/heightCell;
        int y = mParticles[i]->mPositionX/widthCell;
 
        grid[x][y].push_back(i);
 
    }
 
    // trouver les particules voisines
    for(int i = 0; i<nbCellHeight; i++)
        for(int j = 0; j<nbCellWidth; j++)
        {
            while(!grid[i][j].empty()){
                int particule = grid[i][j].back(); // remove particle par la fin pour O(1)
                grid[i][j].pop_back();
 
                for(int x = -1; x<=1; x++){
                    for(int y = -1; y<=1; y++){
                        int a = i + x;
                        int b = j + y;
 
                        if(a < 0 || a >=nbCellHeight || b < 0 || b >= nbCellWidth)
                            continue;
 
                        for(int z = 0; z<grid[a][b].size(); z++)
                        {
                            int voisin = grid[a][b][z];
 
                            float dx = mParticles[voisin]->mPositionX - mParticles[particule]->mPositionX;
                            float dy = mParticles[voisin]->mPositionY - mParticles[particule]->mPositionY;
                            float distance2 = dx * dx + dy * dy;
                            if(distance2 < mParticleWidth * mParticleWidth){
                                float distance = sqrt(distance2);
                                float weight = (1.0f-distance/mParticleWidth);
                                ParticleContact contact(particule, voisin, distance, weight);
                                mParticleContacts.push_back(contact);
 
                            }
                        }
                    }
                }
                mParticles[particule]->mDensity=0.0f;
                mParticles[particule]->mPressure=0.0f;
                mParticles[particule]->mVelocityY -= GRAVITY;
            }
        }
 
    // compute density
    for(int i = 0; i<mParticleContacts.size(); i++){
 
        int p1 = mParticleContacts[i].mParticule1;
        int p2 = mParticleContacts[i].mParticule2;
 
        mParticles[p1]->mDensity += mParticleContacts[i].mWeight;
        mParticles[p2]->mDensity += mParticleContacts[i].mWeight;
 
    }
 
    while(!mParticleContacts.empty()){
 
        ParticleContact contact = mParticleContacts.back();
        mParticleContacts.pop_back();
 
        int p1 = contact.mParticule1;
        int p2 = contact.mParticule2;
 
        //compute pressure
        mParticles[p1]->mPressure = std::max(0.0f, n * (mParticles[p1]->mDensity - w0));
        mParticles[p2]->mPressure = std::max(0.0f, n * (mParticles[p2]->mDensity - w0));
 
        //compute force
        float nX = (mParticles[p2]->mPositionX - mParticles[p1]->mPositionX)/contact.mDistance;
        float nY = (mParticles[p2]->mPositionY - mParticles[p1]->mPositionY)/contact.mDistance;
 
        mParticles[p1]->mVelocityX += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * nX;
        mParticles[p1]->mVelocityY += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * nY;
 
        mParticles[p2]->mVelocityX += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * -nX;
        mParticles[p2]->mVelocityY += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * -nY;
 
        mParticles[p1]->mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p2]->mVelocityX - mParticles[p1]->mVelocityX);
        mParticles[p1]->mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p2]->mVelocityY - mParticles[p1]->mVelocityY);
 
        mParticles[p2]->mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p1]->mVelocityX - mParticles[p2]->mVelocityX);
        mParticles[p2]->mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p1]->mVelocityY - mParticles[p2]->mVelocityY);
 
    }
 
    // update position
    for (int i = 0; i < mParticles.size(); ++i) {
        mParticles[i]->updatePosition(pMap ,TIME_STEP);
    }
}

Les conteneurs dans ParticleSystem.h

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
std::vector<std::shared_ptr<Particle>> mParticles;
        std::vector<ParticleContact> mParticleContacts;

et la classe ParticleContact

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
class ParticleContact
{
    public:
        ParticleContact() {}
        ParticleContact(int p1, int p2, float distance, float weight) {
            mParticule1 = p1;
            mParticule2 = p2;
            mDistance = distance;
            mWeight = weight;
 
        }
        virtual ~ParticleContact() {}
 
        int mParticule1;
        int mParticule2;
        float mDistance;
        float mWeight;
};

[ternel]

Première réaction: Pourquoi un vecteur de shared_ptr? Ca te tue les performances car les données ne sont plus contigües en mémoire.

[oodini]

Je n'ai pas le temps de rentrer dans le détail de ton code, notamment l'algorithmique, mais utiliser des push_back, pour les perfs, ce n'est pas génial.
Pour ce genre de trucs, il vaut mieux coder dans le style C. Cela te permettra en outre de plus facilement faire un portage en CUDA/OpenCL, qui sont des technos bien adaptés à ce genre de problèmes.

[dragonjoker59]

Salut !

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

std::vector<int> grid[nbCellHeight][nbCellWidth];

Ce conteneur n'est pas cache friendly, n'est pas valide (on ne peut pas déclarer un conteneur statique avec des tailles dynamiques), et comme je suppose que ta map ne doit pas changer de taille souvent, tu devrais le mettre en membre de la classe, pour l'initialiser une fois pour toutes, plutôt que de faire les allocations/désallocations qu'il implique, à chaque calcul de tes particules.
Je l'écrirais comme suit (pour l'instant) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
Map const & mMap;
using Particles = std::vector< int >;
std::vector< Particles > mGrid;

et dans le constructeur de ParticleSystem :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
ParticleSystem::ParticleSystem( Map const & map )
  : mGrid( map.WIDTH * map.HEIGHT, Particles() )
  , mMap( map )
{
}

Il faudra donc, pour te faciliter la vie, écrire les accesseurs te permettant d'accéder à l'élément aux coordonnées X et Y :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
Particles & getParticles( size_t x, size_t y )
{
  return mGrid[x * mMap.HEIGHT + y );
}

Ensuite, si tu connais le nombre de particules max que peut contenir Particles, tu peux le définir comme suit, afin de gagner encore au niveau du parcours du tableau :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
size_t constexpr MaxParticles = 12; //par exemple
using Particles = std::array< int, MaxParticles >;

Ca te permettrait de n'avoir qu'une seule allocation en création du ParticleSystem, le reste n'étant que du parcours.
De plus, cette linéarisation des données te permettrait de supprimer la double boucle au niveau de la recherche des voisins (les double boucles c'est mauvais pour les perfs) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
for ( auto & particles : mGrid )
{
  // Recherche voisins...
}

Bon, mais tout ça, ça montre une chose : tu as mal choisi tes structures de données.
Parce que ta recherche de voisins consiste en 6 BOUCLES IMBRIQUEES et si 2 boucles imbriquées c'est mal, je te laisse imaginer avec 6.
Ton soucis de perfs vient exactement de là, car au pire des cas, tu te retrouves avec une complexité de O(N^6) .

[guppyworld]

Première réaction: Pourquoi un vecteur de shared_ptr? Ca te tue les performances car les données ne sont plus contigües en mémoire.

Oui c'est en faisant des test j'ai oublié de remettre comme avant. Donc pour le moment ce sont des unique_ptr.

dragonjoker59:

Donc en faite si j'ai bien compris tu veut que pour la grille contenant l'indice des particules, je passe d'un tableau 2D à un tableau 1D ?

Et qu'il est préférable d'allouer directement un tableau pouvant contenir 10 000 particules par exemple plutôt que de faire un tableau dynamique ?

Parce que ta recherche de voisins consiste en 6 BOUCLES IMBRIQUEES et si 2 boucles imbriquées c'est mal, je te laisse imaginer avec 6.
Ton soucis de perfs vient exactement de là, car au pire des cas, tu te retrouves avec une complexité de O(N^6) .

La j'ai un doute sur le O(N^6) car les 2 premiers for ne dépendent pas du nombre de particules mais plutôt d'un nombre constant. Pareil pour les 2 for avec x et y qui sont également constant.

[Bousk]

Des unique_ptr fragmentent toujours ta mémoire.

[ternel]

En fait, la technique est même plutot recommandée, parce que tu limites les collisions aux collisions d'un même sous-espace, chacun contenant en général beaucoup moins de particules.
Dans un espace quadrillés en B blocs, pour N particules, tu ne calcules que B * (N/B)² = N²/B au lieu de N². À condition d'avoir une répartition uniforme dans la grille.

C'est en partant sur un tel raisonnement que tu arriveras sur le quad-tree (et son cousin en 3D, l'oct-tree), qui a pour vocation d'assurer une bonne distribution dans la grille.

Ta grille peut tout à fait être linéarisée.

D'un point de vue algorithmique, la boucle pourrait être pour toute particule p, pour toute autre particule q de la meme case, calculer l'effet de la collision de q sur la p.

Comme "pour toute particule p" est mathématiquement équivalent à "pour toute case c, pour toute particule p de c", c'est ce que tu as fait.

Par contre, tes calculs m'ont l'air très lourds, et surtout, ton assignation de case.
Je ne pense pas qu'il soit nécessaire de faire des allocations (et réallocation) à chaque pas de temps.

Avec un peu de "mise en cache", ta particule pourrait ressembler à

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
struct particule {
   coord_type coord; //une structure {x, y}
   vector_type vitesse; // une structure {dx, dy}
   vector_type acceleration;
   masse_type masse;
};

La raison pour laquelle je différencie systématiquement coord_type de vector_type, c'est que les opérateurs n'ont pas la même signification. Notamment coord_type + coord_type n'a pas de sens.

Et à chaque pas de temps, ton code devrait globalement être:

pour chaque pas de temps
    pour chaque particule,
        calculer la nouvelle accélération
            remise à 0
            calcul des collisions
        appliquer l'accélération,
        puis la vitesse.

éventuellement, tu peux stocker deux vector/array de particules, au début identiques. A chaque pas de temps, tu utilises l'un comme le temps actuel, et l'autre comme le temps suivant. Et au pas suivant, tu réutilises l'ancien temps actuel comme nouveau temps futur.
Il suffit de prendre les deux listes en références sur ta fonction effectuant un pas de temps.

[ternel]

Par contre, quel symptome qualifies-tu de "ramer"?

Combien as-tu de mémoire vive d'installée sur ton pc?

Normalement, avec un vector sans pointeurs, tu devrais avoir une consommation mémoire de plus ou moins une cacahuète (les quelques constantes de ton programme) + une cacahuète (la pile pour atteindre ta fonction de pas de temps) + sizeof <Particule> * N octets + N² * sizeof<Collision>.

Et cette dernière partie peut couter cher. N² sizeof collision, ca peut rapidement monter si tes particules sont dans un espace petit pour leur vitesse.
Il n'est pas nécessaire de stocker chaque collision, elle peuvent s'additionner au fur et à mesure.

[guppyworld]

Par contre, tes calculs m'ont l'air très lourds, et surtout, ton assignation de case.

Tu veux parler des lignes de code suivante ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
for(int i = 0; i<mParticles.size(); i++){
 
        int x = mParticles[i]->mPositionY/heightCell;
        int y = mParticles[i]->mPositionX/widthCell;
 
        grid[x][y].push_back(i);
 
    }

Sinon pour le moment je ne compte pas utiliser le QuadTree, je reste sur une carte divisé en case de 32 pixels sur 32.

Par contre je sais pas si vous avez compris mais en faite la grille contient les indices des particules pour chaque cellule. Je ne vois donc pas trop comment je peut faire autrement que la structure

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 std::vector<int> mGrid[NB_CELL_HEIGHT][NB_CELL_WIDTH];

. Ou au mieux je peut faire

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

 std::vector<int> mGrid[NB_CELL_HEIGHT * NB_CELL_WIDTH];

.

[moldavi]

Bonjour.

Bonjour,

Je suis actuellement entrain d'implémenter un système de particules dans un environnement en 2D. Pour le moment le programme peut traiter environ 2 000 particules avant de commencer a ramer. Je commence a être a cours d'idée sur les optimisations que je peut faire pour augmenter le nombre de particules.

Je vais juste poser une question : "le programme peut traiter environ 2 000 particules", d'accord, mais quelle est l'objectif ? 2001 particules ? 3000 particules ? 1 millions de particules ? etc...

Je vois plein de messages pour tenter d'optimiser ton code, mais si l'objectif est 1 milliards de particules avec un Pentium 2 et 100 Mo de RAM, c'est peine perdue...

Première étape, définir le nombre de particules à atteindre. Deuxième étape, le matériel est-il capable de prendre en charge l'objectif (processeur/ram/cg/etc...). Dernière étape, optimiser le code.

Bref, vous commencez tous par la dernière étape...

PS : j'abrège le nombre des étapes, il y a en a peut-être plus, selon la réponse à la question.

[ternel]

Théoriquement, un modèle de collision doit au plus gérer N² calculs de collisions.
La mémoire requise est N * sizeof(particules) + sizeof(map)
avec 4Go de ram, on peut raisonnablement avoir 80 octets par particules (10 entiers 64bits) et donc 50000 particules.
Une itération prend alors au plus 2,5 milliards de calculs de collisions, soit à 10000 cycles processeur l'un, 25 secondes à 1Ghz, et 6s à 4Ghz.
C'est long, mais pas tant que ça.

Si on fait une optimisation qui divise par 10 le nombre de collisions à calculer, c'est déjà très intéressant. (moins d'une seconde par itération)
Si on optimise les accès mémoire pour réduire les défauts de cache, on peut gagner beaucoup aussi.

L'objectif semble être l'affichage en temps réel.
Il faut passer en dessous de 50 milli-secondes pour avoir un résultat visuellement satisfaisant.
A noter qu'il faut aussi prendre en compte le temps de dessin, mais là, on est large: 50000 cercles, c'est vite géré.

D'un point de vue technique, il faut absolument séparer la structure Particule de son évolution, ne pas mélanger les variables d'état (position, vitesse, masse) des variables de calcul comme la force appliquée à l'itération actuelle

[guppyworld]

Je pensez atteindre au moins les 20k particules.

Pour le moment, je divise la map en grille de 16 par 16 avec des particules qui ont un diamètre de 10. Et chaque particule test si elle est en collision uniquement avec les particules des cases voisines. Le problème est que cette fonction prend trop de temps.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
void ParticleSystem::findContacts()
{
    for(int i=0; i<mCountParticles; i++)
    {
        int caseX = mParticles[i].mPositionX/16;
        int caseY = mParticles[i].mPositionY/16;
 
 
        // same case
        for(int k=0;k<mGrid[caseY][caseX].size();k++){
            int voisin = mGrid[caseY][caseX][k];
            float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
            float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
            if(dx > 0 || (dx==0 && dy>0)){
                float distance2 = dx * dx + dy * dy;
                if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                    float distance = sqrt(distance2);
                    float weight = (1.0f-distance/WIDTH_PARTICLE);
                    mParticleContacts.emplace_back(i, voisin, distance, weight);
                    mParticles[i].mDensity += weight;
                    mParticles[voisin].mDensity += weight;
                }
            }
        }
 
        // right case
        if(caseX+1<199){
            for(int k=0;k<mGrid[caseY][caseX+1].size();k++){
                int voisin = mGrid[caseY][caseX+1][k];
                float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
                float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
                float distance2 = dx * dx + dy * dy;
                if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                    float distance = sqrt(distance2);
                    float weight = (1.0f-distance/WIDTH_PARTICLE);
                    mParticleContacts.emplace_back(i, voisin, distance, weight);
                    mParticles[i].mDensity += weight;
                    mParticles[voisin].mDensity += weight;
                }
            }
        }
 
        //bottom cases
        for(int j=-1;j<2;j++){
            if(caseY-1>0 && caseX+j>0 && caseX + j<199){
                for(int k=0;k<mGrid[caseY-1][caseX+j].size();k++){
                    int voisin = mGrid[caseY-1][caseX+j][k];
                    float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
                    float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
                    float distance2 = dx * dx + dy * dy;
                    if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                        float distance = sqrt(distance2);
                        float weight = (1.0f-distance/WIDTH_PARTICLE);
                        mParticleContacts.emplace_back(i, voisin, distance, weight);
                        mParticles[i].mDensity += weight;
                        mParticles[voisin].mDensity += weight;
                    }
                }
            }
        }
        //apply gravity
        mParticles[i].mVelocityY -= GRAVITY;
        //compute pressure
        mParticles[i].mPressure = std::max(0.0f, n * (mParticles[i].mDensity - w0));
 
    }
}

[Kannagi]

On considérant que ton algo est correct et optimisé et que t'as mis le -O3.

Quand on vise les performances ,Il faut évité les appel de fonction et méthode , ce qui n'est pas ton cas.
Je pense que si tu arrive a mettre un tableau a la place : mGrid[caseY-1][caseX+j].size() , ça optimiserait pas mal :p
De plus que les tableaux de 'structure' qui l'obligera au compilateur forcément a faire un calcul supplémentaire (pour pointer la bonne adresse) que a un simple tableau.

Sinon tu peux gérer tout cela via ton GPU avec les shader x)

[koala01]

Salut,

As tu remarqué que le gros du temps (plus du tiers, vu qu'on est à 36.45%) se perd dans ParticleSystem::findContacts alors qu'elle n'est appelée que 2800 fois contre plusieurs millions d'appels pour certaines fonctions

Je note aussi que l'on a le même nombre d'appels (pour 1.73% du temps) à ParticlSystem::computeForces et Particle::updatePosition (pour 0.73% du temps).

J'en déduis donc que tu dois avoir "quelque part" une boucle qui contiendra un appel à ces trois fonctions.

Si tu veux améliorer les choses, il va falloir partir de cette fonction "bouffeuse de temps" et voir ce qu'elle fait afin d'améliorer son rendement.

Pourrais tu donc commencer par nous donner le code complet:

• de ta classe ParticleSystem (la définition complète de ta classe, telle qu'elle se trouve dans le fichier hpp qui lui est dévolu)
• l'implémentation de la fonction ParticlSystem::computeForces()
• des fonctions de ParticleSystem qui seraient éventuellement appelées par ParticlSystem::computeForces()

Il nous faudra surement d'autres précisions par la suite, mais ce sera déjà un bon début

[guppyworld]

Je pense que si tu arrive a mettre un tableau a la place : mGrid[caseY-1][caseX+j].size() , ça optimiserait pas mal :p

J'ai le même résultat avec un tableau simple

Mon super code

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#ifndef PARTICLESYSTEM_H
#define PARTICLESYSTEM_H
 
#include <vector>
#include <deque>
#include <algorithm>
#include <iostream>
#include <stdlib.h>
#include <list>
#include <memory>
#include <math.h>
#include "Particle.h"
#include "ParticleContact.h"
#include "Window.h"
#include "Map.h"
#include "Camera.h"
#include "QuadTree.h"
 
class ParticleSystem
{
    public:
        ParticleSystem();
        virtual ~ParticleSystem();
 
        void displayParticle(Window *window, Camera *camera);
        void createParticle(int type, float positionX, float positionY, float velocityX, float velocityY);
        void createParticles(int type, float positionX, float positionY);
 
        void computeParticles(Map *pMap);
        void createGrid();
        void findContacts();
        void computeForces();
        void updatePosition(Map *pMap);
 
 
    protected:
 
    private:
 
        const float PI = 3.14159f;
        const float TIME_STEP = 0.1f;
        const float GRAVITY = 0.5;
        const float n = 5.0f;
        const float w0 = -0.0f;
        const float a = -45.0f;
        const float VISCOSITY = 0.1f;
        const static int MAX_PARTICLES = 8192;
 
        const static int WIDTH_PARTICLE = 10;
        const static int WIDTH_PARTICLE_2 = 100;
        const static int RADIUS_PARTICLE = 5;
 
        Particle mParticles[MAX_PARTICLES];
 
        std::vector<ParticleContact> mParticleContacts;
        std::vector<int> mGrid[96][199];
 
        int mCountParticles;
};
 
#endif // PARTICLESYSTEM_H

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#include "ParticleSystem.h"
 
ParticleSystem::ParticleSystem()
{
    mCountParticles = 0;
}
 
ParticleSystem::~ParticleSystem()
{
 
}
 
void ParticleSystem::displayParticle(Window *window, Camera *camera)
{
    for(int i=0; i<mCountParticles;i++){
        int posXwindow = window->getWidth()/2-(camera->getPosX() - mParticles[i].mPositionX);
        int posYwindow = window->getHeight()/2 + (camera->getPosY() - mParticles[i].mPositionY);
 
 
        window->renderImage("circle", posXwindow-RADIUS_PARTICLE, posYwindow-RADIUS_PARTICLE, WIDTH_PARTICLE, WIDTH_PARTICLE);
    }
 
    window->displayNumber(400, 400, mCountParticles, 0,400,0);
}
 
void ParticleSystem::createParticle(int type, float positionX, float positionY, float velocityX, float velocityY)
{
    mParticles[mCountParticles] = Particle(type, positionX, positionY, velocityX, velocityY);
    mCountParticles ++;
}
 
void ParticleSystem::createParticles(int type, float positionX, float positionY)
{
    float space = 11;
 
    for(int i=0; i< 30; i++)
        for(int j=0; j<30; j++)
            createParticle(type, positionX+space*i, positionY-j*space, 0, 0);
}
 
void ParticleSystem::computeParticles(Map *pMap)
{
    createGrid();
    findContacts();
    computeForces();
    updatePosition(pMap);
}
 
void ParticleSystem::createGrid()
{
    for(int caseX=0; caseX<199; caseX++)
        for(int caseY=0; caseY<96; caseY++)
            mGrid[caseY][caseX].clear();
 
    for(int i = 0; i<mCountParticles; i++){
        int x = mParticles[i].mPositionX/10;
        int y = mParticles[i].mPositionY/10;
        mGrid[y][x].push_back(i);
    }
}
 
void ParticleSystem::findContacts()
{
    for(int i=0; i<mCountParticles; i++)
    {
        int caseX = mParticles[i].mPositionX/10;
        int caseY = mParticles[i].mPositionY/10;
 
 
        // same case
        for(int k=0;k<mGrid[caseY][caseX].size();k++){
            int voisin = mGrid[caseY][caseX][k];
            float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
            float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
            if(dx > 0 || (dx==0 && dy>0)){
                float distance2 = dx * dx + dy * dy;
                if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                    float distance = sqrt(distance2);
                    float weight = (1.0f-distance/WIDTH_PARTICLE);
                    mParticleContacts.emplace_back(i, voisin, distance, weight);
                    mParticles[i].mDensity += weight;
                    mParticles[voisin].mDensity += weight;
                }
            }
        }
 
        // right case
        if(caseX+1<199){
            for(int k=0;k<mGrid[caseY][caseX+1].size();k++){
                int voisin = mGrid[caseY][caseX+1][k];
                float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
                float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
                float distance2 = dx * dx + dy * dy;
                if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                    float distance = sqrt(distance2);
                    float weight = (1.0f-distance/WIDTH_PARTICLE);
                    mParticleContacts.emplace_back(i, voisin, distance, weight);
                    mParticles[i].mDensity += weight;
                    mParticles[voisin].mDensity += weight;
                }
            }
        }
 
        //bottom cases
        for(int j=-1;j<2;j++){
            if(caseY-1>0 && caseX+j>0 && caseX + j<199){
                for(int k=0;k<mGrid[caseY-1][caseX+j].size();k++){
                    int voisin = mGrid[caseY-1][caseX+j][k];
                    float dx = mParticles[voisin].mPositionX - mParticles[i].mPositionX;
                    float dy = mParticles[voisin].mPositionY - mParticles[i].mPositionY;
                    float distance2 = dx * dx + dy * dy;
                    if(distance2 < WIDTH_PARTICLE_2 && distance2 != 0){
                        float distance = sqrt(distance2);
                        float weight = (1.0f-distance/WIDTH_PARTICLE);
                        mParticleContacts.emplace_back(i, voisin, distance, weight);
                        mParticles[i].mDensity += weight;
                        mParticles[voisin].mDensity += weight;
                    }
                }
            }
        }
 
        //apply gravity
        mParticles[i].mVelocityY -= GRAVITY;
        //compute pressure
        mParticles[i].mPressure = std::max(0.0f, n * (mParticles[i].mDensity - w0));
 
    }
 
}
 
void ParticleSystem::computeForces()
{
    while(!mParticleContacts.empty()){
 
        ParticleContact contact = mParticleContacts.back();
        mParticleContacts.pop_back();
 
        int p1 = contact.mParticule1;
        int p2 = contact.mParticule2;
 
        //compute force
 
 
        float nX = (mParticles[p2].mPositionX - mParticles[p1].mPositionX)/contact.mDistance;
        float nY = (mParticles[p2].mPositionY - mParticles[p1].mPositionY)/contact.mDistance;
 
        float forceX = TIME_STEP * a * (mParticles[p1].mPressure + mParticles[p2].mPressure) * contact.mWeight * nX;
        float forceY = TIME_STEP * a * (mParticles[p1].mPressure + mParticles[p2].mPressure) * contact.mWeight * nY;
 
        mParticles[p1].mVelocityX += forceX;
        mParticles[p1].mVelocityY += forceY;
 
        mParticles[p2].mVelocityX -= forceX;
        mParticles[p2].mVelocityY -= forceY;
 
        mParticles[p1].mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p2].mVelocityX - mParticles[p1].mVelocityX);
        mParticles[p1].mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p2].mVelocityY - mParticles[p1].mVelocityY);
 
        mParticles[p2].mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p1].mVelocityX - mParticles[p2].mVelocityX);
        mParticles[p2].mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p1].mVelocityY - mParticles[p2].mVelocityY);
    }
}
 
void ParticleSystem::updatePosition(Map *pMap)
{
    // update position
    for (int i = 0; i < mCountParticles; ++i) {
 
        if(mParticles[i].mVelocityX>50)
            mParticles[i].mVelocityX = 50;
 
        if(mParticles[i].mVelocityX<-50)
            mParticles[i].mVelocityX = -50;
 
        if(mParticles[i].mVelocityY>50)
            mParticles[i].mVelocityY = 50;
 
        if(mParticles[i].mVelocityY<-50)
            mParticles[i].mVelocityY = -50;
 
        mParticles[i].updatePosition(pMap ,TIME_STEP);
 
        //reset
        mParticles[i].mDensity=0.0f;
        mParticles[i].mPressure=0.0f;
    }
}