首页 > > 程序设计 > C/C++ >

Nbody问题 Barnes-Hut 实现

2018-06-17 22:29:10来源：未知阅读 ()

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 }