Lots of cleanup

Cargo.toml

@@ -5,7 +5,7 @@ battleship = []
 
 [package]
 name = "bempp-octree"
-version = "0.0.1-dev"
+version = "0.1.0"
 edition = "2021"
 authors = [
     "Srinath Kailasa <[email protected]>, Timo Betcke <[email protected]>",
@@ -27,15 +27,13 @@ itertools = "0.13.*"
 rand = { version = "0.8.5", features = ["alloc"] }
 rand_chacha = "0.3.*"
 num = "0.4.*"
-vtkio = "0.6.*"
 mpi = { version = "0.8.*", features = ["derive", "user-operations"] }
 
 [profile.release]
 debug = 1
 
 [dev-dependencies]
 rand_distr = "0.4.3"
-#criterion = { version = "0.5.*", features = ["html_reports"]}
 
 [build-dependencies]
 cbindgen = "0.27.0"

examples/mpi_complete_tree.rs

@@ -1,11 +1,8 @@
-//! Test the computation of a complete octree.
+//! Demonstrate the instantiation of a complete octree using MPI.
 
-use bempp_octree::{
-    morton::MortonKey,
-    octree::{is_complete_linear_and_balanced, KeyType, Octree},
-    tools::{gather_to_all, generate_random_points},
-};
-use itertools::Itertools;
+use std::time::Instant;
+
+use bempp_octree::{generate_random_points, Octree};
 use mpi::traits::Communicator;
 use rand::prelude::*;
 use rand_chacha::ChaCha8Rng;
@@ -21,9 +18,9 @@ pub fn main() {
     let mut rng = ChaCha8Rng::seed_from_u64(comm.rank() as u64);
 
     // Create `npoints` per rank.
-    let npoints = 10000;
+    let npoints = 1000000;
 
-    // Generate random points.
+    // Generate random points on the positive octant of the unit sphere.
 
     let mut points = generate_random_points(npoints, &mut rng, &comm);
     // Make sure that the points live on the unit sphere.
@@ -37,101 +34,21 @@ pub fn main() {
         point.coords_mut()[2] /= len;
     }
 
-    let tree = Octree::new(&points, 15, 50, &comm);
-
-    // We now check that each node of the tree has all its neighbors available.
-
-    let leaf_tree = tree.leaf_keys();
-    let all_keys = tree.all_keys();
-
-    assert!(is_complete_linear_and_balanced(leaf_tree, &comm));
-    for &key in leaf_tree {
-        // We only check interior keys. Leaf keys may not have a neighbor
-        // on the same level.
-        let mut parent = key.parent();
-        while parent.level() > 0 {
-            // Check that the key itself is there.
-            assert!(all_keys.contains_key(&key));
-            // Check that all its neighbours are there.
-            for neighbor in parent.neighbours().iter().filter(|&key| key.is_valid()) {
-                assert!(all_keys.contains_key(neighbor));
-            }
-            parent = parent.parent();
-            // Check that the parent is there.
-            assert!(all_keys.contains_key(&parent));
-        }
-    }
-
-    // At the end check that the root of the tree is also contained.
-    assert!(all_keys.contains_key(&MortonKey::root()));
-
-    // Count the number of ghosts on each rank
-    // Count the number of global keys on each rank.
-
-    // Assert that all ghosts are from a different rank and count them.
-
-    let nghosts = all_keys
-        .iter()
-        .filter_map(|(_, &value)| {
-            if let KeyType::Ghost(rank) = value {
-                assert!(rank != comm.size() as usize);
-                Some(rank)
-            } else {
-                None
-            }
-        })
-        .count();
-
-    if comm.size() == 1 {
-        assert_eq!(nghosts, 0);
-    } else {
-        assert!(nghosts > 0);
-    }
+    let start = Instant::now();
+    let octree = Octree::new(&points, 15, 50, &comm);
+    let duration = start.elapsed();
 
-    let nglobal = all_keys
-        .iter()
-        .filter(|(_, &value)| matches!(value, KeyType::Global))
-        .count();
+    let global_number_of_points = octree.global_number_of_points();
+    let global_max_level = octree.global_max_level();
 
-    // Assert that all globals across all ranks have the same count.
-
-    let nglobals = gather_to_all(std::slice::from_ref(&nglobal), &comm);
-
-    assert_eq!(nglobals.iter().unique().count(), 1);
-
-    // Check that the points are associated with the correct leaf keys.
-    let mut npoints = 0;
-    let leaf_point_map = tree.leaf_keys_to_local_point_indices();
-
-    for (key, point_indices) in leaf_point_map {
-        for &index in point_indices {
-            assert!(key.is_ancestor(tree.point_keys()[index]));
-        }
-        npoints += point_indices.len();
-    }
-
-    // Make sure that the number of points and point keys lines up
-    // with the points stored for each leaf key.
-    assert_eq!(npoints, tree.points().len());
-    assert_eq!(npoints, tree.point_keys().len());
-
-    // Check the neighbour relationships.
-
-    let all_neighbours = tree.neighbour_map();
-    let all_keys = tree.all_keys();
-
-    for (key, key_type) in all_keys {
-        // Ghost keys should not be in the neighbour map.
-        match key_type {
-            KeyType::Ghost(_) => assert!(!all_neighbours.contains_key(key)),
-            _ => {
-                // If it is not a ghost the key should be in the neighbour map.
-                assert!(all_neighbours.contains_key(key));
-            }
-        }
-    }
+    // We now check that each node of the tree has all its neighbors available.
 
     if comm.rank() == 0 {
-        println!("No errors were found in setting up tree.");
+        println!(
+            "Setup octree with {} points and maximum level {} in {} ms",
+            global_number_of_points,
+            global_max_level,
+            duration.as_millis()
+        );
     }
 }

examples/mpi_complete_tree_debug.rs

@@ -0,0 +1,137 @@
+//! Test the computation of a complete octree.
+
+use bempp_octree::{
+    morton::MortonKey,
+    octree::{is_complete_linear_and_balanced, KeyType, Octree},
+    tools::{gather_to_all, generate_random_points},
+};
+use itertools::Itertools;
+use mpi::traits::Communicator;
+use rand::prelude::*;
+use rand_chacha::ChaCha8Rng;
+
+pub fn main() {
+    // Initialise MPI
+    let universe = mpi::initialize().unwrap();
+
+    // Get the world communicator
+    let comm = universe.world();
+
+    // Initialise a seeded Rng.
+    let mut rng = ChaCha8Rng::seed_from_u64(comm.rank() as u64);
+
+    // Create `npoints` per rank.
+    let npoints = 10000;
+
+    // Generate random points.
+
+    let mut points = generate_random_points(npoints, &mut rng, &comm);
+    // Make sure that the points live on the unit sphere.
+    for point in points.iter_mut() {
+        let len = point.coords()[0] * point.coords()[0]
+            + point.coords()[1] * point.coords()[1]
+            + point.coords()[2] * point.coords()[2];
+        let len = len.sqrt();
+        point.coords_mut()[0] /= len;
+        point.coords_mut()[1] /= len;
+        point.coords_mut()[2] /= len;
+    }
+
+    let tree = Octree::new(&points, 15, 50, &comm);
+
+    // We now check that each node of the tree has all its neighbors available.
+
+    let leaf_tree = tree.leaf_keys();
+    let all_keys = tree.all_keys();
+
+    assert!(is_complete_linear_and_balanced(leaf_tree, &comm));
+    for &key in leaf_tree {
+        // We only check interior keys. Leaf keys may not have a neighbor
+        // on the same level.
+        let mut parent = key.parent();
+        while parent.level() > 0 {
+            // Check that the key itself is there.
+            assert!(all_keys.contains_key(&key));
+            // Check that all its neighbours are there.
+            for neighbor in parent.neighbours().iter().filter(|&key| key.is_valid()) {
+                assert!(all_keys.contains_key(neighbor));
+            }
+            parent = parent.parent();
+            // Check that the parent is there.
+            assert!(all_keys.contains_key(&parent));
+        }
+    }
+
+    // At the end check that the root of the tree is also contained.
+    assert!(all_keys.contains_key(&MortonKey::root()));
+
+    // Count the number of ghosts on each rank
+    // Count the number of global keys on each rank.
+
+    // Assert that all ghosts are from a different rank and count them.
+
+    let nghosts = all_keys
+        .iter()
+        .filter_map(|(_, &value)| {
+            if let KeyType::Ghost(rank) = value {
+                assert!(rank != comm.size() as usize);
+                Some(rank)
+            } else {
+                None
+            }
+        })
+        .count();
+
+    if comm.size() == 1 {
+        assert_eq!(nghosts, 0);
+    } else {
+        assert!(nghosts > 0);
+    }
+
+    let nglobal = all_keys
+        .iter()
+        .filter(|(_, &value)| matches!(value, KeyType::Global))
+        .count();
+
+    // Assert that all globals across all ranks have the same count.
+
+    let nglobals = gather_to_all(std::slice::from_ref(&nglobal), &comm);
+
+    assert_eq!(nglobals.iter().unique().count(), 1);
+
+    // Check that the points are associated with the correct leaf keys.
+    let mut npoints = 0;
+    let leaf_point_map = tree.leaf_keys_to_local_point_indices();
+
+    for (key, point_indices) in leaf_point_map {
+        for &index in point_indices {
+            assert!(key.is_ancestor(tree.point_keys()[index]));
+        }
+        npoints += point_indices.len();
+    }
+
+    // Make sure that the number of points and point keys lines up
+    // with the points stored for each leaf key.
+    assert_eq!(npoints, tree.points().len());
+    assert_eq!(npoints, tree.point_keys().len());
+
+    // Check the neighbour relationships.
+
+    let all_neighbours = tree.neighbour_map();
+    let all_keys = tree.all_keys();
+
+    for (key, key_type) in all_keys {
+        // Ghost keys should not be in the neighbour map.
+        match key_type {
+            KeyType::Ghost(_) => assert!(!all_neighbours.contains_key(key)),
+            _ => {
+                // If it is not a ghost the key should be in the neighbour map.
+                assert!(all_neighbours.contains_key(key));
+            }
+        }
+    }
+
+    if comm.rank() == 0 {
+        println!("No errors were found in setting up tree.");
+    }
+}

src/geometry.rs

@@ -4,7 +4,11 @@ use mpi::traits::Equivalence;
 
 use crate::constants::DEEPEST_LEVEL;
 
-/// Definition of a point.
+/// Definition of a point in 3d space.
+///
+/// A point consists of three coordinates and a global id.
+/// The global id makes it easier to identify points as they
+/// are distributed across MPI nodes.
 #[derive(Clone, Copy, Equivalence)]
 pub struct Point {
     coords: [f64; 3],
@@ -47,6 +51,9 @@ impl PhysicalBox {
     }
 
     /// Give a slice of points. Compute an associated bounding box.
+    ///
+    /// The created bounding box is slightly bigger than the actual point set.
+    /// This is to make sure that no point lies exactly on the boundary of the box.
     pub fn from_points(points: &[Point]) -> PhysicalBox {
         let mut xmin = f64::MAX;
         let mut xmax = f64::MIN;

src/lib.rs

@@ -70,6 +70,15 @@
 //! - [Ghost(rank)](crate::octree::KeyType::Ghost): A ghost key that is adjacent to the local node and originates on the process with rank `rank`.
 //! - [Global](crate::octree::KeyType::Global): A global key that is stored on all processes. These are the ancestors of the coarse tree leafs.
 //!
+//! An important property in octrees is the neighbour relationship. To get a hash map that stores the neighbours for each non-ghost key
+//! use
+//! ```ignore
+//! let neighbour_map = octree.neighbour_map();
+//! ```
+//! The `neighbour_map` stores for each non-ghost key a vector of keys that are adjacent to the neighbour. Neighbours need not
+//! necesssarily be at the same level in the tree. Only for interior keys neighbours are guaranteed to be at the same level. For leaf keys
+//! neighbours may be one level higher or lower.
+//!
 //! Each key in this library is a [Morton Key](MortonKey). For details of Morton keys see the description in the [morton] module.
 #![cfg_attr(feature = "strict", deny(warnings), deny(unused_crate_dependencies))]
 #![warn(missing_docs)]
@@ -82,6 +91,7 @@ pub mod parsort;
 pub mod tools;
 pub mod types;
 
+pub use crate::geometry::{PhysicalBox, Point};
 pub use crate::octree::Octree;
 pub use crate::tools::generate_random_points;
 pub use morton::MortonKey;

src/octree.rs

@@ -237,6 +237,42 @@ impl<'o, C: CommunicatorCollectives> Octree<'o, C> {
     pub fn neighbour_map(&self) -> &HashMap<MortonKey, Vec<MortonKey>> {
         &self.neighbours
     }
+
+    /// Return the local number of points in the octree.
+    pub fn local_number_of_points(&self) -> usize {
+        self.points.len()
+    }
+
+    /// Return the global number of points in the octree.
+    pub fn global_number_of_points(&self) -> usize {
+        let mut global_num_points = 0;
+        self.comm.all_reduce_into(
+            &self.local_number_of_points(),
+            &mut global_num_points,
+            SystemOperation::sum(),
+        );
+        global_num_points
+    }
+
+    /// Return the local maximum level
+    pub fn local_max_level(&self) -> usize {
+        self.leaf_keys
+            .iter()
+            .map(|key| key.level())
+            .max()
+            .unwrap_or(0)
+    }
+
+    /// Return the global maximum level
+    pub fn global_max_level(&self) -> usize {
+        let mut global_max_level = 0;
+        self.comm.all_reduce_into(
+            &self.local_max_level(),
+            &mut global_max_level,
+            SystemOperation::max(),
+        );
+        global_max_level
+    }
 }
 
 /// Test if an array of keys are the leafs of a complete linear and balanced tree.