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