Commit 9619194c authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Finalize the construction of the Let. Add periodicity in the box and remove...

Finalize the construction of the Let. Add periodicity in the box and remove top_level in the tree constructor
parent ba3b4f25
......@@ -148,6 +148,7 @@ auto run(parallel_manager& para, const std::string& input_file, const int tree_h
use_particle_distribution);
#ifdef SCALFMM_USE_MPI
std::cout << std::flush;
para.get_communicator().barrier();
#endif
......
// See LICENCE file at project root
// See LICENCE file at project root
// File : box.hpp
// --------------------------------
#ifndef SCALFMM_TREE_BOX_HPP
......@@ -39,29 +39,31 @@ namespace scalfmm::component
using space_filling_curve_t = SpaceFillingCurve<dimension>;
private:
position_type _c1; ///< Minimum corner
position_type _c2; ///< Maximum corner
position_type _center; ///< Center
space_filling_curve_t _space_filling_curve;
position_type m_c1; ///< Minimum corner
position_type m_c2; ///< Maximum corner
position_type m_center; ///< Center
space_filling_curve_t m_space_filling_curve;
std::array<bool, dimension> m_periodicity; ///< get the periodicity per direction
bool m_is_periodic; ///< True if a direction is periodic
/** Rearranges the corners to ensure the maximum-minimum predicate */
void rearrange_corners()
{
for(std::size_t i = 0; i < dimension; ++i)
{
if(_c1[i] > _c2[i])
if(m_c1[i] > m_c2[i])
{
std::swap(_c1[i], _c2[i]);
std::swap(m_c1[i], m_c2[i]);
}
}
_center = (_c1 + _c2) / 2;
m_center = (m_c1 + m_c2) / 2;
}
public:
/// Accessor for the minimum corner
[[nodiscard]] auto c1() const noexcept -> position_type const& { return _c1; }
[[nodiscard]] auto c1() const noexcept -> position_type const& { return m_c1; }
/// Accessor for the maximum corner
[[nodiscard]] auto c2() const noexcept -> position_type const& { return _c2; }
[[nodiscard]] auto c2() const noexcept -> position_type const& { return m_c2; }
/** Builds an empty box at the origin */
box() = default;
......@@ -87,21 +89,25 @@ namespace scalfmm::component
* \param side_length The cube's side length
**/
[[deprecated]] box(const position_type& min_corner, value_type side_length)
: _c1(min_corner)
, _c2(min_corner)
, _space_filling_curve()
: m_c1(min_corner)
, m_c2(min_corner)
, m_space_filling_curve()
, m_is_periodic(false)
{
if(side_length < 0)
{
side_length = -side_length;
}
for(auto&& d: _c2)
for(auto& v: m_periodicity)
{
v = false;
}
for(auto&& d: m_c2)
{
d += side_length;
}
_center = (_c1 + _c2) / 2;
m_center = (m_c1 + m_c2) / 2;
}
/** Builds a cube using its center and width
......@@ -111,23 +117,28 @@ namespace scalfmm::component
* \param unused Disambiguisation param from lower corner constructor
*/
box(value_type width, const position_type& box_center)
: _c1(box_center)
, _c2(box_center)
, _center(box_center)
, _space_filling_curve()
: m_c1(box_center)
, m_c2(box_center)
, m_center(box_center)
, m_space_filling_curve()
, m_is_periodic(false)
{
if(width < 0)
{
width = -width;
}
for(auto& v: m_periodicity)
{
v = false;
}
value_type radius = width / 2;
for(auto&& d: _c1)
for(auto&& d: m_c1)
{
d -= radius;
}
for(auto&& d: _c2)
for(auto&& d: m_c2)
{
d += radius;
}
......@@ -141,10 +152,16 @@ namespace scalfmm::component
* \param corner_2 The second corner.
*/
box(const position_type& corner_1, const position_type& corner_2)
: _c1(corner_1)
, _c2(corner_2)
, _space_filling_curve()
: m_c1(corner_1)
, m_c2(corner_2)
, m_space_filling_curve()
, m_is_periodic(false)
{
for(auto& v: m_periodicity)
{
v = false;
}
rearrange_corners();
}
......@@ -157,8 +174,8 @@ namespace scalfmm::component
*/
void set(const position_type& corner_1, const position_type& corner_2)
{
_c1 = corner_1;
_c2 = corner_2;
m_c1 = corner_1;
m_c2 = corner_2;
rearrange_corners();
}
......@@ -171,7 +188,7 @@ namespace scalfmm::component
{
for(std::size_t i = 0; i < dimension; i++)
{
if(p[i] < _c1[i] || p[i] > _c2[i])
if(p[i] < m_c1[i] || p[i] > m_c2[i])
{
return false;
}
......@@ -193,7 +210,7 @@ namespace scalfmm::component
}
/** Accessor for the box center */
[[nodiscard]] auto center() const -> position_type const& { return _center; }
[[nodiscard]] auto center() const -> position_type const& { return m_center; }
/** Accessor for the box corners
*
......@@ -206,9 +223,9 @@ namespace scalfmm::component
{
position_type c;
std::size_t i = 0;
for(bool choice: _space_filling_curve.position(idx))
for(bool choice: m_space_filling_curve.position(idx))
{
c[i] = choice ? _c2[i] : _c1[i];
c[i] = choice ? m_c2[i] : m_c1[i];
++i;
}
return c;
......@@ -225,29 +242,54 @@ namespace scalfmm::component
void corner(std::size_t idx, const position_type& pos)
{
std::size_t i = 0;
for(bool choice: _space_filling_curve.position(idx))
for(bool choice: m_space_filling_curve.position(idx))
{
if(choice)
{
_c2[i] = pos[i];
m_c2[i] = pos[i];
}
else
{
_c1[i] = pos[i];
m_c1[i] = pos[i];
}
++i;
}
rearrange_corners();
}
///
/// \brief set the periodicity in the direction dir
/// \param dir direction
/// \param per true if the direction dir is periodique otherwise false
///
void set_periodicity(int dir, bool per)
{
m_periodicity[dir] = per;
m_is_periodic = m_is_periodic && per;
}
///
/// \brief set the periodicity in the direction dir
/// \param pbl a vector of boolean that specifies if teh direction is periodic or not
///
void set_periodicity(const std::array<bool, dimension> pbl)
{
for(std::size_t d = 0; d < pbl.size(); ++d)
set_periodicity(d, pbl[d]);
}
///
/// \brief is_periodic tell if there is a periodic direction
///
///
auto is_periodic() const noexcept -> bool { return m_is_periodic; }
/** Returns the width for givmusique traditionnelle ameren dimension */
[[nodiscard]] auto width(std::size_t dim) const noexcept -> decltype(std::abs(_c2[dim] - _c1[dim]))
/** Returns the width for given dimension */
[[nodiscard]] auto width(std::size_t dim) const noexcept -> decltype(std::abs(m_c2[dim] - m_c1[dim]))
{
return std::abs(_c2[dim] - _c1[dim]);
return std::abs(m_c2[dim] - m_c1[dim]);
}
/** Sums the corners of two boxes */
auto operator+(const box& other) const -> box { return box(_c1 + other._c1, _c2 + other._c2); }
auto operator+(const box& other) const -> box { return box(m_c1 + other.m_c1, m_c2 + other.m_c2); }
/** Tests two boxes equality */
auto operator==(const box& other) const -> bool { return c1() == other.c1() && c2() == other.c2(); }
......
......@@ -30,26 +30,19 @@ namespace scalfmm::component
template<typename ParticleContainer>
explicit dist_group_tree(std::size_t tree_height, std::size_t order, Box const& box,
std::size_t size_leaf_blocking, std::size_t size_cell_blocking,
ParticleContainer const& particle_container, bool particles_are_sorted = false,
int in_top_level = 2 /*, int in_left_limit = -1*/)
ParticleContainer const& particle_container,
bool particles_are_sorted = false /*, int in_left_limit = -1*/)
: base_type(tree_height, order, box, size_leaf_blocking, size_cell_blocking, particle_container,
particles_are_sorted, in_top_level /*, in_left_limit*/)
particles_are_sorted /*, in_left_limit*/)
{
m_cell_distrib.resize(tree_height);
}
void set_leaf_distribution(const data_distrib_type& in_leaf_distrib)
{
m_leaf_distrib.resize(in_leaf_distrib.size());
std::copy(in_leaf_distrib.begin(), in_leaf_distrib.end(), m_leaf_distrib.begin());
// m_leaf_distrib = in_leaf_distrib;
}
void set_leaf_distribution(const data_distrib_type& in_leaf_distrib) { m_leaf_distrib = in_leaf_distrib; }
void set_cell_distribution(const int in_level, const data_distrib_type& in_cell_distrib)
{
m_cell_distrib[in_level].resize(in_cell_distrib.size());
std::copy(in_cell_distrib.begin(), in_cell_distrib.end(), m_cell_distrib[in_level].begin());
// m_cell_distrib[in_level] = in_cell_distrib;
m_cell_distrib[in_level] = in_cell_distrib;
}
void print_distrib(std::ostream& out, bool verbose = true)
......@@ -58,10 +51,8 @@ namespace scalfmm::component
if(verbose)
{
std::cout << "Tree distribution" << std::endl;
header = "leaf distribution: \n";
}
out::print(out, std::move(header), m_leaf_distrib);
for(int l = base_type::leaf_level(); l >= base_type::top_level(); --l)
for(int l = base_type::top_level(); l < base_type::height(); ++l)
{
if(verbose)
{
......@@ -69,7 +60,13 @@ namespace scalfmm::component
}
out::print(out, std::move(header), m_cell_distrib[l]);
}
if(verbose)
{
header = "leaf distribution: \n";
}
out::print(out, std::move(header), m_leaf_distrib);
}
~dist_group_tree()
{
std::cout << scalfmm::colors::red;
......
......@@ -830,8 +830,10 @@ namespace scalfmm::tree
}
if(need_comm)
{
std::cout << scalfmm::colors::red << "Need to adjust the distribution (not checked)"
<< scalfmm::colors::reset << '\n';
std::cout << scalfmm::colors::red << "Need to adjust the distribution (not checked)\n";
out::print("rank(" + std::to_string(rank) + ") parent_distrib: ", parent_distrib);
std::cout << scalfmm::colors::reset << '\n';
/// Find the good index
///
using element_distrib_type = typename MortonDistribution::value_type;
......@@ -871,13 +873,11 @@ namespace scalfmm::tree
{
mortonCellIndex[i] = mortonCellIndex[i] >> dimension;
}
out::print("rank(" + std::to_string(rank) + ") mortonCellIndex: ", mortonCellIndex);
/// we have to remove some elements at the begining if start != 0
auto last = std::unique(mortonCellIndex.begin(), mortonCellIndex.end());
if(start > 0)
{
VectorMortonIdx new_cell_index;
VectorMortonIdx new_cell_index(std::distance(mortonCellIndex.begin() + start, last));
std::move(mortonCellIndex.begin() + start, last, new_cell_index.begin());
mortonCellIndex = std::move(new_cell_index);
}
......@@ -885,7 +885,7 @@ namespace scalfmm::tree
{
mortonCellIndex.erase(last, mortonCellIndex.end());
}
out::print("rank(" + std::to_string(rank) + ") mortonCellIndex: ", mortonCellIndex);
out::print("rank(" + std::to_string(rank) + ") mortonCellIndex1: ", mortonCellIndex);
std::cout << scalfmm::colors::blue << " --> End distrib::build_upper_distribution at level " << level
<< scalfmm::colors::reset << std::endl;
......@@ -1558,28 +1558,13 @@ namespace scalfmm::tree
////////////////////////////////////////////////////////////////////////////////////////////
///
/// Construct the local tree based on our set of particles
/// and compute the new morton indexes associated to the particles on the process
///
/// Put the particles inside the container or particles (for compatibility)
///
using container_type = scalfmm::container::particle_container<typename Vector_type::value_type>;
container_type container(particle_container.size());
leafMortonIdx.resize(particle_container.size());
#pragma omp parallel for shared(localNumberOfParticles, box, leaf_level)
for(std::size_t i = 0; i < particle_container.size(); ++i)
{
leafMortonIdx[i] = scalfmm::index::get_morton_index(particle_container[i].position(), box, leaf_level);
container.push_particle(i, particle_container[i]);
}
/// Set true because the particles are already sorted
/// In fact we have all the leaves to add in leafMortonIdx - could be used to construct
/// the tree !!!
///
// localGroupTree = new Tree_type(static_cast<std::size_t>(leaf_level + 1), order, box,
// static_cast<std::size_t>(groupSizeCells), container, true);
localGroupTree = new Tree_type(static_cast<std::size_t>(leaf_level + 1), order, box, groupSizeLeaves,
groupSizeCells, /*particle_*/ container, true, 2);
groupSizeCells, particle_container, true);
#ifdef SCALFMM_USE_MPI
// std::cout << scalfmm::colors::red;
......@@ -1588,9 +1573,17 @@ namespace scalfmm::tree
/// End
////////////////////////////////////////////////////////////////////////////////////////////
///
/// Compute the new morton indexes associated to the particles on the process
///
leafMortonIdx.resize(particle_container.size());
#pragma omp parallel for shared(localNumberOfParticles, box, leaf_level)
for(std::size_t i = 0; i < particle_container.size(); ++i)
{
leafMortonIdx[i] = scalfmm::index::get_morton_index(particle_container[i].position(), box, leaf_level);
}
auto last = std::unique(leafMortonIdx.begin(), leafMortonIdx.end());
leafMortonIdx.erase(last, leafMortonIdx.end());
///
if(manager.get_num_processes() > 1)
{
////////////////////////////////////////////////////////////////////////////////////////////
......@@ -1635,9 +1628,11 @@ namespace scalfmm::tree
{
level_dist[l] = std::move(
distrib::build_upper_distribution(manager, dimension, l, leafMortonIdx, level_dist[l + 1]));
// out::print("rank(" + std::to_string(rank) + ") leafMortonIdx: ", leafMortonIdx);
build_let_at_level(manager, *localGroupTree, leafMortonIdx, level_dist[l], l);
// out::print("rank(" + std::to_string(rank) + ") leafMortonIdx: ",
// leafMortonIdx); out::print("xxrank(" + std::to_string(rank) + ") part: ",
// level_dist[l]);
localGroupTree->set_cell_distribution(l, level_dist[l]);
build_let_at_level(manager, *localGroupTree, leafMortonIdx, level_dist[l], l);
}
}
else
......@@ -1645,21 +1640,14 @@ namespace scalfmm::tree
{
localGroupTree->set_leaf_distribution(particles_distrib);
localGroupTree->set_cell_distribution(leaf_level, leaves_distrib);
for(int l = leaf_level - 1; l >= localGroupTree->top_level(); --l)
{
morton_distrib_type local{{leaves_distrib[0][0] >> dimension, leaves_distrib[0][1] >> dimension}};
localGroupTree->set_cell_distribution(l, local);
leaves_distrib[0][0] = leaves_distrib[0][0] >> dimension;
leaves_distrib[0][1] = leaves_distrib[0][1] >> dimension;
localGroupTree->set_cell_distribution(l, leaves_distrib);
}
}
//#else
// localGroupTree->set_leaf_distribution(particles_distrib);
// localGroupTree->set_cell_distribution(leaf_level, leaves_distrib);
// for(int l = leaf_level - 1; l >= localGroupTree->top_level(); --l)
// {
// morton_distrib_type local{{leaves_distrib[0][0] >> dimension, leaves_distrib[0][1] >>
// dimension}}; localGroupTree->set_cell_distribution(l, local);
// }
//#endif
std::cout << scalfmm::colors::green << " --> End let::group_let() " << scalfmm::colors::reset << std::endl;
}
......
......@@ -17,6 +17,7 @@
#include <utility>
#include <vector>
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/operators/generic/tags.hpp"
#include "scalfmm/tree/box.hpp"
#include "scalfmm/tree/for_each.hpp"
......@@ -24,7 +25,6 @@
#include "scalfmm/tree/interaction_list.hpp"
#include "scalfmm/tree/utils.hpp"
#include "scalfmm/utils/sort.hpp"
#include "scalfmm/meta/utils.hpp"
namespace scalfmm::component
{
......@@ -166,10 +166,11 @@ namespace scalfmm::component
/// @param begin_container
///
/// @return
//template<typename ContainerIterator>
//auto create_leaf(std::vector<std::tuple<std::size_t, std::size_t>> const& tuple_of_indexes,
// std::size_t morton_index_of_leaf, std::size_t number_of_particles_in_leaf, box_type const& box,
// std::size_t offset_particle, std::size_t component_id, ContainerIterator begin_container)
// template<typename ContainerIterator>
// auto create_leaf(std::vector<std::tuple<std::size_t, std::size_t>> const& tuple_of_indexes,
// std::size_t morton_index_of_leaf, std::size_t number_of_particles_in_leaf, box_type const&
// box, std::size_t offset_particle, std::size_t component_id, ContainerIterator
// begin_container)
// -> leaf_type
//{
// // get the coordinate of the leaf in the tree
......@@ -293,15 +294,14 @@ namespace scalfmm::component
{
// start and end of the leaves
auto start_index{g * m_number_of_leaves_per_group};
auto end_index{(g * m_number_of_leaves_per_group) + m_number_of_leaves_per_group-1};
auto end_index{(g * m_number_of_leaves_per_group) + m_number_of_leaves_per_group - 1};
// get the starting and ending morton indices i.e the first and last leaves morton indices
auto starting_morton_index{vector_of_mortons.at(start_index)};
auto ending_morton_index{vector_of_mortons.at(end_index)};
// we create the group and stores the pointer in the vector of groups
m_group_of_leaf.at(g) = std::move(
std::make_shared<group_of_leaf_type>(starting_morton_index, ending_morton_index+1,
m_number_of_leaves_per_group, g));
m_group_of_leaf.at(g) = std::move(std::make_shared<group_of_leaf_type>(
starting_morton_index, ending_morton_index + 1, m_number_of_leaves_per_group, g));
}
// here is the residue of the leaves
......@@ -310,16 +310,16 @@ namespace scalfmm::component
// we go at the end of the full groups
auto start_index{number_of_groups_at_leaf_level * m_number_of_leaves_per_group};
auto end_index{(number_of_groups_at_leaf_level * m_number_of_leaves_per_group) +
remain_number_of_leaves-1};
remain_number_of_leaves - 1};
// we get the morton indices
auto starting_morton_index{vector_of_mortons.at(start_index)};
auto ending_morton_index{vector_of_mortons.at(end_index)};
// we create the last group and store the pointer at the end of the vector of groups
m_group_of_leaf.push_back(std::move(std::make_shared<group_of_leaf_type>(
starting_morton_index, ending_morton_index+1, remain_number_of_leaves,
number_of_groups_at_leaf_level + 1)));
m_group_of_leaf.push_back(std::move(
std::make_shared<group_of_leaf_type>(starting_morton_index, ending_morton_index + 1,
remain_number_of_leaves, number_of_groups_at_leaf_level + 1)));
}
}
......@@ -331,8 +331,8 @@ namespace scalfmm::component
///
/// @return void
template<typename MortonIndex>
auto build_groups_of_cells_at_level(std::vector<MortonIndex> const& vector_of_mortons,
std::size_t level) -> void
auto build_groups_of_cells_at_level(std::vector<MortonIndex> const& vector_of_mortons, std::size_t level)
-> void
{
auto number_of_cells{vector_of_mortons.size()};
auto number_of_groups{number_of_cells / m_number_of_cells_per_group};
......@@ -346,28 +346,27 @@ namespace scalfmm::component
{
// start and end of the cells
auto start_index{g * m_number_of_cells_per_group};
auto end_index{(g * m_number_of_cells_per_group) + m_number_of_cells_per_group-1};
auto end_index{(g * m_number_of_cells_per_group) + m_number_of_cells_per_group - 1};
// get the starting and ending morton indices i.e the first and last cells morton indices
auto starting_morton_index{vector_of_mortons.at(start_index)};
auto ending_morton_index{vector_of_mortons.at(end_index)};
// Create a group
m_group_of_cell_per_level.at(level).at(g) = std::move(std::make_shared<group_of_cell_type>(
starting_morton_index, ending_morton_index+1, m_number_of_cells_per_group));
starting_morton_index, ending_morton_index + 1, m_number_of_cells_per_group));
}
if(remain_number_of_cells != 0)
{
// we go at the end of the full groups
auto start_index{number_of_groups * m_number_of_cells_per_group};
auto end_index{(number_of_groups * m_number_of_cells_per_group) +
remain_number_of_cells-1};
auto end_index{(number_of_groups * m_number_of_cells_per_group) + remain_number_of_cells - 1};
// we get the morton indices
auto starting_morton_index{vector_of_mortons.at(start_index)};
auto ending_morton_index{vector_of_mortons.at(end_index)};
m_group_of_cell_per_level.at(level).push_back(std::move(std::make_shared<group_of_cell_type>(
starting_morton_index, ending_morton_index+1, remain_number_of_cells)));
starting_morton_index, ending_morton_index + 1, remain_number_of_cells)));
}
}
......@@ -386,7 +385,7 @@ namespace scalfmm::component
{
std::size_t leaf_index{0};
// loop on groups
for(std::size_t ng{0}; ng<m_group_of_leaf.size(); ++ng)
for(std::size_t ng{0}; ng < m_group_of_leaf.size(); ++ng)
{
auto pg = m_group_of_leaf.at(ng);
auto& group_symbolics{pg->symbolics()};
......@@ -401,7 +400,7 @@ namespace scalfmm::component
index::get_box(box.c1(), box.width(0), coordinate, m_tree_height, m_tree_height - 1)};
leaf = std::move(leaf_type(number_of_particles_per_leaves.at(leaf_index), std::get<1>(width_center),
std::get<0>(width_center), morton_index_of_leaf));
std::get<0>(width_center), morton_index_of_leaf));
// accumulate for set the number of particle in group.
group_symbolics.number_of_particles_in_group += number_of_particles_per_leaves.at(leaf_index);
++leaf_index;
......@@ -418,26 +417,25 @@ namespace scalfmm::component
///
/// @return void
template<typename MortonType>
auto build_cells_in_groups_at_level(std::vector<MortonType> const& vector_of_mortons,
box_type const& box, std::size_t level)
-> void
auto build_cells_in_groups_at_level(std::vector<MortonType> const& vector_of_mortons, box_type const& box,
std::size_t level) -> void
{
std::size_t cell_index{0};
// loop on groups
for(std::size_t ng{0}; ng<m_group_of_cell_per_level.at(level).size(); ++ng)
for(std::size_t ng{0}; ng < m_group_of_cell_per_level.at(level).size(); ++ng)
{
auto pg = m_group_of_cell_per_level.at(level).at(ng);
// loop on leaves
for(auto&& cell : pg->block())
for(auto&& cell: pg->block())
{
auto morton_index_of_cell = vector_of_mortons.at(cell_index);
// get the coordinate of the leaf in the tree
auto coordinate = index::get_coordinate_from_morton_index<dimension>(morton_index_of_cell);
// get the corresponding box to the leaf
auto width_center{
index::get_box(box.c1(), box.width(0), coordinate, m_tree_height, level)};
auto width_center{index::get_box(box.c1(), box.width(0), coordinate, m_tree_height, level)};
cell = std::move(cell_type(std::get<1>(width_center), std::get<0>(width_center), m_order, level, morton_index_of_cell, coordinate));
cell = std::move(cell_type(std::get<1>(width_center), std::get<0>(width_center), m_order, level,
morton_index_of_cell, coordinate));
++cell_index;
}
}
......@@ -493,7 +491,6 @@ namespace scalfmm::component