testUseNewApi.c 4.76 KB
Newer Older
1 2
// See LICENCE file at project root

3 4 5 6 7 8 9 10 11
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


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



12 13 14 15 16 17 18
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];
    }
19 20 21 22 23 24 25 26 27
}



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

    scalfmm_kernel_type myChoice = chebyshev;

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

    //Init our lib
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
33
    scalfmm_handle handle = scalfmm_init(/* TreeHeight,boxWidth,boxCenter, */myChoice,multi_thread); //The tree is built
34 35 36
    struct User_Scalfmm_Cell_Descriptor user_descr;
    user_descr.user_init_cell = NULL;
    user_descr.user_free_cell = NULL;
37 38 39 40 41 42 43 44
    struct User_Scalfmm_Leaf_Descriptor user_descr_leaf;
    user_descr_leaf.user_init_leaf = NULL;
    user_descr_leaf.user_free_leaf = NULL;
    user_descr_leaf.user_get_size = NULL;
    user_descr_leaf.user_copy_leaf = NULL;
    user_descr_leaf.user_restore_leaf = NULL;

    scalfmm_build_tree(handle,TreeHeight,boxWidth,boxCenter,user_descr,user_descr_leaf);
45
    scalfmm_algorithm_config(handle,periodic);
46
    //Creation of an array of particles
47
    int nb_of_parts = 2;
48
    FSize idxPart;
49
    double * positionsXYZ = malloc(sizeof(double)*3*nb_of_parts);
50 51
    memset(positionsXYZ,0,sizeof(double)*3*nb_of_parts);

52 53 54 55 56 57 58 59 60
    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){
61
        array_of_charge[idxPart] = 1.0;
62 63 64
    }

    //Inserting the array in the tree
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
65
    scalfmm_tree_insert_particles_xyz(handle,nb_of_parts,positionsXYZ,BOTH);
66

67
    //Set the charge
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
68
    scalfmm_set_physical_values(handle,nb_of_parts,array_of_charge,BOTH);
69 70 71


    //Computation Part
72
    int nb_iteration = 3;//atoi(av[1]);
73 74
    int curr_iteration = 0;

75
    //Array to store the output
76
    double * array_of_forces = malloc(sizeof(double)*3*nb_of_parts);
77 78 79
    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);
80 81 82

    //Array to store the displacement
    double * array_of_displacement = malloc(sizeof(double)*3*nb_of_parts);
83
    memset(array_of_displacement,0,sizeof(double)*3*nb_of_parts);
84 85 86 87 88 89 90

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

        //Get the resulting forces
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
91
        scalfmm_get_forces_xyz(handle,nb_of_parts,array_of_forces,BOTH);
92
        //Get the resulting potential
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
93
        scalfmm_get_potentials(handle,nb_of_parts,array_of_pot,BOTH);
94

95 96


97 98 99 100
        //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
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
101
        scalfmm_get_positions_xyz(handle,nb_of_parts,positionsXYZ,BOTH);
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118

        //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]);
            }
        }
119
        //Apply displacement computed
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
120
        scalfmm_add_to_positions_xyz(handle,nb_of_parts,array_of_displacement,BOTH);
121

122 123 124
        //Update Consequently the tree
        scalfmm_update_tree(handle);
        //Continue the loop
125 126 127 128 129 130 131
        curr_iteration++;
    }

    //Free memory
    free(positionsXYZ);
    free(array_of_charge);
    free(array_of_forces);
132
    free(array_of_pot);
133 134
    free(array_of_displacement);

135 136 137 138
    //End of Computation, useer get the position after some iterations
    scalfmm_dealloc_handle(handle,NULL);


139 140
    return 0;
}