Formatting and update handling of neighborhood lists in result

Author: w1th0utnam3

splashsurf_lib/src/dense_subdomains.rs

@@ -623,7 +623,7 @@ pub(crate) fn compute_global_densities_and_neighbors<I: Index, R: Real>(
                             .copied()
                             .zip(particle_densities.iter().copied()),
                     )
-                    // Update density values only for particles inside of the subdomain (ghost particles have wrong values)
+                    // Update density values only for particles inside the subdomain (ghost particles have wrong values)
                     .filter(|(is_inside, _)| *is_inside)
                     .for_each(|(_, (particle_idx, density))| {
                         global_particle_densities[particle_idx] = density;
@@ -659,17 +659,6 @@ pub(crate) fn compute_global_densities_and_neighbors<I: Index, R: Real>(
     let global_particle_densities = global_particle_densities.into_inner();
     let global_neighbors = global_neighbors.into_inner();
 
-    /*
-    {
-        let points = PointCloud3d::new(global_particles.clone());
-        let mesh = MeshWithData::new(points).with_cell_data(MeshAttribute::new_real_scalar(
-            "density".to_string(),
-            global_particle_densities.clone(),
-        ));
-        splashsurf_lib::io::vtk_format::write_vtk(&mesh, "out/new_density.vtk", "particles").unwrap();
-    }
-    */
-
     (global_particle_densities, global_neighbors)
 }
 

splashsurf_lib/src/lib.rs

@@ -213,7 +213,7 @@ pub struct SurfaceReconstruction<I: Index, R: Real> {
 }
 
 impl<I: Index, R: Real> Default for SurfaceReconstruction<I, R> {
-    /// Returns an empty [SurfaceReconstruction] to pass into the inplace surface reconstruction
+    /// Returns an empty [SurfaceReconstruction] to pass into the in-place surface reconstruction
     fn default() -> Self {
         Self {
             grid: UniformGrid::new_zero(),
@@ -227,22 +227,22 @@ impl<I: Index, R: Real> Default for SurfaceReconstruction<I, R> {
 }
 
 impl<I: Index, R: Real> SurfaceReconstruction<I, R> {
-    /// Returns a reference to the actual triangulated surface mesh that is the result of the reconstruction
+    /// Returns a reference to the surface mesh that is the result of the reconstruction
     pub fn mesh(&self) -> &TriMesh3d<R> {
         &self.mesh
     }
 
-    /// Returns a reference to the global particle density vector if it was computed during the reconstruction (always `None` when using independent subdomains with domain decomposition)
+    /// Returns a reference to the global particle density vector if computed during the reconstruction (currently, all reconstruction approaches return this)
     pub fn particle_densities(&self) -> Option<&Vec<R>> {
         self.particle_densities.as_ref()
     }
 
-    /// Returns a reference to the global particle density vector if it was computed during the reconstruction (always `None` when using octree based domain decomposition)
+    /// Returns a reference to the global list of per-particle neighborhood lists if computed during the reconstruction (`None` if not specified in the parameters)
     pub fn particle_neighbors(&self) -> Option<&Vec<Vec<usize>>> {
         self.particle_neighbors.as_ref()
     }
 
-    /// Returns a reference to the virtual background grid that was used as a basis for discretization of the density map for marching cubes, can be used to convert the density map to a hex mesh (using [`density_map::sparse_density_map_to_hex_mesh`])
+    /// Returns a reference to the virtual background grid that was used for marching cubes
     pub fn grid(&self) -> &UniformGrid<I, R> {
         &self.grid
     }
@@ -308,7 +308,7 @@ pub fn reconstruct_surface<I: Index, R: Real>(
     Ok(surface)
 }
 
-/// Performs a marching cubes surface construction of the fluid represented by the given particle positions, inplace
+/// Performs a marching cubes surface construction of the fluid represented by the given particle positions, in-place
 pub fn reconstruct_surface_inplace<I: Index, R: Real>(
     particle_positions: &[Vector3<R>],
     parameters: &Parameters<R>,

splashsurf_lib/src/reconstruction.rs

@@ -27,26 +27,9 @@ pub(crate) fn reconstruct_surface_subdomain_grid<I: Index, R: Real>(
         .global_marching_cubes_grid()
         .context("failed to convert global marching cubes grid")?;
 
-    // Filter "narrow band"
-    /*
-    let narrow_band_particles = extract_narrow_band(&parameters, &particles);
-    let particles = narrow_band_particles;
-     */
-
     let subdomains =
         decomposition::<I, R, GhostMarginClassifier<I>>(&internal_parameters, particle_positions)?;
 
-    /*
-    {
-        use super::dense_subdomains::debug::*;
-        subdomain_stats(&parameters, &particle_positions, &subdomains);
-        info!(
-            "Number of subdomains with only ghost particles: {}",
-            count_no_owned_particles_subdomains(&parameters, &particle_positions, &subdomains)
-        );
-    }
-     */
-
     let (particle_densities, particle_neighbors) = compute_global_densities_and_neighbors(
         &internal_parameters,
         particle_positions,
@@ -70,6 +53,7 @@ pub(crate) fn reconstruct_surface_subdomain_grid<I: Index, R: Real>(
     output_surface.mesh = global_mesh;
     output_surface.particle_densities = Some(particle_densities);
     if parameters.global_neighborhood_list {
+        // The global neighborhood list is only non-empty if requested from the user (requires merging)
         output_surface.particle_neighbors = Some(particle_neighbors);
     }
     Ok(())
@@ -92,15 +76,21 @@ pub(crate) fn reconstruct_surface_global<I: Index, R: Real>(
         .borrow_mut();
 
     // Reuse allocated memory: swap particle densities from output object into the workspace if the former has a larger capacity
-    if let Some(output_densities) = output_surface.particle_densities.as_ref() {
-        if output_densities.capacity() > workspace.particle_densities.capacity() {
-            std::mem::swap(
-                output_surface.particle_densities.as_mut().unwrap(),
-                &mut workspace.particle_densities,
-            );
+    fn swap_output_vec<T>(output_vec: &mut Option<Vec<T>>, workspace_vec: &mut Vec<T>) {
+        if output_vec.as_ref().map(Vec::capacity).unwrap_or(0) > workspace_vec.capacity() {
+            std::mem::swap(output_vec.as_mut().unwrap(), workspace_vec);
         }
     }
 
+    swap_output_vec(
+        &mut output_surface.particle_densities,
+        &mut workspace.particle_densities,
+    );
+    swap_output_vec(
+        &mut output_surface.particle_neighbors,
+        &mut workspace.particle_neighbor_lists,
+    );
+
     // Clear the current mesh, as reconstruction will be appended to output
     output_surface.mesh.clear();
     // Perform global reconstruction without domain decomposition
@@ -113,6 +103,8 @@ pub(crate) fn reconstruct_surface_global<I: Index, R: Real>(
     )?;
 
     output_surface.particle_densities = Some(std::mem::take(&mut workspace.particle_densities));
+    output_surface.particle_neighbors =
+        Some(std::mem::take(&mut workspace.particle_neighbor_lists));
 
     Ok(())
 }