Структура, соответствующая мортоновскому дереву More...

#include <TreeBH.h>

Collaboration diagram for BH::MortonTree:

[legend]

Public Member Functions
void	getStatistics (int &numParticles, int &treshold, std::pair< int, int > &partLevel, std::pair< int, int > &leafsLevel, int &lowLevelCells) const

double &	getBinom (int n, int k)

	MortonTree (const params &prm_, int maxTreeLevel_, std::vector< PointsCopy > &points)

void	MakeRootMortonTree ()
	Построение корня дерева и задание его общих параметров

void	BuildMortonInternalTree ()
	Построение внутренних ячеек дерева

void	BuildMortonParticlesTree ()
	Построение верхушек дерева — отдельных частиц

void	FillMortonLowCells (int cell=0)
	Функция заполнения списка нижних ячеек дерева (до которых производится обход)

void	FillMortonLowCellsA ()

void	CalculateMortonTreeParams (int cell, int omplevel)
	Вычисление параметров дерева (циркуляций и высших моментов)

void	CalcLocalCoeffToLowLevel (int lowCell, std::unique_ptr< MortonTree > &treeInf, int fromWho, bool calcCloseTrees)

void	CalcVeloBiotSavart (int lowCell, std::unique_ptr< MortonTree > &treeInf)

void	CalcVeloTaylorExpansion (int lowCell)

void	AddVelo (int lowCell)

void	MinMaxCellsSearch (int numCell)

void	MinMaxCellsCalc (BH::treeCellT &cell)

void	GabsPyr (int numCell)

void	ShareGabs (int numCell)

Public Attributes
std::vector< TParticleCode >	mortonCodes
	Список мортоновских кодов, вычисляемых по положениям частиц

std::vector< double >	iFact
	Спикок обратных факториалов целых чисел

std::vector< double >	binomCft

const int	maxTreeLevel
	Максимальная глубина дерева при его обходе

std::vector< int >	mortonLowCells
	Вектор индесков ячеек нижнего уровня в дереве

std::vector< treeCellT >	mortonTree
	Мортоновское дерево Размер дерева равен сумме числа внутренних вершин (число частиц без единицы) и числа частиц резервируется на 1 штуку больше для удобства работы ячейки [0 ... pos.size()-2] отвечают внутренним узлам, ячейка с идексом pos.size()-1 не используется ячейки [pos.size() ... 2*pos.size()-1] отвечают частицам

Private Member Functions
int	Delta (int i, int j) const

int	PrefixLength (int cell) const

std::pair< unsigned int, int >	Prefix (int cell) const

void	SetCellGeometry (int cell)
	Функция вычисления геометрических параметров внутренней ячейки дерева

void	BuildInternalTreeCell (int i)
	Функция построения i-й внутренней ячейки дерева

void	RSort_Node3 (TParticleCode m, TParticleCode m_temp, unsigned int n)
	Поразрядная сортировка мортоновских кодов

void	RSort_step3 (TParticleCode source, TParticleCode dest, unsigned int n, unsigned int *offset, unsigned char sortable_bit)
	Вспомогательная функция к RSort_Node3.

void	RSort_Parallel (TParticleCode m, TParticleCode m_temp, unsigned int n, unsigned int *s)
	Параллельная версия поразрядной сортировки мортоновских кодов

Static Private Member Functions
static unsigned int	ExpandBits (unsigned int v)
	"Разрежение" двоичного представления беззнакового целого, вставляя по одному нулику между всеми битами

static unsigned int	Morton2D (const Point2D &r)

static std::pair< Point2D, Point2D >	FindEnclosingRectangleOld (const std::vector< PointsCopy > &points)
	Поиск габаритного прямоугольника системы точек (для OpenMP 2.0)

Private Attributes
const params &	prm
	Ссылка на параметры, загружаемые из файла

std::vector< PointsCopy > &	pointsCopy
	Ссылка на список обобщенных вихрей, по которым строится дерево

double	iQuadSideVar

Point2D	lowLeft
	Положение нижнего левого угла дерева

Detailed Description

Структура, соответствующая мортоновскому дереву

Author: Марчевский Илья Константинович; Сокол Ксения Сергеевна; Рятина Евгения Павловна; Колганова Александра Олеговна \Version 1.14

Date: 6 марта 2026 г.

Definition at line 146 of file TreeBH.h.

Constructor & Destructor Documentation

◆ MortonTree()

BH::MortonTree::MortonTree	(	const params &	prm_,
		int	maxTreeLevel_,
		std::vector< PointsCopy > &	points
	)

Definition at line 49 of file TreeBH.cpp.

        : prm(prm_), pointsCopy(points), maxTreeLevel(maxTreeLevel_)
    {
        mortonTree.resize(2 * points.size());
        for (auto& c : mortonTree)
        {
            c.mom.resize(prm.order + 1);
            c.E.resize(prm.order);
        }
 
        mortonCodes.resize(points.size());
        mortonLowCells.reserve(points.size());
        
 
        //Вычисляем факториалы и биномиальные коэффиценты
        iFact.resize(prm.order + 1);
        binomCft.resize((prm.order + 1) * (prm.order + 1));
        
        iFact[0] = 1.0;
        for (int i = 1; i <= prm.order; ++i)
            iFact[i] = iFact[i - 1] / i;
        getBinom(0, 0) = 1.0;
        for (int i = 1; i <= prm.order; ++i) 
        {
            getBinom(i, 0) = 1.0;
            getBinom(i, i) = 1.0;
            for (int j = 1; j < i; ++j)
                getBinom(i, j) = getBinom(i - 1, j - 1) + getBinom(i - 1, j);
        }//for i
    }//MortonTree(...)

