|
2 | 2 | // it under the terms of the GNU General Public License as published by |
3 | 3 | // the Free Software Foundation, either version 3 of the License, or |
4 | 4 | // (at your option) any later version. |
5 | | - |
6 | 5 | // This program is distributed in the hope that it will be useful, |
7 | 6 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
8 | 7 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
9 | 8 | // GNU General Public License for more details. |
10 | | - |
11 | 9 | // You should have received a copy of the GNU General Public License |
12 | 10 | // along with this program. If not, see <http://www.gnu.org/licenses/>. |
13 | 11 |
|
14 | 12 | #pragma once |
15 | 13 | #include "../struct/particle.h" |
16 | | -#include <vector> |
17 | | -#include "octree.h" |
18 | 14 | #include "floatdef.h" |
| 15 | +#include "octree.h" |
| 16 | +#include <vector> |
19 | 17 |
|
20 | | -inline void Step(std::vector<Particle>& p, real dt) { |
21 | | - if (p.empty()) return; |
22 | | - |
23 | | - real theta = 0.5; |
24 | | - real half = dt * real(0.5); |
25 | | - |
26 | | - auto buildTree = [&](Octree*& root) { |
27 | | - // Compute bounding box |
28 | | - real minx=+1e30, miny=+1e30, minz=+1e30; |
29 | | - real maxx=-1e30, maxy=-1e30, maxz=-1e30; |
30 | | - |
31 | | - for (auto& a : p) { |
32 | | - minx = std::min(minx, a.x); |
33 | | - miny = std::min(miny, a.y); |
34 | | - minz = std::min(minz, a.z); |
35 | | - maxx = std::max(maxx, a.x); |
36 | | - maxy = std::max(maxy, a.y); |
37 | | - maxz = std::max(maxz, a.z); |
38 | | - } |
39 | | - |
40 | | - real cx = (minx + maxx) * 0.5; |
41 | | - real cy = (miny + maxy) * 0.5; |
42 | | - real cz = (minz + maxz) * 0.5; |
43 | | - real dx = maxx - minx; |
44 | | - real dy = maxy - miny; |
45 | | - real dz = maxz - minz; |
| 18 | +inline void Step(std::vector<Particle> &p, real dt) { |
| 19 | + if (p.empty()) |
| 20 | + return; |
| 21 | + |
| 22 | + real theta = 0.5; |
| 23 | + real half = dt * real(0.5); |
| 24 | + |
| 25 | + auto buildTree = [&](Octree *&root) { |
| 26 | + // Compute bounding box |
| 27 | + real minx = +1e30, miny = +1e30, minz = +1e30; |
| 28 | + real maxx = -1e30, maxy = -1e30, maxz = -1e30; |
| 29 | + |
| 30 | + for (auto &a : p) { |
| 31 | + minx = std::min(minx, a.x); |
| 32 | + miny = std::min(miny, a.y); |
| 33 | + minz = std::min(minz, a.z); |
| 34 | + maxx = std::max(maxx, a.x); |
| 35 | + maxy = std::max(maxy, a.y); |
| 36 | + maxz = std::max(maxz, a.z); |
| 37 | + } |
46 | 38 |
|
47 | | - real size = std::max(dx, std::max(dy, dz)) * real(0.5); |
| 39 | + real cx = (minx + maxx) * 0.5; |
| 40 | + real cy = (miny + maxy) * 0.5; |
| 41 | + real cz = (minz + maxz) * 0.5; |
| 42 | + real dx = maxx - minx; |
| 43 | + real dy = maxy - miny; |
| 44 | + real dz = maxz - minz; |
48 | 45 |
|
49 | | - if (size <= 0) size = 1; // safety |
| 46 | + real size = std::max(dx, std::max(dy, dz)) * real(0.5); |
50 | 47 |
|
51 | | - root = new Octree(cx, cy, cz, size); |
| 48 | + if (size <= 0) |
| 49 | + size = 1; // safety |
52 | 50 |
|
53 | | - for (auto& a : p) |
54 | | - root->insert(&a); |
| 51 | + root = new Octree(cx, cy, cz, size); |
55 | 52 |
|
56 | | - root->computeMass(); |
57 | | - }; |
| 53 | + for (auto &a : p) |
| 54 | + root->insert(&a); |
58 | 55 |
|
59 | | - // ========================= |
60 | | - // First Kick (dt/2) |
61 | | - // ========================= |
62 | | - { |
63 | | - Octree* root = nullptr; |
64 | | - buildTree(root); |
| 56 | + root->computeMass(); |
| 57 | + }; |
65 | 58 |
|
66 | | - #pragma omp parallel for schedule(static) |
67 | | - for (int i = 0; i < (int)p.size(); i++) { |
68 | | - real ax = 0, ay = 0, az = 0; |
69 | | - bhAccel(root, p[i], theta, ax, ay, az); |
| 59 | + // ========================= |
| 60 | + // First Kick (dt/2) |
| 61 | + // ========================= |
| 62 | + { |
| 63 | + Octree *root = nullptr; |
| 64 | + buildTree(root); |
70 | 65 |
|
71 | | - p[i].vx += ax * half; |
72 | | - p[i].vy += ay * half; |
73 | | - p[i].vz += az * half; |
74 | | - } |
| 66 | +#pragma omp parallel for schedule(static) |
| 67 | + for (int i = 0; i < (int)p.size(); i++) { |
| 68 | + real ax = 0, ay = 0, az = 0; |
| 69 | + bhAccel(root, p[i], theta, ax, ay, az); |
75 | 70 |
|
76 | | - delete root; |
| 71 | + p[i].vx += ax * half; |
| 72 | + p[i].vy += ay * half; |
| 73 | + p[i].vz += az * half; |
77 | 74 | } |
78 | 75 |
|
79 | | - // ========================= |
80 | | - // Drift (dt) |
81 | | - // ========================= |
82 | | - #pragma omp parallel for schedule(static) |
| 76 | + delete root; |
| 77 | + } |
| 78 | + |
| 79 | +// ========================= |
| 80 | +// Drift (dt) |
| 81 | +// ========================= |
| 82 | +#pragma omp parallel for schedule(static) |
| 83 | + for (int i = 0; i < (int)p.size(); i++) { |
| 84 | + p[i].x += p[i].vx * dt; |
| 85 | + p[i].y += p[i].vy * dt; |
| 86 | + p[i].z += p[i].vz * dt; |
| 87 | + } |
| 88 | + |
| 89 | + // ========================= |
| 90 | + // Second Kick (dt/2) |
| 91 | + // ========================= |
| 92 | + { |
| 93 | + Octree *root = nullptr; |
| 94 | + buildTree(root); |
| 95 | + |
| 96 | +#pragma omp parallel for schedule(static) |
83 | 97 | for (int i = 0; i < (int)p.size(); i++) { |
84 | | - p[i].x += p[i].vx * dt; |
85 | | - p[i].y += p[i].vy * dt; |
86 | | - p[i].z += p[i].vz * dt; |
87 | | - } |
| 98 | + real ax = 0, ay = 0, az = 0; |
| 99 | + bhAccel(root, p[i], theta, ax, ay, az); |
88 | 100 |
|
89 | | - // ========================= |
90 | | - // Second Kick (dt/2) |
91 | | - // ========================= |
92 | | - { |
93 | | - Octree* root = nullptr; |
94 | | - buildTree(root); |
95 | | - |
96 | | - #pragma omp parallel for schedule(static) |
97 | | - for (int i = 0; i < (int)p.size(); i++) { |
98 | | - real ax = 0, ay = 0, az = 0; |
99 | | - bhAccel(root, p[i], theta, ax, ay, az); |
100 | | - |
101 | | - p[i].vx += ax * half; |
102 | | - p[i].vy += ay * half; |
103 | | - p[i].vz += az * half; |
104 | | - } |
105 | | - |
106 | | - delete root; |
| 101 | + p[i].vx += ax * half; |
| 102 | + p[i].vy += ay * half; |
| 103 | + p[i].vz += az * half; |
107 | 104 | } |
| 105 | + |
| 106 | + delete root; |
| 107 | + } |
108 | 108 | } |
0 commit comments