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