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