VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
SummScript16.h
Go to the documentation of this file.
1
2 inline void Summ16()
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]; //bounding boxes двух детей
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 mom11;
28 double2 mom12;
29 double2 mom13;
30 double2 mom14;
31 double2 mom15;
32 double2 cen; //центр текущего родительского узла
33 double2 dr; //вектор от центра родителя к центру ребенка
34
35 int m[2]; //для листа это 1, для внутреннего узла это число листьев в нем
36 int cm; //масса текущего узла = сумма масс двух детей
37 const int nnodes = 2 * (int)object.size() - 1;
38 const int nbodies = (int)object.size();
39#pragma omp for schedule(dynamic,5)
40 for (int k = nbodies; k < nnodes; ++k)
41 {
42 //MortonTree:
43 // 0 1 2 ... (nb-2) x (nb+0) (nb+1) (nb+2) ... (nb+(nb-1))
44 // ---------------- -----------------------------------
45 // cells bodies
46
47 //Martin's tree:
48 // 0 1 2 ... (nb-1) x x x x (nn-(nb-1)) ... (nn-2) (nn-1)
49 // ---------------- ----------------------------
50 // bodies sorted and reversed cells
51 //printf("k = %d\n", k);
52
53 j = 0;
54 cm = 0;
55
56 // iterate over all cells assigned to thread
57 while (cm == 0)
58 {
59 j = 2;
60 int srt = indexSort[(nnodes - 1) - k]; //проход снизу вверх - т.е. в обратном порядке
61 int2 chdPair = child[srt];
62
63 for (i = 0; i < 2; i++) {
64 int chd = i * chdPair.second + (1 - i) * chdPair.first; // i==0 => .x; i==1 => .y
65
66 ch = (chd >= nbodies) ? (chd - nbodies) : ((nnodes - 1) - indexSortT[chd]);
67
68 if ((chd >= nbodies) || (mass[nnodes - 1 - ch] >= 0))
69 j--;
70 }
71
72 if (j == 0)
73 {
74 // all children are ready
75 const int kch = ((nnodes - 1) - k) * orderAlignment; //позиция, куда будут записаны мм для текущего внутреннего узла в массиве moms
76 //moms хранится не по индексу srt, а по отсортированному порядку узлов
77
78 //Считаем bounding box текущего узла из двух детей
79 for (i = 0; i < 2; i++)
80 {
81 const int chd = i * chdPair.second + (1 - i) * chdPair.first; // индекс i-го ребенка
82 if (chd >= nbodies)//если ребенок - это лист
83 {
84 ch = chd - nbodies; //номер частицы в отсортированном массиве
85 const int sortedBody = mortonCodesIdx[ch];
86
87 double4 xyAB = gabForLeaves[sortedBody];
88 lu[i] = double4{
89 ::fmin(xyAB[0], xyAB[2]), ::fmin(xyAB[1], xyAB[3]),
90 ::fmax(xyAB[0], xyAB[2]), ::fmax(xyAB[1], xyAB[3])
91 };
92 }
93 else
94 {
95 ch = indexSortT[chd];
96 // если внутренний узел, его bounding box уже должен быть посчитан, тк идем снизу вверх
97 lu[i] = lowerupper[chd];
98 }
99 }//for i
100
101 // Объединяем bounding box двух детей:
102 lowerupper[srt] = double4{
103 ::fmin(lu[0][0], lu[1][0]), ::fmin(lu[0][1], lu[1][1]),
104 ::fmax(lu[0][2], lu[1][2]), ::fmax(lu[0][3], lu[1][3])
105 };
106 // Центр текущего узла
107 cen = center[srt] = double2{ 0.5 * (lowerupper[srt][0] + lowerupper[srt][2]),
108 0.5 * (lowerupper[srt][1] + lowerupper[srt][3]) };
109
110 const double2 zero = { 0.0, 0.0 };
111 double2 momh0 = zero;
112 double2 momh1 = zero;
113 double2 momh2 = zero;
114 double2 momh3 = zero;
115 double2 momh4 = zero;
116 double2 momh5 = zero;
117 double2 momh6 = zero;
118 double2 momh7 = zero;
119 double2 momh8 = zero;
120 double2 momh9 = zero;
121 double2 momh10 = zero;
122 double2 momh11 = zero;
123 double2 momh12 = zero;
124 double2 momh13 = zero;
125 double2 momh14 = zero;
126 double2 momh15 = zero;
127 //переносим мультипольные моменты в центр родительской ячейки
128 for (i = 0; i < 2; i++)
129 {
130 const int chd = i * chdPair.second + (1 - i) * chdPair.first;
131 if (chd >= nbodies) //если ребенок - это лист
132 {
133 ch = chd - nbodies;
134 const int sortedBody = mortonCodesIdx[ch];
135
136 if (objectType == object_T::point4)
137 {
138 mom0 = double2{ gamma[sortedBody], 0.0 };
139 //для вихря все остальные мм нулевые
140mom1 = mom2 = mom3 = mom4 = mom5 = mom6 = mom7 = mom8 = mom9 = mom10 = mom11 = mom12 = mom13 = mom14 = mom15 = double2{ 0.0, 0.0 };
141 double2 pos = object[sortedBody];
142 dr = pos - cen;
143 m[i] = 1;
144 } //objectType==point4
145
146 if (objectType == object_T::panel)
147 {
148 //...
149 } //objectType == panel
150
151 }
152 else // если ребенок - внутренний узел
153 {
154 const int srtT = indexSortT[chd];
155 ch = (nnodes - 1) - srtT;
156
157 const int nch = srtT * orderAlignment;
158
159 mom0 = double2{ moms[nch + 0][0], (double)0 };
160 mom1 = moms[nch + 1];
161 mom2 = moms[nch + 2];
162 mom3 = moms[nch + 3];
163 mom4 = moms[nch + 4];
164 mom5 = moms[nch + 5];
165 mom6 = moms[nch + 6];
166 mom7 = moms[nch + 7];
167 mom8 = moms[nch + 8];
168 mom9 = moms[nch + 9];
169 mom10 = moms[nch + 10];
170 mom11 = moms[nch + 11];
171 mom12 = moms[nch + 12];
172 mom13 = moms[nch + 13];
173 mom14 = moms[nch + 14];
174 mom15 = moms[nch + 15];
175 dr = center[chd] - cen;
176 m[i] = mass[srtT];
177 }
178 // add child's contribution
179
180 momh0 += mom0;
181 momh1 += mom1;
182 momh2 += mom2;
183 momh3 += mom3;
184 momh4 += mom4;
185 momh5 += mom5;
186 momh6 += mom6;
187 momh7 += mom7;
188 momh8 += mom8;
189 momh9 += mom9;
190 momh10 += mom10;
191 momh11 += mom11;
192 momh12 += mom12;
193 momh13 += mom13;
194 momh14 += mom14;
195 momh15 += mom15;
196
197 double2 z = dr;
198
199 momh1 += multz(mom0, z);
200 momh2 += 2 * multz(mom1, z);
201 momh3 += 3 * multz(mom2, z);
202 momh4 += 4 * multz(mom3, z);
203 momh5 += 5 * multz(mom4, z);
204 momh6 += 6 * multz(mom5, z);
205 momh7 += 7 * multz(mom6, z);
206 momh8 += 8 * multz(mom7, z);
207 momh9 += 9 * multz(mom8, z);
208 momh10 += 10 * multz(mom9, z);
209 momh11 += 11 * multz(mom10, z);
210 momh12 += 12 * multz(mom11, z);
211 momh13 += 13 * multz(mom12, z);
212 momh14 += 14 * multz(mom13, z);
213 momh15 += 15 * multz(mom14, z);
214
215 z = multz(z, dr);
216
217 momh2 += multz(mom0, z);
218 momh3 += 3 * multz(mom1, z);
219 momh4 += 6 * multz(mom2, z);
220 momh5 += 10 * multz(mom3, z);
221 momh6 += 15 * multz(mom4, z);
222 momh7 += 21 * multz(mom5, z);
223 momh8 += 28 * multz(mom6, z);
224 momh9 += 36 * multz(mom7, z);
225 momh10 += 45 * multz(mom8, z);
226 momh11 += 55 * multz(mom9, z);
227 momh12 += 66 * multz(mom10, z);
228 momh13 += 78 * multz(mom11, z);
229 momh14 += 91 * multz(mom12, z);
230 momh15 += 105 * multz(mom13, z);
231
232 z = multz(z, dr);
233
234 momh3 += multz(mom0, z);
235 momh4 += 4 * multz(mom1, z);
236 momh5 += 10 * multz(mom2, z);
237 momh6 += 20 * multz(mom3, z);
238 momh7 += 35 * multz(mom4, z);
239 momh8 += 56 * multz(mom5, z);
240 momh9 += 84 * multz(mom6, z);
241 momh10 += 120 * multz(mom7, z);
242 momh11 += 165 * multz(mom8, z);
243 momh12 += 220 * multz(mom9, z);
244 momh13 += 286 * multz(mom10, z);
245 momh14 += 364 * multz(mom11, z);
246 momh15 += 455 * multz(mom12, z);
247
248 z = multz(z, dr);
249
250 momh4 += multz(mom0, z);
251 momh5 += 5 * multz(mom1, z);
252 momh6 += 15 * multz(mom2, z);
253 momh7 += 35 * multz(mom3, z);
254 momh8 += 70 * multz(mom4, z);
255 momh9 += 126 * multz(mom5, z);
256 momh10 += 210 * multz(mom6, z);
257 momh11 += 330 * multz(mom7, z);
258 momh12 += 495 * multz(mom8, z);
259 momh13 += 715 * multz(mom9, z);
260 momh14 += 1001 * multz(mom10, z);
261 momh15 += 1365 * multz(mom11, z);
262
263 z = multz(z, dr);
264
265 momh5 += multz(mom0, z);
266 momh6 += 6 * multz(mom1, z);
267 momh7 += 21 * multz(mom2, z);
268 momh8 += 56 * multz(mom3, z);
269 momh9 += 126 * multz(mom4, z);
270 momh10 += 252 * multz(mom5, z);
271 momh11 += 462 * multz(mom6, z);
272 momh12 += 792 * multz(mom7, z);
273 momh13 += 1287 * multz(mom8, z);
274 momh14 += 2002 * multz(mom9, z);
275 momh15 += 3003 * multz(mom10, z);
276
277 z = multz(z, dr);
278
279 momh6 += multz(mom0, z);
280 momh7 += 7 * multz(mom1, z);
281 momh8 += 28 * multz(mom2, z);
282 momh9 += 84 * multz(mom3, z);
283 momh10 += 210 * multz(mom4, z);
284 momh11 += 462 * multz(mom5, z);
285 momh12 += 924 * multz(mom6, z);
286 momh13 += 1716 * multz(mom7, z);
287 momh14 += 3003 * multz(mom8, z);
288 momh15 += 5005 * multz(mom9, z);
289
290 z = multz(z, dr);
291
292 momh7 += multz(mom0, z);
293 momh8 += 8 * multz(mom1, z);
294 momh9 += 36 * multz(mom2, z);
295 momh10 += 120 * multz(mom3, z);
296 momh11 += 330 * multz(mom4, z);
297 momh12 += 792 * multz(mom5, z);
298 momh13 += 1716 * multz(mom6, z);
299 momh14 += 3432 * multz(mom7, z);
300 momh15 += 6435 * multz(mom8, z);
301
302 z = multz(z, dr);
303
304 momh8 += multz(mom0, z);
305 momh9 += 9 * multz(mom1, z);
306 momh10 += 45 * multz(mom2, z);
307 momh11 += 165 * multz(mom3, z);
308 momh12 += 495 * multz(mom4, z);
309 momh13 += 1287 * multz(mom5, z);
310 momh14 += 3003 * multz(mom6, z);
311 momh15 += 6435 * multz(mom7, z);
312
313 z = multz(z, dr);
314
315 momh9 += multz(mom0, z);
316 momh10 += 10 * multz(mom1, z);
317 momh11 += 55 * multz(mom2, z);
318 momh12 += 220 * multz(mom3, z);
319 momh13 += 715 * multz(mom4, z);
320 momh14 += 2002 * multz(mom5, z);
321 momh15 += 5005 * multz(mom6, z);
322
323 z = multz(z, dr);
324
325 momh10 += multz(mom0, z);
326 momh11 += 11 * multz(mom1, z);
327 momh12 += 66 * multz(mom2, z);
328 momh13 += 286 * multz(mom3, z);
329 momh14 += 1001 * multz(mom4, z);
330 momh15 += 3003 * multz(mom5, z);
331
332 z = multz(z, dr);
333
334 momh11 += multz(mom0, z);
335 momh12 += 12 * multz(mom1, z);
336 momh13 += 78 * multz(mom2, z);
337 momh14 += 364 * multz(mom3, z);
338 momh15 += 1365 * multz(mom4, z);
339
340 z = multz(z, dr);
341
342 momh12 += multz(mom0, z);
343 momh13 += 13 * multz(mom1, z);
344 momh14 += 91 * multz(mom2, z);
345 momh15 += 455 * multz(mom3, z);
346
347 z = multz(z, dr);
348
349 momh13 += multz(mom0, z);
350 momh14 += 14 * multz(mom1, z);
351 momh15 += 105 * multz(mom2, z);
352
353 z = multz(z, dr);
354
355 momh14 += multz(mom0, z);
356 momh15 += 15 * multz(mom1, z);
357
358 z = multz(z, dr);
359
360 momh15 += multz(mom0, z);
361 }
362
363 // Сохраняем итоговые моменты текущего внутреннего узла в общий массив
364 momh0[1] = 0;
365 moms[kch + 0] = momh0;
366 moms[kch + 1] = momh1;
367 moms[kch + 2] = momh2;
368 moms[kch + 3] = momh3;
369 moms[kch + 4] = momh4;
370 moms[kch + 5] = momh5;
371 moms[kch + 6] = momh6;
372 moms[kch + 7] = momh7;
373 moms[kch + 8] = momh8;
374 moms[kch + 9] = momh9;
375 moms[kch + 10] = momh10;
376 moms[kch + 11] = momh11;
377 moms[kch + 12] = momh12;
378 moms[kch + 13] = momh13;
379 moms[kch + 14] = momh14;
380 moms[kch + 15] = momh15;
381 cm = m[0] + m[1];
382 }
383
384#pragma omp flush
385
386 if (cm != 0)
387 {
388 mass[nnodes - 1 - k] = cm;// Записываем массу текущего узла
389 //k += inc;
390 //flag = 0;
391 }
392 }//while flag==0
393 }//for k
394 }
395 }//Summ12()
orderAlignment
#define orderAlignment
Definition Gpudefs.h:94
Summ16
void Summ16()
Definition SummScript16.h:2