3 {
4 using double4 = Point4D;
5 using double2 = Point2D;
6 using int2 = std::pair<int, int>;
7
8
9
10#pragma omp parallel
11 {
12 int i, j, ch;
13
14 double4 lu[2];
15
16 double2 mom0;
17 double2 mom1;
18 double2 mom2;
19 double2 mom3;
20 double2 mom4;
21 double2 mom5;
22 double2 mom6;
23 double2 mom7;
24 double2 mom8;
25 double2 mom9;
26 double2 mom10;
27 double2 cen;
28 double2 dr;
29
30 int m[2];
31 int cm;
32 const int nnodes = 2 * (int)object.size() - 1;
33 const int nbodies = (int)object.size();
34#pragma omp for schedule(dynamic,5)
35 for (int k = nbodies; k < nnodes; ++k)
36 {
37
38
39
40
41
42
43
44
45
46
47
48 j = 0;
49 cm = 0;
50
51
52 while (cm == 0)
53 {
54 j = 2;
55 int srt = indexSort[(nnodes - 1) - k];
56 int2 chdPair = child[srt];
57
58 for (i = 0; i < 2; i++) {
59 int chd = i * chdPair.second + (1 - i) * chdPair.first;
60
61 ch = (chd >= nbodies) ? (chd - nbodies) : ((nnodes - 1) - indexSortT[chd]);
62
63 if ((chd >= nbodies) || (mass[nnodes - 1 - ch] >= 0))
64 j--;
65 }
66
67 if (j == 0)
68 {
69
71
72
73
74 for (i = 0; i < 2; i++)
75 {
76 const int chd = i * chdPair.second + (1 - i) * chdPair.first;
77 if (chd >= nbodies)
78 {
79 ch = chd - nbodies;
80 const int sortedBody = mortonCodesIdx[ch];
81
82 double4 xyAB = gabForLeaves[sortedBody];
83 lu[i] = double4{
84 ::fmin(xyAB[0], xyAB[2]), ::fmin(xyAB[1], xyAB[3]),
85 ::fmax(xyAB[0], xyAB[2]), ::fmax(xyAB[1], xyAB[3])
86 };
87 }
88 else
89 {
90 ch = indexSortT[chd];
91
92 lu[i] = lowerupper[chd];
93 }
94 }
95
96
97 lowerupper[srt] = double4{
98 ::fmin(lu[0][0], lu[1][0]), ::fmin(lu[0][1], lu[1][1]),
99 ::fmax(lu[0][2], lu[1][2]), ::fmax(lu[0][3], lu[1][3])
100 };
101
102 cen = center[srt] = double2{ 0.5 * (lowerupper[srt][0] + lowerupper[srt][2]),
103 0.5 * (lowerupper[srt][1] + lowerupper[srt][3]) };
104
105 const double2 zero = { 0.0, 0.0 };
106 double2 momh0 = zero;
107 double2 momh1 = zero;
108 double2 momh2 = zero;
109 double2 momh3 = zero;
110 double2 momh4 = zero;
111 double2 momh5 = zero;
112 double2 momh6 = zero;
113 double2 momh7 = zero;
114 double2 momh8 = zero;
115 double2 momh9 = zero;
116 double2 momh10 = zero;
117
118 for (i = 0; i < 2; i++)
119 {
120 const int chd = i * chdPair.second + (1 - i) * chdPair.first;
121 if (chd >= nbodies)
122 {
123 ch = chd - nbodies;
124 const int sortedBody = mortonCodesIdx[ch];
125
127 {
128 mom0 = double2{ gamma[sortedBody], 0.0 };
129
130mom1 = mom2 = mom3 = mom4 = mom5 = mom6 = mom7 = mom8 = mom9 = mom10 = double2{ 0.0, 0.0 };
131 double2 pos = object[sortedBody];
132 dr = pos - cen;
133 m[i] = 1;
134 }
135
137 {
138
139 }
140
141 }
142 else
143 {
144 const int srtT = indexSortT[chd];
145 ch = (nnodes - 1) - srtT;
146
148
149 mom0 = double2{ moms[nch + 0][0], (double)0 };
150 mom1 = moms[nch + 1];
151 mom2 = moms[nch + 2];
152 mom3 = moms[nch + 3];
153 mom4 = moms[nch + 4];
154 mom5 = moms[nch + 5];
155 mom6 = moms[nch + 6];
156 mom7 = moms[nch + 7];
157 mom8 = moms[nch + 8];
158 mom9 = moms[nch + 9];
159 mom10 = moms[nch + 10];
160 dr = center[chd] - cen;
161 m[i] = mass[srtT];
162 }
163
164
165 momh0 += mom0;
166 momh1 += mom1;
167 momh2 += mom2;
168 momh3 += mom3;
169 momh4 += mom4;
170 momh5 += mom5;
171 momh6 += mom6;
172 momh7 += mom7;
173 momh8 += mom8;
174 momh9 += mom9;
175 momh10 += mom10;
176
177 double2 z = dr;
178
179 momh1 +=
multz(mom0, z);
180 momh2 += 2 *
multz(mom1, z);
181 momh3 += 3 *
multz(mom2, z);
182 momh4 += 4 *
multz(mom3, z);
183 momh5 += 5 *
multz(mom4, z);
184 momh6 += 6 *
multz(mom5, z);
185 momh7 += 7 *
multz(mom6, z);
186 momh8 += 8 *
multz(mom7, z);
187 momh9 += 9 *
multz(mom8, z);
188 momh10 += 10 *
multz(mom9, z);
189
191
192 momh2 +=
multz(mom0, z);
193 momh3 += 3 *
multz(mom1, z);
194 momh4 += 6 *
multz(mom2, z);
195 momh5 += 10 *
multz(mom3, z);
196 momh6 += 15 *
multz(mom4, z);
197 momh7 += 21 *
multz(mom5, z);
198 momh8 += 28 *
multz(mom6, z);
199 momh9 += 36 *
multz(mom7, z);
200 momh10 += 45 *
multz(mom8, z);
201
203
204 momh3 +=
multz(mom0, z);
205 momh4 += 4 *
multz(mom1, z);
206 momh5 += 10 *
multz(mom2, z);
207 momh6 += 20 *
multz(mom3, z);
208 momh7 += 35 *
multz(mom4, z);
209 momh8 += 56 *
multz(mom5, z);
210 momh9 += 84 *
multz(mom6, z);
211 momh10 += 120 *
multz(mom7, z);
212
214
215 momh4 +=
multz(mom0, z);
216 momh5 += 5 *
multz(mom1, z);
217 momh6 += 15 *
multz(mom2, z);
218 momh7 += 35 *
multz(mom3, z);
219 momh8 += 70 *
multz(mom4, z);
220 momh9 += 126 *
multz(mom5, z);
221 momh10 += 210 *
multz(mom6, z);
222
224
225 momh5 +=
multz(mom0, z);
226 momh6 += 6 *
multz(mom1, z);
227 momh7 += 21 *
multz(mom2, z);
228 momh8 += 56 *
multz(mom3, z);
229 momh9 += 126 *
multz(mom4, z);
230 momh10 += 252 *
multz(mom5, z);
231
233
234 momh6 +=
multz(mom0, z);
235 momh7 += 7 *
multz(mom1, z);
236 momh8 += 28 *
multz(mom2, z);
237 momh9 += 84 *
multz(mom3, z);
238 momh10 += 210 *
multz(mom4, z);
239
241
242 momh7 +=
multz(mom0, z);
243 momh8 += 8 *
multz(mom1, z);
244 momh9 += 36 *
multz(mom2, z);
245 momh10 += 120 *
multz(mom3, z);
246
248
249 momh8 +=
multz(mom0, z);
250 momh9 += 9 *
multz(mom1, z);
251 momh10 += 45 *
multz(mom2, z);
252
254
255 momh9 +=
multz(mom0, z);
256 momh10 += 10 *
multz(mom1, z);
257
259
260 momh10 +=
multz(mom0, z);
261 }
262
263
264 momh0[1] = 0;
265 moms[kch + 0] = momh0;
266 moms[kch + 1] = momh1;
267 moms[kch + 2] = momh2;
268 moms[kch + 3] = momh3;
269 moms[kch + 4] = momh4;
270 moms[kch + 5] = momh5;
271 moms[kch + 6] = momh6;
272 moms[kch + 7] = momh7;
273 moms[kch + 8] = momh8;
274 moms[kch + 9] = momh9;
275 moms[kch + 10] = momh10;
276 cm = m[0] + m[1];
277 }
278
279#pragma omp flush
280
281 if (cm != 0)
282 {
283 mass[nnodes - 1 - k] = cm;
284
285
286 }
287 }
288 }
289 }
290 }
Point2D multz(const Point2D &a, const Point2D &b)
Умножение комплексных чисел
1.9.8