VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
SummScript11.h
Go to the documentation of this file.
1
2 inline void Summ11()
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 cen; //центр текущего родительского узла
28 double2 dr; //вектор от центра родителя к центру ребенка
29
30 int m[2]; //для листа это 1, для внутреннего узла это число листьев в нем
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 //MortonTree:
38 // 0 1 2 ... (nb-2) x (nb+0) (nb+1) (nb+2) ... (nb+(nb-1))
39 // ---------------- -----------------------------------
40 // cells bodies
41
42 //Martin's tree:
43 // 0 1 2 ... (nb-1) x x x x (nn-(nb-1)) ... (nn-2) (nn-1)
44 // ---------------- ----------------------------
45 // bodies sorted and reversed cells
46 //printf("k = %d\n", k);
47
48 j = 0;
49 cm = 0;
50
51 // iterate over all cells assigned to thread
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; // i==0 => .x; i==1 => .y
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 // all children are ready
70 const int kch = ((nnodes - 1) - k) * orderAlignment; //позиция, куда будут записаны мм для текущего внутреннего узла в массиве moms
71 //moms хранится не по индексу srt, а по отсортированному порядку узлов
72
73 //Считаем bounding box текущего узла из двух детей
74 for (i = 0; i < 2; i++)
75 {
76 const int chd = i * chdPair.second + (1 - i) * chdPair.first; // индекс i-го ребенка
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 // если внутренний узел, его bounding box уже должен быть посчитан, тк идем снизу вверх
92 lu[i] = lowerupper[chd];
93 }
94 }//for i
95
96 // Объединяем bounding box двух детей:
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
126 if (objectType == object_T::point4)
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 } //objectType==point4
135
136 if (objectType == object_T::panel)
137 {
138 //...
139 } //objectType == panel
140
141 }
142 else // если ребенок - внутренний узел
143 {
144 const int srtT = indexSortT[chd];
145 ch = (nnodes - 1) - srtT;
146
147 const int nch = srtT * orderAlignment;
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 // add child's contribution
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
190 z = multz(z, dr);
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
202 z = multz(z, dr);
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
213 z = multz(z, dr);
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
223 z = multz(z, dr);
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
232 z = multz(z, dr);
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
240 z = multz(z, dr);
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
247 z = multz(z, dr);
248
249 momh8 += multz(mom0, z);
250 momh9 += 9 * multz(mom1, z);
251 momh10 += 45 * multz(mom2, z);
252
253 z = multz(z, dr);
254
255 momh9 += multz(mom0, z);
256 momh10 += 10 * multz(mom1, z);
257
258 z = multz(z, dr);
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 //k += inc;
285 //flag = 0;
286 }
287 }//while flag==0
288 }//for k
289 }
290 }//Summ12()
orderAlignment
#define orderAlignment
Definition Gpudefs.h:94
Summ11
void Summ11()
Definition SummScript11.h:2