fmm::MortonTree< dim, point, value > Class Template Reference

#include <morton_tree.h>

Collaboration diagram for fmm::MortonTree< dim, point, value >:

Public Types
using	particle_t = particle< point, value >
using	TreeCell_t = TreeCell< point, value >

Public Member Functions
	MortonTree (std::vector< particle_t > &&particles, size_t tree_depth)

Public Attributes
std::vector< std::vector< TreeCell_t > >	levels
std::vector< size_t >	level_sizes
size_t	tree_depth
size_t	targets_num = 0
std::vector< particle_t >	particles
std::vector< size_t >	positions_map

Private Member Functions
void	FillNeighbours ()

Private Attributes
point	shift
double	scale

Detailed Description

template<int dim, typename point, typename value>
class fmm::MortonTree< dim, point, value >

Definition at line 25 of file morton_tree.h.

Member Typedef Documentation

particle_t

template<int dim, typename point , typename value >

using fmm::MortonTree< dim, point, value >::particle_t = particle<point, value>

Definition at line 28 of file morton_tree.h.

TreeCell_t

template<int dim, typename point , typename value >

using fmm::MortonTree< dim, point, value >::TreeCell_t = TreeCell<point, value>

Definition at line 29 of file morton_tree.h.

Constructor & Destructor Documentation

MortonTree()

template<int dim, typename point , typename value >

fmm::MortonTree< dim, point, value >::MortonTree	(	std::vector< particle_t > &&	particles,
		size_t	tree_depth
	)

