From b5c634ab1d913d24f2a7df567e2662bd7732042c Mon Sep 17 00:00:00 2001 From: tmathis720 <34523436+tmathis720@users.noreply.github.com> Date: Fri, 25 Oct 2024 23:34:15 -0700 Subject: [PATCH] Improved the function and concurrency in the domain module by using DashMap in some parts of the code. Improved the geometry module and added calculations for normal vectors for shapes. All tests passing. --- Cargo.lock | 66 ++++++- Cargo.toml | 1 + src/boundary/bc_handler.rs | 212 +++++---------------- src/boundary/dirichlet.rs | 132 +++---------- src/boundary/neumann.rs | 124 +++---------- src/boundary/robin.rs | 112 +++-------- src/domain/entity_fill.rs | 33 ++-- src/domain/mesh/geometry.rs | 97 ++-------- src/domain/overlap.rs | 218 +++++----------------- src/domain/section.rs | 73 ++++---- src/domain/sieve.rs | 276 +++++++++++----------------- src/domain/stratify.rs | 24 ++- src/geometry/docs/about_geometry.md | 78 ++++---- src/geometry/hexahedron.rs | 56 ++---- src/geometry/mod.rs | 206 ++++++--------------- src/geometry/prism.rs | 48 +---- src/geometry/pyramid.rs | 164 +++-------------- src/geometry/quadrilateral.rs | 76 +++++--- src/geometry/tetrahedron.rs | 28 +-- src/geometry/triangle.rs | 68 +++++-- target/.rustc_info.json | 2 +- 21 files changed, 673 insertions(+), 1421 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index fe9326c..63e95c9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -270,6 +270,19 @@ dependencies = [ "typenum", ] +[[package]] +name = "dashmap" +version = "5.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "978747c1d849a7d2ee5e8adc0159961c48fb7e5db2f06af6723b80123bb53856" +dependencies = [ + "cfg-if", + "hashbrown 0.14.5", + "lock_api", + "once_cell", + "parking_lot_core", +] + [[package]] name = "dbgf" version = "0.1.2" @@ -554,6 +567,12 @@ version = "0.12.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "8a9ee70c43aaf417c914396645a0fa852624801b24ebb7ae78fe8272889ac888" +[[package]] +name = "hashbrown" +version = "0.14.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" + [[package]] name = "heck" version = "0.5.0" @@ -575,6 +594,7 @@ version = "0.3.0" dependencies = [ "criterion", "crossbeam", + "dashmap", "faer", "rand", "rayon", @@ -588,7 +608,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "bd070e393353796e801d209ad339e89596eb4c8d430d18ede6a1cced8fafbd99" dependencies = [ "autocfg", - "hashbrown", + "hashbrown 0.12.3", ] [[package]] @@ -633,6 +653,16 @@ version = "0.2.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" +[[package]] +name = "lock_api" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "07af8b9cdd281b7915f413fa73f29ebd5d55d0d3f0155584dade1ff18cea1b17" +dependencies = [ + "autocfg", + "scopeguard", +] + [[package]] name = "log" version = "0.4.22" @@ -800,6 +830,19 @@ version = "6.6.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e2355d85b9a3786f481747ced0e0ff2ba35213a1f9bd406ed906554d7af805a1" +[[package]] +name = "parking_lot_core" +version = "0.9.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e401f977ab385c9e4e3ab30627d6f26d00e2c73eef317493c4ec6d468726cf8" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + [[package]] name = "paste" version = "1.0.15" @@ -1006,6 +1049,15 @@ version = "0.5.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "03251193000f4bd3b042892be858ee50e8b3719f2b08e5833ac4353724632430" +[[package]] +name = "redox_syscall" +version = "0.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b6dfecf2c74bce2466cabf93f6664d6998a69eb21e39f4207930065b27b771f" +dependencies = [ + "bitflags 2.6.0", +] + [[package]] name = "regex" version = "1.10.6" @@ -1056,6 +1108,12 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + [[package]] name = "seq-macro" version = "0.3.5" @@ -1105,6 +1163,12 @@ dependencies = [ "digest", ] +[[package]] +name = "smallvec" +version = "1.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" + [[package]] name = "syn" version = "2.0.77" diff --git a/Cargo.toml b/Cargo.toml index 7bc6829..fd3da1c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,6 +12,7 @@ rustc-hash = "2.0.0" faer = { version = "0.19", features = ["rayon"] } rayon = "1.5" crossbeam = "0.8.4" +dashmap = "5.1" [dev-dependencies] criterion = "0.4" diff --git a/src/boundary/bc_handler.rs b/src/boundary/bc_handler.rs index a1c171b..726746f 100644 --- a/src/boundary/bc_handler.rs +++ b/src/boundary/bc_handler.rs @@ -1,7 +1,7 @@ -use rustc_hash::FxHashMap; +use dashmap::DashMap; +use rayon::prelude::*; use std::sync::Arc; use crate::domain::mesh_entity::MeshEntity; -use crate::domain::section::Section; use crate::boundary::dirichlet::DirichletBC; use crate::boundary::neumann::NeumannBC; use crate::boundary::robin::RobinBC; @@ -10,20 +10,7 @@ use faer::MatMut; pub type BoundaryConditionFn = Arc f64 + Send + Sync>; /// BoundaryCondition represents various types of boundary conditions -/// that can be applied to mesh entities. These include: -/// * Dirichlet: A constant boundary condition where a specific value is enforced. -/// * Neumann: A flux or derivative boundary condition. -/// * Robin: A linear combination of Dirichlet and Neumann boundary conditions. -/// * DirichletFn: A time-dependent Dirichlet boundary condition using a function. -/// * NeumannFn: A time-dependent Neumann boundary condition using a function. -/// Example usage: -/// -/// let bc_handler = BoundaryConditionHandler::new(); -/// let entity = MeshEntity::Vertex(1); -/// let condition = BoundaryCondition::Dirichlet(10.0); -/// bc_handler.set_bc(entity, condition); -/// - +/// that can be applied to mesh entities. #[derive(Clone)] pub enum BoundaryCondition { Dirichlet(f64), @@ -34,150 +21,65 @@ pub enum BoundaryCondition { } /// The BoundaryConditionHandler struct is responsible for managing -/// boundary conditions associated with specific mesh entities. -/// It supports adding, retrieving, and applying boundary conditions -/// for a given set of mesh entities and interacting with the linear system -/// matrices and right-hand side vectors. -/// Example usage: -/// -/// let handler = BoundaryConditionHandler::new(); -/// let entity = MeshEntity::Edge(5); -/// let bc = BoundaryCondition::Neumann(5.0); -/// handler.set_bc(entity, bc); -/// - +/// boundary conditions associated with specific mesh entities. pub struct BoundaryConditionHandler { - conditions: Section, + conditions: DashMap, } impl BoundaryConditionHandler { - /// Creates a new BoundaryConditionHandler with an empty section - /// to store boundary conditions. - /// Example usage: - /// - /// let handler = BoundaryConditionHandler::new(); - /// - + /// Creates a new BoundaryConditionHandler with an empty map to store boundary conditions. pub fn new() -> Self { Self { - conditions: Section::new(), + conditions: DashMap::new(), } } /// Sets a boundary condition for a specific mesh entity. - /// The boundary condition can be Dirichlet, Neumann, Robin, or a functional form. - /// - /// # Arguments: - /// * `entity` - The mesh entity (such as a vertex, edge, or face). - /// * `condition` - The boundary condition to apply to the entity. - /// - /// Example usage: - /// - /// let entity = MeshEntity::Vertex(3); - /// let bc = BoundaryCondition::Dirichlet(1.0); - /// handler.set_bc(entity, bc); - /// pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) { - self.conditions.set_data(entity, condition); + self.conditions.insert(entity, condition); } /// Retrieves the boundary condition applied to a specific mesh entity, if it exists. - /// - /// # Arguments: - /// * `entity` - A reference to the mesh entity. - /// - /// Returns an `Option` indicating the boundary condition - /// if one has been set for the entity. - /// - /// Example usage: - /// - /// let entity = MeshEntity::Edge(2); - /// if let Some(bc) = handler.get_bc(&entity) { - /// println!("Boundary condition found!"); - /// } - /// pub fn get_bc(&self, entity: &MeshEntity) -> Option { - self.conditions.restrict(entity) + self.conditions.get(entity).map(|entry| entry.clone()) } /// Applies the boundary conditions to the system matrices and right-hand side vectors. - /// This modifies the matrix and the right-hand side based on the type of boundary - /// conditions for the provided boundary entities. - /// - /// # Arguments: - /// * `matrix` - The mutable matrix (MatMut) to modify. - /// * `rhs` - The mutable right-hand side vector (MatMut). - /// * `boundary_entities` - A list of mesh entities representing the boundary. - /// * `entity_to_index` - A hash map that associates each mesh entity with an index - /// in the system. - /// * `time` - The current time for time-dependent boundary conditions. - /// - /// Example usage: - /// - /// let mut matrix = MatMut::new(); - /// let mut rhs = MatMut::new(); - /// let boundary_entities = vec![MeshEntity::Vertex(1), MeshEntity::Edge(2)]; - /// let entity_to_index: FxHashMap = FxHashMap::default(); - /// handler.apply_bc(&mut matrix, &mut rhs, &boundary_entities, &entity_to_index, 0.0); - /// pub fn apply_bc( &self, matrix: &mut MatMut, rhs: &mut MatMut, boundary_entities: &[MeshEntity], - entity_to_index: &FxHashMap, + entity_to_index: &DashMap, time: f64, ) { for entity in boundary_entities { if let Some(bc) = self.get_bc(entity) { + let index = *entity_to_index.get(entity).unwrap(); match bc { BoundaryCondition::Dirichlet(value) => { let dirichlet_bc = DirichletBC::new(); - dirichlet_bc.apply_constant_dirichlet( - matrix, - rhs, - *entity_to_index.get(entity).unwrap(), - value, - ); + dirichlet_bc.apply_constant_dirichlet(matrix, rhs, index, value); } BoundaryCondition::Neumann(flux) => { let neumann_bc = NeumannBC::new(); - neumann_bc.apply_constant_neumann( - rhs, - *entity_to_index.get(entity).unwrap(), - flux, - ); + neumann_bc.apply_constant_neumann(rhs, index, flux); } BoundaryCondition::Robin { alpha, beta } => { let robin_bc = RobinBC::new(); - robin_bc.apply_robin( - matrix, - rhs, - *entity_to_index.get(entity).unwrap(), - alpha, - beta, - ); + robin_bc.apply_robin(matrix, rhs, index, alpha, beta); } BoundaryCondition::DirichletFn(fn_bc) => { let coords = [0.0, 0.0, 0.0]; let value = fn_bc(time, &coords); let dirichlet_bc = DirichletBC::new(); - dirichlet_bc.apply_constant_dirichlet( - matrix, - rhs, - *entity_to_index.get(entity).unwrap(), - value, - ); + dirichlet_bc.apply_constant_dirichlet(matrix, rhs, index, value); } BoundaryCondition::NeumannFn(fn_bc) => { let coords = [0.0, 0.0, 0.0]; let value = fn_bc(time, &coords); let neumann_bc = NeumannBC::new(); - neumann_bc.apply_constant_neumann( - rhs, - *entity_to_index.get(entity).unwrap(), - value, - ); + neumann_bc.apply_constant_neumann(rhs, index, value); } } } @@ -186,80 +88,52 @@ impl BoundaryConditionHandler { } /// The BoundaryConditionApply trait defines the `apply` method, which is used to apply -/// a boundary condition to a given mesh entity. It modifies the matrix and right-hand -/// side (rhs) of the system based on the boundary condition type (Dirichlet, Neumann, Robin, etc.). -/// -/// The `apply` method should be implemented by any type that represents boundary conditions -/// and needs to be applied to the system matrix and rhs. -/// -/// Example usage: -/// -/// let bc = BoundaryCondition::Dirichlet(1.0); -/// let entity = MeshEntity::Vertex(1); -/// bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.0); -/// +/// a boundary condition to a given mesh entity. pub trait BoundaryConditionApply { - /// Applies a boundary condition to a specific mesh entity, modifying the system matrix and rhs. - /// - /// # Arguments: - /// * `entity` - The mesh entity to which the boundary condition is applied. - /// * `rhs` - The mutable right-hand side vector that may be modified. - /// * `matrix` - The mutable system matrix that may be modified. - /// * `entity_to_index` - A hash map that associates each mesh entity with its index in the system. - /// * `time` - The current time for applying time-dependent boundary conditions. - /// - /// Example usage: - /// - /// let bc = BoundaryCondition::Neumann(5.0); - /// let entity = MeshEntity::Edge(3); - /// bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.0); - /// - fn apply(&self, entity: &MeshEntity, rhs: &mut MatMut, matrix: &mut MatMut, entity_to_index: &FxHashMap, time: f64); + fn apply( + &self, + entity: &MeshEntity, + rhs: &mut MatMut, + matrix: &mut MatMut, + entity_to_index: &DashMap, + time: f64, + ); } -/// The implementation of BoundaryConditionApply for BoundaryCondition allows different types -/// of boundary conditions (Dirichlet, Neumann, Robin, and their functional forms) to modify -/// the matrix and rhs according to the specific logic of each boundary condition type. -/// -/// Example usage: -/// -/// let bc = BoundaryCondition::DirichletFn(Arc::new(|time, coords| time + coords[0])); -/// let entity = MeshEntity::Vertex(2); -/// bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.5); -/// impl BoundaryConditionApply for BoundaryCondition { - /// Applies the boundary condition to a specific mesh entity. - /// - /// # Arguments: - /// * `entity` - The mesh entity to which the boundary condition is applied. - /// * `rhs` - The mutable right-hand side vector. - /// * `matrix` - The mutable system matrix. - /// * `entity_to_index` - A hash map mapping each mesh entity to its index in the system. - /// * `time` - The current time for time-dependent boundary conditions. - fn apply(&self, entity: &MeshEntity, rhs: &mut MatMut, matrix: &mut MatMut, entity_to_index: &FxHashMap, time: f64) { + fn apply( + &self, + entity: &MeshEntity, + rhs: &mut MatMut, + matrix: &mut MatMut, + entity_to_index: &DashMap, + time: f64, + ) { + let index = *entity_to_index.get(entity).unwrap(); match self { BoundaryCondition::Dirichlet(value) => { let dirichlet_bc = DirichletBC::new(); - dirichlet_bc.apply_constant_dirichlet(matrix, rhs, *entity_to_index.get(entity).unwrap(), *value); + dirichlet_bc.apply_constant_dirichlet(matrix, rhs, index, *value); } BoundaryCondition::Neumann(flux) => { let neumann_bc = NeumannBC::new(); - neumann_bc.apply_constant_neumann(rhs, *entity_to_index.get(entity).unwrap(), *flux); + neumann_bc.apply_constant_neumann(rhs, index, *flux); } - BoundaryCondition::Robin { alpha: _, beta: _ } => { - // Robin boundary condition logic to be implemented here + BoundaryCondition::Robin { alpha, beta } => { + let robin_bc = RobinBC::new(); + robin_bc.apply_robin(matrix, rhs, index, *alpha, *beta); } BoundaryCondition::DirichletFn(fn_bc) => { - let coords = [0.0, 0.0, 0.0]; // Placeholder for mesh entity coordinates + let coords = [0.0, 0.0, 0.0]; let value = fn_bc(time, &coords); let dirichlet_bc = DirichletBC::new(); - dirichlet_bc.apply_constant_dirichlet(matrix, rhs, *entity_to_index.get(entity).unwrap(), value); + dirichlet_bc.apply_constant_dirichlet(matrix, rhs, index, value); } BoundaryCondition::NeumannFn(fn_bc) => { - let coords = [0.0, 0.0, 0.0]; // Placeholder for mesh entity coordinates + let coords = [0.0, 0.0, 0.0]; let value = fn_bc(time, &coords); let neumann_bc = NeumannBC::new(); - neumann_bc.apply_constant_neumann(rhs, *entity_to_index.get(entity).unwrap(), value); + neumann_bc.apply_constant_neumann(rhs, index, value); } } } diff --git a/src/boundary/dirichlet.rs b/src/boundary/dirichlet.rs index ca614d9..e2643d7 100644 --- a/src/boundary/dirichlet.rs +++ b/src/boundary/dirichlet.rs @@ -1,78 +1,43 @@ -// src/boundary/dirichlet.rs - +use dashmap::DashMap; +use std::sync::Arc; use crate::domain::mesh_entity::MeshEntity; -use rustc_hash::FxHashMap; use crate::boundary::bc_handler::{BoundaryCondition, BoundaryConditionApply}; -use crate::domain::section::Section; use faer::MatMut; + /// The `DirichletBC` struct represents a handler for applying Dirichlet boundary conditions -/// to a set of mesh entities. It stores the conditions in a `Section` structure and applies -/// them to modify both the system matrix and the right-hand side (rhs). -/// -/// Example usage: -/// -/// let dirichlet_bc = DirichletBC::new(); -/// let entity = MeshEntity::Vertex(1); -/// dirichlet_bc.set_bc(entity, BoundaryCondition::Dirichlet(10.0)); -/// dirichlet_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); -/// +/// to a set of mesh entities. It stores the conditions in a `DashMap` and applies them to +/// modify both the system matrix and the right-hand side (rhs). pub struct DirichletBC { - conditions: Section, + conditions: DashMap, } impl DirichletBC { - /// Creates a new instance of `DirichletBC` with an empty section to store boundary conditions. - /// - /// Example usage: - /// - /// let dirichlet_bc = DirichletBC::new(); - /// + /// Creates a new instance of `DirichletBC` with an empty `DashMap` to store boundary conditions. pub fn new() -> Self { Self { - conditions: Section::new(), + conditions: DashMap::new(), } } /// Sets a Dirichlet boundary condition for a specific mesh entity. - /// - /// # Arguments: - /// * `entity` - The mesh entity to which the boundary condition will be applied. - /// * `condition` - The boundary condition to set (either a constant value or a functional form). - /// - /// Example usage: - /// - /// let entity = MeshEntity::Vertex(1); - /// dirichlet_bc.set_bc(entity, BoundaryCondition::Dirichlet(5.0)); - /// pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) { - self.conditions.set_data(entity, condition); + self.conditions.insert(entity, condition); } /// Applies the stored Dirichlet boundary conditions to the system matrix and rhs. /// It iterates over the stored conditions and applies either constant or function-based Dirichlet /// boundary conditions to the corresponding entities. - /// - /// # Arguments: - /// * `matrix` - The mutable system matrix. - /// * `rhs` - The mutable right-hand side vector. - /// * `entity_to_index` - A hash map that associates mesh entities with their indices in the system. - /// * `time` - The current time, used for time-dependent boundary conditions. - /// - /// Example usage: - /// - /// let entity_to_index = FxHashMap::default(); - /// dirichlet_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); - /// pub fn apply_bc( &self, matrix: &mut MatMut, rhs: &mut MatMut, - entity_to_index: &FxHashMap, + entity_to_index: &DashMap, time: f64, ) { - let data = self.conditions.data.read().unwrap(); - for (entity, condition) in data.iter() { - if let Some(&index) = entity_to_index.get(entity) { + // Iterate through the conditions and apply each condition accordingly. + self.conditions.iter().for_each(|entry| { + let (entity, condition) = entry.pair(); + if let Some(index) = entity_to_index.get(entity).map(|i| *i) { match condition { BoundaryCondition::Dirichlet(value) => { self.apply_constant_dirichlet(matrix, rhs, index, *value); @@ -85,21 +50,10 @@ impl DirichletBC { _ => {} } } - } + }); } /// Applies a constant Dirichlet boundary condition to the matrix and rhs for a specific index. - /// - /// # Arguments: - /// * `matrix` - The mutable system matrix. - /// * `rhs` - The mutable right-hand side vector. - /// * `index` - The index of the matrix row and rhs corresponding to the mesh entity. - /// * `value` - The Dirichlet condition value to be applied. - /// - /// Example usage: - /// - /// dirichlet_bc.apply_constant_dirichlet(&mut matrix, &mut rhs, 1, 5.0); - /// pub fn apply_constant_dirichlet( &self, matrix: &mut MatMut, @@ -115,20 +69,7 @@ impl DirichletBC { rhs.write(index, 0, value); } - /// Retrieves the coordinates of the mesh entity. - /// - /// This method currently returns a default placeholder value, but it can be expanded - /// to extract real entity coordinates if needed. - /// - /// # Arguments: - /// * `_entity` - The mesh entity for which coordinates are being requested. - /// - /// Returns an array of default coordinates `[0.0, 0.0, 0.0]`. - /// - /// Example usage: - /// - /// let coords = dirichlet_bc.get_coordinates(&entity); - /// + /// Retrieves the coordinates of the mesh entity (placeholder for real coordinates). fn get_coordinates(&self, _entity: &MeshEntity) -> [f64; 3] { [0.0, 0.0, 0.0] } @@ -136,21 +77,14 @@ impl DirichletBC { impl BoundaryConditionApply for DirichletBC { /// Applies the stored Dirichlet boundary conditions for a specific mesh entity. - /// - /// This implementation utilizes the general `apply_bc` method to modify the matrix and rhs. - /// - /// # Arguments: - /// * `_entity` - The mesh entity to which the boundary condition applies. - /// * `rhs` - The mutable right-hand side vector. - /// * `matrix` - The mutable system matrix. - /// * `entity_to_index` - A hash map that associates mesh entities with their indices. - /// * `time` - The current time, used for time-dependent boundary conditions. - /// - /// Example usage: - /// - /// dirichlet_bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.0); - /// - fn apply(&self, _entity: &MeshEntity, rhs: &mut MatMut, matrix: &mut MatMut, entity_to_index: &FxHashMap, time: f64) { + fn apply( + &self, + _entity: &MeshEntity, + rhs: &mut MatMut, + matrix: &mut MatMut, + entity_to_index: &DashMap, + time: f64, + ) { self.apply_bc(matrix, rhs, entity_to_index, time); } } @@ -158,7 +92,6 @@ impl BoundaryConditionApply for DirichletBC { #[cfg(test)] mod tests { use super::*; - use rustc_hash::FxHashMap; use faer::Mat; use crate::domain::mesh_entity::MeshEntity; use std::sync::Arc; @@ -178,7 +111,7 @@ mod tests { dirichlet_bc.set_bc(entity, BoundaryCondition::Dirichlet(10.0)); // Verify that the condition was set correctly - let condition = dirichlet_bc.conditions.restrict(&entity); + let condition = dirichlet_bc.conditions.get(&entity).map(|entry| entry.clone()); assert!(matches!(condition, Some(BoundaryCondition::Dirichlet(10.0)))); } @@ -186,30 +119,24 @@ mod tests { fn test_apply_constant_dirichlet() { let dirichlet_bc = DirichletBC::new(); let entity = MeshEntity::Vertex(1); - let mut entity_to_index = FxHashMap::default(); + let entity_to_index = DashMap::new(); entity_to_index.insert(entity, 1); - // Set a Dirichlet boundary condition dirichlet_bc.set_bc(entity, BoundaryCondition::Dirichlet(5.0)); - // Create a test matrix and RHS vector let (mut matrix, mut rhs) = create_test_matrix_and_rhs(); let mut matrix_mut = matrix.as_mut(); let mut rhs_mut = rhs.as_mut(); - // Apply the Dirichlet condition dirichlet_bc.apply_bc(&mut matrix_mut, &mut rhs_mut, &entity_to_index, 0.0); - // Check that the row in the matrix corresponding to the entity index is zeroed out for col in 0..matrix_mut.ncols() { if col == 1 { - assert_eq!(matrix_mut[(1, col)], 1.0); // Diagonal should be 1 + assert_eq!(matrix_mut[(1, col)], 1.0); } else { - assert_eq!(matrix_mut[(1, col)], 0.0); // Off-diagonal should be 0 + assert_eq!(matrix_mut[(1, col)], 0.0); } } - - // Check that the RHS has the correct value for the Dirichlet condition assert_eq!(rhs_mut[(1, 0)], 5.0); } @@ -217,7 +144,7 @@ mod tests { fn test_apply_function_based_dirichlet() { let dirichlet_bc = DirichletBC::new(); let entity = MeshEntity::Vertex(2); - let mut entity_to_index = FxHashMap::default(); + let entity_to_index = DashMap::new(); entity_to_index.insert(entity, 2); dirichlet_bc.set_bc( @@ -238,7 +165,6 @@ mod tests { assert_eq!(matrix_mut[(2, col)], 0.0); } } - assert_eq!(rhs_mut[(2, 0)], 7.0); } } diff --git a/src/boundary/neumann.rs b/src/boundary/neumann.rs index 04e3d42..9bc9dbd 100644 --- a/src/boundary/neumann.rs +++ b/src/boundary/neumann.rs @@ -1,80 +1,43 @@ -// src/boundary/neumann.rs - +use dashmap::DashMap; +use std::sync::Arc; use crate::domain::mesh_entity::MeshEntity; -use rustc_hash::FxHashMap; use crate::boundary::bc_handler::{BoundaryCondition, BoundaryConditionApply}; -use crate::domain::section::Section; use faer::MatMut; /// The `NeumannBC` struct represents a handler for applying Neumann boundary conditions /// to a set of mesh entities. Neumann boundary conditions involve specifying the flux across /// a boundary, and they modify only the right-hand side (RHS) of the system without modifying /// the system matrix. -/// -/// Example usage: -/// -/// let neumann_bc = NeumannBC::new(); -/// let entity = MeshEntity::Vertex(1); -/// neumann_bc.set_bc(entity, BoundaryCondition::Neumann(10.0)); -/// neumann_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); -/// pub struct NeumannBC { - conditions: Section, + conditions: DashMap, } impl NeumannBC { - /// Creates a new instance of `NeumannBC` with an empty section to store boundary conditions. - /// - /// Example usage: - /// - /// let neumann_bc = NeumannBC::new(); - /// + /// Creates a new instance of `NeumannBC` with an empty `DashMap` to store boundary conditions. pub fn new() -> Self { Self { - conditions: Section::new(), + conditions: DashMap::new(), } } /// Sets a Neumann boundary condition for a specific mesh entity. - /// - /// # Arguments: - /// * `entity` - The mesh entity to which the boundary condition will be applied. - /// * `condition` - The boundary condition to set (either a constant flux or a functional form). - /// - /// Example usage: - /// - /// let entity = MeshEntity::Vertex(1); - /// neumann_bc.set_bc(entity, BoundaryCondition::Neumann(5.0)); - /// pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) { - self.conditions.set_data(entity, condition); + self.conditions.insert(entity, condition); } /// Applies the stored Neumann boundary conditions to the right-hand side (RHS) of the system. /// It iterates over the stored conditions and applies either constant or function-based Neumann /// boundary conditions to the corresponding entities. - /// - /// # Arguments: - /// * `_matrix` - The mutable system matrix (unused in Neumann BC). - /// * `rhs` - The mutable right-hand side vector. - /// * `entity_to_index` - A hash map that associates mesh entities with their indices in the system. - /// * `time` - The current time, used for time-dependent boundary conditions. - /// - /// Example usage: - /// - /// let entity_to_index = FxHashMap::default(); - /// neumann_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); - /// pub fn apply_bc( &self, _matrix: &mut MatMut, rhs: &mut MatMut, - entity_to_index: &FxHashMap, + entity_to_index: &DashMap, time: f64, ) { - let data = self.conditions.data.read().unwrap(); - for (entity, condition) in data.iter() { - if let Some(&index) = entity_to_index.get(entity) { + self.conditions.iter().for_each(|entry| { + let (entity, condition) = entry.pair(); + if let Some(index) = entity_to_index.get(entity).map(|i| *i) { match condition { BoundaryCondition::Neumann(value) => { self.apply_constant_neumann(rhs, index, *value); @@ -87,38 +50,15 @@ impl NeumannBC { _ => {} } } - } + }); } /// Applies a constant Neumann boundary condition to the right-hand side (RHS) for a specific index. - /// - /// # Arguments: - /// * `rhs` - The mutable right-hand side vector. - /// * `index` - The index of the rhs corresponding to the mesh entity. - /// * `value` - The Neumann flux to be applied. - /// - /// Example usage: - /// - /// neumann_bc.apply_constant_neumann(&mut rhs, 1, 5.0); - /// pub fn apply_constant_neumann(&self, rhs: &mut MatMut, index: usize, value: f64) { rhs.write(index, 0, rhs.read(index, 0) + value); } - /// Retrieves the coordinates of the mesh entity. - /// - /// This method currently returns a default placeholder value, but it can be expanded - /// to extract real entity coordinates if needed. - /// - /// # Arguments: - /// * `_entity` - The mesh entity for which coordinates are being requested. - /// - /// Returns an array of default coordinates `[0.0, 0.0, 0.0]`. - /// - /// Example usage: - /// - /// let coords = neumann_bc.get_coordinates(&entity); - /// + /// Retrieves the coordinates of the mesh entity (currently a placeholder). fn get_coordinates(&self, _entity: &MeshEntity) -> [f64; 3] { [0.0, 0.0, 0.0] } @@ -126,21 +66,14 @@ impl NeumannBC { impl BoundaryConditionApply for NeumannBC { /// Applies the stored Neumann boundary conditions for a specific mesh entity. - /// - /// This implementation utilizes the general `apply_bc` method to modify the rhs. - /// - /// # Arguments: - /// * `_entity` - The mesh entity to which the boundary condition applies. - /// * `rhs` - The mutable right-hand side vector. - /// * `_matrix` - The mutable system matrix (unused in Neumann BC). - /// * `entity_to_index` - A hash map that associates mesh entities with their indices. - /// * `time` - The current time, used for time-dependent boundary conditions. - /// - /// Example usage: - /// - /// neumann_bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.0); - /// - fn apply(&self, _entity: &MeshEntity, rhs: &mut MatMut, _matrix: &mut MatMut, entity_to_index: &FxHashMap, time: f64) { + fn apply( + &self, + _entity: &MeshEntity, + rhs: &mut MatMut, + _matrix: &mut MatMut, + entity_to_index: &DashMap, + time: f64, + ) { self.apply_bc(_matrix, rhs, entity_to_index, time); } } @@ -148,18 +81,13 @@ impl BoundaryConditionApply for NeumannBC { #[cfg(test)] mod tests { use super::*; - use rustc_hash::FxHashMap; use faer::Mat; use crate::domain::mesh_entity::MeshEntity; use std::sync::Arc; fn create_test_matrix_and_rhs() -> (Mat, Mat) { - // Create a 3x3 test matrix initialized to identity matrix (though unused for NeumannBC) let matrix = Mat::from_fn(3, 3, |i, j| if i == j { 1.0 } else { 0.0 }); - - // Create a 3x1 test RHS vector initialized to zero let rhs = Mat::zeros(3, 1); - (matrix, rhs) } @@ -168,11 +96,9 @@ mod tests { let neumann_bc = NeumannBC::new(); let entity = MeshEntity::Vertex(1); - // Set a Neumann boundary condition neumann_bc.set_bc(entity, BoundaryCondition::Neumann(10.0)); - // Verify that the condition was set correctly - let condition = neumann_bc.conditions.restrict(&entity); + let condition = neumann_bc.conditions.get(&entity).map(|entry| entry.clone()); assert!(matches!(condition, Some(BoundaryCondition::Neumann(10.0)))); } @@ -180,20 +106,16 @@ mod tests { fn test_apply_constant_neumann() { let neumann_bc = NeumannBC::new(); let entity = MeshEntity::Vertex(1); - let mut entity_to_index = FxHashMap::default(); + let entity_to_index = DashMap::new(); entity_to_index.insert(entity, 1); - // Set a Neumann boundary condition neumann_bc.set_bc(entity, BoundaryCondition::Neumann(5.0)); - // Create a test matrix and RHS vector let (mut matrix, mut rhs) = create_test_matrix_and_rhs(); let mut rhs_mut = rhs.as_mut(); - // Apply the Neumann condition neumann_bc.apply_bc(&mut matrix.as_mut(), &mut rhs_mut, &entity_to_index, 0.0); - // Check that the RHS has been updated with the Neumann flux assert_eq!(rhs_mut[(1, 0)], 5.0); } @@ -201,7 +123,7 @@ mod tests { fn test_apply_function_based_neumann() { let neumann_bc = NeumannBC::new(); let entity = MeshEntity::Vertex(2); - let mut entity_to_index = FxHashMap::default(); + let entity_to_index = DashMap::new(); entity_to_index.insert(entity, 2); neumann_bc.set_bc( diff --git a/src/boundary/robin.rs b/src/boundary/robin.rs index 174d3dd..fcc1e0b 100644 --- a/src/boundary/robin.rs +++ b/src/boundary/robin.rs @@ -1,80 +1,42 @@ -// src/boundary/robin.rs - +use dashmap::DashMap; use crate::domain::mesh_entity::MeshEntity; -use rustc_hash::FxHashMap; use crate::boundary::bc_handler::{BoundaryCondition, BoundaryConditionApply}; -use crate::domain::section::Section; use faer::MatMut; /// The `RobinBC` struct represents a handler for applying Robin boundary conditions /// to a set of mesh entities. Robin boundary conditions involve a linear combination /// of Dirichlet and Neumann boundary conditions, and they modify both the system matrix /// and the right-hand side (RHS). -/// -/// Example usage: -/// -/// let robin_bc = RobinBC::new(); -/// let entity = MeshEntity::Vertex(1); -/// robin_bc.set_bc(entity, BoundaryCondition::Robin { alpha: 2.0, beta: 3.0 }); -/// robin_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); -/// pub struct RobinBC { - conditions: Section, + conditions: DashMap, } impl RobinBC { - /// Creates a new instance of `RobinBC` with an empty section to store boundary conditions. - /// - /// Example usage: - /// - /// let robin_bc = RobinBC::new(); - /// + /// Creates a new instance of `RobinBC` with an empty `DashMap` to store boundary conditions. pub fn new() -> Self { Self { - conditions: Section::new(), + conditions: DashMap::new(), } } /// Sets a Robin boundary condition for a specific mesh entity. - /// - /// # Arguments: - /// * `entity` - The mesh entity to which the boundary condition will be applied. - /// * `condition` - The boundary condition to set (Robin condition with alpha and beta). - /// - /// Example usage: - /// - /// let entity = MeshEntity::Vertex(1); - /// robin_bc.set_bc(entity, BoundaryCondition::Robin { alpha: 2.0, beta: 3.0 }); - /// pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) { - self.conditions.set_data(entity, condition); + self.conditions.insert(entity, condition); } /// Applies the stored Robin boundary conditions to both the system matrix and rhs. /// It iterates over the stored conditions and applies the Robin boundary condition /// to the corresponding entities. - /// - /// # Arguments: - /// * `matrix` - The mutable system matrix. - /// * `rhs` - The mutable right-hand side vector. - /// * `entity_to_index` - A hash map that associates mesh entities with their indices in the system. - /// * `_time` - The current time (unused in Robin BC). - /// - /// Example usage: - /// - /// let entity_to_index = FxHashMap::default(); - /// robin_bc.apply_bc(&mut matrix, &mut rhs, &entity_to_index, 0.0); - /// pub fn apply_bc( &self, matrix: &mut MatMut, rhs: &mut MatMut, - entity_to_index: &FxHashMap, + entity_to_index: &DashMap, _time: f64, ) { - let data = self.conditions.data.read().unwrap(); - for (entity, condition) in data.iter() { - if let Some(&index) = entity_to_index.get(entity) { + self.conditions.iter().for_each(|entry| { + let (entity, condition) = entry.pair(); + if let Some(index) = entity_to_index.get(entity).map(|i| *i) { match condition { BoundaryCondition::Robin { alpha, beta } => { self.apply_robin(matrix, rhs, index, *alpha, *beta); @@ -82,22 +44,10 @@ impl RobinBC { _ => {} } } - } + }); } /// Applies a Robin boundary condition to the system matrix and rhs for a specific index. - /// - /// # Arguments: - /// * `matrix` - The mutable system matrix. - /// * `rhs` - The mutable right-hand side vector. - /// * `index` - The index of the matrix row and rhs corresponding to the mesh entity. - /// * `alpha` - The coefficient for the Robin condition in the matrix (modifying the diagonal). - /// * `beta` - The constant value for the Robin condition (added to rhs). - /// - /// Example usage: - /// - /// robin_bc.apply_robin(&mut matrix, &mut rhs, 1, 2.0, 3.0); - /// pub fn apply_robin( &self, matrix: &mut MatMut, @@ -113,21 +63,14 @@ impl RobinBC { impl BoundaryConditionApply for RobinBC { /// Applies the stored Robin boundary conditions for a specific mesh entity. - /// - /// This implementation utilizes the general `apply_bc` method to modify the matrix and rhs. - /// - /// # Arguments: - /// * `_entity` - The mesh entity to which the boundary condition applies. - /// * `rhs` - The mutable right-hand side vector. - /// * `matrix` - The mutable system matrix. - /// * `entity_to_index` - A hash map that associates mesh entities with their indices. - /// * `time` - The current time (unused in Robin BC). - /// - /// Example usage: - /// - /// robin_bc.apply(&entity, &mut rhs, &mut matrix, &entity_to_index, 0.0); - /// - fn apply(&self, _entity: &MeshEntity, rhs: &mut MatMut, matrix: &mut MatMut, entity_to_index: &FxHashMap, time: f64) { + fn apply( + &self, + _entity: &MeshEntity, + rhs: &mut MatMut, + matrix: &mut MatMut, + entity_to_index: &DashMap, + time: f64, + ) { self.apply_bc(matrix, rhs, entity_to_index, time); } } @@ -135,17 +78,12 @@ impl BoundaryConditionApply for RobinBC { #[cfg(test)] mod tests { use super::*; - use rustc_hash::FxHashMap; use faer::Mat; use crate::domain::mesh_entity::MeshEntity; fn create_test_matrix_and_rhs() -> (Mat, Mat) { - // Create a 3x3 test matrix initialized to identity matrix let matrix = Mat::from_fn(3, 3, |i, j| if i == j { 1.0 } else { 0.0 }); - - // Create a 3x1 test RHS vector initialized to zero let rhs = Mat::zeros(3, 1); - (matrix, rhs) } @@ -154,11 +92,9 @@ mod tests { let robin_bc = RobinBC::new(); let entity = MeshEntity::Vertex(1); - // Set a Robin boundary condition robin_bc.set_bc(entity, BoundaryCondition::Robin { alpha: 2.0, beta: 3.0 }); - // Verify that the condition was set correctly - let condition = robin_bc.conditions.restrict(&entity); + let condition = robin_bc.conditions.get(&entity).map(|entry| entry.clone()); assert!(matches!(condition, Some(BoundaryCondition::Robin { alpha: 2.0, beta: 3.0 }))); } @@ -166,24 +102,18 @@ mod tests { fn test_apply_robin_bc() { let robin_bc = RobinBC::new(); let entity = MeshEntity::Vertex(1); - let mut entity_to_index = FxHashMap::default(); + let entity_to_index = DashMap::new(); entity_to_index.insert(entity, 1); - // Set a Robin boundary condition robin_bc.set_bc(entity, BoundaryCondition::Robin { alpha: 2.0, beta: 3.0 }); - // Create a test matrix and RHS vector let (mut matrix, mut rhs) = create_test_matrix_and_rhs(); let mut matrix_mut = matrix.as_mut(); let mut rhs_mut = rhs.as_mut(); - // Apply the Robin condition robin_bc.apply_bc(&mut matrix_mut, &mut rhs_mut, &entity_to_index, 0.0); - // Check that the matrix diagonal has been updated (alpha term) assert_eq!(matrix_mut[(1, 1)], 3.0); // Initial value 1.0 + alpha 2.0 - - // Check that the RHS has been updated (beta term) - assert_eq!(rhs_mut[(1, 0)], 3.0); // Beta term applied + assert_eq!(rhs_mut[(1, 0)], 3.0); // Beta term applied } } diff --git a/src/domain/entity_fill.rs b/src/domain/entity_fill.rs index 1563cb9..851173b 100644 --- a/src/domain/entity_fill.rs +++ b/src/domain/entity_fill.rs @@ -1,6 +1,6 @@ use crate::domain::mesh_entity::MeshEntity; use crate::domain::sieve::Sieve; -use rustc_hash::FxHashSet; +use dashmap::DashMap; impl Sieve { /// Infers and adds missing edges (in 2D) or faces (in 3D) based on existing cells and vertices. @@ -13,32 +13,31 @@ impl Sieve { /// sieve.fill_missing_entities(); /// pub fn fill_missing_entities(&self) { - let mut edge_set: FxHashSet<(MeshEntity, MeshEntity)> = FxHashSet::default(); - - // Acquire a read lock to access the adjacency data. - let adjacency = self.adjacency.read().unwrap(); + // Use DashMap instead of FxHashSet for concurrent access. + let edge_set: DashMap<(MeshEntity, MeshEntity), ()> = DashMap::new(); // Loop through each cell and infer its edges (for 2D meshes) - for (cell, vertices) in adjacency.iter() { + self.adjacency.iter().for_each(|entry| { + let cell = entry.key(); if let MeshEntity::Cell(_) = cell { - let vertices: Vec<_> = vertices.iter().collect(); + let vertices: Vec<_> = entry.value().iter().map(|v| v.key().clone()).collect(); // Connect each vertex with its neighboring vertex to form edges. for i in 0..vertices.len() { - let v1 = vertices[i]; - let v2 = vertices[(i + 1) % vertices.len()]; - let edge = if v1 < v2 { (*v1, *v2) } else { (*v2, *v1) }; - edge_set.insert(edge); + let v1 = vertices[i].clone(); + let v2 = vertices[(i + 1) % vertices.len()].clone(); + let edge = if v1 < v2 { (v1, v2) } else { (v2, v1) }; + edge_set.insert(edge, ()); } } - } + }); // Add the deduced edges to the sieve. - let adjacency = self.adjacency.write().unwrap(); - for (v1, v2) in edge_set { + let edge_count = self.adjacency.len(); + edge_set.into_iter().enumerate().for_each(|(index, ((v1, v2), _))| { // Generate a unique ID for the new edge. - let edge = MeshEntity::Edge(adjacency.len()); - self.add_arrow(v1, edge); + let edge = MeshEntity::Edge(edge_count + index); + self.add_arrow(v1, edge.clone()); self.add_arrow(v2, edge); - } + }); } } diff --git a/src/domain/mesh/geometry.rs b/src/domain/mesh/geometry.rs index 6a3d2d5..3110db6 100644 --- a/src/domain/mesh/geometry.rs +++ b/src/domain/mesh/geometry.rs @@ -1,7 +1,7 @@ use super::Mesh; use crate::domain::mesh_entity::MeshEntity; use crate::geometry::{Geometry, CellShape, FaceShape}; -use rustc_hash::FxHashSet; +use dashmap::DashMap; impl Mesh { /// Retrieves all the faces of a given cell. @@ -12,12 +12,12 @@ impl Mesh { /// Returns a set of `MeshEntity` representing the faces of the cell, or /// `None` if the cell has no connected faces. /// - /// Example usage: - /// - /// let faces = mesh.get_faces_of_cell(&cell); - /// - pub fn get_faces_of_cell(&self, cell: &MeshEntity) -> Option> { - self.sieve.cone(cell).map(|set| set.clone()) + pub fn get_faces_of_cell(&self, cell: &MeshEntity) -> Option> { + self.sieve.cone(cell).map(|set| { + let faces = DashMap::new(); + set.into_iter().for_each(|face| { faces.insert(face, ()); }); + faces + }) } /// Retrieves all the cells that share the given face. @@ -27,12 +27,10 @@ impl Mesh { /// /// Returns a set of `MeshEntity` representing the neighboring cells. /// - /// Example usage: - /// - /// let cells = mesh.get_cells_sharing_face(&face); - /// - pub fn get_cells_sharing_face(&self, face: &MeshEntity) -> FxHashSet { - self.sieve.support(face) + pub fn get_cells_sharing_face(&self, face: &MeshEntity) -> DashMap { + let cells = DashMap::new(); + self.sieve.support(face).into_iter().for_each(|cell| { cells.insert(cell, ()); }); + cells } /// Computes the Euclidean distance between two cells based on their centroids. @@ -40,10 +38,6 @@ impl Mesh { /// This method calculates the centroids of both cells and then uses the `Geometry` /// module to compute the distance between these centroids. /// - /// Example usage: - /// - /// let distance = mesh.get_distance_between_cells(&cell1, &cell2); - /// pub fn get_distance_between_cells(&self, cell_i: &MeshEntity, cell_j: &MeshEntity) -> f64 { let centroid_i = self.get_cell_centroid(cell_i); let centroid_j = self.get_cell_centroid(cell_j); @@ -55,12 +49,6 @@ impl Mesh { /// This method determines the face shape (triangle or quadrilateral) and /// uses the `Geometry` module to compute the area. /// - /// Panics if the face has an unsupported number of vertices. - /// - /// Example usage: - /// - /// let area = mesh.get_face_area(&face); - /// pub fn get_face_area(&self, face: &MeshEntity) -> f64 { let face_vertices = self.get_face_vertices(face); let face_shape = match face_vertices.len() { @@ -69,22 +57,14 @@ impl Mesh { _ => panic!("Unsupported face shape with {} vertices", face_vertices.len()), }; - // Use a new or existing geometry instance, providing a unique ID for caching. let mut geometry = Geometry::new(); - let face_id = face.id(); // Assuming `id` method provides a unique ID for the face. + let face_id = face.id(); geometry.compute_face_area(face_id, face_shape, &face_vertices) } /// Computes the centroid of a cell based on its vertices. /// - /// This method determines the cell shape (tetrahedron, pyramid, prism, or hexahedron) - /// and uses the `Geometry` module to compute the centroid. - /// - /// Panics if the cell has an unsupported number of vertices. - /// - /// Example usage: - /// - /// let centroid = mesh.get_cell_centroid(&cell); + /// This method determines the cell shape and uses the `Geometry` module to compute the centroid. /// pub fn get_cell_centroid(&self, cell: &MeshEntity) -> [f64; 3] { let cell_vertices = self.get_cell_vertices(cell); @@ -96,7 +76,6 @@ impl Mesh { _ => panic!("Unsupported cell shape with {} vertices", cell_vertices.len()), }; - // Use a new or existing Geometry instance, providing a unique ID for caching. let mut geometry = Geometry::new(); geometry.compute_cell_centroid(self, cell) } @@ -106,49 +85,32 @@ impl Mesh { /// This method uses the `support` function of the sieve to find cells that /// contain the given vertex and then retrieves all other vertices in those cells. /// - /// Example usage: - /// - /// let neighbors = mesh.get_neighboring_vertices(&vertex); - /// pub fn get_neighboring_vertices(&self, vertex: &MeshEntity) -> Vec { - let mut neighbors = FxHashSet::default(); + let neighbors = DashMap::new(); let connected_cells = self.sieve.support(vertex); - for cell in &connected_cells { - if let Some(cell_vertices) = self.sieve.cone(cell).as_ref() { + connected_cells.into_iter().for_each(|cell| { + if let Some(cell_vertices) = self.sieve.cone(&cell).as_ref() { for v in cell_vertices { if v != vertex && matches!(v, MeshEntity::Vertex(_)) { - neighbors.insert(*v); + neighbors.insert(v.clone(), ()); } } } else { panic!("Cell {:?} has no connected vertices", cell); } - } - neighbors.into_iter().collect() + }); + neighbors.into_iter().map(|(vertex, _)| vertex).collect() } /// Returns an iterator over the IDs of all vertices in the mesh. /// - /// Example usage: - /// - /// for vertex_id in mesh.iter_vertices() { - /// println!("Vertex ID: {}", vertex_id); - /// } - /// pub fn iter_vertices(&self) -> impl Iterator { self.vertex_coordinates.keys() } /// Determines the shape of a cell based on the number of vertices it has. /// - /// This method supports tetrahedrons, pyramids, prisms, and hexahedrons. - /// Returns an error if the cell has an unsupported number of vertices. - /// - /// Example usage: - /// - /// let cell_shape = mesh.get_cell_shape(&cell).unwrap(); - /// pub fn get_cell_shape(&self, cell: &MeshEntity) -> Result { let cell_vertices = self.get_cell_vertices(cell); match cell_vertices.len() { @@ -165,22 +127,11 @@ impl Mesh { /// Retrieves the vertices of a cell and their coordinates. /// - /// This method uses the sieve to find all vertices connected to the cell, - /// and retrieves their 3D coordinates from the `vertex_coordinates` map. - /// - /// Panics if the coordinates for a vertex are not found. - /// - /// Example usage: - /// - /// let vertices = mesh.get_cell_vertices(&cell); - /// pub fn get_cell_vertices(&self, cell: &MeshEntity) -> Vec<[f64; 3]> { let mut vertices = Vec::new(); - // Retrieve entities connected to the cell. if let Some(connected_entities) = self.sieve.cone(cell) { for entity in connected_entities { - // Check if the entity is a vertex and extract its coordinates. if let MeshEntity::Vertex(vertex_id) = entity { if let Some(coords) = self.get_vertex_coordinates(vertex_id) { vertices.push(coords); @@ -190,8 +141,7 @@ impl Mesh { } } } - - // Sort and deduplicate vertices to ensure unique entries. + vertices.sort_by(|a, b| a.partial_cmp(b).unwrap()); vertices.dedup(); vertices @@ -199,13 +149,6 @@ impl Mesh { /// Retrieves the vertices of a face and their coordinates. /// - /// This method uses the sieve to find all vertices connected to the face, - /// and retrieves their 3D coordinates from the `vertex_coordinates` map. - /// - /// Example usage: - /// - /// let face_vertices = mesh.get_face_vertices(&face); - /// pub fn get_face_vertices(&self, face: &MeshEntity) -> Vec<[f64; 3]> { let mut vertices = Vec::new(); if let Some(connected_vertices) = self.sieve.cone(face) { diff --git a/src/domain/overlap.rs b/src/domain/overlap.rs index 9840409..037ee8e 100644 --- a/src/domain/overlap.rs +++ b/src/domain/overlap.rs @@ -1,241 +1,124 @@ -use rustc_hash::{FxHashMap, FxHashSet}; -use std::sync::{Arc, RwLock}; +use dashmap::DashMap; +use std::sync::Arc; use crate::domain::mesh_entity::MeshEntity; /// The `Overlap` struct manages two sets of `MeshEntity` elements: /// - `local_entities`: Entities that are local to the current partition. -/// - `ghost_entities`: Entities that are shared with other partitions. -/// -/// It supports adding entities to these sets, checking whether an entity -/// is local or ghost, and merging overlaps from different partitions. -/// -/// Example usage: -/// -/// let overlap = Overlap::new(); -/// overlap.add_local_entity(MeshEntity::Vertex(1)); -/// assert!(overlap.is_local(&MeshEntity::Vertex(1))); -/// +/// - `ghost_entities`: Entities that are shared with other partitions. pub struct Overlap { /// A thread-safe set of local entities. - pub local_entities: Arc>>, + pub local_entities: Arc>, /// A thread-safe set of ghost entities. - pub ghost_entities: Arc>>, + pub ghost_entities: Arc>, } impl Overlap { - /// Creates a new `Overlap` with empty sets for local and ghost entities. - /// - /// Example usage: - /// - /// let overlap = Overlap::new(); - /// + /// Creates a new `Overlap` with empty sets for local and ghost entities. pub fn new() -> Self { Overlap { - local_entities: Arc::new(RwLock::new(FxHashSet::default())), - ghost_entities: Arc::new(RwLock::new(FxHashSet::default())), + local_entities: Arc::new(DashMap::new()), + ghost_entities: Arc::new(DashMap::new()), } } - /// Adds a `MeshEntity` to the set of local entities. - /// - /// Example usage: - /// - /// overlap.add_local_entity(MeshEntity::Vertex(1)); - /// + /// Adds a `MeshEntity` to the set of local entities. pub fn add_local_entity(&self, entity: MeshEntity) { - self.local_entities.write().unwrap().insert(entity); + self.local_entities.insert(entity, ()); } - /// Adds a `MeshEntity` to the set of ghost entities. - /// - /// Example usage: - /// - /// overlap.add_ghost_entity(MeshEntity::Vertex(2)); - /// + /// Adds a `MeshEntity` to the set of ghost entities. pub fn add_ghost_entity(&self, entity: MeshEntity) { - self.ghost_entities.write().unwrap().insert(entity); + self.ghost_entities.insert(entity, ()); } - /// Checks if a `MeshEntity` is a local entity. - /// - /// Returns `true` if the entity is in the local entities set, otherwise `false`. - /// - /// Example usage: - /// - /// let is_local = overlap.is_local(&MeshEntity::Vertex(1)); - /// + /// Checks if a `MeshEntity` is a local entity. pub fn is_local(&self, entity: &MeshEntity) -> bool { - let local = self.local_entities.read().unwrap(); - local.contains(entity) + self.local_entities.contains_key(entity) } - /// Checks if a `MeshEntity` is a ghost entity. - /// - /// Returns `true` if the entity is in the ghost entities set, otherwise `false`. - /// - /// Example usage: - /// - /// let is_ghost = overlap.is_ghost(&MeshEntity::Vertex(2)); - /// + /// Checks if a `MeshEntity` is a ghost entity. pub fn is_ghost(&self, entity: &MeshEntity) -> bool { - let ghost = self.ghost_entities.read().unwrap(); - ghost.contains(entity) + self.ghost_entities.contains_key(entity) } - /// Retrieves a clone of all local entities. - /// - /// Example usage: - /// - /// let local_entities = overlap.local_entities(); - /// - pub fn local_entities(&self) -> FxHashSet { - let local = self.local_entities.read().unwrap(); - local.clone() + /// Retrieves a clone of all local entities. + pub fn local_entities(&self) -> Vec { + self.local_entities.iter().map(|entry| entry.key().clone()).collect() } - /// Retrieves a clone of all ghost entities. - /// - /// Example usage: - /// - /// let ghost_entities = overlap.ghost_entities(); - /// - pub fn ghost_entities(&self) -> FxHashSet { - let ghost = self.ghost_entities.read().unwrap(); - ghost.clone() + /// Retrieves a clone of all ghost entities. + pub fn ghost_entities(&self) -> Vec { + self.ghost_entities.iter().map(|entry| entry.key().clone()).collect() } /// Merges another `Overlap` instance into this one, combining local - /// and ghost entities from both overlaps. - /// - /// Example usage: - /// - /// overlap1.merge(&overlap2); - /// + /// and ghost entities from both overlaps. pub fn merge(&self, other: &Overlap) { - let mut local = self.local_entities.write().unwrap(); - let other_local = other.local_entities.read().unwrap(); - local.extend(other_local.iter().cloned()); + other.local_entities.iter().for_each(|entry| { + self.local_entities.insert(entry.key().clone(), ()); + }); - let mut ghost = self.ghost_entities.write().unwrap(); - let other_ghost = other.ghost_entities.read().unwrap(); - ghost.extend(other_ghost.iter().cloned()); + other.ghost_entities.iter().for_each(|entry| { + self.ghost_entities.insert(entry.key().clone(), ()); + }); } } /// The `Delta` struct manages transformation data for `MeshEntity` elements /// in overlapping regions. It is used to store and apply data transformations -/// across entities in distributed environments. -/// -/// Example usage: -/// -/// let delta = Delta::new(); -/// delta.set_data(MeshEntity::Vertex(1), 42); -/// assert_eq!(delta.get_data(&MeshEntity::Vertex(1)), Some(42)); -/// +/// across entities in distributed environments. pub struct Delta { /// A thread-safe map storing transformation data associated with `MeshEntity` objects. - pub data: Arc>>, // Transformation data over overlapping regions + pub data: Arc>, // Transformation data over overlapping regions } impl Delta { - /// Creates a new, empty `Delta`. - /// - /// Example usage: - /// - /// let delta = Delta::new(); - /// + /// Creates a new, empty `Delta`. pub fn new() -> Self { Delta { - data: Arc::new(RwLock::new(FxHashMap::default())), + data: Arc::new(DashMap::new()), } } - /// Sets the transformation data for a specific `MeshEntity`. - /// - /// Example usage: - /// - /// delta.set_data(MeshEntity::Vertex(1), 42); - /// + /// Sets the transformation data for a specific `MeshEntity`. pub fn set_data(&self, entity: MeshEntity, value: T) { - let mut data = self.data.write().unwrap(); - data.insert(entity, value); + self.data.insert(entity, value); } - /// Retrieves the transformation data associated with a specific `MeshEntity`. - /// - /// Returns `None` if no data is found for the entity. - /// - /// Example usage: - /// - /// let value = delta.get_data(&MeshEntity::Vertex(1)); - /// + /// Retrieves the transformation data associated with a specific `MeshEntity`. pub fn get_data(&self, entity: &MeshEntity) -> Option where T: Clone, { - let data = self.data.read().unwrap(); - data.get(entity).cloned() + self.data.get(entity).map(|entry| entry.clone()) } - /// Removes the transformation data associated with a specific `MeshEntity`. - /// - /// Returns the removed data if it exists, otherwise `None`. - /// - /// Example usage: - /// - /// let removed_value = delta.remove_data(&MeshEntity::Vertex(1)); - /// + /// Removes the transformation data associated with a specific `MeshEntity`. pub fn remove_data(&self, entity: &MeshEntity) -> Option { - let mut data = self.data.write().unwrap(); - data.remove(entity) + self.data.remove(entity).map(|(_, value)| value) } - /// Checks if there is transformation data for a specific `MeshEntity`. - /// - /// Returns `true` if the entity has associated data, otherwise `false`. - /// - /// Example usage: - /// - /// let has_data = delta.has_data(&MeshEntity::Vertex(1)); - /// + /// Checks if there is transformation data for a specific `MeshEntity`. pub fn has_data(&self, entity: &MeshEntity) -> bool { - let data = self.data.read().unwrap(); - data.contains_key(entity) + self.data.contains_key(entity) } - /// Applies a function to all entities in the delta. - /// - /// Example usage: - /// - /// delta.apply(|entity, value| { - /// println!("{:?}: {:?}", entity, value); - /// }); - /// + /// Applies a function to all entities in the delta. pub fn apply(&self, mut func: F) where F: FnMut(&MeshEntity, &T), { - let data = self.data.read().unwrap(); - for (entity, value) in data.iter() { - func(entity, value); - } + self.data.iter().for_each(|entry| func(entry.key(), entry.value())); } - /// Merges another `Delta` instance into this one, combining data from both deltas. - /// - /// Example usage: - /// - /// delta1.merge(&delta2); - /// + /// Merges another `Delta` instance into this one, combining data from both deltas. pub fn merge(&self, other: &Delta) where T: Clone, { - let mut data = self.data.write().unwrap(); - let other_data = other.data.read().unwrap(); - for (entity, value) in other_data.iter() { - data.insert(entity.clone(), value.clone()); - } + other.data.iter().for_each(|entry| { + self.data.insert(entry.key().clone(), entry.value().clone()); + }); } } @@ -245,8 +128,6 @@ mod tests { use crate::domain::mesh_entity::MeshEntity; #[test] - /// Test that verifies the addition and retrieval of local and ghost entities - /// in the `Overlap` structure. fn test_overlap_local_and_ghost_entities() { let overlap = Overlap::new(); let vertex_local = MeshEntity::Vertex(1); @@ -258,7 +139,6 @@ mod tests { } #[test] - /// Test that verifies merging two `Overlap` structures works as expected. fn test_overlap_merge() { let overlap1 = Overlap::new(); let overlap2 = Overlap::new(); @@ -280,8 +160,6 @@ mod tests { } #[test] - /// Test that verifies setting and retrieving data in the `Delta` structure - /// works as expected. fn test_delta_set_and_get_data() { let delta = Delta::new(); let vertex = MeshEntity::Vertex(1); @@ -293,7 +171,6 @@ mod tests { } #[test] - /// Test that verifies removing data from the `Delta` structure works as expected. fn test_delta_remove_data() { let delta = Delta::new(); let vertex = MeshEntity::Vertex(1); @@ -304,7 +181,6 @@ mod tests { } #[test] - /// Test that verifies merging two `Delta` structures works as expected. fn test_delta_merge() { let delta1 = Delta::new(); let delta2 = Delta::new(); diff --git a/src/domain/section.rs b/src/domain/section.rs index 1666d33..1d05666 100644 --- a/src/domain/section.rs +++ b/src/domain/section.rs @@ -1,14 +1,13 @@ -use rustc_hash::FxHashMap; -use std::sync::{Arc, RwLock}; +use dashmap::DashMap; use rayon::prelude::*; use crate::domain::mesh_entity::MeshEntity; /// A generic `Section` struct that associates data of type `T` with `MeshEntity` elements. /// It provides methods for setting, updating, and retrieving data, and supports /// parallel updates for performance improvements. -/// +/// /// Example usage: -/// +/// /// let section = Section::new(); /// let vertex = MeshEntity::Vertex(1); /// section.set_data(vertex, 42); @@ -16,20 +15,20 @@ use crate::domain::mesh_entity::MeshEntity; /// pub struct Section { /// A thread-safe map storing data of type `T` associated with `MeshEntity` objects. - pub data: Arc>>, + pub data: DashMap, } impl Section { /// Creates a new `Section` with an empty data map. /// /// Example usage: - /// + /// /// let section = Section::new(); - /// assert!(section.data.read().unwrap().is_empty()); + /// assert!(section.data.is_empty()); /// pub fn new() -> Self { Section { - data: Arc::new(RwLock::new(FxHashMap::default())), + data: DashMap::new(), } } @@ -37,13 +36,12 @@ impl Section { /// This method inserts the `entity` and its corresponding `value` into the data map. /// /// Example usage: - /// + /// /// let section = Section::new(); /// section.set_data(MeshEntity::Vertex(1), 10); /// pub fn set_data(&self, entity: MeshEntity, value: T) { - let mut data = self.data.write().unwrap(); - data.insert(entity, value); + self.data.insert(entity, value); } /// Restricts the data for a given `MeshEntity` by returning an immutable copy of the data @@ -52,7 +50,7 @@ impl Section { /// Returns `None` if no data is found for the entity. /// /// Example usage: - /// + /// /// let section = Section::new(); /// let vertex = MeshEntity::Vertex(1); /// section.set_data(vertex, 42); @@ -62,14 +60,13 @@ impl Section { where T: Clone, { - self.data.read().unwrap().get(entity).cloned() + self.data.get(entity).map(|v| v.clone()) } - /// Applies the given function in parallel to update all data values in the section. - /// This method uses Rayon to iterate over the data map concurrently for performance. + /// Applies the given function in parallel to update all data values in the section. /// /// Example usage: - /// + /// /// section.parallel_update(|v| *v += 1); /// pub fn parallel_update(&self, update_fn: F) @@ -77,8 +74,15 @@ impl Section { F: Fn(&mut T) + Sync + Send, T: Send + Sync, { - let mut data = self.data.write().unwrap(); - data.par_iter_mut().for_each(|(_, v)| update_fn(v)); + // Clone the keys to ensure safe access to each mutable entry in parallel. + let keys: Vec = self.data.iter().map(|entry| entry.key().clone()).collect(); + + // Apply the update function to each entry in parallel. + keys.into_par_iter().for_each(|key| { + if let Some(mut entry) = self.data.get_mut(&key) { + update_fn(entry.value_mut()); + } + }); } /// Restricts the data for a given `MeshEntity` by returning a mutable copy of the data @@ -87,7 +91,7 @@ impl Section { /// Returns `None` if no data is found for the entity. /// /// Example usage: - /// + /// /// let section = Section::new(); /// let vertex = MeshEntity::Vertex(1); /// section.set_data(vertex, 5); @@ -99,32 +103,29 @@ impl Section { where T: Clone, { - let data = self.data.write().unwrap(); - data.get(entity).cloned() + self.data.get(entity).map(|v| v.clone()) } /// Updates the data for a specific `MeshEntity` by replacing the existing value /// with the new value. /// /// Example usage: - /// + /// /// section.update_data(&MeshEntity::Vertex(1), 15); /// pub fn update_data(&self, entity: &MeshEntity, new_value: T) { - let mut data = self.data.write().unwrap(); - data.insert(*entity, new_value); + self.data.insert(*entity, new_value); } /// Clears all data from the section, removing all entity associations. /// /// Example usage: - /// + /// /// section.clear(); - /// assert!(section.data.read().unwrap().is_empty()); + /// assert!(section.data.is_empty()); /// pub fn clear(&self) { - let mut data = self.data.write().unwrap(); - data.clear(); + self.data.clear(); } /// Retrieves all `MeshEntity` objects associated with the section. @@ -132,12 +133,11 @@ impl Section { /// Returns a vector containing all mesh entities currently stored in the section. /// /// Example usage: - /// + /// /// let entities = section.entities(); /// pub fn entities(&self) -> Vec { - let data = self.data.read().unwrap(); - data.keys().cloned().collect() + self.data.iter().map(|entry| entry.key().clone()).collect() } /// Retrieves all data stored in the section as immutable copies. @@ -145,15 +145,14 @@ impl Section { /// Returns a vector of data values. /// /// Example usage: - /// + /// /// let all_data = section.all_data(); /// pub fn all_data(&self) -> Vec where T: Clone, { - let data = self.data.read().unwrap(); - data.values().cloned().collect() + self.data.iter().map(|entry| entry.value().clone()).collect() } /// Retrieves all data stored in the section with mutable access. @@ -161,18 +160,18 @@ impl Section { /// Returns a vector of data values that can be modified. /// /// Example usage: - /// + /// /// let all_data_mut = section.all_data_mut(); /// pub fn all_data_mut(&self) -> Vec where T: Clone, { - let data = self.data.write().unwrap(); - data.values().cloned().collect() + self.data.iter_mut().map(|entry| entry.value().clone()).collect() } } + #[cfg(test)] mod tests { use super::*; diff --git a/src/domain/sieve.rs b/src/domain/sieve.rs index c472f92..281453d 100644 --- a/src/domain/sieve.rs +++ b/src/domain/sieve.rs @@ -1,211 +1,151 @@ -use rustc_hash::{FxHashMap, FxHashSet}; -use std::sync::{Arc, RwLock}; +use dashmap::DashMap; use rayon::prelude::*; use crate::domain::mesh_entity::MeshEntity; /// A `Sieve` struct that manages the relationships (arrows) between `MeshEntity` -/// elements, organized in an adjacency map. +/// elements, organized in an adjacency map. /// /// The adjacency map tracks directed relations between entities in the mesh. /// It supports operations such as adding relationships, querying direct -/// relations (cones), and computing closure and star sets for entities. -/// -/// Example usage: -/// -/// let sieve = Sieve::new(); -/// let vertex = MeshEntity::Vertex(1); -/// let edge = MeshEntity::Edge(2); -/// sieve.add_arrow(vertex, edge); -/// let cone = sieve.cone(&vertex).unwrap(); -/// +/// relations (cones), and computing closure and star sets for entities. #[derive(Clone, Debug)] pub struct Sieve { /// A thread-safe adjacency map where each key is a `MeshEntity`, /// and the value is a set of `MeshEntity` objects related to the key. - pub adjacency: Arc>>>, + pub adjacency: DashMap>, } impl Sieve { - /// Creates a new empty `Sieve` instance with an empty adjacency map. - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// assert!(sieve.adjacency.read().unwrap().is_empty()); - /// + /// Creates a new empty `Sieve` instance with an empty adjacency map. pub fn new() -> Self { Sieve { - adjacency: Arc::new(RwLock::new(FxHashMap::default())), + adjacency: DashMap::new(), } } /// Adds a directed relationship (arrow) between two `MeshEntity` elements. /// The relationship is stored in the adjacency map from the `from` entity - /// to the `to` entity. - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex = MeshEntity::Vertex(1); - /// let edge = MeshEntity::Edge(2); - /// sieve.add_arrow(vertex, edge); - /// + /// to the `to` entity. pub fn add_arrow(&self, from: MeshEntity, to: MeshEntity) { - let mut adjacency = self.adjacency.write().unwrap(); - adjacency.entry(from).or_insert_with(FxHashSet::default).insert(to); + self.adjacency + .entry(from) + .or_insert_with(DashMap::new) + .insert(to, ()); } /// Retrieves all entities directly related to the given entity (`point`). /// This operation is referred to as retrieving the cone of the entity. - /// Returns `None` if there are no related entities. - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex = MeshEntity::Vertex(1); - /// sieve.add_arrow(vertex, MeshEntity::Edge(1)); - /// let cone = sieve.cone(&vertex); - /// - pub fn cone(&self, point: &MeshEntity) -> Option> { - self.adjacency.read().unwrap().get(point).cloned() + /// Returns `None` if there are no related entities. + pub fn cone(&self, point: &MeshEntity) -> Option> { + self.adjacency.get(point).map(|cone| { + cone.iter().map(|entry| entry.key().clone()).collect() + }) } - /// Computes the transitive closure for a given `MeshEntity`. - /// The closure includes the entity itself and all entities reachable - /// through arrows from the entity. - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex = MeshEntity::Vertex(1); - /// sieve.add_arrow(vertex, MeshEntity::Edge(1)); - /// let closure_set = sieve.closure(&vertex); - /// - pub fn closure(&self, point: &MeshEntity) -> FxHashSet { - let mut result = FxHashSet::default(); - let mut stack = vec![point.clone()]; - while let Some(p) = stack.pop() { - if let Some(cones) = self.cone(&p) { - for q in cones { - if result.insert(q.clone()) { - stack.push(q.clone()); + /// Computes the closure of a given `MeshEntity`. + /// The closure includes the entity itself and all entities it covers (cones) recursively. + pub fn closure(&self, point: &MeshEntity) -> DashMap { + let result = DashMap::new(); + let stack = DashMap::new(); + stack.insert(point.clone(), ()); + + while !stack.is_empty() { + let keys: Vec = stack.iter().map(|entry| entry.key().clone()).collect(); + for p in keys { + if result.insert(p.clone(), ()).is_none() { + if let Some(cones) = self.cone(&p) { + for q in cones { + stack.insert(q, ()); + } } } + stack.remove(&p); } } result } /// Computes the star of a given `MeshEntity`. - /// The star includes all entities directly related to the entity - /// (cone and support sets). - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex = MeshEntity::Vertex(1); - /// let edge = MeshEntity::Edge(1); - /// sieve.add_arrow(vertex, edge); - /// let star_set = sieve.star(&vertex); - /// - pub fn star(&self, point: &MeshEntity) -> FxHashSet { - let mut result = FxHashSet::default(); - let mut stack = vec![point.clone()]; - - while let Some(p) = stack.pop() { - if result.insert(p.clone()) { - if let Some(cones) = self.cone(&p) { - for q in cones { - stack.push(q.clone()); - } - } - let supports = self.support(&p); - for q in supports { - stack.push(q.clone()); - } - } + /// The star includes the entity itself and all entities that directly cover it (supports). + pub fn star(&self, point: &MeshEntity) -> DashMap { + let result = DashMap::new(); + result.insert(point.clone(), ()); + let supports = self.support(point); + for support in supports { + result.insert(support, ()); } result } /// Retrieves all entities that support the given entity (`point`). /// These are the entities that have an arrow pointing to `point`. - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex = MeshEntity::Vertex(1); - /// let edge = MeshEntity::Edge(1); - /// sieve.add_arrow(vertex, edge); - /// let support_set = sieve.support(&edge); - /// - pub fn support(&self, point: &MeshEntity) -> FxHashSet { - let adjacency = self.adjacency.read().unwrap(); - adjacency - .iter() - .filter_map(|(from, to_set)| if to_set.contains(point) { Some(from.clone()) } else { None }) - .collect() + pub fn support(&self, point: &MeshEntity) -> Vec { + let mut supports = Vec::new(); + self.adjacency.iter().for_each(|entry| { + let from = entry.key(); + if entry.value().contains_key(point) { + supports.push(from.clone()); + } + }); + supports } /// Computes the meet operation for two entities, `p` and `q`. - /// This is the intersection of their closures (minimal separator). - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex1 = MeshEntity::Vertex(1); - /// let vertex2 = MeshEntity::Vertex(2); - /// let edge = MeshEntity::Edge(1); - /// sieve.add_arrow(vertex1, edge); - /// sieve.add_arrow(vertex2, edge); - /// let meet_set = sieve.meet(&vertex1, &vertex2); - /// - pub fn meet(&self, p: &MeshEntity, q: &MeshEntity) -> FxHashSet { + /// This is the intersection of their closures. + pub fn meet(&self, p: &MeshEntity, q: &MeshEntity) -> DashMap { let closure_p = self.closure(p); let closure_q = self.closure(q); - closure_p.intersection(&closure_q).cloned().collect() + let result = DashMap::new(); + + closure_p.iter().for_each(|entry| { + let key = entry.key(); + if closure_q.contains_key(key) { + result.insert(key.clone(), ()); + } + }); + + result } /// Computes the join operation for two entities, `p` and `q`. - /// This is the union of their stars (minimal separator). - /// - /// Example usage: - /// - /// let sieve = Sieve::new(); - /// let vertex1 = MeshEntity::Vertex(1); - /// let vertex2 = MeshEntity::Vertex(2); - /// let edge = MeshEntity::Edge(1); - /// sieve.add_arrow(vertex1, edge); - /// sieve.add_arrow(vertex2, edge); - /// let join_set = sieve.join(&vertex1, &vertex2); - /// - pub fn join(&self, p: &MeshEntity, q: &MeshEntity) -> FxHashSet { + /// This is the union of their stars. + pub fn join(&self, p: &MeshEntity, q: &MeshEntity) -> DashMap { let star_p = self.star(p); let star_q = self.star(q); + let result = DashMap::new(); + + star_p.iter().for_each(|entry| { + result.insert(entry.key().clone(), ()); + }); + star_q.iter().for_each(|entry| { + result.insert(entry.key().clone(), ()); + }); - let join_result: FxHashSet = star_p.union(&star_q).cloned().collect(); - join_result + result } /// Applies a given function in parallel to all adjacency map entries. /// This function is executed concurrently over each entity and its /// corresponding set of related entities. - /// - /// Example usage: - /// - /// sieve.par_for_each_adjacent(|(entity, adj_set)| { - /// println!("{:?} -> {:?}", entity, adj_set); - /// }); - /// pub fn par_for_each_adjacent(&self, func: F) where - F: Fn((&MeshEntity, &FxHashSet)) + Sync + Send, + F: Fn((&MeshEntity, Vec)) + Sync + Send, { - let adjacency = self.adjacency.read().unwrap(); - adjacency.par_iter().for_each(func); + // Collect entries from DashMap to avoid borrow conflicts + let entries: Vec<_> = self.adjacency.iter().map(|entry| { + let key = entry.key().clone(); + let values: Vec = entry.value().iter().map(|e| e.key().clone()).collect(); + (key, values) + }).collect(); + + // Execute in parallel over collected entries + entries.par_iter().for_each(|entry| { + func((&entry.0, entry.1.clone())); + }); } } + #[cfg(test)] mod tests { use super::*; @@ -213,7 +153,7 @@ mod tests { #[test] /// Test that verifies adding an arrow between two entities and querying - /// the cone of an entity works as expected. + /// the cone of an entity works as expected. fn test_add_arrow_and_cone() { let sieve = Sieve::new(); let vertex = MeshEntity::Vertex(1); @@ -225,19 +165,24 @@ mod tests { #[test] /// Test that verifies the closure of a vertex correctly includes - /// all transitive relationships. + /// all transitive relationships and the entity itself. fn test_closure() { let sieve = Sieve::new(); let vertex = MeshEntity::Vertex(1); let edge = MeshEntity::Edge(1); + let face = MeshEntity::Face(1); sieve.add_arrow(vertex, edge); + sieve.add_arrow(edge, face); let closure_result = sieve.closure(&vertex); - assert!(closure_result.contains(&edge)); + assert!(closure_result.contains_key(&vertex)); + assert!(closure_result.contains_key(&edge)); + assert!(closure_result.contains_key(&face)); + assert_eq!(closure_result.len(), 3); } #[test] /// Test that verifies the support of an entity includes the - /// correct supporting entities. + /// correct supporting entities. fn test_support() { let sieve = Sieve::new(); let vertex = MeshEntity::Vertex(1); @@ -247,64 +192,57 @@ mod tests { let support_result = sieve.support(&edge); assert!(support_result.contains(&vertex)); + assert_eq!(support_result.len(), 1); } #[test] - /// Test that verifies the star of an entity includes both cone and - /// support sets of the entity. + /// Test that verifies the star of an entity includes both the entity itself and + /// its immediate supports. fn test_star() { let sieve = Sieve::new(); - let vertex = MeshEntity::Vertex(1); let edge = MeshEntity::Edge(1); let face = MeshEntity::Face(1); - sieve.add_arrow(vertex, edge); sieve.add_arrow(edge, face); let star_result = sieve.star(&face); - assert!(star_result.contains(&edge)); - assert!(star_result.contains(&vertex)); + assert!(star_result.contains_key(&face)); + assert!(star_result.contains_key(&edge)); + assert_eq!(star_result.len(), 2); } #[test] /// Test that verifies the meet operation between two entities returns - /// the correct minimal separator (intersection of closures). + /// the correct intersection of their closures. fn test_meet() { let sieve = Sieve::new(); let vertex1 = MeshEntity::Vertex(1); let vertex2 = MeshEntity::Vertex(2); let edge = MeshEntity::Edge(1); - let face = MeshEntity::Face(1); sieve.add_arrow(vertex1, edge); sieve.add_arrow(vertex2, edge); - sieve.add_arrow(edge, face); let meet_result = sieve.meet(&vertex1, &vertex2); - assert!(meet_result.contains(&edge)); + assert!(meet_result.contains_key(&edge)); + assert_eq!(meet_result.len(), 1); } #[test] /// Test that verifies the join operation between two entities returns - /// the correct minimal separator (union of stars). + /// the correct union of their stars. fn test_join() { let sieve = Sieve::new(); let vertex1 = MeshEntity::Vertex(1); let vertex2 = MeshEntity::Vertex(2); - let edge = MeshEntity::Edge(1); - let face = MeshEntity::Face(1); - - sieve.add_arrow(vertex1, edge); - sieve.add_arrow(vertex2, edge); - sieve.add_arrow(edge, face); let join_result = sieve.join(&vertex1, &vertex2); - assert!(join_result.contains(&vertex1), "Join result should contain vertex1"); - assert!(join_result.contains(&vertex2), "Join result should contain vertex2"); - assert!(join_result.contains(&edge), "Join result should contain the edge"); - assert!(join_result.contains(&face), "Join result should contain the face"); + assert!(join_result.contains_key(&vertex1), "Join result should contain vertex1"); + assert!(join_result.contains_key(&vertex2), "Join result should contain vertex2"); + assert_eq!(join_result.len(), 2); } } + diff --git a/src/domain/stratify.rs b/src/domain/stratify.rs index 11d0474..a1f1cfd 100644 --- a/src/domain/stratify.rs +++ b/src/domain/stratify.rs @@ -1,6 +1,6 @@ use crate::domain::mesh_entity::MeshEntity; use crate::domain::sieve::Sieve; -use rustc_hash::FxHashMap; +use dashmap::DashMap; /// Implements a stratification method for the `Sieve` structure. /// Stratification organizes the mesh entities into different strata based on @@ -11,7 +11,7 @@ use rustc_hash::FxHashMap; /// - Stratum 3: Cells /// /// This method categorizes each `MeshEntity` into its corresponding stratum and -/// returns a map where the keys are the dimension (stratum) and the values +/// returns a `DashMap` where the keys are the dimension (stratum) and the values /// are vectors of mesh entities in that stratum. /// /// Example usage: @@ -34,14 +34,12 @@ impl Sieve { /// sieve.add_arrow(MeshEntity::Vertex(1), MeshEntity::Edge(1)); /// let strata = sieve.stratify(); /// - pub fn stratify(&self) -> FxHashMap> { - let mut strata: FxHashMap> = FxHashMap::default(); - - // Acquire a read lock to access the adjacency data. - let adjacency = self.adjacency.read().unwrap(); - + pub fn stratify(&self) -> DashMap> { + let strata: DashMap> = DashMap::new(); + // Iterate over the adjacency map to classify entities by their dimension. - for (entity, _) in adjacency.iter() { + self.adjacency.iter().for_each(|entry| { + let entity = entry.key(); // Determine the dimension of the current entity. let dimension = match entity { MeshEntity::Vertex(_) => 0, // Stratum 0 for vertices @@ -50,10 +48,10 @@ impl Sieve { MeshEntity::Cell(_) => 3, // Stratum 3 for cells }; - // Add the entity to the appropriate stratum in the map. - strata.entry(dimension).or_insert_with(Vec::new).push(*entity); - } - + // Insert entity into the appropriate stratum in a thread-safe manner. + strata.entry(dimension).or_insert_with(Vec::new).push(entity.clone()); + }); + strata } } diff --git a/src/geometry/docs/about_geometry.md b/src/geometry/docs/about_geometry.md index 86426c2..21c0969 100644 --- a/src/geometry/docs/about_geometry.md +++ b/src/geometry/docs/about_geometry.md @@ -1,22 +1,25 @@ ### Overview of the `src/geometry/` Module #### Overview -The `src/geometry/` module in the Hydra project is responsible for performing geometric calculations on 3D shapes commonly used in geophysical fluid dynamics simulations, such as prisms, pyramids, hexahedrons, and tetrahedrons. This module provides essential methods to compute properties like the **centroid** (the geometric center) and **volume** of these cells, which are crucial for finite volume methods (FVM) that solve partial differential equations over complex domains. The geometry module is built with extensibility in mind, supporting various 3D cell shapes through modular functions and decomposition techniques. +The `src/geometry/` module in the Hydra project provides essential geometric calculations for 3D shapes used in geophysical fluid dynamics simulations, including prisms, pyramids, hexahedrons, tetrahedrons, and more. This module underpins the finite volume methods (FVM) used to solve partial differential equations across complex domains. The module emphasizes computational efficiency, mathematical rigor, and robustness to handle both regular and degenerate shapes that may arise in realistic simulations. Key operations in this module include: -- Centroid calculation for common geometric shapes (hexahedrons, prisms, pyramids, and tetrahedrons) -- Volume computation for both simple and complex shapes (e.g., through tetrahedral decomposition) -- Support for degenerate cases where cells collapse into lower-dimensional shapes, ensuring robustness in numerical simulations +- **Centroid Calculation**: Determines the geometric center of various shapes. +- **Volume Computation**: Uses techniques like tetrahedral decomposition for complex shapes. +- **Face Normal Calculation**: Computes outward-pointing normal vectors for face flux calculations. +- **Handling Degenerate Cases**: Manages cases where shapes collapse into lower dimensions. -#### Key Classes and Functions +#### Key Structures and Functions -##### 1. **Centroid Calculation** +The primary struct, `Geometry`, encapsulates the geometric data and methods for a mesh, including vertex positions, centroids, volumes, and cached computed properties to avoid redundant calculations. -Centroid calculation determines the average position of all vertices in a geometric shape. This is critical in mesh-based simulations where centroids often serve as reference points for finite volume methods or flux calculations. +##### 1. **Centroid Calculation** +Centroid calculations are vital for determining the center of mass of each control volume. These centroids are used as reference points in numerical methods for flux integration, and they enable efficient computation by avoiding repetitive geometric queries. - **`compute_hexahedron_centroid`**: - - Calculates the centroid of a hexahedron (cube or cuboid) by averaging the coordinates of its 8 vertices. - - Example: + - Calculates the centroid of a hexahedron (e.g., cube or cuboid) by averaging the coordinates of its 8 vertices. + - This method is optimized for stability and handles both regular and degenerate cases where vertices may collapse onto a plane. + - **Example Usage**: ```rust,ignore let geometry = Geometry::new(); let hexahedron_vertices = vec![ @@ -28,8 +31,8 @@ Centroid calculation determines the average position of all vertices in a geomet ``` - **`compute_prism_centroid`**: - - Computes the centroid of a triangular prism by calculating the centroids of the top and bottom triangles and averaging them. - - Example: + - Calculates the centroid of a triangular prism by averaging the centroids of the top and bottom triangular faces. + - **Example Usage**: ```rust,ignore let geometry = Geometry::new(); let prism_vertices = vec![ @@ -41,12 +44,11 @@ Centroid calculation determines the average position of all vertices in a geomet ``` ##### 2. **Volume Calculation** - -The volume of a 3D geometric cell is a fundamental quantity in finite volume methods, determining the amount of flux passing through a volume. This module supports volume computation for common geometric shapes. +Volume computation enables the determination of flux through a control volume, which is central to FVM methods. Complex shapes are decomposed into tetrahedrons for efficient and accurate volume computation. - **`compute_hexahedron_volume`**: - Computes the volume of a hexahedron by decomposing it into 5 tetrahedrons and summing their volumes. - - Example: + - **Example Usage**: ```rust,ignore let geometry = Geometry::new(); let hexahedron_vertices = vec![ @@ -54,34 +56,33 @@ The volume of a 3D geometric cell is a fundamental quantity in finite volume met [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0], [0.0, 1.0, 1.0], ]; let volume = geometry.compute_hexahedron_volume(&hexahedron_vertices); - assert!((volume - 1.0).abs() < 1e-10); // Volume of a unit cube + assert!((volume - 1.0).abs() < 1e-10); // Volume of a unit cube ``` -- **`compute_prism_volume`**: - - Computes the volume of a triangular prism by multiplying the area of the base (bottom triangle) by the height (distance between the centroids of the top and bottom triangles). - - Example: +- **`compute_pyramid_volume`**: + - Calculates the volume of a square or triangular-based pyramid by breaking down the geometry into tetrahedrons and summing their volumes. + - **Example Usage**: ```rust,ignore let geometry = Geometry::new(); - let prism_vertices = vec![ - [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], - [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], + let pyramid_vertices = vec![ + [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], [0.5, 0.5, 1.0], ]; - let volume = geometry.compute_prism_volume(&prism_vertices); - assert!((volume - 0.5).abs() < 1e-10); + let volume = geometry.compute_pyramid_volume(&pyramid_vertices); + assert!((volume - 1.0 / 3.0).abs() < 1e-10); ``` -##### 3. **Tetrahedral Decomposition** +##### 3. **Face Normal Calculation** +Face normal vectors are crucial for FVM, as they define the direction and magnitude of flux through faces. Normal calculations are available for triangular and quadrilateral faces. -For more complex shapes, the geometry module uses **tetrahedral decomposition** to compute volumes. This technique is applied in the volume calculation of hexahedrons, where the shape is decomposed into 5 tetrahedrons. The volume of each tetrahedron is calculated and summed to get the total volume of the hexahedron. - -- **`compute_tetrahedron_volume`**: - - This function is used internally to calculate the volume of individual tetrahedrons, which is then used in the decomposition process. +- **`compute_triangle_normal`**: + - Computes the normal vector of a triangle face using the cross product of two edge vectors. +- **`compute_quadrilateral_normal`**: + - Computes the normal vector of a quadrilateral by dividing it into two triangles and averaging their normals. ##### 4. **Handling Degenerate Cases** +For robustness, the module handles degenerate cases where cells collapse into lower-dimensional shapes. For instance, if a hexahedron collapses into a plane, the volume will be computed as zero, and centroid calculations will default to the average position of the vertices, ensuring stable simulations under various conditions. -The module is designed to handle **degenerate cases**, where cells collapse into lower-dimensional shapes (e.g., all vertices of a hexahedron lie on the same plane). In such cases, the volume should be zero, and the module appropriately returns this value. This ensures that numerical simulations remain stable even when encountering degenerate geometries. - -- Example: +- **Example of Degenerate Handling**: ```rust,ignore let geometry = Geometry::new(); let degenerate_hexahedron_vertices = vec![ @@ -89,16 +90,17 @@ The module is designed to handle **degenerate cases**, where cells collapse into [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0], ]; let volume = geometry.compute_hexahedron_volume(°enerate_hexahedron_vertices); - assert_eq!(volume, 0.0); // Degenerate hexahedron should have zero volume + assert_eq!(volume, 0.0); ``` -#### Test Coverage +### Test Coverage -The module includes comprehensive unit tests for both centroid and volume calculations. Each geometric shape has multiple test cases, covering both regular and degenerate scenarios. These tests ensure that the module behaves correctly across a wide range of input cases and helps prevent errors when adding new features or modifications. +The module includes rigorous unit testing for all core functions: +- **Regular Shape Tests**: Verifies centroid, volume, and normal calculations for standard shapes like cubes and prisms. +- **Degenerate Shape Tests**: Ensures that the module returns zero volumes for degenerate shapes, handles edge cases in centroid calculations, and tests normals for collapsed faces. -- **Regular Case Tests**: Tests for normal shapes, such as cubes and prisms, to verify the correctness of centroid and volume calculations. -- **Degenerate Case Tests**: Tests for degenerate shapes, such as a collapsed hexahedron or prism, ensuring that the computed volume is zero and that centroids are handled appropriately. +These tests are integral for ensuring numerical stability across a wide range of shapes and configurations. -#### Summary +### Summary -The `src/geometry/` module in Hydra provides core geometric utilities for mesh-based simulations. It offers highly efficient and mathematically robust methods for calculating centroids and volumes for common 3D cells like hexahedrons, prisms, and tetrahedrons. The ability to handle degenerate cases ensures the reliability and stability of numerical simulations, while the use of tetrahedral decomposition enables complex shape handling. Overall, this module serves as a fundamental building block in Hydra’s mesh handling and simulation processes. \ No newline at end of file +The `src/geometry/` module in Hydra provides fundamental tools for calculating geometric properties within FVM-based simulations. Its capabilities in handling centroids, volumes, and normals for a variety of 3D shapes ensure it can meet the needs of complex geophysical simulations. The module’s flexibility and robust handling of degenerate cases make it a reliable foundation for advancing Hydra's mesh handling and simulation processes. \ No newline at end of file diff --git a/src/geometry/hexahedron.rs b/src/geometry/hexahedron.rs index e271051..2766683 100644 --- a/src/geometry/hexahedron.rs +++ b/src/geometry/hexahedron.rs @@ -14,36 +14,22 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not exactly 8. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let hexahedron_vertices = vec![ - /// [0.0, 0.0, 0.0], // vertex 1 - /// [1.0, 0.0, 0.0], // vertex 2 - /// [1.0, 1.0, 0.0], // vertex 3 - /// [0.0, 1.0, 0.0], // vertex 4 - /// [0.0, 0.0, 1.0], // vertex 5 - /// [1.0, 0.0, 1.0], // vertex 6 - /// [1.0, 1.0, 1.0], // vertex 7 - /// [0.0, 1.0, 1.0], // vertex 8 - /// ]; - /// let centroid = geometry.compute_hexahedron_centroid(&hexahedron_vertices); - /// assert_eq!(centroid, [0.5, 0.5, 0.5]); // Centroid of a unit cube - /// ``` pub fn compute_hexahedron_centroid(&self, cell_vertices: &Vec<[f64; 3]>) -> [f64; 3] { + assert!(cell_vertices.len() == 8, "Hexahedron must have exactly 8 vertices"); + let mut centroid = [0.0, 0.0, 0.0]; - for v in cell_vertices { - centroid[0] += v[0]; - centroid[1] += v[1]; - centroid[2] += v[2]; + for vertex in cell_vertices { + centroid[0] += vertex[0]; + centroid[1] += vertex[1]; + centroid[2] += vertex[2]; } + let num_vertices = cell_vertices.len() as f64; - centroid[0] /= num_vertices; - centroid[1] /= num_vertices; - centroid[2] /= num_vertices; - centroid + [ + centroid[0] / num_vertices, + centroid[1] / num_vertices, + centroid[2] / num_vertices, + ] } /// Computes the volume of a hexahedral cell using tetrahedral decomposition. @@ -61,24 +47,6 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not exactly 8. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let hexahedron_vertices = vec![ - /// [0.0, 0.0, 0.0], // vertex 1 - /// [1.0, 0.0, 0.0], // vertex 2 - /// [1.0, 1.0, 0.0], // vertex 3 - /// [0.0, 1.0, 0.0], // vertex 4 - /// [0.0, 0.0, 1.0], // vertex 5 - /// [1.0, 0.0, 1.0], // vertex 6 - /// [1.0, 1.0, 1.0], // vertex 7 - /// [0.0, 1.0, 1.0], // vertex 8 - /// ]; - /// let volume = geometry.compute_hexahedron_volume(&hexahedron_vertices); - /// assert!((volume - 1.0).abs() < 1e-10); // Volume of a unit cube - /// ``` pub fn compute_hexahedron_volume(&self, cell_vertices: &Vec<[f64; 3]>) -> f64 { assert!(cell_vertices.len() == 8, "Hexahedron must have exactly 8 vertices"); diff --git a/src/geometry/mod.rs b/src/geometry/mod.rs index 68177dd..52f1a14 100644 --- a/src/geometry/mod.rs +++ b/src/geometry/mod.rs @@ -1,6 +1,8 @@ use rayon::prelude::*; use rustc_hash::FxHashMap; use crate::domain::{mesh::Mesh, MeshEntity}; +use faer::{Mat, zipped, unzipped}; +use std::sync::Mutex; // Module for handling geometric data and computations // 2D Shape Modules @@ -14,28 +16,22 @@ pub mod pyramid; /// The `Geometry` struct stores geometric data for a mesh, including vertex coordinates, /// cell centroids, and volumes. It also maintains a cache of computed properties such as -/// volume and centroid for reuse, optimizing performance. -/// -/// Example usage: -/// -/// let mut geometry = Geometry::new(); -/// geometry.set_vertex(0, [1.0, 2.0, 3.0]); -/// let centroid = geometry.compute_cell_centroid(&mesh, &cell); -/// +/// volume and centroid for reuse, optimizing performance by avoiding redundant calculations. pub struct Geometry { - pub vertices: Vec<[f64; 3]>, // 3D coordinates for each vertex - pub cell_centroids: Vec<[f64; 3]>, // Centroid positions for each cell - pub cell_volumes: Vec, // Volumes of each cell - pub cache: FxHashMap, // Cache for computed properties + pub vertices: Vec<[f64; 3]>, // 3D coordinates for each vertex + pub cell_centroids: Vec<[f64; 3]>, // Centroid positions for each cell + pub cell_volumes: Vec, // Volumes of each cell + pub cache: Mutex>, // Cache for computed properties, with thread safety } /// The `GeometryCache` struct stores computed properties of geometric entities, -/// including volume, centroid, and area. It helps avoid recomputation and thus improves efficiency. +/// including volume, centroid, and area, with an optional "dirty" flag for lazy evaluation. #[derive(Default)] pub struct GeometryCache { pub volume: Option, pub centroid: Option<[f64; 3]>, pub area: Option, + pub normal: Option<[f64; 3]>, // Stores a precomputed normal vector for a face } /// `CellShape` enumerates the different cell shapes in a mesh, including: @@ -60,153 +56,86 @@ pub enum FaceShape { Quadrilateral, } -// Initial implementation of the Geometry module impl Geometry { /// Initializes a new `Geometry` instance with empty data. - /// - /// Example usage: - /// - /// let geometry = Geometry::new(); - /// pub fn new() -> Geometry { Geometry { vertices: Vec::new(), cell_centroids: Vec::new(), cell_volumes: Vec::new(), - cache: FxHashMap::default(), + cache: Mutex::new(FxHashMap::default()), } } - /// Adds or updates a vertex in the geometry. If the vertex already exists (based on ID or index), + /// Adds or updates a vertex in the geometry. If the vertex already exists, /// it updates its coordinates; otherwise, it adds a new vertex. /// /// # Arguments - /// * `vertex_index` - The index or ID of the vertex to be set or updated. - /// * `coords` - The 3D coordinates of the vertex as an array of `[f64; 3]`. - /// - /// Example usage: - /// - /// geometry.set_vertex(0, [1.0, 2.0, 3.0]); - /// + /// * `vertex_index` - The index of the vertex. + /// * `coords` - The 3D coordinates of the vertex. pub fn set_vertex(&mut self, vertex_index: usize, coords: [f64; 3]) { - // Ensure the vector has enough capacity to accommodate the vertex at the given index if vertex_index >= self.vertices.len() { - // Resize the vertices vector to hold up to `vertex_index` and beyond self.vertices.resize(vertex_index + 1, [0.0, 0.0, 0.0]); } - - // Update the vertex at the specified index self.vertices[vertex_index] = coords; - self.invalidate_cache(); // Invalidate cache when vertices change + self.invalidate_cache(); } - - /// Computes the centroid of a given cell based on its shape and vertices. - /// - /// # Arguments - /// * `mesh` - A reference to the mesh structure containing the cells and vertices. - /// * `cell` - A reference to the `MeshEntity` representing the cell for which the centroid is being computed. - /// - /// # Returns - /// * `[f64; 3]` - The 3D coordinates of the cell centroid. - /// - /// Example usage: - /// - /// let centroid = geometry.compute_cell_centroid(&mesh, &cell); - /// + + /// Computes and returns the centroid of a specified cell using the cell's shape and vertices. + /// Caches the result for reuse. pub fn compute_cell_centroid(&mut self, mesh: &Mesh, cell: &MeshEntity) -> [f64; 3] { let cell_id = cell.id(); - if let Some(cached) = self.cache.get(&cell_id).and_then(|c| c.centroid) { + if let Some(cached) = self.cache.lock().unwrap().get(&cell_id).and_then(|c| c.centroid) { return cached; } - let cell_shape = mesh.get_cell_shape(cell); + let cell_shape = mesh.get_cell_shape(cell).expect("Cell shape not found"); let cell_vertices = mesh.get_cell_vertices(cell); let centroid = match cell_shape { - Ok(CellShape::Tetrahedron) => self.compute_tetrahedron_centroid(&cell_vertices), - Ok(CellShape::Hexahedron) => self.compute_hexahedron_centroid(&cell_vertices), - Ok(CellShape::Prism) => self.compute_prism_centroid(&cell_vertices), - Ok(CellShape::Pyramid) => self.compute_pyramid_centroid(&cell_vertices), - Err(_) => todo!(), + CellShape::Tetrahedron => self.compute_tetrahedron_centroid(&cell_vertices), + CellShape::Hexahedron => self.compute_hexahedron_centroid(&cell_vertices), + CellShape::Prism => self.compute_prism_centroid(&cell_vertices), + CellShape::Pyramid => self.compute_pyramid_centroid(&cell_vertices), }; - self.cache.entry(cell_id).or_default().centroid = Some(centroid); + self.cache.lock().unwrap().entry(cell_id).or_default().centroid = Some(centroid); centroid } - /// Computes the volume of a given cell based on its shape and vertex coordinates. - /// - /// # Arguments - /// * `mesh` - A reference to the mesh structure containing the cells and vertices. - /// * `cell` - A reference to the `MeshEntity` representing the cell for which the volume is being computed. - /// - /// # Returns - /// * `f64` - The volume of the cell. - /// - /// Example usage: - /// - /// let volume = geometry.compute_cell_volume(&mesh, &cell); - /// + /// Computes the volume of a given cell using its shape and vertex coordinates. + /// The computed volume is cached for efficiency. pub fn compute_cell_volume(&mut self, mesh: &Mesh, cell: &MeshEntity) -> f64 { let cell_id = cell.id(); - if let Some(cached) = self.cache.get(&cell_id).and_then(|c| c.volume) { + if let Some(cached) = self.cache.lock().unwrap().get(&cell_id).and_then(|c| c.volume) { return cached; } - let cell_shape = mesh.get_cell_shape(cell); + let cell_shape = mesh.get_cell_shape(cell).expect("Cell shape not found"); let cell_vertices = mesh.get_cell_vertices(cell); let volume = match cell_shape { - Ok(CellShape::Tetrahedron) => self.compute_tetrahedron_volume(&cell_vertices), - Ok(CellShape::Hexahedron) => self.compute_hexahedron_volume(&cell_vertices), - Ok(CellShape::Prism) => self.compute_prism_volume(&cell_vertices), - Ok(CellShape::Pyramid) => self.compute_pyramid_volume(&cell_vertices), - Err(_) => todo!(), + CellShape::Tetrahedron => self.compute_tetrahedron_volume(&cell_vertices), + CellShape::Hexahedron => self.compute_hexahedron_volume(&cell_vertices), + CellShape::Prism => self.compute_prism_volume(&cell_vertices), + CellShape::Pyramid => self.compute_pyramid_volume(&cell_vertices), }; - self.cache.entry(cell_id).or_default().volume = Some(volume); + self.cache.lock().unwrap().entry(cell_id).or_default().volume = Some(volume); volume } - /// Computes the Euclidean distance between two points in 3D space. - /// - /// # Arguments - /// * `p1` - The first point as an array of 3D coordinates `[f64; 3]`. - /// * `p2` - The second point as an array of 3D coordinates `[f64; 3]`. - /// - /// # Returns - /// * `f64` - The Euclidean distance between the two points. - /// - /// Example usage: - /// - /// let distance = Geometry::compute_distance(&p1, &p2); - /// + /// Calculates Euclidean distance between two points in 3D space. pub fn compute_distance(p1: &[f64; 3], p2: &[f64; 3]) -> f64 { - // Compute the squared differences for each coordinate (x, y, z) let dx = p1[0] - p2[0]; let dy = p1[1] - p2[1]; let dz = p1[2] - p2[2]; - - // Apply the Euclidean distance formula: sqrt(dx^2 + dy^2 + dz^2) (dx.powi(2) + dy.powi(2) + dz.powi(2)).sqrt() } - /// Computes the area of a 2D face based on its shape. - /// - /// # Arguments - /// * `face_id` - The ID of the face to cache the result. - /// * `face_shape` - Enum defining the shape of the face (e.g., Triangle, Quadrilateral). - /// * `face_vertices` - A vector of vertices representing the 3D coordinates of the face's vertices. - /// - /// # Returns - /// * `f64` - The area of the face. - /// - /// Example usage: - /// - /// let area = geometry.compute_face_area(1, FaceShape::Triangle, &triangle_vertices); - /// + /// Computes the area of a 2D face based on its shape, caching the result. pub fn compute_face_area(&mut self, face_id: usize, face_shape: FaceShape, face_vertices: &Vec<[f64; 3]>) -> f64 { - if let Some(cached) = self.cache.get(&face_id).and_then(|c| c.area) { + if let Some(cached) = self.cache.lock().unwrap().get(&face_id).and_then(|c| c.area) { return cached; } @@ -215,7 +144,7 @@ impl Geometry { FaceShape::Quadrilateral => self.compute_quadrilateral_area(face_vertices), }; - self.cache.entry(face_id).or_default().area = Some(area); + self.cache.lock().unwrap().entry(face_id).or_default().area = Some(area); area } @@ -234,63 +163,31 @@ impl Geometry { } } - /// Clears the cache when geometry changes, such as when vertices are updated. - /// - /// Example usage: - /// - /// geometry.invalidate_cache(); - /// + /// Invalidate the cache when geometry changes (e.g., vertex updates). fn invalidate_cache(&mut self) { - self.cache.clear(); + self.cache.lock().unwrap().clear(); } - /// Computes the total volume of all cells in the geometry. - /// - /// # Returns - /// * `f64` - The total volume of all cells. - /// - /// Example usage: - /// - /// let total_volume = geometry.compute_total_volume(); - /// + /// Computes the total volume of all cells. pub fn compute_total_volume(&self) -> f64 { self.cell_volumes.par_iter().sum() } - /// Updates all cell volumes in parallel by computing them for each cell in the mesh. - /// - /// # Arguments - /// * `mesh` - A reference to the `Mesh` structure containing the cells and vertices. - /// - /// Example usage: - /// - /// geometry.update_all_cell_volumes(&mesh); - /// + /// Updates all cell volumes in parallel using mesh information. pub fn update_all_cell_volumes(&mut self, mesh: &Mesh) { - // Using a mutable reference inside the parallel iteration requires `Mutex` or similar for thread safety. let new_volumes: Vec = mesh .get_cells() .par_iter() .map(|cell| { - // Use a temporary Geometry instance to compute the volume. let mut temp_geometry = Geometry::new(); temp_geometry.compute_cell_volume(mesh, cell) }) .collect(); - - // Update the cell volumes after the parallel iteration. + self.cell_volumes = new_volumes; } /// Computes the total centroid of all cells. - /// - /// # Returns - /// * `[f64; 3]` - The total centroid of all cells. - /// - /// Example usage: - /// - /// let total_centroid = geometry.compute_total_centroid(); - /// pub fn compute_total_centroid(&self) -> [f64; 3] { let total_centroid: [f64; 3] = self.cell_centroids .par_iter() @@ -311,9 +208,9 @@ impl Geometry { total_centroid[2] / num_centroids, ] } - } + #[cfg(test)] mod tests { use crate::geometry::{Geometry, CellShape, FaceShape}; @@ -391,12 +288,8 @@ mod tests { mesh.add_arrow(cell, vertex3); mesh.add_arrow(cell, vertex4); - // Debug: Print the sieve structure. - println!("Sieve contents: {:?}", mesh.sieve); - // Verify that `get_cell_vertices` retrieves the correct vertices. let cell_vertices = mesh.get_cell_vertices(&cell); - println!("Cell vertices: {:?}", cell_vertices); assert_eq!(cell_vertices.len(), 4, "Expected 4 vertices for a tetrahedron cell."); // Validate the shape before computing. @@ -444,4 +337,17 @@ mod tests { // Expected centroid is the geometric center: (0.5, 0.5, 0.0) assert_eq!(centroid, [0.5, 0.5, 0.0]); } + + #[test] + fn test_compute_total_volume() { + let mut geometry = Geometry::new(); + let mut mesh = Mesh::new(); + + // Example setup: Define cells with known volumes + // Here, you would typically define several cells and their volumes for the test + geometry.cell_volumes = vec![1.0, 2.0, 3.0]; + + // Expected total volume is the sum of individual cell volumes: 1.0 + 2.0 + 3.0 = 6.0 + assert_eq!(geometry.compute_total_volume(), 6.0); + } } diff --git a/src/geometry/prism.rs b/src/geometry/prism.rs index 81cd933..bcbf48e 100644 --- a/src/geometry/prism.rs +++ b/src/geometry/prism.rs @@ -1,7 +1,6 @@ use crate::geometry::Geometry; impl Geometry { - /// Computes the centroid of a triangular prism. /// /// The prism is assumed to have two parallel triangular faces (top and bottom), @@ -17,22 +16,6 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not exactly 6. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let prism_vertices = vec![ - /// [0.0, 0.0, 0.0], // top triangle vertex 1 - /// [1.0, 0.0, 0.0], // top triangle vertex 2 - /// [0.0, 1.0, 0.0], // top triangle vertex 3 - /// [0.0, 0.0, 1.0], // bottom triangle vertex 1 - /// [1.0, 0.0, 1.0], // bottom triangle vertex 2 - /// [0.0, 1.0, 1.0], // bottom triangle vertex 3 - /// ]; - /// let centroid = geometry.compute_prism_centroid(&prism_vertices); - /// assert_eq!(centroid, [1.0 / 3.0, 1.0 / 3.0, 0.5]); - /// ``` pub fn compute_prism_centroid(&self, cell_vertices: &Vec<[f64; 3]>) -> [f64; 3] { assert!(cell_vertices.len() == 6, "Triangular prism must have exactly 6 vertices"); @@ -41,17 +24,15 @@ impl Geometry { let bottom_triangle = vec![cell_vertices[3], cell_vertices[4], cell_vertices[5]]; // Compute the centroids of both triangles - let top_centroid = self.compute_tetrahedron_centroid(&top_triangle); - let bottom_centroid = self.compute_tetrahedron_centroid(&bottom_triangle); + let top_centroid = self.compute_triangle_centroid(&top_triangle); + let bottom_centroid = self.compute_triangle_centroid(&bottom_triangle); // Compute the centroid of the prism by averaging the top and bottom centroids - let prism_centroid = [ + [ (top_centroid[0] + bottom_centroid[0]) / 2.0, (top_centroid[1] + bottom_centroid[1]) / 2.0, (top_centroid[2] + bottom_centroid[2]) / 2.0, - ]; - - prism_centroid + ] } /// Computes the volume of a triangular prism. @@ -68,35 +49,18 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not exactly 6. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let prism_vertices = vec![ - /// [0.0, 0.0, 0.0], // top triangle vertex 1 - /// [1.0, 0.0, 0.0], // top triangle vertex 2 - /// [0.0, 1.0, 0.0], // top triangle vertex 3 - /// [0.0, 0.0, 1.0], // bottom triangle vertex 1 - /// [1.0, 0.0, 1.0], // bottom triangle vertex 2 - /// [0.0, 1.0, 1.0], // bottom triangle vertex 3 - /// ]; - /// let volume = geometry.compute_prism_volume(&prism_vertices); - /// assert!((volume - 0.5).abs() < 1e-10); - /// ``` pub fn compute_prism_volume(&self, cell_vertices: &Vec<[f64; 3]>) -> f64 { assert!(cell_vertices.len() == 6, "Triangular prism must have exactly 6 vertices"); // Split into top and bottom triangles - let top_triangle = vec![cell_vertices[0], cell_vertices[1], cell_vertices[2]]; let bottom_triangle = vec![cell_vertices[3], cell_vertices[4], cell_vertices[5]]; // Compute the area of the base triangle (bottom triangle) let base_area = self.compute_triangle_area(&bottom_triangle); // Compute the height of the prism as the distance between the top and bottom triangle centroids - let top_centroid = self.compute_tetrahedron_centroid(&top_triangle); - let bottom_centroid = self.compute_tetrahedron_centroid(&bottom_triangle); + let top_centroid = self.compute_triangle_centroid(&cell_vertices[0..3].to_vec()); + let bottom_centroid = self.compute_triangle_centroid(&bottom_triangle); let height = Geometry::compute_distance(&top_centroid, &bottom_centroid); // Volume of the prism = base area * height diff --git a/src/geometry/pyramid.rs b/src/geometry/pyramid.rs index 9d340b5..1d9e98e 100644 --- a/src/geometry/pyramid.rs +++ b/src/geometry/pyramid.rs @@ -6,8 +6,8 @@ impl Geometry { /// The centroid is computed for a pyramid with either a triangular or square base. /// The pyramid is represented by 4 vertices (triangular base) or 5 vertices (square base). /// - /// - For a triangular base, the function treats the pyramid as a tetrahedron. - /// - For a square base, the function splits the pyramid into two tetrahedrons and combines their centroids. + /// - For a triangular base, the pyramid is treated as a tetrahedron. + /// - For a square base, the pyramid is split into two tetrahedrons, and the weighted centroid of both is calculated. /// /// # Arguments /// * `cell_vertices` - A vector of 4 vertices for a triangular-based pyramid, or 5 vertices for a square-based pyramid. @@ -17,76 +17,27 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not 4 (for a triangular base) or 5 (for a square base). - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let pyramid_vertices = vec![ - /// [0.0, 0.0, 0.0], // base vertex 1 - /// [1.0, 0.0, 0.0], // base vertex 2 - /// [1.0, 1.0, 0.0], // base vertex 3 - /// [0.0, 1.0, 0.0], // base vertex 4 - /// [0.5, 0.5, 1.0], // apex - /// ]; - /// let centroid = geometry.compute_pyramid_centroid(&pyramid_vertices); - /// assert_eq!(centroid, [0.5, 0.5, 0.25]); - /// ``` pub fn compute_pyramid_centroid(&self, cell_vertices: &Vec<[f64; 3]>) -> [f64; 3] { - let mut _total_volume = 0.0; - let mut weighted_centroid = [0.0, 0.0, 0.0]; - - let num_vertices = cell_vertices.len(); - assert!(num_vertices == 4 || num_vertices == 5, "Pyramid must have 4 (triangular) or 5 (square) vertices"); + assert!(cell_vertices.len() == 4 || cell_vertices.len() == 5, "Pyramid must have 4 or 5 vertices"); - // Define the apex and base vertices - let apex = if num_vertices == 4 { cell_vertices[3] } else { cell_vertices[4] }; - let base_vertices = if num_vertices == 4 { - vec![cell_vertices[0], cell_vertices[1], cell_vertices[2]] // Triangular base + let apex = if cell_vertices.len() == 4 { cell_vertices[3] } else { cell_vertices[4] }; + let base_vertices = if cell_vertices.len() == 4 { + vec![cell_vertices[0], cell_vertices[1], cell_vertices[2]] } else { - vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], cell_vertices[3]] // Square base + vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], cell_vertices[3]] }; - // For a square-based pyramid, split into two tetrahedrons - if num_vertices == 5 { - let tetra1 = vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], apex]; - let tetra2 = vec![cell_vertices[0], cell_vertices[2], cell_vertices[3], apex]; - - // Compute volumes and centroids - let volume1 = self.compute_tetrahedron_volume(&tetra1); - let volume2 = self.compute_tetrahedron_volume(&tetra2); - - let centroid1 = self.compute_tetrahedron_centroid(&tetra1); - let centroid2 = self.compute_tetrahedron_centroid(&tetra2); - - _total_volume = volume1 + volume2; - - // Compute the weighted centroid - weighted_centroid[0] = (centroid1[0] * volume1 + centroid2[0] * volume2) / _total_volume; - weighted_centroid[1] = (centroid1[1] * volume1 + centroid2[1] * volume2) / _total_volume; - weighted_centroid[2] = (centroid1[2] * volume1 + centroid2[2] * volume2) / _total_volume; - } else { - // For triangular base pyramid (tetrahedron) - let tetra = vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], apex]; - let volume = self.compute_tetrahedron_volume(&tetra); - let centroid = self.compute_tetrahedron_centroid(&tetra); - - _total_volume = volume; - weighted_centroid = centroid; - } - - // Adjust the centroid between the base centroid and the apex let base_centroid = self.compute_face_centroid( - if num_vertices == 4 { FaceShape::Triangle } else { FaceShape::Quadrilateral }, + if cell_vertices.len() == 4 { FaceShape::Triangle } else { FaceShape::Quadrilateral }, &base_vertices, ); - // Apply pyramid centroid formula: 3/4 base_centroid + 1/4 apex - weighted_centroid[0] = (3.0 * base_centroid[0] + apex[0]) / 4.0; - weighted_centroid[1] = (3.0 * base_centroid[1] + apex[1]) / 4.0; - weighted_centroid[2] = (3.0 * base_centroid[2] + apex[2]) / 4.0; - - weighted_centroid + // Apply pyramid centroid formula: 3/4 * base_centroid + 1/4 * apex + [ + (3.0 * base_centroid[0] + apex[0]) / 4.0, + (3.0 * base_centroid[1] + apex[1]) / 4.0, + (3.0 * base_centroid[2] + apex[2]) / 4.0, + ] } /// Computes the volume of a pyramid cell (triangular or square base). @@ -101,83 +52,26 @@ impl Geometry { /// /// # Panics /// This function will panic if the number of vertices is not 4 (triangular) or 5 (square). - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let pyramid_vertices = vec![ - /// [0.0, 0.0, 0.0], // base vertex 1 - /// [1.0, 0.0, 0.0], // base vertex 2 - /// [1.0, 1.0, 0.0], // base vertex 3 - /// [0.0, 1.0, 0.0], // base vertex 4 - /// [0.5, 0.5, 1.0], // apex - /// ]; - /// let volume = geometry.compute_pyramid_volume(&pyramid_vertices); - /// assert!((volume - (1.0 / 3.0)).abs() < 1e-10); // Volume of a square pyramid - /// ``` pub fn compute_pyramid_volume(&self, cell_vertices: &Vec<[f64; 3]>) -> f64 { - let mut _total_volume = 0.0; - - let num_vertices = cell_vertices.len(); - assert!(num_vertices == 4 || num_vertices == 5, "Pyramid must have 4 (triangular) or 5 (square) vertices"); - - // Define the apex and base vertices - let apex = if num_vertices == 4 { cell_vertices[3] } else { cell_vertices[4] }; - let _base_vertices = if num_vertices == 4 { - vec![cell_vertices[0], cell_vertices[1], cell_vertices[2]] // Triangular base - } else { - vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], cell_vertices[3]] // Square base - }; + assert!(cell_vertices.len() == 4 || cell_vertices.len() == 5, "Pyramid must have 4 or 5 vertices"); - // Split into two tetrahedrons if square base - if num_vertices == 5 { + let apex = if cell_vertices.len() == 4 { cell_vertices[3] } else { cell_vertices[4] }; + + if cell_vertices.len() == 5 { let tetra1 = vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], apex]; let tetra2 = vec![cell_vertices[0], cell_vertices[2], cell_vertices[3], apex]; - let volume1 = self.compute_tetrahedron_volume_local(&tetra1); - let volume2 = self.compute_tetrahedron_volume_local(&tetra2); + let volume1 = self.compute_tetrahedron_volume(&tetra1); + let volume2 = self.compute_tetrahedron_volume(&tetra2); - _total_volume = volume1 + volume2; + volume1 + volume2 } else { let tetra = vec![cell_vertices[0], cell_vertices[1], cell_vertices[2], apex]; - _total_volume = self.compute_tetrahedron_volume_local(&tetra); + self.compute_tetrahedron_volume(&tetra) } - - _total_volume - } - - /// Computes the volume of a tetrahedron using its four vertices. - /// - /// # Arguments - /// * `vertices` - A vector of 4 vertices representing the tetrahedron. - /// - /// # Returns - /// * `f64` - The volume of the tetrahedron. - fn compute_tetrahedron_volume_local(&self, vertices: &Vec<[f64; 3]>) -> f64 { - assert!(vertices.len() == 4, "Tetrahedron must have 4 vertices"); - - let a = vertices[0]; - let b = vertices[1]; - let c = vertices[2]; - let d = vertices[3]; - - let v1 = [b[0] - d[0], b[1] - d[1], b[2] - d[2]]; - let v2 = [c[0] - d[0], c[1] - d[1], c[2] - d[2]]; - let v3 = [a[0] - d[0], a[1] - d[1], a[2] - d[2]]; - - let cross_product = [ - v2[1] * v3[2] - v2[2] * v3[1], - v2[2] * v3[0] - v2[0] * v3[2], - v2[0] * v3[1] - v2[1] * v3[0], - ]; - - let dot_product = v1[0] * cross_product[0] + v1[1] * cross_product[1] + v1[2] * cross_product[2]; - (dot_product.abs() / 6.0).abs() } } - #[cfg(test)] mod tests { use crate::geometry::Geometry; @@ -197,8 +91,7 @@ mod tests { let volume = geometry.compute_pyramid_volume(&pyramid_vertices); - // The expected volume of the square pyramid is (1/3) * base area * height - // Base area = 1.0, height = 1.0 -> Volume = (1/3) * 1.0 * 1.0 = 1/3 + // Expected volume: (1/3) * base area * height, where base area = 1.0, height = 1.0 assert!((volume - (1.0 / 3.0)).abs() < 1e-10, "Volume of the square pyramid is incorrect"); } @@ -216,8 +109,7 @@ mod tests { let volume = geometry.compute_pyramid_volume(&pyramid_vertices); - // The expected volume of the triangular pyramid (tetrahedron) is (1/3) * base area * height - // Base area = 0.5, height = 1.0 -> Volume = (1/3) * 0.5 * 1.0 = 1/6 + // Expected volume: (1/3) * base area * height, where base area = 0.5, height = 1.0 assert!((volume - (1.0 / 6.0)).abs() < 1e-10, "Volume of the triangular pyramid is incorrect"); } @@ -236,7 +128,7 @@ mod tests { let centroid = geometry.compute_pyramid_centroid(&pyramid_vertices); - // The correct centroid is at (0.5, 0.5, 0.25) + // Expected centroid: (0.5, 0.5, 0.25) assert_eq!(centroid, [0.5, 0.5, 0.25]); } @@ -254,7 +146,7 @@ mod tests { let centroid = geometry.compute_pyramid_centroid(&pyramid_vertices); - // The correct centroid is at (0.375, 0.375, 0.25) + // Expected centroid: (0.375, 0.375, 0.25) assert_eq!(centroid, [0.375, 0.375, 0.25]); } @@ -272,7 +164,7 @@ mod tests { let volume = geometry.compute_pyramid_volume(°enerate_pyramid_vertices); - // The volume of a degenerate pyramid should be zero + // Volume of a degenerate pyramid should be zero assert_eq!(volume, 0.0); } @@ -290,7 +182,7 @@ mod tests { let centroid = geometry.compute_pyramid_centroid(°enerate_pyramid_vertices); - // The centroid of this degenerate pyramid should still be computed correctly + // Expected centroid in the plane: (0.375, 0.375, 0.0) let expected_centroid = [0.375, 0.375, 0.0]; assert_eq!(centroid, expected_centroid); } diff --git a/src/geometry/quadrilateral.rs b/src/geometry/quadrilateral.rs index 74c8279..a080610 100644 --- a/src/geometry/quadrilateral.rs +++ b/src/geometry/quadrilateral.rs @@ -18,20 +18,6 @@ impl Geometry { /// # Panics /// /// This function will panic if the number of vertices provided is not exactly four. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let quad_vertices = vec![ - /// [0.0, 0.0, 0.0], - /// [1.0, 0.0, 0.0], - /// [1.0, 1.0, 0.0], - /// [0.0, 1.0, 0.0] - /// ]; - /// let area = geometry.compute_quadrilateral_area(&quad_vertices); - /// assert_eq!(area, 1.0); // The area of a square with side length 1 is 1.0 - /// ``` pub fn compute_quadrilateral_area(&self, quad_vertices: &Vec<[f64; 3]>) -> f64 { assert!(quad_vertices.len() == 4, "Quadrilateral must have exactly 4 vertices"); @@ -63,20 +49,6 @@ impl Geometry { /// # Panics /// /// This function will panic if the number of vertices provided is not exactly four. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let quad_vertices = vec![ - /// [0.0, 0.0, 0.0], - /// [1.0, 0.0, 0.0], - /// [1.0, 1.0, 0.0], - /// [0.0, 1.0, 0.0] - /// ]; - /// let centroid = geometry.compute_quadrilateral_centroid(&quad_vertices); - /// assert_eq!(centroid, [0.5, 0.5, 0.0]); // The centroid of the square is its center - /// ``` pub fn compute_quadrilateral_centroid(&self, quad_vertices: &Vec<[f64; 3]>) -> [f64; 3] { assert!(quad_vertices.len() == 4, "Quadrilateral must have exactly 4 vertices"); @@ -94,12 +66,58 @@ impl Geometry { centroid } + + /// Computes the normal vector for a quadrilateral face in 3D space. + /// + /// The normal vector is approximated by dividing the quadrilateral into two triangles, + /// calculating each triangle's normal, and averaging them. + /// + /// # Arguments + /// * `quad_vertices` - A vector of 3D coordinates representing the vertices of the quadrilateral. + /// + /// # Returns + /// * `[f64; 3]` - The normal vector for the quadrilateral face. + pub fn compute_quadrilateral_normal(&self, quad_vertices: &Vec<[f64; 3]>) -> [f64; 3] { + assert!(quad_vertices.len() == 4, "Quadrilateral must have exactly 4 vertices"); + + // Divide the quadrilateral into two triangles + let triangle1 = vec![quad_vertices[0], quad_vertices[1], quad_vertices[2]]; + let triangle2 = vec![quad_vertices[2], quad_vertices[3], quad_vertices[0]]; + + // Compute normals for each triangle + let normal1 = self.compute_triangle_normal(&triangle1); + let normal2 = self.compute_triangle_normal(&triangle2); + + // Average the two normals + [ + (normal1[0] + normal2[0]) / 2.0, + (normal1[1] + normal2[1]) / 2.0, + (normal1[2] + normal2[2]) / 2.0, + ] + } } #[cfg(test)] mod tests { use crate::geometry::Geometry; + #[test] + fn test_quadrilateral_normal() { + let geometry = Geometry::new(); + + let quad_vertices = vec![ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0], + ]; + + let normal = geometry.compute_quadrilateral_normal(&quad_vertices); + + // Expected normal should be [0.0, 0.0, 1.0] for this quadrilateral in the XY plane + assert!((normal[2] - 1.0).abs() < 1e-10); + } + #[test] fn test_quadrilateral_area_square() { let geometry = Geometry::new(); @@ -150,7 +168,7 @@ mod tests { let centroid = geometry.compute_quadrilateral_centroid(&quad_vertices); - // The centroid of a square centered at (0.5, 0.5, 0.0) + // The centroid of a square is at (0.5, 0.5, 0.0) assert_eq!(centroid, [0.5, 0.5, 0.0]); } diff --git a/src/geometry/tetrahedron.rs b/src/geometry/tetrahedron.rs index 589d2c9..64d2916 100644 --- a/src/geometry/tetrahedron.rs +++ b/src/geometry/tetrahedron.rs @@ -5,7 +5,7 @@ impl Geometry { /// /// This function calculates the centroid (geometric center) of a tetrahedron using /// the 3D coordinates of its four vertices. The centroid is the average of the - /// positions of all the vertices. + /// positions of all vertices. /// /// # Arguments /// @@ -20,25 +20,29 @@ impl Geometry { /// ```rust,ignore /// let geometry = Geometry::new(); /// let vertices = vec![ - /// [0.0, 0.0, 0.0], - /// [1.0, 0.0, 0.0], - /// [0.0, 1.0, 0.0], + /// [0.0, 0.0, 0.0], + /// [1.0, 0.0, 0.0], + /// [0.0, 1.0, 0.0], /// [0.0, 0.0, 1.0] /// ]; /// let centroid = geometry.compute_tetrahedron_centroid(&vertices); /// assert_eq!(centroid, [0.25, 0.25, 0.25]); /// ``` pub fn compute_tetrahedron_centroid(&self, cell_vertices: &Vec<[f64; 3]>) -> [f64; 3] { + assert_eq!(cell_vertices.len(), 4, "Tetrahedron must have exactly 4 vertices"); + let mut centroid = [0.0, 0.0, 0.0]; for v in cell_vertices { centroid[0] += v[0]; centroid[1] += v[1]; centroid[2] += v[2]; } - let num_vertices = cell_vertices.len() as f64; - centroid[0] /= num_vertices; - centroid[1] /= num_vertices; - centroid[2] /= num_vertices; + + // Average the vertex coordinates to get the centroid + for i in 0..3 { + centroid[i] /= 4.0; + } + centroid } @@ -65,16 +69,16 @@ impl Geometry { /// ```rust,ignore /// let geometry = Geometry::new(); /// let vertices = vec![ - /// [0.0, 0.0, 0.0], - /// [1.0, 0.0, 0.0], - /// [0.0, 1.0, 0.0], + /// [0.0, 0.0, 0.0], + /// [1.0, 0.0, 0.0], + /// [0.0, 1.0, 0.0], /// [0.0, 0.0, 1.0] /// ]; /// let volume = geometry.compute_tetrahedron_volume(&vertices); /// assert!((volume - 1.0 / 6.0).abs() < 1e-10); /// ``` pub fn compute_tetrahedron_volume(&self, tet_vertices: &Vec<[f64; 3]>) -> f64 { - assert!(tet_vertices.len() == 4, "Tetrahedron must have exactly 4 vertices"); + assert_eq!(tet_vertices.len(), 4, "Tetrahedron must have exactly 4 vertices"); let v0 = tet_vertices[0]; let v1 = tet_vertices[1]; diff --git a/src/geometry/triangle.rs b/src/geometry/triangle.rs index 84dd1f0..31c22b1 100644 --- a/src/geometry/triangle.rs +++ b/src/geometry/triangle.rs @@ -17,19 +17,9 @@ impl Geometry { /// # Panics /// /// This function will panic if the number of vertices provided is not exactly three. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let vertices = vec![[0.0, 0.0, 0.0], [3.0, 0.0, 0.0], [0.0, 4.0, 0.0]]; - /// let centroid = geometry.compute_triangle_centroid(&vertices); - /// assert_eq!(centroid, [1.0, 4.0 / 3.0, 0.0]); - /// ``` pub fn compute_triangle_centroid(&self, triangle_vertices: &Vec<[f64; 3]>) -> [f64; 3] { assert!(triangle_vertices.len() == 3, "Triangle must have exactly 3 vertices"); - // Sequential centroid calculation let mut centroid = [0.0, 0.0, 0.0]; for vertex in triangle_vertices { centroid[0] += vertex[0]; @@ -58,15 +48,6 @@ impl Geometry { /// # Panics /// /// This function will panic if the number of vertices provided is not exactly three. - /// - /// # Example - /// - /// ```rust,ignore - /// let geometry = Geometry::new(); - /// let vertices = vec![[0.0, 0.0, 0.0], [3.0, 0.0, 0.0], [0.0, 4.0, 0.0]]; - /// let area = geometry.compute_triangle_area(&vertices); - /// assert!((area - 6.0).abs() < 1e-10); // Area should be approximately 6.0 - /// ``` pub fn compute_triangle_area(&self, triangle_vertices: &Vec<[f64; 3]>) -> f64 { assert!( triangle_vertices.len() == 3, @@ -97,12 +78,59 @@ impl Geometry { // Area is half the magnitude of the cross product 0.5 * cross_product_magnitude } + + /// Computes the normal vector for a triangular face in 3D space. + /// + /// The normal vector is computed using the cross product of two edges from the triangle. + /// The length of the normal vector will be proportional to the area of the triangle. + /// + /// # Arguments + /// * `triangle_vertices` - A vector of 3D coordinates representing the vertices of the triangle. + /// + /// # Returns + /// * `[f64; 3]` - The normal vector for the triangular face. + pub fn compute_triangle_normal(&self, triangle_vertices: &Vec<[f64; 3]>) -> [f64; 3] { + assert!(triangle_vertices.len() == 3, "Triangle must have exactly 3 vertices"); + + let v0 = triangle_vertices[0]; + let v1 = triangle_vertices[1]; + let v2 = triangle_vertices[2]; + + // Compute edge vectors + let e1 = [v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]]; + let e2 = [v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]]; + + // Compute the cross product to get the normal vector + let normal = [ + e1[1] * e2[2] - e1[2] * e2[1], + e1[2] * e2[0] - e1[0] * e2[2], + e1[0] * e2[1] - e1[1] * e2[0], + ]; + + normal + } } #[cfg(test)] mod tests { use crate::geometry::Geometry; + #[test] + fn test_triangle_normal() { + let geometry = Geometry::new(); + + let triangle_vertices = vec![ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + ]; + + let normal = geometry.compute_triangle_normal(&triangle_vertices); + + // Expected normal should be [0.0, 0.0, 1.0] for this triangle in the XY plane + assert!((normal[2] - 1.0).abs() < 1e-10); + } + #[test] fn test_triangle_area() { let geometry = Geometry::new(); @@ -134,7 +162,7 @@ mod tests { let centroid = geometry.compute_triangle_centroid(&triangle_vertices); // The centroid of this triangle is the average of the vertices: - // ([0.0, 0.0, 0.0] + [3.0, 0.0, 0.0] + [0.0, 4.0, 0.0]) / 3 = [1.0, 1.3333, 0.0] + // ([0.0, 0.0, 0.0] + [3.0, 0.0, 0.0] + [0.0, 4.0, 0.0]) / 3 = [1.0, 4.0 / 3.0, 0.0] assert_eq!(centroid, [1.0, 4.0 / 3.0, 0.0]); } diff --git a/target/.rustc_info.json b/target/.rustc_info.json index 4f47a19..e380f6f 100644 --- a/target/.rustc_info.json +++ b/target/.rustc_info.json @@ -1 +1 @@ -{"rustc_fingerprint":18432001303766494279,"outputs":{"14305839045545480886":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\n","stderr":""},"15729799797837862367":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\nC:\\Users\\Tea\\.rustup\\toolchains\\stable-x86_64-pc-windows-msvc\npacked\n___\ndebug_assertions\noverflow_checks\npanic=\"unwind\"\nproc_macro\nrelocation_model=\"pic\"\ntarget_abi=\"\"\ntarget_arch=\"x86_64\"\ntarget_endian=\"little\"\ntarget_env=\"msvc\"\ntarget_family=\"windows\"\ntarget_feature=\"cmpxchg16b\"\ntarget_feature=\"fxsr\"\ntarget_feature=\"lahfsahf\"\ntarget_feature=\"sse\"\ntarget_feature=\"sse2\"\ntarget_feature=\"sse3\"\ntarget_has_atomic\ntarget_has_atomic=\"128\"\ntarget_has_atomic=\"16\"\ntarget_has_atomic=\"32\"\ntarget_has_atomic=\"64\"\ntarget_has_atomic=\"8\"\ntarget_has_atomic=\"ptr\"\ntarget_has_atomic_equal_alignment=\"128\"\ntarget_has_atomic_equal_alignment=\"16\"\ntarget_has_atomic_equal_alignment=\"32\"\ntarget_has_atomic_equal_alignment=\"64\"\ntarget_has_atomic_equal_alignment=\"8\"\ntarget_has_atomic_equal_alignment=\"ptr\"\ntarget_has_atomic_load_store\ntarget_has_atomic_load_store=\"128\"\ntarget_has_atomic_load_store=\"16\"\ntarget_has_atomic_load_store=\"32\"\ntarget_has_atomic_load_store=\"64\"\ntarget_has_atomic_load_store=\"8\"\ntarget_has_atomic_load_store=\"ptr\"\ntarget_os=\"windows\"\ntarget_pointer_width=\"64\"\ntarget_thread_local\ntarget_vendor=\"pc\"\nub_checks\nwindows\n","stderr":""},"12744816824612481171":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\nC:\\Users\\Tea\\.rustup\\toolchains\\stable-x86_64-pc-windows-msvc\npacked\n___\ndebug_assertions\npanic=\"unwind\"\nproc_macro\ntarget_abi=\"\"\ntarget_arch=\"x86_64\"\ntarget_endian=\"little\"\ntarget_env=\"msvc\"\ntarget_family=\"windows\"\ntarget_feature=\"cmpxchg16b\"\ntarget_feature=\"fxsr\"\ntarget_feature=\"sse\"\ntarget_feature=\"sse2\"\ntarget_feature=\"sse3\"\ntarget_has_atomic=\"128\"\ntarget_has_atomic=\"16\"\ntarget_has_atomic=\"32\"\ntarget_has_atomic=\"64\"\ntarget_has_atomic=\"8\"\ntarget_has_atomic=\"ptr\"\ntarget_os=\"windows\"\ntarget_pointer_width=\"64\"\ntarget_vendor=\"pc\"\nwindows\n","stderr":""},"4614504638168534921":{"success":true,"status":"","code":0,"stdout":"rustc 1.80.1 (3f5fd8dd4 2024-08-06)\nbinary: rustc\ncommit-hash: 3f5fd8dd41153bc5fdca9427e9e05be2c767ba23\ncommit-date: 2024-08-06\nhost: x86_64-pc-windows-msvc\nrelease: 1.80.1\nLLVM version: 18.1.7\n","stderr":""},"16495917692426387086":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\n","stderr":""}},"successes":{}} \ No newline at end of file +{"rustc_fingerprint":15376570777996937782,"outputs":{"4614504638168534921":{"success":true,"status":"","code":0,"stdout":"rustc 1.80.1 (3f5fd8dd4 2024-08-06)\nbinary: rustc\ncommit-hash: 3f5fd8dd41153bc5fdca9427e9e05be2c767ba23\ncommit-date: 2024-08-06\nhost: x86_64-pc-windows-msvc\nrelease: 1.80.1\nLLVM version: 18.1.7\n","stderr":""},"15729799797837862367":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\nC:\\Users\\Tea\\.rustup\\toolchains\\stable-x86_64-pc-windows-msvc\npacked\n___\ndebug_assertions\npanic=\"unwind\"\nproc_macro\ntarget_abi=\"\"\ntarget_arch=\"x86_64\"\ntarget_endian=\"little\"\ntarget_env=\"msvc\"\ntarget_family=\"windows\"\ntarget_feature=\"cmpxchg16b\"\ntarget_feature=\"fxsr\"\ntarget_feature=\"sse\"\ntarget_feature=\"sse2\"\ntarget_feature=\"sse3\"\ntarget_has_atomic=\"128\"\ntarget_has_atomic=\"16\"\ntarget_has_atomic=\"32\"\ntarget_has_atomic=\"64\"\ntarget_has_atomic=\"8\"\ntarget_has_atomic=\"ptr\"\ntarget_os=\"windows\"\ntarget_pointer_width=\"64\"\ntarget_vendor=\"pc\"\nwindows\n","stderr":""},"16495917692426387086":{"success":true,"status":"","code":0,"stdout":"___.exe\nlib___.rlib\n___.dll\n___.dll\n___.lib\n___.dll\n","stderr":""}},"successes":{}} \ No newline at end of file