Nbody问题 Barnes-Hut 实现

 

Barnes-Hut  算法

参考 http://arborjs.org/docs/barnes-hut

该算法对区域进行4分割。直到区域中只包含1个或者0个元素。

如下图

通过分割构造出如下树。

 

 

 

递归构造树的算法

 1 bool Tree::buildTree(NbodyNode *tree, complex<double> start, complex<double> end, vector<Plant> &plants) {
 2     if (tree == nullptr)
 3         return false;
 4     complex<double> check = end - start;
 5     if (check.real() <= 0 || check.imag() <= 0) {
 6         printf("check failed\n");
 7         return false;
 8     }
 9     if (plants.size() == 1) {
10         tree->body() = plants.front();
11         tree->isPlant() = true;
12         this->plants.push_back(&tree->body());
13         return true;
14     }
15     vector<Plant> wrapers[4];
16     int centerX = (start.real() + end.real()) / 2;
17     int centerY = (start.imag() + end.imag()) / 2;
18     complex<double> center = complex<double>(centerX, centerY);
19     complex<double> sub = complex<double>();
20     for (vector<Plant>::iterator i = plants.begin(); i != plants.end(); i++) {
21         sub = i->location() - center;
22         if (sub.real() <= 0 && sub.imag() <= 0) {
23             wrapers[0].push_back(*i);
24         } else if (sub.real() < 0 && sub.imag() > 0) {
25             wrapers[2].push_back(*i);
26         } else if (sub.real() > 0 && sub.imag() < 0) {
27             wrapers[1].push_back(*i);
28         } else if (sub.real() >= 0 && sub.imag() >= 0) {
29             wrapers[3].push_back(*i);
30         }
31     }
32     int width = tree->width() / 4;
33     tree->body() = Plant();
34     bool ret = true;
35     if (wrapers[0].size() > 0) {
36         tree->leftTop() = new NbodyNode(width);
37         ret = ret && buildTree(tree->leftTop(), start, center, wrapers[0]);
38         tree->body() = tree->body() + tree->leftTop()->body();
39     }
40     if (wrapers[1].size() > 0) {
41         tree->rightTop() = new NbodyNode(width);
42         ret = ret && buildTree(tree->rightTop(), complex<double>(start.real() + centerX, start.imag()),
43                                complex<double>(end.real(), centerY), wrapers[1]);
44         tree->body() = tree->body() + tree->rightTop()->body();
45     }
46     if (wrapers[2].size() > 0) {
47         tree->leftButtom() = new NbodyNode(width);
48         ret = ret && buildTree(tree->leftButtom(), complex<double>(start.real(), centerY),
49                                complex<double>(centerX, end.imag()), wrapers[2]);
50         tree->body() = tree->body() + tree->leftButtom()->body();
51     }
52     if (wrapers[3].size() > 0) {
53         tree->rightButtom() = new NbodyNode(width);
54         ret = ret && buildTree(tree->rightButtom(), center, end, wrapers[3]);
55         tree->body() = tree->body() + tree->rightButtom()->body();
56     }
57     return ret;
58 }

树中每一个非NULL节点保存该区域中星体的等效值。

若是星体，保存本身。若不是，保存该区域中的等效星体。

即

星体1 质量M1 位置（x1,y1）星体2 质量M2 位置（x2,y2）

等效星体 质量M = M1+M2 位置（x = (x1*M1+x*M2)/M, y = (y1*M1+y2*M2)/M）;

如下图

s 为该区域的宽度

d 为A星体到蓝色区域等效星体的距离

若 d/s < θ

则该区域可以被等效，否则计算该区域的子区域。

若区域本身是一个星体，则直接计算该星体对A的万有引力。不用计算 d/s

github 实现

转载于:https://www.cnblogs.com/like1/p/6933952.html