diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 00000000..fda74983 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,9 @@ +{ + "rust-analyzer.linkedProjects": [ + "./tree/Cargo.toml", + "./tree/Cargo.toml", + "./tree/Cargo.toml", + "./tree/Cargo.toml", + "./tree/Cargo.toml" + ] +} \ No newline at end of file diff --git a/traits/src/tree.rs b/traits/src/tree.rs index ec6e180e..cb737bd0 100644 --- a/traits/src/tree.rs +++ b/traits/src/tree.rs @@ -44,14 +44,6 @@ pub trait Tree { where Self: 'a; - fn new( - points: Self::PointDataSlice<'_>, - adaptive: bool, - n_crit: Option, - depth: Option, - global_idxs: Self::GlobalIndexSlice<'_> - ) -> Self; - // Get depth of tree. fn get_depth(&self) -> u64; diff --git a/tree/examples/test_uniform.rs b/tree/examples/test_uniform.rs index bd75e406..ab29191c 100644 --- a/tree/examples/test_uniform.rs +++ b/tree/examples/test_uniform.rs @@ -1,31 +1,49 @@ -// //? mpirun -n {{NPROCESSES}} --features "mpi" +// ? mpirun -n {{NPROCESSES}} --features "mpi" -// use rand::prelude::*; -// use rand::SeedableRng; +use rand::prelude::*; +use rand::SeedableRng; -// use mpi::{environment::Universe, topology::UserCommunicator, traits::*}; +use mpi::{environment::Universe, topology::UserCommunicator, traits::*}; -// use bempp_traits::tree::Tree; +use bempp_traits::tree::Tree; -// use bempp_tree::types::{ -// domain::Domain, morton::MortonKey, multi_node::MultiNodeTree, point::PointType, -// }; +use rlst::{ + dense::{ + base_matrix::BaseMatrix, data_container::VectorContainer, matrix::Matrix, rlst_col_vec, + rlst_mat, rlst_pointer_mat, traits::*, Dot, global, + }, +}; -// pub fn points_fixture(npoints: i32) -> Vec<[f64; 3]> { -// let mut range = StdRng::seed_from_u64(0); -// let between = rand::distributions::Uniform::from(0.0..1.0); -// let mut points: Vec<[PointType; 3]> = Vec::new(); +use bempp_tree::types::{ + domain::Domain, morton::MortonKey, multi_node::MultiNodeTree, point::PointType, +}; -// for _ in 0..npoints { -// points.push([ -// between.sample(&mut range), -// between.sample(&mut range), -// between.sample(&mut range), -// ]) -// } +fn points_fixture( + npoints: usize, + min: Option, + max: Option, +) -> Matrix, Dynamic, Dynamic>, Dynamic, Dynamic> +{ + // Generate a set of randomly distributed points + let mut range = StdRng::seed_from_u64(0); -// points -// } + let between; + if let (Some(min), Some(max)) = (min, max) { + between = rand::distributions::Uniform::from(min..max); + } else { + between = rand::distributions::Uniform::from(0.0_f64..1.0_f64); + } + + let mut points = rlst_mat![f64, (npoints, 3)]; + + for i in 0..npoints { + points[[i, 0]] = between.sample(&mut range); + points[[i, 1]] = between.sample(&mut range); + points[[i, 2]] = between.sample(&mut range); + } + + points +} // /// Test that the leaves on separate nodes do not overlap. // fn test_no_overlaps(world: &UserCommunicator, tree: &MultiNodeTree) { @@ -82,34 +100,41 @@ // } // } -// fn main() { -// // Setup an MPI environment -// let universe: Universe = mpi::initialize().unwrap(); -// let world = universe.world(); -// let comm = world.duplicate(); - -// // Setup tree parameters -// let adaptive = false; -// let n_crit: Option<_> = None; -// let k: Option<_> = None; -// let depth = Some(4); -// let n_points = 10000; - -// let points = points_fixture(n_points); - -// let tree = MultiNodeTree::new(&comm, k, &points, adaptive, n_crit, depth); - -// test_global_bounds(&comm); -// if world.rank() == 0 { -// println!("\t ... test_global_bounds passed on uniform tree"); -// } -// test_uniform(&tree); -// if world.rank() == 0 { -// println!("\t ... test_uniform passed on uniform tree"); -// } -// test_no_overlaps(&comm, &tree); -// if world.rank() == 0 { -// println!("\t ... test_no_overlaps passed on uniform tree"); -// } -// } -fn main() {} +fn main() { + // Setup an MPI environment + let universe: Universe = mpi::initialize().unwrap(); + let world = universe.world(); + let comm = world.duplicate(); + + // Setup tree parameters + let adaptive = false; + let n_crit: Option<_> = None; + let k = 2; + let depth = Some(3); + let n_points = 10000; + + // Generate some random test data local to each process + let points = points_fixture(n_points, None, None); + let global_idxs: Vec<_> = (0..n_points).collect(); + + // Calculate the global domain + let domain = Domain::from_global_points(points.data(), &comm); + + // Create a uniform tree + let tree = MultiNodeTree::new(&comm, &points.data(), adaptive, n_crit, depth, k, &global_idxs); + + println!("Rank {:?}, leaves {:?}", world.rank(), tree.leaves.len()) + // test_global_bounds(&comm); + // if world.rank() == 0 { + // println!("\t ... test_global_bounds passed on uniform tree"); + // } + // test_uniform(&tree); + // if world.rank() == 0 { + // println!("\t ... test_uniform passed on uniform tree"); + // } + // test_no_overlaps(&comm, &tree); + // if world.rank() == 0 { + // println!("\t ... test_no_overlaps passed on uniform tree"); + // } +} +// fn main() {} diff --git a/tree/src/implementations/impl_domain_mpi.rs b/tree/src/implementations/impl_domain_mpi.rs index 3713f883..9388fa60 100644 --- a/tree/src/implementations/impl_domain_mpi.rs +++ b/tree/src/implementations/impl_domain_mpi.rs @@ -1,82 +1,82 @@ -// use memoffset::offset_of; -// use mpi::{ -// datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, -// topology::UserCommunicator, -// traits::*, -// Address, -// }; +use memoffset::offset_of; +use mpi::{ + datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, + topology::UserCommunicator, + traits::*, + Address, +}; -// use crate::types::{domain::Domain, point::PointType}; +use crate::types::{domain::Domain, point::PointType}; -// unsafe impl Equivalence for Domain { -// type Out = UserDatatype; -// fn equivalent_datatype() -> Self::Out { -// UserDatatype::structured( -// &[1, 1], -// &[ -// offset_of!(Domain, origin) as Address, -// offset_of!(Domain, diameter) as Address, -// ], -// &[ -// UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), -// UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), -// ], -// ) -// } -// } +unsafe impl Equivalence for Domain { + type Out = UserDatatype; + fn equivalent_datatype() -> Self::Out { + UserDatatype::structured( + &[1, 1], + &[ + offset_of!(Domain, origin) as Address, + offset_of!(Domain, diameter) as Address, + ], + &[ + UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), + UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), + ], + ) + } +} -// impl Domain { -// /// Compute the points domain over all nodes. -// pub fn from_global_points(local_points: &[[PointType; 3]], comm: &UserCommunicator) -> Domain { -// let size = comm.size(); +impl Domain { + /// Compute the points domain over all nodes. + pub fn from_global_points(local_points: &[PointType], comm: &UserCommunicator) -> Domain { + let size = comm.size(); -// let local_domain = Domain::from_local_points(local_points); -// let local_bounds: Vec = vec![local_domain; size as usize]; -// let mut buffer = vec![Domain::default(); size as usize]; + let local_domain = Domain::from_local_points(local_points); + let local_bounds: Vec = vec![local_domain; size as usize]; + let mut buffer = vec![Domain::default(); size as usize]; -// comm.all_to_all_into(&local_bounds, &mut buffer[..]); + comm.all_to_all_into(&local_bounds, &mut buffer[..]); -// // Find minimum origin -// let min_x = buffer -// .iter() -// .min_by(|a, b| a.origin[0].partial_cmp(&b.origin[0]).unwrap()) -// .unwrap() -// .origin[0]; -// let min_y = buffer -// .iter() -// .min_by(|a, b| a.origin[1].partial_cmp(&b.origin[1]).unwrap()) -// .unwrap() -// .origin[1]; -// let min_z = buffer -// .iter() -// .min_by(|a, b| a.origin[2].partial_cmp(&b.origin[2]).unwrap()) -// .unwrap() -// .origin[2]; + // Find minimum origin + let min_x = buffer + .iter() + .min_by(|a, b| a.origin[0].partial_cmp(&b.origin[0]).unwrap()) + .unwrap() + .origin[0]; + let min_y = buffer + .iter() + .min_by(|a, b| a.origin[1].partial_cmp(&b.origin[1]).unwrap()) + .unwrap() + .origin[1]; + let min_z = buffer + .iter() + .min_by(|a, b| a.origin[2].partial_cmp(&b.origin[2]).unwrap()) + .unwrap() + .origin[2]; -// let min_origin = [min_x, min_y, min_z]; + let min_origin = [min_x, min_y, min_z]; -// // Find maximum diameter (+max origin) -// let max_x = buffer -// .iter() -// .max_by(|a, b| a.diameter[0].partial_cmp(&b.diameter[0]).unwrap()) -// .unwrap() -// .diameter[0]; -// let max_y = buffer -// .iter() -// .max_by(|a, b| a.diameter[1].partial_cmp(&b.diameter[1]).unwrap()) -// .unwrap() -// .diameter[1]; -// let max_z = buffer -// .iter() -// .max_by(|a, b| a.diameter[2].partial_cmp(&b.diameter[2]).unwrap()) -// .unwrap() -// .diameter[2]; + // Find maximum diameter (+max origin) + let max_x = buffer + .iter() + .max_by(|a, b| a.diameter[0].partial_cmp(&b.diameter[0]).unwrap()) + .unwrap() + .diameter[0]; + let max_y = buffer + .iter() + .max_by(|a, b| a.diameter[1].partial_cmp(&b.diameter[1]).unwrap()) + .unwrap() + .diameter[1]; + let max_z = buffer + .iter() + .max_by(|a, b| a.diameter[2].partial_cmp(&b.diameter[2]).unwrap()) + .unwrap() + .diameter[2]; -// let max_diameter = [max_x, max_y, max_z]; + let max_diameter = [max_x, max_y, max_z]; -// Domain { -// origin: min_origin, -// diameter: max_diameter, -// } -// } -// } + Domain { + origin: min_origin, + diameter: max_diameter, + } + } +} diff --git a/tree/src/implementations/impl_morton_mpi.rs b/tree/src/implementations/impl_morton_mpi.rs index 1374ea60..8118c042 100644 --- a/tree/src/implementations/impl_morton_mpi.rs +++ b/tree/src/implementations/impl_morton_mpi.rs @@ -1,23 +1,23 @@ -// use crate::types::morton::{KeyType, MortonKey}; -// use memoffset::offset_of; -// use mpi::{ -// datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, -// Address, -// }; +use crate::types::morton::{KeyType, MortonKey}; +use memoffset::offset_of; +use mpi::{ + datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, + Address, +}; -// unsafe impl Equivalence for MortonKey { -// type Out = UserDatatype; -// fn equivalent_datatype() -> Self::Out { -// UserDatatype::structured( -// &[1, 1], -// &[ -// offset_of!(MortonKey, anchor) as Address, -// offset_of!(MortonKey, morton) as Address, -// ], -// &[ -// UncommittedUserDatatype::contiguous(3, &KeyType::equivalent_datatype()).as_ref(), -// UncommittedUserDatatype::contiguous(1, &KeyType::equivalent_datatype()).as_ref(), -// ], -// ) -// } -// } +unsafe impl Equivalence for MortonKey { + type Out = UserDatatype; + fn equivalent_datatype() -> Self::Out { + UserDatatype::structured( + &[1, 1], + &[ + offset_of!(MortonKey, anchor) as Address, + offset_of!(MortonKey, morton) as Address, + ], + &[ + UncommittedUserDatatype::contiguous(3, &KeyType::equivalent_datatype()).as_ref(), + UncommittedUserDatatype::contiguous(1, &KeyType::equivalent_datatype()).as_ref(), + ], + ) + } +} diff --git a/tree/src/implementations/impl_multi_node.rs b/tree/src/implementations/impl_multi_node.rs index dbc4ab3d..08c18a99 100644 --- a/tree/src/implementations/impl_multi_node.rs +++ b/tree/src/implementations/impl_multi_node.rs @@ -1,417 +1,554 @@ -// use itertools::Itertools; -// use std::collections::{HashMap, HashSet}; - -// use mpi::{collective::SystemOperation, topology::UserCommunicator, traits::*, Count, Rank}; - -// use hyksort::hyksort; - -// use bempp_traits::tree::{FmmData, FmmTree, Tree}; - -// use crate::{ -// constants::{DEEPEST_LEVEL, K, LEVEL_SIZE, NCRIT, ROOT}, -// implementations::{ -// impl_morton::{complete_region, encode_anchor}, -// impl_single_node::{ -// assign_nodes_to_points, assign_points_to_nodes, find_seeds, split_blocks, -// }, -// mpi_helpers::all_to_allv_sparse, -// }, -// types::{ -// data::NodeData, -// domain::Domain, -// morton::{KeyType, MortonKey, MortonKeys}, -// multi_node::MultiNodeTree, -// point::{Point, PointType, Points}, -// }, -// }; - -// impl MultiNodeTree { -// /// Create a new MultiNodeTree from a set of distributed points which define a domain. -// pub fn new( -// world: &UserCommunicator, -// k: Option, -// points: &[[PointType; 3]], -// adaptive: bool, -// n_crit: Option, -// depth: Option, -// ) -> MultiNodeTree { -// let domain = Domain::from_global_points(points, world); - -// let n_crit = n_crit.unwrap_or(NCRIT); -// let depth = depth.unwrap_or(DEEPEST_LEVEL); -// let k = k.unwrap_or(K); - -// if adaptive { -// MultiNodeTree::adaptive_tree(world, k, points, &domain, n_crit) -// } else { -// MultiNodeTree::uniform_tree(world, k, points, &domain, depth) -// } -// } - -// /// Specialization for uniform trees -// pub fn uniform_tree( -// world: &UserCommunicator, -// k: i32, -// points: &[[PointType; 3]], -// domain: &Domain, -// depth: u64, -// ) -> MultiNodeTree { -// // Encode points at specified depth -// let mut points = points -// .iter() -// .enumerate() -// .map(|(i, p)| Point { -// coordinate: *p, -// key: MortonKey::from_point(p, &domain, depth), -// global_idx: i, -// }) -// .collect(); - -// // 2.i Perform parallel Morton sort over encoded points -// let comm = world.duplicate(); -// hyksort(&mut points, k, comm); - -// // 2.ii Find leaf keys on each processor -// let min = points.iter().min().unwrap().key; -// let max = points.iter().max().unwrap().key; - -// let diameter = 1 << (DEEPEST_LEVEL - depth as u64); - -// let mut leaves = MortonKeys { -// keys: (0..LEVEL_SIZE) -// .step_by(diameter) -// .flat_map(|i| (0..LEVEL_SIZE).step_by(diameter).map(move |j| (i, j))) -// .flat_map(|(i, j)| (0..LEVEL_SIZE).step_by(diameter).map(move |k| [i, j, k])) -// .map(|anchor| { -// let morton = encode_anchor(&anchor, depth as u64); -// MortonKey { anchor, morton } -// }) -// .filter(|k| (k >= &min) && (k <= &max)) -// .collect(), -// index: 0, -// }; - -// // 3. Create bi-directional maps between keys and points -// let points_to_leaves = assign_points_to_nodes(&points, &leaves); -// let leaves_to_points = assign_nodes_to_points(&leaves, &points); - -// // Only retain keys that contain points -// leaves = MortonKeys { -// keys: leaves_to_points.keys().cloned().collect(), -// index: 0, -// }; - -// let leaves_set: HashSet = leaves.iter().cloned().collect(); - -// let mut keys_set: HashSet = HashSet::new(); -// for key in leaves.iter() { -// let ancestors = key.ancestors(); -// keys_set.extend(&ancestors); -// } - -// let min = leaves.iter().min().unwrap(); -// let max = leaves.iter().max().unwrap(); -// let range = [world.rank() as KeyType, min.morton, max.morton]; - -// let keys_to_data: HashMap = HashMap::new(); - -// MultiNodeTree { -// world: world.duplicate(), -// adaptive: false, -// points, -// keys_set, -// leaves, -// leaves_set, -// domain: *domain, -// points_to_leaves, -// leaves_to_points, -// keys_to_data, -// range, -// } -// } - -// /// Specialization for adaptive tree. -// pub fn adaptive_tree( -// world: &UserCommunicator, -// k: i32, -// points: &[[PointType; 3]], -// domain: &Domain, -// n_crit: usize, -// ) -> MultiNodeTree { -// // 1. Encode Points to Leaf Morton Keys, add a global index related to the processor -// let mut points: Points = points -// .iter() -// .enumerate() -// .map(|(i, p)| Point { -// coordinate: *p, -// global_idx: i, -// key: MortonKey::from_point(p, domain, DEEPEST_LEVEL), -// }) -// .collect(); - -// // 2.i Perform parallel Morton sort over encoded points -// let comm = world.duplicate(); - -// hyksort(&mut points, k, comm); - -// // 2.ii Find unique leaf keys on each processor -// let mut local = MortonKeys { -// keys: points.iter().map(|p| p.key).collect(), -// index: 0, -// }; - -// // 3. Linearise received keys (remove overlaps if they exist). -// local.linearize(); - -// // 4. Complete region spanned by node. -// local.complete(); - -// // 5.i Find seeds and compute the coarse blocktree -// let mut seeds = find_seeds(&local); - -// let blocktree = MultiNodeTree::complete_blocktree(world, &mut seeds); - -// // 5.ii any data below the min seed sent to partner process -// let points = MultiNodeTree::transfer_points_to_blocktree(world, &points, &blocktree); - -// // 6. Split blocks based on ncrit constraint -// let mut locally_balanced = split_blocks(&points, blocktree, n_crit); - -// locally_balanced.sort(); - -// // 7. Create a minimal balanced octree for local octants spanning their domain and linearize -// locally_balanced.balance(); -// locally_balanced.linearize(); - -// // 8. Find new maps between points and locally balanced tree - -// let points_to_locally_balanced = assign_points_to_nodes(&points, &locally_balanced); -// let mut points: Points = points -// .iter() -// .map(|p| Point { -// coordinate: p.coordinate, -// global_idx: p.global_idx, -// key: *points_to_locally_balanced.get(p).unwrap(), -// }) -// .collect(); - -// // 9. Perform another distributed sort and remove overlaps locally -// let comm = world.duplicate(); - -// hyksort(&mut points, k, comm); - -// let mut globally_balanced = MortonKeys { -// keys: points.iter().map(|p| p.key).collect(), -// index: 0, -// }; -// globally_balanced.linearize(); - -// // 10. Find final bidirectional maps to non-overlapping tree -// let points_to_globally_balanced = assign_points_to_nodes(&points, &globally_balanced); -// let globally_balanced_to_points = assign_nodes_to_points(&globally_balanced, &points); -// let points: Points = points -// .iter() -// .map(|p| Point { -// coordinate: p.coordinate, -// global_idx: p.global_idx, -// key: *points_to_globally_balanced.get(p).unwrap(), -// }) -// .collect(); - -// let leaves_set: HashSet = globally_balanced.iter().cloned().collect(); -// let mut keys_set: HashSet = HashSet::new(); -// for key in globally_balanced.iter() { -// let ancestors = key.ancestors(); -// keys_set.extend(&ancestors); -// } - -// let min = globally_balanced.iter().min().unwrap(); -// let max = globally_balanced.iter().max().unwrap(); -// let range = [world.rank() as KeyType, min.morton, max.morton]; - -// let keys_to_data: HashMap = HashMap::new(); - -// MultiNodeTree { -// world: world.duplicate(), -// adaptive: true, -// points, -// keys_set, -// leaves: globally_balanced, -// leaves_set, -// domain: *domain, -// points_to_leaves: points_to_globally_balanced, -// leaves_to_points: globally_balanced_to_points, -// keys_to_data, -// range, -// } -// } - -// /// Complete a distributed block tree from the seed octants, algorithm 4 in [1] (parallel). -// fn complete_blocktree(world: &UserCommunicator, seeds: &mut MortonKeys) -> MortonKeys { -// let rank = world.rank(); -// let size = world.size(); - -// // Define the tree's global domain with the finest first/last descendants -// if rank == 0 { -// let ffc_root = ROOT.finest_first_child(); -// let min = seeds.iter().min().unwrap(); -// let fa = ffc_root.finest_ancestor(min); -// let first_child = fa.children().into_iter().min().unwrap(); -// // Check for overlap -// if first_child < *min { -// seeds.push(first_child) -// } -// seeds.sort(); -// } - -// if rank == (size - 1) { -// let flc_root = ROOT.finest_last_child(); -// let max = seeds.iter().max().unwrap(); -// let fa = flc_root.finest_ancestor(max); -// let last_child = fa.children().into_iter().max().unwrap(); - -// if last_child > *max -// && !max.ancestors().contains(&last_child) -// && !last_child.ancestors().contains(max) -// { -// seeds.push(last_child); -// } -// } - -// let next_rank = if rank + 1 < size { rank + 1 } else { 0 }; -// let previous_rank = if rank > 0 { rank - 1 } else { size - 1 }; - -// let previous_process = world.process_at_rank(previous_rank); -// let next_process = world.process_at_rank(next_rank); - -// // Send required data to partner process. -// if rank > 0 { -// let min = *seeds.iter().min().unwrap(); -// previous_process.send(&min); -// } - -// let mut boundary = MortonKey::default(); - -// if rank < (size - 1) { -// next_process.receive_into(&mut boundary); -// seeds.push(boundary); -// } - -// // Complete region between seeds at each process -// let mut complete = MortonKeys::new(); +use itertools::Itertools; +use std::collections::{HashMap, HashSet}; + +use mpi::{collective::SystemOperation, topology::UserCommunicator, traits::*, Count, Rank}; + +use hyksort::hyksort; + +use bempp_traits::tree::Tree; + +use crate::{ + constants::{DEEPEST_LEVEL, K, LEVEL_SIZE, NCRIT, ROOT, LEVEL_DISPLACEMENT}, + implementations::{ + impl_morton::{complete_region, encode_anchor}, + mpi_helpers::all_to_allv_sparse, + }, + types::{ + domain::Domain, + morton::{KeyType, MortonKey, MortonKeys}, + multi_node::MultiNodeTree, + single_node::SingleNodeTree, + point::{Point, PointType, Points, self}, + }, +}; + +impl MultiNodeTree { + /// Constructor for uniform trees + pub fn uniform_tree( + world: &UserCommunicator, + k: i32, + points: &[PointType], + domain: &Domain, + depth: u64, + global_idxs: &[usize] + ) -> MultiNodeTree { + // Encode points at deepest level, and map to specified depth. + let dim = 3; + let npoints = points.len() / dim; + + let mut tmp = Points::default(); + for i in 0..npoints { + let point = [points[i], points[i + npoints], points[i + 2 * npoints]]; + let base_key = MortonKey::from_point(&point, &domain, DEEPEST_LEVEL); + let encoded_key = MortonKey::from_point(&point, &domain, depth); + tmp.points.push(Point { + coordinate: point, + base_key, + encoded_key, + global_idx: global_idxs[i], + }) + } + let mut points = tmp; + + // 2.i Perform parallel Morton sort over encoded points + let comm = world.duplicate(); + hyksort(&mut points.points, k, comm); + + // 2.ii Find leaf keys on each processor + let min = points.points.iter().min().unwrap().encoded_key; + let max = points.points.iter().max().unwrap().encoded_key; + + let diameter = 1 << (DEEPEST_LEVEL - depth as u64); + + // Find leaves within ths processor's range + let leaves = MortonKeys { + keys: (min.anchor[0]..max.anchor[0]) + .step_by(diameter) + .flat_map(|i| (min.anchor[1]..max.anchor[1]).step_by(diameter).map(move |j| (i, j))) + .flat_map(|(i, j)| (min.anchor[2]..max.anchor[2]).step_by(diameter).map(move |k| [i, j, k])) + .map(|anchor| { + let morton = encode_anchor(&anchor, depth); + MortonKey { anchor, morton } + }) + .collect(), + index: 0, + }; + + // 3. Assign keys to points + let unmapped = SingleNodeTree::assign_nodes_to_points(&leaves, &mut points); + + // Group points by leaves + points.sort(); + + let mut leaves_to_points = HashMap::new(); + let mut curr = points.points[0].clone(); + let mut curr_idx = 0; + + for (i, point) in points.points.iter().enumerate() { + if point.encoded_key != curr.encoded_key { + leaves_to_points.insert(curr.encoded_key, (curr_idx, i)); + curr_idx = i; + curr = point.clone(); + } + } + leaves_to_points.insert(curr.encoded_key, (curr_idx, points.points.len())); + + // Add unmapped leaves + let leaves = MortonKeys { + keys: leaves_to_points + .keys() + .cloned() + .chain(unmapped.iter().copied()) + .collect_vec(), + index: 0, + }; + + // Find all keys in tree + let tmp: HashSet = leaves + .iter() + .flat_map(|leaf| leaf.ancestors().into_iter()) + .collect(); + + let mut keys = MortonKeys { + keys: tmp.into_iter().collect_vec(), + index: 0, + }; + + let leaves_set: HashSet = leaves.iter().cloned().collect(); + let keys_set: HashSet = keys.iter().cloned().collect(); + + let min = leaves.iter().min().unwrap(); + let max = leaves.iter().max().unwrap(); + let range = [world.rank() as KeyType, min.morton, max.morton]; + + // Group by level to perform efficient lookup of nodes + keys.sort_by_key(|a| a.level()); + + let mut levels_to_keys = HashMap::new(); + let mut curr = keys[0]; + let mut curr_idx = 0; + for (i, key) in keys.iter().enumerate() { + if key.level() != curr.level() { + levels_to_keys.insert(curr.level(), (curr_idx, i)); + curr_idx = i; + curr = *key; + } + } + levels_to_keys.insert(curr.level(), (curr_idx, keys.len())); + + // Return tree in sorted order + for l in 0..=depth { + let &(l, r) = levels_to_keys.get(&l).unwrap(); + let subset = &mut keys[l..r]; + subset.sort(); + } + + MultiNodeTree { + world: world.duplicate(), + depth, + domain: *domain, + points, + leaves, + keys, + leaves_to_points, + levels_to_keys, + leaves_set, + keys_set, + range, + } + } + + /// Constructor for adaptive tree. + pub fn adaptive_tree( + world: &UserCommunicator, + k: i32, + points: &[PointType], + domain: &Domain, + n_crit: u64, + global_idxs: &[usize] + ) -> MultiNodeTree + { + // 1. Encode Points to Leaf Morton Keys, add a global index related to the processor + let dim = 3; + let npoints = points.len() / dim; + + let mut tmp = Points::default(); + for i in 0..npoints { + let point = [points[i], points[i + npoints], points[i + 2 * npoints]]; + let key = MortonKey::from_point(&point, &domain, DEEPEST_LEVEL); + tmp.points.push(Point { + coordinate: point, + base_key: key, + encoded_key: key, + global_idx: global_idxs[i], + }) + } + let mut points = tmp; + + // 2.i Perform parallel Morton sort over encoded points + let comm = world.duplicate(); + + hyksort(&mut points.points, k, comm); + + // 2.ii Find unique leaf keys on each processor + let mut local = MortonKeys { + keys: points.points.iter().map(|p| p.encoded_key).collect(), + index: 0, + }; + + // 3. Linearise received keys (remove overlaps if they exist). + local.linearize(); + + // 4. Complete region spanned by node. + local.complete(); + + // 5.i Find seeds and compute the coarse blocktree + let mut seeds = SingleNodeTree::find_seeds(&local); + + let blocktree = MultiNodeTree::complete_blocktree(world, &mut seeds); + + // 5.ii any data below the min seed sent to partner process + let mut points = MultiNodeTree::transfer_points_to_blocktree(world, &points.points[..], &blocktree); + + // 6. Split blocks based on ncrit constraint + let mut locally_balanced = SingleNodeTree::split_blocks(&mut points, blocktree, n_crit as usize); + + // 7. Create a minimal balanced octree for local octants spanning their domain and linearize + locally_balanced.sort(); + locally_balanced.balance(); + locally_balanced.linearize(); + + // 8. Find new maps between points and locally balanced tree + let unmapped = SingleNodeTree::assign_nodes_to_points(&locally_balanced, &mut points); + + + // 9. Perform another distributed sort and remove overlaps locally + let comm = world.duplicate(); + + hyksort(&mut points.points, k, comm); + + let mut globally_balanced = MortonKeys { + keys: points.points.iter().map(|p| p.encoded_key).collect(), + index: 0, + }; + globally_balanced.linearize(); + + // Group points by leaves + points.sort(); + + let mut leaves_to_points = HashMap::new(); + let mut curr = points.points[0].clone(); + let mut curr_idx = 0; + + for (i, point) in points.points.iter().enumerate() { + if point.encoded_key != curr.encoded_key { + leaves_to_points.insert(curr.encoded_key, (curr_idx, i)); + curr_idx = i; + curr = point.clone(); + } + } + leaves_to_points.insert(curr.encoded_key, (curr_idx, points.points.len())); + + // 10. Find final maps to non-overlapping tree + let unmapped = SingleNodeTree::assign_nodes_to_points(&globally_balanced, &mut points); + + // Add unmapped leaves + let globaly_balanced = MortonKeys { + keys: leaves_to_points + .keys() + .cloned() + .chain(unmapped.iter().copied()) + .collect_vec(), + index: 0, + }; + + // Find all keys in tree + let tmp: HashSet = globally_balanced + .iter() + .flat_map(|leaf| leaf.ancestors().into_iter()) + .collect(); + + let mut keys = MortonKeys { + keys: tmp.into_iter().collect_vec(), + index: 0, + }; + + let leaves_set: HashSet = globally_balanced.iter().cloned().collect(); + let mut keys_set: HashSet = HashSet::new(); + + // Create maps between points and globally balanced tree + let mut leaves_to_points = HashMap::new(); + let mut curr = points.points[0].clone(); + let mut curr_idx = 0; + + for (i, point) in points.points.iter().enumerate() { + if point.encoded_key != curr.encoded_key { + leaves_to_points.insert(curr.encoded_key, (curr_idx, i)); + curr_idx = i; + curr = point.clone(); + } + } + leaves_to_points.insert(curr.encoded_key, (curr_idx, points.points.len())); + + for key in globally_balanced.iter() { + let ancestors = key.ancestors(); + keys_set.extend(&ancestors); + } + + let min = globally_balanced.iter().min().unwrap(); + let max = globally_balanced.iter().max().unwrap(); + let range = [world.rank() as KeyType, min.morton, max.morton]; + + // Group by level to perform efficient lookup of nodes + keys.sort_by_key(|a| a.level()); + + let mut levels_to_keys = HashMap::new(); + let mut curr = keys[0]; + let mut curr_idx = 0; + for (i, key) in keys.iter().enumerate() { + if key.level() != curr.level() { + levels_to_keys.insert(curr.level(), (curr_idx, i)); + curr_idx = i; + curr = *key; + } + } + levels_to_keys.insert(curr.level(), (curr_idx, keys.len())); + + // Find depth + let depth = points + .points + .iter() + .map(|p| p.encoded_key.level()) + .max() + .unwrap(); + + // Return tree in sorted order + for l in 0..=depth { + let &(l, r) = levels_to_keys.get(&l).unwrap(); + let subset = &mut keys[l..r]; + subset.sort(); + } + + MultiNodeTree { + world: world.duplicate(), + depth, + domain: *domain, + points, + leaves: globally_balanced, + keys, + leaves_to_points, + levels_to_keys, + leaves_set, + keys_set, + range + } + } + + + /// Create a new multi-node tree. If non-adaptive (uniform) trees are created, they are specified + /// by a user defined maximum depth, if an adaptive tree is created it is specified by only by the + /// user defined maximum leaf maximum occupancy n_crit. + pub fn new( + world: &UserCommunicator, + points: &[PointType], + adaptive: bool, + n_crit: Option, + depth: Option, + k: i32, + global_idxs: &[usize] + ) -> MultiNodeTree { + // TODO: Come back and reconcile a runtime point dimension detector + + let domain = Domain::from_global_points(points, &world); + + let n_crit = n_crit.unwrap_or(NCRIT); + let depth = depth.unwrap_or(DEEPEST_LEVEL); + + if adaptive { + MultiNodeTree::adaptive_tree(world, k, points, &domain, n_crit, global_idxs) + } else { + MultiNodeTree::uniform_tree(world, k, points, &domain, depth, &global_idxs) + } + } + + /// Complete a distributed block tree from the seed octants, algorithm 4 in [1] (parallel). + fn complete_blocktree(world: &UserCommunicator, seeds: &mut MortonKeys) -> MortonKeys { + let rank = world.rank(); + let size = world.size(); + + // Define the tree's global domain with the finest first/last descendants + if rank == 0 { + let ffc_root = ROOT.finest_first_child(); + let min = seeds.iter().min().unwrap(); + let fa = ffc_root.finest_ancestor(min); + let first_child = fa.children().into_iter().min().unwrap(); + // Check for overlap + if first_child < *min { + seeds.push(first_child) + } + seeds.sort(); + } + + if rank == (size - 1) { + let flc_root = ROOT.finest_last_child(); + let max = seeds.iter().max().unwrap(); + let fa = flc_root.finest_ancestor(max); + let last_child = fa.children().into_iter().max().unwrap(); + + if last_child > *max + && !max.ancestors().contains(&last_child) + && !last_child.ancestors().contains(max) + { + seeds.push(last_child); + } + } + + let next_rank = if rank + 1 < size { rank + 1 } else { 0 }; + let previous_rank = if rank > 0 { rank - 1 } else { size - 1 }; + + let previous_process = world.process_at_rank(previous_rank); + let next_process = world.process_at_rank(next_rank); + + // Send required data to partner process. + if rank > 0 { + let min = *seeds.iter().min().unwrap(); + previous_process.send(&min); + } + + let mut boundary = MortonKey::default(); + + if rank < (size - 1) { + next_process.receive_into(&mut boundary); + seeds.push(boundary); + } + + // Complete region between seeds at each process + let mut complete = MortonKeys::new(); + + for i in 0..(seeds.iter().len() - 1) { + let a = seeds[i]; + let b = seeds[i + 1]; + + let mut tmp: Vec = complete_region(&a, &b); + complete.keys.push(a); + complete.keys.append(&mut tmp); + } + + if rank == (size - 1) { + complete.keys.push(seeds.last().unwrap()); + } + + complete.sort(); + complete + } + + // Transfer points to correct processor based on the coarse distributed blocktree. + fn transfer_points_to_blocktree( + world: &UserCommunicator, + points: &[Point], + blocktree: &[MortonKey], + ) -> Points { + let rank = world.rank(); + let size = world.size(); + + let mut received_points = Vec::new(); + + let min = blocktree.iter().min().unwrap(); + + let prev_rank = if rank > 0 { rank - 1 } else { size - 1 }; + let next_rank = if rank + 1 < size { rank + 1 } else { 0 }; + + if rank > 0 { + let msg: Vec<_> = points.iter().filter(|&p| p.encoded_key < *min).cloned().collect_vec(); + + let msg_size: Rank = msg.len() as Rank; + world.process_at_rank(prev_rank).send(&msg_size); + world.process_at_rank(prev_rank).send(&msg[..]); + } + + if rank < (size - 1) { + let mut bufsize = 0; + world.process_at_rank(next_rank).receive_into(&mut bufsize); + let mut buffer = vec![Point::default(); bufsize as usize]; + world + .process_at_rank(next_rank) + .receive_into(&mut buffer[..]); + received_points.append(&mut buffer); + } + + // Filter out local points that's been sent to partner + received_points = points.iter().filter(|&p| p.encoded_key >= *min).cloned().collect(); + + received_points.sort(); + + Points { + points: received_points, + index: 0 + } + } +} + +impl Tree for MultiNodeTree { + type Domain = Domain; + type NodeIndex = MortonKey; + type NodeIndexSlice<'a> = &'a [MortonKey]; + type NodeIndices = MortonKeys; + type Point = Point; + type PointSlice<'a> = &'a [Point]; + type PointData = f64; + type PointDataSlice<'a> = &'a [f64]; + type GlobalIndex = usize; + type GlobalIndexSlice<'a> = &'a [usize]; + + fn get_depth(&self) -> u64 { + self.depth + } + + fn get_domain(&self) -> &'_ Self::Domain { + &self.domain + } + + fn get_keys(&self, level: u64) -> Option> { + if let Some(&(l, r)) = self.levels_to_keys.get(&level) { + Some(&self.keys[l..r]) + } else { + None + } + } + + fn get_all_keys(&self) -> Option> { + Some(&self.keys) + } + + fn get_all_keys_set(&self) -> &'_ HashSet { + &self.keys_set + } + + fn get_all_leaves_set(&self) -> &'_ HashSet { + &self.leaves_set + } + + fn get_leaves(&self) -> Option> { + Some(&self.leaves) + } + + fn get_points<'a>(&'a self, key: &Self::NodeIndex) -> Option> { + if let Some(&(l, r)) = self.leaves_to_points.get(key) { + Some(&self.points.points[l..r]) + } else { + None + } + } + + fn is_leaf(&self, key: &Self::NodeIndex) -> bool { + self.leaves_set.contains(key) + } + + fn is_node(&self, key: &Self::NodeIndex) -> bool { + self.keys_set.contains(key) + } +} -// for i in 0..(seeds.iter().len() - 1) { -// let a = seeds[i]; -// let b = seeds[i + 1]; - -// let mut tmp: Vec = complete_region(&a, &b); -// complete.keys.push(a); -// complete.keys.append(&mut tmp); -// } - -// if rank == (size - 1) { -// complete.keys.push(seeds.last().unwrap()); -// } - -// complete.sort(); -// complete -// } - -// // Transfer points to correct processor based on the coarse distributed blocktree. -// fn transfer_points_to_blocktree( -// world: &UserCommunicator, -// points: &[Point], -// blocktree: &[MortonKey], -// ) -> Points { -// let rank = world.rank(); -// let size = world.size(); - -// let mut received_points = Points::new(); - -// let min = blocktree.iter().min().unwrap(); - -// let prev_rank = if rank > 0 { rank - 1 } else { size - 1 }; -// let next_rank = if rank + 1 < size { rank + 1 } else { 0 }; - -// if rank > 0 { -// let msg: Points = points.iter().filter(|&p| p.key < *min).cloned().collect(); - -// let msg_size: Rank = msg.len() as Rank; -// world.process_at_rank(prev_rank).send(&msg_size); -// world.process_at_rank(prev_rank).send(&msg[..]); -// } - -// if rank < (size - 1) { -// let mut bufsize = 0; -// world.process_at_rank(next_rank).receive_into(&mut bufsize); -// let mut buffer = vec![Point::default(); bufsize as usize]; -// world -// .process_at_rank(next_rank) -// .receive_into(&mut buffer[..]); -// received_points.append(&mut buffer); -// } - -// // Filter out local points that's been sent to partner -// received_points = points.iter().filter(|&p| p.key >= *min).cloned().collect(); - -// received_points.sort(); -// received_points -// } -// } - -// impl Tree for MultiNodeTree { -// type Domain = Domain; -// type Point = Point; -// type Points = Points; -// type NodeIndex = MortonKey; -// type NodeIndices = MortonKeys; -// type NodeIndicesSet = HashSet; -// type NodeDataType = NodeData; - -// // Get adaptivity information -// fn get_adaptive(&self) -> bool { -// self.adaptive -// } - -// // Get all keys, gets local keys in multi-node setting -// fn get_keys(&self) -> &MortonKeys { -// &self.leaves -// } - -// fn get_keys_set(&self) -> &HashSet { -// &self.leaves_set -// } - -// // Get all points, gets local keys in multi-node setting -// fn get_all_points(&self) -> &Points { -// &self.points -// } - -// // Get domain, gets global domain in multi-node setting -// fn get_domain(&self) -> &Domain { -// &self.domain -// } - -// // Get leaf key associated with a point -// fn get_leaf(&self, point: &Point) -> Option<&MortonKey> { -// self.points_to_leaves.get(point) -// } - -// // Get points associated with a tree leaf -// fn get_points(&self, leaf: &MortonKey) -> Option<&Points> { -// self.leaves_to_points.get(leaf) -// } - -// // Set data associated with a tree node key -// fn set_data(&mut self, node_index: &Self::NodeIndex, data: Self::NodeDataType) { -// self.keys_to_data.insert(*node_index, data); -// } - -// // Get data associated with a tree node key -// fn get_data(&self, key: &Self::NodeIndex) -> Option<&Self::NodeDataType> { -// self.keys_to_data.get(key) -// } -// } // // Helper function for creating locally essential trees. // fn let_helper(tree: &mut MultiNodeTree) { @@ -755,151 +892,3 @@ // tree.points_to_leaves = assign_points_to_nodes(&tree.points, &tree.leaves); // tree.leaves_to_points = assign_nodes_to_points(&tree.leaves, &tree.points); // } - -// impl FmmTree for MultiNodeTree { -// type NodeIndex = MortonKey; -// type NodeIndices = MortonKeys; -// type NodeData = NodeData; -// type NodeDataType = Vec; - -// fn create_let(&mut self) { -// // Create an LET -// let_helper(self); -// // Load balance the LET -// load_balance_let(self); -// // Reform LET based, now load balanced. -// let_helper(self); -// } - -// // Calculate near field interaction list of keys. -// fn get_near_field(&self, leaf: &Self::NodeIndex) -> Option { -// let mut keys = Vec::::new(); -// let neighbours = leaf.neighbors(); - -// // Child level -// let mut neighbors_children_adj: Vec = neighbours -// .iter() -// .flat_map(|n| n.children()) -// .filter(|nc| self.keys_set.contains(nc) && leaf.is_adjacent(nc)) -// .collect(); - -// // Key level -// let mut neighbors_adj: Vec = neighbours -// .iter() -// .filter(|n| self.keys_set.contains(n) && leaf.is_adjacent(n)) -// .cloned() -// .collect(); - -// // Parent level -// let mut neighbors_parents_adj: Vec = neighbours -// .iter() -// .map(|n| n.parent()) -// .filter(|np| self.leaves_set.contains(np) && leaf.is_adjacent(np)) -// .collect(); - -// keys.append(&mut neighbors_children_adj); -// keys.append(&mut neighbors_adj); -// keys.append(&mut neighbors_parents_adj); - -// if keys.len() > 0 { -// Some(MortonKeys { keys, index: 0 }) -// } else { -// None -// } -// } - -// // Calculate compressible far field interactions of leaf & other keys. -// fn get_interaction_list(&self, key: &Self::NodeIndex) -> Option { -// if key.level() >= 2 { -// let keys = key -// .parent() -// .neighbors() -// .iter() -// .flat_map(|pn| pn.children()) -// .filter(|pnc| self.keys_set.contains(pnc) && key.is_adjacent(pnc)) -// .collect_vec(); - -// if keys.len() > 0 { -// return Some(MortonKeys { keys, index: 0 }); -// } else { -// return None; -// } -// } -// { -// None -// } -// } - -// // Calculate M2P interactions of leaf key. -// fn get_w_list(&self, leaf: &Self::NodeIndex) -> Option { -// // Child level -// let keys = leaf -// .neighbors() -// .iter() -// .flat_map(|n| n.children()) -// .filter(|nc| self.keys_set.contains(nc) && !leaf.is_adjacent(nc)) -// .collect_vec(); - -// if keys.len() > 0 { -// Some(MortonKeys { keys, index: 0 }) -// } else { -// None -// } -// } - -// // Calculate P2L interactions of leaf key. -// fn get_x_list(&self, leaf: &Self::NodeIndex) -> Option { -// let keys = leaf -// .parent() -// .neighbors() -// .into_iter() -// .filter(|pn| self.keys_set.contains(pn) && !leaf.is_adjacent(pn)) -// .collect_vec(); - -// if keys.len() > 0 { -// Some(MortonKeys { keys, index: 0 }) -// } else { -// None -// } -// } - -// // Set data associated with a tree node key -// fn set_multipole_expansion(&mut self, node_index: &Self::NodeIndex, data: &Self::NodeDataType) { -// if let Some(x) = self.keys_to_data.get_mut(node_index) { -// x.set_multipole_expansion(data); -// } -// } - -// // Get data associated with a tree node key -// fn get_multipole_expansion(&self, node_index: &Self::NodeIndex) -> Option { -// if let Some(x) = self.keys_to_data.get(node_index) { -// Some(x.get_multipole_expansion()) -// } else { -// None -// } -// } - -// fn set_local_expansion(&mut self, node_index: &Self::NodeIndex, data: &Self::NodeDataType) { -// if let Some(x) = self.keys_to_data.get_mut(node_index) { -// x.set_local_expansion(data); -// } -// } - -// // Get data associated with a tree node key -// fn get_local_expansion(&self, node_index: &Self::NodeIndex) -> Option { -// if let Some(x) = self.keys_to_data.get(node_index) { -// Some(x.get_local_expansion()) -// } else { -// None -// } -// } - -// // TODO: Not implemented -// fn downward_pass(&self) {} - -// // TODO: Not implemented -// fn upward_pass(&self) {} - -// // TODO: Not implemented -// fn run(&self, _expansion_order: usize) {} -// } diff --git a/tree/src/implementations/impl_point.rs b/tree/src/implementations/impl_point.rs index 6664cec2..5b2b7ac5 100644 --- a/tree/src/implementations/impl_point.rs +++ b/tree/src/implementations/impl_point.rs @@ -71,6 +71,7 @@ impl FromIterator for Points { } } + #[cfg(test)] mod test { use rand::prelude::*; diff --git a/tree/src/implementations/impl_point_mpi.rs b/tree/src/implementations/impl_point_mpi.rs index e96b4a5e..483749f5 100644 --- a/tree/src/implementations/impl_point_mpi.rs +++ b/tree/src/implementations/impl_point_mpi.rs @@ -1,44 +1,57 @@ -// use crate::types::{ -// morton::{KeyType, MortonKey}, -// point::{Point, PointType}, -// }; -// use memoffset::offset_of; -// use mpi::{ -// datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, -// Address, -// }; +use crate::types::{ + morton::{KeyType, MortonKey}, + point::{Point, PointType}, +}; +use memoffset::offset_of; +use mpi::{ + datatype::{Equivalence, UncommittedUserDatatype, UserDatatype}, + Address, +}; -// unsafe impl Equivalence for Point { -// type Out = UserDatatype; -// fn equivalent_datatype() -> Self::Out { -// UserDatatype::structured( -// &[1, 1, 1], -// &[ -// offset_of!(Point, coordinate) as Address, -// offset_of!(Point, global_idx) as Address, -// offset_of!(Point, key) as Address, -// offset_of!(Point, data) as Address, -// ], -// &[ -// UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), -// UncommittedUserDatatype::contiguous(1, &usize::equivalent_datatype()).as_ref(), -// UncommittedUserDatatype::structured( -// &[1, 1], -// &[ -// offset_of!(MortonKey, anchor) as Address, -// offset_of!(MortonKey, morton) as Address, -// ], -// &[ -// UncommittedUserDatatype::contiguous(3, &KeyType::equivalent_datatype()) -// .as_ref(), -// UncommittedUserDatatype::contiguous(1, &KeyType::equivalent_datatype()) -// .as_ref(), -// ], -// ) -// .as_ref(), -// UncommittedUserDatatype::contiguous(1, &::equivalent_datatype()).as_ref(), -// ], -// UncommittedUserDatatype::contiguous(1, &PointType::equivalent_datatype()).as_ref(), -// ) -// } -// } +unsafe impl Equivalence for Point { + type Out = UserDatatype; + fn equivalent_datatype() -> Self::Out { + UserDatatype::structured( + &[1, 1, 1, 1], + &[ + offset_of!(Point, coordinate) as Address, + offset_of!(Point, global_idx) as Address, + offset_of!(Point, base_key) as Address, + offset_of!(Point, encoded_key) as Address, + + ], + &[ + UncommittedUserDatatype::contiguous(3, &PointType::equivalent_datatype()).as_ref(), + UncommittedUserDatatype::contiguous(1, &usize::equivalent_datatype()).as_ref(), + UncommittedUserDatatype::structured( + &[1, 1], + &[ + offset_of!(MortonKey, anchor) as Address, + offset_of!(MortonKey, morton) as Address, + ], + &[ + UncommittedUserDatatype::contiguous(3, &KeyType::equivalent_datatype()) + .as_ref(), + UncommittedUserDatatype::contiguous(1, &KeyType::equivalent_datatype()) + .as_ref(), + ], + ) + .as_ref(), + UncommittedUserDatatype::structured( + &[1, 1], + &[ + offset_of!(MortonKey, anchor) as Address, + offset_of!(MortonKey, morton) as Address, + ], + &[ + UncommittedUserDatatype::contiguous(3, &KeyType::equivalent_datatype()) + .as_ref(), + UncommittedUserDatatype::contiguous(1, &KeyType::equivalent_datatype()) + .as_ref(), + ], + ) + .as_ref(), + ], + ) + } +} diff --git a/tree/src/implementations/impl_single_node.rs b/tree/src/implementations/impl_single_node.rs index 997fb055..0bacc1b8 100644 --- a/tree/src/implementations/impl_single_node.rs +++ b/tree/src/implementations/impl_single_node.rs @@ -241,6 +241,7 @@ impl SingleNodeTree { } levels_to_keys.insert(curr.level(), (curr_idx, keys.len())); + // Return tree in sorted order for l in 0..=depth { let &(l, r) = levels_to_keys.get(&l).unwrap(); let subset = &mut keys[l..r]; @@ -419,29 +420,16 @@ impl SingleNodeTree { unmapped } -} - -impl Tree for SingleNodeTree { - type Domain = Domain; - type NodeIndex = MortonKey; - type NodeIndexSlice<'a> = &'a [MortonKey]; - type NodeIndices = MortonKeys; - type Point = Point; - type PointSlice<'a> = &'a [Point]; - type PointData = f64; - type PointDataSlice<'a> = &'a [f64]; - type GlobalIndex = usize; - type GlobalIndexSlice<'a> = &'a [usize]; /// Create a new single-node tree. If non-adaptive (uniform) trees are created, they are specified /// by a user defined maximum depth, if an adaptive tree is created it is specified by only by the /// user defined maximum leaf maximum occupancy n_crit. - fn new( - points: Self::PointDataSlice<'_>, + pub fn new( + points: &[PointType], adaptive: bool, n_crit: Option, depth: Option, - global_idxs: Self::GlobalIndexSlice<'_> + global_idxs: &[usize] ) -> SingleNodeTree { // TODO: Come back and reconcile a runtime point dimension detector @@ -456,6 +444,19 @@ impl Tree for SingleNodeTree { SingleNodeTree::uniform_tree(points, &domain, depth, &global_idxs) } } +} + +impl Tree for SingleNodeTree { + type Domain = Domain; + type NodeIndex = MortonKey; + type NodeIndexSlice<'a> = &'a [MortonKey]; + type NodeIndices = MortonKeys; + type Point = Point; + type PointSlice<'a> = &'a [Point]; + type PointData = f64; + type PointDataSlice<'a> = &'a [f64]; + type GlobalIndex = usize; + type GlobalIndexSlice<'a> = &'a [usize]; fn get_depth(&self) -> u64 { self.depth diff --git a/tree/src/types/multi_node.rs b/tree/src/types/multi_node.rs index 25442f40..5ad8ad35 100644 --- a/tree/src/types/multi_node.rs +++ b/tree/src/types/multi_node.rs @@ -1,47 +1,46 @@ -// //! Data structures and methods to create distributed octrees with MPI. -// use mpi::topology::UserCommunicator; +//! Data structures and methods to create distributed octrees with MPI. +use mpi::topology::UserCommunicator; -// use std::collections::{HashMap, HashSet}; +use std::collections::{HashMap, HashSet}; -// use crate::types::{ -// data::NodeData, -// domain::Domain, -// morton::{KeyType, MortonKey, MortonKeys}, -// point::{Point, Points}, -// }; +use crate::types::{ + domain::Domain, + morton::{KeyType, MortonKey, MortonKeys}, + point::{Point, Points}, +}; -// /// Concrete distributed multi-node tree. -// pub struct MultiNodeTree { -// /// Global communicator for this Tree -// pub world: UserCommunicator, +/// Concrete distributed multi-node tree. +pub struct MultiNodeTree { + /// Global communicator for this Tree + pub world: UserCommunicator, -// /// Adaptivity is optional. -// pub adaptive: bool, + // Depth of the tree + pub depth: u64, -// /// A vector of Cartesian points. -// pub points: Points, + /// Domain spanned by the points. + pub domain: Domain, -// /// All ancestors of leaves in tree, as a set. -// pub keys_set: HashSet, + /// A vector of Cartesian points. + pub points: Points, -// /// The leaf nodes that span the tree, defined by its leaf nodes. -// pub leaves: MortonKeys, + /// The leaves that span the tree. + pub leaves: MortonKeys, -// /// The leaf nodes that span the tree, defined by its leaf nodes, as a set. -// pub leaves_set: HashSet, + /// All nodes in tree. + pub keys: MortonKeys, -// /// Domain spanned by the points in the tree. -// pub domain: Domain, + /// Associate leaves with point indices. + pub leaves_to_points: HashMap, -// /// Map between the points and the leaves in the tree. -// pub points_to_leaves: HashMap, + /// Associate levels with key indices. + pub levels_to_keys: HashMap, -// /// Map between the nodes in the tree and the points they contain. -// pub leaves_to_points: HashMap, + /// All leaves, returned as a set. + pub leaves_set: HashSet, -// // Map between keys and data -// pub keys_to_data: HashMap, + /// All keys, returned as a set. + pub keys_set: HashSet, -// /// Range of Morton keys at this processor, and their current rank [rank, min, max] -// pub range: [KeyType; 3], -// } + /// Range of Morton keys at this processor, and their current rank [rank, min, max] + pub range: [KeyType; 3], +} diff --git a/tree/src/types/point.rs b/tree/src/types/point.rs index 75fb3b25..c2303c91 100644 --- a/tree/src/types/point.rs +++ b/tree/src/types/point.rs @@ -8,7 +8,7 @@ pub type PointType = f64; /// Morton encoding at the lowest possible level of discretization (DEEPEST_LEVEL), and an 'encoded key' /// specifiying its encoding at a given level of discretization. Points also have associated data #[repr(C)] -#[derive(Clone, Debug, Default)] +#[derive(Clone, Debug, Default, Copy)] pub struct Point { /// Physical coordinate in Cartesian space. pub coordinate: [PointType; 3],