diff --git a/Cargo.lock b/Cargo.lock index 6f90181..653089d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -173,7 +173,7 @@ checksum = "fdde5c9cd29ebd706ce1b35600920a33550e402fc998a2e53ad3b42c3c47a192" dependencies = [ "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", ] [[package]] @@ -275,9 +275,9 @@ dependencies = [ [[package]] name = "clap" -version = "4.3.9" +version = "4.3.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bba77a07e4489fb41bd90e8d4201c3eb246b3c2c9ea2ba0bddd6c1d1df87db7d" +checksum = "384e169cc618c613d5e3ca6404dda77a8685a63e08660dcc64abaf7da7cb0c7a" dependencies = [ "clap_builder", "clap_derive", @@ -286,13 +286,12 @@ dependencies = [ [[package]] name = "clap_builder" -version = "4.3.9" +version = "4.3.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2c9b4a88bb4bc35d3d6f65a21b0f0bafe9c894fa00978de242c555ec28bea1c0" +checksum = "ef137bbe35aab78bdb468ccfba75a5f4d8321ae011d34063770780545176af2d" dependencies = [ "anstream", "anstyle", - "bitflags 1.3.2", "clap_lex", "strsim", ] @@ -306,7 +305,7 @@ dependencies = [ "heck", "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", ] [[package]] @@ -585,9 +584,9 @@ checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" [[package]] name = "hermit-abi" -version = "0.3.1" +version = "0.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fed44880c466736ef9a5c5b5facefb5ed0785676d0c02d612db14e54f0d84286" +checksum = "443144c8cdadd93ebf52ddb4056d257f5b52c04d3c804e657d19eb73fc33668b" [[package]] name = "iana-time-zone" @@ -647,13 +646,12 @@ dependencies = [ [[package]] name = "is-terminal" -version = "0.4.7" +version = "0.4.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "adcf93614601c8129ddf72e2d5633df827ba6551541c6d8c59520a371475be1f" +checksum = "24fddda5af7e54bf7da53067d6e802dbcc381d0a8eef629df528e3ebf68755cb" dependencies = [ "hermit-abi", - "io-lifetimes", - "rustix", + "rustix 0.38.2", "windows-sys 0.48.0", ] @@ -677,9 +675,9 @@ dependencies = [ [[package]] name = "itoa" -version = "1.0.6" +version = "1.0.8" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "453ad9f582a441959e5f0d088b02ce04cfe8d51a8eaf077f12ac6d3e94164ca6" +checksum = "62b02a5381cc465bd3041d84623d0fa3b66738b52b8e2fc3bab8ad63ab032f4a" [[package]] name = "js-sys" @@ -729,6 +727,12 @@ version = "0.3.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ef53942eb7bf7ff43a617b3e2c1c4a5ecf5944a7c1bc12d7ee39bbb15e5c1519" +[[package]] +name = "linux-raw-sys" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09fc20d2ca12cb9f044c93e3bd6d32d523e6e2ec3db4f7b2939cd99026ecd3f0" + [[package]] name = "lock_api" version = "0.4.10" @@ -974,9 +978,9 @@ dependencies = [ [[package]] name = "paste" -version = "1.0.12" +version = "1.0.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9f746c4065a8fa3fe23974dd82f15431cc8d40779821001404d10d2e79ca7d79" +checksum = "b4b27ab7be369122c218afc2079489cdcb4b517c0a3fc386ff11e1fedfcc2b35" [[package]] name = "peg" @@ -1181,9 +1185,21 @@ dependencies = [ [[package]] name = "regex" -version = "1.8.4" +version = "1.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d0ab3ca65655bb1e41f2a8c8cd662eb4fb035e67c3f78da1d61dffe89d07300f" +checksum = "89089e897c013b3deb627116ae56a6955a72b8bed395c9526af31c9fe528b484" +dependencies = [ + "aho-corasick", + "memchr 2.5.0", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fa250384981ea14565685dea16a9ccc4d1c541a13f82b9c168572264d1df8c56" dependencies = [ "aho-corasick", "memchr 2.5.0", @@ -1192,9 +1208,9 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.7.2" +version = "0.7.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "436b050e76ed2903236f032a59761c1eb99e1b0aead2c257922771dab1fc8c78" +checksum = "2ab07dc67230e4a4718e70fd5c20055a4334b121f1f9db8fe63ef39ce9b8c846" [[package]] name = "rstar" @@ -1218,23 +1234,36 @@ dependencies = [ [[package]] name = "rustix" -version = "0.37.20" +version = "0.37.22" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b96e891d04aa506a6d1f318d2771bcb1c7dfda84e126660ace067c9b474bb2c0" +checksum = "8818fa822adcc98b18fedbb3632a6a33213c070556b5aa7c4c8cc21cff565c4c" dependencies = [ "bitflags 1.3.2", "errno", "io-lifetimes", "libc", - "linux-raw-sys", + "linux-raw-sys 0.3.8", + "windows-sys 0.48.0", +] + +[[package]] +name = "rustix" +version = "0.38.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aabcb0461ebd01d6b79945797c27f8529082226cb630a9865a71870ff63532a4" +dependencies = [ + "bitflags 2.3.3", + "errno", + "libc", + "linux-raw-sys 0.4.3", "windows-sys 0.48.0", ] [[package]] name = "ryu" -version = "1.0.13" +version = "1.0.14" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f91339c0467de62360649f8d3e185ca8de4224ff281f66000de5eb2a77a79041" +checksum = "fe232bdf6be8c8de797b22184ee71118d63780ea42ac85b61d1baa6d3b782ae9" [[package]] name = "safe_arch" @@ -1279,29 +1308,29 @@ dependencies = [ [[package]] name = "serde" -version = "1.0.164" +version = "1.0.166" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9e8c8cf938e98f769bc164923b06dce91cea1751522f46f8466461af04c9027d" +checksum = "d01b7404f9d441d3ad40e6a636a7782c377d2abdbe4fa2440e2edcc2f4f10db8" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.164" +version = "1.0.166" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d9735b638ccc51c28bf6914d90a2e9725b377144fc612c49a611fddd1b631d68" +checksum = "5dd83d6dde2b6b2d466e14d9d1acce8816dedee94f735eac6395808b3483c6d6" dependencies = [ "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", ] [[package]] name = "serde_json" -version = "1.0.99" +version = "1.0.100" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "46266871c240a00b8f503b877622fe33430b3c7d963bdc0f2adc511e54a1eae3" +checksum = "0f1e14e89be7aa4c4b78bdbdc9eb5bf8517829a600ae8eaa39a6e1d960b5185c" dependencies = [ "itoa", "ryu", @@ -1353,7 +1382,7 @@ dependencies = [ [[package]] name = "splashsurf" -version = "0.9.3" +version = "0.10.0" dependencies = [ "anyhow", "bytemuck", @@ -1373,7 +1402,7 @@ dependencies = [ [[package]] name = "splashsurf_lib" -version = "0.9.2" +version = "0.10.0" dependencies = [ "anyhow", "arrayvec", @@ -1430,9 +1459,9 @@ dependencies = [ [[package]] name = "syn" -version = "2.0.22" +version = "2.0.23" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2efbeae7acf4eabd6bcdcbd11c92f45231ddda7539edc7806bd1a04a03b24616" +checksum = "59fb7d6d8281a51045d62b8eb3a7d1ce347b76f312af50cd3dc0af39c87c1737" dependencies = [ "proc-macro2", "quote", @@ -1449,28 +1478,28 @@ dependencies = [ "cfg-if", "fastrand", "redox_syscall", - "rustix", + "rustix 0.37.22", "windows-sys 0.48.0", ] [[package]] name = "thiserror" -version = "1.0.40" +version = "1.0.41" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "978c9a314bd8dc99be594bc3c175faaa9794be04a5a5e153caba6915336cebac" +checksum = "c16a64ba9387ef3fdae4f9c1a7f07a0997fce91985c0336f1ddc1822b3b37802" dependencies = [ "thiserror-impl", ] [[package]] name = "thiserror-impl" -version = "1.0.40" +version = "1.0.41" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f9456a42c5b0d803c8cd86e73dd7cc9edd429499f37a3550d286d5e86720569f" +checksum = "d14928354b01c4d6a4f0e549069adef399a284e7995c7ccca94e8a07a5346c59" dependencies = [ "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", ] [[package]] @@ -1530,9 +1559,9 @@ dependencies = [ [[package]] name = "unicode-ident" -version = "1.0.9" +version = "1.0.10" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b15811caf2415fb889178633e7724bad2509101cde276048e013b9def5e51fa0" +checksum = "22049a19f4a68748a168c0fc439f9516686aa045927ff767eca0a85101fb6e73" [[package]] name = "unicode-width" @@ -1614,7 +1643,7 @@ dependencies = [ "once_cell", "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", "wasm-bindgen-shared", ] @@ -1636,7 +1665,7 @@ checksum = "54681b18a46765f095758388f2d0cf16eb8d4169b639ab575a8f5693af210c7b" dependencies = [ "proc-macro2", "quote", - "syn 2.0.22", + "syn 2.0.23", "wasm-bindgen-backend", "wasm-bindgen-shared", ] diff --git a/splashsurf/Cargo.toml b/splashsurf/Cargo.toml index ad240a2..033b064 100644 --- a/splashsurf/Cargo.toml +++ b/splashsurf/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "splashsurf" -version = "0.9.3" +version = "0.10.0" authors = ["Fabian Löschner "] license = "MIT" description = "Command-line tool for surface reconstruction of SPH particle data" @@ -13,7 +13,7 @@ homepage = "http://splashsurf.physics-simulation.org" repository = "https://github.com/InteractiveComputerGraphics/splashsurf" [dependencies] -splashsurf_lib = { path = "../splashsurf_lib", version = "0.9.2", features = ["vtk_extras", "profiling", "io"] } +splashsurf_lib = { path = "../splashsurf_lib", version = "0.10.0", features = ["vtk_extras", "profiling", "io"] } clap = { version = "4.3", features = ["derive"] } log = "0.4" fern = "0.6" diff --git a/splashsurf/src/reconstruction.rs b/splashsurf/src/reconstruction.rs index 9de8811..60f0fe2 100644 --- a/splashsurf/src/reconstruction.rs +++ b/splashsurf/src/reconstruction.rs @@ -118,6 +118,20 @@ pub struct ReconstructSubcommandArgs { #[arg(help_heading = ARGS_ADV, long, short = 'n')] pub num_threads: Option, + /// Whether to enable spatial decomposition using a regular grid-based approach + #[arg( + help_heading = ARGS_OCTREE, + long, + default_value = "off", + value_name = "off|on", + ignore_case = true, + require_equals = true + )] + pub subdomain_grid: Switch, + /// Each subdomain will be a cube consisting of this number of MC cube cells along each coordinate axis + #[arg(help_heading = ARGS_OCTREE, long, default_value="64")] + pub subdomain_cubes: u32, + /// Whether to enable spatial decomposition using an octree (faster) instead of a global approach #[arg( help_heading = ARGS_OCTREE, @@ -229,6 +243,8 @@ impl Switch { /// Executes the `reconstruct` subcommand pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<(), anyhow::Error> { + profile!("reconstruct CLI"); + let paths = ReconstructionRunnerPathCollection::try_from(cmd_args) .context("Failed parsing input file path(s) from command line")? .collect(); @@ -388,6 +404,10 @@ mod arguments { cube_size, iso_surface_threshold: args.surface_threshold, domain_aabb, + subdomain_num_cubes_per_dim: args + .subdomain_grid + .into_bool() + .then_some(args.subdomain_cubes), enable_multi_threading: args.parallelize_over_particles.into_bool(), spatial_decomposition, }; diff --git a/splashsurf_lib/Cargo.toml b/splashsurf_lib/Cargo.toml index 2bb4ea6..1ab3f6f 100644 --- a/splashsurf_lib/Cargo.toml +++ b/splashsurf_lib/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "splashsurf_lib" -version = "0.9.2" +version = "0.10.0" authors = ["Fabian Löschner "] license = "MIT" description = "Library for surface reconstruction of SPH particle data" @@ -67,7 +67,7 @@ serde_json = { version = "1.0", optional = true } lazy_static = { version = "1.4", optional = true } [dev-dependencies] -criterion = "0.5" +criterion = "0.5.1" ultraviolet = "0.9" sdfu = { git = "https://github.com/w1th0utnam3/sdfu", features = ["ultraviolet"], rev = "e39a4a8685a56a3430218b9f2dfd546ab2dbe2d6" } diff --git a/splashsurf_lib/benches/benches/bench_full.rs b/splashsurf_lib/benches/benches/bench_full.rs index f558a17..4d5efbb 100644 --- a/splashsurf_lib/benches/benches/bench_full.rs +++ b/splashsurf_lib/benches/benches/bench_full.rs @@ -22,7 +22,7 @@ pub fn surface_reconstruction_canyon(c: &mut Criterion) { let parameters = Parameters { particle_radius, rest_density: 1000.0, - compact_support_radius: compact_support_radius, + compact_support_radius, cube_size, iso_surface_threshold: 0.6, domain_aabb: None, @@ -98,11 +98,12 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) { let parameters = Parameters { particle_radius, rest_density: 1000.0, - compact_support_radius: compact_support_radius, + compact_support_radius, cube_size, iso_surface_threshold: 0.6, domain_aabb: None, enable_multi_threading: true, + subdomain_num_cubes_per_dim: None, spatial_decomposition: None, }; @@ -156,6 +157,15 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) { }, ); + group.bench_function("surface_reconstruction_dam_break_par_grid_64", |b| { + b.iter(|| { + let mut parameters = parameters.clone(); + parameters.subdomain_num_cubes_per_dim = Some(64); + reconstruction = + reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap() + }) + }); + group.finish(); /* @@ -179,11 +189,12 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) { let parameters = Parameters { particle_radius, rest_density: 1000.0, - compact_support_radius: compact_support_radius, + compact_support_radius, cube_size, iso_surface_threshold: 0.6, domain_aabb: None, enable_multi_threading: true, + subdomain_num_cubes_per_dim: None, spatial_decomposition: None, }; @@ -237,6 +248,15 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) { }, ); + group.bench_function("surface_reconstruction_double_dam_break_par_grid_64", |b| { + b.iter(|| { + let mut parameters = parameters.clone(); + parameters.subdomain_num_cubes_per_dim = Some(64); + reconstruction = + reconstruct_surface::(particle_positions.as_slice(), ¶meters).unwrap() + }) + }); + group.finish(); /* @@ -260,11 +280,12 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) { let parameters = Parameters { particle_radius, rest_density: 1000.0, - compact_support_radius: compact_support_radius, + compact_support_radius, cube_size, iso_surface_threshold: 0.6, domain_aabb: None, enable_multi_threading: true, + subdomain_num_cubes_per_dim: None, spatial_decomposition: None, }; @@ -335,6 +356,22 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) { }, ); + group.bench_function( + "surface_reconstruction_double_dam_break_inplace_par_grid_64", + |b| { + b.iter(|| { + let mut parameters = parameters.clone(); + parameters.subdomain_num_cubes_per_dim = Some(64); + reconstruct_surface_inplace::( + particle_positions.as_slice(), + ¶meters, + &mut reconstruction, + ) + .unwrap() + }) + }, + ); + group.finish(); write_vtk( diff --git a/splashsurf_lib/benches/benches/bench_mesh.rs b/splashsurf_lib/benches/benches/bench_mesh.rs index feb507e..e04abac 100644 --- a/splashsurf_lib/benches/benches/bench_mesh.rs +++ b/splashsurf_lib/benches/benches/bench_mesh.rs @@ -18,11 +18,12 @@ fn reconstruct_particles>(particle_file: P) -> SurfaceReconstruct let parameters = Parameters { particle_radius, rest_density: 1000.0, - compact_support_radius: compact_support_radius, + compact_support_radius, cube_size, iso_surface_threshold: 0.6, domain_aabb: None, enable_multi_threading: true, + subdomain_num_cubes_per_dim: None, spatial_decomposition: Some(SpatialDecompositionParameters { subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto, ghost_particle_safety_factor: None, diff --git a/splashsurf_lib/src/dense_subdomains.rs b/splashsurf_lib/src/dense_subdomains.rs new file mode 100644 index 0000000..fa69131 --- /dev/null +++ b/splashsurf_lib/src/dense_subdomains.rs @@ -0,0 +1,1543 @@ +use anyhow::Context; +use arrayvec::ArrayVec; +use log::{info, trace}; +use nalgebra::Vector3; +use num_integer::Integer; +use num_traits::{FromPrimitive, NumCast}; +use parking_lot::Mutex; +use rayon::prelude::*; +use std::cell::RefCell; +use std::sync::atomic::{AtomicUsize, Ordering}; +use thread_local::ThreadLocal; + +use crate::density_map::sequential_compute_particle_densities; +use crate::kernel::{CubicSplineKernel, SymmetricKernel3d}; +use crate::marching_cubes::marching_cubes_lut::marching_cubes_triangulation_iter; +use crate::mesh::{HexMesh3d, TriMesh3d}; +use crate::neighborhood_search::{ + neighborhood_search_spatial_hashing, neighborhood_search_spatial_hashing_parallel, +}; +use crate::uniform_grid::{EdgeIndex, PointIndex, UniformCartesianCubeGrid3d}; +use crate::{ + new_map, new_parallel_map, profile, Aabb3d, MapType, Parameters, SurfaceReconstruction, +}; +use crate::{Index, Real}; + +type GlobalIndex = u64; + +/// Converts any literal or expression to the Index type I (panics if value does not fit) +macro_rules! to_index { + ($value:literal) => { + ::from($value).expect("literal has to fit in index type") + }; + ($value:expr) => { + ::from($value).expect("value has to fit in index type") + }; +} + +/// Converts any literal or expression to the Real type R (panics if value does not fit) +macro_rules! to_real { + ($value:literal) => { + ::from($value).expect("literal has to fit in real type") + }; + ($value:expr) => { + ::from($value).expect("value has to fit in real type") + }; +} + +pub(crate) struct ParametersSubdomainGrid { + /// SPH particle radius (in simulation units) + #[allow(unused)] + particle_radius: R, + /// Rest mass of each particle + particle_rest_mass: R, + /// SPH kernel compact support radius (in simulation units) + compact_support_radius: R, + /// Density value for the iso-surface + surface_threshold: R, + /// MC cube size (in simulation units) + cube_size: R, + /// Size of a subdomain in multiplies of MC cubes + subdomain_cubes: I, + /// Margin for ghost particles around each subdomain + ghost_particle_margin: R, + /// Implicit global MC background grid (required to compute consistent float coordinates at domain boundaries) + global_marching_cubes_grid: UniformCartesianCubeGrid3d, + /// Implicit subdomain grid + subdomain_grid: UniformCartesianCubeGrid3d, + /// Chunk size for chunked parallel processing + chunk_size: usize, +} + +/// Result of the subdomain decomposition procedure +pub(crate) struct Subdomains { + // Flat subdomain coordinate indices (same order as the particle list) + flat_subdomain_indices: Vec, + // Particles of each subdomain (including ghost particles) + per_subdomain_particles: Vec>, +} + +pub(crate) fn initialize_parameters<'a, I: Index, R: Real>( + parameters: &Parameters, + _particles: &[Vector3], + output_surface: &'a SurfaceReconstruction, +) -> Result, anyhow::Error> { + let chunk_size = 500; + + // A subdomain will be a cube consisting of this number of MC cubes along each coordinate axis + let subdomain_cubes_in = parameters.subdomain_num_cubes_per_dim.unwrap_or(64); + let subdomain_cubes = I::from_u32(subdomain_cubes_in) + .expect("number of subdomain cubes has to fit in index type"); + let subdomain_cubes_global = GlobalIndex::from_u32(subdomain_cubes_in) + .expect("number of subdomain cubes has to fit in global index type"); + + // Physical particle properties + let particle_radius = parameters.particle_radius; + let particle_rest_density = parameters.rest_density; + let compact_support_radius = parameters.compact_support_radius; + let cube_size = parameters.cube_size; + let surface_threshold = parameters.iso_surface_threshold; + + let particle_rest_volume = to_real!(4) * R::frac_pi_3() * particle_radius.powi(3); + let particle_rest_mass = particle_rest_volume * particle_rest_density; + + let ghost_particle_margin = + (compact_support_radius / cube_size).ceil() * cube_size * to_real!(1.01); + + // Compute information of ghost margin volume for debugging + { + let ghost_margin_cubes = I::from((ghost_particle_margin / cube_size).ceil()) + .expect("ghost margin cube count has to fit in index type"); + + let vol_subdomain = subdomain_cubes + .checked_cubed() + .expect("number of cubes per subdomain has to be representable in index type"); + let vol_margin = (ghost_margin_cubes * to_index!(2) + subdomain_cubes) + .checked_cubed() + .expect( + "number of cubes per subdomain with margin has to be representable in index type", + ) + - vol_subdomain; + + info!( + "The ghost margin volume is {:.2}% of the subdomain volume", + (to_real!(vol_margin) / to_real!(vol_subdomain)) * to_real!(100.0) + ); + info!( + "The ghost margin is {:.2} MC cells or {:.2} subdomains thick", + ghost_particle_margin / cube_size, + ghost_particle_margin / (cube_size * subdomain_cubes.to_real_unchecked()) + ); + + if ghost_margin_cubes > subdomain_cubes / to_index!(2) { + panic!("The ghost margin is {ghost_margin_cubes} cubes thick (rounded up), while the subdomain only has an extent of {subdomain_cubes} cubes. The subdomain has to have at least twice the number of cubes ({})!", ghost_margin_cubes.times(2)); + } + } + + // AABB of the particles + let aabb = output_surface.grid.aabb(); + + let global_mc_grid = UniformCartesianCubeGrid3d::::new( + aabb.min(), + &output_surface + .grid + .cells_per_dim() + .map(|c| ::from(c).unwrap()), + cube_size, + ) + .context("construct initial global marching cubes cell grid")?; + trace!("Initial global MC Grid: {:?}", global_mc_grid); + + // MC cubes along each coordinate axis of the entire global MC background grid + let cells_per_dim = global_mc_grid.cells_per_dim(); + // Compute the number of subdomains along each coordinate axis + let num_subdomains = [ + int_ceil_div(cells_per_dim[0], subdomain_cubes_global), + int_ceil_div(cells_per_dim[1], subdomain_cubes_global), + int_ceil_div(cells_per_dim[2], subdomain_cubes_global), + ]; + + let num_global_mc_cells = (|| -> Option<_> { + Some([ + num_subdomains[0].checked_mul(subdomain_cubes_global)?, + num_subdomains[1].checked_mul(subdomain_cubes_global)?, + num_subdomains[2].checked_mul(subdomain_cubes_global)?, + ]) + })() + .context("compute global number of marching cubes cells per dimension")?; + + let global_mc_grid = UniformCartesianCubeGrid3d::::new( + &global_mc_grid.aabb().min(), + &num_global_mc_cells, + cube_size, + ) + .context("construct final global marching cubes cell grid")?; + trace!("Global MC Grid: {:?}", global_mc_grid); + + // Convert number of subdomains back to local index type + let num_subdomains = (|| -> Option<_> { + Some([ + I::from(num_subdomains[0])?, + I::from(num_subdomains[1])?, + I::from(num_subdomains[2])?, + ]) + })() + .context("convert number of subdomains per dimension to local index type")?; + + // Compute total number of subdomains + let subdomain_count: I = (|| -> Option<_> { + num_subdomains[0].checked_mul(&num_subdomains[1].checked_mul(&num_subdomains[2])?) + })() + .context("compute total number of subdomains")?; + // Edge length of a subdomain in absolute units + let subdomain_size = cube_size * subdomain_cubes.to_real_unchecked(); + // Background grid of the subdomains + let subdomain_grid = UniformCartesianCubeGrid3d::::new( + &global_mc_grid.aabb().min(), + &num_subdomains, + subdomain_size, + )?; + + { + let nc = subdomain_cubes; + let [nx, ny, nz] = subdomain_grid.cells_per_dim(); + let [mc_x, mc_y, mc_z] = global_mc_grid.cells_per_dim(); + info!("Number of subdomains: {subdomain_count} ({nx}x{ny}x{nz})"); + info!( + "Number of MC cells per subdomain: {} ({nc}x{nc}x{nc})", + nc.cubed() + ); + info!( + "Number of MC cells globally: {} ({mc_x}x{mc_y}x{mc_z})", + mc_x * mc_y * mc_z + ); + trace!("Subdomain grid: {:?}", subdomain_grid); + } + + Ok(ParametersSubdomainGrid { + particle_radius, + particle_rest_mass, + compact_support_radius, + surface_threshold, + cube_size, + subdomain_cubes, + ghost_particle_margin, + global_marching_cubes_grid: global_mc_grid, + subdomain_grid, + chunk_size, + }) +} + +#[allow(unused)] +pub(crate) fn extract_narrow_band( + parameters: &ParametersSubdomainGrid, + particles: &[Vector3], +) -> Vec> { + profile!("Filter narrow band"); + + let compact_support_radius = parameters.compact_support_radius; + let ghost_particle_margin = (compact_support_radius / parameters.cube_size).ceil() + * parameters.cube_size + * to_real!(1.01); + + // AABB of the particles + let aabb = { + let mut aabb = Aabb3d::::par_from_points(&particles); + // Add some safety margin, this should be large enough such that all mesh vertices are guaranteed to be inside of it + aabb.grow_uniformly(ghost_particle_margin); + info!("Enlarged AABB: {:?}", aabb); + aabb + }; + + let neighbor_lists = { + profile!("Global neighborhood search"); + let mut neighbors = Vec::new(); + neighborhood_search_spatial_hashing_parallel::( + &aabb, + &particles, + compact_support_radius, + &mut neighbors, + ); + neighbors + }; + + let narrow_band = { + profile!("Identify narrow band"); + let surface_particle_neighbors = 30; + let surface_particles = neighbor_lists + .iter() + .enumerate() + .filter_map(|(i, nl)| (nl.len() < surface_particle_neighbors).then_some(i)) + .collect::>(); + + info!(""); + info!( + "Number of pure \"surface particles\": {}", + surface_particles.len() + ); + + let mut first_ring = surface_particles + .iter() + .copied() + .map(|i| neighbor_lists[i].iter().copied()) + .flatten() + .collect::>(); + info!("First ring before dedup: {}", first_ring.len()); + { + profile!("Dedup first ring"); + first_ring.sort_unstable(); + first_ring.dedup(); + } + info!("First ring after dedup: {}", first_ring.len()); + + let mut second_ring = first_ring + .iter() + .copied() + .map(|i| neighbor_lists[i].iter().copied()) + .flatten() + .collect::>(); + info!("Second ring before dedup: {}", second_ring.len()); + { + profile!("Dedup second ring"); + second_ring.sort_unstable(); + second_ring.dedup(); + } + info!("Second ring after dedup: {}", second_ring.len()); + + let mut narrow_band = surface_particles; + narrow_band.append(&mut first_ring); + narrow_band.append(&mut second_ring); + info!("All before dedup: {}", narrow_band.len()); + { + profile!("Dedup entire narrow band"); + narrow_band.sort_unstable(); + narrow_band.dedup(); + } + info!( + "All after dedup: {} ({:.3}%), interior particles: {}", + narrow_band.len(), + (narrow_band.len() as f64 / particles.len() as f64) * 100.0, + particles.len() - narrow_band.len() + ); + narrow_band + }; + + { + profile!("Collect narrow band positions"); + let narrow_band_positions = narrow_band + .into_iter() + .map(|i| particles[i].clone()) + .collect::>(); + narrow_band_positions + } +} + +/// Performs classification and decomposition of particles into a regular grid of subdomains +pub(crate) fn decomposition< + I: Index, + R: Real, + C: subdomain_classification::ParticleToSubdomainClassifier, +>( + parameters: &ParametersSubdomainGrid, + particles: &[Vector3], +) -> Result, anyhow::Error> { + profile!("decomposition"); + info!("Starting classification of particles into subdomains."); + + // Count the number of particles and ghost particles per subdomain (with thread local counters) + let per_subdomain_counter_tls = ThreadLocal::>>::new(); + { + profile!("classifying particles"); + + particles + .par_chunks(parameters.chunk_size) + .for_each(|particle_chunk| { + let mut per_subdomain_counter = per_subdomain_counter_tls + .get_or(|| RefCell::new(new_map())) + .borrow_mut(); + let mut classifier = C::new(); + + for particle in particle_chunk.iter() { + classifier.classify_particle( + particle, + ¶meters.subdomain_grid, + parameters.ghost_particle_margin, + ); + for i in 0..classifier.len() { + let flat_subdomain_idx = classifier.get(i); + *per_subdomain_counter.entry(flat_subdomain_idx).or_insert(0) += 1; + } + } + }); + } + + // Merge all thread local subdomain particle counters + let global_per_subdomain_counter = new_parallel_map(); + { + profile!("merging TL per cell particle counters"); + + let per_subdomain_counter_tls = per_subdomain_counter_tls + .into_iter() + .map(RefCell::into_inner) + .collect::>(); + + per_subdomain_counter_tls + .into_par_iter() + .for_each(|per_cell_counter| { + for (flat_cell_index, count) in per_cell_counter { + *global_per_subdomain_counter + .entry(flat_cell_index) + .or_insert(0) += count; + } + }); + } + + // Mapping from flat subdomain coordinate index to offset into contiguous subdomain storage + let mut subdomain_compressed_indices = new_map(); + // Inverse mapping (offset to flat subdomain coordinate) + let mut flat_subdomain_indices = vec![I::zero(); global_per_subdomain_counter.len()]; + + let per_subdomain_particle_count: Vec = { + profile!("initializing flat subdomain data and index mapping"); + + global_per_subdomain_counter + .into_iter() + .enumerate() + .map(|(i, (flat_cell_index, particle_count))| { + subdomain_compressed_indices.insert(flat_cell_index, i); + flat_subdomain_indices[i] = flat_cell_index; + particle_count + }) + .collect() + }; + + let mut per_subdomain_particles = Vec::with_capacity(per_subdomain_particle_count.len()); + per_subdomain_particles.resize_with(per_subdomain_particle_count.len(), || { + Mutex::new(Vec::new()) + }); + + { + profile!("copying particles to subdomains"); + + particles + .par_chunks(parameters.chunk_size) + .enumerate() + .for_each(|(chunk_idx, particle_chunk)| { + let chunk_offset = chunk_idx * parameters.chunk_size; + let mut classifier = C::new(); + for (particle_idx, particle) in particle_chunk.iter().enumerate() { + let particle_idx = chunk_offset + particle_idx; + classifier.classify_particle( + particle, + ¶meters.subdomain_grid, + parameters.ghost_particle_margin, + ); + for i in 0..classifier.len() { + let flat_subdomain_idx = classifier.get(i); + + let compressed_subdomain_idx = + subdomain_compressed_indices[&flat_subdomain_idx]; + + // Lock the subdomain for writing + let mut subdomain_particles = + per_subdomain_particles[compressed_subdomain_idx].lock(); + // Reserve full size of subdomain if it's still empty + if subdomain_particles.is_empty() { + let particle_count = + per_subdomain_particle_count[compressed_subdomain_idx]; + subdomain_particles.reserve(particle_count); + } + // Add the particle to the subdomain + subdomain_particles.push(particle_idx); + } + } + }); + } + + // Remove mutexes + let mut per_subdomain_particles = per_subdomain_particles + .into_iter() + .map(Mutex::into_inner) + .collect::>(); + + // Sort subdomain particle sets by index so that overlapping particle regions of subdomains + // will be in the same order + { + profile!("sort subdomain particles"); + per_subdomain_particles + .par_iter_mut() + .for_each(|particles| { + //use rand::prelude::SliceRandom; + //let mut rng = rand::thread_rng(); + //particles.shuffle(&mut rng) + particles.sort_unstable(); + }); + } + + Ok(Subdomains { + flat_subdomain_indices, + per_subdomain_particles, + }) +} + +pub(crate) fn compute_global_density_vector( + parameters: &ParametersSubdomainGrid, + global_particles: &[Vector3], + subdomains: &Subdomains, +) -> Vec { + profile!(parent, "compute_global_density_vector"); + info!("Starting computation of global density vector."); + + let global_particle_densities = Mutex::new(vec![R::zero(); global_particles.len()]); + + #[derive(Default)] + struct SubdomainWorkspace { + // Particle positions of this subdomain + subdomain_particles: Vec>, + // Per particle neighborhood lists + neighborhood_lists: Vec>, + // Per particle density values of this subdomain + particle_densities: Vec, + } + + let workspace_tls = ThreadLocal::>>::new(); + + subdomains + .flat_subdomain_indices + .par_iter() + .copied() + .zip(subdomains.per_subdomain_particles.par_iter()) + .for_each(|(flat_subdomain_idx, subdomain_particle_indices)| { + profile!("inner subdomain_tasks", parent = parent); + + // Obtain thread local workspace and clear it + let mut workspace = workspace_tls.get_or_default().borrow_mut(); + + let SubdomainWorkspace { + subdomain_particles, + neighborhood_lists, + particle_densities, + } = &mut *workspace; + + let flat_subdomain_idx: I = flat_subdomain_idx; + let subdomain_particle_indices: &[usize] = subdomain_particle_indices.as_slice(); + + // Collect all particle positions of this subdomain + { + profile!("collect subdomain data"); + gather_subdomain_data( + global_particles, + subdomain_particle_indices, + subdomain_particles, + ); + } + + // Get the cell index and AABB of the subdomain + let subdomain_idx = parameters + .subdomain_grid + .try_unflatten_cell_index(flat_subdomain_idx) + .expect("Subdomain cell does not exist"); + let subdomain_aabb = parameters.subdomain_grid.cell_aabb(&subdomain_idx); + + let margin_aabb = { + let mut margin_aabb = subdomain_aabb.clone(); + // TODO: Verify if we can omit this extra margin? + margin_aabb.grow_uniformly(parameters.ghost_particle_margin * to_real!(1.5)); + margin_aabb + }; + + neighborhood_search_spatial_hashing::( + &margin_aabb, + &subdomain_particles, + parameters.compact_support_radius, + neighborhood_lists, + ); + + sequential_compute_particle_densities::( + &subdomain_particles, + neighborhood_lists, + parameters.compact_support_radius, + parameters.particle_rest_mass, + particle_densities, + ); + + // Write particle densities into global storage + { + profile!("update global density values"); + // Lock global vector while this subdomain writes into it + let mut global_particle_densities = global_particle_densities.lock(); + subdomain_particle_indices + .iter() + .copied() + .zip(particle_densities.iter().copied()) + // Update density values only for particles inside of the subdomain (ghost particles have wrong values) + .filter(|(particle_idx, _)| { + subdomain_aabb.contains_point(&global_particles[*particle_idx]) + }) + .for_each(|(particle_idx, density)| { + global_particle_densities[particle_idx] = density; + }); + } + }); + + let global_particle_densities = global_particle_densities.into_inner(); + + /* + { + let points = PointCloud3d::new(global_particles.clone()); + let mesh = MeshWithData::new(points).with_cell_data(MeshAttribute::new_real_scalar( + "density".to_string(), + global_particle_densities.clone(), + )); + splashsurf_lib::io::vtk_format::write_vtk(&mesh, "out/new_density.vtk", "particles").unwrap(); + } + */ + + global_particle_densities +} + +pub(crate) struct SurfacePatch { + pub vertices: Vec>, + pub triangles: Vec<[usize; 3]>, + pub vertex_inside_count: usize, + pub triangle_inside_count: usize, + pub vertex_inside_flags: Vec, + pub triangle_inside_flags: Vec, + pub exterior_vertex_edge_indices: Vec<(I, EdgeIndex)>, +} + +pub(crate) fn reconstruction( + parameters: &ParametersSubdomainGrid, + global_particles: &[Vector3], + global_particle_densities: &[R], + subdomains: &Subdomains, +) -> Vec> { + profile!(parent, "reconstruction"); + info!("Starting reconstruction (level-set evaluation and local triangulation)."); + + let squared_support = parameters.compact_support_radius * parameters.compact_support_radius; + // Add 1% so that we don't exclude grid points that are just on the kernel boundary + let squared_support_with_margin = squared_support * to_real!(1.01); + // Compute radial distance in terms of grid points we have to evaluate for each particle + let cube_radius = I::from((parameters.compact_support_radius / parameters.cube_size).ceil()) + .expect("kernel radius in cubes has to fit in index type"); + // Kernel + let kernel = CubicSplineKernel::new(parameters.compact_support_radius); + //let kernel = DiscreteSquaredDistanceCubicKernel::new::(1000, parameters.compact_support_radius); + + let mc_total_cells = parameters.subdomain_cubes.cubed(); + let mc_total_points = (parameters.subdomain_cubes + I::one()).cubed(); + + assert!( + mc_total_points.to_usize().is_some(), + "number of mc cubes per subdomain must be fit into usize" + ); + + let globalize_local_edge = |mc_grid: &UniformCartesianCubeGrid3d, + subdomain_grid: &UniformCartesianCubeGrid3d, + subdomain_index: I, + local_edge: &EdgeIndex| + -> (I, EdgeIndex) { + // We globalize the boundary edge index by translating the local edge index to the subdomain + // where it lies on the lower boundary of that domain. + + let max_mc_point_index = mc_grid.points_per_dim().map(|i| i - I::one()); + let max_subdomain_index = subdomain_grid + .cells_per_dim() + .map(|i| i.saturating_sub(&I::one()).max(I::zero())); + + // Check along which axes this edge is on the max boundary + let is_max = local_edge.axis().orthogonal_axes().map(|orth_axis| { + if local_edge.origin().index()[orth_axis.dim()] == max_mc_point_index[orth_axis.dim()] { + // We are on the max side of this domain along the axis + true + } else { + // We are either + // - On the min side of this domain along the axis + // - Somewhere in the middle (in this case this axis is irrelevant) + false + } + }); + + if !is_max[0] && !is_max[1] { + // Edge is already in the correct subdomain + (subdomain_index, local_edge.clone()) + } else { + // We have to translate to the neighboring subdomain (+1 in all directions where is_max == true) + let subdomain_cell = subdomain_grid + .try_unflatten_cell_index(subdomain_index) + .expect("invalid subdomain index"); + + let mut target_subdomain_ijk = subdomain_cell.index().clone(); + let mut target_local_origin_ijk = local_edge.origin().index().clone(); + + // Obtain index of new subdomain and new origin point + for (&orth_axis, &is_max) in local_edge + .axis() + .orthogonal_axes() + .iter() + .zip(is_max.iter()) + { + if is_max { + // Clamp the step to the subdomain grid because we are not interested in subdomains outside the grid + // (globalization is not needed on the outermost boundary of the entire problem domain) + target_subdomain_ijk[orth_axis.dim()] = (target_subdomain_ijk[orth_axis.dim()] + + I::one()) + .min(max_subdomain_index[orth_axis.dim()]); + // Move origin point from max boundary to min boundary + target_local_origin_ijk[orth_axis.dim()] = I::zero(); + } + } + + let target_subdomain = subdomain_grid + .get_cell(target_subdomain_ijk) + .expect("target subdomain has to exist"); + let flat_target_subdomain = subdomain_grid.flatten_cell_index(&target_subdomain); + + // We re-use the same marching cubes domain here because the domain is anyway rectangular, + // therefore this shift gives the same result + let new_local_edge = mc_grid + .get_edge(target_local_origin_ijk, local_edge.axis()) + .expect("failed to translate edge"); + + (flat_target_subdomain, new_local_edge) + } + }; + + #[derive(Default)] + struct SubdomainWorkspace { + // Particle positions of this subdomain + subdomain_particles: Vec>, + // Per particle density values of this subdomain + subdomain_particle_densities: Vec, + // Cache for the level-set values + levelset_grid: Vec, + // Cache for indices + index_cache: Vec>, + } + + let workspace_tls = ThreadLocal::>>::new(); + + let reconstruct_dense = |flat_subdomain_idx: I, subdomain_particle_indices: &Vec| { + profile!("inner subdomain_tasks (dense)", parent = parent); + + // Obtain thread local workspace and clear it + let mut workspace = workspace_tls.get_or_default().borrow_mut(); + + let SubdomainWorkspace { + subdomain_particles, + subdomain_particle_densities, + levelset_grid, + index_cache: _index_cache, + } = &mut *workspace; + + let flat_subdomain_idx: I = flat_subdomain_idx; + let subdomain_particle_indices: &[usize] = subdomain_particle_indices.as_slice(); + + // Collect all particle positions and densities of this subdomain + { + //profile!("collect subdomain data"); + gather_subdomain_data( + global_particles, + subdomain_particle_indices, + subdomain_particles, + ); + gather_subdomain_data( + global_particle_densities, + subdomain_particle_indices, + subdomain_particle_densities, + ); + } + + // Get the cell index and AABB of the subdomain + let subdomain_idx = parameters + .subdomain_grid + .try_unflatten_cell_index(flat_subdomain_idx) + .expect("Subdomain cell does not exist"); + let subdomain_aabb = parameters.subdomain_grid.cell_aabb(&subdomain_idx); + + let mc_grid = UniformCartesianCubeGrid3d::new( + subdomain_aabb.min(), + &[parameters.subdomain_cubes; 3], + parameters.cube_size, + ) + .unwrap(); + + levelset_grid.fill(R::zero()); + levelset_grid.resize(mc_total_points.to_usize().unwrap(), R::zero()); + + { + //profile!("density grid loop"); + + let extents = mc_grid.points_per_dim(); + + for (p_i, rho_i) in subdomain_particles + .iter() + .copied() + .zip(subdomain_particle_densities.iter().copied()) + { + // Get grid cell containing particle + let particle_cell = mc_grid.enclosing_cell(&p_i); + + // Compute lower and upper bounds of the grid points possibly affected by the particle + // We want to loop over the vertices of the enclosing cells plus all points in `cube_radius` distance from the cell + + let lower = [ + (particle_cell[0] - cube_radius).max(I::zero()), + (particle_cell[1] - cube_radius).max(I::zero()), + (particle_cell[2] - cube_radius).max(I::zero()), + ]; + + let upper = [ + // We add 2 because + // - we want to loop over all grid points of the cell (+1 for upper points) + the radius + // - the upper range limit is exclusive (+1) + (particle_cell[0] + cube_radius + I::two()).min(extents[0]), + (particle_cell[1] + cube_radius + I::two()).min(extents[1]), + (particle_cell[2] + cube_radius + I::two()).min(extents[2]), + ]; + + // Loop over all grid points around the enclosing cell + for i in I::range(lower[0], upper[0]).iter() { + for j in I::range(lower[1], upper[1]).iter() { + for k in I::range(lower[2], upper[2]).iter() { + let point_ijk = [i, j, k]; + let local_point = mc_grid + .get_point(point_ijk) + .expect("point has to be part of the subdomain grid"); + //let point_coordinates = mc_grid.point_coordinates(&point); + + let subdomain_ijk = subdomain_idx.index(); + let mc_cells_per_subdomain = mc_grid.cells_per_dim(); + + fn local_to_global_point_ijk( + local_point_ijk: [I; 3], + subdomain_ijk: [I; 3], + cells_per_subdomain: [I; 3], + ) -> [GlobalIndex; 3] { + let local_point_ijk = local_point_ijk + .map(|i| ::from(i).unwrap()); + let subdomain_ijk = subdomain_ijk + .map(|i| ::from(i).unwrap()); + let cells_per_subdomain = cells_per_subdomain + .map(|i| ::from(i).unwrap()); + let [i, j, k] = local_point_ijk; + + [ + subdomain_ijk[0] * cells_per_subdomain[0] + i, + subdomain_ijk[1] * cells_per_subdomain[1] + j, + subdomain_ijk[2] * cells_per_subdomain[2] + k, + ] + } + + // Use global coordinate calculation for consistency with neighboring domains + let global_point_ijk = local_to_global_point_ijk( + point_ijk, + subdomain_ijk.clone(), + mc_cells_per_subdomain.clone(), + ); + let global_point = parameters + .global_marching_cubes_grid + .get_point(global_point_ijk) + .expect("point has to be part of the global mc grid"); + let point_coordinates = parameters + .global_marching_cubes_grid + .point_coordinates(&global_point); + + let dx = p_i - point_coordinates; + let dx_norm_sq = dx.norm_squared(); + + if dx_norm_sq < squared_support_with_margin { + let v_i = parameters.particle_rest_mass / rho_i; + let r = dx_norm_sq.sqrt(); + let w_ij = kernel.evaluate(r); + //let w_ij = kernel.evaluate(dx_norm_sq); + + let interpolated_value = v_i * w_ij; + + let flat_point_idx = mc_grid.flatten_point_index(&local_point); + let flat_point_idx = flat_point_idx.to_usize().unwrap(); + levelset_grid[flat_point_idx] += interpolated_value; + } + } + } + } + } + } + + let mut vertices = Vec::new(); + let mut triangles = Vec::new(); + + let mut vertex_inside_count = 0; + let mut triangle_inside_count = 0; + + let mut vertex_inside_flags = Vec::new(); + let mut triangle_inside_flags = Vec::new(); + + let mut exterior_vertex_edge_indices = Vec::new(); + + let mut edge_to_vertex = new_map(); + + { + //profile!("mc triangulation loop"); + + for flat_cell_idx in I::range(I::zero(), mc_total_cells).iter() { + let cell = mc_grid.try_unflatten_cell_index(flat_cell_idx).unwrap(); + + let mut vertices_inside = [true; 8]; + for local_point_index in 0..8 { + let point = cell.global_point_index_of(local_point_index).unwrap(); + let flat_point_idx = mc_grid.flatten_point_index(&point); + let flat_point_idx = flat_point_idx.to_usize().unwrap(); + // Get value of density map + let density_value = levelset_grid[flat_point_idx]; + // Update inside/outside surface flag + vertices_inside[local_point_index] = + density_value > parameters.surface_threshold; + } + + for triangle in marching_cubes_triangulation_iter(&vertices_inside) { + let mut global_triangle = [0; 3]; + for (v_idx, local_edge_index) in triangle.iter().copied().enumerate() { + let edge = cell + .global_edge_index_of(local_edge_index as usize) + .unwrap(); + let vertex_index = *edge_to_vertex.entry(edge).or_insert_with(|| { + // TODO: Nonlinear interpolation + + let origin_coords = mc_grid.point_coordinates(&edge.origin()); + let target_coords = mc_grid.point_coordinates(&edge.target()); + + let flat_origin_idx = mc_grid + .flatten_point_index(&edge.origin()) + .to_usize() + .unwrap(); + let flat_target_idx = mc_grid + .flatten_point_index(&edge.target()) + .to_usize() + .unwrap(); + + let origin_value = levelset_grid[flat_origin_idx]; + let target_value = levelset_grid[flat_target_idx]; + + let alpha = (parameters.surface_threshold - origin_value) + / (target_value - origin_value); + let interpolated_coords = + origin_coords * (R::one() - alpha) + target_coords * alpha; + let vertex_coords = interpolated_coords; + + vertices.push(vertex_coords); + let vertex_index = vertices.len() - 1; + + let is_interior_vertex = !mc_grid.is_boundary_edge(&edge); + vertex_inside_count += is_interior_vertex as usize; + vertex_inside_flags.push(is_interior_vertex); + + if !is_interior_vertex { + exterior_vertex_edge_indices.push(globalize_local_edge( + &mc_grid, + ¶meters.subdomain_grid, + flat_subdomain_idx, + &edge, + )); + } + + vertex_index + }); + + global_triangle[v_idx] = vertex_index; + } + + let all_tri_vertices_inside = global_triangle + .iter() + .copied() + .all(|v_idx| vertex_inside_flags[v_idx]); + + triangles.push(global_triangle); + triangle_inside_count += all_tri_vertices_inside as usize; + triangle_inside_flags.push(all_tri_vertices_inside); + } + } + } + + SurfacePatch { + vertices, + triangles, + vertex_inside_count, + triangle_inside_count, + vertex_inside_flags, + triangle_inside_flags, + exterior_vertex_edge_indices, + } + }; + + let mut surface_patches = Vec::with_capacity(subdomains.flat_subdomain_indices.len()); + subdomains + .flat_subdomain_indices + .par_iter() + .copied() + .zip(subdomains.per_subdomain_particles.par_iter()) + .map(|(flat_subdomain_idx, subdomain_particle_indices)| { + reconstruct_dense(flat_subdomain_idx, subdomain_particle_indices) + }) + .collect_into_vec(&mut surface_patches); + + surface_patches +} + +pub(crate) fn stitching( + surface_patches: Vec>, +) -> TriMesh3d { + profile!("stitching"); + info!("Starting stitching of subdomains to global mesh."); + + // Calculate offsets of interior vertices and triangles + let vert_and_tri_offsets = { + profile!("surface patch offset scan"); + + std::iter::once((0, 0)) + .chain(surface_patches.iter().scan((0, 0), |offsets, patch| { + let (vert_offset, tri_offset) = offsets; + *vert_offset += patch.vertex_inside_count; + *tri_offset += patch.triangle_inside_count; + Some(*offsets) + })) + .collect::>() + }; + + let (total_interior_vert_count, total_interior_tri_count) = vert_and_tri_offsets + .last() + .copied() + .expect("there has to be at least one entry in the offset list"); + + let mut interior_vertices = vec![Vector3::::zeros(); total_interior_vert_count]; + let mut interior_triangles = vec![[0, 0, 0]; total_interior_tri_count]; + + let mut exterior_vertices = Vec::new(); + let mut exterior_triangles = Vec::new(); + let mut exterior_vertex_mapping = new_map(); + + { + profile!("copy interior verts/tris and deduplicate exterior verts"); + + if vert_and_tri_offsets.len() > 1 { + vert_and_tri_offsets + .windows(2) + .zip(surface_patches.iter()) + .for_each(|(offsets, patch)| { + if let [start_offsets, end_offsets] = offsets { + let (start_verts, start_tris) = *start_offsets; + let (end_verts, end_tris) = *end_offsets; + + let global_vertex_offset = start_verts; + + let mut local_to_global_vertex_mapping = vec![0; patch.vertices.len()]; + + // Copy interior vertices + { + let out_verts = &mut interior_vertices[start_verts..end_verts]; + + // Copy all interior vertices into global storage + patch + .vertices + .iter() + .zip(patch.vertex_inside_flags.iter()) + .enumerate() + // Skip all exterior vertices + .filter_map(|(i, (v, is_interior))| is_interior.then_some((i, v))) + .enumerate() + .for_each(|(new_local_idx, (old_local_idx, vert))| { + out_verts[new_local_idx] = *vert; + local_to_global_vertex_mapping[old_local_idx] = + global_vertex_offset + new_local_idx; + }); + } + + // Copy interior triangles + { + let out_tris = &mut interior_triangles[start_tris..end_tris]; + + // Copy all interior triangle into global storage + patch + .triangles + .iter() + .zip(patch.triangle_inside_flags.iter()) + // Skip all exterior triangles + .filter_map(|(tri, is_interior)| is_interior.then_some(tri)) + // Count only the interior triangles + .enumerate() + .for_each(|(tri_idx, tri)| { + out_tris[tri_idx] = [ + local_to_global_vertex_mapping[tri[0]], + local_to_global_vertex_mapping[tri[1]], + local_to_global_vertex_mapping[tri[2]], + ]; + }); + } + + // Insert & deduplicate exterior vertices + { + patch + .vertices + .iter() + .zip(patch.vertex_inside_flags.iter()) + .enumerate() + // Skip all interior vertices + .filter_map(|(i, (v, is_interior))| { + (!is_interior).then_some((i, v)) + }) + // For each exterior vertex there is a corresponding globalized edge index + .zip(patch.exterior_vertex_edge_indices.iter()) + .enumerate() + .for_each( + |(_exterior_local_idx, ((old_local_idx, vert), edge_index))| { + let global_index = *exterior_vertex_mapping + .entry(edge_index.clone()) + .or_insert_with(|| { + // Exterior vertices will come after all interior vertices in the mesh + let global_index = total_interior_vert_count + + exterior_vertices.len(); + exterior_vertices.push(*vert); + global_index + }); + local_to_global_vertex_mapping[old_local_idx] = + global_index; + }, + ); + } + + // Insert exterior triangles + { + patch + .triangles + .iter() + .zip(patch.triangle_inside_flags.iter()) + // Skip all exterior triangles + .filter_map(|(tri, is_interior)| (!is_interior).then_some(tri)) + .for_each(|tri| { + exterior_triangles.push(tri.map(|local_vert| { + local_to_global_vertex_mapping[local_vert] + })) + }); + } + } + }); + } + } + + let mut vertices = interior_vertices; + vertices.append(&mut exterior_vertices); + + let mut triangles = interior_triangles; + triangles.append(&mut exterior_triangles); + + TriMesh3d { + vertices, + triangles, + } +} + +pub(crate) mod subdomain_classification { + use super::*; + + /// Trait for assignment of particles to their subdomains + pub trait ParticleToSubdomainClassifier { + /// Constructs a new classifier of this type + fn new() -> Self; + /// Classifies a particle into all subdomains it belongs to + fn classify_particle( + &mut self, + particle: &Vector3, + subdomain_grid: &UniformCartesianCubeGrid3d, + ghost_particle_margin: R, + ); + /// Returns the number of subdomains that were assigned to a particle in the last call to `classify_particle` + fn len(&self) -> usize; + /// Returns the `i`-th subdomain assigned to a particle in the last call to `classify_particle` + fn get(&self, i: usize) -> I; + } + + /// Classifier that assigns only the owning subdomain to a particle + pub struct NoMarginClassifier { + subdomain: I, + } + + impl ParticleToSubdomainClassifier for NoMarginClassifier { + #[inline(always)] + fn new() -> Self { + NoMarginClassifier { + subdomain: I::zero(), + } + } + + #[inline(always)] + fn classify_particle( + &mut self, + particle: &Vector3, + subdomain_grid: &UniformCartesianCubeGrid3d, + _ghost_particle_margin: R, + ) { + // Find the owning subdomain of the particle + let subdomain_ijk = subdomain_grid.enclosing_cell(particle); + // And store its flattened index + let flat_subdomain_idx = subdomain_grid.flatten_cell_index_array(&subdomain_ijk); + self.subdomain = flat_subdomain_idx; + } + + #[inline(always)] + fn len(&self) -> usize { + 1 + } + + #[inline(always)] + fn get(&self, _i: usize) -> I { + self.subdomain + } + } + + /// Classifier that assign a particle to its owning subdomain and all subdomains where it's a ghost particle + pub struct GhostMarginClassifier { + subdomains: ArrayVec, + } + + impl ParticleToSubdomainClassifier for GhostMarginClassifier { + fn new() -> Self { + GhostMarginClassifier { + subdomains: ArrayVec::new(), + } + } + + fn classify_particle( + &mut self, + particle: &Vector3, + subdomain_grid: &UniformCartesianCubeGrid3d, + ghost_particle_margin: R, + ) { + self.subdomains.clear(); + ghost_particle_classification( + particle, + subdomain_grid, + ghost_particle_margin, + &mut self.subdomains, + ) + } + + #[inline(always)] + fn len(&self) -> usize { + self.subdomains.len() + } + + #[inline(always)] + fn get(&self, i: usize) -> I { + self.subdomains[i] + } + } + + /// Assigns a particle into the subdomain that contains it and all subdomains where it's a ghost particle. + fn ghost_particle_classification( + particle: &Vector3, + subdomain_grid: &UniformCartesianCubeGrid3d, + ghost_particle_margin: R, + subdomains: &mut ArrayVec, + ) { + // Find the owning subdomain of the particle + let subdomain_ijk = subdomain_grid.enclosing_cell(particle); + // Make sure particle is part of computational domain + if subdomain_grid.get_cell(subdomain_ijk).is_none() { + return; + } + + // Get corner points spanning the owning subdomain + let subdomain_aabb = subdomain_grid.cell_aabb( + &subdomain_grid + .get_cell(subdomain_ijk) + .expect("Subdomain has to be part of grid"), + ); + let min_corner = subdomain_aabb.min(); + let max_corner = subdomain_aabb.max(); + + // Checks whether the current particle is within the ghost particle margin of the half plane reached by the given step and along an axis + let is_in_ghost_margin_single_dim = |step: i8, axis: usize| -> bool { + match step { + -1 => particle[axis] - min_corner[axis] < ghost_particle_margin, + 0 => true, + 1 => max_corner[axis] - particle[axis] < ghost_particle_margin, + _ => unsafe { std::hint::unreachable_unchecked() }, + } + }; + + // Checks whether the current particle is within the ghost particle margin of the neighbor subdomain reached by the given steps + let is_in_ghost_margin = |x_step: i8, y_step: i8, z_step: i8| -> bool { + is_in_ghost_margin_single_dim(x_step, 0) + && is_in_ghost_margin_single_dim(y_step, 1) + && is_in_ghost_margin_single_dim(z_step, 2) + }; + + // Loop over all 27 subdomains around and including the owning subdomain + for &i in &[-1, 0, 1] { + for &j in &[-1, 0, 1] { + for &k in &[-1, 0, 1] { + // Check if the particle is in the ghost particle margin of the current subdomain + let in_ghost_margin = is_in_ghost_margin(i, j, k); + + if in_ghost_margin { + let neighbor_subdomain_ijk = [ + subdomain_ijk[0] + to_index!(i), + subdomain_ijk[1] + to_index!(j), + subdomain_ijk[2] + to_index!(k), + ]; + // The potential neighbor subdomain might not even be part of our computation domain + if let Some(cell) = subdomain_grid.get_cell(neighbor_subdomain_ijk) { + // If it is, it can be added as a subdomain of the particle + subdomains.push(subdomain_grid.flatten_cell_index(&cell)); + } + } + } + } + } + } +} + +pub(crate) mod debug { + use super::*; + + /// Prints statistics of the given list of subdomains + #[allow(unused)] + pub(crate) fn subdomain_stats( + parameters: &ParametersSubdomainGrid, + particles: &[Vector3], + subdomains: &Subdomains, + ) { + profile!("subdomain stats"); + info!("Statistics"); + + let per_subdomain_particle_count = subdomains + .per_subdomain_particles + .iter() + .map(|p| p.len()) + .collect::>(); + + info!( + "Occupied subdomains: {}", + per_subdomain_particle_count.len() + ); + + /* + info!("Printing subdomain particle counts:"); + for subdomain in subdomains.per_subdomain_particles { + info!("{}", subdomain.lock().len()); + } + */ + + info!("Smallest Subdomains:"); + for i in 0..11 { + let c = per_subdomain_particle_count + .iter() + .copied() + .filter(|c| *c == i) + .count(); + + info!( + "Number of subdomains with {} particles: {} ({:.2}% of number of subdomains)", + i, + c, + (c as f64 / per_subdomain_particle_count.len() as f64) * 100.0 + ); + } + + info!("Other stats:"); + for n in [10, 50, 100, 500, 1000, 2000, 10000] { + let c = per_subdomain_particle_count + .iter() + .copied() + .filter(|c| *c <= n) + .count(); + + info!( + "Number of subdomains with {} or fewer particles: {} ({:.2}% of number of subdomains)", + n, + c, + (c as f64 / per_subdomain_particle_count.len() as f64) * 100.0 + ); + } + + { + let mut per_subdomain_particle_count = per_subdomain_particle_count.clone(); + per_subdomain_particle_count.sort(); + + { + let largest_subdomain_particle_count = + per_subdomain_particle_count[per_subdomain_particle_count.len() - 1]; + info!("Largest subdomain has {largest_subdomain_particle_count} particles"); + + for f in [0.95, 0.9, 0.8, 0.75, 0.5, 0.1] { + let c = per_subdomain_particle_count + .iter() + .copied() + .filter(|c| *c >= (f * largest_subdomain_particle_count as f64) as usize) + .count(); + + let n = per_subdomain_particle_count + .iter() + .copied() + .filter(|c| *c >= (f * largest_subdomain_particle_count as f64) as usize) + .sum::(); + + info!( + "Number of subdomains with {} or more particles ({}% of largest subdomain): {} ({:.2}% of number of subdomains), in sum {} particles ({:.2}% of all particles)", + (f * largest_subdomain_particle_count as f64) as usize, + f * 100.0, + c, + (c as f64 / per_subdomain_particle_count.len() as f64) * 100.0, + n, + 100.0 * (n as f64 / particles.len() as f64) + ); + } + } + + info!("Largest subdomains:"); + for i in 0..10 { + if let Some(&count) = + per_subdomain_particle_count.get(per_subdomain_particle_count.len() - 1 - i) + { + info!( + "{} particles ({:.2}% of all particles)", + count, + 100.0 * (count as f64 / particles.len() as f64) + ); + } + } + } + } + + #[allow(unused)] + pub(crate) fn subdomains_to_hexmesh( + parameters: &ParametersSubdomainGrid, + subdomains: &Subdomains, + ) -> HexMesh3d { + let mut hexmesh = HexMesh3d::default(); + let subdomain_grid = ¶meters.subdomain_grid; + + // Loop over all non-empty subdomains + for &flat_subdomain_idx in &subdomains.flat_subdomain_indices { + let subdomain_ijk = subdomain_grid + .try_unflatten_cell_index(flat_subdomain_idx as I) + .unwrap(); + let [i, j, k] = subdomain_ijk.index().clone(); + + let vertex_offset = hexmesh.vertices.len(); + + { + let mut push_vertex = |a: i32, b: i32, c: i32| { + hexmesh.vertices.push( + subdomain_grid.point_coordinates( + &subdomain_grid + .get_point([i + to_index!(a), j + to_index!(b), k + to_index!(c)]) + .unwrap(), + ), + ); + }; + + push_vertex(0, 0, 0); + push_vertex(1, 0, 0); + push_vertex(1, 1, 0); + push_vertex(0, 1, 0); + push_vertex(0, 0, 1); + push_vertex(1, 0, 1); + push_vertex(1, 1, 1); + push_vertex(0, 1, 1); + } + + hexmesh.cells.push([ + vertex_offset + 0, + vertex_offset + 1, + vertex_offset + 2, + vertex_offset + 3, + vertex_offset + 4, + vertex_offset + 5, + vertex_offset + 6, + vertex_offset + 7, + ]); + } + + hexmesh + } + + /// Counts subdomains with only ghost particles and no particles in interior + #[allow(unused)] + pub(crate) fn count_no_owned_particles_subdomains( + parameters: &ParametersSubdomainGrid, + particles: &[Vector3], + subdomains: &Subdomains, + ) -> usize { + profile!(parent, "count_no_owned_particles_subdomains"); + + let no_owned_particles_counter = AtomicUsize::new(0); + + subdomains + .flat_subdomain_indices + .par_iter() + .copied() + .zip(subdomains.per_subdomain_particles.par_iter()) + .for_each(|(flat_subdomain_idx, subdomain_particle_indices)| { + profile!("inner subdomain_tasks", parent = parent); + + let flat_subdomain_idx: I = flat_subdomain_idx; + let subdomain_particle_indices: &[usize] = subdomain_particle_indices.as_slice(); + + // Collect all particle positions of this subdomain + let subdomain_particles = subdomain_particle_indices + .iter() + .copied() + .map(|idx| particles[idx]) + .collect::>(); + + // Get the cell index and AABB of the subdomain + let subdomain_idx = parameters + .subdomain_grid + .try_unflatten_cell_index(flat_subdomain_idx) + .expect("Subdomain cell does not exist"); + let subdomain_aabb = parameters.subdomain_grid.cell_aabb(&subdomain_idx); + + // Count the number of owned (non-ghost particles) of this domain + let non_ghost_particle_count = subdomain_particles + .iter() + .filter(|p| subdomain_aabb.contains_point(*p)) + .count(); + if non_ghost_particle_count == 0 { + no_owned_particles_counter.fetch_add(1, Ordering::AcqRel); + } + }); + + let no_owned_particles_counter = no_owned_particles_counter.into_inner(); + no_owned_particles_counter + } +} + +/// Performs integer division and rounds the result up if there is a remainder +fn int_ceil_div(numerator: T, denominator: T) -> T { + numerator / denominator + (numerator % denominator).min(T::one()) +} + +/// Ensures that at least the specified total capacity is reserved for the given vector +fn reserve_total(vec: &mut Vec, total_capacity: usize) { + if total_capacity > vec.capacity() { + vec.reserve(total_capacity - vec.capacity()); + } +} + +/// Gathers particle related data from global storage to subdomain storage +fn gather_subdomain_data( + global_data: &[T], + subdomain_particle_indices: &[usize], + subdomain_data: &mut Vec, +) { + subdomain_data.clear(); + reserve_total(subdomain_data, subdomain_particle_indices.len()); + subdomain_data.extend( + subdomain_particle_indices + .iter() + .copied() + .map(|idx| global_data[idx]), + ); +} diff --git a/splashsurf_lib/src/lib.rs b/splashsurf_lib/src/lib.rs index 78424e7..b35769d 100644 --- a/splashsurf_lib/src/lib.rs +++ b/splashsurf_lib/src/lib.rs @@ -27,6 +27,7 @@ use log::info; /// Re-export the version of `nalgebra` used by this crate pub use nalgebra; use nalgebra::Vector3; +use std::hash::Hash; use thiserror::Error as ThisError; /// Re-export the version of `vtkio` used by this crate, if vtk support is enabled #[cfg(feature = "vtk_extras")] @@ -52,6 +53,7 @@ pub mod profiling; pub mod profiling_macro; mod aabb; +pub(crate) mod dense_subdomains; pub mod density_map; pub mod generic_tree; #[cfg(feature = "io")] @@ -62,7 +64,8 @@ pub mod marching_cubes; pub mod mesh; pub mod neighborhood_search; pub mod octree; -mod reconstruction; +pub mod reconstruction; +mod reconstruction_octree; pub mod sph_interpolation; pub mod topology; mod traits; @@ -107,6 +110,9 @@ pub(crate) fn new_map() -> MapType { */ pub(crate) type ParallelMapType = dashmap::DashMap; +fn new_parallel_map() -> ParallelMapType { + ParallelMapType::with_hasher(HashState::default()) +} /// Parameters for the spatial decomposition #[derive(Clone, Debug)] @@ -196,6 +202,8 @@ pub struct Parameters { pub domain_aabb: Option>, /// Whether to allow multi threading within the surface reconstruction procedure pub enable_multi_threading: bool, + /// Each subdomain will be a cube consisting of this number of MC cube cells along each coordinate axis + pub subdomain_num_cubes_per_dim: Option, /// Parameters for the spatial decomposition (octree subdivision) of the particles. /// If not provided, no octree is generated and a global approach is used instead. pub spatial_decomposition: Option>, @@ -212,6 +220,7 @@ impl Parameters { iso_surface_threshold: self.iso_surface_threshold.try_convert()?, domain_aabb: map_option!(&self.domain_aabb, aabb => aabb.try_convert()?), enable_multi_threading: self.enable_multi_threading, + subdomain_num_cubes_per_dim: self.subdomain_num_cubes_per_dim, spatial_decomposition: map_option!(&self.spatial_decomposition, sd => sd.try_convert()?), }) } @@ -356,14 +365,24 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>( output_surface.grid.log_grid_info(); - if parameters.spatial_decomposition.is_some() { - reconstruction::reconstruct_surface_domain_decomposition( + if parameters.subdomain_num_cubes_per_dim.is_some() { + reconstruction::reconstruct_surface_subdomain_grid::( + particle_positions, + parameters, + output_surface, + )?; + } else if parameters.spatial_decomposition.is_some() { + reconstruction_octree::reconstruct_surface_domain_decomposition( particle_positions, parameters, output_surface, )?; } else { - reconstruction::reconstruct_surface_global(particle_positions, parameters, output_surface)?; + reconstruction_octree::reconstruct_surface_global( + particle_positions, + parameters, + output_surface, + )?; } Ok(()) diff --git a/splashsurf_lib/src/marching_cubes.rs b/splashsurf_lib/src/marching_cubes.rs index 1b91981..0c4c8d2 100644 --- a/splashsurf_lib/src/marching_cubes.rs +++ b/splashsurf_lib/src/marching_cubes.rs @@ -200,6 +200,7 @@ pub fn check_mesh_consistency( grid: &UniformGrid, mesh: &TriMesh3d, ) -> Result<(), String> { + profile!("check_mesh_consistency"); let boundary_edges = mesh.find_boundary_edges(); if boundary_edges.is_empty() { diff --git a/splashsurf_lib/src/reconstruction.rs b/splashsurf_lib/src/reconstruction.rs index 4bfe66d..ac96b91 100644 --- a/splashsurf_lib/src/reconstruction.rs +++ b/splashsurf_lib/src/reconstruction.rs @@ -1,733 +1,60 @@ -//! 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, Parameters, - ParticleDensityComputationStrategy, Real, ReconstructionError, SpatialDecompositionParameters, - SurfaceReconstruction, -}; -use log::{debug, info, trace}; +use log::info; 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(()) -} +use crate::dense_subdomains::{ + compute_global_density_vector, decomposition, initialize_parameters, reconstruction, stitching, + subdomain_classification::GhostMarginClassifier, +}; +use crate::{profile, Index, Parameters, Real, SurfaceReconstruction}; -/// Performs a surface reconstruction with an octree for domain decomposition -pub(crate) fn reconstruct_surface_domain_decomposition<'a, I: Index, R: Real>( +/// Performs a surface reconstruction with a regular grid for domain decomposition +pub(crate) fn reconstruct_surface_subdomain_grid<'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: SpatialDecompositionParameters, - /// 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(); - - // Construct the octree - let octree = if let Some(decomposition_parameters) = ¶meters.spatial_decomposition { - 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, - ) - } else { - // TODO: Use default values instead? - - // If there are no decomposition parameters, we cannot construct an octree. - return None; - }; - - // 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: parameters.spatial_decomposition.as_ref().unwrap().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, +) -> Result<(), anyhow::Error> { + let mesh = { + profile!("surface reconstruction prototype"); + + let parameters = initialize_parameters(parameters, &particle_positions, output_surface)?; + + // Filter "narrow band" + /* + let narrow_band_particles = extract_narrow_band(¶meters, &particles); + let particles = narrow_band_particles; + */ + + let subdomains = + decomposition::>(¶meters, &particle_positions)?; + + /* + subdomain_stats(¶meters, &particle_positions, &subdomains); + info!( + "Number of subdomains with only ghost particles: {}", + count_no_owned_particles_subdomains(¶meters, &particle_positions, &subdomains) ); - 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 - }; + let particle_densities = + compute_global_density_vector(¶meters, &particle_positions, &subdomains); - // 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]), + let surface_patches = reconstruction( + ¶meters, + &particle_positions, + &particle_densities, + &subdomains, ); - } - /// 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]), + let global_mesh = stitching(surface_patches); + info!( + "Global mesh has {} vertices and {} triangles.", + global_mesh.vertices.len(), + global_mesh.triangles.len() ); - } -} - -/// 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() + global_mesh }; - // 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, - )?; - + let _ = std::mem::replace(&mut output_surface.mesh, 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/reconstruction_octree.rs b/splashsurf_lib/src/reconstruction_octree.rs new file mode 100644 index 0000000..4bfe66d --- /dev/null +++ b/splashsurf_lib/src/reconstruction_octree.rs @@ -0,0 +1,733 @@ +//! 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, Parameters, + ParticleDensityComputationStrategy, Real, ReconstructionError, SpatialDecompositionParameters, + 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: SpatialDecompositionParameters, + /// 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(); + + // Construct the octree + let octree = if let Some(decomposition_parameters) = ¶meters.spatial_decomposition { + 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, + ) + } else { + // TODO: Use default values instead? + + // If there are no decomposition parameters, we cannot construct an octree. + return None; + }; + + // 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: parameters.spatial_decomposition.as_ref().unwrap().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/traits.rs b/splashsurf_lib/src/traits.rs index bdcd31b..8bdb185 100644 --- a/splashsurf_lib/src/traits.rs +++ b/splashsurf_lib/src/traits.rs @@ -5,12 +5,40 @@ use std::ops::{AddAssign, MulAssign, SubAssign}; use bytemuck::Pod; use nalgebra::{RealField, SVector}; use num_integer::Integer; -use num_traits::{Bounded, CheckedAdd, CheckedMul, CheckedSub, FromPrimitive, ToPrimitive}; +use num_traits::{ + Bounded, CheckedAdd, CheckedMul, CheckedSub, FromPrimitive, NumCast, SaturatingSub, ToPrimitive, +}; /// Convenience trait that combines `Send` and `Sync` pub trait ThreadSafe: Sync + Send {} impl ThreadSafe for T where T: Sync + Send {} +pub struct IndexRange { + start: I, + end: I, +} + +impl IndexRange { + fn new(start: I, end: I) -> Self { + assert!(start <= end, "start must be less or equal to end"); + Self { start, end } + } + + pub fn iter(self) -> impl Iterator { + let end = self.end; + let mut counter = self.start; + std::iter::from_fn(move || { + let current = counter; + if current < end { + counter += I::one(); + Some(current) + } else { + None + } + }) + } +} + /// Trait that has to be implemented for types to be used as background grid cell indices in the context of the library pub trait Index: Copy @@ -20,11 +48,13 @@ pub trait Index: + CheckedAdd + CheckedSub + CheckedMul + + SaturatingSub + AddAssign + SubAssign + MulAssign + FromPrimitive + ToPrimitive + + NumCast + Default + Debug + Display @@ -32,6 +62,14 @@ pub trait Index: + ThreadSafe + 'static { + fn range(start: Self, end: Self) -> IndexRange { + IndexRange::new(start, end) + } + + fn two() -> Self { + Self::one() + Self::one() + } + /// Converts this value to the specified [`Real`] type `T` by converting first to `f64` followed by `T::from_f64`. If the value cannot be represented by the target type, `None` is returned. fn to_real(self) -> Option { R::from_f64(self.to_f64()?) @@ -46,21 +84,37 @@ pub trait Index: fn times(self, n: i32) -> Self { self.mul(Self::from_i32(n).unwrap()) } + + /// Returns the squared value of this value. + fn squared(self) -> Self { + self * self + } + + /// Returns the cubed value of this value. + fn cubed(self) -> Self { + self * self * self + } + + fn checked_cubed(self) -> Option { + self.checked_mul(&self) + .and_then(|val| val.checked_mul(&self)) + } } /// Trait that has to be implemented for types to be used as floating points values in the context of the library (e.g. for coordinates, density values) pub trait Real: - RealField - // Required by RStar and not part of RealFied anymore - + Bounded - // Not part of RealField anymore - + Copy - + FromPrimitive - + ToPrimitive - + Debug - + Default - + Pod - + ThreadSafe +RealField +// Required by RStar and not part of RealFied anymore ++ Bounded +// Not part of RealField anymore ++ Copy ++ FromPrimitive ++ ToPrimitive ++ NumCast ++ Debug ++ Default ++ Pod ++ ThreadSafe { /// Tries to convert this value to another [`Real`] type `T` by converting first to `f64` followed by `T::from_f64`. If the value cannot be represented by the target type, `None` is returned. fn try_convert(self) -> Option { @@ -69,8 +123,8 @@ pub trait Real: /// Tries to convert the values of a statically sized `nalgebra::SVector` to another type, same behavior as [`Real::try_convert`] fn try_convert_vec_from(vec: &SVector) -> Option> - where - R: Real, + where + R: Real, { let mut converted = SVector::::zeros(); for i in 0..D { @@ -108,11 +162,13 @@ impl Index for T where + CheckedAdd + CheckedSub + CheckedMul + + SaturatingSub + AddAssign + SubAssign + MulAssign + FromPrimitive + ToPrimitive + + NumCast + Debug + Default + Display @@ -128,6 +184,7 @@ impl< + Copy + FromPrimitive + ToPrimitive + + NumCast + Debug + Default + Pod diff --git a/splashsurf_lib/src/uniform_grid.rs b/splashsurf_lib/src/uniform_grid.rs index ec51dc7..998e996 100644 --- a/splashsurf_lib/src/uniform_grid.rs +++ b/splashsurf_lib/src/uniform_grid.rs @@ -382,6 +382,19 @@ impl UniformCartesianCubeGrid3d { } } + pub fn get_edge(&self, origin_ijk: [I; 3], axis: Axis) -> Option> { + let mut target_ijk = origin_ijk.clone(); + target_ijk[axis.dim()] += I::one(); + if self.point_exists(&origin_ijk) && self.point_exists(&target_ijk) { + Some(EdgeIndex { + origin: PointIndex::from_ijk(origin_ijk), + axis, + }) + } else { + None + } + } + /// Returns whether a point exists in the grid #[inline(always)] pub fn point_exists(&self, point_ijk: &[I; 3]) -> bool { @@ -412,6 +425,15 @@ impl UniformCartesianCubeGrid3d { || cell_index.index[2] + I::one() == self.n_cells_per_dim[2]) } + /// Returns whether the edge is between two points on the boundary of the grid + pub fn is_boundary_edge(&self, edge_index: &EdgeIndex) -> bool { + let origin_index = edge_index.origin().index(); + edge_index.axis.orthogonal_axes().iter().any(|orth_ax| { + origin_index[orth_ax.dim()] == I::zero() + || origin_index[orth_ax.dim()] + I::one() == self.n_points_per_dim[orth_ax.dim()] + }) + } + /// Flattens the grid point index triplet to a single index #[inline(always)] pub fn flatten_point_indices(&self, i: I, j: I, k: I) -> I { @@ -769,7 +791,7 @@ impl UniformCartesianCubeGrid3d { impl OwningSubdomainGrid { /// Creates a new subdomain grid - pub(crate) fn new( + pub fn new( global_grid: UniformGrid, subdomain_grid: UniformGrid, subdomain_offset: [I; 3], @@ -964,6 +986,11 @@ impl EdgeIndex { .expect("Index type overflow"); PointIndex::from_ijk(new_index) } + + /// The axis this edge is parallel to + pub fn axis(&self) -> Axis { + self.axis + } } #[test] diff --git a/splashsurf_lib/tests/integration_tests/test_full.rs b/splashsurf_lib/tests/integration_tests/test_full.rs index 519ad67..0d7d7dc 100644 --- a/splashsurf_lib/tests/integration_tests/test_full.rs +++ b/splashsurf_lib/tests/integration_tests/test_full.rs @@ -14,6 +14,7 @@ enum Strategy { Global, Octree, OctreeStitching, + SubdomainGrid, } fn params_with_aabb( @@ -35,6 +36,7 @@ fn params_with_aabb( iso_surface_threshold, domain_aabb, enable_multi_threading: false, + subdomain_num_cubes_per_dim: None, spatial_decomposition: None, }; @@ -58,6 +60,9 @@ fn params_with_aabb( ParticleDensityComputationStrategy::SynchronizeSubdomains, }); } + Strategy::SubdomainGrid => { + parameters.subdomain_num_cubes_per_dim = Some(64); + } } parameters @@ -156,15 +161,22 @@ 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.vtk" => "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.vtk" => "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.vtk" => "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.vtk" => "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.vtk" => "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.vtk" => "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); generate_test!(f32, surface_reconstruction_free_particles_02, "free_particles_125_particles.vtk" => "reconstruct_surface_free_particles_02_global.vtk", params_with_aabb(0.5, 4.0, 1.5, 0.45, Some(Aabb3d::new(Vector3::new(-10.0, -10.0, -10.0), Vector3::new(210.0, 210.0, 210.0))), Strategy::Global), 1500, 1600); diff --git a/splashsurf_lib/tests/integration_tests/test_octree.rs b/splashsurf_lib/tests/integration_tests/test_octree.rs index 5eebee0..933e50c 100644 --- a/splashsurf_lib/tests/integration_tests/test_octree.rs +++ b/splashsurf_lib/tests/integration_tests/test_octree.rs @@ -175,6 +175,7 @@ impl Default for TestParameters { } } +#[allow(unused)] impl TestParameters { fn new(particle_radius: f64, compact_support_factor: f64, cube_size_factor: f64) -> Self { let params = Self::default();