First version of function evaluator

bempp · Jan 25, 2025 · 0a30d4a · 0a30d4a
1 parent 03caaef
commit 0a30d4a
Show file tree

Hide file tree

Showing 11 changed files with 446 additions and 335 deletions.
diff --git a/Cargo.toml b/Cargo.toml
@@ -27,7 +27,7 @@ crate-type = ["lib", "cdylib"]
 [dependencies]
 bempp-quadrature = { version = "0.1.0" }
 itertools = "0.13.*"
-mpi = { version = "0.8.*"}
+mpi = { version = "0.8.*", features = ["complex"]}
 num = "0.4"
 ndelement = { git="https://github.com/bempp/ndelement.git", features = ["mpi"]}
 # ndelement = { path = "../ndelement", features = ["mpi"] }

diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs
@@ -1,5 +1,6 @@
 use bempp::boundary_assemblers::BoundaryAssemblerOptions;
 use bempp::function::FunctionSpace;
+use bempp::function::FunctionSpaceTrait;
 use bempp::function::LocalFunctionSpaceTrait;
 use bempp::laplace::assembler::single_layer;
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
@@ -31,8 +32,8 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
         group.bench_function(
             format!(
                 "Assembly of singular terms of {}x{} matrix",
-                space.global_size(),
-                space.global_size()
+                space.local_space().global_size(),
+                space.local_space().global_size()
             ),
             |b| b.iter(|| black_box(assembler.assemble_singular(&space, &space))),
         );

diff --git a/examples/laplace_evaluator.rs b/examples/laplace_evaluator.rs
@@ -3,14 +3,15 @@
 use std::rc::Rc;
 
 use bempp::{
-    boundary_assemblers::BoundaryAssemblerOptions, evaluator_tools::NeighbourEvaluator,
-    function::LocalFunctionSpaceTrait,
+    boundary_assemblers::BoundaryAssemblerOptions,
+    evaluator_tools::NeighbourEvaluator,
+    function::{FunctionSpaceTrait, LocalFunctionSpaceTrait},
 };
 use green_kernels::laplace_3d::Laplace3dKernel;
 use mpi::traits::Communicator;
 use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
 use ndgrid::{
-    traits::{Entity, GeometryMap, Grid},
+    traits::{Entity, GeometryMap, Grid, ParallelGrid},
     types::Ownership,
 };
 use rand::SeedableRng;
