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 實現