Definition at line 287 of file morton_tree.cpp.

                                                                                                : tree_depth(tree_depth_), particles(particles_)
{
    double T = omp_get_wtime();
    //if (IAmRoot()) std::cout << "\n************ Start Building Tree ************" << std::endl;
    const int dim2 = MyPow(2, dim);
    double t = omp_get_wtime();
    {
        auto be = box_extent(particles);
        shift = be.first; scale = be.second;
 
        while (scale  < FORCE_EPS * (1 << tree_depth))
        {
            --tree_depth;
        }
 
        if constexpr (dim == 2)
            shift -= std::complex<double> {1.e-7, 1.e-7}; // min particles shouldnt have coordinates equal zero
        else
            shift -= Vector3d {1.e-7, 1.e-7, 1.e-7}; // min particles shouldnt have coordinates equal zero
        scale += 1.e-6; // max particles shouldnt have coordinate equal 1
    }
    //if (IAmRoot()) std::cout << "compute extent time: " << omp_get_wtime() - t << std::endl;
    //if (IAmRoot()) std::cout << "scale factors: shift = " << shift << ", extent = " << scale << std::endl;
 
    size_t num_particles = particles.size();
    // if (IAmRoot()) std::cout << "particles num = " << num_particles << std::endl;
 
    levels.resize(tree_depth);
    tbb::parallel_for(size_t(0), tree_depth, [&](size_t i) {
        levels[i].resize(std::min(num_particles, (size_t)MyPow(size_t(dim2), (unsigned int)i)));
        });
    level_sizes.resize(tree_depth);
    positions_map.resize(num_particles);
 
    t = omp_get_wtime();
    tbb::parallel_for_each(particles.begin(), particles.end(), [&](particle_t& x) {
        x.morton_code = MortonCode<dim, point>(x.center, (int)(tree_depth - 1), shift, scale);
        assert(x.morton_code < MyPow(dim2, tree_depth - 1));
        });
    //if (IAmRoot()) std::cout << "compute morton codes time: " << omp_get_wtime() - t << std::endl;
 
    assert(MyPow(dim2, tree_depth) > std::max_element(particles.begin(), particles.end(), [&](const particle_t& x, const particle_t& y) {return x.morton_code < y.morton_code; })->morton_code);
    t = omp_get_wtime();
    {
#ifdef deterministic_behaviour
        tbb::parallel_sort(particles.begin(), particles.end(), [](const auto& x, const auto& y) {return x.morton_code < y.morton_code; });
        //boost::sort::block_indirect_sort(particles.begin(), particles.end(), [](const auto& x, const auto& y) {return x.morton_code < y.morton_code; }, 1);
#else
        if (IAmRoot())
        {
            // counting sort, using known number of unique numbers
            tbb::concurrent_vector<std::atomic_size_t> counts(MyPow(dim2, (unsigned int)(tree_depth - 1)) + 1);
            tbb::parallel_for_each(particles.begin(), particles.end(), [&](const auto& x) {counts[1 + x.morton_code]++; });
            for (int i = 1; i < counts.size(); ++i)
                counts[i] += counts[i - 1];
            assert(counts.back() == num_particles);
 
            tbb::scalable_allocator<small_particle<dim, value>> temp;
            constexpr size_t memsize = sizeof(small_particle<dim, value>);
            auto ptr = temp.allocate(num_particles);
 
            //t = omp_get_wtime();
            tbb::parallel_for(size_t(0), num_particles, [&](size_t i) {
                const auto& x = particles[i];
                size_t pos = counts[x.morton_code]++;
                assert(pos < num_particles);
                positions_map[i] = pos;
                memcpy(&ptr[pos], &x, memsize);
                });
 
            tbb::parallel_for(size_t(0), num_particles, [&](size_t i) {
                memcpy(&particles[i], &ptr[i], memsize);
                });
 
            temp.deallocate(ptr, num_particles);
 
            assert(std::is_sorted(particles.begin(), particles.end(), [](const auto& x, const auto& y) {return x.morton_code < y.morton_code; }));
        }
#ifdef FMM_MPI
        MPI_Bcast(particles.data(), num_particles * sizeof(particle_t), MPI_BYTE, RootID(), MPI_COMM_WORLD);
        MPI_Bcast(positions_map.data(), num_particles * sizeof(size_t), MPI_BYTE, RootID(), MPI_COMM_WORLD);
#endif
 
#endif
    }
    //if (IAmRoot()) std::cout << "sort time: " << omp_get_wtime() - t << std::endl;
    t = omp_get_wtime();
#ifdef deterministic_behaviour
    tbb::parallel_for_each(particles.begin(), particles.end(), [&](auto& x) {x.morton_code >>= dim * (30 - tree_depth); });
#endif
    {
        auto& top_level = levels.back();
        auto ptr = std::unique_copy(exec_policy, particles.begin(), particles.end(), top_level.begin(), [&](const auto& x, const auto& y) {return (x.morton_code) == (y.morton_code); });
        level_sizes.back() = ptr - top_level.begin();
 
        int idx_bottom = 0, idx_top = 0;
        for (int j = 0; j < num_particles; ++j)
        {
            if ((particles[j].morton_code) != (top_level[idx_top].morton_code))
            {
                auto& cell = top_level[idx_top];
                cell.source_range = { idx_bottom, j };
                cell.center = PointFromMortonCode(cell.morton_code, (int)(tree_depth - 1), shift, scale);
                ++idx_top;
                idx_bottom = j;
            }
        }
        auto& cell = top_level[idx_top];
        cell.source_range = { idx_bottom, num_particles };
        cell.center = PointFromMortonCode(cell.morton_code, (int)(tree_depth - 1), shift, scale);
 
        //if (IAmRoot()) std::cout << "level " << tree_depth - 1 << " size: " <<  level_sizes.back() << std::endl;
    }
 
    for (int i = (int)(tree_depth - 2); i >= 0; --i)
    {
        auto& top_level = levels[i];
        auto& bottom_level = levels[i + 1];
 
        auto ptr = std::unique_copy(exec_policy, bottom_level.begin(), bottom_level.begin() + level_sizes[i + 1], top_level.begin(), [&](const auto& x, const auto& y) {return (x.morton_code >> dim) == (y.morton_code >> dim); });
        level_sizes[i] = ptr - top_level.begin();
 
        //if (IAmRoot()) std::cout << "level " << i << " size: " << level_sizes[i] << std::endl;
 
        int idx_bottom = 0, idx_top = 0;
        for (int j = 0; j < level_sizes[i + 1]; ++j)
        {
            if ((bottom_level[j].morton_code >> dim) != (top_level[idx_top].morton_code >> dim))
            {
                auto& cell = top_level[idx_top];
                cell.source_range = { idx_bottom, j };
                cell.morton_code >>= dim;
                cell.center = PointFromMortonCode(cell.morton_code, i, shift, scale);
                ++idx_top;
                idx_bottom = j;
            }
            bottom_level[j].parent = idx_top;
        }
        auto& cell = top_level[idx_top];
        cell.source_range = { idx_bottom, level_sizes[i + 1] };
        cell.morton_code >>= dim;
        cell.center = PointFromMortonCode(cell.morton_code, i, shift, scale);
    }
    //if (IAmRoot()) std::cout << "create levels time: " << omp_get_wtime() - t << std::endl;
 
    t = omp_get_wtime();
    {
        std::vector<size_t> inverse_positions_map(num_particles);
        tbb::parallel_for(size_t(0), num_particles, [&](size_t i){
            inverse_positions_map[positions_map[i]] = i;
        });
 
        value Zero;
        if constexpr (std::is_same_v<value, double>)
            Zero = 0;
        if constexpr (std::is_same_v<value, Vector3d>)
            Zero = { 0,0,0 };
        auto& leaves = levels.back();
        tbb::enumerable_thread_specific<size_t> local_num(0);
        tbb::parallel_for(size_t(0), level_sizes.back(), [&](size_t k){
            auto& leaf = leaves[k];
            auto& [s_begin, s_end] = leaf.source_range;
            auto& [t_begin, t_end] = leaf.target_range;
            t_end = s_end;
            t_begin = s_begin;
            for (size_t i = s_begin; i < s_end; ++i)
            {
                if (particles[i].q != Zero)
                {
                    std::swap(particles[t_begin], particles[i]);
                    std::swap(inverse_positions_map[t_begin++], inverse_positions_map[i]);
                }
            }
            s_end = t_begin;
            local_num.local() += t_end - t_begin;
 
        });
 
        tbb::parallel_for(size_t(0), num_particles, [&](size_t i){
                positions_map[inverse_positions_map[i]] = i;
            });
 
        for (const auto& x : local_num)
            targets_num += x;
        //if (IAmRoot()) std::cout << "targets num: " << targets_num << std::endl;
 
    }
    //if (IAmRoot()) std::cout << "source/target partition time: " << omp_get_wtime() - t << std::endl;
 
    t = omp_get_wtime();
    FillNeighbours();
    //if (IAmRoot()) std::cout << "fill neighbours time: " << omp_get_wtime() - t << std::endl;
 
    //������ ������ ������� ������������
    /*
    t = omp_get_wtime();
    if (targets_num == 0)
    {
        auto& leaves = levels.back();
        for (int i = 0; i < level_sizes.back(); ++i)
        {
            const auto& cell = leaves[i];
            for (const auto& idx : cell.closeneighbours)
            {
                auto& neighbour_cell = leaves[idx];
                auto ptr = std::find(neighbour_cell.closeneighbours.begin(), neighbour_cell.closeneighbours.end(), i);
                neighbour_cell.closeneighbours.erase(ptr);
            }
        }
    }
    if (IAmRoot()) std::cout << "remove duplicate nbr indices time:" << omp_get_wtime() - t << std::endl;
    */
 
    //if (IAmRoot()) std::cout << "total build time: " << omp_get_wtime() - T << std::endl;
    //if (IAmRoot()) std::cout << "************ End Building Tree ************\n\n";
}

