VM2D 1.14
Vortex methods for 2D flows simulation
Loading...
Searching...
No Matches
Functions
SummScript6.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void Summ6 ()
 

Function Documentation

◆ Summ6()

void Summ6 ( )
inline

Definition at line 2 of file SummScript6.h.

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