Here is the call graph for this function:

Member Function Documentation

◆ AddVelo()

void BH::MortonTree::AddVelo ( int lowCell )

Добавление влияния дальней зоны для правой части

Parameters

[in] lowCell — индекс ячейки нижнего уровня, в которой рассчитывается влияние

◆ BuildInternalTreeCell()

void BH::MortonTree::BuildInternalTreeCell ( int i )

private

Функция построения i-й внутренней ячейки дерева

Definition at line 264 of file TreeBH.cpp.

    {
        const int n = (int)pointsCopy.size();
        int d = sign(Delta(i, i + 1) - Delta(i, i - 1));
        int delta_min = Delta(i, i - d);
        
        int Lmax = 2;
        while (Delta(i, i + Lmax * d) > delta_min)
            Lmax *= 2;
 
        int L = 0;
        for (int t = (Lmax >> 1); t >= 1; t >>= 1)      
            if (Delta(i, i + (L + t) * d) > delta_min)
                L += t;
        
        int j = i + L * d;
        int delta_node = Delta(i, j);
        int s = 0;
        for (int p = 1, t = ceilhalf(L); L > (1 << (p - 1)); ++p, t = ceilpow2(L, p))
        {
            int dl = Delta(i, i + (s + t) * d);
            if (dl > delta_node)
                s += t;
        }//for p
        int gamma = i + s * d + std::min(d, 0);
 
        auto min_max = std::minmax(i, j);
 
        const int& left = gamma;
        const int& right = gamma + 1;
 
        treeCellT& treeCell = mortonTree[i];
 
        // Левый потомок - лист или внутренний узел
        bool ifLeftParticle = (min_max.first == gamma);
        treeCell.child[0] = ifLeftParticle ? n + left : left;
        mortonTree[treeCell.child[0]].range = { min_max.first, gamma };
        mortonTree[treeCell.child[0]].particle = ifLeftParticle;
        mortonTree[treeCell.child[0]].parent = i;
 
        // Правый потомок - лист или внутренний узел
        bool ifRightParticle = (min_max.second == gamma + 1);
        treeCell.child[1] = ifRightParticle ? n + right : right;
        mortonTree[treeCell.child[1]].range = { gamma + 1, min_max.second };
        mortonTree[treeCell.child[1]].particle = ifRightParticle;
        mortonTree[treeCell.child[1]].parent = i;
    }//BuildInternalTreeCell(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ BuildMortonInternalTree()

void BH::MortonTree::BuildMortonInternalTree ( )

Построение внутренних ячеек дерева

Definition at line 314 of file TreeBH.cpp.

    {
#pragma omp parallel for schedule(dynamic, 1000)
        for (int q = 0; q < (int)pointsCopy.size() - 1; ++q)
            BuildInternalTreeCell(q);
        
        mortonTree[pointsCopy.size() - 1].leaf = false;
        mortonTree[pointsCopy.size() - 1].particle = false;
 
#pragma omp parallel for schedule(dynamic, 1000)
        for (int q = 0; q < (int)pointsCopy.size() - 1; ++q)
            SetCellGeometry(q);
    }//BuildMortonInternalTree()

Here is the call graph for this function:

◆ BuildMortonParticlesTree()

void BH::MortonTree::BuildMortonParticlesTree ( )

Построение верхушек дерева — отдельных частиц

Definition at line 330 of file TreeBH.cpp.

    {
#pragma omp parallel for        
        for (int q = 0; q < (int)pointsCopy.size(); ++q)
        {
            auto& pt = mortonCodes[q];
 
            auto& cell = mortonTree[pointsCopy.size() + q];
            cell.center = pt.r;
 
            cell.leaf = true;
            cell.size = { 0.0, 0.0 };
        }//for q
    }//BuildMortonParticlesTree(...)

Here is the call graph for this function:

◆ CalcLocalCoeffToLowLevel()

void BH::MortonTree::CalcLocalCoeffToLowLevel	(	int	lowCell,
		std::unique_ptr< MortonTree > &	treeInf,
		int	fromWho,
		bool	calcCloseTrees
	)

Расчет коэффицентов локального разложения для нижних ячеек

Parameters

[in]	lowCell	индекс ячейки нижнего уровня, в которой будет рассчитываться влияние
[in]	treeInf	ссылка на влияющее дерево
[in]	fromWho	индекс влияющей ячейки во влияющем дереве
[in]	calcCloseTrees	признак заполнения списка ближних ячеек (нижнего уровня) к ячейке lowCell

