BarnesHut_MPI/update_node.cc at master · Andrei-EPFL/BarnesHut_MPI · GitHub

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
#include "update_node.h"
#include <iostream>

MyNode* initialize_node(MyParticle particle, double bound_min_x, double bound_max_x, double bound_min_y, double bound_max_y)
{
    MyNode *node = new MyNode;
    node->elements = 1;
    node->totalmass = particle.mass;
    node->COM_x = particle.x;
    node->COM_y = particle.y;

    node->COM_vx = particle.vx;
    node->COM_vy = particle.vy;

    node->bound_min_x = bound_min_x;
    node->bound_max_x = bound_max_x;
    node->bound_min_y = bound_min_y;
    node->bound_max_y = bound_max_y;

    node->nw = NULL;
    node->ne = NULL;
    node->sw = NULL;
    node->se = NULL;
    return node;
}

void free_node(MyNode *node)
{
    if(node->nw){free_node(node->nw);}
    if(node->ne){free_node(node->ne);}
    if(node->sw){free_node(node->sw);}
    if(node->se){free_node(node->se);}
    if(node){delete(node);}
}

void update_child_node(MyNode *node, MyParticle particle, double bound_min_x, double bound_max_x, double bound_min_y, double bound_max_y)
{
    double x = particle.x;
    double y = particle.y;
    double x_center = (bound_min_x + bound_max_x)/2.;
    double y_center = (bound_min_y + bound_max_y)/2.;

    if(x<=x_center && y>y_center)
    {
        if(node->nw) {add_particle(node->nw, particle, bound_min_x, x_center, y_center, bound_max_y);}
        else {node->nw = initialize_node(particle, bound_min_x, x_center, y_center, bound_max_y);}
    }
    if(x<=x_center && y<=y_center)
    {
        if(node->sw) {add_particle(node->sw, particle, bound_min_x, x_center, bound_min_y, y_center );}
        else {node->sw = initialize_node(particle, bound_min_x, x_center, bound_min_y, y_center);}
    }
    if(x>x_center && y>y_center)
    {
        if(node->ne) {add_particle(node->ne, particle, x_center, bound_max_x, y_center, bound_max_y);}
        else {node->ne = initialize_node(particle, x_center, bound_max_x, y_center, bound_max_y);}
    }
    if(x>x_center && y<=y_center)
    {
        if(node->se) {add_particle(node->se, particle, x_center, bound_max_x, bound_min_y, y_center );}
        else {node->se = initialize_node(particle, x_center, bound_max_x, bound_min_y, y_center );}
    }
}

void update_current_node(MyNode *node, MyParticle particle)
{
    node->elements++;
    node->COM_x = (node->totalmass * node->COM_x + particle.mass * particle.x)/(node->totalmass+particle.mass);
    node->COM_y = (node->totalmass * node->COM_y + particle.mass * particle.y)/(node->totalmass+particle.mass);
    node->COM_vx = (node->totalmass * node->COM_vx + particle.mass * particle.vx)/(node->totalmass+particle.mass);
    node->COM_vy = (node->totalmass * node->COM_vy + particle.mass * particle.vy)/(node->totalmass+particle.mass);
    node->totalmass+=particle.mass;
}

void add_particle(MyNode *node, MyParticle newparticle, double bound_min_x, double bound_max_x, double bound_min_y, double bound_max_y)
{
    double x = newparticle.x;
    double y = newparticle.y;

    if(!node)
    {
        if(y>=bound_min_y && y<=bound_max_y && x>=bound_min_x && x<=bound_max_x)
        {
            std::cout << "This is the !node branch of the if clauses\n";
        }
        else
        {
            std::cout << "The particle is placed outside the boundaries(Node\n";
        }
    }
    else
    {
        //std::cout<<"This is the other branch of the if clause\n";
        if(y>=bound_min_y && y<=bound_max_y && x>=bound_min_x && x<=bound_max_x)
        {
            if(node->elements == 1)
            {
                MyParticle oldparticle;
                //std::cout<<"This node has one element\n";
                oldparticle.x = node->COM_x;
                oldparticle.y = node->COM_y;
                oldparticle.vx = node->COM_vx;
                oldparticle.vy = node->COM_vy;
                oldparticle.mass = node->totalmass;

                update_child_node(node, oldparticle, bound_min_x, bound_max_x, bound_min_y, bound_max_y);
            }
            update_current_node(node, newparticle);
            update_child_node(node, newparticle, bound_min_x, bound_max_x, bound_min_y, bound_max_y);
        }
        else
        {
            std::cout << "The particle is placed outside the boundaries\n";
        }
    }
}