Here is the call graph for this function:

Member Function Documentation

FillNeighbours()

template<int dim, typename point , typename value >

void fmm::MortonTree< dim, point, value >::FillNeighbours ( )

private

Definition at line 148 of file morton_tree.cpp.

{
    if (tree_depth > 1)
    for (int i = 0; i < level_sizes[1]; ++i)
    {
        auto& cell = levels[1][i];
        auto& [s, e] = levels[0][0].source_range;
        std::vector<size_t> indexes(e - s);
        std::iota(indexes.begin(), indexes.end(), s);
        indexes.erase(indexes.begin() + i);
        cell.closeneighbours = indexes;
    }
    for (int i = 2; i < tree_depth; ++i)
    {
        double dw = sqrt(dim + 1.e-5) * scale / (1 << i);
        dw *= dw;
        const auto& parent_level = levels[i - 1];
        auto& current_level = levels[i];
 
        tbb::parallel_for_each(current_level.begin(), current_level.begin() + level_sizes[i], [&](auto& cell)
            {
                if constexpr (dim == 2)
                {
                    cell.closeneighbours.reserve(9);
                    cell.farneighbours.reserve(27);
                }
                else
                {
                    cell.closeneighbours.reserve(27);
                    cell.farneighbours.reserve(189);
                }
 
                const auto& parent_cell = parent_level[cell.parent];
                for (const auto& parent_neighbours : parent_cell.closeneighbours)
                {
                    const auto& [s, e] = parent_level[parent_neighbours].source_range;
                    for (auto k = s; k < e; ++k)
                    {
                        if (norm(current_level[k].center - cell.center) < dw)
                            cell.closeneighbours.push_back(k);
                        else
                            cell.farneighbours.push_back(k);
                    }
                }
 
                const auto& [s, e] = parent_level[cell.parent].source_range;
                for (auto k = s; k < e; ++k)
                {
                    if (&cell != &(current_level[k]))
                    {
                        if (norm(current_level[k].center - cell.center) < dw)
                            cell.closeneighbours.push_back(k);
                        else
                            cell.farneighbours.push_back(k);
                    }
                }
                cell.farneighbours.shrink_to_fit();
                cell.closeneighbours.shrink_to_fit();
            });
    }
}

Here is the call graph for this function:

Member Data Documentation

level_sizes

template<int dim, typename point , typename value >

std::vector<size_t> fmm::MortonTree< dim, point, value >::level_sizes

Definition at line 34 of file morton_tree.h.

levels

template<int dim, typename point , typename value >

std::vector<std::vector<TreeCell_t> > fmm::MortonTree< dim, point, value >::levels

Definition at line 33 of file morton_tree.h.

particles

template<int dim, typename point , typename value >

std::vector<particle_t> fmm::MortonTree< dim, point, value >::particles

Definition at line 37 of file morton_tree.h.

positions_map

template<int dim, typename point , typename value >

std::vector<size_t> fmm::MortonTree< dim, point, value >::positions_map

Definition at line 38 of file morton_tree.h.

scale

template<int dim, typename point , typename value >

double fmm::MortonTree< dim, point, value >::scale

private

Definition at line 43 of file morton_tree.h.

shift

template<int dim, typename point , typename value >

point fmm::MortonTree< dim, point, value >::shift

private

Definition at line 42 of file morton_tree.h.

targets_num

template<int dim, typename point , typename value >

size_t fmm::MortonTree< dim, point, value >::targets_num = 0

Definition at line 36 of file morton_tree.h.

tree_depth

template<int dim, typename point , typename value >

size_t fmm::MortonTree< dim, point, value >::tree_depth

Definition at line 35 of file morton_tree.h.

---

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

• VMfmm/cpu/morton_tree.h
• VMfmm/cpu/morton_tree.cpp