@@ -31,7 +32,7 @@ fn main() {
     let quad_degree = 6;
     // Get the number of cells in the grid.
 
-    let n_cells = grid.entity_iter(2).count();
+    let n_cells = grid.local_grid().entity_iter(2).count();
 
     println!("Number of cells: {}", n_cells);
 
@@ -56,7 +57,7 @@ fn main() {
     //First initialise the index layouts.
 
     let space_layout = Rc::new(bempp_distributed_tools::IndexLayout::from_local_counts(
-        space.global_size(),
+        space.local_space().global_size(),
         &world,
     ));
 
@@ -84,9 +85,12 @@ fn main() {
 
     let mut points = vec![0 as f64; 3 * quad_degree * n_cells];
 
-    let geometry_map = grid.geometry_map(ReferenceCellType::Triangle, &qrule.points);
+    let geometry_map = grid
+        .local_grid()
+        .geometry_map(ReferenceCellType::Triangle, &qrule.points);
 
     for cell in grid
+        .local_grid()
         .entity_iter(2)
         .filter(|e| matches!(e.ownership(), Ownership::Owned))
     {
@@ -100,24 +104,24 @@ fn main() {
         );
     }
 
-    let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
-        &points,
-        &points,
-        green_kernels::types::GreenKernelEvalType::Value,
-        true,
-        Laplace3dKernel::default(),
-        &world,
-    );
-
-    // let kernel_evaluator = bempp::greens_function_evaluators::kifmm_evaluator::KiFmmEvaluator::new(
-    //     &points, &points, 1, 3, 5, &world,
+    // let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
+    //     &points,
+    //     &points,
+    //     green_kernels::types::GreenKernelEvalType::Value,
+    //     true,
+    //     Laplace3dKernel::default(),
+    //     &world,
     // );
 
+    let kernel_evaluator = bempp::greens_function_evaluators::kifmm_evaluator::KiFmmEvaluator::new(
+        &points, &points, 1, 3, 5, &world,
+    );
+
     let correction = NeighbourEvaluator::new(
         &qrule.points,
         Laplace3dKernel::default(),
         green_kernels::types::GreenKernelEvalType::Value,
-        space.support_cells(),
+        space.local_space().support_cells(),
         &grid,
     );
 
@@ -138,7 +142,7 @@ fn main() {
 
     let res_local = res.view().local();
 
-    let mut expected = rlst_dynamic_array1!(f64, [space.global_size()]);
+    let mut expected = rlst_dynamic_array1!(f64, [space.local_space().global_size()]);
 
     expected
         .r_mut()

diff --git a/examples/neighbour_evaluator.rs b/examples/neighbour_evaluator.rs
@@ -8,7 +8,7 @@ use itertools::Itertools;
 use mpi::traits::Communicator;
 use ndelement::types::ReferenceCellType;
 use ndgrid::{
-    traits::{Entity, GeometryMap, Grid},
+    traits::{Entity, GeometryMap, Grid, ParallelGrid},
     types::Ownership,
 };
 use rand::SeedableRng;
@@ -34,6 +34,7 @@ fn main() {
     // Now get the active cells on the current process.
 
     let n_cells = grid
+        .local_grid()
         .entity_iter(2)
         .filter(|e| matches!(e.ownership(), Ownership::Owned))
         .count();
@@ -46,7 +47,7 @@ fn main() {
         points.data(),
         Laplace3dKernel::default(),
         green_kernels::types::GreenKernelEvalType::Value,
-        &(0..grid.entity_count(ReferenceCellType::Triangle)).collect_vec(),
+        &(0..grid.local_grid().entity_count(ReferenceCellType::Triangle)).collect_vec(),
         &grid,
     );
 
@@ -56,9 +57,12 @@ fn main() {
 
     let mut physical_points = vec![0 as f64; 3 * n_points * n_cells];
 
-    let geometry_map = grid.geometry_map(ReferenceCellType::Triangle, points.data());
+    let geometry_map = grid
+        .local_grid()
+        .geometry_map(ReferenceCellType::Triangle, points.data());
 
     for cell in grid
+        .local_grid()
         .entity_iter(2)
         .filter(|e| matches!(e.ownership(), Ownership::Owned))
     {

diff --git a/examples/parallel_laplace_evaluator.rs b/examples/parallel_laplace_evaluator.rs
@@ -3,15 +3,16 @@
 use std::rc::Rc;
 
 use bempp::{
-    boundary_assemblers::BoundaryAssemblerOptions, evaluator_tools::NeighbourEvaluator,
-    function::LocalFunctionSpaceTrait,
+    boundary_assemblers::BoundaryAssemblerOptions,
+    evaluator_tools::NeighbourEvaluator,
+    function::{FunctionSpaceTrait, LocalFunctionSpaceTrait},
 };
 use green_kernels::laplace_3d::Laplace3dKernel;
 use itertools::Itertools;
 use mpi::traits::Communicator;
 use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
 use ndgrid::{
-    traits::{Entity, GeometryMap, Grid, Topology},
+    traits::{Entity, GeometryMap, Grid, ParallelGrid, Topology},
     types::Ownership,
 };
 use rand::SeedableRng;
@@ -32,7 +33,7 @@ fn main() {
     let quad_degree = 6;
     // Get the number of cells in the grid.
 
-    let n_cells = grid.entity_iter(2).count();
+    let n_cells = grid.local_grid().entity_iter(2).count();
 
     println!("Number of cells: {}", n_cells);
 
@@ -41,8 +42,9 @@ fn main() {
         &LagrangeElementFamily::<f64>::new(1, ndelement::types::Continuity::Discontinuous),
     );
 
-    let dofs = space.cell_dofs(50).unwrap().to_vec();
+    let dofs = space.local_space().cell_dofs(50).unwrap().to_vec();
     let vertex_ids = grid
+        .local_grid()
         .entity(2, 50)
         .unwrap()
         .topology()
@@ -68,7 +70,7 @@ fn main() {
     //First initialise the index layouts.
 
     let space_layout = Rc::new(bempp_distributed_tools::IndexLayout::from_local_counts(
-        space.global_size(),
+        space.local_space().global_size(),
         &world,
     ));
 
@@ -96,9 +98,12 @@ fn main() {
 
     let mut points = vec![0 as f64; 3 * quad_degree * n_cells];
 
-    let geometry_map = grid.geometry_map(ReferenceCellType::Triangle, &qrule.points);
+    let geometry_map = grid
+        .local_grid()
+        .geometry_map(ReferenceCellType::Triangle, &qrule.points);
 
     for cell in grid
+        .local_grid()
         .entity_iter(2)
         .filter(|e| matches!(e.ownership(), Ownership::Owned))
     {
@@ -129,7 +134,7 @@ fn main() {
         &qrule.points,
         Laplace3dKernel::default(),
         green_kernels::types::GreenKernelEvalType::Value,
-        space.support_cells(),
+        space.local_space().support_cells(),
         &grid,
     );
 
@@ -150,7 +155,7 @@ fn main() {
 
     let res_local = res.view().local();
 
-    let mut expected = rlst_dynamic_array1!(f64, [space.global_size()]);
+    let mut expected = rlst_dynamic_array1!(f64, [space.local_space().global_size()]);
 
     expected
         .r_mut()

diff --git a/src/boundary_assemblers.rs b/src/boundary_assemblers.rs
@@ -129,7 +129,10 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, K:
         Space::LocalFunctionSpace: Sync,
         T: Equivalence,
     {
-        let shape = [test_space.global_size(), trial_space.global_size()];
+        let shape = [
+            test_space.local_space().global_size(),
+            trial_space.local_space().global_size(),
+        ];
         let sparse_matrix =
             self.assemble_singular_part(shape, trial_space.local_space(), test_space.local_space());
 
@@ -157,8 +160,13 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, K:
             panic!("Dense assembly can only be used for function spaces stored in serial");
         }
 
-        let mut output =
-            rlst_dynamic_array2!(T, [test_space.global_size(), trial_space.global_size()]);
+        let mut output = rlst_dynamic_array2!(
+            T,
+            [
+                test_space.local_space().global_size(),
+                trial_space.local_space().global_size()
+            ]
+        );
 
         self.assemble_into_memory(trial_space, test_space, output.data_mut());
 
@@ -176,15 +184,18 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, K:
     {
         assert_eq!(
             output.len(),
-            test_space.global_size() * trial_space.global_size()
+            test_space.local_space().global_size() * trial_space.local_space().global_size()
         );
         if trial_space.comm().size() > 1 || test_space.comm().size() > 1 {
             panic!("Dense assembly can only be used for function spaces stored in serial");
         }
 
-        let test_colouring = test_space.cell_colouring();
-        let trial_colouring = trial_space.cell_colouring();
-        let shape = [test_space.global_size(), trial_space.global_size()];
+        let test_colouring = test_space.local_space().cell_colouring();
+        let trial_colouring = trial_space.local_space().cell_colouring();
+        let shape = [
+            test_space.local_space().global_size(),
+            trial_space.local_space().global_size(),
+        ];
         let output_raw = RawData2D {
             data: output.as_mut_ptr(),
             shape,