diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ff87ed..edf5d37 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,10 @@ The following changes are present in the `main` branch of the repository and are not yet part of a release: - - N/A + - Lib: Remove octree-based domain decomposition (superseded by regular grid subdivision approach `--subdomain-grid`) + - CLI: Set `--subdomain-grid=on` by default + - CLI: Remove all arguments for octree-based domain decomposition + - CLI: Remove options to output some debug files (octree grid, density map, etc.) ## Version 0.10.0 diff --git a/splashsurf/src/reconstruction.rs b/splashsurf/src/reconstruction.rs index 6b0ca4c..cc7129b 100644 --- a/splashsurf/src/reconstruction.rs +++ b/splashsurf/src/reconstruction.rs @@ -4,10 +4,10 @@ use clap::value_parser; use indicatif::{ProgressBar, ProgressStyle}; use log::info; use rayon::prelude::*; -use splashsurf_lib::mesh::{AttributeData, Mesh3d, MeshAttribute, MeshWithData, PointCloud3d}; +use splashsurf_lib::mesh::{AttributeData, Mesh3d, MeshAttribute, MeshWithData}; use splashsurf_lib::nalgebra::{Unit, Vector3}; use splashsurf_lib::sph_interpolation::SphInterpolator; -use splashsurf_lib::{density_map, profile, Aabb3d, Index, Real}; +use splashsurf_lib::{profile, Aabb3d, Index, Real}; use std::borrow::Cow; use std::convert::TryFrom; use std::path::PathBuf; @@ -122,7 +122,7 @@ pub struct ReconstructSubcommandArgs { #[arg( help_heading = ARGS_OCTREE, long, - default_value = "off", + default_value = "on", value_name = "off|on", ignore_case = true, require_equals = true @@ -132,55 +132,6 @@ pub struct ReconstructSubcommandArgs { #[arg(help_heading = ARGS_OCTREE, long, default_value="64")] pub subdomain_cubes: u32, - /// Enable spatial decomposition using an octree (faster) instead of a global approach - #[arg( - help_heading = ARGS_OCTREE, - long, - default_value = "on", - value_name = "off|on", - ignore_case = true, - require_equals = true - )] - pub octree_decomposition: Switch, - /// Enable stitching of the disconnected local meshes resulting from the reconstruction when spatial decomposition is enabled (slower, but without stitching meshes will not be closed) - #[arg( - help_heading = ARGS_OCTREE, - long, - default_value = "on", - value_name = "off|on", - ignore_case = true, - require_equals = true - )] - pub octree_stitch_subdomains: Switch, - /// The maximum number of particles for leaf nodes of the octree, default is to compute it based on the number of threads and particles - #[arg(help_heading = ARGS_OCTREE, long)] - pub octree_max_particles: Option, - /// Safety factor applied to the kernel compact support radius when it's used as a margin to collect ghost particles in the leaf nodes when performing the spatial decomposition - #[arg(help_heading = ARGS_OCTREE, long)] - pub octree_ghost_margin_factor: Option, - /// Enable computing particle densities in a global step before domain decomposition (slower) - #[arg( - help_heading = ARGS_OCTREE, - long, - default_value = "off", - value_name = "off|on", - ignore_case = true, - require_equals = true - )] - pub octree_global_density: Switch, - /// Enable computing particle densities per subdomain but synchronize densities for ghost-particles (faster, recommended). - /// Note: if both this and global particle density computation is disabled the ghost particle margin has to be increased to at least 2.0 - /// to compute correct density values for ghost particles. - #[arg( - help_heading = ARGS_OCTREE, - long, - default_value = "on", - value_name = "off|on", - ignore_case = true, - require_equals = true - )] - pub octree_sync_local_density: Switch, - /// Enable omputing surface normals at the mesh vertices and write them to the output file #[arg( help_heading = ARGS_INTERP, @@ -337,15 +288,6 @@ pub struct ReconstructSubcommandArgs { )] pub output_raw_mesh: Switch, - /// Optional filename for writing the point cloud representation of the intermediate density map to disk - #[arg(help_heading = ARGS_DEBUG, long, value_parser = value_parser!(PathBuf))] - pub output_dm_points: Option, - /// Optional filename for writing the grid representation of the intermediate density map to disk - #[arg(help_heading = ARGS_DEBUG, long, value_parser = value_parser!(PathBuf))] - pub output_dm_grid: Option, - /// Optional filename for writing the octree used to partition the particles to disk - #[arg(help_heading = ARGS_DEBUG, long, value_parser = value_parser!(PathBuf))] - pub output_octree: Option, /// Enable checking the final mesh for holes and non-manifold edges and vertices #[arg( help_heading = ARGS_DEBUG, @@ -472,7 +414,7 @@ mod arguments { use log::info; use regex::{escape, Regex}; use splashsurf_lib::nalgebra::Vector3; - use splashsurf_lib::{Aabb3d, ParticleDensityComputationStrategy}; + use splashsurf_lib::Aabb3d; use std::convert::TryFrom; use std::fs; use std::path::{Path, PathBuf}; @@ -577,37 +519,7 @@ mod arguments { Some(splashsurf_lib::SpatialDecomposition::UniformGrid( splashsurf_lib::GridDecompositionParameters { subdomain_num_cubes_per_dim: args.subdomain_cubes, - }, - )) - } else if args.octree_decomposition.into_bool() { - let subdivision_criterion = if let Some(max_particles) = args.octree_max_particles { - splashsurf_lib::SubdivisionCriterion::MaxParticleCount(max_particles) - } else { - splashsurf_lib::SubdivisionCriterion::MaxParticleCountAuto - }; - let ghost_particle_safety_factor = args.octree_ghost_margin_factor; - let enable_stitching = args.octree_stitch_subdomains.into_bool(); - - let particle_density_computation = if args.octree_global_density.into_bool() - && args.octree_sync_local_density.into_bool() - { - return Err(anyhow!("Cannot enable both global and merged local particle density computation at the same time. Switch off at least one.")); - } else { - if args.octree_global_density.into_bool() { - ParticleDensityComputationStrategy::Global - } else if args.octree_sync_local_density.into_bool() { - ParticleDensityComputationStrategy::SynchronizeSubdomains - } else { - ParticleDensityComputationStrategy::IndependentSubdomains - } - }; - - Some(splashsurf_lib::SpatialDecomposition::Octree( - splashsurf_lib::OctreeDecompositionParameters { - subdivision_criterion, - ghost_particle_safety_factor, - enable_stitching, - particle_density_computation, + ..Default::default() }, )) } else { @@ -673,9 +585,6 @@ mod arguments { is_sequence: bool, input_file: PathBuf, output_file: PathBuf, - output_density_map_points_file: Option, - output_density_map_grid_file: Option, - output_octree_file: Option, sequence_range: (Option, Option), } @@ -685,17 +594,11 @@ mod arguments { input_file: P, output_base_path: Option

, output_file: P, - output_density_map_points_file: Option

, - output_density_map_grid_file: Option

, - output_octree_file: Option