Definition at line 436 of file TreeBH.cpp.

    {       
        auto& lt = mortonTree[lowCell];
        //if (calcCloseTrees)
        //  lt.closeCells.resize(0);
 
        if (treeInf.get() != this || lowCell != fromWho)
        {
            auto& wh = treeInf->mortonTree[fromWho];
 
            double h, h0;
            h = wh.size[0] + wh.size[1];
            h0 = lt.size[0] + lt.size[1];
 
            Point2D& r0 = lt.center;
            Point2D& r1 = wh.center;
 
            Point2D rr = r0 - r1;
 
            double crit2 = rr.length2();
            // если выполнен критерий дальности => считаем коэффициенты
            if ((crit2 >= sqr((h0 + h + 2.0 * prm.eps) / prm.theta)) && (wh.level <= prm.NumOfLevelsVortex))
            {
                if (wh.range[1] == wh.range[0])
                    lt.closeCells.push_back(fromWho);
                else
                {
                    Point2D  rr = r0 - r1;
                    rr *= 1.0 / rr.length2();
 
                    Point2D varTheta = rr;
                    for (int diag = 0; diag < prm.order; ++diag)
                    {
                        for (int q = 0; q <= diag; ++q)
                            lt.E[q] += ((q % 2 ? -1.0 : 1.0) * iFact[diag - q]) * multzA(varTheta, wh.mom[diag - q]);
                        
                        varTheta = (diag + 1) * multz(varTheta, rr);
                    }
                    lt.E[0] += iFact[prm.order] * multzA(varTheta, wh.mom[prm.order]);
                }
            }//if crit2
            else // если не выполнен критерий, то рекурсия 
            {
                if ((!wh.leaf) && (wh.level < prm.NumOfLevelsVortex))
                {
                    CalcLocalCoeffToLowLevel(lowCell, treeInf, wh.child[0], calcCloseTrees);
                    CalcLocalCoeffToLowLevel(lowCell, treeInf, wh.child[1], calcCloseTrees);
                }
                else if (calcCloseTrees)
                {   
                    lt.closeCells.push_back(fromWho);                   
                }
            }//else if
        }//if (lowTree != this)
        else if (calcCloseTrees)
        {
            lt.closeCells.push_back(lowCell); //себя тоже добавляем в ближнюю зону          
        }
    }//CalcLocalCoeffToLowLevel(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ CalculateMortonTreeParams()

void BH::MortonTree::CalculateMortonTreeParams	(	int	cell,
		int	omplevel
	)

Вычисление параметров дерева (циркуляций и высших моментов)

Расчет мультипольных моментов всех ячеек дерева

Parameters

[in]	cell	ячейка для вычисления мультипольных моментов
[in]	omplevel	текущий уровень вложенности OpenMP

Definition at line 378 of file TreeBH.cpp.

    {
        const int maxLevelOmp = int(log2(omp_get_max_threads()));
 
        auto& cl = mortonTree[cell];
        if (!cl.particle)
        {
#pragma omp parallel /*default(none)*/ shared(cl, omplevel) num_threads(2) if (omplevel < maxLevelOmp + 1)
            {               
                std::vector<Point2D> h(prm.order);
#pragma omp for 
                for (int s = 0; s < 2; ++s)
                    CalculateMortonTreeParams(cl.child[s], omplevel + 1);
 
#pragma omp single
                {                   
                    for (auto& m : cl.mom)
                        m.toZero();
                                        
                    for (int s = 0; s < 2; ++s)//цикл по двум потомкам
                    {
                        auto& chld = mortonTree[cl.child[s]];
                        const Point2D& g = chld.mom[0];
                        cl.mom[0] += g;
 
                        
                        h[0] = (chld.center) - (cl.center);
                        for (int i = 1; i < prm.order; ++i)
                            h[i] = multz(h[i - 1], h[0]);
                        
                        for (int p = 1; p <= prm.order; ++p) 
                        {
                            Point2D shiftMom = chld.mom[p];
                            for (int k = 1; k <= p; ++k)
                            {
                                shiftMom += getBinom(p, k) *  multz(chld.mom[p - k], h[k - 1]);
                                //ADDOP(2);
                            }//for k
                            cl.mom[p] += shiftMom;
                        }//for m
                        
                    }//for s
                }//single               
            }//parallel
        }//if(!particle)
        else
        {
            //расчет моментов ячейки нижнего уровня для вихрей или для панелей
                //C учетом того, что дерево строится до частиц --- только монопольный момент отличен от нуля
                cl.mom[0][0] = mortonCodes[cl.range[0]].g;
                cl.mom[0][1] = 0.0;
 
                for (int p = 1; p <= prm.order; ++p)
                    cl.mom[p].toZero();
        }//else
    }//CalculateMortonTreeParams(...)       

Here is the call graph for this function:

Here is the caller graph for this function:

◆ CalcVeloBiotSavart()

void BH::MortonTree::CalcVeloBiotSavart	(	int	lowCell,
		std::unique_ptr< MortonTree > &	treeInf
	)

Расчет влияния от ближних ячеек по формуле Био - Савара

Parameters

[in]	lowCell	индекс ячейки нижнего уровня, в которой рассчитывается влияние
[in]	treeInf	влияющее дерево (чтобы брать ближние ячейки по индексу из нужного массива)

Definition at line 499 of file TreeBH.cpp.

    {
        auto& lc = mortonTree[lowCell];
        for (int i = lc.range[0]; i <= lc.range[1]; ++i)        
        {
            //Локальные переменные для цикла
            Point2D velI, velLinI;
            double dst2eps;
            double r2;          
 
            for (size_t k = 0; k < lc.closeCells.size(); ++k)
            {
                velI.toZero();
                velLinI.toZero();           
 
                const Point2D& posI = mortonCodes[i].r; 
 
                auto& rg = treeInf->mortonTree[lc.closeCells[k]].range;
                for (int j = rg[0]; j <= rg[1]; ++j)
                {
                    auto& pt = treeInf->mortonCodes[j];
 
                    const Point2D& posJ = pt.r;
                    const double& gamJ = pt.g;
                    r2 = (posI - posJ).length2();
                    dst2eps = std::max(r2, prm.eps2);
                    velI += (gamJ / dst2eps) * (posI - posJ).kcross();          
                }//for j
 
                pointsCopy[mortonCodes[i].originNumber].veloCopy += velI;       
            }//for k 
 
        }//for i
    }//CalcVeloBiotSavart(...)

Here is the call graph for this function:

◆ CalcVeloTaylorExpansion()

void BH::MortonTree::CalcVeloTaylorExpansion ( int lowCell )

Расчет влияния от дальней зоны при помощи Тейлоровского разложения

Parameters

[in] lowCell — индекс ячейки нижнего уровня, в которой рассчитывается влияние

