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