Integrating first Laplace evaluator

benches/assembly_benchmark.rs

@@ -2,6 +2,7 @@ use bempp::boundary_assemblers::BoundaryAssemblerOptions;
 use bempp::function::FunctionSpace;
 use bempp::function::LocalFunctionSpaceTrait;
 use bempp::laplace::assembler::single_layer;
+use bempp_distributed_tools::SingleProcessIndexLayout;
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use ndelement::ciarlet::LagrangeElementFamily;
 use ndelement::types::{Continuity, ReferenceCellType};
@@ -27,14 +28,24 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
         options.set_batch_size(128);
 
         let assembler = single_layer(&options);
+        let index_layout = SingleProcessIndexLayout::new(0, space.global_size(), &comm);
 
         group.bench_function(
             format!(
                 "Assembly of singular terms of {}x{} matrix",
                 space.global_size(),
                 space.global_size()
             ),
-            |b| b.iter(|| black_box(assembler.assemble_singular(&space, &space))),
+            |b| {
+                b.iter(|| {
+                    black_box(assembler.assemble_singular(
+                        &space,
+                        &index_layout,
+                        &space,
+                        &index_layout,
+                    ))
+                })
+            },
         );
     }
     group.finish();

examples/laplace_evaluator.rs

@@ -0,0 +1,64 @@
+//! This file implements an example Laplace evaluator test the different involved operators.
+
+use bempp::{boundary_assemblers::BoundaryAssemblerOptions, function::LocalFunctionSpaceTrait};
+use bempp_quadrature::types::NumericalQuadratureGenerator;
+use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
+use ndgrid::traits::Grid;
+use rlst::operator::interface::DistributedArrayVectorSpace;
+
+fn main() {
+    let universe = mpi::initialize().unwrap();
+    let world = universe.world();
+
+    let grid = bempp::shapes::regular_sphere::<f64, _>(3, 1, &world);
+
+    // Get the number of cells in the grid.
+
+    let n_cells = grid.entity_iter(2).count();
+
+    let space = bempp::function::FunctionSpace::new(
+        &grid,
+        &LagrangeElementFamily::<f64>::new(1, ndelement::types::Continuity::Standard),
+    );
+
+    let mut options = BoundaryAssemblerOptions::default();
+    options.set_regular_quadrature_degree(ReferenceCellType::Triangle, 6);
+
+    let quad_degree = options
+        .get_regular_quadrature_degree(ReferenceCellType::Triangle)
+        .unwrap();
+
+    let assembler = bempp::laplace::assembler::single_layer::<f64>(&options);
+
+    let dense_matrix = assembler.assemble(&space, &space);
+
+    // Now let's build an evaluator.
+
+    //First initialise the index layouts.
+
+    let space_layout =
+        bempp_distributed_tools::SingleProcessIndexLayout::new(0, space.global_size(), &world);
+
+    let point_layout = bempp_distributed_tools::SingleProcessIndexLayout::new(
+        0,
+        quad_degree * space.global_size(),
+        &world,
+    );
+
+    // Instantiate function spaces.
+
+    let array_function_space = DistributedArrayVectorSpace::<_, f64>::new(&space_layout);
+    let point_function_space = DistributedArrayVectorSpace::<_, f64>::new(&point_layout);
+
+    let qrule = bempp_quadrature::simplex_rules::simplex_rule_triangle(6).unwrap();
+
+    println!("Here");
+    let space_to_point = bempp::evaluator_tools::basis_to_point_map(
+        &space,
+        &array_function_space,
+        &point_function_space,
+        &qrule.points,
+        &qrule.weights,
+        false,
+    );
+}

src/boundary_assemblers.rs

@@ -9,14 +9,15 @@ use crate::boundary_assemblers::cell_pair_assemblers::{
 use crate::boundary_assemblers::helpers::KernelEvaluator;
 use crate::boundary_assemblers::helpers::{equal_grids, RawData2D, RlstArray, SparseMatrixData};
 use crate::function::{FunctionSpaceTrait, LocalFunctionSpaceTrait};
+use bempp_distributed_tools::index_layout;
 use bempp_quadrature::duffy::{
     quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy,
 };
 use bempp_quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition};
 use green_kernels::traits::Kernel;
 use integrands::BoundaryIntegrand;
 use itertools::izip;
-use mpi::traits::Communicator;
+use mpi::traits::{Communicator, Equivalence};
 use ndelement::quadrature::simplex_rule;
 use ndelement::reference_cell;
 use ndelement::traits::FiniteElement;
@@ -25,8 +26,9 @@ use ndgrid::traits::{Entity, Grid, Topology};
 use ndgrid::types::Ownership;
 use rayon::prelude::*;
 use rlst::{
-    rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, DynamicArray,
-    MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape,
+    rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, DistributedCsrMatrix,
+    DynamicArray, IndexLayout, MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar,
+    Shape,
 };
 use std::collections::HashMap;
 
@@ -120,43 +122,46 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, K:
     BoundaryAssembler<'o, T, Integrand, K>
 {
     /// Assemble the singular part into a CSR matrix.
-    pub fn assemble_singular<Space: FunctionSpaceTrait<T = T>>(
+    pub fn assemble_singular<
+        'a,
+        C: Communicator,
+        TrialLayout: IndexLayout<Comm = C>,
+        TestLayout: IndexLayout<Comm = C>,
+        Space: FunctionSpaceTrait<T = T>,
+    >(
         &self,
         trial_space: &Space,
+        trial_index_layout: &'a TrialLayout,
         test_space: &Space,
-    ) -> CsrMatrix<T>
+        test_index_layout: &'a TestLayout,
+    ) -> DistributedCsrMatrix<'a, TrialLayout, TestLayout, T, C>
     where
         Space::LocalFunctionSpace: Sync,
+        T: Equivalence,
     {
+        assert_eq!(
+            trial_space.global_size(),
+            trial_index_layout.number_of_global_indices()
+        );
+
+        assert_eq!(
+            test_space.global_size(),
+            test_index_layout.number_of_global_indices()
+        );
+
         let shape = [test_space.global_size(), trial_space.global_size()];
         let sparse_matrix =
             self.assemble_singular_part(shape, trial_space.local_space(), test_space.local_space());
 
-        if sparse_matrix.data.is_empty()
-            || sparse_matrix
-                .data
-                .iter()
-                .map(|i| i.abs())
-                .filter(|i| *i > T::from(1e-10).unwrap().re())
-                .count()
-                == 0
-        {
-            // TODO: remove this hack once https://github.com/linalg-rs/rlst/pull/100 is merged and there a new release of RLST
-            CsrMatrix::<T>::new(
-                sparse_matrix.shape,
-                vec![],
-                vec![0; sparse_matrix.shape[0] + 1],
-                vec![],
-            )
-        } else {
-            CsrMatrix::<T>::from_aij(
-                sparse_matrix.shape,
-                &sparse_matrix.rows,
-                &sparse_matrix.cols,
-                &sparse_matrix.data,
-            )
-            .unwrap()
-        }
+        // Instantiate the CSR matrix.
+
+        DistributedCsrMatrix::from_aij(
+            trial_index_layout,
+            test_index_layout,
+            &sparse_matrix.rows,
+            &sparse_matrix.cols,
+            &sparse_matrix.data,
+        )
     }
 
     /// Assemble into a dense matrix.