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