Definition at line 536 of file TreeBH.cpp.

    {
        auto& lc = mortonTree[lowCell];
        for (int i = lc.range[0]; i <= lc.range[1]; ++i)
        {
            Point2D deltaPos = mortonCodes[i].r - lc.center;
            Point2D v = lc.E[0];
 
            Point2D hh = deltaPos;
            for (int k = 1; k < prm.order - 1; ++k) 
            {
                v += iFact[k] * multzA(lc.E[k], hh);
                hh = multz(hh, deltaPos);
            }//for k
            v += iFact[prm.order-1] * multzA(lc.E[prm.order-1], hh);
 
            pointsCopy[mortonCodes[i].originNumber].veloCopy += Point2D({ -v[1], v[0] });   
 
        }//for i
    }//CalcVeloTaylorExpansion(...)

◆ Delta()

int BH::MortonTree::Delta	(	int	i,
		int	j
	)		const

private

Функция вычисления длины общей части префиксов двух (а значит - и диапазона) частиц

Parameters

[in]	i	индекс одной из частиц в отсортированном по мортоновским кодам списке
[in]	j	индекс второй из частиц в отсортированном по мортоновским кодам списке

Returns: длину общей части мортоновского кода

Definition at line 187 of file TreeBH.cpp.

    {
        if ((j < 0) || (j > (int)pointsCopy.size() - 1))
            return -1;
 
        if (i > j)
            std::swap(i, j);
 
        //if ((i < 0) || (j > n-1))
        //    exit(111);
 
        const unsigned int& ki = mortonCodes[i].key;
        const unsigned int& kj = mortonCodes[j].key;
        
        //Поиск номера самого старшего ненулевого бита в числе c 
        int count = 0;
        for (unsigned int c = (ki ^ kj); c; c >>= 1, ++count);
 
        if ((!count) && (i != j))
        {
            int addCount = 0;
            //единички к номерам i и j добавлены для совместимости с Wolfram Mathematica, 
            //для кода они не важны, но дерево без них почти наверняка построится по-другому        
            for (unsigned int add = ((i + 1) ^ (j + 1)); add; add >>= 1, ++addCount);
            return 2 * codeLength + (2 * codeLength - addCount);
        }//if ((!count) && (i != j))
 
        return (2 * codeLength - count);
    }//Delta(...)

Here is the caller graph for this function:

◆ ExpandBits()

unsigned int BH::MortonTree::ExpandBits ( unsigned int v )

staticprivate

"Разрежение" двоичного представления беззнакового целого, вставляя по одному нулику между всеми битами

Definition at line 559 of file TreeBH.cpp.

    {
        // вставит 1 нуль
        v = (v | (v << 8)) & 0x00FF00FF;      //  00000000`00000000`abcdefgh`ijklmnop 
        //                                      | 00000000`abcdefgh`ijklmnop`00000000
        //                                      = 00000000`abcdefgh`XXXXXXXX`ijklmnop
        //                                      & 00000000`11111111`00000000`11111111
        //                                      = 00000000`abcdefgh`00000000`ijklmnop
 
        v = (v | (v << 4)) & 0x0F0F0F0F;      //  00000000`abcdefgh`00000000`ijklmnop 
        //                                      | 0000abcd`efgh0000`0000ijkl`mnop0000
        //                                      = 0000abcd`XXXXefgh`0000ijkl`XXXXmnop
        //                                      & 00001111`00001111`00001111`00001111
        //                                      = 0000abcd`0000efgh`0000ijkl`0000mnop
 
        v = (v | (v << 2)) & 0x33333333;      //  0000abcd`0000efgh`0000ijkl`0000mnop 
        //                                      | 00abcd00`00efgh00`00ijkl00`00mnop00
        //                                      = 00abXXcd`00efXXgh`00ijXXkl`00mnXXop
        //                                      & 00110011`00110011`00110011`00110011
        //                                      = 00ab00cd`00ef00gh`00ij00kl`00mn00op
 
        v = (v | (v << 1)) & 0x55555555;      //  00ab00cd`00ef00gh`00ij00kl`00mn00op 
        //                                      | 0ab00cd0`0ef00gh0`0ij00kl0`0mn00op0
        //                                      = 0aXb0cXd`0eXf0gXh`0iXj0kXl`0mXn0oXp
        //                                      & 01010101`01010101`01010101`01010101
        //                                      = 0a0b0c0d`0e0f0g0h`0i0j0k0l`0m0n0o0p
        return v;
    }

Here is the caller graph for this function:

◆ FillMortonLowCells()

void BH::MortonTree::FillMortonLowCells ( int cell = 0 )

Функция заполнения списка нижних ячеек дерева (до которых производится обход)

Definition at line 348 of file TreeBH.cpp.

    {
        if (mortonTree[cell].leaf)
            mortonLowCells.push_back(cell);
        else
        {
            FillMortonLowCells(mortonTree[cell].child[0]);
            FillMortonLowCells(mortonTree[cell].child[1]);
        }//else
    }//FillMortonLowCells(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ FillMortonLowCellsA()

void BH::MortonTree::FillMortonLowCellsA ( )

Definition at line 362 of file TreeBH.cpp.

    {
#pragma omp parallel
        {
            std::vector<int> good_matches_private;
#pragma omp for nowait
            for (int cell = 0; cell < (int)mortonTree.size(); ++cell)           
                if (!mortonTree[mortonTree[cell].parent].leaf && mortonTree[cell].leaf)             
                    good_matches_private.push_back(cell);               
#pragma omp critical
            mortonLowCells.insert(mortonLowCells.end(), good_matches_private.begin(), good_matches_private.end());
        }//parallel
    }//FillMortonLowCellsA()

◆ FindEnclosingRectangleOld()

