Implemented parallel communication with rust concurrency primitives a…

…nd crossbeam for domain module.

Cargo.toml

@@ -11,6 +11,7 @@ path = "src/lib.rs"
 rustc-hash = "2.0.0"
 faer = { version = "0.19", features = ["rayon"] }
 rayon = "1.5"
+crossbeam = "0.8.4"
 
 [dev-dependencies]
 criterion = "0.4"

src/boundary/bc_handler.rs

@@ -1,21 +1,23 @@
 // src/boundary/bc_handler.rs
 
 use rustc_hash::FxHashMap;
+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;  // Import NeumannBC
+use crate::boundary::neumann::NeumannBC;
 use crate::boundary::robin::RobinBC;
 use faer::MatMut;
 
-pub type BoundaryConditionFn = Box<dyn Fn(f64, &[f64]) -> f64 + Send + Sync>;
+pub type BoundaryConditionFn = Arc<dyn Fn(f64, &[f64]) -> f64 + Send + Sync>;
 
+#[derive(Clone)]
 pub enum BoundaryCondition {
     Dirichlet(f64),
     Neumann(f64),
     Robin { alpha: f64, beta: f64 },
-    DirichletFn(BoundaryConditionFn), // Dirichlet with a function
-    NeumannFn(BoundaryConditionFn),   // Neumann with a function
+    DirichletFn(BoundaryConditionFn),
+    NeumannFn(BoundaryConditionFn),
 }
 
 pub struct BoundaryConditionHandler {
@@ -29,43 +31,72 @@ impl BoundaryConditionHandler {
         }
     }
 
-    pub fn set_bc(&mut self, entity: MeshEntity, condition: BoundaryCondition) {
+    pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) {
         self.conditions.set_data(entity, condition);
     }
 