, sequence_range: (Option, Option), ) -> Result { let input_file = input_file.into(); let output_base_path = output_base_path.map(|p| p.into()); let output_file = output_file.into(); - let output_density_map_points_file = output_density_map_points_file.map(|p| p.into()); - let output_density_map_grid_file = output_density_map_grid_file.map(|p| p.into()); - let output_octree_file = output_octree_file.map(|p| p.into()); if let (Some(start), Some(end)) = sequence_range { if start > end { @@ -727,11 +630,6 @@ mod arguments { is_sequence, input_file, output_file, - output_density_map_points_file: output_density_map_points_file - .map(|f| output_base_path.join(f)), - output_density_map_grid_file: output_density_map_grid_file - .map(|f| output_base_path.join(f)), - output_octree_file: output_octree_file.map(|f| output_base_path.join(f)), sequence_range, }) } else { @@ -739,9 +637,6 @@ mod arguments { is_sequence, input_file, output_file, - output_density_map_points_file, - output_density_map_grid_file, - output_octree_file, sequence_range, }) } @@ -821,14 +716,7 @@ mod arguments { let output_filename_i = output_pattern.replace("{}", index); let output_file_i = output_dir.join(output_filename_i); - paths.push(ReconstructionRunnerPaths::new( - input_file_i, - output_file_i, - // Don't write density maps etc. when processing a sequence of files - None, - None, - None, - )); + paths.push(ReconstructionRunnerPaths::new(input_file_i, output_file_i)); } } @@ -851,9 +739,6 @@ mod arguments { ReconstructionRunnerPaths::new( self.input_file.clone(), self.output_file.clone(), - self.output_density_map_points_file.clone(), - self.output_density_map_grid_file.clone(), - self.output_octree_file.clone(), ); 1 ] @@ -950,9 +835,6 @@ mod arguments { args.input_file_or_sequence.clone(), args.output_dir.clone(), output_filename, - args.output_dm_points.clone(), - args.output_dm_grid.clone(), - args.output_octree.clone(), (args.start_index, args.end_index), ) } @@ -963,25 +845,13 @@ mod arguments { pub(crate) struct ReconstructionRunnerPaths { pub input_file: PathBuf, pub output_file: PathBuf, - pub output_density_map_points_file: Option, - pub output_density_map_grid_file: Option, - pub output_octree_file: Option, } impl ReconstructionRunnerPaths { - fn new( - input_file: PathBuf, - output_file: PathBuf, - output_density_map_points_file: Option, - output_density_map_grid_file: Option, - output_octree_file: Option, - ) -> Self { + fn new(input_file: PathBuf, output_file: PathBuf) -> Self { ReconstructionRunnerPaths { input_file, output_file, - output_density_map_points_file, - output_density_map_grid_file, - output_octree_file, } } } @@ -1423,74 +1293,6 @@ pub(crate) fn reconstruction_pipeline_generic( info!("Done."); } - // Store octree leaf nodes as hex cells - if let Some(output_octree_file) = &paths.output_octree_file { - info!("Writing octree to \"{}\"...", output_octree_file.display()); - io::vtk_format::write_vtk( - reconstruction.octree().unwrap().hexmesh(grid, true), - output_octree_file, - "mesh", - ) - .with_context(|| { - format!( - "Failed to write octree to output file \"{}\"", - output_octree_file.display() - ) - })?; - info!("Done."); - } - - // Store point cloud density map - if let Some(output_density_map_points_file) = &paths.output_density_map_points_file { - info!("Constructing density map point cloud..."); - let density_map = reconstruction - .density_map() - .ok_or_else(|| anyhow::anyhow!("No density map was created during reconstruction"))?; - - let point_cloud: PointCloud3d = { - let mut points = Vec::with_capacity(density_map.len()); - density_map.for_each(|flat_point_index, _| { - let point = grid.try_unflatten_point_index(flat_point_index).unwrap(); - points.push(grid.point_coordinates(&point)); - }); - - PointCloud3d::new(points) - }; - - info!( - "Saving density map point cloud to \"{}\"...", - output_density_map_points_file.display() - ); - - io::vtk_format::write_vtk( - &point_cloud, - output_density_map_points_file, - "density_map_points", - )?; - - info!("Done."); - } - - // Store hex-mesh density map - if let Some(output_density_map_grid_file) = &paths.output_density_map_grid_file { - info!("Constructing density map hex mesh..."); - let density_map = reconstruction - .density_map() - .ok_or_else(|| anyhow::anyhow!("No density map was created during reconstruction"))?; - - let density_mesh = - density_map::sparse_density_map_to_hex_mesh(&density_map, &grid, R::zero()); - - info!( - "Saving density map hex mesh to \"{}\"...", - output_density_map_grid_file.display() - ); - - io::vtk_format::write_vtk(density_mesh, output_density_map_grid_file, "density_map")?; - - info!("Done."); - } - if postprocessing.check_mesh_closed || postprocessing.check_mesh_manifold || postprocessing.check_mesh_debug diff --git a/splashsurf_lib/README.md b/splashsurf_lib/README.md index 13aaf09..5aae6a5 100644 --- a/splashsurf_lib/README.md +++ b/splashsurf_lib/README.md @@ -37,13 +37,15 @@ For each of the features, `splashsurf_lib` re-exports the corresponding dependen Currently, only one method based on a "spatial hashing" strategy is implemented. +`TODO: This section is missing a description of the domain decomposition for more efficient parallelization` + **Short summary**: The fluid density is evaluated or mapped onto a sparse grid using spatial hashing in the support radius of each fluid particle. This implies that memory is only allocated in areas where the fluid density is non-zero. This is in contrast to a naive approach where the marching cubes background grid is allocated for the whole domain. Finally, the marching cubes reconstruction is performed only in those grid cells where an edge crosses the surface threshold. Cells completely in the interior of the fluid are skipped in the marching cubes phase. **Individual steps**: 1. Construct a "virtual background grid" with the desired resolution of the marching cubes algorithm. In the end, the procedure will place a single surface mesh vertex on each edge of this virtual grid, where the fluid surface crosses the edge (or rather, where the fluid density crosses the specified threshold). Virtual means that no storage is actually allocated for this grid yet; only its topology is used implicitly later. 2. Compute the density of each fluid particle - - Perform a neighborhood search - - Per particle, evaluate an SPH sum over the neighbors to compute its density (based on input parameters of kernel radius and particle rest mass) + - Perform a neighborhood search + - Per particle, evaluate an SPH sum over the neighbors to compute its density (based on input parameters of kernel radius and particle rest mass) 3. Optional: filter out (or rather mask as inactive) single particles if the user provided a "splash detection radius". This is done by performing an additional neighborhood search using this splash detection radius instead of the kernel radius. 4. Compute a "sparse density map": a map from the index of a vertex of the virtual background grid to the corresponding fluid density value. The map will only contain entries for vertices where the fluid density is non-zero. Construction of the map: - Iterate over all active particles diff --git a/splashsurf_lib/benches/benches/bench_full.rs b/splashsurf_lib/benches/benches/bench_full.rs index 8b8cf2d..def7307 100644 --- a/splashsurf_lib/benches/benches/bench_full.rs +++ b/splashsurf_lib/benches/benches/bench_full.rs @@ -4,9 +4,8 @@ use splashsurf_lib::io::particles_from_file; #[allow(dead_code)] use splashsurf_lib::io::vtk_format::write_vtk; use splashsurf_lib::{ - reconstruct_surface, reconstruct_surface_inplace, GridDecompositionParameters, - OctreeDecompositionParameters, Parameters, ParticleDensityComputationStrategy, - SpatialDecomposition, SubdivisionCriterion, SurfaceReconstruction, + reconstruct_surface, reconstruct_surface_inplace, GridDecompositionParameters, Parameters, + SpatialDecomposition, SurfaceReconstruction, }; use std::time::Duration; @@ -121,46 +120,6 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) { }) }); - group.bench_function("surface_reconstruction_dam_break_par_octree", |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: false, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruction = - reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap() - }) - }); - - group.bench_function( - "surface_reconstruction_dam_break_par_octree_stitching", - |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: true, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruction = - reconstruct_surface::(particle_positions.as_slice(), ¶meters) - .unwrap() - }) - }, - ); - group.bench_function("surface_reconstruction_dam_break_par_grid_64", |b| { b.iter(|| { let mut parameters = parameters.clone(); @@ -220,46 +179,6 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) { }) }); - group.bench_function("surface_reconstruction_double_dam_break_par_octree", |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: false, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruction = - reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap() - }) - }); - - group.bench_function( - "surface_reconstruction_double_dam_break_par_octree_stitching", - |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: true, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruction = - reconstruct_surface::(particle_positions.as_slice(), ¶meters) - .unwrap() - }) - }, - ); - group.bench_function("surface_reconstruction_double_dam_break_par_grid_64", |b| { b.iter(|| { let mut parameters = parameters.clone(); @@ -326,56 +245,6 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) { }, ); - group.bench_function( - "surface_reconstruction_double_dam_break_inplace_par_octree", - |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: false, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruct_surface_inplace::( - particle_positions.as_slice(), - ¶meters, - &mut reconstruction, - ) - .unwrap() - }) - }, - ); - - group.bench_function( - "surface_reconstruction_double_dam_break_inplace_par_octree_stitching", - |b| { - b.iter(|| { - let mut parameters = parameters.clone(); - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(1.0), - enable_stitching: true, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - - reconstruct_surface_inplace::( - particle_positions.as_slice(), - ¶meters, - &mut reconstruction, - ) - .unwrap() - }) - }, - ); - group.bench_function( "surface_reconstruction_double_dam_break_inplace_par_grid_64", |b| { diff --git a/splashsurf_lib/benches/benches/bench_mesh.rs b/splashsurf_lib/benches/benches/bench_mesh.rs index 6df2448..7ee2464 100644 --- a/splashsurf_lib/benches/benches/bench_mesh.rs +++ b/splashsurf_lib/benches/benches/bench_mesh.rs @@ -2,8 +2,7 @@ use criterion::{criterion_group, Criterion}; use splashsurf_lib::io::particles_from_file; use splashsurf_lib::nalgebra::Vector3; use splashsurf_lib::{ - reconstruct_surface, OctreeDecompositionParameters, Parameters, - ParticleDensityComputationStrategy, SpatialDecomposition, SubdivisionCriterion, + reconstruct_surface, GridDecompositionParameters, Parameters, SpatialDecomposition, SurfaceReconstruction, }; use std::path::Path; @@ -24,14 +23,8 @@ fn reconstruct_particles>(particle_file: P) -> SurfaceReconstruct iso_surface_threshold: 0.6, particle_aabb: None, enable_multi_threading: true, - spatial_decomposition: Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: None, - enable_stitching: true, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, + spatial_decomposition: Some(SpatialDecomposition::UniformGrid( + GridDecompositionParameters::default(), )), global_neighborhood_list: false, }; diff --git a/splashsurf_lib/benches/benches/bench_octree.rs b/splashsurf_lib/benches/benches/bench_octree.rs deleted file mode 100644 index ff573fd..0000000 --- a/splashsurf_lib/benches/benches/bench_octree.rs +++ /dev/null @@ -1,53 +0,0 @@ -use criterion::{criterion_group, BatchSize, Criterion}; -use nalgebra::Vector3; -use splashsurf_lib::io::particles_from_file; -use splashsurf_lib::octree::Octree; -use splashsurf_lib::{grid_for_reconstruction, SubdivisionCriterion, UniformGrid}; -use std::time::Duration; - -pub fn subdivide_recursively_benchmark(c: &mut Criterion) { - let particle_positions: &Vec> = - &particles_from_file("../data/hilbert_46843_particles.bgeo").unwrap(); - let particles_per_cell = 20000; - - let particle_radius = 0.011; - let compact_support_radius = 4.0 * particle_radius; - - let cube_size = 1.5 * particle_radius; - let grid: &UniformGrid = &grid_for_reconstruction( - particle_positions.as_slice(), - particle_radius, - compact_support_radius, - cube_size, - None, - true, - ) - .unwrap(); - - let mut group = c.benchmark_group("octree subdivision"); - group.sample_size(80); - group.warm_up_time(Duration::from_secs(5)); - group.measurement_time(Duration::from_secs(25)); - - let get_tree = || Octree::new(&grid, particle_positions.len()); - - group.bench_function("subdivide_recursively_margin_par", move |b| { - b.iter_batched( - get_tree, - |mut tree| { - tree.par_subdivide_recursively_margin( - grid, - particle_positions.as_slice(), - SubdivisionCriterion::MaxParticleCount(particles_per_cell), - compact_support_radius, - false, - ) - }, - BatchSize::LargeInput, - ) - }); - - group.finish(); -} - -criterion_group!(bench_octree, subdivide_recursively_benchmark); diff --git a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs index c853d69..7b0180e 100644 --- a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs +++ b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs @@ -23,6 +23,7 @@ fn parameters_canyon() -> Parameters { spatial_decomposition: Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { subdomain_num_cubes_per_dim: 32, + ..Default::default() }, )), global_neighborhood_list: false, @@ -50,6 +51,7 @@ pub fn grid_canyon(c: &mut Criterion) { parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { subdomain_num_cubes_per_dim: 32, + ..Default::default() }, )); reconstruction = @@ -64,6 +66,7 @@ pub fn grid_canyon(c: &mut Criterion) { parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { subdomain_num_cubes_per_dim: 48, + ..Default::default() }, )); reconstruction = @@ -78,6 +81,7 @@ pub fn grid_canyon(c: &mut Criterion) { parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { subdomain_num_cubes_per_dim: 64, + ..Default::default() }, )); reconstruction = @@ -111,6 +115,7 @@ pub fn grid_optimal_num_cubes_canyon(c: &mut Criterion) { parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { subdomain_num_cubes_per_dim: num_cubes, + ..Default::default() }, )); reconstruction = diff --git a/splashsurf_lib/benches/benches/mod.rs b/splashsurf_lib/benches/benches/mod.rs index 52c1c76..cd2402c 100644 --- a/splashsurf_lib/benches/benches/mod.rs +++ b/splashsurf_lib/benches/benches/mod.rs @@ -2,5 +2,4 @@ pub mod bench_aabb; pub mod bench_full; pub mod bench_mesh; pub mod bench_neighborhood; -pub mod bench_octree; pub mod bench_subdomain_grid; diff --git a/splashsurf_lib/benches/splashsurf_lib_benches.rs b/splashsurf_lib/benches/splashsurf_lib_benches.rs index c397e1f..e9be90b 100644 --- a/splashsurf_lib/benches/splashsurf_lib_benches.rs +++ b/splashsurf_lib/benches/splashsurf_lib_benches.rs @@ -6,13 +6,11 @@ use benches::bench_aabb::bench_aabb; use benches::bench_full::bench_full; use benches::bench_mesh::bench_mesh; use benches::bench_neighborhood::bench_neighborhood; -use benches::bench_octree::bench_octree; use benches::bench_subdomain_grid::bench_subdomain_grid; criterion_main!( bench_aabb, bench_mesh, - bench_octree, bench_full, bench_neighborhood, bench_subdomain_grid, diff --git a/splashsurf_lib/src/density_map.rs b/splashsurf_lib/src/density_map.rs index 9aebac8..0b8802b 100644 --- a/splashsurf_lib/src/density_map.rs +++ b/splashsurf_lib/src/density_map.rs @@ -22,7 +22,7 @@ use crate::aabb::Aabb3d; use crate::kernel::DiscreteSquaredDistanceCubicKernel; use crate::mesh::{HexMesh3d, MeshAttribute, MeshWithData}; use crate::neighborhood_search::NeighborhoodList; -use crate::uniform_grid::{OwningSubdomainGrid, Subdomain, UniformGrid}; +use crate::uniform_grid::UniformGrid; use crate::utils::{ChunkSize, ParallelPolicy}; use crate::{new_map, profile, HashState, Index, MapType, ParallelMapType, Real}; use dashmap::ReadOnlyView as ReadDashMap; @@ -224,6 +224,12 @@ pub enum DensityMap { DashMap(ReadDashMap), } +impl Default for DensityMap { + fn default() -> Self { + DensityMap::Standard(MapType::default()) + } +} + impl From> for DensityMap { fn from(map: MapType) -> Self { Self::Standard(map) @@ -261,17 +267,6 @@ impl DensityMap { } } - /// Returns a mutable reference to the contained standard map, replaces itself if not of standard type - fn standard_or_insert_mut(&mut self) -> &mut MapType { - match self { - DensityMap::Standard(map) => return map, - _ => {} - } - - *self = new_map().into(); - self.standard_or_insert_mut() - } - /// Calls a closure for each `(flat_point_index, density_value)` tuple in the map pub fn for_each(&self, f: F) { let mut f = f; @@ -286,7 +281,6 @@ impl DensityMap { #[inline(never)] pub fn generate_sparse_density_map( grid: &UniformGrid, - subdomain: Option<&OwningSubdomainGrid>, particle_positions: &[Vector3], particle_densities: &[R], active_particles: Option<&[usize]>, @@ -305,44 +299,27 @@ pub fn generate_sparse_density_map( } ); - if let Some(subdomain) = subdomain { - if allow_threading { - panic!("Multi threading not implemented for density map with subdomain"); - } else { - sequential_generate_sparse_density_map_subdomain( - subdomain, - particle_positions, - particle_densities, - active_particles, - particle_rest_mass, - compact_support_radius, - cube_size, - density_map, - )?; - } + if allow_threading { + *density_map = parallel_generate_sparse_density_map( + grid, + particle_positions, + particle_densities, + active_particles, + particle_rest_mass, + compact_support_radius, + cube_size, + )? } else { - if allow_threading { - *density_map = parallel_generate_sparse_density_map( - grid, - particle_positions, - particle_densities, - active_particles, - particle_rest_mass, - compact_support_radius, - cube_size, - )? - } else { - *density_map = sequential_generate_sparse_density_map( - grid, - particle_positions, - particle_densities, - active_particles, - particle_rest_mass, - compact_support_radius, - cube_size, - )? - } - }; + *density_map = sequential_generate_sparse_density_map( + grid, + particle_positions, + particle_densities, + active_particles, + particle_rest_mass, + compact_support_radius, + cube_size, + )? + } trace!( "Sparse density map was constructed. (Output: density map with {} grid point data entries)", @@ -399,55 +376,6 @@ pub fn sequential_generate_sparse_density_map( Ok(sparse_densities.into()) } -/// Computes a sparse density map for the fluid restricted to the specified subdomain -#[inline(never)] -pub fn sequential_generate_sparse_density_map_subdomain( - subdomain: &OwningSubdomainGrid, - particle_positions: &[Vector3], - particle_densities: &[R], - active_particles: Option<&[usize]>, - particle_rest_mass: R, - compact_support_radius: R, - cube_size: R, - density_map: &mut DensityMap, -) -> Result<(), DensityMapError> { - profile!("sequential_generate_sparse_density_map_subdomain"); - - let mut sparse_densities = density_map.standard_or_insert_mut(); - sparse_densities.clear(); - - let density_map_generator = SparseDensityMapGenerator::try_new( - &subdomain.global_grid(), - compact_support_radius, - cube_size, - particle_rest_mass, - )?; - - let process_particle = |particle_data: (&Vector3, R)| { - let (particle, particle_density) = particle_data; - density_map_generator.compute_particle_density_contribution_subdomain( - subdomain, - &mut sparse_densities, - particle, - particle_density, - ); - }; - - match active_particles { - None => particle_positions - .iter() - .zip(particle_densities.iter().copied()) - .for_each(process_particle), - Some(indices) => indices - .iter() - .map(|&i| &particle_positions[i]) - .zip(indices.iter().map(|&i| particle_densities[i])) - .for_each(process_particle), - } - - Ok(()) -} - /// Computes a sparse density map for the fluid based on the specified background grid, multi-threaded implementation #[inline(never)] pub fn parallel_generate_sparse_density_map( @@ -702,84 +630,6 @@ impl SparseDensityMapGenerator { ); } - /// Computes all density contributions of a particle to a subdomain of the background grid into the given map - fn compute_particle_density_contribution_subdomain( - &self, - subdomain: &OwningSubdomainGrid, - sparse_densities: &mut MapType, - particle: &Vector3, - particle_density: R, - ) { - let grid = subdomain.global_grid(); - let subdomain_grid = subdomain.subdomain_grid(); - let subdomain_offset = subdomain.subdomain_offset(); - - // Skip particles outside of allowed domain - if !self.allowed_domain.contains_point(particle) { - return; - } - - let global_subdomain_min_point = subdomain_offset; - // Note: max supported point is one past the actual max point index - let global_subdomain_max_point = [ - global_subdomain_min_point[0] + subdomain_grid.points_per_dim()[0], - global_subdomain_min_point[1] + subdomain_grid.points_per_dim()[1], - global_subdomain_min_point[2] + subdomain_grid.points_per_dim()[2], - ]; - - // Compute cuboid region of grid points that may be affected by the particle - // This excludes grid points outside of the current subdomain - - let min_supported_point_ijk = { - let cell_ijk = grid.enclosing_cell(particle); - [ - (cell_ijk[0] - self.half_supported_cells).max(global_subdomain_min_point[0]), - (cell_ijk[1] - self.half_supported_cells).max(global_subdomain_min_point[1]), - (cell_ijk[2] - self.half_supported_cells).max(global_subdomain_min_point[2]), - ] - }; - - let max_supported_point_ijk = { - [ - (min_supported_point_ijk[0] + self.supported_points) - .min(global_subdomain_max_point[0]), - (min_supported_point_ijk[1] + self.supported_points) - .min(global_subdomain_max_point[1]), - (min_supported_point_ijk[2] + self.supported_points) - .min(global_subdomain_max_point[2]), - ] - }; - - // Check if lower corner of the supported domain is above of subdomain - if min_supported_point_ijk - .iter() - .copied() - .zip(global_subdomain_max_point.iter().copied()) - .any(|(min_support, subdomain_max)| min_support > subdomain_max) - { - return; - } - - // Check if upper corner of the supported domain is below of subdomain - if max_supported_point_ijk - .iter() - .copied() - .zip(global_subdomain_min_point.iter().copied()) - .any(|(max_support, subdomain_min)| max_support < subdomain_min) - { - return; - } - - self.particle_support_loop( - sparse_densities, - grid, - &min_supported_point_ijk, - &max_supported_point_ijk, - particle, - particle_density, - ); - } - /// Loops over a cube of background grid points that are potentially in the support radius of the particle and evaluates density contributions #[inline(always)] fn particle_support_loop( diff --git a/splashsurf_lib/src/lib.rs b/splashsurf_lib/src/lib.rs index 20ac58e..ddd8b18 100644 --- a/splashsurf_lib/src/lib.rs +++ b/splashsurf_lib/src/lib.rs @@ -36,14 +36,12 @@ pub use vtkio; pub use crate::aabb::{Aabb2d, Aabb3d, AxisAlignedBoundingBox}; pub use crate::density_map::DensityMap; -pub use crate::octree::SubdivisionCriterion; pub use crate::traits::{Index, Real, RealConvert, ThreadSafe}; pub use crate::uniform_grid::UniformGrid; use crate::density_map::DensityMapError; use crate::marching_cubes::MarchingCubesError; use crate::mesh::TriMesh3d; -use crate::octree::Octree; use crate::uniform_grid::GridConstructionError; use crate::workspace::ReconstructionWorkspace; @@ -65,10 +63,8 @@ pub mod kernel; pub mod marching_cubes; pub mod mesh; pub mod neighborhood_search; -pub mod octree; pub mod postprocessing; pub(crate) mod reconstruction; -mod reconstruction_octree; pub mod sph_interpolation; pub mod topology; mod traits; @@ -121,33 +117,21 @@ fn new_parallel_map() -> ParallelMapType { /// Approach used for spatial decomposition of the surface reconstruction and its parameters #[derive(Clone, Debug)] -pub enum SpatialDecomposition { +pub enum SpatialDecomposition { /// Use a uniform grid of subdomains with contiguous (dense) marching cubes grids per subdomain /// /// Only subdomains containing at least one particle will be processed. /// The small contiguous grid per subdomain make this approach very cache efficient. UniformGrid(GridDecompositionParameters), - /// Use an octree for decomposition with hashing based (sparse) marching cubes grids per subdomain - Octree(OctreeDecompositionParameters), } /// Default parameters for the spatial decomposition use the uniform grid based decomposition approach -impl Default for SpatialDecomposition { +impl Default for SpatialDecomposition { fn default() -> Self { Self::UniformGrid(GridDecompositionParameters::default()) } } -impl SpatialDecomposition { - /// Tries to convert the parameters from one [`Real`] type to another [`Real`] type, returns `None` if conversion fails - pub fn try_convert(&self) -> Option> { - match &self { - Self::UniformGrid(g) => Some(SpatialDecomposition::UniformGrid(g.clone())), - Self::Octree(o) => o.try_convert().map(|o| SpatialDecomposition::Octree(o)), - } - } -} - /// Parameters for the uniform grid-based spatial decomposition #[derive(Clone, Debug)] pub struct GridDecompositionParameters { @@ -163,87 +147,6 @@ impl Default for GridDecompositionParameters { } } -/// Parameters for the octree-based spatial decomposition -#[derive(Clone, Debug)] -pub struct OctreeDecompositionParameters { - /// Criterion used for subdivision of the octree cells - pub subdivision_criterion: SubdivisionCriterion, - /// Safety factor applied to the kernel radius when it's used as a margin to collect ghost particles in the leaf nodes - pub ghost_particle_safety_factor: Option, - /// Whether to enable stitching of all disjoint subdomain meshes to a global manifold mesh - pub enable_stitching: bool, - /// Which method to use for computing the densities of the particles - pub particle_density_computation: ParticleDensityComputationStrategy, -} - -/// Available strategies for the computation of the particle densities -#[derive(Copy, Clone, Debug)] -pub enum ParticleDensityComputationStrategy { - /// Compute the particle densities globally before performing domain decomposition. - /// - /// With this approach the particle densities are computed globally on all particles before any - /// domain decomposition is performed. - /// - /// This approach is guaranteed to lead to consistent results and does not depend on the following - /// decomposition. However, it is also by far the *slowest method* as global operations (especially - /// the global neighborhood search) are much slower. - Global, - /// Compute particle densities for all particles locally followed by a synchronization step. - /// - /// **This is the recommended approach.** - /// The particle densities will be evaluated for all particles per subdomain, possibly in parallel. - /// Afterwards, the values for all non-ghost particles are written to a global array. - /// This happens in a separate step before performing any reconstructions - /// For the following reconstruction procedure, each subdomain will update the densities of its ghost particles - /// from this global array. This ensures that all ghost-particles receive correct density values - /// without requiring to double the width of the ghost-particle margin just to ensure correct values - /// for the actual inner ghost-particles (i.e. in contrast to the completely local approach). - /// - /// The actual synchronization overhead is relatively low and this approach is often the fastest method. - /// - /// This approach should always lead consistent results. Only in very rare cases when a particle is not - /// uniquely assigned during domain decomposition this might lead to problems. If you encounter such - /// problems with this approach please report it as a bug. - SynchronizeSubdomains, - /// Compute densities locally per subdomain without global synchronization. - /// - /// The particle densities will be evaluated per subdomain on-the-fly just before the reconstruction - /// of the subdomain happens. In order to compute correct densities for the ghost particles of each - /// subdomain it is required that the ghost-particle margin is at least two times the kernel compact - /// support radius. This may add a lot of additional ghost-particles to each subdomain. - /// - /// If the ghost-particle margin is not set wide enough, this may lead to density differences on subdomain - /// boundaries. Otherwise this approach robust with respect to the classification of particles into the - /// subdomains. - IndependentSubdomains, -} - -impl Default for OctreeDecompositionParameters { - fn default() -> Self { - Self { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: None, - enable_stitching: true, - particle_density_computation: ParticleDensityComputationStrategy::SynchronizeSubdomains, - } - } -} - -impl OctreeDecompositionParameters { - /// Tries to convert the parameters from one [`Real`] type to another [`Real`] type, returns `None` if conversion fails - pub fn try_convert(&self) -> Option> { - Some(OctreeDecompositionParameters { - subdivision_criterion: self.subdivision_criterion.clone(), - ghost_particle_safety_factor: map_option!( - &self.ghost_particle_safety_factor, - r => r.try_convert()? - ), - enable_stitching: self.enable_stitching, - particle_density_computation: self.particle_density_computation, - }) - } -} - /// Parameters for the surface reconstruction #[derive(Clone, Debug)] pub struct Parameters { @@ -268,7 +171,7 @@ pub struct Parameters { pub enable_multi_threading: bool, /// Parameters for the spatial decomposition of the surface reconstruction /// If not provided, no spatial decomposition is performed and a global approach is used instead. - pub spatial_decomposition: Option>, + pub spatial_decomposition: Option, /// Whether to return the global particle neighborhood list from the reconstruction. /// Depending on the settings of the reconstruction, neighborhood lists are only computed locally /// in subdomains. Enabling this flag joins this data over all particles which can add a small overhead. @@ -286,7 +189,7 @@ impl Parameters { iso_surface_threshold: self.iso_surface_threshold.try_convert()?, particle_aabb: map_option!(&self.particle_aabb, aabb => aabb.try_convert()?), enable_multi_threading: self.enable_multi_threading, - spatial_decomposition: map_option!(&self.spatial_decomposition, sd => sd.try_convert()?), + spatial_decomposition: self.spatial_decomposition.clone(), global_neighborhood_list: self.global_neighborhood_list, }) } @@ -297,10 +200,6 @@ impl Parameters { pub struct SurfaceReconstruction { /// Background grid that was used as a basis for generating the density map for marching cubes grid: UniformGrid, - /// Octree constructed for domain decomposition (if octree-based spatial decomposition was enabled) - octree: Option>, - /// Point-based density map generated from the particles that was used as input to marching cubes - density_map: Option>, /// Per particle densities (contains only data of particles inside the domain) particle_densities: Option>, /// If an AABB was specified to restrict the reconstruction, this stores per input particle whether they were inside @@ -310,7 +209,7 @@ pub struct SurfaceReconstruction { /// Surface mesh that is the result of the surface reconstruction mesh: TriMesh3d, /// Workspace with allocated memory for subsequent surface reconstructions - workspace: ReconstructionWorkspace, + workspace: ReconstructionWorkspace, } impl Default for SurfaceReconstruction { @@ -318,8 +217,6 @@ impl Default for SurfaceReconstruction { fn default() -> Self { Self { grid: UniformGrid::new_zero(), - octree: None, - density_map: None, particle_densities: None, particle_neighbors: None, particle_inside_aabb: None, @@ -335,16 +232,6 @@ impl SurfaceReconstruction { &self.mesh } - /// Returns a reference to the octree generated for spatial decomposition of the input particles (mostly useful for debugging visualization) - pub fn octree(&self) -> Option<&Octree> { - self.octree.as_ref() - } - - /// Returns a reference to the sparse density map (discretized on the vertices of the background grid) that is used as input for marching cubes (always `None` when using domain decomposition) - pub fn density_map(&self) -> Option<&DensityMap> { - self.density_map.as_ref() - } - /// 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) pub fn particle_densities(&self) -> Option<&Vec> { self.particle_densities.as_ref() @@ -489,14 +376,7 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>( output_surface, )? } - Some(SpatialDecomposition::Octree(_)) => { - reconstruction_octree::reconstruct_surface_domain_decomposition( - particle_positions, - parameters, - output_surface, - )? - } - None => reconstruction_octree::reconstruct_surface_global( + None => reconstruction::reconstruct_surface_global( particle_positions, parameters, output_surface, @@ -512,12 +392,6 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>( Ok(()) } -/* -fn filter_particles(particle_positions: &[Vector3], particle_aabb: &Aabb3d, enable_multi_threading: bool) -> Vec> { - use rayon::prelude::*; - let is_inside = particle_positions.par_iter().map(|p| particle_aabb.contains_point(p)).collect::>(); -}*/ - /// Constructs the background grid for marching cubes based on the parameters supplied to the surface reconstruction pub fn grid_for_reconstruction( particle_positions: &[Vector3], diff --git a/splashsurf_lib/src/marching_cubes.rs b/splashsurf_lib/src/marching_cubes.rs index 81d0007..8af6794 100644 --- a/splashsurf_lib/src/marching_cubes.rs +++ b/splashsurf_lib/src/marching_cubes.rs @@ -1,23 +1,16 @@ //! Triangulation of [`DensityMap`]s using marching cubes -use crate::marching_cubes::narrow_band_extraction::{ - construct_mc_input, construct_mc_input_with_stitching_data, -}; -use crate::marching_cubes::triangulation::{ - triangulate, triangulate_with_criterion, DebugTriangleGenerator, TriangulationSkipBoundaryCells, -}; +use crate::marching_cubes::narrow_band_extraction::construct_mc_input; +use crate::marching_cubes::triangulation::triangulate; use crate::mesh::TriMesh3d; -use crate::uniform_grid::{DummySubdomain, OwningSubdomainGrid, Subdomain}; use crate::{new_map, profile, DensityMap, Index, MapType, Real, UniformGrid}; use nalgebra::Vector3; use thiserror::Error as ThisError; pub mod marching_cubes_lut; mod narrow_band_extraction; -mod stitching; mod triangulation; -pub(crate) use stitching::{stitch_surface_patches, SurfacePatch}; pub use triangulation::TriangulationError; /// Error enum for the marching cubes functions @@ -59,15 +52,6 @@ impl RelativeToThreshold { } } } - - /// Returns if the value is above the iso-surface threshold or `None` if the value is indeterminate - fn is_indeterminate(&self) -> bool { - if let RelativeToThreshold::Indeterminate = self { - true - } else { - false - } - } } /// Data for a single cell required by marching cubes @@ -121,78 +105,28 @@ pub fn triangulate_density_map( profile!("triangulate_density_map"); let mut mesh = TriMesh3d::default(); - triangulate_density_map_append(grid, None, density_map, iso_surface_threshold, &mut mesh)?; + triangulate_density_map_append(grid, density_map, iso_surface_threshold, &mut mesh)?; Ok(mesh) } /// Performs a marching cubes triangulation of a density map on the given background grid, appends triangles to the given mesh pub fn triangulate_density_map_append( grid: &UniformGrid, - subdomain: Option<&OwningSubdomainGrid>, density_map: &DensityMap, iso_surface_threshold: R, mesh: &mut TriMesh3d, ) -> Result<(), MarchingCubesError> { profile!("triangulate_density_map_append"); - let marching_cubes_data = if let Some(subdomain) = subdomain { - construct_mc_input( - subdomain, - &density_map, - iso_surface_threshold, - &mut mesh.vertices, - ) - } else { - let subdomain = DummySubdomain::new(grid); - construct_mc_input( - &subdomain, - &density_map, - iso_surface_threshold, - &mut mesh.vertices, - ) - }; - - triangulate(marching_cubes_data, mesh)?; - Ok(()) -} - -/// Performs triangulation of the given density map to a surface patch -pub(crate) fn triangulate_density_map_to_surface_patch( - subdomain: &OwningSubdomainGrid, - density_map: &DensityMap, - iso_surface_threshold: R, -) -> Result, MarchingCubesError> { - profile!("triangulate_density_map_append"); - - let mut mesh = TriMesh3d::default(); - let subdomain = subdomain.clone(); - - assert!( - subdomain.subdomain_grid().cells_per_dim().iter().all(|&n_cells| n_cells > I::one() + I::one()), - "Interpolation procedure with stitching support only works on grids & subdomains with more than 2 cells in each dimension!" - ); - - let (marching_cubes_data, boundary_data) = construct_mc_input_with_stitching_data( - &subdomain, + let marching_cubes_data = construct_mc_input( + grid, &density_map, iso_surface_threshold, &mut mesh.vertices, ); - triangulate_with_criterion( - &subdomain, - marching_cubes_data, - &mut mesh, - TriangulationSkipBoundaryCells, - DebugTriangleGenerator, - )?; - - Ok(SurfacePatch { - mesh, - subdomain, - data: boundary_data, - stitching_level: 0, - }) + triangulate(marching_cubes_data, mesh)?; + Ok(()) } /// Checks the consistency of the mesh (currently checks for holes, non-manifold edges and vertices) and returns a string with debug information in case of problems @@ -397,15 +331,12 @@ fn test_interpolate_cell_data() { let mut sparse_data = new_map(); - let marching_cubes_data = { - let subdomain = DummySubdomain::new(&grid); - construct_mc_input( - &subdomain, - &sparse_data.clone().into(), - iso_surface_threshold, - &mut trimesh.vertices, - ) - }; + let marching_cubes_data = construct_mc_input( + &grid, + &sparse_data.clone().into(), + iso_surface_threshold, + &mut trimesh.vertices, + ); assert_eq!(trimesh.vertices.len(), 0); assert_eq!(marching_cubes_data.cell_data.len(), 0); @@ -425,15 +356,12 @@ fn test_interpolate_cell_data() { sparse_data.insert(grid.flatten_point_index_array(&ijk), val); } - let marching_cubes_data = { - let subdomain = DummySubdomain::new(&grid); - construct_mc_input( - &subdomain, - &sparse_data.clone().into(), - iso_surface_threshold, - &mut trimesh.vertices, - ) - }; + let marching_cubes_data = construct_mc_input( + &grid, + &sparse_data.clone().into(), + iso_surface_threshold, + &mut trimesh.vertices, + ); assert_eq!(marching_cubes_data.cell_data.len(), 1); // Check that the correct number of vertices was created diff --git a/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs b/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs index 67357bc..863da93 100644 --- a/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs +++ b/splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs @@ -1,67 +1,27 @@ -use crate::marching_cubes::stitching::{collect_boundary_cell_data, BoundaryData}; use crate::marching_cubes::{CellData, MarchingCubesInput, RelativeToThreshold}; -use crate::topology::{Axis, DirectedAxisArray}; -use crate::uniform_grid::{CellIndex, GridBoundaryFaceFlags, PointIndex, Subdomain}; -use crate::{profile, DensityMap, Index, MapType, Real}; +use crate::uniform_grid::CellIndex; +use crate::{profile, DensityMap, Index, MapType, Real, UniformGrid}; use log::trace; use nalgebra::Vector3; -/// Trait used by [interpolate_points_to_cell_data_generic] to filter out points and edges during construction of iso-surface vertices in narrow-band cells -trait DensityMapFilter> { - /// Returns whether the given point should be considered for the density map to cell data conversion - fn process_point( - &mut self, - density_map: &DensityMap, - subdomain: &S, - flat_point_index: I, - subdomain_point: &PointIndex, - point_value: R, - ) -> bool; - - /// Returns whether the given edge should be considered for the density map to cell data conversion - fn process_edge( - &mut self, - density_map: &DensityMap, - subdomain: &S, - flat_point_index: I, - subdomain_point: &PointIndex, - flat_neighbor_index: I, - subdomain_neighbor: &PointIndex, - ) -> bool; -} - -/// Cell data interpolation filter that accepts all points and edges -struct IdentityDensityMapFilter; - -/// Cell data interpolation filter that skips one boundary layer of cells of the domain and builds the boundary density maps -struct SkipBoundaryLayerFilter { - boundary_density_maps: DirectedAxisArray>, -} - -/// Cell data interpolation filter for processing the stitching domain between two reconstructed surface patches -struct StitchingDomainNarrowBandFilter { - stitching_axis: Axis, -} - /// Returns the marching cubes input data for the narrow band of a single contiguous domain without support for stitching -pub(crate) fn construct_mc_input>( - subdomain: &S, +pub(crate) fn construct_mc_input( + grid: &UniformGrid, density_map: &DensityMap, iso_surface_threshold: R, vertices: &mut Vec>, ) -> MarchingCubesInput { let mut marching_cubes_data = MarchingCubesInput::default(); - let _ = interpolate_points_to_cell_data_generic::( - subdomain, + let _ = interpolate_points_to_cell_data_generic::( + grid, density_map, iso_surface_threshold, vertices, &mut marching_cubes_data, - IdentityDensityMapFilter, ); update_cell_data_threshold_flags( - subdomain, + grid, density_map, iso_surface_threshold, true, @@ -71,69 +31,6 @@ pub(crate) fn construct_mc_input>( marching_cubes_data } -/// Returns marching cubes input data for the narrow band of the subdomain and collects all data on the boundaries required for stitching of this domain to neighboring domains -pub(crate) fn construct_mc_input_with_stitching_data>( - subdomain: &S, - density_map: &DensityMap, - iso_surface_threshold: R, - vertices: &mut Vec>, -) -> (MarchingCubesInput, DirectedAxisArray>) { - let mut marching_cubes_data = MarchingCubesInput::default(); - let boundary_filter = interpolate_points_to_cell_data_generic( - subdomain, - density_map, - iso_surface_threshold, - vertices, - &mut marching_cubes_data, - SkipBoundaryLayerFilter::new(), - ); - - update_cell_data_threshold_flags( - subdomain, - density_map, - iso_surface_threshold, - true, - &mut marching_cubes_data, - ); - - let mut boundary_density_maps = boundary_filter.into_inner(); - let mut boundary_cell_data = collect_boundary_cell_data(subdomain, &mut marching_cubes_data); - - let boundary_data = DirectedAxisArray::new_with(|axis| BoundaryData { - boundary_density_map: std::mem::take(boundary_density_maps.get_mut(axis)), - boundary_cell_data: std::mem::take(boundary_cell_data.get_mut(axis)), - }); - - (marching_cubes_data, boundary_data) -} - -/// Updates the marching cubes input data in the narrow band of a stitching domain between two reconstructed patches -pub(crate) fn update_mc_input_for_stitching_domain>( - subdomain: &S, - density_map: &DensityMap, - iso_surface_threshold: R, - stitching_axis: Axis, - vertices: &mut Vec>, - marching_cubes_input: &mut MarchingCubesInput, -) { - let _ = interpolate_points_to_cell_data_generic( - subdomain, - density_map, - iso_surface_threshold, - vertices, - marching_cubes_input, - StitchingDomainNarrowBandFilter::new(stitching_axis), - ); - - update_cell_data_threshold_flags( - subdomain, - density_map, - iso_surface_threshold, - false, - marching_cubes_input, - ); -} - /// Generates input data for performing the actual marching cubes triangulation /// /// This function interpolates iso-surface vertices for cells in the narrow band around the iso-surface. @@ -150,23 +47,14 @@ pub(crate) fn update_mc_input_for_stitching_domain, - F: DensityMapFilter, ->( - subdomain: &S, +fn interpolate_points_to_cell_data_generic( + grid: &UniformGrid, density_map: &DensityMap, iso_surface_threshold: R, vertices: &mut Vec>, marching_cubes_data: &mut MarchingCubesInput, - mut filter: F, -) -> F { - let grid = subdomain.global_grid(); - let subdomain_grid = subdomain.subdomain_grid(); - - profile!("interpolate_points_to_cell_data_skip_boundary"); +) { + profile!("interpolate_points_to_cell_data_generic"); trace!("Starting interpolation of cell data for marching cubes (excluding boundary layer)... (Input: {} existing vertices)", vertices.len()); // Map from flat cell index to all data that is required per cell for the marching cubes triangulation @@ -176,23 +64,7 @@ fn interpolate_points_to_cell_data_generic< { profile!("generate_iso_surface_vertices"); density_map.for_each(|flat_point_index, point_value| { - let global_point = grid.try_unflatten_point_index(flat_point_index).unwrap(); - let point = subdomain - .map_point(&global_point) - .expect("Point cannot be mapped into subdomain."); - - // First check if the point should be processed, this can also be used to store certain - // values in a separate storage, e.g. to filter out points on the boundary but store them - // in a separate list to process them later. - if !filter.process_point( - density_map, - subdomain, - flat_point_index, - &point, - point_value, - ) { - return; - } + 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. @@ -206,14 +78,12 @@ fn interpolate_points_to_cell_data_generic< return; } - let neighborhood = subdomain_grid.get_point_neighborhood(&point); + let neighborhood = grid.get_point_neighborhood(&point); // Iterate over all neighbors of the point to find edges crossing the iso-surface for neighbor_edge in neighborhood.neighbor_edge_iter() { let neighbor = neighbor_edge.neighbor_index(); - // Get flat index of neighbor on global grid - let global_neighbor = subdomain.inv_map_point(neighbor).unwrap(); - let flat_neighbor_index = grid.flatten_point_index(&global_neighbor); + let flat_neighbor_index = grid.flatten_point_index(&neighbor); // Try to read out the function value at the neighboring point let neighbor_value = if let Some(v) = density_map.get(flat_neighbor_index) { @@ -231,21 +101,10 @@ fn interpolate_points_to_cell_data_generic< continue; } - if !filter.process_edge( - density_map, - subdomain, - flat_point_index, - &point, - flat_neighbor_index, - neighbor, - ) { - continue; - } - // Interpolate iso-surface vertex on the edge let alpha = (iso_surface_threshold - point_value) / (neighbor_value - point_value); - let point_coords = subdomain_grid.point_coordinates(&point); - let neighbor_coords = subdomain_grid.point_coordinates(neighbor); + let point_coords = grid.point_coordinates(&point); + let neighbor_coords = grid.point_coordinates(neighbor); let interpolated_coords = (point_coords) * (R::one() - alpha) + neighbor_coords * alpha; @@ -256,13 +115,8 @@ fn interpolate_points_to_cell_data_generic< // Store the data required for the marching cubes triangulation for // each cell adjacent to the edge crossing the iso-surface. // This includes the above/below iso-surface flags and the interpolated vertex index. - for cell in subdomain_grid - .cells_adjacent_to_edge(&neighbor_edge) - .iter() - .flatten() - { - let global_cell = subdomain.inv_map_cell(cell).unwrap(); - let flat_cell_index = grid.flatten_cell_index(&global_cell); + for cell in grid.cells_adjacent_to_edge(&neighbor_edge).iter().flatten() { + let flat_cell_index = grid.flatten_cell_index(&cell); let cell_data_entry = cell_data .entry(flat_cell_index) @@ -292,13 +146,11 @@ fn interpolate_points_to_cell_data_generic< cell_data.len(), vertices.len() ); - - filter } /// Loops through all corner vertices in the given marching cubes input and updates the above/below threshold flags -fn update_cell_data_threshold_flags>( - subdomain: &S, +fn update_cell_data_threshold_flags( + grid: &UniformGrid, density_map: &DensityMap, iso_surface_threshold: R, skip_points_above_threshold: bool, @@ -319,8 +171,6 @@ fn update_cell_data_threshold_flags>( // consisting of missing points and points below the surface. profile!("relative_to_threshold_postprocessing"); - let grid = subdomain.global_grid(); - let update_corner_flags = |cell: &CellIndex, local_point_index: usize, flag: &mut RelativeToThreshold| { // Otherwise try to look up its value and potentially mark it as above the threshold @@ -369,204 +219,3 @@ fn update_cell_data_threshold_flags>( } } } - -impl> DensityMapFilter for IdentityDensityMapFilter { - #[inline(always)] - fn process_point( - &mut self, - _density_map: &DensityMap, - _subdomain: &S, - _flat_point_index: I, - _subdomain_point: &PointIndex, - _point_value: R, - ) -> bool { - return true; - } - - #[inline(always)] - fn process_edge( - &mut self, - _density_map: &DensityMap, - _subdomain: &S, - _flat_point_index: I, - _subdomain_point: &PointIndex, - _flat_neighbor_index: I, - _subdomain_neighbor: &PointIndex, - ) -> bool { - return true; - } -} - -impl SkipBoundaryLayerFilter { - /// Constructs a new empty filter - fn new() -> Self { - Self { - boundary_density_maps: Default::default(), - } - } - - /// Consumes self and returns the collected boundary density maps - fn into_inner(self) -> DirectedAxisArray> { - self.boundary_density_maps - } -} - -impl> DensityMapFilter - for SkipBoundaryLayerFilter -{ - #[inline(always)] - fn process_point( - &mut self, - density_map: &DensityMap, - subdomain: &S, - flat_point_index: I, - subdomain_point: &PointIndex, - point_value: R, - ) -> bool { - let grid = subdomain.global_grid(); - let subdomain_grid = subdomain.subdomain_grid(); - - // Skip points directly at the boundary but add them to the respective boundary density map - { - let point_boundary_flags = - GridBoundaryFaceFlags::classify_point(subdomain_grid, &subdomain_point); - if !point_boundary_flags.is_empty() { - // Insert the point into each boundary density map it belongs to - for boundary in point_boundary_flags.iter_individual() { - let boundary_map = self.boundary_density_maps.get_mut(&boundary); - boundary_map.insert(flat_point_index, point_value); - - // Also insert second row neighbor, if present - if let Some(flat_neighbor_index) = subdomain_grid - // Get neighbor in subdomain - .get_point_neighbor(&subdomain_point, boundary.opposite()) - // Map neighbor from subdomain to global grid - .and_then(|neighbor| subdomain.inv_map_point(&neighbor)) - // Flatten on global grid - .map(|global_neighbor| grid.flatten_point_index(&global_neighbor)) - { - if let Some(density_value) = density_map.get(flat_neighbor_index) { - boundary_map.insert(flat_neighbor_index, density_value); - } - } - } - // Skip this point for interpolation - return false; - } - } - - // Otherwise the point can be processed - return true; - } - - #[inline(always)] - fn process_edge( - &mut self, - _density_map: &DensityMap, - subdomain: &S, - _flat_point_index: I, - _subdomain_point: &PointIndex, - _flat_neighbor_index: I, - subdomain_neighbor: &PointIndex, - ) -> bool { - let subdomain_grid = subdomain.subdomain_grid(); - - // Check if the neighbor is on the boundary of the subdomain - let point_boundary_flags = - GridBoundaryFaceFlags::classify_point(subdomain_grid, subdomain_neighbor); - let point_is_on_outer_boundary = !point_boundary_flags.is_empty(); - - // Skip edges that go into the boundary layer - if point_is_on_outer_boundary { - return false; - } - - return true; - } -} - -impl StitchingDomainNarrowBandFilter { - /// Constructs a new filter for the narrow band of a stitching domain - fn new(stitching_axis: Axis) -> Self { - StitchingDomainNarrowBandFilter { stitching_axis } - } - - /// Detects points that are on the positive/negative side of the stitching domain, along the stitching axis - fn point_is_on_stitching_surface>( - &self, - subdomain: &S, - p: &PointIndex, - ) -> bool { - let subdomain_grid = subdomain.subdomain_grid(); - - let index = p.index(); - index[self.stitching_axis.dim()] == I::zero() - || index[self.stitching_axis.dim()] - == subdomain_grid.points_per_dim()[self.stitching_axis.dim()] - I::one() - } - - /// Detects points that are on a boundary other than the stitching surfaces - fn point_is_outside_stitching_volume>( - &self, - subdomain: &S, - p: &PointIndex, - ) -> bool { - let subdomain_grid = subdomain.subdomain_grid(); - - let index = p.index(); - self.stitching_axis - .orthogonal_axes() - .iter() - .copied() - .any(|axis| { - index[axis.dim()] == I::zero() - || index[axis.dim()] == subdomain_grid.points_per_dim()[axis.dim()] - I::one() - }) - } -} - -impl> DensityMapFilter - for StitchingDomainNarrowBandFilter -{ - #[inline(always)] - fn process_point( - &mut self, - _density_map: &DensityMap, - subdomain: &S, - _flat_point_index: I, - subdomain_point: &PointIndex, - _point_value: R, - ) -> bool { - // Skip points on the outside of the stitching domain (except if they are on the stitching surface) - if self.point_is_outside_stitching_volume(subdomain, subdomain_point) { - return false; - } - - true - } - - #[inline(always)] - fn process_edge( - &mut self, - _density_map: &DensityMap, - subdomain: &S, - _flat_point_index: I, - subdomain_point: &PointIndex, - _flat_neighbor_index: I, - subdomain_neighbor: &PointIndex, - ) -> bool { - // Skip edges that are on the stitching surface (were already triangulated by the patches) - if self.point_is_on_stitching_surface(subdomain, subdomain_point) - && self.point_is_on_stitching_surface(subdomain, subdomain_neighbor) - { - return false; - } - - // Skip edges that go out of the stitching domain - if self.point_is_outside_stitching_volume(subdomain, subdomain_neighbor) { - return false; - } - - true - } -} diff --git a/splashsurf_lib/src/marching_cubes/stitching.rs b/splashsurf_lib/src/marching_cubes/stitching.rs deleted file mode 100644 index 14962b0..0000000 --- a/splashsurf_lib/src/marching_cubes/stitching.rs +++ /dev/null @@ -1,594 +0,0 @@ -use crate::marching_cubes::narrow_band_extraction::update_mc_input_for_stitching_domain; -use crate::marching_cubes::triangulation::{ - triangulate_with_criterion, DebugTriangleGenerator, TriangulationStitchingInterior, -}; -use crate::marching_cubes::{CellData, MarchingCubesError, MarchingCubesInput}; -use crate::mesh::TriMesh3d; -use crate::topology::{Axis, DirectedAxis, DirectedAxisArray, Direction}; -use crate::uniform_grid::{GridBoundaryFaceFlags, OwningSubdomainGrid, Subdomain, UniformGrid}; -use crate::{profile, Index, MapType, Real, ReconstructionError}; -use log::{debug, trace}; - -/// Stitches the two given surface patches by triangulating the domain between them -pub(crate) fn stitch_surface_patches( - iso_surface_threshold: R, - stitching_axis: Axis, - mut negative_side: SurfacePatch, - mut positive_side: SurfacePatch, -) -> Result, ReconstructionError> { - profile!("stitch_surface_patches"); - - assert_eq!( - negative_side.subdomain.global_grid(), - positive_side.subdomain.global_grid(), - "The global grid of the two subdomains that should be stitched is not identical!" - ); - - // Take out boundary data to satisfy borrow checker - let mut negative_data = std::mem::take(&mut negative_side.data); - let mut positive_data = std::mem::take(&mut positive_side.data); - - let global_grid = negative_side.subdomain.global_grid(); - - debug!( - "Stitching patches orthogonal to {:?}-axis. (-) side (offset: {:?}, cells_per_dim: {:?}, stitching_level: {:?}), (+) side (offset: {:?}, cells_per_dim: {:?}, stitching_level: {:?})", - stitching_axis, - negative_side.subdomain.subdomain_offset(), - negative_side.subdomain.subdomain_grid().cells_per_dim(), - negative_side.stitching_level, - positive_side.subdomain.subdomain_offset(), - positive_side.subdomain.subdomain_grid().cells_per_dim(), - positive_side.stitching_level, - ); - - // Construct domain for the triangulation of the boundary layer between the sides - let stitching_subdomain = compute_stitching_domain( - stitching_axis, - global_grid, - &negative_side.subdomain, - &positive_side.subdomain, - ); - - // Merge the two input meshes structures and get vertex offset for all vertices of the positive side - let (mut output_mesh, negative_vertex_offset, positive_vertex_offset) = { - let mut negative_mesh = std::mem::take(&mut negative_side.mesh); - let mut positive_mesh = std::mem::take(&mut positive_side.mesh); - - if negative_mesh.vertices.len() > positive_mesh.vertices.len() { - let positive_vertex_offset = negative_mesh.vertices.len(); - negative_mesh.append(&mut positive_mesh); - (negative_mesh, None, Some(positive_vertex_offset)) - } else { - let negative_vertex_offset = positive_mesh.vertices.len(); - positive_mesh.append(&mut negative_mesh); - (positive_mesh, Some(negative_vertex_offset), None) - } - }; - - // Merge the boundary data at the stitching boundaries of the two patches - let merged_boundary_data = { - // On the negative side we need the data of its positive boundary and vice versa - let negative_data = std::mem::take( - negative_data.get_mut(&DirectedAxis::new(stitching_axis, Direction::Positive)), - ); - let positive_data = std::mem::take( - positive_data.get_mut(&DirectedAxis::new(stitching_axis, Direction::Negative)), - ); - - // Merge the boundary layer density and cell data maps of the two sides - merge_boundary_data( - &stitching_subdomain, - &negative_side.subdomain, - negative_data, - negative_vertex_offset, - &positive_side.subdomain, - positive_data, - positive_vertex_offset, - ) - }; - - let BoundaryData { - boundary_density_map, - boundary_cell_data, - } = merged_boundary_data; - - let mut marching_cubes_input = MarchingCubesInput { - cell_data: boundary_cell_data, - }; - - // Perform marching cubes on the stitching domain - let mut boundary_cell_data = { - update_mc_input_for_stitching_domain( - &stitching_subdomain, - &boundary_density_map.into(), - iso_surface_threshold, - stitching_axis, - &mut output_mesh.vertices, - &mut marching_cubes_input, - ); - - // Collect the boundary cell data of the stitching domain - let boundary_cell_data = - collect_boundary_cell_data(&stitching_subdomain, &marching_cubes_input); - - triangulate_with_criterion( - &stitching_subdomain, - marching_cubes_input, - &mut output_mesh, - TriangulationStitchingInterior::new(stitching_axis), - DebugTriangleGenerator, - ) - .map_err(MarchingCubesError::from)?; - - boundary_cell_data - }; - - // Get domain for the whole stitched domain - let output_subdomain_grid = compute_stitching_result_domain( - stitching_axis, - global_grid, - &negative_side.subdomain, - &positive_side.subdomain, - ); - - // Merge all remaining boundary data - let output_boundary_data = DirectedAxisArray::new_with(|&directed_axis| { - // The positive and negative sides of the result domain can be taken directly from the inputs - if directed_axis == stitching_axis.with_direction(Direction::Negative) { - let data = std::mem::take(negative_data.get_mut(&directed_axis)); - data.map_to_domain( - &output_subdomain_grid, - &negative_side.subdomain, - negative_vertex_offset, - ) - } else if directed_axis == stitching_axis.with_direction(Direction::Positive) { - let data = std::mem::take(positive_data.get_mut(&directed_axis)); - data.map_to_domain( - &output_subdomain_grid, - &positive_side.subdomain, - positive_vertex_offset, - ) - } else { - // Otherwise, they have to be merged first - let mut merged_data = merge_boundary_data( - &output_subdomain_grid, - &negative_side.subdomain, - std::mem::take(negative_data.get_mut(&directed_axis)), - negative_vertex_offset, - &positive_side.subdomain, - std::mem::take(positive_data.get_mut(&directed_axis)), - positive_vertex_offset, - ); - - // Map cell indices from stitching domain to result domain and append to cell data map - let current_boundary_cell_data = - std::mem::take(boundary_cell_data.get_mut(&directed_axis)); - for (flat_cell_index, cell_data) in current_boundary_cell_data { - let global_cell = global_grid - .try_unflatten_cell_index(flat_cell_index) - .unwrap(); - - // Skip cells not part of the output grid - if output_subdomain_grid.map_cell(&global_cell).is_none() { - continue; - } - - merged_data - .boundary_cell_data - .entry(flat_cell_index) - // Merge with existing cell data - .and_modify(|existing_cell_data| { - // Should be fine to just replace these values as they will be overwritten anyway in the next stitching process - existing_cell_data.corner_above_threshold = - cell_data.corner_above_threshold; - // For the iso-surface vertices we need the union - for (existing_vertex, new_vertex) in existing_cell_data - .iso_surface_vertices - .iter_mut() - .zip(cell_data.iso_surface_vertices.iter()) - { - if existing_vertex != new_vertex { - assert!( - existing_vertex.is_none(), - "Overwriting already existing vertex. This is a bug." - ); - *existing_vertex = *new_vertex - } - } - }) - // Or insert new cell data - .or_insert_with(|| cell_data.clone()); - } - - merged_data - } - }); - - Ok(SurfacePatch { - subdomain: output_subdomain_grid, - mesh: output_mesh, - data: output_boundary_data, - stitching_level: negative_side - .stitching_level - .max(positive_side.stitching_level), - }) -} - -/// A surface patch representing a local part of a larger surface reconstruction -#[derive(Clone, Debug)] -pub(crate) struct SurfacePatch { - /// The local surface mesh of this side - pub(crate) mesh: TriMesh3d, - /// The subdomain of this local mesh - pub(crate) subdomain: OwningSubdomainGrid, - /// All additional data required for stitching - pub(crate) data: DirectedAxisArray>, - /// The maximum number of times parts of this patch where stitched together - pub(crate) stitching_level: usize, -} - -impl SurfacePatch { - /// Creates an empty surface patch for the given subdomain - pub(crate) fn new_empty(subdomain: OwningSubdomainGrid) -> Self { - Self { - mesh: Default::default(), - subdomain, - data: Default::default(), - stitching_level: 0, - } - } -} - -/// Stitching data per boundary -#[derive(Clone, Default, Debug)] -pub(crate) struct BoundaryData { - /// The density map for all vertices of this boundary - pub(crate) boundary_density_map: MapType, - /// The cell data for all cells of this boundary - pub(crate) boundary_cell_data: MapType, -} - -/// Extracts the cell data of all cells on the boundary of the subdomain -#[inline(never)] -pub(crate) fn collect_boundary_cell_data>( - subdomain: &S, - input: &MarchingCubesInput, -) -> DirectedAxisArray> { - let mut boundary_cell_data: DirectedAxisArray> = Default::default(); - - for (&flat_cell_index, cell_data) in &input.cell_data { - let global_cell = subdomain - .global_grid() - .try_unflatten_cell_index(flat_cell_index) - .unwrap(); - let local_cell = subdomain.map_cell(&global_cell).unwrap(); - - // Check which grid boundary faces this cell is part of - let cell_grid_boundaries = - GridBoundaryFaceFlags::classify_cell(subdomain.subdomain_grid(), &local_cell); - // Only process cells that are part of some boundary - if !cell_grid_boundaries.is_empty() { - for boundary in cell_grid_boundaries.iter_individual() { - boundary_cell_data - .get_mut(&boundary) - .insert(flat_cell_index, cell_data.clone()); - } - } - } - - boundary_cell_data -} - -impl BoundaryData { - /// Maps this boundary data to another domain by converting all indices to the new subdomain - fn map_to_domain( - mut self, - target_domain: &OwningSubdomainGrid, - source_domain: &OwningSubdomainGrid, - source_offset: Option, - ) -> Self { - assert_eq!( - target_domain.global_grid(), - source_domain.global_grid(), - "The global grid of target and source domain has to match!" - ); - let grid = target_domain.global_grid(); - - // Process density map - { - let mut points_to_remove = Vec::new(); - for flat_point_index in self.boundary_density_map.keys().copied() { - let global_point = grid.try_unflatten_point_index(flat_point_index).unwrap(); - - // If point is not part of the target domain, it is marked to be removed - if target_domain.map_point(&global_point).is_none() { - points_to_remove.push(flat_point_index) - } - } - - // Remove all points not part of the target domain - for flat_point_index in points_to_remove { - self.boundary_density_map.remove(&flat_point_index); - } - } - - // Process cell data - { - let mut cells_to_remove = Vec::new(); - for flat_cell_index in self.boundary_cell_data.keys().copied() { - let global_cell = grid.try_unflatten_cell_index(flat_cell_index).unwrap(); - - // If cell is not part of the target domain, it is marked to be removed - if target_domain.map_cell(&global_cell).is_none() { - cells_to_remove.push(flat_cell_index); - } - } - - // Remove all cells not part of the target domain - for flat_cell_index in cells_to_remove { - self.boundary_cell_data.remove(&flat_cell_index); - } - } - - // Apply vertex offset if required - self.apply_cell_data_vertex_offset(source_offset); - - self - } - - /// Merges another boundary data object into this object - fn merge_with( - mut self, - target_domain: &OwningSubdomainGrid, - mut other: BoundaryData, - other_domain: &OwningSubdomainGrid, - other_vertex_offset: Option, - ) -> Self { - assert_eq!( - target_domain.global_grid(), - other_domain.global_grid(), - "The global grid of target and source domain has to match!" - ); - let grid = target_domain.global_grid(); - - // Apply vertex offset if required - other.apply_cell_data_vertex_offset(other_vertex_offset); - - let BoundaryData { - boundary_density_map: other_boundary_density_map, - boundary_cell_data: other_boundary_cell_data, - } = other; - - // Process density map - for (flat_point_index, density_contribution) in other_boundary_density_map { - let global_point = grid.try_unflatten_point_index(flat_point_index).unwrap(); - - // Skip points that are not part of the target domain - if target_domain.map_point(&global_point).is_none() { - continue; - } - - // TODO: For a proper average we would have to keep track how often this point is merged with other domains - // and use an incremental average over this number of contributions - self.boundary_density_map - .entry(flat_point_index) - // Compute average with existing value - .and_modify(|density| { - *density += density_contribution; - *density /= R::one() + R::one(); - }) - // Or just insert the new value - .or_insert(density_contribution); - } - - // Process cell data - for (flat_cell_index, cell_data) in other_boundary_cell_data { - let global_cell = grid.try_unflatten_cell_index(flat_cell_index).unwrap(); - - // Skip points that are not part of the target domain - if target_domain.map_cell(&global_cell).is_none() { - continue; - } - - self.boundary_cell_data - .entry(flat_cell_index) - // The cell data interpolation function should only populate cells that are part of their subdomain - .and_modify(|_| { - panic!("Merge conflict: there is duplicate cell data for this cell index") - }) - // Otherwise insert the additional cell data - .or_insert(cell_data); - } - - self - } - - /// Adds an offset to all iso-surface vertex indices stored in the cell data map - fn apply_cell_data_vertex_offset(&mut self, vertex_offset: Option) { - // Apply the vertex offset - if let Some(vertex_offset) = vertex_offset { - for v in self - .boundary_cell_data - .iter_mut() - // Iterate over all vertices in the cell data - .flat_map(|(_, cell_data)| cell_data.iso_surface_vertices.iter_mut()) - // Skip missing vertices - .flatten() - { - *v += vertex_offset; - } - } - } -} - -/// Merges boundary such that only density values and cell data in the result subdomain are part of the result -fn merge_boundary_data( - target_subdomain: &OwningSubdomainGrid, - negative_subdomain: &OwningSubdomainGrid, - negative_data: BoundaryData, - negative_vertex_offset: Option, - positive_subdomain: &OwningSubdomainGrid, - positive_data: BoundaryData, - positive_vertex_offset: Option, -) -> BoundaryData { - trace!("Merging boundary data. Size of containers: (-) side (density map: {}, cell data map: {}), (+) side (density map: {}, cell data map: {})", negative_data.boundary_density_map.len(), negative_data.boundary_cell_data.len(), positive_data.boundary_density_map.len(), positive_data.boundary_cell_data.len()); - - let negative_len = - negative_data.boundary_density_map.len() + negative_data.boundary_cell_data.len(); - let positive_len = - positive_data.boundary_density_map.len() + positive_data.boundary_cell_data.len(); - - let merged_boundary_data = if negative_len > positive_len { - let merged_boundary_data = negative_data.map_to_domain( - target_subdomain, - negative_subdomain, - negative_vertex_offset, - ); - merged_boundary_data.merge_with( - target_subdomain, - positive_data, - positive_subdomain, - positive_vertex_offset, - ) - } else { - let merged_boundary_data = positive_data.map_to_domain( - target_subdomain, - positive_subdomain, - positive_vertex_offset, - ); - merged_boundary_data.merge_with( - target_subdomain, - negative_data, - negative_subdomain, - negative_vertex_offset, - ) - }; - - trace!("Finished merging boundary data. Size of containers: result (density map: {}, cell data map: {})", merged_boundary_data.boundary_density_map.len(), merged_boundary_data.boundary_cell_data.len()); - - merged_boundary_data -} - -/// Computes the [OwningSubdomainGrid] for stitching region between the two sides that has to be triangulated -fn compute_stitching_domain( - stitching_axis: Axis, - global_grid: &UniformGrid, - negative_subdomain: &OwningSubdomainGrid, - positive_subdomain: &OwningSubdomainGrid, -) -> OwningSubdomainGrid { - // Ensure that global grids are equivalent - assert_eq!( - negative_subdomain.global_grid(), - global_grid, - "The global grid of the two subdomains that should be stitched is not identical!" - ); - assert_eq!( - positive_subdomain.global_grid(), - global_grid, - "The global grid of the two subdomains that should be stitched is not identical!" - ); - - // Check that the two domains actually meet - { - // Starting at the offset of the negative subdomain and going along the stitching axis... - let lower_corner_end = stitching_axis - .with_direction(Direction::Positive) - .checked_apply_step_ijk( - negative_subdomain.subdomain_offset(), - negative_subdomain.subdomain_grid().cells_per_dim(), - ) - .expect("Index type out of range?"); - - // ...we should arrive at the lower corner of the positive side - assert_eq!( - lower_corner_end, - *positive_subdomain.subdomain_offset(), - "The two subdomains that should be stitched do not meet directly!" - ); - } - - // Get the number of cells of the stitching domain - let n_cells_per_dim = { - let mut n_cells_per_dim_neg = negative_subdomain.subdomain_grid().cells_per_dim().clone(); - let mut n_cells_per_dim_pos = positive_subdomain.subdomain_grid().cells_per_dim().clone(); - - // Between the two subdomains are only two layers of cells - n_cells_per_dim_neg[stitching_axis.dim()] = I::one() + I::one(); - n_cells_per_dim_pos[stitching_axis.dim()] = I::one() + I::one(); - - // Ensure that the stitching domain is identical from both sides - assert_eq!( - n_cells_per_dim_neg, n_cells_per_dim_pos, - "The cross sections of the two subdomains that should be stitched is not identical!" - ); - - let n_cells_per_dim = n_cells_per_dim_neg; - n_cells_per_dim - }; - - // Obtain the index of the lower corner of the stitching domain - let stitching_grid_offset = { - let axis_index = stitching_axis.dim(); - - // Start at offset of negative domain - let mut stitching_grid_offset = negative_subdomain.subdomain_offset().clone(); - - // Go to the end of the negative domain along the stitching axis - stitching_grid_offset[axis_index] += - negative_subdomain.subdomain_grid().cells_per_dim()[axis_index]; - // Subtract one because the stitching domain starts in the boundary layer - stitching_grid_offset[axis_index] -= I::one(); - stitching_grid_offset - }; - // Get coordinates of offset point - let lower_corner_coords = global_grid.point_coordinates_array(&stitching_grid_offset); - - // Build the grid for the stitching domain - let stitching_grid = UniformGrid::new( - &lower_corner_coords, - &n_cells_per_dim, - global_grid.cell_size(), - ) - .expect("Unable to construct stitching domain grid"); - - trace!( - "Constructed domain for stitching. offset: {:?}, cells_per_dim: {:?}", - stitching_grid_offset, - n_cells_per_dim - ); - - OwningSubdomainGrid::new(global_grid.clone(), stitching_grid, stitching_grid_offset) -} - -/// Computes the [OwningSubdomainGrid] for the final combined domain of the two sides -fn compute_stitching_result_domain( - stitching_axis: Axis, - global_grid: &UniformGrid, - negative_subdomain: &OwningSubdomainGrid, - positive_subdomain: &OwningSubdomainGrid, -) -> OwningSubdomainGrid { - // Get the number of cells of the result domain by adding all cells in stitching direction - let n_cells_per_dim = { - let length_neg = negative_subdomain.subdomain_grid().cells_per_dim()[stitching_axis.dim()]; - let length_pos = positive_subdomain.subdomain_grid().cells_per_dim()[stitching_axis.dim()]; - - let mut n_cells_per_dim = negative_subdomain.subdomain_grid().cells_per_dim().clone(); - n_cells_per_dim[stitching_axis.dim()] = length_neg + length_pos; - - n_cells_per_dim - }; - - // Construct the grid - let subdomain_grid = UniformGrid::new( - &negative_subdomain.subdomain_grid().aabb().min(), - &n_cells_per_dim, - global_grid.cell_size(), - ) - .expect("Unable to construct stitching domain grid"); - - OwningSubdomainGrid::new( - global_grid.clone(), - subdomain_grid, - negative_subdomain.subdomain_offset().clone(), - ) -} diff --git a/splashsurf_lib/src/marching_cubes/triangulation.rs b/splashsurf_lib/src/marching_cubes/triangulation.rs index 056a32d..9e39597 100644 --- a/splashsurf_lib/src/marching_cubes/triangulation.rs +++ b/splashsurf_lib/src/marching_cubes/triangulation.rs @@ -1,13 +1,9 @@ use crate::marching_cubes::marching_cubes_lut::marching_cubes_triangulation_iter; use crate::marching_cubes::{CellData, MarchingCubesInput}; use crate::mesh::TriMesh3d; -use crate::topology::Axis; -use crate::uniform_grid::{DummySubdomain, GridBoundaryFaceFlags, Subdomain, UniformGrid}; use crate::{profile, Index, Real}; use anyhow::Context; use log::trace; -use nalgebra::Vector3; -use std::marker::PhantomData; use thiserror::Error as ThisError; /// Error enum for the marching cubes triangulation stage @@ -22,68 +18,13 @@ pub enum TriangulationError { ), } -/// Trait that is used by the marching cubes [triangulate_with_criterion] function to query whether a cell should be triangulated -pub(crate) trait TriangulationCriterion> { - /// Returns whether the given cell should be triangulated - fn triangulate_cell(&self, subdomain: &S, flat_cell_index: I, cell_data: &CellData) -> bool; -} - -/// An identity triangulation criterion that accepts all cells -pub(crate) struct TriangulationIdentityCriterion; -/// A triangulation criterion that ensures that every cell is part of the subdomain but skips one layer of boundary cells -pub(crate) struct TriangulationSkipBoundaryCells; -/// A triangulation criterion that ensures that only the interior of the stitching domain is triangulated (boundary layer except in stitching direction is skipped) -pub(crate) struct TriangulationStitchingInterior { - stitching_axis: Axis, -} - -/// Trait that is used by the marching cubes [triangulate_with_criterion] function to convert a marching cubes triangulation to actual triangle-vertex connectivity -pub(crate) trait TriangleGenerator> { - fn triangle_connectivity( - &self, - subdomain: &S, - flat_cell_index: I, - cell_data: &CellData, - edge_indices: [i32; 3], - ) -> Result<[usize; 3], anyhow::Error>; -} - -/// Maps the edges indices directly to the vertex indices in the cell data -pub(crate) struct DefaultTriangleGenerator; -/// Tries to map the edge indices to the vertex indices in the cell data, returns an error with debug information if vertices are missing -pub(crate) struct DebugTriangleGenerator; - /// Converts the marching cubes input cell data into a triangle surface mesh, appends triangles to existing mesh #[inline(never)] pub(crate) fn triangulate( input: MarchingCubesInput, mesh: &mut TriMesh3d, ) -> Result<(), TriangulationError> { - triangulate_with_criterion( - &DummySubdomain::new(&UniformGrid::new_zero()), - input, - mesh, - TriangulationIdentityCriterion, - DebugTriangleGenerator, - ) -} - -/// Converts the marching cubes input cell data into a triangle surface mesh, appends triangles to existing mesh with custom criterion to filter out cells during triangulation -#[inline(never)] -pub(crate) fn triangulate_with_criterion< - I: Index, - R: Real, - S: Subdomain, - C: TriangulationCriterion, - G: TriangleGenerator, ->( - subdomain: &S, - input: MarchingCubesInput, - mesh: &mut TriMesh3d, - triangulation_criterion: C, - triangle_generator: G, -) -> Result<(), TriangulationError> { - profile!("triangulate_with_criterion"); + profile!("triangulate"); let MarchingCubesInput { cell_data } = input; @@ -95,18 +36,12 @@ pub(crate) fn triangulate_with_criterion< ); // Triangulate affected cells - for (&flat_cell_index, cell_data) in &cell_data { - // Skip cells that don't fulfill triangulation criterion - if !triangulation_criterion.triangulate_cell(subdomain, flat_cell_index, cell_data) { - continue; - } - + for (&_flat_cell_index, cell_data) in &cell_data { // TODO: Replace `are_vertices_above_unchecked` with something that can return an error for triangle in marching_cubes_triangulation_iter(&cell_data.are_vertices_above_unchecked()) { // TODO: Allow user to set option to skip invalid triangles? - let global_triangle = triangle_generator - .triangle_connectivity(subdomain, flat_cell_index, cell_data, triangle) + let global_triangle = get_triangle(cell_data, triangle) .map_err(|e| TriangulationError::TriangleConnectivityError(e))?; mesh.triangles.push(global_triangle); } @@ -121,143 +56,6 @@ pub(crate) fn triangulate_with_criterion< Ok(()) } -/// Forwards to the wrapped triangulation criterion but first makes some assertions on the cell data -struct DebugTriangulationCriterion< - I: Index, - R: Real, - S: Subdomain, - C: TriangulationCriterion, -> { - triangulation_criterion: C, - phantom: PhantomData<(I, R, S)>, -} - -impl> TriangulationCriterion - for TriangulationIdentityCriterion -{ - #[inline(always)] - fn triangulate_cell(&self, _: &S, _: I, _: &CellData) -> bool { - true - } -} - -impl> TriangulationCriterion - for TriangulationSkipBoundaryCells -{ - #[inline(always)] - fn triangulate_cell(&self, subdomain: &S, flat_cell_index: I, _: &CellData) -> bool { - let global_cell = subdomain - .global_grid() - .try_unflatten_cell_index(flat_cell_index) - .unwrap(); - let local_cell = subdomain - .map_cell(&global_cell) - .expect("Cell is not part of the subdomain"); - let cell_grid_boundaries = - GridBoundaryFaceFlags::classify_cell(subdomain.subdomain_grid(), &local_cell); - - return cell_grid_boundaries.is_empty(); - } -} - -impl TriangulationStitchingInterior { - pub(crate) fn new(stitching_axis: Axis) -> Self { - Self { stitching_axis } - } -} - -impl> TriangulationCriterion - for TriangulationStitchingInterior -{ - #[inline(always)] - fn triangulate_cell(&self, subdomain: &S, flat_cell_index: I, _: &CellData) -> bool { - let global_cell = subdomain - .global_grid() - .try_unflatten_cell_index(flat_cell_index) - .unwrap(); - - let local_cell = subdomain - .map_cell(&global_cell) - .expect("Cell is not part of the subdomain"); - - let subdomain_grid = subdomain.subdomain_grid(); - let index = local_cell.index(); - - // Skip only boundary cells in directions orthogonal to the stitching axis - !self.stitching_axis.orthogonal_axes().iter().any(|&axis| { - index[axis.dim()] == I::zero() - || index[axis.dim()] == subdomain_grid.cells_per_dim()[axis.dim()] - I::one() - }) - } -} - -impl, C: TriangulationCriterion> - DebugTriangulationCriterion -{ - #[allow(unused)] - fn new(triangulation_criterion: C) -> Self { - Self { - triangulation_criterion, - phantom: Default::default(), - } - } -} - -impl, C: TriangulationCriterion> - TriangulationCriterion for DebugTriangulationCriterion -{ - #[inline(always)] - fn triangulate_cell(&self, subdomain: &S, flat_cell_index: I, cell_data: &CellData) -> bool { - assert!( - !cell_data - .corner_above_threshold - .iter() - .any(|c| c.is_indeterminate()), - "{}", - cell_debug_string(subdomain, flat_cell_index, cell_data) - ); - - self.triangulation_criterion - .triangulate_cell(subdomain, flat_cell_index, cell_data) - } -} - -impl> TriangleGenerator - for DefaultTriangleGenerator -{ - #[inline(always)] - fn triangle_connectivity( - &self, - _subdomain: &S, - _flat_cell_index: I, - cell_data: &CellData, - edge_indices: [i32; 3], - ) -> Result<[usize; 3], anyhow::Error> { - // Note: If the one of the following expect calls causes a panic, it is probably because - // a cell was added improperly to the marching cubes input, e.g. a cell was added to the - // cell data map that is not part of the domain. This results in only those edges of the cell - // that are neighboring to the domain having correct iso surface vertices. The remaining - // edges would have missing iso-surface vertices and overall this results in an invalid triangulation - // - // If this happens, it's a bug in the cell data map generation. - get_triangle(cell_data, edge_indices) - } -} - -impl> TriangleGenerator for DebugTriangleGenerator { - #[inline(always)] - fn triangle_connectivity( - &self, - subdomain: &S, - flat_cell_index: I, - cell_data: &CellData, - edge_indices: [i32; 3], - ) -> Result<[usize; 3], anyhow::Error> { - get_triangle(cell_data, edge_indices) - .with_context(|| cell_debug_string(subdomain, flat_cell_index, cell_data)) - } -} - /// Helper function that extracts vertex indices from the [`CellData`] for the triangle with the given edge indices fn get_triangle(cell_data: &CellData, edge_indices: [i32; 3]) -> Result<[usize; 3], anyhow::Error> { let [edge_idx_0, edge_idx_1, edge_idx_2] = edge_indices; @@ -295,32 +93,3 @@ fn get_triangle(cell_data: &CellData, edge_indices: [i32; 3]) -> Result<[usize; })?, ]) } - -/// Helper function that returns a formatted string to debug triangulation failures -fn cell_debug_string>( - subdomain: &S, - flat_cell_index: I, - cell_data: &CellData, -) -> String { - let global_cell_index = subdomain - .global_grid() - .try_unflatten_cell_index(flat_cell_index) - .expect("Failed to get cell index"); - - let point_index = subdomain - .global_grid() - .get_point(*global_cell_index.index()) - .expect("Unable to get point index of cell"); - let cell_center = subdomain.global_grid().point_coordinates(&point_index) - + &Vector3::repeat(subdomain.global_grid().cell_size().times_f64(0.5)); - - format!( - "Unable to construct triangle for cell {:?}, with center coordinates {:?} and edge length {}.\n{:?}\nStitching domain: (offset: {:?}, cells_per_dim: {:?})", - global_cell_index.index(), - cell_center, - subdomain.global_grid().cell_size(), - cell_data, - subdomain.subdomain_offset(), - subdomain.subdomain_grid().cells_per_dim(), - ) -} diff --git a/splashsurf_lib/src/octree.rs b/splashsurf_lib/src/octree.rs deleted file mode 100644 index b8679cf..0000000 --- a/splashsurf_lib/src/octree.rs +++ /dev/null @@ -1,1252 +0,0 @@ -//! Octree for spatially partitioning particle sets - -use crate::generic_tree::*; -use crate::marching_cubes::SurfacePatch; -use crate::mesh::{HexMesh3d, MeshAttribute, MeshWithData, TriMesh3d}; -use crate::topology::{Axis, Direction}; -use crate::uniform_grid::{PointIndex, UniformGrid}; -use crate::utils::{ChunkSize, ParallelPolicy}; -use crate::{ - marching_cubes, new_map, profile, Aabb3d, GridConstructionError, Index, MapType, Real, - ReconstructionError, -}; -use arrayvec::ArrayVec; -use log::info; -use nalgebra::Vector3; -use octant_helper::{HalfspaceFlags, Octant, OctantAxisDirections}; -use rayon::prelude::*; -use smallvec::SmallVec; -use split_criterion::{default_split_criterion, LeafSplitCriterion}; -use std::cell::RefCell; -use std::sync::atomic::{AtomicUsize, Ordering}; -use thread_local::ThreadLocal; - -// TODO: Make margin an Option - -/// Criterion used for the subdivision of the spatial decomposition of the particle collection -#[derive(Clone, Debug)] -pub enum SubdivisionCriterion { - /// Perform octree subdivision until an upper limit of particles is reached per chunk, automatically chosen based on number of threads - MaxParticleCountAuto, - /// Perform octree subdivision until an upper limit of particles is reached per chunk, based on the given fixed number of particles - MaxParticleCount(usize), -} - -/// Data structure for octree based spatial subdivision of particles sets, for tree iteration/visitation use the [`root`](Self::root) [`OctreeNode`] -#[derive(Clone, Debug)] -pub struct Octree { - /// Root node of the tree - root: OctreeNode, - /// Counter for assigning ids to subdivided nodes - next_id: usize, -} - -/// Represents a node in the octree hierarchy and stores child nodes, implements tree iteration/visitation from the [`generic_tree`](crate::generic_tree) module -#[derive(Clone, Debug)] -pub struct OctreeNode { - /// Id of the node used to identify it for debugging - id: usize, - /// All child nodes of this octree node - children: ArrayVec, 8>, - /// Lower corner point of the octree node on the background grid - min_corner: PointIndex, - /// Upper corner point of the octree node on the background grid - max_corner: PointIndex, - /// AABB of the octree node - aabb: Aabb3d, - /// Additional data associated to this octree node - data: NodeData, -} - -impl TreeNode for OctreeNode { - /// Returns a slice of all child nodes - fn children(&self) -> &[Box] { - self.children.as_slice() - } -} - -impl TreeNodeMut for OctreeNode { - /// Returns a mutable slice of all child nodes - fn children_mut(&mut self) -> &mut [Box] { - self.children.as_mut_slice() - } -} - -/// Optional data that may be stored in [`OctreeNode`]s -#[derive(Clone, Debug)] -pub enum NodeData { - /// Empty variant - None, - /// Storage for a set of SPH particles - ParticleSet(ParticleSet), - /// A patch that was already meshed - SurfacePatch(SurfacePatchWrapper), -} - -impl Default for NodeData { - /// Returns an empty data instance - fn default() -> Self { - Self::None - } -} - -/// Stores the particle ids and the number of ghost particles inside an octree leaf -#[derive(Clone, Debug)] -pub struct ParticleSet { - // The particles belonging to this set - pub particles: OctreeNodeParticleStorage, - // Number of ghost particles in this particle set - pub ghost_particle_count: usize, -} - -/// Wrapper for an internal `SurfacePatch` to avoid leaking too much implementation details -#[derive(Clone, Debug)] -pub struct SurfacePatchWrapper { - pub(crate) patch: SurfacePatch, -} - -impl From> for SurfacePatchWrapper { - fn from(patch: SurfacePatch) -> Self { - Self { patch } - } -} - -impl SurfacePatchWrapper { - pub fn mesh(&self) -> &TriMesh3d { - &self.patch.mesh - } -} - -type OctreeNodeParticleStorage = SmallVec<[usize; 6]>; - -impl Octree { - /// Creates a new octree with a single leaf node containing all vertices - pub fn new(grid: &UniformGrid, n_particles: usize) -> Self { - Self { - root: OctreeNode::new_root(grid, n_particles), - next_id: 0, - } - } - - /// Create a new octree and perform subdivision with the specified margin - /// - /// The margin is used to assign ghost particles to octree nodes. Each octant resulting - /// from the subdivision gets assigned all particles that are directly inside it plus all - /// particles from its parent that are within the given margin around the octant. - pub fn new_subdivided( - grid: &UniformGrid, - particle_positions: &[Vector3], - subdivision_criterion: SubdivisionCriterion, - margin: R, - enable_multi_threading: bool, - enable_stitching: bool, - ) -> Self { - let mut tree = Octree::new(&grid, particle_positions.len()); - - if enable_multi_threading { - tree.par_subdivide_recursively_margin( - grid, - particle_positions, - subdivision_criterion, - margin, - enable_stitching, - ); - } else { - tree.subdivide_recursively_margin( - grid, - particle_positions, - subdivision_criterion, - margin, - enable_stitching, - ); - } - - tree - } - - /// Returns a reference to the root node of the octree - pub fn root(&self) -> &OctreeNode { - &self.root - } - - /// Returns a mutable reference to the root node of the octree - pub fn root_mut(&mut self) -> &mut OctreeNode { - &mut self.root - } - - /// Subdivide the octree recursively using the given splitting criterion and a margin to add ghost particles - pub fn subdivide_recursively_margin( - &mut self, - grid: &UniformGrid, - particle_positions: &[Vector3], - subdivision_criterion: SubdivisionCriterion, - margin: R, - enable_stitching: bool, - ) { - profile!("octree subdivide_recursively_margin"); - - let split_criterion = default_split_criterion( - subdivision_criterion, - particle_positions.len(), - enable_stitching, - ); - - let next_id = AtomicUsize::new(0); - self.root.visit_mut_bfs(|node| { - // Stop recursion if split criterion is not fulfilled - if !split_criterion.split_leaf(node) { - return; - } - - // Perform one octree split on the node - node.subdivide_with_margin(grid, particle_positions, margin, &next_id); - }); - self.next_id = next_id.into_inner(); - } - - /// Subdivide the octree recursively and in parallel using the given splitting criterion and a margin to add ghost particles - pub fn par_subdivide_recursively_margin( - &mut self, - grid: &UniformGrid, - particle_positions: &[Vector3], - subdivision_criterion: SubdivisionCriterion, - margin: R, - enable_stitching: bool, - ) { - profile!("octree subdivide_recursively_margin_par"); - - let split_criterion = default_split_criterion( - subdivision_criterion, - particle_positions.len(), - enable_stitching, - ); - let parallel_policy = ParallelPolicy::default(); - - let next_id = AtomicUsize::new(0); - let visitor = { - let next_id = &next_id; - move |node: &mut OctreeNode| { - // Stop recursion if split criterion is not fulfilled - if !split_criterion.split_leaf(node) { - return; - } - - // Perform one octree split on the leaf - if node - .data - .particle_set() - .expect("Node is not a leaf") - .particles - .len() - < parallel_policy.min_task_size - { - node.subdivide_with_margin(grid, particle_positions, margin, &next_id); - } else { - node.par_subdivide_with_margin( - grid, - particle_positions, - margin, - ¶llel_policy, - &next_id, - ); - } - } - }; - - self.root.par_visit_mut_bfs(visitor); - self.next_id = next_id.into_inner(); - } - - /// Constructs a hex mesh visualizing the cells of the octree, may contain hanging and duplicate vertices as cells are not connected - pub fn hexmesh( - &self, - grid: &UniformGrid, - only_non_empty: bool, - ) -> MeshWithData> { - profile!("convert octree into hexmesh"); - - let mut mesh = HexMesh3d { - vertices: Vec::new(), - cells: Vec::new(), - }; - - let mut ids = Vec::new(); - self.root.dfs_iter().for_each(|node| { - if node.children().is_empty() { - if only_non_empty - && node - .data() - .particle_set() - .map(|ps| ps.particles.is_empty()) - .unwrap_or(true) - { - return; - } - - let lower_coords = grid.point_coordinates(&node.min_corner); - let upper_coords = grid.point_coordinates(&node.max_corner); - - let vertices = vec![ - lower_coords, - Vector3::new(upper_coords[0], lower_coords[1], lower_coords[2]), - Vector3::new(upper_coords[0], upper_coords[1], lower_coords[2]), - Vector3::new(lower_coords[0], upper_coords[1], lower_coords[2]), - Vector3::new(lower_coords[0], lower_coords[1], upper_coords[2]), - Vector3::new(upper_coords[0], lower_coords[1], upper_coords[2]), - upper_coords, - Vector3::new(lower_coords[0], upper_coords[1], upper_coords[2]), - ]; - - let offset = mesh.vertices.len(); - let cell = [ - offset + 0, - offset + 1, - offset + 2, - offset + 3, - offset + 4, - offset + 5, - offset + 6, - offset + 7, - ]; - - mesh.vertices.extend(vertices); - mesh.cells.push(cell); - ids.push(node.id as u64); - } - }); - - assert_eq!(mesh.cells.len(), ids.len()); - MeshWithData::new(mesh).with_cell_data(MeshAttribute::new("node_id".to_string(), ids)) - } -} - -impl OctreeNode { - pub fn new( - id: usize, - min_corner: PointIndex, - max_corner: PointIndex, - aabb: Aabb3d, - ) -> Self { - Self::with_data(id, min_corner, max_corner, aabb, NodeData::None) - } - - fn new_root(grid: &UniformGrid, n_particles: usize) -> Self { - let n_points = grid.points_per_dim(); - let min_point = [I::zero(), I::zero(), I::zero()]; - let max_point = [ - n_points[0] - I::one(), - n_points[1] - I::one(), - n_points[2] - I::one(), - ]; - - Self::with_data( - 0, - grid.get_point(min_point) - .expect("Cannot get lower corner of grid"), - grid.get_point(max_point) - .expect("Cannot get upper corner of grid"), - grid.aabb().clone(), - NodeData::new_particle_set((0..n_particles).collect::>(), 0), - ) - } - - fn with_data( - id: usize, - min_corner: PointIndex, - max_corner: PointIndex, - aabb: Aabb3d, - data: NodeData, - ) -> Self { - Self { - id, - children: Default::default(), - min_corner, - max_corner, - aabb, - data, - } - } - - /// Returns the id of the node - pub fn id(&self) -> usize { - self.id - } - - /// Returns a reference to the data stored in the node - pub fn data(&self) -> &NodeData { - &self.data - } - - /// Returns a mutable reference to the data stored in the node - pub(crate) fn data_mut(&mut self) -> &mut NodeData { - &mut self.data - } - - /// Returns the [`PointIndex`] of the lower corner of the octree node - pub fn min_corner(&self) -> &PointIndex { - &self.min_corner - } - - /// Returns the [`PointIndex`] of the upper corner of the octree node - pub fn max_corner(&self) -> &PointIndex { - &self.max_corner - } - - /// Returns the AABB represented by this octree node - pub fn aabb(&self) -> &Aabb3d { - &self.aabb - } - - /// Constructs a [`UniformGrid`] that represents the domain of this octree node - pub fn grid( - &self, - min: &Vector3, - cell_size: R, - ) -> Result, GridConstructionError> { - let min_corner = self.min_corner.index(); - let max_corner = self.max_corner.index(); - - let n_cells_per_dim = [ - max_corner[0] - min_corner[0], - max_corner[1] - min_corner[1], - max_corner[2] - min_corner[2], - ]; - - UniformGrid::new(min, &n_cells_per_dim, cell_size) - } - - /// Performs a subdivision of this node while considering a margin for "ghost particles" around each octant - pub fn subdivide_with_margin( - &mut self, - grid: &UniformGrid, - particle_positions: &[Vector3], - margin: R, - next_id: &AtomicUsize, - ) { - // Convert node body from Leaf to Children - if let NodeData::ParticleSet(particle_set) = self.data.take() { - let particles = particle_set.particles; - - // Obtain the point used as the octree split/pivot point - let split_point = get_split_point(grid, &self.min_corner, &self.max_corner) - .expect("Failed to get split point of octree node"); - let split_coordinates = grid.point_coordinates(&split_point); - - let mut halfspace_flags = vec![HalfspaceFlags::empty(); particles.len()]; - let mut counters: [usize; 8] = [0, 0, 0, 0, 0, 0, 0, 0]; - let mut non_ghost_counters: [usize; 8] = [0, 0, 0, 0, 0, 0, 0, 0]; - - // Classify all particles of this leaf into the halfspaces relative to the split point - assert_eq!(particles.len(), halfspace_flags.len()); - for (particle_idx, particle_halfspace_flags) in - particles.iter().copied().zip(halfspace_flags.iter_mut()) - { - let pos = particle_positions[particle_idx]; - let relative_pos = pos - split_coordinates; - - let is_ghost_particle = !self.aabb.contains_point(&pos); - - // Check what the main octant (without margin) of the particle is to count ghost particles later - if !is_ghost_particle { - let main_octant: Octant = OctantAxisDirections::classify(&relative_pos).into(); - non_ghost_counters[main_octant as usize] += 1; - } - - // Classify into all halfspaces with margin - { - *particle_halfspace_flags = - HalfspaceFlags::classify_with_margin(&relative_pos, margin); - - // Increase the counter of each octant that contains the current particle - HalfspaceFlags::all_unique_octants() - .iter() - .zip(counters.iter_mut()) - .filter(|(octant, _)| particle_halfspace_flags.contains(**octant)) - .for_each(|(_, counter)| { - *counter += 1; - }); - } - } - - // Construct the node for each octant - let mut children = ArrayVec::new(); - for (¤t_octant, (&octant_particle_count, &octant_non_ghost_count)) in - Octant::all() - .iter() - .zip(counters.iter().zip(non_ghost_counters.iter())) - { - let current_octant_dir = OctantAxisDirections::from(current_octant); - let current_octant_flags = HalfspaceFlags::from(current_octant); - - let min_corner = current_octant_dir - .combine_point_index(grid, &self.min_corner, &split_point) - .expect("Failed to get corner point of octree subcell"); - let max_corner = current_octant_dir - .combine_point_index(grid, &split_point, &self.max_corner) - .expect("Failed to get corner point of octree subcell"); - - let child_aabb = Aabb3d::new( - grid.point_coordinates(&min_corner), - grid.point_coordinates(&max_corner), - ); - - let mut octant_particles = SmallVec::with_capacity(octant_particle_count); - octant_particles.extend( - particles - .iter() - .copied() - .zip(halfspace_flags.iter()) - // Skip particles from other octants - .filter(|(_, &particle_i_halfspaces)| { - particle_i_halfspaces.contains(current_octant_flags) - }) - .map(|(particle_i, _)| particle_i), - ); - assert_eq!(octant_particles.len(), octant_particle_count); - - let child = Box::new(OctreeNode::with_data( - next_id.fetch_add(1, Ordering::SeqCst), - min_corner, - max_corner, - child_aabb, - NodeData::new_particle_set( - octant_particles, - octant_particle_count - octant_non_ghost_count, - ), - )); - - children.push(child); - } - - // Assign new children to the current node - self.children = children; - } else { - panic!("Only nodes with ParticleSet data can be subdivided"); - }; - } - - /// Parallel subdivision of this node while considering a margin for "ghost particles" around each octant - pub fn par_subdivide_with_margin( - &mut self, - grid: &UniformGrid, - particle_positions: &[Vector3], - margin: R, - parallel_policy: &ParallelPolicy, - next_id: &AtomicUsize, - ) { - // Convert node body from Leaf to Children - if let NodeData::ParticleSet(particle_set) = self.data.take() { - let particles = particle_set.particles; - - // Obtain the point used as the octree split/pivot point - let split_point = get_split_point(grid, &self.min_corner, &self.max_corner) - .expect("Failed to get split point of octree node"); - let split_coordinates = grid.point_coordinates(&split_point); - - let mut octant_flags = vec![HalfspaceFlags::empty(); particles.len()]; - - // Initial values for the thread local counters - let zeros = || ([0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0]); - let zeros_cell = || RefCell::new(zeros()); - - let tl_counters = ThreadLocal::new(); - let chunk_size = ChunkSize::new(parallel_policy, particles.len()).chunk_size; - - // Classify all particles of this leaf into its octants - assert_eq!(particles.len(), octant_flags.len()); - particles - .par_chunks(chunk_size) - .zip(octant_flags.par_chunks_mut(chunk_size)) - .for_each(|(idx_chunk, flags_chunk)| { - // Obtain references to the thread-local counters - let mut counters_borrow_mut = tl_counters.get_or(zeros_cell).borrow_mut(); - let counters_ref_mut = &mut *counters_borrow_mut; - let (counters, non_ghost_counters) = - (&mut counters_ref_mut.0, &mut counters_ref_mut.1); - - idx_chunk - .iter() - .copied() - .zip(flags_chunk.iter_mut()) - .for_each(|(particle_idx, particle_octant_flags)| { - let pos = particle_positions[particle_idx]; - let relative_pos = particle_positions[particle_idx] - split_coordinates; - - let is_ghost_particle = !self.aabb.contains_point(&pos); - - // Check what the main octant of the particle is (to count ghost particles) - if !is_ghost_particle { - let main_octant: Octant = - OctantAxisDirections::classify(&relative_pos).into(); - non_ghost_counters[main_octant as usize] += 1; - } - - // Classify into all octants with margin - { - *particle_octant_flags = - HalfspaceFlags::classify_with_margin(&relative_pos, margin); - - // Increase the counter of each octant that contains the current particle - HalfspaceFlags::all_unique_octants() - .iter() - .zip(counters.iter_mut()) - .filter(|(octant, _)| particle_octant_flags.contains(**octant)) - .for_each(|(_, counter)| { - *counter += 1; - }); - } - }) - }); - - // Sum up all thread local counter arrays - let (counters, non_ghost_counters) = tl_counters.into_iter().fold( - zeros(), - |(mut counters_acc, mut non_ghost_counters_acc), counter_cell| { - let (counters, non_ghost_counters) = counter_cell.into_inner(); - for i in 0..8 { - counters_acc[i] += counters[i]; - non_ghost_counters_acc[i] += non_ghost_counters[i]; - } - (counters_acc, non_ghost_counters_acc) - }, - ); - - // TODO: Would be nice to collect directly into a ArrayVec but that doesn't seem to be possible - // (at least without some unsafe magic with uninit) - let mut children = Vec::with_capacity(8); - // Construct the octree node for each octant - Octant::all() - .par_iter() - .zip(counters.par_iter().zip(non_ghost_counters.par_iter())) - .map( - |(¤t_octant, (&octant_particle_count, &octant_non_ghost_count))| { - let current_octant_dir = OctantAxisDirections::from(current_octant); - let current_octant_flags = HalfspaceFlags::from(current_octant); - - let min_corner = current_octant_dir - .combine_point_index(grid, &self.min_corner, &split_point) - .expect("Failed to get corner point of octree subcell"); - let max_corner = current_octant_dir - .combine_point_index(grid, &split_point, &self.max_corner) - .expect("Failed to get corner point of octree subcell"); - - let child_aabb = Aabb3d::new( - grid.point_coordinates(&min_corner), - grid.point_coordinates(&max_corner), - ); - - let mut octant_particles = SmallVec::with_capacity(octant_particle_count); - octant_particles.extend( - particles - .iter() - .copied() - .zip(octant_flags.iter()) - // Skip particles from other octants - .filter(|(_, &particle_i_octant)| { - particle_i_octant.contains(current_octant_flags) - }) - .map(|(particle_i, _)| particle_i), - ); - assert_eq!(octant_particles.len(), octant_particle_count); - - let child = Box::new(OctreeNode::with_data( - next_id.fetch_add(1, Ordering::SeqCst), - min_corner, - max_corner, - child_aabb, - NodeData::new_particle_set( - octant_particles, - octant_particle_count - octant_non_ghost_count, - ), - )); - - child - }, - ) - .collect_into_vec(&mut children); - - // Assign children to this node - self.children = children.into_iter().collect::>(); - } else { - panic!("Only nodes with ParticleSet data can be subdivided"); - }; - } - - fn stitch_children_orthogonal_to( - &mut self, - children_map: &mut MapType>, - stitching_axis: Axis, - iso_surface_threshold: R, - ) -> Result<(), ReconstructionError> { - profile!("stitch_children_orthogonal_to"); - - for mut octant in OctantAxisDirections::all().iter().copied() { - // Iterate over every octant pair only once - if octant.direction(stitching_axis).is_positive() { - continue; - } - - // First try to get negative side, it might not exist because children were already merged before to another octant in the map - let negative_side = if let Some(negative_patch) = children_map.remove(&octant) { - negative_patch - } else { - continue; - }; - - // If the negative side on the stitching axis exists, the positive side must also exist - octant.set_direction(stitching_axis, Direction::Positive); - let positive_side = children_map.remove(&octant).expect("Child node missing!"); - - let stitched_patch = marching_cubes::stitch_surface_patches( - iso_surface_threshold, - stitching_axis, - negative_side, - positive_side, - )?; - - // Add stitched surface back to map, setting the direction of the octant of the stitched patch to positive - children_map.insert(octant, stitched_patch); - } - - Ok(()) - } - - /// Stitches together the [`SurfacePatch`]es stored in the children of this node if this is the direct parent of only leaf nodes - pub(crate) fn stitch_surface_patches( - &mut self, - iso_surface_threshold: R, - ) -> Result<(), ReconstructionError> { - profile!("stitch_surface_patches"); - - // If this node has no children there is nothing to stitch - if self.children.is_empty() { - panic!("A node can only be stitched if it has children!"); - } - - // Don't try to stitch if there are children that still have children - for child in self.children.iter() { - if !child.children.is_empty() { - panic!("A node can only be stitched if all children are leaf nodes!"); - } - } - - let mut children_map: MapType<_, SurfacePatch> = { - let old_children = std::mem::take(&mut self.children); - - let mut children_map = new_map(); - for (child, octant) in old_children.into_iter().zip(Octant::all().iter().copied()) { - let octant_directions = OctantAxisDirections::from(octant); - children_map.insert( - octant_directions, - child - .data - .into_surface_patch() - .expect("Surface patch missing!") - .patch, - ); - } - - children_map - }; - - self.stitch_children_orthogonal_to(&mut children_map, Axis::X, iso_surface_threshold)?; - self.stitch_children_orthogonal_to(&mut children_map, Axis::Y, iso_surface_threshold)?; - self.stitch_children_orthogonal_to(&mut children_map, Axis::Z, iso_surface_threshold)?; - - assert_eq!( - children_map.len(), - 1, - "After stitching, there should be only one child left." - ); - - for (_, mut stitched_patch) in children_map.into_iter() { - stitched_patch.stitching_level += 1; - self.data = NodeData::SurfacePatch(stitched_patch.into()); - break; - } - - assert!( - self.children.is_empty(), - "After stitching, the node should not have any children." - ); - - Ok(()) - } -} - -impl NodeData { - fn new_particle_set>( - particles: P, - ghost_particle_count: usize, - ) -> Self { - let particles = particles.into(); - NodeData::ParticleSet(ParticleSet { - particles, - ghost_particle_count, - }) - } - - /// Returns a reference to the contained particle set if it contains one - pub fn particle_set(&self) -> Option<&ParticleSet> { - if let Self::ParticleSet(particle_set) = self { - Some(particle_set) - } else { - None - } - } - - /// Returns a reference to the contained surface patch if it contains one - pub fn surface_patch(&self) -> Option<&SurfacePatchWrapper> { - if let Self::SurfacePatch(surface_patch) = self { - Some(surface_patch) - } else { - None - } - } - - /// Consumes self and returns the ParticleSet if it contained one - pub fn into_particle_set(self) -> Option { - if let Self::ParticleSet(particle_set) = self { - Some(particle_set) - } else { - None - } - } - - /// Consumes self and returns the SurfacePatch if it contained one - pub fn into_surface_patch(self) -> Option> { - if let Self::SurfacePatch(surface_patch) = self { - Some(surface_patch) - } else { - None - } - } - - /// Returns the stored data and leaves `None` in its place - pub fn take(&mut self) -> Self { - std::mem::take(self) - } - - /// Returns the stored data and leaves `None` in its place - pub fn replace(&mut self, new_data: Self) { - *self = new_data; - } -} - -/// Returns the [`PointIndex`] of the octree subdivision point for an [`OctreeNode`] with the given lower and upper points -fn get_split_point( - grid: &UniformGrid, - lower: &PointIndex, - upper: &PointIndex, -) -> Option> { - let two = I::one() + I::one(); - - let lower = lower.index(); - let upper = upper.index(); - - let mid_indices = [ - (lower[0] + upper[0]) / two, - (lower[1] + upper[1]) / two, - (lower[2] + upper[2]) / two, - ]; - - grid.get_point(mid_indices) -} - -mod split_criterion { - use super::*; - - /// Trait that is used by an octree to decide whether an octree node should be further split or subdivided - pub(super) trait LeafSplitCriterion { - /// Returns whether the specified node should be split - fn split_leaf(&self, node: &OctreeNode) -> bool; - } - - /// Split criterion that decides based on whether the number of non-ghost particles in a node is above a limit - pub(super) struct MaxNonGhostParticleLeafSplitCriterion { - max_particles: usize, - } - - impl MaxNonGhostParticleLeafSplitCriterion { - fn new(max_particles: usize) -> Self { - Self { max_particles } - } - } - - impl LeafSplitCriterion for MaxNonGhostParticleLeafSplitCriterion { - /// Returns true if the number of non-ghost particles in a node is above a limit - fn split_leaf(&self, node: &OctreeNode) -> bool { - match &node.data { - NodeData::ParticleSet(particle_set) => { - // Check if this leaf is already below the limit of particles per cell - return particle_set.particles.len() - particle_set.ghost_particle_count - > self.max_particles; - } - // Early out if called on a non-leaf node - _ => return false, - }; - } - } - - /// Split criterion that decides based on whether the node's extents are larger than 1 cell in all dimensions - pub(super) struct MinimumExtentSplitCriterion { - minimum_extent: I, - } - - impl MinimumExtentSplitCriterion { - fn new(minimum_extent: I) -> Self { - Self { minimum_extent } - } - } - - impl LeafSplitCriterion for MinimumExtentSplitCriterion { - /// Only returns true if a splitting of the node does not result in a node that is smaller than the allowed minimum extent - fn split_leaf(&self, node: &OctreeNode) -> bool { - let lower = node.min_corner.index(); - let upper = node.max_corner.index(); - - upper[0] - lower[0] >= self.minimum_extent + self.minimum_extent - && upper[1] - lower[1] >= self.minimum_extent + self.minimum_extent - && upper[2] - lower[2] >= self.minimum_extent + self.minimum_extent - } - } - - impl LeafSplitCriterion for (A, B) - where - A: LeafSplitCriterion, - B: LeafSplitCriterion, - { - fn split_leaf(&self, node: &OctreeNode) -> bool { - self.0.split_leaf(node) && self.1.split_leaf(node) - } - } - - pub(super) fn default_split_criterion( - subdivision_criterion: SubdivisionCriterion, - num_particles: usize, - enable_stitching: bool, - ) -> ( - MaxNonGhostParticleLeafSplitCriterion, - MinimumExtentSplitCriterion, - ) { - let particles_per_cell = match subdivision_criterion { - SubdivisionCriterion::MaxParticleCount(count) => count, - SubdivisionCriterion::MaxParticleCountAuto => { - ChunkSize::new(&ParallelPolicy::default(), num_particles) - .with_log("particles", "octree generation") - .chunk_size - } - }; - - info!( - "Building octree with at most {} particles per leaf", - particles_per_cell - ); - - ( - MaxNonGhostParticleLeafSplitCriterion::new(particles_per_cell), - MinimumExtentSplitCriterion::new(if enable_stitching { - I::one() + I::one() + I::one() - } else { - I::one() - }), - ) - } -} - -mod octant_helper { - use bitflags::bitflags; - use nalgebra::Vector3; - - use crate::topology::{Axis, Direction}; - use crate::uniform_grid::{PointIndex, UniformGrid}; - use crate::{Index, Real}; - - /// All octants of a 3D cartesian coordinate system - #[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Debug)] - #[repr(u8)] - pub enum Octant { - NegNegNeg = 0, - PosNegNeg = 1, - NegPosNeg = 2, - PosPosNeg = 3, - NegNegPos = 4, - PosNegPos = 5, - NegPosPos = 6, - PosPosPos = 7, - } - - /// Representation of a cartesian coordinate system octant using a direction along each coordinate axis - #[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Debug)] - pub struct OctantAxisDirections { - pub x_axis: Direction, - pub y_axis: Direction, - pub z_axis: Direction, - } - - bitflags! { - #[derive(Copy, Clone, Debug)] - pub struct HalfspaceFlags: u8 { - const X_NEG = 0b00000001; - const X_POS = 0b00000010; - const Y_NEG = 0b00000100; - const Y_POS = 0b00001000; - const Z_NEG = 0b00010000; - const Z_POS = 0b00100000; - - const NEG_NEG_NEG = Self::X_NEG.bits() | Self::Y_NEG.bits() | Self::Z_NEG.bits(); - const POS_NEG_NEG = Self::X_POS.bits() | Self::Y_NEG.bits() | Self::Z_NEG.bits(); - const NEG_POS_NEG = Self::X_NEG.bits() | Self::Y_POS.bits() | Self::Z_NEG.bits(); - const POS_POS_NEG = Self::X_POS.bits() | Self::Y_POS.bits() | Self::Z_NEG.bits(); - const NEG_NEG_POS = Self::X_NEG.bits() | Self::Y_NEG.bits() | Self::Z_POS.bits(); - const POS_NEG_POS = Self::X_POS.bits() | Self::Y_NEG.bits() | Self::Z_POS.bits(); - const NEG_POS_POS = Self::X_NEG.bits() | Self::Y_POS.bits() | Self::Z_POS.bits(); - const POS_POS_POS = Self::X_POS.bits() | Self::Y_POS.bits() | Self::Z_POS.bits(); - } - } - - impl HalfspaceFlags { - #[inline(always)] - pub const fn all_unique_octants() -> &'static [HalfspaceFlags] { - &ALL_UNIQUE_OCTANT_DIRECTION_FLAGS - } - - /// Classifies a point relative to zero into all halfspaces it belongs including a margin around the halfspace boundary - #[inline(always)] - pub fn classify_with_margin(point: &Vector3, margin: R) -> Self { - let mut flags = HalfspaceFlags::empty(); - flags.set(HalfspaceFlags::X_NEG, point.x < margin); - flags.set(HalfspaceFlags::X_POS, point.x > -margin); - flags.set(HalfspaceFlags::Y_NEG, point.y < margin); - flags.set(HalfspaceFlags::Y_POS, point.y > -margin); - flags.set(HalfspaceFlags::Z_NEG, point.z < margin); - flags.set(HalfspaceFlags::Z_POS, point.z > -margin); - flags - } - - #[inline(always)] - pub fn from_octant(octant: Octant) -> Self { - match octant { - Octant::NegNegNeg => Self::NEG_NEG_NEG, - Octant::PosNegNeg => Self::POS_NEG_NEG, - Octant::NegPosNeg => Self::NEG_POS_NEG, - Octant::PosPosNeg => Self::POS_POS_NEG, - Octant::NegNegPos => Self::NEG_NEG_POS, - Octant::PosNegPos => Self::POS_NEG_POS, - Octant::NegPosPos => Self::NEG_POS_POS, - Octant::PosPosPos => Self::POS_POS_POS, - } - } - - #[inline(always)] - pub fn from_directions(directions: OctantAxisDirections) -> Self { - let mut flags = HalfspaceFlags::empty(); - flags.set(HalfspaceFlags::X_NEG, directions.x_axis.is_negative()); - flags.set(HalfspaceFlags::X_POS, directions.x_axis.is_positive()); - flags.set(HalfspaceFlags::Y_NEG, directions.y_axis.is_negative()); - flags.set(HalfspaceFlags::Y_POS, directions.y_axis.is_positive()); - flags.set(HalfspaceFlags::Z_NEG, directions.z_axis.is_negative()); - flags.set(HalfspaceFlags::Z_POS, directions.z_axis.is_positive()); - flags - } - } - - impl From for HalfspaceFlags { - fn from(octant: Octant) -> Self { - Self::from_octant(octant) - } - } - - impl From for HalfspaceFlags { - fn from(directions: OctantAxisDirections) -> Self { - Self::from_directions(directions) - } - } - - const ALL_UNIQUE_OCTANT_DIRECTION_FLAGS: [HalfspaceFlags; 8] = [ - HalfspaceFlags::NEG_NEG_NEG, - HalfspaceFlags::POS_NEG_NEG, - HalfspaceFlags::NEG_POS_NEG, - HalfspaceFlags::POS_POS_NEG, - HalfspaceFlags::NEG_NEG_POS, - HalfspaceFlags::POS_NEG_POS, - HalfspaceFlags::NEG_POS_POS, - HalfspaceFlags::POS_POS_POS, - ]; - - impl OctantAxisDirections { - #[allow(dead_code)] - #[inline(always)] - pub const fn all() -> &'static [OctantAxisDirections; 8] { - &ALL_OCTANT_DIRECTIONS - } - - #[inline(always)] - pub const fn from_bool(x_positive: bool, y_positive: bool, z_positive: bool) -> Self { - Self { - x_axis: Direction::new_positive(x_positive), - y_axis: Direction::new_positive(y_positive), - z_axis: Direction::new_positive(z_positive), - } - } - - pub const fn from_octant(octant: Octant) -> Self { - match octant { - Octant::NegNegNeg => Self::from_bool(false, false, false), - Octant::PosNegNeg => Self::from_bool(true, false, false), - Octant::NegPosNeg => Self::from_bool(false, true, false), - Octant::PosPosNeg => Self::from_bool(true, true, false), - Octant::NegNegPos => Self::from_bool(false, false, true), - Octant::PosNegPos => Self::from_bool(true, false, true), - Octant::NegPosPos => Self::from_bool(false, true, true), - Octant::PosPosPos => Self::from_bool(true, true, true), - } - } - - pub fn direction(&self, axis: Axis) -> Direction { - match axis { - Axis::X => self.x_axis, - Axis::Y => self.y_axis, - Axis::Z => self.z_axis, - } - } - - pub fn set_direction(&mut self, axis: Axis, direction: Direction) { - match axis { - Axis::X => self.x_axis = direction, - Axis::Y => self.y_axis = direction, - Axis::Z => self.z_axis = direction, - } - } - - /// Classifies a point relative to zero into the corresponding octant - #[inline(always)] - pub fn classify(point: &Vector3) -> Self { - Self::from_bool( - point[0].is_positive(), - point[1].is_positive(), - point[2].is_positive(), - ) - } - - /// Combines two vectors by choosing between their components depending on the octant - pub fn combine_point_index( - &self, - grid: &UniformGrid, - lower: &PointIndex, - upper: &PointIndex, - ) -> Option> { - let lower = lower.index(); - let upper = upper.index(); - - let combined_index = [ - if self.x_axis.is_positive() { - upper[0] - } else { - lower[0] - }, - if self.y_axis.is_positive() { - upper[1] - } else { - lower[1] - }, - if self.z_axis.is_positive() { - upper[2] - } else { - lower[2] - }, - ]; - - grid.get_point(combined_index) - } - } - - impl From for OctantAxisDirections { - fn from(octant: Octant) -> Self { - Self::from_octant(octant) - } - } - - const ALL_OCTANT_DIRECTIONS: [OctantAxisDirections; 8] = [ - OctantAxisDirections::from_octant(Octant::NegNegNeg), - OctantAxisDirections::from_octant(Octant::PosNegNeg), - OctantAxisDirections::from_octant(Octant::NegPosNeg), - OctantAxisDirections::from_octant(Octant::PosPosNeg), - OctantAxisDirections::from_octant(Octant::NegNegPos), - OctantAxisDirections::from_octant(Octant::PosNegPos), - OctantAxisDirections::from_octant(Octant::NegPosPos), - OctantAxisDirections::from_octant(Octant::PosPosPos), - ]; - - impl Octant { - #[inline(always)] - pub const fn all() -> &'static [Octant; 8] { - &ALL_OCTANTS - } - - #[inline(always)] - pub const fn from_directions(directions: OctantAxisDirections) -> Self { - use Direction::*; - let OctantAxisDirections { - x_axis, - y_axis, - z_axis, - } = directions; - match (x_axis, y_axis, z_axis) { - (Negative, Negative, Negative) => Octant::NegNegNeg, - (Positive, Negative, Negative) => Octant::PosNegNeg, - (Negative, Positive, Negative) => Octant::NegPosNeg, - (Positive, Positive, Negative) => Octant::PosPosNeg, - (Negative, Negative, Positive) => Octant::NegNegPos, - (Positive, Negative, Positive) => Octant::PosNegPos, - (Negative, Positive, Positive) => Octant::NegPosPos, - (Positive, Positive, Positive) => Octant::PosPosPos, - } - } - } - - impl From for Octant { - fn from(directions: OctantAxisDirections) -> Self { - Self::from_directions(directions) - } - } - - const ALL_OCTANTS: [Octant; 8] = [ - Octant::NegNegNeg, - Octant::PosNegNeg, - Octant::NegPosNeg, - Octant::PosPosNeg, - Octant::NegNegPos, - Octant::PosNegPos, - Octant::NegPosPos, - Octant::PosPosPos, - ]; - - #[cfg(test)] - mod test_octant { - use super::*; - - #[test] - fn test_octant_iter_all_consistency() { - for (i, octant) in Octant::all().iter().copied().enumerate() { - assert_eq!(octant as usize, i); - assert_eq!(octant, unsafe { - std::mem::transmute::(i as u8) - }); - } - } - - #[test] - fn test_octant_directions_iter_all_consistency() { - assert_eq!(Octant::all().len(), OctantAxisDirections::all().len()); - for (octant, octant_directions) in Octant::all() - .iter() - .copied() - .zip(OctantAxisDirections::all().iter().copied()) - { - assert_eq!(octant, Octant::from(octant_directions)); - assert_eq!(octant_directions, OctantAxisDirections::from(octant)); - } - } - } -} diff --git a/splashsurf_lib/src/reconstruction.rs b/splashsurf_lib/src/reconstruction.rs index 56b5764..ff91006 100644 --- a/splashsurf_lib/src/reconstruction.rs +++ b/splashsurf_lib/src/reconstruction.rs @@ -1,12 +1,17 @@ -use anyhow::Context; -use log::info; -use nalgebra::Vector3; - use crate::dense_subdomains::{ compute_global_densities_and_neighbors, decomposition, initialize_parameters, reconstruction, stitching, subdomain_classification::GhostMarginClassifier, }; -use crate::{profile, Index, Parameters, Real, SurfaceReconstruction}; +use crate::mesh::TriMesh3d; +use crate::uniform_grid::UniformGrid; +use crate::workspace::LocalReconstructionWorkspace; +use crate::{ + density_map, marching_cubes, neighborhood_search, profile, Index, Parameters, Real, + ReconstructionError, SurfaceReconstruction, +}; +use anyhow::Context; +use log::{info, trace}; +use nalgebra::Vector3; /// Performs a surface reconstruction with a regular grid for domain decomposition pub(crate) fn reconstruct_surface_subdomain_grid<'a, I: Index, R: Real>( @@ -69,3 +74,129 @@ pub(crate) fn reconstruct_surface_subdomain_grid<'a, I: Index, R: Real>( } Ok(()) } + +/// Performs a global surface reconstruction without domain decomposition +pub(crate) fn reconstruct_surface_global<'a, I: Index, R: Real>( + particle_positions: &[Vector3], + parameters: &Parameters, + output_surface: &'a mut SurfaceReconstruction, +) -> Result<(), ReconstructionError> { + profile!("reconstruct_surface_global"); + + // Multiple local workspaces are only needed for processing different subdomains in parallel. + // However, in this global surface reconstruction without domain decomposition, each step in the + // reconstruction pipeline manages its memory on its own. + let mut workspace = output_surface + .workspace + .get_local_with_capacity(particle_positions.len()) + .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, + ); + } + } + + // Clear the current mesh, as reconstruction will be appended to output + output_surface.mesh.clear(); + // Perform global reconstruction without domain decomposition + reconstruct_single_surface_append( + &mut *workspace, + &output_surface.grid, + particle_positions, + parameters, + &mut output_surface.mesh, + )?; + + output_surface.particle_densities = Some(std::mem::take(&mut workspace.particle_densities)); + + Ok(()) +} + +/// Performs global neighborhood search and computes per particle densities +pub(crate) fn compute_particle_densities_and_neighbors( + grid: &UniformGrid, + particle_positions: &[Vector3], + parameters: &Parameters, + particle_neighbor_lists: &mut Vec>, + densities: &mut Vec, +) { + profile!("compute_particle_densities_and_neighbors"); + + let particle_rest_density = parameters.rest_density; + let particle_rest_volume = R::from_f64((4.0 / 3.0) * std::f64::consts::PI).unwrap() + * parameters.particle_radius.powi(3); + let particle_rest_mass = particle_rest_volume * particle_rest_density; + + trace!("Starting neighborhood search..."); + neighborhood_search::search_inplace::( + &grid.aabb(), + particle_positions, + parameters.compact_support_radius, + parameters.enable_multi_threading, + particle_neighbor_lists, + ); + + trace!("Computing particle densities..."); + density_map::compute_particle_densities_inplace::( + particle_positions, + particle_neighbor_lists.as_slice(), + parameters.compact_support_radius, + particle_rest_mass, + parameters.enable_multi_threading, + densities, + ); +} + +/// Reconstruct a surface, appends triangulation to the given mesh +pub(crate) fn reconstruct_single_surface_append<'a, I: Index, R: Real>( + workspace: &mut LocalReconstructionWorkspace, + grid: &UniformGrid, + particle_positions: &[Vector3], + parameters: &Parameters, + output_mesh: &'a mut TriMesh3d, +) -> Result<(), ReconstructionError> { + let particle_rest_density = parameters.rest_density; + let particle_rest_volume = R::from_f64((4.0 / 3.0) * std::f64::consts::PI).unwrap() + * parameters.particle_radius.powi(3); + let particle_rest_mass = particle_rest_volume * particle_rest_density; + + let particle_densities = { + compute_particle_densities_and_neighbors( + grid, + particle_positions, + parameters, + &mut workspace.particle_neighbor_lists, + &mut workspace.particle_densities, + ); + workspace.particle_densities.as_slice() + }; + + // Create a new density map, reusing memory with the workspace is bad for cache efficiency + // Alternatively one could reuse memory with a custom caching allocator + let mut density_map = Default::default(); + density_map::generate_sparse_density_map( + grid, + particle_positions, + particle_densities, + None, + particle_rest_mass, + parameters.compact_support_radius, + parameters.cube_size, + parameters.enable_multi_threading, + &mut density_map, + )?; + + marching_cubes::triangulate_density_map_append( + grid, + &density_map, + parameters.iso_surface_threshold, + output_mesh, + )?; + + Ok(()) +} diff --git a/splashsurf_lib/src/reconstruction_octree.rs b/splashsurf_lib/src/reconstruction_octree.rs deleted file mode 100644 index 37422da..0000000 --- a/splashsurf_lib/src/reconstruction_octree.rs +++ /dev/null @@ -1,737 +0,0 @@ -//! Helper functions calling the individual steps of the reconstruction pipeline - -use crate::generic_tree::*; -use crate::marching_cubes::SurfacePatch; -use crate::mesh::TriMesh3d; -use crate::octree::{NodeData, Octree, OctreeNode}; -use crate::uniform_grid::{OwningSubdomainGrid, Subdomain, UniformGrid}; -use crate::workspace::LocalReconstructionWorkspace; -use crate::{ - density_map, marching_cubes, neighborhood_search, new_map, profile, utils, Index, - OctreeDecompositionParameters, Parameters, ParticleDensityComputationStrategy, Real, - ReconstructionError, SpatialDecomposition, SurfaceReconstruction, -}; -use log::{debug, info, trace}; -use nalgebra::Vector3; -use num_traits::Bounded; -use parking_lot::Mutex; - -/// Performs a global surface reconstruction without domain decomposition -pub(crate) fn reconstruct_surface_global<'a, I: Index, R: Real>( - particle_positions: &[Vector3], - parameters: &Parameters, - output_surface: &'a mut SurfaceReconstruction, -) -> Result<(), ReconstructionError> { - profile!("reconstruct_surface_global"); - - // Multiple local workspaces are only needed for processing different subdomains in parallel. - // However, in this global surface reconstruction without domain decomposition, each step in the - // reconstruction pipeline manages its memory on its own. - let mut workspace = output_surface - .workspace - .get_local_with_capacity(particle_positions.len()) - .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() > output_surface.workspace.densities().capacity() { - std::mem::swap( - output_surface.particle_densities.as_mut().unwrap(), - &mut workspace.particle_densities, - ); - } - } - - // Clear the current mesh, as reconstruction will be appended to output - output_surface.mesh.clear(); - // Perform global reconstruction without octree - reconstruct_single_surface_append( - &mut *workspace, - &output_surface.grid, - None, - particle_positions, - None, - parameters, - &mut output_surface.mesh, - )?; - - // TODO: Set this correctly - output_surface.density_map = None; - output_surface.particle_densities = Some(std::mem::take(&mut workspace.particle_densities)); - - Ok(()) -} - -/// Performs a surface reconstruction with an octree for domain decomposition -pub(crate) fn reconstruct_surface_domain_decomposition<'a, I: Index, R: Real>( - particle_positions: &[Vector3], - parameters: &Parameters, - output_surface: &'a mut SurfaceReconstruction, -) -> Result<(), ReconstructionError> { - profile!("reconstruct_surface_domain_decomposition"); - - OctreeBasedSurfaceReconstruction::new(particle_positions, parameters, output_surface) - .expect("Unable to construct octree. Missing/invalid decomposition parameters?") - .run(particle_positions, output_surface)?; - - Ok(()) -} - -/// Helper type for performing an octree based surface reconstruction -struct OctreeBasedSurfaceReconstruction { - /// General parameters for the surface reconstruction - parameters: Parameters, - /// Spatial decomposition specific parameters extracted from the general parameters - spatial_decomposition: OctreeDecompositionParameters, - /// The implicit global grid for the entire surface reconstruction - grid: UniformGrid, - /// Octree containing the individual particle subdomains, built during the construction of this helper type - octree: Octree, -} - -// TODO: Make this less object oriented? -impl OctreeBasedSurfaceReconstruction { - /// Initializes the octree based surface reconstruction by constructing the corresponding octree - fn new( - global_particle_positions: &[Vector3], - parameters: &Parameters, - output_surface: &SurfaceReconstruction, - ) -> Option { - // The grid was already generated by the calling public function - let grid = output_surface.grid.clone(); - - let Some(SpatialDecomposition::Octree(decomposition_parameters)) = - ¶meters.spatial_decomposition - else { - // TODO: Use default values instead? - - // If there are no decomposition parameters, we cannot construct an octree. - return None; - }; - - // Construct the octree - let octree = { - let margin_factor = decomposition_parameters - .ghost_particle_safety_factor - .unwrap_or(R::one()); - - Octree::new_subdivided( - &grid, - global_particle_positions, - decomposition_parameters.subdivision_criterion.clone(), - parameters.compact_support_radius * margin_factor, - parameters.enable_multi_threading, - decomposition_parameters.enable_stitching, - ) - }; - - // Disable all multi-threading in sub-tasks for now (instead, entire sub-tasks are processed in parallel) - let parameters = { - let mut p = parameters.clone(); - p.enable_multi_threading = false; - p - }; - - Some(Self { - octree, - spatial_decomposition: decomposition_parameters.clone(), - grid, - parameters, - }) - } - - /// Runs the surface reconstruction on the initialized octree - fn run( - self, - global_particle_positions: &[Vector3], - output_surface: &mut SurfaceReconstruction, - ) -> Result<(), ReconstructionError> { - // 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() > output_surface.workspace.densities().capacity() { - std::mem::swap( - output_surface.particle_densities.as_mut().unwrap(), - output_surface.workspace.densities_mut(), - ); - } - } - - // Compute particle densities depending on the selected strategy - let global_particle_densities_vec = - match self.spatial_decomposition.particle_density_computation { - // Strategy 1: compute particle densities globally - ParticleDensityComputationStrategy::Global => { - Some(Self::compute_particle_densities_global( - global_particle_positions, - &self.grid, - &self.parameters, - output_surface, - )); - Some(std::mem::take(output_surface.workspace.densities_mut())) - } - // Strategy 2: compute and merge particle densities per subdomain - ParticleDensityComputationStrategy::SynchronizeSubdomains => { - Some(Self::compute_particle_densities_local( - global_particle_positions, - &self.grid, - &self.octree, - &self.parameters, - output_surface, - )); - Some(std::mem::take(output_surface.workspace.densities_mut())) - } - // Strategy 3: each subdomain will compute densities later on its own - // (can only work correctly if margin is large enough) - ParticleDensityComputationStrategy::IndependentSubdomains => None, - }; - - { - let global_particle_densities = - global_particle_densities_vec.as_ref().map(|v| v.as_slice()); - - // Run surface reconstruction - if self.spatial_decomposition.enable_stitching { - self.run_with_stitching( - global_particle_positions, - global_particle_densities, - output_surface, - )?; - } else { - self.run_without_stitching( - global_particle_positions, - global_particle_densities, - output_surface, - )?; - } - - info!( - "Global mesh has {} triangles and {} vertices.", - output_surface.mesh.triangles.len(), - output_surface.mesh.vertices.len() - ); - } - - output_surface.octree = Some(self.octree); - output_surface.density_map = None; - output_surface.particle_densities = global_particle_densities_vec; - - Ok(()) - } - - /// Computes the particle densities globally on all particles without any domain decomposition - fn compute_particle_densities_global( - global_particle_positions: &[Vector3], - grid: &UniformGrid, - parameters: &Parameters, - output_surface: &mut SurfaceReconstruction, - ) { - let mut densities = std::mem::take(output_surface.workspace.densities_mut()); - - { - let mut workspace = output_surface.workspace.get_local().borrow_mut(); - compute_particle_densities_and_neighbors( - grid, - global_particle_positions, - parameters, - &mut workspace.particle_neighbor_lists, - &mut densities, - ); - } - - *output_surface.workspace.densities_mut() = densities; - } - - /// Computes the particles densities per subdomain followed by merging them into a global vector - fn compute_particle_densities_local( - global_particle_positions: &[Vector3], - grid: &UniformGrid, - octree: &Octree, - parameters: &Parameters, - output_surface: &mut SurfaceReconstruction, - ) { - profile!( - parent_scope, - "parallel subdomain particle density computation" - ); - info!("Starting computation of particle densities."); - - // Take the global density storage from workspace to move it behind a mutex - let mut global_densities = std::mem::take(output_surface.workspace.densities_mut()); - utils::resize_and_fill( - &mut global_densities, - global_particle_positions.len(), - ::min_value(), - parameters.enable_multi_threading, - ); - let global_densities = Mutex::new(global_densities); - - let tl_workspaces = &output_surface.workspace; - - octree - .root() - .par_visit_bfs(|octree_node: &OctreeNode| { - profile!( - "visit octree node for density computation", - parent = parent_scope - ); - - let node_particles = if let Some(particle_set) = octree_node.data().particle_set() { - &particle_set.particles - } else { - // Skip non-leaf nodes - return; - }; - - let mut tl_workspace_ref_mut = tl_workspaces - .get_local_with_capacity(node_particles.len()) - .borrow_mut(); - let tl_workspace = &mut *tl_workspace_ref_mut; - - Self::collect_node_particle_positions( - node_particles, - global_particle_positions, - &mut tl_workspace.particle_positions, - ); - - compute_particle_densities_and_neighbors( - grid, - tl_workspace.particle_positions.as_slice(), - parameters, - &mut tl_workspace.particle_neighbor_lists, - &mut tl_workspace.particle_densities, - ); - - { - profile!("update global density values"); - - let mut global_densities = global_densities.lock(); - for (&global_idx, (&density, position)) in node_particles.iter().zip( - tl_workspace - .particle_densities - .iter() - .zip(tl_workspace.particle_positions.iter()), - ) { - // Check if the particle is actually inside of the cell and not a ghost particle - if octree_node.aabb().contains_point(position) { - global_densities[global_idx] = density; - } - } - } - }); - - // Unpack densities from mutex and move back into workspace - *output_surface.workspace.densities_mut() = global_densities.into_inner(); - } - - /// Performs surface reconstruction without stitching by visiting all octree leaf nodes - fn run_without_stitching( - &self, - global_particle_positions: &[Vector3], - global_particle_densities: Option<&[R]>, - output_surface: &mut SurfaceReconstruction, - ) -> Result<(), ReconstructionError> { - // Clear all local meshes - { - let tl_workspaces = &mut output_surface.workspace; - // Clear all thread local meshes - tl_workspaces - .local_workspaces_mut() - .iter_mut() - .for_each(|local_workspace| { - local_workspace.borrow_mut().mesh.clear(); - }); - } - - // Perform individual surface reconstructions on all non-empty leaves of the octree - { - let tl_workspaces = &output_surface.workspace; - - profile!(parent_scope, "parallel subdomain surf. rec."); - info!("Starting triangulation of surface patches."); - - self.octree - .root() - .try_par_visit_bfs(|octree_node: &OctreeNode| -> Result<(), ReconstructionError> { - let particles = if let Some(particle_set) = octree_node.data().particle_set() { - &particle_set.particles - } else { - // Skip non-leaf nodes: as soon as all leaves are processed all work is done - return Ok(()); - }; - - profile!("visit octree node for reconstruction", parent = parent_scope); - trace!("Processing octree leaf with {} particles", particles.len()); - - if particles.is_empty() { - return Ok(()); - } else { - let subdomain_grid = self.extract_node_subdomain(octree_node); - - debug!( - "Surface reconstruction of local patch with {} particles. (offset: {:?}, cells_per_dim: {:?})", - particles.len(), - subdomain_grid.subdomain_offset(), - subdomain_grid.subdomain_grid().cells_per_dim()); - - let mut tl_workspace = tl_workspaces - .get_local_with_capacity(particles.len()) - .borrow_mut(); - - // Take particle position storage from workspace and fill it with positions of the leaf - let mut node_particle_positions = std::mem::take(&mut tl_workspace.particle_positions); - Self::collect_node_particle_positions(particles, global_particle_positions, &mut node_particle_positions); - - // Take particle density storage from workspace and fill it with densities of the leaf - let node_particle_densities = if let Some(global_particle_densities) = global_particle_densities { - let mut node_particle_densities = std::mem::take(&mut tl_workspace.particle_densities); - Self::collect_node_particle_densities(particles, global_particle_densities, &mut node_particle_densities); - Some(node_particle_densities) - } else { - None - }; - - // Take the thread local mesh and append to it without clearing - let mut node_mesh = std::mem::take(&mut tl_workspace.mesh); - - reconstruct_single_surface_append( - &mut *tl_workspace, - &self.grid, - Some(&subdomain_grid), - node_particle_positions.as_slice(), - node_particle_densities.as_ref().map(|v| v.as_slice()), - &self.parameters, - &mut node_mesh, - )?; - - trace!("Surface patch successfully processed."); - - // Put back everything taken from the workspace - tl_workspace.particle_positions = node_particle_positions; - tl_workspace.mesh = node_mesh; - if let Some(node_particle_densities) = node_particle_densities { - tl_workspace.particle_densities = node_particle_densities; - } - - Ok(()) - } - })?; - }; - - // Append all thread local meshes to global mesh - { - let tl_workspaces = &mut output_surface.workspace; - tl_workspaces.local_workspaces_mut().iter_mut().fold( - &mut output_surface.mesh, - |global_mesh, local_workspace| { - global_mesh.append(&mut local_workspace.borrow_mut().mesh); - global_mesh - }, - ); - } - - Ok(()) - } - - /// Performs surface reconstruction with stitching on octree using DFS visitation: reconstruct leaf nodes first, then stitch the parent node as soon as all children are processed - fn run_with_stitching( - &self, - global_particle_positions: &[Vector3], - global_particle_densities: Option<&[R]>, - output_surface: &mut SurfaceReconstruction, - ) -> Result<(), ReconstructionError> { - let mut octree = self.octree.clone(); - - // Perform individual surface reconstructions on all non-empty leaves of the octree - { - let tl_workspaces = &output_surface.workspace; - - profile!( - parent_scope, - "parallel domain decomposed surf. rec. with stitching" - ); - info!("Starting triangulation of surface patches."); - - octree - .root_mut() - // Use DFS visitation as we can only start stitching after all child nodes of one node are reconstructed/stitched. - .try_par_visit_mut_dfs_post(|octree_node: &mut OctreeNode| -> Result<(), ReconstructionError> { - profile!("visit octree node (reconstruct or stitch)", parent = parent_scope); - - // Extract the set of particles of the current node - let particles = if let Some(particle_set) = octree_node.data().particle_set() { - &particle_set.particles - } else { - // If node has no particle set, its children were already processed so it can be stitched - octree_node.stitch_surface_patches(self.parameters.iso_surface_threshold)?; - // After stitching we can directly continue visting the next node - return Ok(()); - }; - - trace!("Processing octree leaf with {} particles", particles.len()); - - let subdomain_grid = self.extract_node_subdomain(octree_node); - let surface_patch = if particles.is_empty() { - SurfacePatch::new_empty(subdomain_grid) - } else { - debug!( - "Reconstructing surface of local patch with {} particles. (offset: {:?}, cells_per_dim: {:?})", - particles.len(), - subdomain_grid.subdomain_offset(), - subdomain_grid.subdomain_grid().cells_per_dim() - ); - - let mut tl_workspace = tl_workspaces - .get_local_with_capacity(particles.len()) - .borrow_mut(); - - // Take particle position storage from workspace and fill it with positions of the leaf - let mut node_particle_positions = std::mem::take(&mut tl_workspace.particle_positions); - Self::collect_node_particle_positions(particles, global_particle_positions, &mut node_particle_positions); - - // Take particle density storage from workspace and fill it with densities of the leaf - let node_particle_densities = if let Some(global_particle_densities) = global_particle_densities { - let mut node_particle_densities = std::mem::take(&mut tl_workspace.particle_densities); - Self::collect_node_particle_densities(particles, global_particle_densities, &mut node_particle_densities); - Some(node_particle_densities) - } else { - None - }; - - let surface_patch = reconstruct_surface_patch( - &mut *tl_workspace, - &subdomain_grid, - node_particle_positions.as_slice(), - node_particle_densities.as_ref().map(|v| v.as_slice()), - &self.parameters, - ); - - // Put back everything taken from the workspace - tl_workspace.particle_positions = node_particle_positions; - if let Some(node_particle_densities) = node_particle_densities { - tl_workspace.particle_densities = node_particle_densities; - } - - surface_patch? - }; - - trace!("Surface patch successfully processed."); - - // Store triangulation in the leaf - octree_node - .data_mut() - .replace(NodeData::SurfacePatch(surface_patch.into())); - - Ok(()) - })?; - - info!("Generation of surface patches is done."); - }; - - // Move stitched mesh out of octree - { - let surface_path = octree - .root_mut() - .data_mut() - .take() - .into_surface_patch() - .expect("Cannot extract stitched mesh from root node") - .patch; - output_surface.mesh = surface_path.mesh; - } - - Ok(()) - } - - /// Computes the subdomain grid for the given octree node - fn extract_node_subdomain(&self, octree_node: &OctreeNode) -> OwningSubdomainGrid { - let grid = &self.grid; - - let subdomain_grid = octree_node - .grid(octree_node.aabb().min(), grid.cell_size()) - .expect("Unable to construct Octree node grid"); - let subdomain_offset = octree_node.min_corner(); - subdomain_grid.log_grid_info(); - - OwningSubdomainGrid::new(grid.clone(), subdomain_grid, *subdomain_offset.index()) - } - - /// Collects the particle positions of all particles in the node - fn collect_node_particle_positions( - node_particles: &[usize], - global_particle_positions: &[Vector3], - node_particle_positions: &mut Vec>, - ) { - node_particle_positions.clear(); - utils::reserve_total(node_particle_positions, node_particles.len()); - - // Extract the particle positions of the leaf - node_particle_positions.extend( - node_particles - .iter() - .map(|&idx| global_particle_positions[idx]), - ); - } - - /// Collects the density values of all particles in the node - fn collect_node_particle_densities( - node_particles: &[usize], - global_particle_densities: &[R], - node_particle_densities: &mut Vec, - ) { - node_particle_densities.clear(); - utils::reserve_total(node_particle_densities, node_particles.len()); - - // Extract the particle densities of the leaf - node_particle_densities.extend( - node_particles - .iter() - .map(|&idx| global_particle_densities[idx]), - ); - } -} - -/// Computes per particle densities into the workspace, also performs the required neighborhood search -pub(crate) fn compute_particle_densities_and_neighbors( - grid: &UniformGrid, - particle_positions: &[Vector3], - parameters: &Parameters, - particle_neighbor_lists: &mut Vec>, - densities: &mut Vec, -) { - profile!("compute_particle_densities_and_neighbors"); - - let particle_rest_density = parameters.rest_density; - let particle_rest_volume = R::from_f64((4.0 / 3.0) * std::f64::consts::PI).unwrap() - * parameters.particle_radius.powi(3); - let particle_rest_mass = particle_rest_volume * particle_rest_density; - - trace!("Starting neighborhood search..."); - neighborhood_search::search_inplace::( - &grid.aabb(), - particle_positions, - parameters.compact_support_radius, - parameters.enable_multi_threading, - particle_neighbor_lists, - ); - - trace!("Computing particle densities..."); - density_map::compute_particle_densities_inplace::( - particle_positions, - particle_neighbor_lists.as_slice(), - parameters.compact_support_radius, - particle_rest_mass, - parameters.enable_multi_threading, - densities, - ); -} - -/// Reconstruct a surface, appends triangulation to the given mesh -pub(crate) fn reconstruct_single_surface_append<'a, I: Index, R: Real>( - workspace: &mut LocalReconstructionWorkspace, - grid: &UniformGrid, - subdomain_grid: Option<&OwningSubdomainGrid>, - particle_positions: &[Vector3], - particle_densities: Option<&[R]>, - parameters: &Parameters, - output_mesh: &'a mut TriMesh3d, -) -> Result<(), ReconstructionError> { - let particle_rest_density = parameters.rest_density; - let particle_rest_volume = R::from_f64((4.0 / 3.0) * std::f64::consts::PI).unwrap() - * parameters.particle_radius.powi(3); - let particle_rest_mass = particle_rest_volume * particle_rest_density; - - let particle_densities = if let Some(particle_densities) = particle_densities { - assert_eq!(particle_densities.len(), particle_positions.len()); - particle_densities - } else { - compute_particle_densities_and_neighbors( - grid, - particle_positions, - parameters, - &mut workspace.particle_neighbor_lists, - &mut workspace.particle_densities, - ); - workspace.particle_densities.as_slice() - }; - - // Create a new density map, reusing memory with the workspace is bad for cache efficiency - // Alternatively one could reuse memory with a custom caching allocator - let mut density_map = new_map().into(); - density_map::generate_sparse_density_map( - grid, - subdomain_grid, - particle_positions, - particle_densities, - None, - particle_rest_mass, - parameters.compact_support_radius, - parameters.cube_size, - parameters.enable_multi_threading, - &mut density_map, - )?; - - marching_cubes::triangulate_density_map_append( - grid, - subdomain_grid, - &density_map, - parameters.iso_surface_threshold, - output_mesh, - )?; - - Ok(()) -} - -/// Reconstruct a surface, appends triangulation to the given mesh -pub(crate) fn reconstruct_surface_patch( - workspace: &mut LocalReconstructionWorkspace, - subdomain_grid: &OwningSubdomainGrid, - particle_positions: &[Vector3], - particle_densities: Option<&[R]>, - parameters: &Parameters, -) -> Result, ReconstructionError> { - profile!("reconstruct_surface_patch"); - - let particle_rest_density = parameters.rest_density; - let particle_rest_volume = R::from_f64((4.0 / 3.0) * std::f64::consts::PI).unwrap() - * parameters.particle_radius.powi(3); - let particle_rest_mass = particle_rest_volume * particle_rest_density; - - let particle_densities = if let Some(particle_densities) = particle_densities { - assert_eq!(particle_densities.len(), particle_positions.len()); - particle_densities - } else { - compute_particle_densities_and_neighbors( - subdomain_grid.global_grid(), - particle_positions, - parameters, - &mut workspace.particle_neighbor_lists, - &mut workspace.particle_densities, - ); - workspace.particle_densities.as_slice() - }; - - // Create a new density map, reusing memory with the workspace is bad for cache efficiency - // Alternatively, one could reuse memory with a custom caching allocator - let mut density_map = new_map().into(); - density_map::generate_sparse_density_map( - subdomain_grid.global_grid(), - Some(subdomain_grid), - particle_positions, - particle_densities, - None, - particle_rest_mass, - parameters.compact_support_radius, - parameters.cube_size, - parameters.enable_multi_threading, - &mut density_map, - )?; - - // Run marching cubes and get boundary data - let patch = marching_cubes::triangulate_density_map_to_surface_patch::( - subdomain_grid, - &density_map, - parameters.iso_surface_threshold, - )?; - - Ok(patch) -} diff --git a/splashsurf_lib/src/uniform_grid.rs b/splashsurf_lib/src/uniform_grid.rs index 998e996..d9e2696 100644 --- a/splashsurf_lib/src/uniform_grid.rs +++ b/splashsurf_lib/src/uniform_grid.rs @@ -82,105 +82,6 @@ pub struct NeighborEdge<'a, 'b: 'a, I: Index> { connectivity: DirectedAxis, } -/// A [`UniformGrid`] that represents a subdomain with an offset inside of a global grid, can be used for mapping between the grids using the [`Subdomain`] trait -#[derive(Clone, Debug)] -pub struct OwningSubdomainGrid { - /// The global or outer grid - global_grid: UniformGrid, - /// The smaller subdomain grid inside of the global grid - subdomain_grid: UniformGrid, - /// The offset of the subdomain grid relative to the global grid - subdomain_offset: [I; 3], -} - -/// A dummy subdomain grid that is actually identical to its global grid, note that this type assumes that any point and cell indices are valid for this grid -#[derive(Clone, Debug)] -pub struct DummySubdomain<'a, I: Index, R: Real> { - /// The global or outer grid - global_grid: &'a UniformGrid, - /// The zero offset of the grid - subdomain_offset: [I; 3], -} - -/// Trait for subdomains of a global grid, providing mappings of indices between the grids -pub trait Subdomain { - /// Returns a reference to the global grid corresponding to the subdomain - fn global_grid(&self) -> &UniformGrid; - /// Returns a reference to the subdomain grid - fn subdomain_grid(&self) -> &UniformGrid; - /// Returns the offset of the subdomain grid relative to the global grid - fn subdomain_offset(&self) -> &[I; 3]; - - /// Returns the lower corner point index of the subdomain in the global grid - fn min_point(&self) -> PointIndex { - self.global_grid() - .get_point(self.subdomain_offset().clone()) - .expect("Invalid subdomain") - } - - /// Returns the upper corner point index of the subdomain in the global grid - fn max_point(&self) -> PointIndex { - let max_point = Direction::Positive - // Take `cells_per_dim` steps as this is equal to `points_per_dim - 1` - .checked_apply_step_ijk( - &self.subdomain_offset(), - self.subdomain_grid().cells_per_dim(), - ) - .expect("Invalid subdomain"); - self.global_grid() - .get_point(max_point) - .expect("Invalid subdomain") - } - - /// Maps the point index from the global grid into the subdomain grid - fn map_point(&self, global_point: &PointIndex) -> Option> { - let new_point = Direction::Negative - .checked_apply_step_ijk(global_point.index(), &self.subdomain_offset())?; - self.subdomain_grid().get_point(new_point) - } - - /// Maps the point index from the subdomain grid to the global grid - fn inv_map_point(&self, subdomain_point: &PointIndex) -> Option> { - let new_point = Direction::Positive - .checked_apply_step_ijk(subdomain_point.index(), &self.subdomain_offset())?; - self.global_grid().get_point(new_point) - } - - /// Maps a point from this subdomain to the other - fn map_point_to>( - &self, - other_subdomain: &S, - subdomain_point: &PointIndex, - ) -> Option> { - let global_point = self.inv_map_point(subdomain_point)?; - other_subdomain.map_point(&global_point) - } - - /// Maps the cell index from the global grid into the subdomain grid - fn map_cell(&self, global_cell: &CellIndex) -> Option> { - let new_cell = Direction::Negative - .checked_apply_step_ijk(global_cell.index(), self.subdomain_offset())?; - self.subdomain_grid().get_cell(new_cell) - } - - /// Maps the cell index from the subdomain grid to the global grid - fn inv_map_cell(&self, subdomain_cell: &CellIndex) -> Option> { - let new_cell = Direction::Positive - .checked_apply_step_ijk(subdomain_cell.index(), self.subdomain_offset())?; - self.global_grid().get_cell(new_cell) - } - - /// Maps a cell from this subdomain to the other - fn map_cell_to>( - &self, - other_subdomain: &S, - subdomain_cell: &CellIndex, - ) -> Option> { - let global_point = self.inv_map_cell(subdomain_cell)?; - other_subdomain.map_cell(&global_point) - } -} - bitflags! { /// Flags naming the outer faces of a grid cell or an entire grid, can be used to select multiple faces at once #[derive(Copy, Clone, Debug)] @@ -789,103 +690,6 @@ impl UniformCartesianCubeGrid3d { } } -impl OwningSubdomainGrid { - /// Creates a new subdomain grid - pub fn new( - global_grid: UniformGrid, - subdomain_grid: UniformGrid, - subdomain_offset: [I; 3], - ) -> Self { - OwningSubdomainGrid { - global_grid, - subdomain_grid, - subdomain_offset, - } - } -} - -impl Subdomain for OwningSubdomainGrid { - /// Returns a reference to the global grid corresponding to the subdomain - #[inline(always)] - fn global_grid(&self) -> &UniformGrid { - &self.global_grid - } - /// Returns a reference to the subdomain grid - #[inline(always)] - fn subdomain_grid(&self) -> &UniformGrid { - &self.subdomain_grid - } - /// Returns the offset of the subdomain grid relative to the global grid - #[inline(always)] - fn subdomain_offset(&self) -> &[I; 3] { - &self.subdomain_offset - } -} - -impl<'a, I: Index, R: Real> DummySubdomain<'a, I, R> { - /// Construct a dummy subdomain without any offset from the given global grid - pub fn new(global_grid: &'a UniformGrid) -> Self { - Self { - global_grid, - subdomain_offset: [I::zero(), I::zero(), I::zero()], - } - } -} - -/// Provide more efficient implementations specific for a dummy grid, note that this type assumes that any point and cell indices provided for mapping are valid for this grid -impl<'a, I: Index, R: Real> Subdomain for DummySubdomain<'a, I, R> { - #[inline(always)] - fn global_grid(&self) -> &UniformGrid { - &self.global_grid - } - - #[inline(always)] - fn subdomain_grid(&self) -> &UniformGrid { - &self.global_grid - } - - #[inline(always)] - fn subdomain_offset(&self) -> &[I; 3] { - &self.subdomain_offset - } - - /// Returns the lower corner of the global domain - #[inline(always)] - fn min_point(&self) -> PointIndex { - PointIndex::from_ijk(self.subdomain_offset) - } - - /// Returns the upper corner of the global domain - #[inline(always)] - fn max_point(&self) -> PointIndex { - PointIndex::from_ijk(self.global_grid.n_cells_per_dim) - } - - /// Just clones the given point index, does not check if it is part of the grid - #[inline(always)] - fn map_point(&self, global_point: &PointIndex) -> Option> { - Some(global_point.clone()) - } - - /// Just clones the given point index, does not check if it is part of the grid - #[inline(always)] - fn inv_map_point(&self, subdomain_point: &PointIndex) -> Option> { - Some(subdomain_point.clone()) - } - - /// Just clones the given cell index, does not check if it is part of the grid - #[inline(always)] - fn map_cell(&self, global_cell: &CellIndex) -> Option> { - Some(global_cell.clone()) - } - - /// Just clones the given cell index, does not check if it is part of the grid - #[inline(always)] - fn inv_map_cell(&self, subdomain_cell: &CellIndex) -> Option> { - Some(subdomain_cell.clone()) - } -} - impl PointIndex { #[inline(always)] fn from_ijk(point_ijk: [I; 3]) -> Self { diff --git a/splashsurf_lib/src/utils.rs b/splashsurf_lib/src/utils.rs index fdcd96a..5675ad5 100644 --- a/splashsurf_lib/src/utils.rs +++ b/splashsurf_lib/src/utils.rs @@ -94,6 +94,7 @@ pub(crate) fn reserve_total(vec: &mut Vec, total_capacity: usize) { } /// Resizes the given vector to the given length and fills new entries with `value.clone()`, parallel or sequential depending on runtime parameter +#[allow(unused)] pub(crate) fn resize_and_fill( vec: &mut Vec, new_len: usize, @@ -108,6 +109,7 @@ pub(crate) fn resize_and_fill( } /// Resizes the given vector to the given length and fills new entries with `value.clone()`, sequential version +#[allow(unused)] pub(crate) fn seq_resize_and_fill(vec: &mut Vec, new_len: usize, value: T) { let old_len = vec.len(); vec.iter_mut() @@ -117,6 +119,7 @@ pub(crate) fn seq_resize_and_fill(vec: &mut Vec, new_len: usize, va } /// Resizes the given vector to the given length and fills new entries with `value.clone()`, parallel version +#[allow(unused)] pub(crate) fn par_resize_and_fill( vec: &mut Vec, new_len: usize, diff --git a/splashsurf_lib/src/workspace.rs b/splashsurf_lib/src/workspace.rs index eac3a2e..e239051 100644 --- a/splashsurf_lib/src/workspace.rs +++ b/splashsurf_lib/src/workspace.rs @@ -1,7 +1,6 @@ //! Workspace for reusing allocated memory between multiple surface reconstructions -use crate::mesh::TriMesh3d; -use crate::{new_map, DensityMap, Index, Real}; +use crate::Real; use nalgebra::Vector3; use std::cell::RefCell; use std::fmt; @@ -10,59 +9,36 @@ use thread_local::ThreadLocal; /// Collection of all thread local workspaces used to reduce allocations on subsequent surface reconstructions #[derive(Default)] -pub struct ReconstructionWorkspace { - global_densities: Vec, +pub struct ReconstructionWorkspace { /// Temporary storage for storing a filtered set of the user provided particles filtered_particles: Vec>, - local_workspaces: ThreadLocal>>, + local_workspaces: ThreadLocal>>, } -impl ReconstructionWorkspace { - /// Returns a reference to the global particle density vector - pub(crate) fn densities(&self) -> &Vec { - &self.global_densities - } - - /// Returns a mutable reference to the global particle density vector - pub(crate) fn densities_mut(&mut self) -> &mut Vec { - &mut self.global_densities - } - +impl ReconstructionWorkspace { /// Returns a mutable reference to the global filtered particles vector pub(crate) fn filtered_particles_mut(&mut self) -> &mut Vec> { &mut self.filtered_particles } - /// Returns a reference to a thread local workspace - pub(crate) fn get_local(&self) -> &RefCell> { - self.local_workspaces.get_or_default() - } - /// Returns a reference to a thread local workspace, initializes it with the given capacity if not already initialized pub(crate) fn get_local_with_capacity( &self, capacity: usize, - ) -> &RefCell> { + ) -> &RefCell> { self.local_workspaces .get_or(|| RefCell::new(LocalReconstructionWorkspace::with_capacity(capacity))) } - - /// Returns a mutable reference to the thread local workspaces - pub(crate) fn local_workspaces_mut( - &mut self, - ) -> &mut ThreadLocal>> { - &mut self.local_workspaces - } } -impl Clone for ReconstructionWorkspace { +impl Clone for ReconstructionWorkspace { /// Returns a new default workspace without any allocated memory fn clone(&self) -> Self { ReconstructionWorkspace::default() } } -impl Debug for ReconstructionWorkspace { +impl Debug for ReconstructionWorkspace { /// Only print the name of type to the formatter fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { f.debug_struct("ReconstructionWorkspace").finish() @@ -70,47 +46,34 @@ impl Debug for ReconstructionWorkspace { } /// Workspace used by [`reconstruct_surface_inplace`] internally to re-use allocated memory -pub(crate) struct LocalReconstructionWorkspace { - /// Storage for the particle positions (only used in octree based approach) - pub particle_positions: Vec>, +pub(crate) struct LocalReconstructionWorkspace { /// Storage for per particle neighbor lists pub particle_neighbor_lists: Vec>, /// Storage for per particle densities pub particle_densities: Vec, - /// Storage for the final surface mesh - pub mesh: TriMesh3d, - /// Storage for the density level-set - #[allow(unused)] - pub density_map: DensityMap, } -impl Default for LocalReconstructionWorkspace { +impl Default for LocalReconstructionWorkspace { /// Constructs a workspace without allocating additional memory fn default() -> Self { Self::new() } } -impl LocalReconstructionWorkspace { +impl LocalReconstructionWorkspace { /// Constructs a workspace without allocating additional memory pub(crate) fn new() -> Self { Self { - particle_positions: Default::default(), particle_neighbor_lists: Default::default(), particle_densities: Default::default(), - mesh: Default::default(), - density_map: new_map().into(), } } /// Constructs a workspace with capacity for the given number of particles pub(crate) fn with_capacity(capacity: usize) -> Self { Self { - particle_positions: Vec::with_capacity(capacity), particle_neighbor_lists: Vec::with_capacity(capacity), particle_densities: Vec::with_capacity(capacity), - mesh: Default::default(), - density_map: new_map().into(), } } } diff --git a/splashsurf_lib/tests/integration_tests/mod.rs b/splashsurf_lib/tests/integration_tests/mod.rs index 9b76c4c..2a21be6 100644 --- a/splashsurf_lib/tests/integration_tests/mod.rs +++ b/splashsurf_lib/tests/integration_tests/mod.rs @@ -3,5 +3,3 @@ pub mod test_full; #[cfg(feature = "io")] pub mod test_mesh; pub mod test_neighborhood_search; -#[cfg(feature = "io")] -pub mod test_octree; diff --git a/splashsurf_lib/tests/integration_tests/test_full.rs b/splashsurf_lib/tests/integration_tests/test_full.rs index c46d9b1..ff4b13e 100644 --- a/splashsurf_lib/tests/integration_tests/test_full.rs +++ b/splashsurf_lib/tests/integration_tests/test_full.rs @@ -3,9 +3,8 @@ use splashsurf_lib::io::particles_from_file; use splashsurf_lib::io::vtk_format::write_vtk; use splashsurf_lib::marching_cubes::check_mesh_consistency; use splashsurf_lib::{ - reconstruct_surface, Aabb3d, GridDecompositionParameters, OctreeDecompositionParameters, - Parameters, ParticleDensityComputationStrategy, Real, SpatialDecomposition, - SubdivisionCriterion, + reconstruct_surface, Aabb3d, GridDecompositionParameters, Parameters, Real, + SpatialDecomposition, }; use std::path::Path; @@ -14,8 +13,6 @@ use std::path::Path; enum Strategy { Global, - Octree, - OctreeStitching, SubdomainGrid, } @@ -44,28 +41,6 @@ fn params_with_aabb( match strategy { Strategy::Global => {} - Strategy::Octree => { - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(R::one()), - enable_stitching: false, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - } - Strategy::OctreeStitching => { - parameters.spatial_decomposition = Some(SpatialDecomposition::Octree( - OctreeDecompositionParameters { - subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, - ghost_particle_safety_factor: Some(R::one() + R::one()), - enable_stitching: true, - particle_density_computation: - ParticleDensityComputationStrategy::SynchronizeSubdomains, - }, - )); - } Strategy::SubdomainGrid => { parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { @@ -109,12 +84,8 @@ fn default_params() -> Parameters { default_params_with(Strategy::Global) } -fn test_for_boundary(params: &Parameters) -> bool { - match ¶ms.spatial_decomposition { - Some(SpatialDecomposition::Octree(p)) => p.enable_stitching, - Some(SpatialDecomposition::UniformGrid(_)) => true, - None => true, - } +fn test_for_boundary(_params: &Parameters) -> bool { + true } macro_rules! generate_test { @@ -170,23 +141,15 @@ macro_rules! generate_test { } generate_test!(f32, surface_reconstruction_bunny_global, "bunny_frame_14_7705_particles.vtk" => "reconstruct_surface_bunny_par_global.vtk", default_params(), 60000, 80000, cfg_attr(debug_assertions, ignore)); -generate_test!(f32, surface_reconstruction_bunny_no_stitching, "bunny_frame_14_7705_particles.vtk" => "reconstruct_surface_bunny_par_no_stitching.vtk", default_params_with(Strategy::Octree), 60000, 80000); -generate_test!(f32, surface_reconstruction_bunny_stitching, "bunny_frame_14_7705_particles.vtk" => "reconstruct_surface_bunny_par_stitching.vtk", default_params_with(Strategy::OctreeStitching), 60000, 80000, cfg_attr(debug_assertions, ignore)); generate_test!(f32, surface_reconstruction_bunny_grid, "bunny_frame_14_7705_particles.vtk" => "reconstruct_surface_bunny_par_grid.vtk", default_params_with(Strategy::SubdomainGrid), 60000, 80000, cfg_attr(debug_assertions, ignore)); -generate_test!(f32, surface_reconstruction_hexecontahedron_stitching, "pentagonal_hexecontahedron_32286_particles.bgeo" => "reconstruct_surface_pentagonal_hexecontahedron_par_stitching.vtk", default_params_with(Strategy::OctreeStitching), 550000, 650000, cfg_attr(debug_assertions, ignore)); generate_test!(f32, surface_reconstruction_hexecontahedron_grid, "pentagonal_hexecontahedron_32286_particles.bgeo" => "reconstruct_surface_pentagonal_hexecontahedron_par_grid.vtk", default_params_with(Strategy::SubdomainGrid), 550000, 650000, cfg_attr(debug_assertions, ignore)); -generate_test!(f32, surface_reconstruction_hilbert_stitching, "hilbert_46843_particles.bgeo" => "reconstruct_surface_hilbert_par_stitching.vtk", default_params_with(Strategy::OctreeStitching), 360000, 400000, cfg_attr(debug_assertions, ignore)); -generate_test!(f32, surface_reconstruction_hilbert2_stitching, "hilbert2_7954_particles.vtk" => "reconstruct_surface_hilbert2_par_stitching.vtk", params(0.025, 4.0, 1.1, 0.6, Strategy::OctreeStitching), 90000, 100000); -generate_test!(f32, surface_reconstruction_octocat_stitching, "octocat_32614_particles.bgeo" => "reconstruct_surface_octocat_par_stitching.vtk", params(0.025, 4.0, 0.75, 0.6, Strategy::OctreeStitching), 140000, 180000, cfg_attr(debug_assertions, ignore)); - generate_test!(f32, surface_reconstruction_hilbert_grid, "hilbert_46843_particles.bgeo" => "reconstruct_surface_hilbert_par_grid.vtk", default_params_with(Strategy::SubdomainGrid), 360000, 400000, cfg_attr(debug_assertions, ignore)); generate_test!(f32, surface_reconstruction_hilbert2_grid, "hilbert2_7954_particles.vtk" => "reconstruct_surface_hilbert2_par_grid.vtk", params(0.025, 4.0, 1.1, 0.6, Strategy::SubdomainGrid), 90000, 100000); generate_test!(f32, surface_reconstruction_octocat_grid, "octocat_32614_particles.bgeo" => "reconstruct_surface_octocat_par_grid.vtk", params(0.025, 4.0, 0.75, 0.6, Strategy::SubdomainGrid), 140000, 180000, cfg_attr(debug_assertions, ignore)); generate_test!(f32, surface_reconstruction_knot_global, "sailors_knot_19539_particles.vtk" => "reconstruct_surface_knot_par_global.vtk", params(0.025, 4.0, 1.1, 0.6, Strategy::Global), 40000, 70000, cfg_attr(debug_assertions, ignore)); -generate_test!(f32, surface_reconstruction_knot_stitching, "sailors_knot_19539_particles.vtk" => "reconstruct_surface_knot_par_stitching.vtk", params(0.025, 4.0, 1.1, 0.6, Strategy::OctreeStitching), 40000, 70000); generate_test!(f32, surface_reconstruction_knot_grid, "sailors_knot_19539_particles.vtk" => "reconstruct_surface_knot_par_grid.vtk", params(0.025, 4.0, 1.1, 0.6, Strategy::SubdomainGrid), 40000, 70000); generate_test!(f32, surface_reconstruction_free_particles_01, "free_particles_1000_particles.vtk" => "reconstruct_surface_free_particles_01_global.vtk", params(0.5, 4.0, 1.5, 0.45, Strategy::Global), 21000, 25000); diff --git a/splashsurf_lib/tests/integration_tests/test_mesh.rs b/splashsurf_lib/tests/integration_tests/test_mesh.rs index fcfba89..6fcfe09 100644 --- a/splashsurf_lib/tests/integration_tests/test_mesh.rs +++ b/splashsurf_lib/tests/integration_tests/test_mesh.rs @@ -21,7 +21,7 @@ fn test_halfedge_ico() -> Result<(), anyhow::Error> { let (tri_mesh, _vertex_map) = he_mesh.into_parts(true); let _tri_vertex_map = tri_mesh.vertex_vertex_connectivity(); - io::obj_format::mesh_to_obj(&MeshWithData::new(tri_mesh), "../icosphere_new.obj")?; + io::obj_format::mesh_to_obj(&MeshWithData::new(tri_mesh), "../out/icosphere_new.obj")?; Ok(()) } @@ -34,7 +34,7 @@ fn test_halfedge_plane() -> Result<(), anyhow::Error> { v.y = -0.1 * (v.x * v.x + v.z * v.z); } - io::obj_format::mesh_to_obj(&MeshWithData::new(mesh.clone()), "../plane_new_0.obj")?; + io::obj_format::mesh_to_obj(&MeshWithData::new(mesh.clone()), "../out/plane_new_0.obj")?; let mut he_mesh = HalfEdgeTriMesh::from(mesh); @@ -61,6 +61,6 @@ fn test_halfedge_plane() -> Result<(), anyhow::Error> { let (tri_mesh, _vertex_map) = he_mesh.into_parts(true); let _tri_vertex_map = tri_mesh.vertex_vertex_connectivity(); - io::obj_format::mesh_to_obj(&MeshWithData::new(tri_mesh), "../plane_new_1.obj")?; + io::obj_format::mesh_to_obj(&MeshWithData::new(tri_mesh), "../out/plane_new_1.obj")?; Ok(()) } diff --git a/splashsurf_lib/tests/integration_tests/test_octree.rs b/splashsurf_lib/tests/integration_tests/test_octree.rs deleted file mode 100644 index 3773964..0000000 --- a/splashsurf_lib/tests/integration_tests/test_octree.rs +++ /dev/null @@ -1,411 +0,0 @@ -use splashsurf_lib::generic_tree::VisitableTree; -use splashsurf_lib::io; -use splashsurf_lib::nalgebra::Vector3; -use splashsurf_lib::octree::Octree; -use splashsurf_lib::{grid_for_reconstruction, Index, Real, SubdivisionCriterion, UniformGrid}; -use std::path::Path; - -/* -#[allow(dead_code)] -fn particles_to_file, R: Real>(particles: Vec>, path: P) { - let points = PointCloud3d { points: particles }; - io::vtk_format::write_vtk( - UnstructuredGridPiece::from(&points), - path.as_ref(), - "particles", - ) - .unwrap(); -} - */ - -#[allow(dead_code)] -fn octree_to_file, I: Index, R: Real>( - octree: &Octree, - grid: &UniformGrid, - path: P, -) { - let mesh = octree.hexmesh(&grid, false); - io::vtk_format::write_vtk(mesh, path.as_ref(), "octree").unwrap(); -} - -#[test] -fn build_octree() { - let particles_per_dim = Vector3::new(20, 10, 6); - let distance = 0.05; - let start = Vector3::new(-0.5, -0.5, -0.5); - - let mut particles = - Vec::with_capacity(particles_per_dim.x * particles_per_dim.y * particles_per_dim.z); - - for i in 0..particles_per_dim.x { - for j in 0..particles_per_dim.y { - if ((i as f64) > 0.25 * particles_per_dim.x as f64) - && ((i as f64) < 0.75 * particles_per_dim.x as f64) - && ((j as f64) > 0.5 * particles_per_dim.y as f64) - { - continue; - } - - for k in 0..particles_per_dim.z { - let particle = Vector3::new( - start.x + (i + 1) as f64 * distance, - start.y + (j + 1) as f64 * distance, - start.z + (k + 1) as f64 * distance, - ); - particles.push(particle); - } - } - } - - let extents = Vector3::new( - (particles_per_dim.x + 2) as f64, - (particles_per_dim.y + 2) as f64, - (particles_per_dim.z + 2) as f64, - ) * distance; - - let lower_corner = start - Vector3::repeat(distance * 0.11); - - let cube_size = distance * 0.3; - let n_cubes = extents * (2.0 / cube_size); - let n_cubes = [ - n_cubes.x.ceil() as i64, - n_cubes.y.ceil() as i64, - n_cubes.z.ceil() as i64, - ]; - - let grid = UniformGrid::new(&lower_corner, &n_cubes, cube_size).unwrap(); - //println!("{:?}", grid); - - let mut octree = Octree::new(&grid, particles.as_slice().len()); - octree.subdivide_recursively_margin( - &grid, - particles.as_slice(), - SubdivisionCriterion::MaxParticleCount(30), - 0.0, - false, - ); - - /* - let root = octree.root(); - - println!("min: {:?}, max: {:?}", root.min_corner(), root.max_corner()); - println!( - "min coord: {:?}, max coord: {:?}", - grid.point_coordinates(root.min_corner()), - grid.point_coordinates(root.max_corner()) - ); - */ - - //particles_to_file(particles, "U:\\particles.vtk"); - //octree_to_file(&octree, &grid, "U:\\octree.vtk"); -} - -/// Loads particles from a VTK file and checks that octree contains all particles and that each particle is in the correct leaf -#[test] -fn build_octree_from_vtk() { - let file = "../data/double_dam_break_frame_26_4732_particles.vtk"; - let particles = io::vtk_format::particles_from_vtk::(file).unwrap(); - //println!("Loaded {} particles from {}", particles.len(), file); - - let grid = grid_for_reconstruction::( - particles.as_slice(), - 0.025, - 4.0 * 0.025, - 0.2, - None, - true, - ) - .unwrap(); - //println!("{:?}", grid); - - let mut octree = Octree::new(&grid, particles.as_slice().len()); - octree.subdivide_recursively_margin( - &grid, - particles.as_slice(), - SubdivisionCriterion::MaxParticleCount(60), - 0.0, - false, - ); - - // Sum the number of particles per leaf - let mut particle_count = 0; - for node in octree.root().dfs_iter() { - if let Some(particle_set) = node.data().particle_set() { - let node_particles = &particle_set.particles; - //println!("Leaf with: {} particles", node_particles.len()); - particle_count += node_particles.len(); - - // Ensure that all particles are within extents of octree cell - let aabb = node.aabb(); - for &idx in node_particles.iter() { - let particle = particles[idx]; - assert!(aabb.contains_point(&particle)); - } - } - } - // Ensure that every particle was exactly in one cell - assert_eq!(particle_count, particles.len()); - - //octree_to_file(&octree, &grid, "U:\\double_dam_break_frame_26_4732_particles_octree.vtk"); -} - -struct TestParameters { - particle_radius: R, - compact_support_radius: R, - cube_size: R, - margin: Option, - max_particles_per_cell: Option, - compare_seq_par: bool, -} - -impl Default for TestParameters { - fn default() -> Self { - let particle_radius = R::from_f64(0.025).unwrap(); - let compact_support_radius = particle_radius.times_f64(4.0); - let cube_size = particle_radius.times_f64(0.5); - - Self { - particle_radius, - compact_support_radius, - cube_size, - margin: None, - max_particles_per_cell: None, - compare_seq_par: true, - } - } -} - -#[allow(unused)] -impl TestParameters { - fn new(particle_radius: f64, compact_support_factor: f64, cube_size_factor: f64) -> Self { - let params = Self::default(); - params.with_parameters(particle_radius, compact_support_factor, cube_size_factor) - } - - fn with_margin(mut self, margin: Option) -> Self { - self.margin = margin.map(|m| self.particle_radius.times_f64(m)); - self - } - fn with_max_particles_per_cell(mut self, max_particles_per_cell: Option) -> Self { - self.max_particles_per_cell = max_particles_per_cell; - self - } - - fn with_compare_seq_par(mut self, compare_seq_par: bool) -> Self { - self.compare_seq_par = compare_seq_par; - self - } - - fn with_parameters( - mut self, - particle_radius: f64, - compact_support_factor: f64, - cube_size_factor: f64, - ) -> Self { - self.particle_radius = R::from_f64(particle_radius).unwrap(); - self.compact_support_radius = self.particle_radius.times_f64(compact_support_factor); - self.cube_size = self.particle_radius.times_f64(cube_size_factor); - self - } - - fn build_grid(&self, particle_positions: &[Vector3]) -> UniformGrid { - grid_for_reconstruction( - particle_positions, - self.particle_radius, - self.compact_support_radius, - self.cube_size, - None, - true, - ) - .unwrap() - } -} - -/// Returns a vector containing per particle how often it is a non-ghost particle in the octree -fn count_non_ghost_particles( - particle_positions: &[Vector3], - octree: &Octree, -) -> Vec { - let mut non_ghost_particles = vec![0; particle_positions.len()]; - for node in octree.root().dfs_iter() { - if let Some(particle_set) = node.data().particle_set() { - for &idx in particle_set.particles.iter() { - if node.aabb().contains_point(&particle_positions[idx]) { - non_ghost_particles[idx] += 1; - } - } - } - } - - non_ghost_particles -} - -/// Asserts whether each particle has a unique octree node where it is not a ghost particle -fn assert_unique_node_per_particle( - particle_positions: &[Vector3], - octree: &Octree, -) { - let non_ghost_particles_counts = count_non_ghost_particles(particle_positions, octree); - let no_unique_node_particles: Vec<_> = non_ghost_particles_counts - .into_iter() - .enumerate() - .filter(|&(_, count)| count != 1) - .collect(); - - assert_eq!( - no_unique_node_particles, - vec![], - "There are nodes that don't have a unique octree node where they are not ghost particles!" - ); -} - -/// Asserts whether both trees have equivalent particle sets in the same octree nodes -fn assert_tree_equivalence(left_tree: &Octree, right_tree: &Octree) { - for (node_seq, node_par) in left_tree - .root() - .dfs_iter() - .zip(right_tree.root().dfs_iter()) - { - match ( - node_seq.data().particle_set(), - node_par.data().particle_set(), - ) { - (Some(particles_seq), Some(particles_par)) => { - let particles_seq = &particles_seq.particles; - let particles_par = &particles_par.particles; - - // Ensure that we have the same number of particles for sequential and parallel result - assert_eq!(particles_seq.len(), particles_par.len(), "Found corresponding leaves in the trees that don't have the same number of particles!"); - - // Sort the particle lists - let mut particles_seq = particles_seq.to_vec(); - let mut particles_par = particles_par.to_vec(); - particles_seq.sort_unstable(); - particles_par.sort_unstable(); - - // Ensure that the particle lists are identical - assert_eq!(particles_seq, particles_par, "Found corresponding leaves in the trees that don't have the same sorted particle lists!"); - } - (None, None) => { - // Both nodes are leaves - continue; - } - _ => { - // Parallel and sequential nodes do not match (one is leaf, one has children) - panic!("Encountered a node where one octree has a particle set but the other does not!"); - } - } - } -} - -fn build_octree_par_consistency>( - file: P, - parameters: TestParameters, -) { - let particles = io::particles_from_file::(file).unwrap(); - - let grid = parameters.build_grid::(particles.as_slice()); - - let octree_seq = if parameters.compare_seq_par { - let mut octree_seq = Octree::new(&grid, particles.as_slice().len()); - octree_seq.subdivide_recursively_margin( - &grid, - particles.as_slice(), - parameters - .max_particles_per_cell - .map(|n| SubdivisionCriterion::MaxParticleCount(n)) - .unwrap_or(SubdivisionCriterion::MaxParticleCountAuto), - parameters.margin.unwrap_or(R::zero()), - false, - ); - - assert_unique_node_per_particle(particles.as_slice(), &octree_seq); - - Some(octree_seq) - } else { - None - }; - - let mut octree_par = Octree::new(&grid, particles.as_slice().len()); - octree_par.par_subdivide_recursively_margin( - &grid, - particles.as_slice(), - parameters - .max_particles_per_cell - .map(|n| SubdivisionCriterion::MaxParticleCount(n)) - .unwrap_or(SubdivisionCriterion::MaxParticleCountAuto), - parameters.margin.unwrap_or(R::zero()), - false, - ); - - assert_unique_node_per_particle(particles.as_slice(), &octree_par); - - if let Some(octree_seq) = octree_seq { - assert_tree_equivalence(&octree_seq, &octree_par) - } -} - -#[test] -fn build_octree_cube() { - build_octree_par_consistency::( - "../data/cube_2366_particles.vtk", - TestParameters::default() - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(70)), - ); -} - -#[test] -fn build_octree_double_dam_break() { - build_octree_par_consistency::( - "../data/double_dam_break_frame_26_4732_particles.vtk", - TestParameters::default() - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(200)), - ); -} - -#[test] -#[cfg_attr(debug_assertions, ignore)] -fn build_octree_dam_break() { - build_octree_par_consistency::( - "../data/dam_break_frame_23_24389_particles.bgeo", - TestParameters::default() - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(1000)), - ); -} - -#[test] -fn build_octree_bunny() { - build_octree_par_consistency::( - "../data/bunny_frame_14_7705_particles.vtk", - TestParameters::default() - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(200)), - ); -} - -#[test] -#[cfg_attr(debug_assertions, ignore)] -fn build_octree_hilbert() { - build_octree_par_consistency::( - "../data/hilbert_46843_particles.bgeo", - TestParameters::default() - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(1000)), - ); -} - -/* -#[test] -fn build_octree_canyon() { - build_octree_par_consistency::( - "../../canyon_13353401_particles.vtk", - TestParameters::new(0.011, 4.0, 1.5) - .with_margin(Some(1.0)) - .with_max_particles_per_cell(Some(52161)) - .with_compare_seq_par(false), - ); -} -*/