std::pair< Point2D, Point2D > BH::MortonTree::FindEnclosingRectangleOld ( const std::vector< PointsCopy > & points )

staticprivate

Поиск габаритного прямоугольника системы точек (для OpenMP 2.0)

Definition at line 82 of file TreeBH.cpp.

    {
        double shared_minx = 1e+10, shared_maxx = -1e+10, shared_miny = 1e+10, shared_maxy = -1e+10;
#pragma omp parallel 
        {
            double minx = 1e+10, maxx = -1e+10, miny = 1e+10, maxy = -1e+10;
#pragma omp for nowait
            for (int ii = 0; ii < (int)points.size(); ++ii)
            {
                const Point2D& r = points[ii].r();
                minx = std::min(r[0], minx);
                maxx = std::max(r[0], maxx);
                miny = std::min(r[1], miny);                
                maxy = std::max(r[1], maxy);
            }//for ii
#pragma omp critical 
            {
                shared_minx = std::min(shared_minx, minx);
                shared_maxx = std::max(shared_maxx, maxx);
                shared_miny = std::min(shared_miny, miny);
                shared_maxy = std::max(shared_maxy, maxy);
            }//critical
        }//parallel
        
        return { {shared_minx, shared_miny}, {shared_maxx, shared_maxy} };
    }//FindEnclosingRectangleOld(...)

Here is the caller graph for this function:

◆ GabsPyr()

void BH::MortonTree::GabsPyr ( int numCell )

Definition at line 870 of file TreeBH.cpp.

    {
        auto& cell = mortonTree[numCell];
        if (cell.hasGabs)
            return;
 
        int ch0 = cell.child[0];
        int ch1 = cell.child[1];
        if (!mortonTree[ch0].hasGabs)
            GabsPyr(ch0);
        if (!mortonTree[ch1].hasGabs)
            GabsPyr(ch1);
        
        ShareGabs(numCell/*, cell.minx_miny, cell.maxx_maxy*/);
    }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getBinom()

double & BH::MortonTree::getBinom	(	int	n,
		int	k
	)

inline

Извлечение биномиального коэффициента из массива binomCft

Returns: ссылку на C_n^k

Definition at line 238 of file TreeBH.h.

        {
            return binomCft[n * (prm.order + 1) + k];
        }

Here is the caller graph for this function:

◆ getStatistics()

void BH::MortonTree::getStatistics	(	int &	numParticles,
		int &	treshold,
		std::pair< int, int > &	partLevel,
		std::pair< int, int > &	leafsLevel,
		int &	lowLevelCells
	)		const

Сбор статистики по дереву

Parameters

[out]	numParticles	число частиц в дереве
[out]	treshold	уровень, до которого производится обход
[out]	partLevel	наименьший и наибольший уровни, на которых находятся частицы
[out]	leafsLevel	наименьший и наибольший уровни, на которых ячейки, до которых производится обход
[out]	lowLevelCells	число ячеек "нижнего" уровня (до которых производится обход)

Definition at line 749 of file TreeBH.cpp.

    {
        numParticles = (int)pointsCopy.size();
        treshold = maxTreeLevel;
 
        int leafsLevelMin = (int)1e+9;
        int leafsLevelMax = 0;
        for (int q = 0; q < numParticles - 1; ++q)
        {
            if (!mortonTree[mortonTree[q].parent].leaf)
            {
                if (mortonTree[q].leaf)
                {
                    int level = mortonTree[q].level;
 
                    if (level > leafsLevelMax)
                        leafsLevelMax = level;
 
                    if (level < leafsLevelMin)
                        leafsLevelMin = level;
 
                    continue;
                }
 
                if ((mortonTree[mortonTree[q].child[0]].particle) || (mortonTree[mortonTree[q].child[1]].particle))
                {
                    int level = mortonTree[q].level + 1;
                    if (level > leafsLevelMax)
                        leafsLevelMax = level;
 
                    if (level < leafsLevelMin)
                        leafsLevelMin = level;
 
                    continue;
                }
            }
        }
        leafsLevel = { leafsLevelMin, leafsLevelMax };
 
 
        int partLevelMin = (int)1e+9;
        int partLevelMax = 0;
        for (int q = 0; q < numParticles - 1; ++q)
        {
            if ((mortonTree[mortonTree[q].child[0]].particle) || (mortonTree[mortonTree[q].child[1]].particle))
            {
                int level = mortonTree[q].level;
                if (level > partLevelMax)
                    partLevelMax = level;
 
                if (level < partLevelMin)
                    partLevelMin = level;
            }
        }
        partLevel = { partLevelMin + 1, partLevelMax + 1 };
        lowLevelCells = (int)mortonLowCells.size();
    }

◆ MakeRootMortonTree()

void BH::MortonTree::MakeRootMortonTree ( )

Построение корня дерева и задание его общих параметров