-    pub fn get_bc(&self, entity: &MeshEntity) -> Option<&BoundaryCondition> {
+    pub fn get_bc(&self, entity: &MeshEntity) -> Option<BoundaryCondition> {
         self.conditions.restrict(entity)
     }
 
-    pub fn apply_bc(&self, matrix: &mut MatMut<f64>, rhs: &mut MatMut<f64>, boundary_entities: &[MeshEntity], entity_to_index: &FxHashMap<MeshEntity, usize>, time: f64) {
+    pub fn apply_bc(
+        &self,
+        matrix: &mut MatMut<f64>,
+        rhs: &mut MatMut<f64>,
+        boundary_entities: &[MeshEntity],
+        entity_to_index: &FxHashMap<MeshEntity, usize>,
+        time: f64,
+    ) {
         for entity in boundary_entities {
             if let Some(bc) = self.get_bc(entity) {
                 match bc {
                     BoundaryCondition::Dirichlet(value) => {
-                        // Apply Dirichlet boundary condition
                         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,
+                            *entity_to_index.get(entity).unwrap(),
+                            value,
+                        );
                     }
                     BoundaryCondition::Neumann(flux) => {
-                        // Apply Neumann boundary condition
                         let neumann_bc = NeumannBC::new();
-                        neumann_bc.apply_constant_neumann(rhs, *entity_to_index.get(entity).unwrap(), *flux);
+                        neumann_bc.apply_constant_neumann(
+                            rhs,
+                            *entity_to_index.get(entity).unwrap(),
+                            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,
+                            *entity_to_index.get(entity).unwrap(),
+                            alpha,
+                            beta,
+                        );
                     }
                     BoundaryCondition::DirichletFn(fn_bc) => {
-                        let coords = [0.0, 0.0, 0.0];  // Placeholder: 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,
+                            *entity_to_index.get(entity).unwrap(),
+                            value,
+                        );
                     }
                     BoundaryCondition::NeumannFn(fn_bc) => {
-                        let coords = [0.0, 0.0, 0.0];  // Placeholder: 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,
+                            *entity_to_index.get(entity).unwrap(),
+                            value,
+                        );
                     }
                 }
             }

src/boundary/dirichlet.rs

@@ -7,7 +7,7 @@ use crate::domain::section::Section;
 use faer::MatMut;
 
 pub struct DirichletBC {
-    conditions: Section<BoundaryCondition>,  // Section to hold Dirichlet conditions
+    conditions: Section<BoundaryCondition>,
 }
 
 impl DirichletBC {
@@ -17,22 +17,26 @@ impl DirichletBC {
         }
     }
 
-    // Set a Dirichlet boundary condition for a specific entity
-    // Modify set_bc to accept BoundaryCondition instead of just f64
-    pub fn set_bc(&mut self, entity: MeshEntity, condition: BoundaryCondition) {
+    pub fn set_bc(&self, entity: MeshEntity, condition: BoundaryCondition) {
         self.conditions.set_data(entity, condition);
     }
 
-    // Apply the Dirichlet boundary condition during the system matrix assembly
-    pub fn apply_bc(&self, matrix: &mut MatMut<f64>, rhs: &mut MatMut<f64>, entity_to_index: &FxHashMap<MeshEntity, usize>, time: f64) {
-        for (entity, condition) in self.conditions.data.iter() {
+    pub fn apply_bc(
+        &self,
+        matrix: &mut MatMut<f64>,
+        rhs: &mut MatMut<f64>,
+        entity_to_index: &FxHashMap<MeshEntity, usize>,
+        time: f64,
+    ) {
+        let data = self.conditions.data.read().unwrap();
+        for (entity, condition) in data.iter() {
             if let Some(&index) = entity_to_index.get(entity) {
-                match condition {  // Dereference the condition
+                match condition {
                     BoundaryCondition::Dirichlet(value) => {
                         self.apply_constant_dirichlet(matrix, rhs, index, *value);
                     }
                     BoundaryCondition::DirichletFn(fn_bc) => {
-                        let coords = self.get_coordinates(entity);  // Placeholder for actual method
+                        let coords = self.get_coordinates(entity);
                         let value = fn_bc(time, &coords);
                         self.apply_constant_dirichlet(matrix, rhs, index, value);
                     }
@@ -42,28 +46,27 @@ impl DirichletBC {
         }
     }
 
-    pub fn apply_constant_dirichlet(&self, matrix: &mut MatMut<f64>, rhs: &mut MatMut<f64>, index: usize, value: f64) {
+    pub fn apply_constant_dirichlet(
+        &self,
+        matrix: &mut MatMut<f64>,
+        rhs: &mut MatMut<f64>,
+        index: usize,
+        value: f64,
+    ) {
         let ncols = matrix.ncols();
-
-        // Zero out the entire row for the Dirichlet condition
         for col in 0..ncols {
-            matrix.write(index, col, 0.0);  // Set each element in the row to 0
+            matrix.write(index, col, 0.0);
         }
-
-        // Set the diagonal element to 1 to maintain the boundary condition in the system
         matrix.write(index, index, 1.0);
-
-        // Set the corresponding value in the RHS vector
         rhs.write(index, 0, value);
     }
 
     fn get_coordinates(&self, _entity: &MeshEntity) -> [f64; 3] {
-        // Placeholder: Implement logic to retrieve the coordinates of the mesh entity.
-        // This method would retrieve the spatial coordinates associated with a vertex or entity.
-        [0.0, 0.0, 0.0]  // Example placeholder return value
+        [0.0, 0.0, 0.0]
     }
 }
 
+
 impl BoundaryConditionApply for DirichletBC {
     fn apply(&self, _entity: &MeshEntity, rhs: &mut MatMut<f64>, matrix: &mut MatMut<f64>, entity_to_index: &FxHashMap<MeshEntity, usize>, time: f64) {
         // Dirichlet-specific logic
@@ -77,14 +80,11 @@ mod tests {
     use rustc_hash::FxHashMap;
     use faer::Mat;
     use crate::domain::mesh_entity::MeshEntity;
+    use std::sync::Arc;
 
     fn create_test_matrix_and_rhs() -> (Mat<f64>, Mat<f64>) {
-        // 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)
     }
 
@@ -134,32 +134,30 @@ mod tests {
 
     #[test]
     fn test_apply_function_based_dirichlet() {
-        let mut dirichlet_bc = DirichletBC::new();
+        let dirichlet_bc = DirichletBC::new();
         let entity = MeshEntity::Vertex(2);
         let mut entity_to_index = FxHashMap::default();
         entity_to_index.insert(entity, 2);
 
-        // Set a function-based Dirichlet boundary condition
-        dirichlet_bc.set_bc(entity, BoundaryCondition::DirichletFn(Box::new(|_time: f64, _coords: &[f64]| 7.0)));
+        dirichlet_bc.set_bc(
+            entity,
+            BoundaryCondition::DirichletFn(Arc::new(|_time: f64, _coords: &[f64]| 7.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 function-based Dirichlet condition
         dirichlet_bc.apply_bc(&mut matrix_mut, &mut rhs_mut, &entity_to_index, 1.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 == 2 {
-                assert_eq!(matrix_mut[(2, col)], 1.0); // Diagonal should be 1
+                assert_eq!(matrix_mut[(2, col)], 1.0);
             } else {
-                assert_eq!(matrix_mut[(2, col)], 0.0); // Off-diagonal should be 0
+                assert_eq!(matrix_mut[(2, col)], 0.0);
             }
         }
 
-        // Check that the RHS has the correct value based on the function
         assert_eq!(rhs_mut[(2, 0)], 7.0);
     }
 }