From 4c0ea2772a74df25eaf2b12ee286746caafae371 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20L=C3=B6schner?= Date: Wed, 26 Jul 2023 15:32:21 +0200 Subject: [PATCH] Documentation, small refactorings --- splashsurf/src/main.rs | 5 + splashsurf_lib/src/dense_subdomains.rs | 4 +- splashsurf_lib/src/density_map.rs | 44 ++-- splashsurf_lib/src/lib.rs | 4 +- splashsurf_lib/src/neighborhood_search.rs | 234 ++++++---------------- 5 files changed, 86 insertions(+), 205 deletions(-) diff --git a/splashsurf/src/main.rs b/splashsurf/src/main.rs index 17e5636..c38b4c3 100644 --- a/splashsurf/src/main.rs +++ b/splashsurf/src/main.rs @@ -1,3 +1,8 @@ +//! The `splashsurf` surface reconstruction CLI. +//! +//! For documentation of the CLI see the [README](https://github.com/InteractiveComputerGraphics/splashsurf) in the project repository. +//! The reconstruction procedure and other internals of the CLI are provided by the [`splashsurf_lib`] crate. + mod convert; mod io; mod reconstruction; diff --git a/splashsurf_lib/src/dense_subdomains.rs b/splashsurf_lib/src/dense_subdomains.rs index ed6ab17..5ce3c7d 100644 --- a/splashsurf_lib/src/dense_subdomains.rs +++ b/splashsurf_lib/src/dense_subdomains.rs @@ -567,7 +567,7 @@ pub(crate) fn compute_global_density_vector( &subdomain_particles, parameters.compact_support_radius, neighborhood_lists, - is_inside, + |i| is_inside[i], ); sequential_compute_particle_densities_filtered::( @@ -576,7 +576,7 @@ pub(crate) fn compute_global_density_vector( parameters.compact_support_radius, parameters.particle_rest_mass, particle_densities, - is_inside, + |i| is_inside[i], ); // Write particle densities into global storage diff --git a/splashsurf_lib/src/density_map.rs b/splashsurf_lib/src/density_map.rs index d96e3b4..b136b07 100644 --- a/splashsurf_lib/src/density_map.rs +++ b/splashsurf_lib/src/density_map.rs @@ -72,7 +72,7 @@ pub fn compute_particle_densities( &mut densities, ) } else { - sequential_compute_particle_densities::( + sequential_compute_particle_densities::( particle_positions, particle_neighbor_lists, compact_support_radius, @@ -102,7 +102,7 @@ pub fn compute_particle_densities_inplace( densities, ) } else { - sequential_compute_particle_densities::( + sequential_compute_particle_densities::( particle_positions, particle_neighbor_lists, compact_support_radius, @@ -120,43 +120,37 @@ fn init_density_storage(densities: &mut Vec, new_len: usize) { /// Computes the individual densities of particles using a standard SPH sum, sequential implementation #[inline(never)] -pub fn sequential_compute_particle_densities( +pub fn sequential_compute_particle_densities( particle_positions: &[Vector3], - particle_neighbor_lists: &[Vec], + particle_neighbor_lists: &Nl, compact_support_radius: R, particle_rest_mass: R, particle_densities: &mut Vec, ) { profile!("sequential_compute_particle_densities"); - init_density_storage(particle_densities, particle_positions.len()); - - // Pre-compute the kernel which can be queried using squared distances - let kernel = DiscreteSquaredDistanceCubicKernel::new::(1000, compact_support_radius); - - for (i, (particle_i_position, particle_i_neighbors)) in particle_positions - .iter() - .zip(particle_neighbor_lists.iter()) - .enumerate() - { - let mut particle_i_density = kernel.evaluate(R::zero()); - for particle_j_position in particle_i_neighbors.iter().map(|&j| &particle_positions[j]) { - let r_squared = (particle_j_position - particle_i_position).norm_squared(); - particle_i_density += kernel.evaluate(r_squared); - } - particle_i_density *= particle_rest_mass; - particle_densities[i] = particle_i_density; - } + sequential_compute_particle_densities_filtered::( + particle_positions, + particle_neighbor_lists, + compact_support_radius, + particle_rest_mass, + particle_densities, + |_| true, + ) } #[inline(never)] -pub fn sequential_compute_particle_densities_filtered( +pub fn sequential_compute_particle_densities_filtered< + I: Index, + R: Real, + Nl: NeighborhoodList + ?Sized, +>( particle_positions: &[Vector3], particle_neighbor_lists: &Nl, compact_support_radius: R, particle_rest_mass: R, particle_densities: &mut Vec, - filter: &[bool], + filter: impl Fn(usize) -> bool, ) { profile!("sequential_compute_particle_densities_filtered"); @@ -168,7 +162,7 @@ pub fn sequential_compute_particle_densities_filtered>, new_len: neighborhood_list.resize_with(new_len, || Vec::with_capacity(15)); } -/// Performs a neighborhood search, returning the indices of all neighboring particles in the given search radius per particle, sequential implementation +/// Performs a neighborhood search (sequential implementation) +/// +/// Returns the indices of all neighboring particles in the given search radius per particle as a `Vec>`. #[inline(never)] pub fn neighborhood_search_spatial_hashing( domain: &Aabb3d, @@ -128,9 +130,30 @@ pub fn neighborhood_search_spatial_hashing( search_radius: R, neighborhood_list: &mut Vec>, ) { - // TODO: Use ArrayStorage from femproto instead of Vec of Vecs? + neighborhood_search_spatial_hashing_filtered::( + domain, + particle_positions, + search_radius, + neighborhood_list, + |_| true, + ) +} + +/// Performs a neighborhood search (sequential implementation, with filter) +/// +/// Returns the indices of all neighboring particles in the given search radius per particle as a `Vec>`. +/// The filter specifies which particles the neighbor lists should be computed for (`true`: compute neighbors). +/// Note that the particles that were filtered out will still appear in the neighbor lists of the particles that were not filtered out. +#[inline(never)] +pub fn neighborhood_search_spatial_hashing_filtered( + domain: &Aabb3d, + particle_positions: &[Vector3], + search_radius: R, + neighborhood_list: &mut Vec>, + filter: impl Fn(usize) -> bool, +) { // FIXME: Replace unwraps? - profile!("neighborhood_search_spatial_hashing"); + profile!("neighborhood_search_spatial_hashing_filtered"); assert!( search_radius > R::zero(), @@ -182,6 +205,10 @@ pub fn neighborhood_search_spatial_hashing( // Iterate over all particles of the current cell for &particle_i in particles { + if !filter(particle_i) { + continue; + } + let pos_i = &particle_positions[particle_i]; let particle_i_neighbors = &mut neighborhood_list[particle_i]; @@ -202,6 +229,7 @@ pub fn neighborhood_search_spatial_hashing( } } +/// Stores particle neighborhood lists contiguously in memory using a second offset array pub struct FlatNeighborhoodList { /// Offsets to the start of the neighborhood list of the given particle (very last entry contains total number of neighbor particles) pub neighbor_ptr: Vec, @@ -224,6 +252,7 @@ impl FlatNeighborhoodList { self.neighbor_ptr.len() - 1 } + /// Returns a slice containing the neighborhood list of the given particle pub fn get_neighbors(&self, particle_i: usize) -> Option<&[usize]> { let range = self .neighbor_ptr @@ -233,6 +262,7 @@ impl FlatNeighborhoodList { range.and_then(|(start, end)| self.neighbors.get(start..end)) } + /// Returns a mutable slice containing the neighborhood list of the given particle pub fn get_neighbors_mut(&mut self, particle_i: usize) -> Option<&mut [usize]> { let range = self .neighbor_ptr @@ -272,16 +302,20 @@ impl NeighborhoodList for FlatNeighborhoodList { } } -impl NeighborhoodList for Vec> { +/// Implementation for `Vec>` and `[Vec]` +impl]> + ?Sized> NeighborhoodList for T { fn len(&self) -> usize { - self.len() + self.as_ref().len() } fn neighbors(&self, particle_i: usize) -> &[usize] { - &self[particle_i] + &self.as_ref()[particle_i] } } +/// Performs a neighborhood search (sequential implementation, returning a [`FlatNeighborhoodList`]) +/// +/// Returns the indices of all neighboring particles in the given search radius per particle as a [`FlatNeighborhoodList`]. pub fn neighborhood_search_spatial_hashing_flat( domain: &Aabb3d, particle_positions: &[Vector3], @@ -290,90 +324,26 @@ pub fn neighborhood_search_spatial_hashing_flat( ) { profile!("neighborhood_search_spatial_hashing_flat"); - assert!( - search_radius > R::zero(), - "Search radius for neighborhood search has to be positive!" - ); - assert!( - domain.is_consistent(), - "Domain for neighborhood search has to be consistent!" - ); - assert!( - !domain.is_degenerate(), - "Domain for neighborhood search cannot be degenerate!" - ); - - let search_radius_squared = search_radius * search_radius; - - // Create a new grid for neighborhood search - let grid = UniformGrid::from_aabb(&domain, search_radius) - .expect("Failed to construct grid for neighborhood search!"); - // Map for spatially hashed storage of all particles (map from cell -> enclosed particles) - let particles_per_cell = - sequential_generate_cell_to_particle_map::(&grid, particle_positions); - - { - neighborhood_list.neighbor_ptr.clear(); - neighborhood_list - .neighbor_ptr - .resize(particle_positions.len() + 1, 0); - neighborhood_list.neighbors.clear(); - } - - { - profile!("write particle neighbors"); - let mut counter = 0; - for (particle_i, pos_i) in particle_positions.iter().enumerate() { - // Store start of current particle neighbor list - neighborhood_list.neighbor_ptr[particle_i] = counter; - - // Cell of the current particle - let current_cell = grid.get_cell(grid.enclosing_cell(pos_i)).unwrap(); - - let mut neighbor_test = |particle_j| { - let pos_j: &Vector3 = &particle_positions[particle_j]; - if (pos_j - pos_i).norm_squared() < search_radius_squared { - // A neighbor was found - neighborhood_list.neighbors.push(particle_j); - counter += 1; - } - }; - - // Loop over adjacent cells - grid.cells_adjacent_to_cell(¤t_cell) - .filter_map(|c| { - let flat_cell_index = grid.flatten_cell_index(&c); - particles_per_cell.get(&flat_cell_index) - }) - .flatten() - .copied() - .for_each(|particle_j| neighbor_test(particle_j)); - - // Loop over current cell - std::iter::once(current_cell) - .filter_map(|c| { - let flat_cell_index = grid.flatten_cell_index(&c); - particles_per_cell.get(&flat_cell_index) - }) - .flatten() - .copied() - .for_each(|particle_j| { - if particle_i != particle_j { - neighbor_test(particle_j); - } - }); - } - // Store end of the last neighbor list - *neighborhood_list.neighbor_ptr.last_mut().unwrap() = counter; - } + neighborhood_search_spatial_hashing_flat_filtered::( + domain, + particle_positions, + search_radius, + neighborhood_list, + |_| true, + ) } +/// Performs a neighborhood search (sequential implementation, with filter, returning a [`FlatNeighborhoodList`]) +/// +/// Returns the indices of all neighboring particles in the given search radius per particle as a [`FlatNeighborhoodList`]. +/// The filter specifies which particles the neighbor lists should be computed for (`true`: compute neighbors). +/// Note that the particles that were filtered out will still appear in the neighbor lists of the particles that were not filtered out. pub fn neighborhood_search_spatial_hashing_flat_filtered( domain: &Aabb3d, particle_positions: &[Vector3], search_radius: R, neighborhood_list: &mut FlatNeighborhoodList, - filter: &[bool], + filter: impl Fn(usize) -> bool, ) { profile!("neighborhood_search_spatial_hashing_flat_filtered"); @@ -410,15 +380,11 @@ pub fn neighborhood_search_spatial_hashing_flat_filtered( { profile!("write particle neighbors"); let mut counter = 0; - for (particle_i, (pos_i, filter)) in particle_positions - .iter() - .zip(filter.iter().copied()) - .enumerate() - { + for (particle_i, pos_i) in particle_positions.iter().enumerate() { // Store start of current particle neighbor list neighborhood_list.neighbor_ptr[particle_i] = counter; - if !filter { + if !filter(particle_i) { continue; } @@ -463,93 +429,9 @@ pub fn neighborhood_search_spatial_hashing_flat_filtered( } } -#[inline(never)] -pub fn neighborhood_search_spatial_hashing_filtered( - domain: &Aabb3d, - particle_positions: &[Vector3], - search_radius: R, - neighborhood_list: &mut Vec>, - filter: &[bool], -) { - // TODO: Use ArrayStorage from femproto instead of Vec of Vecs? - // FIXME: Replace unwraps? - profile!("neighborhood_search_spatial_hashing_filtered"); - - assert!( - search_radius > R::zero(), - "Search radius for neighborhood search has to be positive!" - ); - assert!( - domain.is_consistent(), - "Domain for neighborhood search has to be consistent!" - ); - assert!( - !domain.is_degenerate(), - "Domain for neighborhood search cannot be degenerate!" - ); - - let search_radius_squared = search_radius * search_radius; - - // Create a new grid for neighborhood search - let grid = UniformGrid::from_aabb(&domain, search_radius) - .expect("Failed to construct grid for neighborhood search!"); - // Map for spatially hashed storage of all particles (map from cell -> enclosed particles) - let particles_per_cell = - sequential_generate_cell_to_particle_map::(&grid, particle_positions); - - // Build neighborhood lists cell by cell - init_neighborhood_list(neighborhood_list, particle_positions.len()); - { - profile!("calculate_particle_neighbors_seq"); - let mut potential_neighbor_particle_vecs = Vec::new(); - for (&flat_cell_index, particles) in &particles_per_cell { - let current_cell = grid.try_unflatten_cell_index(flat_cell_index).unwrap(); - - // Collect references to the particle lists of all existing adjacent cells and the cell itself - potential_neighbor_particle_vecs.clear(); - potential_neighbor_particle_vecs.extend( - grid.cells_adjacent_to_cell(¤t_cell) - .chain(std::iter::once(current_cell)) - .filter_map(|c| { - let flat_cell_index = grid.flatten_cell_index(&c); - particles_per_cell.get(&flat_cell_index) - }), - ); - - // Returns an iterator over all particles of all adjacent cells and the cell itself - let potential_neighbor_particle_iter = || { - potential_neighbor_particle_vecs - .iter() - .flat_map(|v| v.iter()) - }; - - // Iterate over all particles of the current cell - for &particle_i in particles { - if !filter[particle_i] { - continue; - } - - let pos_i = &particle_positions[particle_i]; - let particle_i_neighbors = &mut neighborhood_list[particle_i]; - - // Check for neighborhood with all neighboring cells - for &particle_j in potential_neighbor_particle_iter() { - if particle_i == particle_j { - continue; - } - - let pos_j = &particle_positions[particle_j]; - if (pos_j - pos_i).norm_squared() < search_radius_squared { - // A neighbor was found - particle_i_neighbors.push(particle_j); - } - } - } - } - } -} - -/// Performs a neighborhood search, returning the indices of all neighboring particles in the given search radius per particle, multi-threaded implementation +/// Performs a neighborhood search (multi-threaded implementation) +/// +/// Returns the indices of all neighboring particles in the given search radius per particle as a `Vec>`. #[inline(never)] pub fn neighborhood_search_spatial_hashing_parallel( domain: &Aabb3d,