Definition at line 140 of file TreeBH.cpp.

    {
        auto gabs = FindEnclosingRectangleOld(pointsCopy);
 
        double iLeft = gabs.first[0];
        double iBottom = gabs.first[1];
        
        double iRight = gabs.second[0];
        double iTop = gabs.second[1];
 
        Point2D posCentre = { 0.5 * (iLeft + iRight), 0.5 * (iBottom + iTop) };
 
        double quadSide = std::max(iRight - iLeft, iTop - iBottom);
        iQuadSideVar = 1.0 / (quadSide * (1.0 + 1.0 / ((1 << codeLength) - 1)));
        lowLeft = posCentre - 0.5*quadSide*Point2D{ 1.0, 1.0 };
                
        //iQuadSideVar = 1.0;
        //lowLeft = Point2D{ 0.0, 0.0 };
                
#pragma omp parallel for schedule(dynamic, 1000)
        for (int i = 0; i < (int)pointsCopy.size(); ++i)
        {
            Point2D locR = (pointsCopy[i].r() - lowLeft) * iQuadSideVar;
            unsigned int mortonCode = Morton2D(locR);
                
            auto& mc = mortonCodes[i];
            mc.key = mortonCode;
            mc.originNumber = i;
            mc.r = pointsCopy[i].r();
            mc.g = pointsCopy[i].g();
        }//for i
 
        //std::cout << "omp_threads = " << omp_get_max_threads() << std::endl;
 
        //RSort_Node3(mortonCodes.data(), mortonCodes_temp.data(), (int)pointsCopy.size());             
        
        //Временные массивы для сортировки      
        std::vector<unsigned int> s(256 * omp_get_max_threads());
        std::vector<TParticleCode> mortonCodes_temp(pointsCopy.size());
        RSort_Parallel(mortonCodes.data(), mortonCodes_temp.data(), (int)pointsCopy.size(), s.data());
 
        mortonTree[0].range = { 0, (int)pointsCopy.size() - 1 };
        mortonTree[0].particle = false;
    }//MakeRootMortonTree()

Here is the call graph for this function:

◆ MinMaxCellsCalc()

void BH::MortonTree::MinMaxCellsCalc ( BH::treeCellT & cell )

Definition at line 829 of file TreeBH.cpp.

    {
        int spoint = cell.range[0];
        int epoint = cell.range[1];
        
        double minx = 1e+10, maxx = -1e+10, miny = 1e+10, maxy = -1e+10;
 
        for (int ii = spoint; ii <= epoint; ++ii)
        {
            const Point2D& r = mortonCodes[ii].r;
            minx = std::min(r[0], minx);
            maxx = std::max(r[0], maxx);
            miny = std::min(r[1], miny);
            maxy = std::max(r[1], maxy);
        }
        cell.minx_miny[0] = minx;
        cell.minx_miny[1] = miny;
        cell.maxx_maxy[0] = maxx;
        cell.maxx_maxy[1] = maxy;
 
        cell.center = (cell.minx_miny + cell.maxx_maxy) * 0.5;
        cell.size = cell.maxx_maxy - cell.minx_miny;
 
    }

Here is the caller graph for this function:

◆ MinMaxCellsSearch()

void BH::MortonTree::MinMaxCellsSearch ( int numCell )

Definition at line 809 of file TreeBH.cpp.

    {
        auto& cell = mortonTree[numCell];
        if (cell.leaf)
        {
            MinMaxCellsCalc(cell);
            cell.hasGabs = true;
        }
        else
        {
            std::cout << "!!!!!!!!!!!!" << std::endl;
 
            int ch0 = cell.child[0];
            int ch1 = cell.child[1];
            MinMaxCellsSearch(ch0);
            MinMaxCellsSearch(ch1);
            cell.hasGabs = false;
        }
    }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Morton2D()

unsigned int BH::MortonTree::Morton2D ( const Point2D & r )

staticprivate

Мортоновский код для пары из чисел типа double Исходные числа (r[0], r[1]) - строго в диапазоне [0, 1 / (1.0 + 1/(2^codeLength - 1) ))

Definition at line 591 of file TreeBH.cpp.

    {
        /*
        if (x < 0) std::cout << "x*... < 0" << std::endl;
        if (y < 0) std::cout << "y*... < 0" << std::endl;
 
        if (x * twoPowCodeLength > twoPowCodeLengthMinus)
            std::cout << "x*... > ..." << std::endl;
 
        if (y * twoPowCodeLength > twoPowCodeLengthMinus)
            std::cout << "y*... > ..." << std::endl;
 
        x = std::min(std::max(x * twoPowCodeLength, 0.0), twoPowCodeLengthMinus);
        y = std::min(std::max(y * twoPowCodeLength, 0.0), twoPowCodeLengthMinus);
 
        unsigned int xx = expandBits((unsigned int)x);
        unsigned int yy = expandBits((unsigned int)y);
        */
        const Point2D& rscale = twoPowCodeLength * r;
        const unsigned int& xx = ExpandBits((unsigned int)(rscale[0]));
        const unsigned int& yy = ExpandBits((unsigned int)(rscale[1]));
        return yy | (xx << 1);
    }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Prefix()

std::pair< unsigned int, int > BH::MortonTree::Prefix ( int cell ) const

private

Функция вычисления общего префикса и его длины для всех частиц в конкретной ячейке

Parameters

[in] cell индекс ячейки дерева в линейном массиве

Returns: пару из общей части мортоновских кодов всех частиц в данной ячейке и ее длины

Definition at line 227 of file TreeBH.cpp.

    {
        int length = PrefixLength(cell);
        unsigned int el = mortonCodes[mortonTree[cell].range[0]].key;
        return { el >> (2 * codeLength - length), length };
    }//Prefix(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PrefixLength()

int BH::MortonTree::PrefixLength ( int cell ) const

private

Функция вычисления длины общей части префиксов всех частиц в конкретной ячейке

Parameters

[in] cell индекс ячейки дерева в линейном массиве

Returns: длину общей части мортоновских кодов всех частиц в данной ячейке

Definition at line 219 of file TreeBH.cpp.

    {
        const auto& range = mortonTree[cell].range;
        return std::min(Delta(range[0], range[1]), 2 * codeLength);
    }//PrefixLength(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ RSort_Node3()

void BH::MortonTree::RSort_Node3	(	TParticleCode *	m,
		TParticleCode *	m_temp,
		unsigned int	n
	)

private

Поразрядная сортировка мортоновских кодов

