Update computation of particle rest volume

Author: w1th0utnam3

splashsurf_lib/examples/minimal_sph_levelset.rs

@@ -244,7 +244,7 @@ fn reconstruct() -> Result<(), anyhow::Error> {
     let grid = UniformGrid::<I, R>::from_aabb(&particle_aabb, cube_size)?;
     let kernel = kernel::CubicSplineKernel::new(compact_support_radius);
 
-    let particle_rest_volume = (2.0 * args.particle_radius).powi(3);
+    let particle_rest_volume = kernel::Volume::cube_particle(args.particle_radius);
     let particle_rest_mass = particle_rest_volume * args.rest_density;
 
     info!(

splashsurf_lib/src/dense_subdomains.rs

@@ -22,8 +22,8 @@ use crate::neighborhood_search::{
 use crate::topology::Direction;
 use crate::uniform_grid::{EdgeIndex, GridConstructionError, UniformCartesianCubeGrid3d};
 use crate::{
-    Aabb3d, MapType, Parameters, RealConvert, SpatialDecomposition, SurfaceReconstruction, new_map,
-    new_parallel_map, profile,
+    Aabb3d, MapType, Parameters, RealConvert, SpatialDecomposition, SurfaceReconstruction, kernel,
+    new_map, new_parallel_map, profile,
 };
 use crate::{Index, Real};
 
@@ -110,7 +110,7 @@ pub(crate) fn initialize_parameters<I: Index, R: Real>(
     let cube_size = parameters.cube_size;
     let surface_threshold = parameters.iso_surface_threshold;
 
-    let particle_rest_volume = R::from_float(4.0) * R::frac_pi_3() * particle_radius.powi(3);
+    let particle_rest_volume = kernel::Volume::cube_particle(particle_radius);
     let particle_rest_mass = particle_rest_volume * particle_rest_density;
 
     let ghost_particle_margin =

splashsurf_lib/src/kernel.rs

@@ -6,6 +6,21 @@ use numeric_literals::replace_float_literals;
 
 // TODO: Add reference for the kernel function, document formula
 
+/// Utility functions for computing the volume of fluid particles
+pub struct Volume;
+
+impl Volume {
+    /// Returns the volume of a particle representing a cube of fluid
+    pub fn cube_particle<R: Real>(particle_radius: R) -> R {
+        (particle_radius + particle_radius).powi(3)
+    }
+
+    /// Returns the volume of a particle representing a sphere of fluid
+    pub fn sphere_particle<R: Real>(particle_radius: R) -> R {
+        R::from_float(4.0) * R::frac_pi_3() * particle_radius.powi(3)
+    }
+}
+
 /// Trait for symmetric kernel functions in three dimensions
 pub trait SymmetricKernel3d<R: Real> {
     /// Evaluates the kernel at the radial distance `r` relative to the origin

splashsurf_lib/src/reconstruction.rs

@@ -6,7 +6,7 @@ use crate::mesh::TriMesh3d;
 use crate::uniform_grid::UniformGrid;
 use crate::workspace::LocalReconstructionWorkspace;
 use crate::{
-    Index, Parameters, Real, ReconstructionError, SurfaceReconstruction, density_map,
+    Index, Parameters, Real, ReconstructionError, SurfaceReconstruction, density_map, kernel,
     marching_cubes, neighborhood_search, profile,
 };
 use anyhow::Context;
@@ -120,8 +120,7 @@ pub(crate) fn compute_particle_densities_and_neighbors<I: Index, R: Real>(
     profile!("compute_particle_densities_and_neighbors");
 
     let particle_rest_density = parameters.rest_density;
-    let particle_rest_volume =
-        R::from_float(4.0) * R::frac_pi_3() * parameters.particle_radius.powi(3);
+    let particle_rest_volume = kernel::Volume::cube_particle(parameters.particle_radius);
     let particle_rest_mass = particle_rest_volume * particle_rest_density;
 
     trace!("Starting neighborhood search...");
@@ -153,8 +152,7 @@ pub(crate) fn reconstruct_single_surface_append<I: Index, R: Real>(
     output_mesh: &mut TriMesh3d<R>,
 ) -> Result<(), ReconstructionError<I, R>> {
     let particle_rest_density = parameters.rest_density;
-    let particle_rest_volume =
-        R::from_float(4.0) * R::frac_pi_3() * parameters.particle_radius.powi(3);
+    let particle_rest_volume = kernel::Volume::cube_particle(parameters.particle_radius);
     let particle_rest_mass = particle_rest_volume * particle_rest_density;
 
     let particle_densities = {