Discussion :

[dfshr8]

Bon je me lance en C++:
- car après l'objectif de mon stage est de modifier une librairie écrite en C++ donc autant m'y habituer (même si Python est très attractif)
- même si ce n'est pas un langage de script mais de compilation, en lançant le build and Run je peux voir immédiatement ce qui se passe en affichant au fur et à mesure du dévellopement les variables que j'utilise.

Je souhaite dévelloper l' algorithme FMM pour calculer le potentiel en un point généré par des charhes ponctuels. En gros je veux refaire ceci: www.umiacs.umd.edu/~ramani/fmm/files/fmmdemo.jar mais sans la partie graphique (pas besoin de l'animation je veux juste faire les calculs). Pour ceux qui voudraient une idée plus précise de l'algo, je m'appuie sur le document suivant (je garderas toutes les notations de ce doc): http://www.umiacs.umd.edu/labs/cvl/p...M_tutorial.pdf.
Je procéderait par essai erreur donc toutes vos remarques me seront utiles au cours du dévellopement!



Voici le début de mon code:

1) Je crée une distribution uniforme de N Particules dans le carré
2) On veut calculer le potentiel en la cible y_j par calcul directe. On considère tout le long de l'algo un potentiel en 1/r (on parle du Kernel) (ie) le potentiel en un point y arbitraire généré par toutes les sources x_i est $\phi(y_j)=\sum_{i=1}^n 1/\|y-x_i\|$ (peut on mettre du latex dans le forum ?)

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
#include<iostream>
#include <ctime>
#include<cstdlib>
#include<cmath>
using namespace std;
 
 
int main()
{
    srand(time(0));
 
    int N = 10; //Nombre particules
    int threshold = 5;
    int P=5; //ordre du dévelopement
 
    double x[N],y[N],m[N],p_dir[N];
 
 
// Création d'une distribution de 10 Particules au sein du cube [0,1]x[0,1]
    for( int i=0; i<N; i++ )
    {
        x[i] = rand()/(1.+RAND_MAX);
        y[i] = rand()/(1.+RAND_MAX);
        m[i] = 1.0/N;
        cout<<x[i]<<" "<<m[i]<<endl;
    }
 
//Calcul directe du potentiel en la cible y_j: $\phi(y_j)=\sup_{i=1}^n 1/\|y-x_i\|$
    double dx,dy,dz;//Représente la distance entre le centre et la particule concernée
    double r;
    for( int j=0; j<N; j++ )
    {
        for( int i=0; i<N; i++ )
        {
            if (i==j) continue;
            dx = x[i]-x[j];
            dy = y[i]-y[j];
 
            r = sqrt(dx*dx+dy*dy+dz*dz);
            p_dir[j] += m[i] / r;
        }
        cout<<"Le potentiel par calcul directe est p["<<j<<"]: "<<p_dir[j]<<endl;
    }

Remarques: On aura différent moyen (de qualité croissante: rapidité/précision) de calculer le potentiel en chacune des charges:
*calcul directe
* Middleman
* PreFmm/ Middleman amélioré
* SLFMM
*MLFMM
et pour chacunes de ces méthodes on utilise les mêmes types de calculs:
* calcul directe, Expansion multipolaire (dévelopement en série entière), opérateur de translation... Je pensais donc mettre ces dernières en fonctions, car réutilisés souvents.
Donc je pensais mettre chaque méthode dans une fonction.

[dfshr8]

Pour les méthodes de calculs suivantes (Middleman,...) c'est plus complexe. il faut que je partitionne/subdivise mon espace (mon carré) en niveaux (découpages en carrés cf. animation) et que je puisse détecter si une particule appartient à quel niveau/feuille. Pour ça je ne sais pas comment procéder.

PS: pour ceux qui comme moi ne savent pas comment voir l'animation, il suffit de faire java -jar fmm_demo2.jar.
PS: Dans l'animation, pour comprendre le partitionnement de l' espace, on choisit le nombre de niveaux et on a un découpage

[dfshr8]

J'ai modifié le code original pour le calcul directe (en m'inspirant du code Java de l'animation):

Classe point (car j'ai vu que c'est fait systèmatiquement, et ça me permet de tester l' égalite de deux points)

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
#include<iostream>
#include <ctime>
#include<cstdlib>
#include<cmath>
#include "Point.h"
using namespace std;
 
 
int main()
{
    srand(time(0));
 
    int N= 10; //Nombre particules sources et cibles
    //int M=12;
    //int threshold = 5;
    //int P=5; //ordre du dévelopement
 
    Point x[N]; //particules sources
    Point y[N];//particules cibles
    double m[N];// "Potential Weight": ici masses des particules
    double potential[N];//le " Potential Kernel" sera en 1/r
 
 
// Création d'une distribution de 10 Particules au sein du carré [0,1]x[0,1]
    for( int i=0; i<N; i++ )
    {
        x[i].setX(rand()/(1.+RAND_MAX));
        x[i].setY(rand()/(1.+RAND_MAX));
        y[i].setX(rand()/(1.+RAND_MAX));
        y[i].setY(rand()/(1.+RAND_MAX));
        m[i] = 1.0/N;
        cout<<x[i].getX()<<" "<<m[i]<<endl;
    }
 
//Calcul directe du potentiel en la cible y_j avec un " Potential Kernel" en 1/r: $\phi(y_j)=\sup_{i=1}^n 1/\|y-x_i\|$
    double dX,dY;//Représente la distance entre le centre et la particule concernée
    double r;
    for( int j=0; j<N; j++ )
    {
        for( int i=0; i<N; i++ )
        {
            if (x[i]==y[j]) continue;
            dX = x[i].getX()-y[j].getX();
            dY = x[i].getY()-y[j].getY();
            r = sqrt(dX*dX+dY*dY);
            potential[j] += m[i] / r;
        }
        cout<<"Le potentiel par calcul directe est p["<<j<<"]: "<<potential[j]<<endl;
    }