Documentation improvements

particular/README.md

@@ -1,7 +1,7 @@
 # Particular
 
 <div align="center">
-    <img src="https://github.com/Canleskis/particular/blob/main/particular/particular-showcase.gif?raw=true">
+    <img src="https://github.com/Canleskis/particular/blob/main/particular/particular-showcase.gif?raw=true" alt="showcase gif">
 </div>
 
 [![MIT/Apache 2.0](https://img.shields.io/badge/license-MIT%2FApache-blue.svg)](https://github.com/canleskis/particular#license)
@@ -19,21 +19,27 @@ Particular is also built with performance in mind and provides multiple ways of
 
 ### Computation algorithms
 
-There are currently 2 algorithms used by the available compute methods: [BruteForce](https://en.wikipedia.org/wiki/N-body_problem#Simulation) and [BarnesHut](https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation).
+There are currently 2 algorithms used by the available compute methods: [Brute-force](https://en.wikipedia.org/wiki/N-body_problem#Simulation) and [Barnes-Hut](https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation).
 
-Generally speaking, the BruteForce algorithm is more accurate, but slower. The BarnesHut algorithm allows trading accuracy for speed by increasing the `theta` parameter.  
-You can read more about their relative performance [here](#notes-on-performance).
+Generally speaking, the Brute-force algorithm is more accurate, but slower. The Barnes-Hut algorithm allows trading accuracy for speed by increasing the `theta` parameter.  
+You can see more about their relative performance [here](https://particular.rs/benchmarks/).
 
-Particular uses [rayon](https://github.com/rayon-rs/rayon) for parallelization. Enable the "parallel" feature to access the available compute methods.  
-Particular uses [wgpu](https://github.com/gfx-rs/wgpu) for GPU computation. Enable the "gpu" feature to access the available compute methods.
+Particular uses [rayon](https://github.com/rayon-rs/rayon) for parallelization and [wgpu](https://github.com/gfx-rs/wgpu) for GPU computation. Enable the respective `parallel` and `gpu` features to access the available compute methods.
 
 ## Using Particular
 
-### Implementing the [`Particle`] trait
+Particular consists of two "modules", a "backend" that takes care of the abstraction of the computation of the gravitational forces between bodies for different floating-point types and dimensions, and a "frontend" that facilitates usage of that abstraction for user-defined and non-user-defined types. For most simple use cases, the latter is all that you need to know about.
+
+### Simple usage
+
+The [`Particle`] trait provides the main abstraction layer between the internal representation of the position and mass of an object in N-dimensional space and external types, by defining methods to retrieve the position and the gravitational parameter of a type.  
+The return types of these methods, respectively an array of scalars and a scalar, are converted using the [point_mass] method to interface with the backend of Particular.
+
+#### Implementing the [`Particle`] trait
 
 When possible, it can be useful to implement [`Particle`] on a type.
 
-#### Deriving
+##### Deriving
 
 Used when the type has fields named `position` and `mu`:
 
@@ -47,7 +53,7 @@ struct Body {
 }
 ```
 
-#### Manual implementation
+##### Manual implementation
 
 Used when the type does not directly provide a position and a gravitational parameter.
 
@@ -64,7 +70,7 @@ impl Particle for Body {
     fn position(&self) -> [f32; 3] {
         self.position.into()
     }
-    
+
     fn mu(&self) -> f32 {
         self.mass * G
     }
@@ -80,15 +86,26 @@ assert_eq!(particle.position(), [1.0, 1.0, 0.0]);
 assert_eq!(particle.mu(), 5.0);
 ```
 
-### Computing and using the gravitational acceleration
+#### Computing and using the gravitational acceleration
+
+In order to compute the accelerations of your particles, you can use the [accelerations] method on iterators, passing in a mutable reference to a [`ComputeMethod`] of your choice.  
+This method collects the mapped particles in a [`ParticleReordered`] in order to optimise the computation of forces of massless particles, which requires one additional allocation. See the [advanced usage](#advanced-usage) section for information on how to opt out.
 
-In order to compute the accelerations of your particles, you can use the [`accelerations`] method on iterators, passing in a mutable reference to a [`ComputeMethod`] of your choice.
+##### When the iterated type implements [`Particle`]
+
+```rust
+for (acceleration, body) in bodies.iter().accelerations(&mut cm).zip(&mut bodies) {
+    body.velocity += Vec3::from(acceleration) * DT;
+    body.position += body.velocity * DT;
+}
+```
 
-#### When the iterated type doesn't implement [`Particle`]
+##### When the iterated type doesn't implement [`Particle`]
 
 ```rust
 // Items are a tuple of a velocity, a position and a mass.
-// We map them to a tuple of the positions as an array and the mu, since this implements `Particle`.
+// We map them to a tuple of the positions as an array and the mu, 
+// since this implements `Particle`.
 let accelerations = items
     .iter()
     .map(|(_, position, mass)| (*position.as_array(), *mass * G))
@@ -100,25 +117,72 @@ for (acceleration, (velocity, position, _)) in accelerations.zip(&mut items) {
 }
 ```
 
-#### When the iterated type implements [`Particle`]
+### Advanced usage
+
+The "frontend" is built on top of the "backend" but in some instances the abstraction provided by the frontend might not be flexible enough. For example, you might need to access the tree built from the particles for the Barnes-Hut algorithm, want to compute the gravitational forces between two distinct collections of particles, or both at the same time.  
+In that case, you can use the backend directly by calling [compute] on a struct implementing [`ComputeMethod`], passing in an appropriate storage.
+
+#### Storages and built-in [`ComputeMethod`] implementations
+
+Storages are containers that make it easy to apply certain optimisation or algorithms on collections of particles when computing their gravitational acceleration.
+
+The [`ParticleSystem`] storage defines an `affected` slice of particles and a `massive` storage, allowing algorithms to compute gravitational forces the particles in the `massive` storage exert on the `affected` particles. It is used to implement most compute methods, and blanket implementations with the other storages allow a [`ComputeMethod`] implemented with [`ParticleSliceSystem`] or [`ParticleTreeSystem`] to also be implemented with the other storages.
+
+The [`ParticleReordered`] similarly defines a slice of particles, but stores a copy of them in a [`ParticleOrdered`]. These two storages make it easy for algorithms to skip particles with no mass when computing the gravitational forces of particles.
+
+##### Example
 
 ```rust
-for (acceleration, body) in bodies.iter().accelerations(&mut cm).zip(&mut bodies) {
-    body.velocity += Vec3::from(acceleration) * DT;
-    body.position += body.velocity * DT;
+use storage::{ParticleOrdered, ParticleSystem, ParticleTree};
+
+let particles = vec![
+    // ...
+];
+
+// Create a `ParticleOrdered` to split massive and massless particles.
+let ordered = ParticleOrdered::from(&*particles);
+
+// Build a `ParticleTree` from the massive particles.
+let tree = ParticleTree::from(ordered.massive());
+
+// Do something with the tree.
+for (node, data) in std::iter::zip(&tree.get().nodes, &tree.get().data) {
+    // ...
 }
+
+let bh = &mut sequential::BarnesHut { theta: 0.5 };
+// The implementation computes the acceleration exerted on the particles in
+// the `affected` slice.
+// As such, this only computes the acceleration of the massless particles.
+let accelerations = bh.compute(ParticleSystem {
+    affected: ordered.massless(),
+    massive: &tree,
+});
 ```
 
-## Notes on performance
+#### Custom [`ComputeMethod`] implementations
 
-Here is a comparison between 7 available compute methods using an i9 9900KF and an RTX 3080:
+In order to work with the highest number of cases, built-in compute method implementations may not be the most appropriate for your specific use case. You can implement the [`ComputeMethod`] trait on your own type to satisfy your specific requirements but also if you want to implement other algorithms.
 
-<div align="center">
-    <img src="https://github.com/Canleskis/particular/blob/main/particular/particular-comparison.png?raw=true" alt="Performance chart" />
-</div>
+##### Example
+
+```rust
+use storage::{ParticleReordered, ParticleSystem};
+
+struct MyComputeMethod;
 
-Depending on your needs and platform, you may opt for one compute method or another.
-You can also implement the trait on your own type to use other algorithms or combine multiple compute methods and switch between them depending on certain conditions (e.g. the particle count).
+impl ComputeMethod<ParticleReordered<'_, Vec3, f32>> for MyComputeMethod {
+    type Output = Vec<Vec3>;
+
+    fn compute(&mut self, storage: ParticleReordered<Vec3, f32>) -> Self::Output {
+        // Only return the accelerations of the massless particles.
+        sequential::BruteForceScalar.compute(ParticleSliceSystem {
+            affected: storage.massless(),
+            massive: storage.massive(),
+        })
+    }
+}
+```
 
 ## License
 
@@ -128,6 +192,13 @@ This project is licensed under either of [Apache License, Version 2.0](https://g
 
 Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in this project by you, as defined in the Apache 2.0 license, shall be dual licensed as above, without any additional terms or conditions.
 
+[accelerations]: https://docs.rs/particular/latest/particular/particle/trait.Accelerations.html#method.accelerations
+[point_mass]: https://docs.rs/particular/latest/particular/particle/trait.IntoPointMass.html#method.point_mass
+[compute]: https://docs.rs/particular/latest/particular/compute_method/trait.ComputeMethod.html#tymethod.compute
 [`Particle`]: https://docs.rs/particular/latest/particular/particle/trait.Particle.html
 [`ComputeMethod`]: https://docs.rs/particular/latest/particular/compute_method/trait.ComputeMethod.html
-[`accelerations`]: https://docs.rs/particular/latest/particular/particle/trait.Accelerations.html#method.accelerations
+[`ParticleReordered`]: https://docs.rs/particular/latest/particular/compute_method/storage/struct.ParticleReordered.html
+[`ParticleOrdered`]: https://docs.rs/particular/latest/particular/compute_method/storage/struct.ParticleOrdered.html
+[`ParticleSystem`]: https://docs.rs/particular/latest/particular/compute_method/storage/struct.ParticleSystem.html
+[`ParticleSliceSystem`]: https://docs.rs/particular/latest/particular/compute_method/storage/struct.ParticleSliceSystem.html
+[`ParticleTreeSystem`]: https://docs.rs/particular/latest/particular/compute_method/storage/struct.ParticleTreeSystem.html

particular/src/compute_method/gpu_compute/mod.rs

@@ -13,7 +13,7 @@ pub struct WgpuResources {
 }
 
 impl WgpuResources {
-    /// Creates a new [`WgpuResources`] with the given [`wgpu::Device`] and size in bytes of affected and massive particles.
+    /// Creates a new [`WgpuResources`] with the given [`wgpu::Device`] and count of affected and massive particles.
     #[inline]
     pub fn init(device: &wgpu::Device, affected_count: usize, massive_count: usize) -> Self {
         let affected_size = (std::mem::size_of::<PointMass>() * affected_count) as u64;

particular/src/compute_method/mod.rs

@@ -3,8 +3,8 @@
 pub mod gpu_compute;
 /// Trait abstractions for generic vectors and associated floating-point numbers.
 pub mod math;
-/// Representation of the position and mass of an object in N-dimensional space, operations to compute
-/// gravitational forces and the storages used by built-in [`ComputeMethods`](crate::compute_method::ComputeMethod).
+/// Representation of the position and mass of an object in N-dimensional space and collections used by
+/// built-in [`ComputeMethods`](crate::compute_method::ComputeMethod).
 pub mod storage;
 /// Tree, bounding box and BarnesHut implementation details.
 pub mod tree;

particular/src/compute_method/parallel.rs

@@ -77,11 +77,11 @@ where
 
     #[inline]
     fn compute(&mut self, system: ParticleTreeSystem<'_, X, D, V, S>) -> Self::Output {
-        let storage = system.massive;
+        let tree = system.massive;
         system
             .affected
             .par_iter()
-            .map(|p| p.acceleration_tree(storage.tree(), storage.root(), self.theta))
+            .map(|p| p.acceleration_tree(tree.get(), tree.root(), self.theta))
             .collect()
     }
 }

particular/src/compute_method/sequential.rs

@@ -169,11 +169,11 @@ where
 
     #[inline]
     fn compute(&mut self, system: ParticleTreeSystem<X, D, V, S>) -> Self::Output {
-        let storage = system.massive;
+        let tree = system.massive;
         system
             .affected
             .iter()
-            .map(|p| p.acceleration_tree(storage.tree(), storage.root(), self.theta))
+            .map(|p| p.acceleration_tree(tree.get(), tree.root(), self.theta))
             .collect()
     }
 }

particular/src/compute_method/storage.rs

@@ -257,7 +257,7 @@ impl<V, S> PointMass<V, S> {
 pub struct ParticleSystem<'p, V, S, T: ?Sized> {
     /// Particles for which the acceleration is computed.
     pub affected: &'p [PointMass<V, S>],
-    /// Particle storage responsible for the acceleration exerted on the `affected` particles.
+    /// Particles responsible for the acceleration exerted on the `affected` particles, in a storage `S`.
     pub massive: &'p T,
 }
 
@@ -280,7 +280,7 @@ impl<'p, V, S, T: ?Sized> ParticleSystem<'p, V, S, T> {
 /// [`ParticleSystem`] with a slice of particles for the massive storage.
 pub type ParticleSliceSystem<'p, V, S> = ParticleSystem<'p, V, S, [PointMass<V, S>]>;
 
-/// Storage with particles in a [`Orthtree`].
+/// Storage with particles in an [`Orthtree`] and its root.
 #[derive(Clone, Debug)]
 pub struct ParticleTree<const X: usize, const D: usize, V, S> {
     root: Option<NodeID>,
@@ -296,7 +296,7 @@ impl<const X: usize, const D: usize, V, S> ParticleTree<X, D, V, S> {
 
     /// Returns a reference to the [`Orthtree`].
     #[inline]
-    pub const fn tree(&self) -> &Orthtree<X, D, S, PointMass<V, S>> {
+    pub const fn get(&self) -> &Orthtree<X, D, S, PointMass<V, S>> {
         &self.tree
     }
 }
@@ -453,25 +453,6 @@ where
     }
 }
 
-impl<const X: usize, const D: usize, V, S, C, O> ComputeMethod<ParticleSliceSystem<'_, V, S>> for C
-where
-    O: IntoIterator,
-    for<'a> C: ComputeMethod<ParticleTreeSystem<'a, X, D, V, S>, Output = O>,
-    V: Copy + FloatVector<Float = S, Array = [S; D]>,
-    S: Copy + Float + Sum + PartialOrd + FromPrimitive<usize>,
-    BoundingBox<[S; D]>: SubDivide<Division = [BoundingBox<[S; D]>; X]>,
-{
-    type Output = O;
-
-    #[inline]
-    fn compute(&mut self, system: ParticleSliceSystem<V, S>) -> Self::Output {
-        self.compute(ParticleTreeSystem {
-            affected: system.affected,
-            massive: &ParticleTree::from(system.massive),
-        })
-    }
-}
-
 impl<V, S, C, O> ComputeMethod<&[PointMass<V, S>]> for C
 where
     O: IntoIterator,
@@ -519,3 +500,22 @@ where
         })
     }
 }
+
+impl<const X: usize, const D: usize, V, S, C, O> ComputeMethod<ParticleSliceSystem<'_, V, S>> for C
+where
+    O: IntoIterator,
+    for<'a> C: ComputeMethod<ParticleTreeSystem<'a, X, D, V, S>, Output = O>,
+    V: Copy + FloatVector<Float = S, Array = [S; D]>,
+    S: Copy + Float + Sum + PartialOrd + FromPrimitive<usize>,
+    BoundingBox<[S; D]>: SubDivide<Division = [BoundingBox<[S; D]>; X]>,
+{
+    type Output = O;
+
+    #[inline]
+    fn compute(&mut self, system: ParticleSliceSystem<V, S>) -> Self::Output {
+        self.compute(ParticleTreeSystem {
+            affected: system.affected,
+            massive: &ParticleTree::from(system.massive),
+        })
+    }
+}