Definition at line 630 of file TreeBH.cpp.

    {
        // Заводим массив корзин
        unsigned int s[sizeof(m->key) * 256] = { 0 };
        // Заполняем массив корзин для всех разрядов
        for (unsigned int i = 0; i < n; ++i)
        {
            unsigned int key = m[i].key;
            for (unsigned int digit = 0; digit < sizeof(m->key); digit++)
                ++s[((key >> (digit << 3)) & 0xFF) + (digit << 8)];
        }
 
        // Пересчитываем смещения для корзин
        for (unsigned int digit = 0; digit < sizeof(m->key); digit++)
        {
            unsigned int off = 0;
            for (unsigned int i = 0; i < 256; i++)
            {
                auto& rs = s[i + (digit << 8)];
                unsigned int value = rs;
                rs = off;
                off += value;
            }
        }
 
        // Вызов сортировки по битам от младших к старшим (LSD)
        for (unsigned int digit = 0; digit < sizeof(m->key); digit++)
        {
            RSort_step3(m, m_temp, n, &s[digit << 8], digit);
            std::swap(m, m_temp);
            //TParticleCode* temp = m;
            //m = m_temp;
            //m_temp = temp;
        }
 
        // Если ключ структуры однобайтовый, копируем отсортированное в исходный массив
        if (sizeof(m->key) == 1)
        {
            std::swap(m, m_temp);
            //TParticleCode* temp = m;
            //m = m_temp;
            //m_temp = temp;
            memcpy(m, m_temp, n * sizeof(TParticleCode));
        }
    }

Here is the call graph for this function:

◆ RSort_Parallel()

void BH::MortonTree::RSort_Parallel	(	TParticleCode *	m,
		TParticleCode *	m_temp,
		unsigned int	n,
		unsigned int *	s
	)

inlineprivate

Параллельная версия поразрядной сортировки мортоновских кодов

Definition at line 677 of file TreeBH.cpp.

    {
        if (n == 0)
            return;
        // Количество задействованных потоков
        unsigned char threads = omp_get_max_threads();
        //std::cout << "threads = " << (int)threads << std::endl;
 
#pragma omp parallel num_threads(threads)
        {
            TParticleCode* source = m;
            TParticleCode* dest = m_temp;
            unsigned int l = omp_get_thread_num();
            unsigned int div = n / omp_get_num_threads();
            unsigned int mod = n % omp_get_num_threads();
            unsigned int left_index = l < mod ? (div + (mod == 0 ? 0 : 1)) * l : n - (omp_get_num_threads() - l) * div;
            unsigned int right_index = left_index + div - (mod > l ? 0 : 1);
 
            for (unsigned int digit = 0; digit < sizeof(m->key); ++digit)
            {
                unsigned int s_sum[256] = { 0 };
                unsigned int s0[256] = { 0 };
                unsigned char* b1 = (unsigned char*)&source[right_index].key;
                unsigned char* b2 = (unsigned char*)&source[left_index].key;
                while (b1 >= b2)
                {
                    ++s0[*(b1 + digit)];
                    b1 -= sizeof(TParticleCode);
                }
                for (unsigned int i = 0; i < 256; i++)
                {
                    s[i + 256 * l] = s0[i];
                }
 
#pragma omp barrier
                for (unsigned int j = 0; j < threads; j++)
                {
                    for (unsigned int i = 0; i < 256; i++)
                    {
                        s_sum[i] += s[i + 256 * j];
                        if (j < l)
                        {
                            s0[i] += s[i + 256 * j];
                        }
                    }
                }
 
                for (unsigned int i = 1; i < 256; ++i)
                {
                    s_sum[i] += s_sum[i - 1];
                    s0[i] += s_sum[i - 1];
                }
                unsigned char* b = (unsigned char*)&source[right_index].key + digit;
                TParticleCode* v1 = &source[right_index];
                TParticleCode* v2 = &source[left_index];
                while (v1 >= v2)
                {
                    dest[--s0[*b]] = *v1--;
                    b -= sizeof(TParticleCode);
                }
#pragma omp barrier
                std::swap(source, dest);
            }
        }
 
        // Если ключ структуры однобайтовый, просто копируем в исходный массив
        if (sizeof(m->key) == 1)
        {
            memcpy(m, m_temp, n * sizeof(TParticleCode));
        }
    }

Here is the caller graph for this function:

◆ RSort_step3()

void BH::MortonTree::RSort_step3	(	TParticleCode *	source,
		TParticleCode *	dest,
		unsigned int	n,
		unsigned int *	offset,
		unsigned char	sortable_bit
	)

inlineprivate

Вспомогательная функция к RSort_Node3.

Definition at line 618 of file TreeBH.cpp.

    {
        //unsigned char* b = (unsigned char*)&source[n].key + sortable_bit;
 
        for (unsigned int i = 0; i < n; ++i)
        {
            TParticleCode* src = &source[i];
            unsigned int off = (src->key >> (sortable_bit * 8)) & 0xFF;
            dest[offset[off]++] = *src;
        }
    }

Here is the caller graph for this function:

◆ SetCellGeometry()

void BH::MortonTree::SetCellGeometry ( int cell )

private

Функция вычисления геометрических параметров внутренней ячейки дерева

