diff --git a/splashsurf_lib/src/density_map.rs b/splashsurf_lib/src/density_map.rs index 0b8802b..763c3b6 100644 --- a/splashsurf_lib/src/density_map.rs +++ b/splashsurf_lib/src/density_map.rs @@ -528,7 +528,7 @@ pub(crate) fn compute_kernel_evaluation_radius( // The number of points corresponding to the number of supported cells let supported_points: I = I::one() + supported_cells; - // Evaluate kernel in a smaller domain, points outside of this radius have to be assumed to be outside of the iso-surface + // Evaluate kernel in a smaller domain, points outside of this radius have to be assumed to be outside the iso-surface let kernel_evaluation_radius = cube_size * half_supported_cells_real * (R::one() + R::default_epsilon().sqrt()); @@ -599,7 +599,7 @@ impl SparseDensityMapGenerator { particle: &Vector3, particle_density: R, ) { - // Skip particles outside of allowed domain + // Skip particles outside allowed domain if !self.allowed_domain.contains_point(particle) { return; } diff --git a/splashsurf_lib/src/lib.rs b/splashsurf_lib/src/lib.rs index 375db6e..f803fd3 100644 --- a/splashsurf_lib/src/lib.rs +++ b/splashsurf_lib/src/lib.rs @@ -412,7 +412,7 @@ pub fn grid_for_reconstruction( } else { Aabb3d::from_points(particle_positions) }; - // TODO: Is this really necessary? + // TODO: Is this really necessary? This seems unnecessary purely for the density map... aabb.grow_uniformly(particle_radius); aabb }; diff --git a/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs b/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs index 863da93..9057aad 100644 --- a/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs +++ b/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs @@ -45,7 +45,7 @@ pub(crate) fn construct_mc_input( /// Note: The threshold flags in the resulting cell data are not complete and still have to be updated after /// this procedure using the [update_cell_data_threshold_flags] function. /// -/// Note: This functions assumes that the default value for missing point data is below the iso-surface threshold. +/// Note: This functions assumes that the default value for missing point data is zero. #[inline(never)] fn interpolate_points_to_cell_data_generic( grid: &UniformGrid, @@ -66,15 +66,13 @@ fn interpolate_points_to_cell_data_generic( density_map.for_each(|flat_point_index, point_value| { let point = grid.try_unflatten_point_index(flat_point_index).unwrap(); - // We want to find edges that cross the iso-surface, - // therefore we can choose to either skip all points above or below the threshold. - // - // In most scenes, the sparse density map should contain more entries above than - // below the threshold, as it contains the whole fluid interior, whereas areas completely - // devoid of fluid are not part of the density map. - // - // Skip points with densities above the threshold to improve efficiency - if point_value > iso_surface_threshold { + // We want to find edges that cross the iso-surface, i.e. edges where one point is above + // and the other below the threshold. To not process an edge twice, we skip all points + // that are already below the iso-surface threshold. Note that we cannot do it the other + // way around, as some points below the threshold might not be part of the density map + // at all (e.g. points outside the kernel evaluation radius). This could lead to missing + // edges that go directly from above the threshold to e.g. zero. + if point_value < iso_surface_threshold { return; } @@ -89,15 +87,13 @@ fn interpolate_points_to_cell_data_generic( let neighbor_value = if let Some(v) = density_map.get(flat_neighbor_index) { v } else { - // Neighbors that are not in the point-value map were outside of the kernel evaluation radius. - // This should only happen for cells that are completely outside of the compact support of a particle. - // The point-value map has to be consistent such that for each cell, where at least one point-value - // is missing like this, the cell has to be completely below the iso-surface threshold. - continue; + // Neighbors that are not in the point-value map were outside the kernel evaluation radius. + // Assume zero density for these points. + R::zero() }; // Skip edges that don't cross the iso-surface - if !(neighbor_value > iso_surface_threshold) { + if !(neighbor_value < iso_surface_threshold) { continue; } @@ -132,8 +128,8 @@ fn interpolate_points_to_cell_data_generic( ); cell_data_entry.iso_surface_vertices[local_edge_index] = Some(vertex_index); - // Mark the neighbor as above the iso-surface threshold - let local_vertex_index = cell.local_point_index_of(neighbor.index()).unwrap(); + // Mark the corner of the current point as above the iso-surface threshold + let local_vertex_index = cell.local_point_index_of(point.index()).unwrap(); cell_data_entry.corner_above_threshold[local_vertex_index] = RelativeToThreshold::Above; } diff --git a/splashsurf_lib/tests/integration_tests/mod.rs b/splashsurf_lib/tests/integration_tests/mod.rs index 2a21be6..f52c33f 100644 --- a/splashsurf_lib/tests/integration_tests/mod.rs +++ b/splashsurf_lib/tests/integration_tests/mod.rs @@ -3,3 +3,4 @@ pub mod test_full; #[cfg(feature = "io")] pub mod test_mesh; pub mod test_neighborhood_search; +pub mod test_simple; diff --git a/splashsurf_lib/tests/integration_tests/test_simple.rs b/splashsurf_lib/tests/integration_tests/test_simple.rs new file mode 100644 index 0000000..7aa18b6 --- /dev/null +++ b/splashsurf_lib/tests/integration_tests/test_simple.rs @@ -0,0 +1,141 @@ +use nalgebra::Vector3; +use splashsurf_lib::marching_cubes::check_mesh_consistency; +use splashsurf_lib::{ + reconstruct_surface, Aabb3d, GridDecompositionParameters, Parameters, Real, + SpatialDecomposition, +}; + +enum Strategy { + Global, + SubdomainGrid, +} + +fn params_with_aabb( + particle_radius: R, + compact_support_radius: R, + cube_size: R, + iso_surface_threshold: R, + domain_aabb: Option>, + strategy: Strategy, +) -> Parameters { + let compact_support_radius = particle_radius * compact_support_radius; + let cube_size = particle_radius * cube_size; + + let mut parameters = Parameters { + particle_radius, + rest_density: R::from_f64(1000.0).unwrap(), + compact_support_radius, + cube_size, + iso_surface_threshold, + particle_aabb: domain_aabb, + enable_multi_threading: false, + spatial_decomposition: None, + global_neighborhood_list: false, + }; + + match strategy { + Strategy::Global => {} + Strategy::SubdomainGrid => { + parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( + GridDecompositionParameters { + subdomain_num_cubes_per_dim: 64, + }, + )) + } + } + + parameters +} + +fn params( + particle_radius: R, + compact_support_radius: R, + cube_size: R, + iso_surface_threshold: R, + strategy: Strategy, +) -> Parameters { + params_with_aabb( + particle_radius, + compact_support_radius, + cube_size, + iso_surface_threshold, + None, + strategy, + ) +} + +/// This tests ensures that a surface is reconstructed properly if a single edge is above the threshold +/// on one side but outside the compact support on the other side. +#[test] +fn test_edge_above_threshold_to_outside_of_compact_support_global() { + let particle_positions: Vec> = vec![Vector3::new(0.01, 0.0, 0.0)]; + let parameters = params(1.0, 1.0, 1.0, 0.1, Strategy::Global); + let reconstruction = + reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap(); + + println!( + "Reconstructed mesh from {} particles has {} triangles.", + particle_positions.len(), + reconstruction.mesh().triangles.len() + ); + + assert_eq!( + reconstruction.mesh().vertices.len(), + 6, + "Number of vertices" + ); + assert_eq!( + reconstruction.mesh().triangles.len(), + 8, + "Number of triangles" + ); + + // Ensure that the mesh does not have a boundary + if let Err(e) = check_mesh_consistency( + reconstruction.grid(), + reconstruction.mesh(), + true, + true, + true, + ) { + eprintln!("{}", e); + panic!("Mesh contains topological/manifold errors"); + } +} + +#[test] +fn test_edge_above_threshold_to_outside_of_compact_support_subdomains() { + let particle_positions: Vec> = vec![Vector3::new(0.01, 0.0, 0.0)]; + let parameters = params(1.0, 1.0, 1.0, 0.1, Strategy::SubdomainGrid); + let reconstruction = + reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap(); + + println!( + "Reconstructed mesh from {} particles has {} triangles.", + particle_positions.len(), + reconstruction.mesh().triangles.len() + ); + + assert_eq!( + reconstruction.mesh().vertices.len(), + 6, + "Number of vertices" + ); + assert_eq!( + reconstruction.mesh().triangles.len(), + 8, + "Number of triangles" + ); + + // Ensure that the mesh does not have a boundary + if let Err(e) = check_mesh_consistency( + reconstruction.grid(), + reconstruction.mesh(), + true, + true, + true, + ) { + eprintln!("{}", e); + panic!("Mesh contains topological/manifold errors"); + } +}