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