Definition at line 236 of file TreeBH.cpp.

    {
        int prLength;
        unsigned int pr;
        std::tie<unsigned int, int>(pr, prLength) = Prefix(cell);
        prLength -= PrefixLength(0);
 
        Point2D sz = { 1.0 / (double)(1 << ceilhalf(prLength)), 1.0 / (1 << (prLength / 2)) };
 
        Point2D pos = 0.5 * sz;
 
        for (int i = 0; i < prLength; i += 2)
            if (pr & (1 << (prLength - i - 1)))
                pos[0] += 1.0 / (1 << (1 + i / 2));
 
        for (int i = 1; i < prLength; i += 2)
            if (pr & (1 << (prLength - i - 1)))
                pos[1] += 1.0 / (1 << (1 + i / 2));
 
        mortonTree[cell].center = pos * (1.0 / iQuadSideVar) + lowLeft;
        mortonTree[cell].size = sz * (1.0 / iQuadSideVar);
        mortonTree[cell].level = prLength;
 
        mortonTree[cell].leaf = (prLength >= maxTreeLevel);
    }//SetCellGeometry(...)

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ShareGabs()

void BH::MortonTree::ShareGabs ( int numCell )

Definition at line 854 of file TreeBH.cpp.

    {
        auto& lcell = mortonTree[mortonTree[numCell].child[0]]; //левый ребенок     
        auto& rcell = mortonTree[mortonTree[numCell].child[1]]; //правый ребенок 
        auto& pcell = mortonTree[numCell];
        {
            pcell.minx_miny[0] = std::min(lcell.minx_miny[0], rcell.minx_miny[0]);
            pcell.minx_miny[1] = std::min(lcell.minx_miny[1], rcell.minx_miny[1]);
            pcell.maxx_maxy[0] = std::max(lcell.maxx_maxy[0], rcell.maxx_maxy[0]);
            pcell.maxx_maxy[1] = std::max(lcell.maxx_maxy[1], rcell.maxx_maxy[1]);
        }
 
        pcell.center = (pcell.minx_miny + pcell.maxx_maxy) * 0.5;
        pcell.size = pcell.maxx_maxy - pcell.minx_miny;
    }

Here is the caller graph for this function:

Member Data Documentation

◆ binomCft

std::vector<double> BH::MortonTree::binomCft

список биномиальных коэффициентов (имеет смысл только нижняя треугольная часть матрицы, остальные = 0)

Definition at line 223 of file TreeBH.h.

◆ iFact

std::vector<double> BH::MortonTree::iFact

Спикок обратных факториалов целых чисел

Definition at line 219 of file TreeBH.h.

◆ iQuadSideVar

double BH::MortonTree::iQuadSideVar

private

Масштабный фактор дерева (на сколько нужно умножить сторону квадрата содержащего все частицы, чтобы новый квадрат имел размер [0, U]x[0, U], где U выбирается так, чтобы все мортоновские коды заданной длины были в диапазоне [0, 1), т.е. положение частицы должно быть в диапазоне [0, 1 / (1.0 + 1/(2^codeLength - 1) ))

Definition at line 165 of file TreeBH.h.

◆ lowLeft

Point2D BH::MortonTree::lowLeft

private

Положение нижнего левого угла дерева

Definition at line 168 of file TreeBH.h.

◆ maxTreeLevel

const int BH::MortonTree::maxTreeLevel

Максимальная глубина дерева при его обходе

Definition at line 226 of file TreeBH.h.

◆ mortonCodes

std::vector<TParticleCode> BH::MortonTree::mortonCodes

Список мортоновских кодов, вычисляемых по положениям частиц

Definition at line 157 of file TreeBH.h.

◆ mortonLowCells

std::vector<int> BH::MortonTree::mortonLowCells

Вектор индесков ячеек нижнего уровня в дереве

Definition at line 247 of file TreeBH.h.

◆ mortonTree

std::vector<treeCellT> BH::MortonTree::mortonTree

Мортоновское дерево Размер дерева равен сумме числа внутренних вершин (число частиц без единицы) и числа частиц резервируется на 1 штуку больше для удобства работы ячейки [0 ... pos.size()-2] отвечают внутренним узлам, ячейка с идексом pos.size()-1 не используется ячейки [pos.size() ... 2*pos.size()-1] отвечают частицам

Definition at line 257 of file TreeBH.h.

◆ pointsCopy

std::vector<PointsCopy>& BH::MortonTree::pointsCopy

private

Ссылка на список обобщенных вихрей, по которым строится дерево

Definition at line 153 of file TreeBH.h.

◆ prm

const params& BH::MortonTree::prm

private

Ссылка на параметры, загружаемые из файла

Definition at line 150 of file TreeBH.h.

The documentation for this class was generated from the following files:

VM2D/Velocity2D/TreeBH.h
VM2D/Velocity2D/TreeBH.cpp

Public Member Functions

Public Attributes

Private Member Functions

Static Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MortonTree()

Member Function Documentation

◆ AddVelo()

◆ BuildInternalTreeCell()

◆ BuildMortonInternalTree()

◆ BuildMortonParticlesTree()

◆ CalcLocalCoeffToLowLevel()

◆ CalculateMortonTreeParams()

◆ CalcVeloBiotSavart()

◆ CalcVeloTaylorExpansion()

◆ Delta()

◆ ExpandBits()

◆ FillMortonLowCells()

◆ FillMortonLowCellsA()

◆ FindEnclosingRectangleOld()

◆ GabsPyr()

◆ getBinom()

◆ getStatistics()

◆ MakeRootMortonTree()

◆ MinMaxCellsCalc()

◆ MinMaxCellsSearch()

◆ Morton2D()

◆ Prefix()

◆ PrefixLength()

◆ RSort_Node3()

◆ RSort_Parallel()

◆ RSort_step3()

◆ SetCellGeometry()

◆ ShareGabs()

Member Data Documentation

◆ binomCft

◆ iFact

◆ iQuadSideVar

◆ lowLeft

◆ maxTreeLevel

◆ mortonCodes

◆ mortonLowCells

◆ mortonTree

◆ pointsCopy

◆ prm