310{
312
313 double preTime = -omp_get_wtime();
314
315 const size_t n = vtx.size();
316
318 auto xx = std::minmax_element(vtx.begin(), vtx.end(), [](
const Vortex2D& P1,
const Vortex2D& P2) {return P1.r()[0] < P2.r()[0]; });
319 auto yy = std::minmax_element(vtx.begin(), vtx.end(), [](
const Vortex2D& P1,
const Vortex2D& P2) {return P1.r()[1] < P2.r()[1]; });
320
321 double scale = std::max(xx.second->r()[0] - xx.first->r()[0], yy.second->r()[1] - yy.first->r()[1]);
322
323
324 std::vector<unsigned int> s(256 * omp_get_max_threads());
325
326 std::vector<TParticleCode> mcdata(vtx.size());
327 std::vector<TParticleCode> mcdata_temp(vtx.size());
328
329
330
331
332 const int nSdvig = 5;
333 double tStart[nSdvig], tFinish[nSdvig];
334 double tm[nSdvig][5];
335
336 preTime += omp_get_wtime();
337
338 std::pair<double, size_t> zeroPair = std::make_pair<double, size_t>(0.0, 0);
339 std::pair<double, size_t> minusOnePair = std::make_pair<double, size_t>(-1.0, 0);
340
341 std::vector<std::pair<double, size_t>>
dist(2 * k, { -1.0, 0 });
342 std::vector<size_t> loc(2 * k);
343 std::vector<size_t> counter(2 * k, 0);
344 std::vector<size_t> offset(2 * k, 0);
345 std::vector<size_t> counterScan(2 * k, 0);
346 std::vector<std::pair<double, size_t>> dstKeys1(2 * k, zeroPair);
347
348 std::vector<std::pair<double, size_t>> updateNN(k, zeroPair);
349 std::vector<std::pair<double, size_t>> dstKeys(k, zeroPair);
350
351
352 for (size_t sdvig = 0; sdvig < nSdvig ; ++sdvig)
353 {
354 tStart[sdvig] = omp_get_wtime();
355
356 const Point2D& lowLeft = { xx.first->r()[0], yy.first->r()[1] };
357 double* time = tm[sdvig];
358
359
360
361 time[0] = omp_get_wtime();
362#pragma omp parallel for
363 for (int i = 0; i < vtx.size(); ++i)
364 {
365 Point2D sh = (vtx[i].r() - lowLeft) * (0.75 / scale) + sdvig *
Point2D{ 0.05, 0.05 };
367 mcdata[i].originNumber = i;
368 }
369
370 time[1] = omp_get_wtime();
371
372
373 RSort_Parallel(mcdata.data(), mcdata_temp.data(), (
int)mcdata.size(), s.data());
374
375 time[2] = omp_get_wtime();
376
377 for (size_t idx = 0; idx < 2 * k; ++idx)
378 {
379 dist[idx] = minusOnePair;
380 dstKeys1[idx] = zeroPair;
381 loc[idx] = counter[idx] = offset[idx] = counterScan[idx] = 0;
382 }
383
384 for (size_t idx = 0; idx < k; ++idx)
385 updateNN[idx] = dstKeys[idx] = zeroPair;
386
387 time[3] = omp_get_wtime();
388
389#pragma omp parallel for firstprivate(dist, loc, counter, offset, counterScan, updateNN, dstKeys, dstKeys1)
390 for (int i = 0; i < initdist.size(); ++i)
391 {
392 const Vortex2D& vtxi = vtx[mcdata[i].originNumber];
393
394 int cntr = 0;
395 int search = i - 1;
396 std::vector<std::pair<double, size_t>>& fillPosition = ((sdvig == 0) ? initdist[mcdata[i].originNumber] :
dist);
397
398 while ((cntr < k) && (search >= 0))
399 {
400 const Vortex2D& vtxk = vtx[mcdata[search].originNumber];
401 bool flagExit, check;
402 double d2;
403 std::tie(flagExit, check) =
calcCheck(vtxi, vtxk, maxG, cSP, cRBP, epsCol, type, d2);
404
405 if (flagExit)
406 break;
407
408 if ((mcdata[search].originNumber > mcdata[i].originNumber) && check)
409 fillPosition[cntr++] = { d2, mcdata[search].originNumber };
410
411 --search;
412 }
413
414 for (int w = cntr; w < k; ++w)
415 fillPosition[cntr++] = { 100000000.0, 0 };
416
417 search = i + 1;
418
419 while ((cntr < 2 * k) && (search < mcdata.size()))
420 {
421 const Vortex2D& vtxk = vtx[mcdata[search].originNumber];
422 bool flagExit, check;
423 double d2;
424 std::tie(flagExit, check) =
calcCheck(vtxi, vtxk, maxG, cSP, cRBP, epsCol, type, d2);
425
426 if (flagExit)
427 break;
428
429 if ((mcdata[search].originNumber > mcdata[i].originNumber) && check)
430 fillPosition[cntr++] = { d2, mcdata[search].originNumber };
431
432 ++search;
433 }
434 for (int w = cntr; w < 2 * k; ++w)
435 fillPosition[cntr++] = { 100000000.0, 0 };
436
437 if (sdvig == 0)
438 {
439 newSort(initdist[mcdata[i].originNumber], dstKeys1);
440 initdist[mcdata[i].originNumber].resize(k);
441 }
442 else
443 {
445 newMerge(mcdata[i].originNumber, initdist[mcdata[i].originNumber], dist, loc, counter, offset, counterScan, updateNN);
446 newSort(initdist[mcdata[i].originNumber], dstKeys);
447 }
448 }
449
450 time[4] = omp_get_wtime();
451
452 tFinish[sdvig] = omp_get_wtime();
453
454 }
455}
Класс, опеделяющий двумерный вихревой элемент
R dist(const numvector< T, n > &x, const numvector< P, n > &y)
Вычисление расстояния между двумя точками
void newMerge(int iii, std::vector< std::pair< double, size_t > > ¤tNN, const std::vector< std::pair< double, size_t > > &candidateNN, std::vector< size_t > &loc, std::vector< size_t > &counter, std::vector< size_t > &offset, std::vector< size_t > &counterScan, std::vector< std::pair< double, size_t > > &updateNN)
void RSort_Parallel(TParticleCode *m, TParticleCode *m_temp, unsigned int n, unsigned int *s)
Сортировка массива из мортоновский кодов
unsigned int Morton2D(const Point2D &r)
std::pair< bool, bool > calcCheck(const Vortex2D &vtxi, const Vortex2D &vtxk, double maxG, double cSP, double cRBP, double epsCol, int type, double &d2)
void newSort(std::vector< std::pair< double, size_t > > &mass, std::vector< std::pair< double, size_t > > &dstKeys)