testUseNewApi.c 4.2 KB
Newer Older
1 2 3 4 5 6 7 8 9
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


#include "../Src/CScalfmmApi.h"



10 11 12 13 14 15 16
void compute_displacement_from_forces(double * fXYZ, double * charge,double* dXYZ, int nb_xyz){
    int idPart;
    for(idPart = 0 ; idPart < nb_xyz ; ++idPart){
        dXYZ[idPart*3+0] = fXYZ[idPart*3+0]*0.001/charge[idPart];
        dXYZ[idPart*3+1] = fXYZ[idPart*3+1]*0.001/charge[idPart];
        dXYZ[idPart*3+2] = fXYZ[idPart*3+2]*0.001/charge[idPart];
    }
17 18 19 20 21 22 23 24 25
}



int main(int argc, char ** av){

    scalfmm_kernel_type myChoice = chebyshev;

    //Octree configuration
26
    int TreeHeight = 7;
27
    double boxWidth = 2.0;
28 29 30 31
    double boxCenter[3] = {0.0,0.0,0.0};

    //Init our lib
    scalfmm_handle handle = scalfmm_init(TreeHeight,boxWidth,boxCenter,myChoice); //The tree is built
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
32
    scalfmm_algorithm_config(handle,periodic);
33
    //Creation of an array of particles
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
34
    int nb_of_parts = 2;
35 36
    int idxPart;
    double * positionsXYZ = malloc(sizeof(double)*3*nb_of_parts);
37 38
    memset(positionsXYZ,0,sizeof(double)*3*nb_of_parts);

39 40 41 42 43 44 45 46 47
    for(idxPart = 0; idxPart<nb_of_parts ; ++idxPart){
        positionsXYZ[idxPart*3+0] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[0];
        positionsXYZ[idxPart*3+1] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[1];
        positionsXYZ[idxPart*3+2] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[2];
    }

    //Creation of charge for each part
    double * array_of_charge = malloc(sizeof(double)*nb_of_parts);
    for(idxPart = 0; idxPart<nb_of_parts ; ++idxPart){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
48
        array_of_charge[idxPart] = 1.0;
49 50 51 52
    }

    //Inserting the array in the tree
    scalfmm_tree_insert_particles_xyz(handle,nb_of_parts,positionsXYZ);
53

54 55 56 57 58
    //Set the charge
    scalfmm_set_physical_values(handle,nb_of_parts,array_of_charge);


    //Computation Part
59
    int nb_iteration = 10;//atoi(av[1]);
60 61
    int curr_iteration = 0;

62
    //Array to store the output
63
    double * array_of_forces = malloc(sizeof(double)*3*nb_of_parts);
64 65 66
    memset(array_of_forces,0,sizeof(double)*3*nb_of_parts);
    double * array_of_pot = malloc(sizeof(double)*nb_of_parts);
    memset(array_of_pot,0,sizeof(double)*nb_of_parts);
67 68 69

    //Array to store the displacement
    double * array_of_displacement = malloc(sizeof(double)*3*nb_of_parts);
70
    memset(array_of_displacement,0,sizeof(double)*3*nb_of_parts);
71 72 73 74 75 76 77 78

    //Start of an iteration loop
    while(curr_iteration < nb_iteration){
        //Execute
        scalfmm_execute_fmm(handle);

        //Get the resulting forces
        scalfmm_get_forces_xyz(handle,nb_of_parts,array_of_forces);
79 80 81
        //Get the resulting potential
        scalfmm_get_potentials(handle,nb_of_parts,array_of_pot);

82 83


84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
        //User function to compute the movement of each part
        compute_displacement_from_forces(array_of_forces,array_of_charge,array_of_displacement,nb_of_parts);

        //get position in order to display
        scalfmm_get_positions_xyz(handle,nb_of_parts,positionsXYZ);

        //Display forces :
        {
            printf("Iteration n° %d\n",curr_iteration);
            for(idxPart = 0 ; idxPart< nb_of_parts ; ++idxPart){
                printf("Pos : %e, %e, %e, Fxyz %e %e %e, \n \t displ : %e, %e, %e \n",
                       positionsXYZ[idxPart*3+0],
                       positionsXYZ[idxPart*3+1],
                       positionsXYZ[idxPart*3+2],
                       array_of_forces[idxPart*3],
                       array_of_forces[idxPart*3+1],
                       array_of_forces[idxPart*3+2],
                       array_of_displacement[idxPart*3+0],
                       array_of_displacement[idxPart*3+1],
                       array_of_displacement[idxPart*3+2]);
            }
        }
106 107 108
        //Apply displacement computed
        scalfmm_add_to_positions_xyz(handle,nb_of_parts,array_of_displacement);

109 110 111
        //Update Consequently the tree
        scalfmm_update_tree(handle);
        //Continue the loop
112 113 114 115 116 117 118
        curr_iteration++;
    }

    //Free memory
    free(positionsXYZ);
    free(array_of_charge);
    free(array_of_forces);
119
    free(array_of_pot);
120 121
    free(array_of_displacement);

122 123 124 125
    //End of Computation, useer get the position after some iterations
    scalfmm_dealloc_handle(handle,NULL);


126 127
    return 0;
}