WIP: Refactor

src/boundary_assemblers.rs

@@ -121,36 +121,24 @@ 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<'a, C: Communicator, Space: FunctionSpaceTrait<T = T>>(
+    pub fn assemble_singular<'a, C: Communicator, Space: FunctionSpaceTrait<T = T, C = C>>(
         &self,
-        trial_space: &Space,
-        trial_index_layout: &'a IndexLayout<'a, C>,
-        test_space: &Space,
-        test_index_layout: &'a IndexLayout<'a, C>,
+        trial_space: &'a Space,
+        test_space: &'a Space,
     ) -> DistributedCsrMatrix<'a, 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());
 
         // Instantiate the CSR matrix.
 
         DistributedCsrMatrix::from_aij(
-            trial_index_layout,
-            test_index_layout,
+            trial_space.index_layout(),
+            test_space.index_layout(),
             &sparse_matrix.rows,
             &sparse_matrix.cols,
             &sparse_matrix.data,

src/evaluator_tools.rs

@@ -37,13 +37,8 @@ use crate::function::FunctionSpaceTrait;
 
 /// Create a linear operator from the map of a basis to points. The points are sorted by global
 /// index of the corresponding cells.
-pub fn basis_to_point_map<
-    'a,
-    C: Communicator,
-    T: RlstScalar + Equivalence,
-    Space: FunctionSpaceTrait<T = T>,
->(
-    function_space: &Space,
+pub fn basis_to_point_map<'a, T: RlstScalar + Equivalence, Space: FunctionSpaceTrait<T = T>>(
+    function_space: &'a Space,
     quadrature_points: &[T::Real],
     quadrature_weights: &[T::Real],
     return_transpose: bool,
@@ -64,7 +59,7 @@ where
 
     let reference_cell = grid.entity_types(tdim)[0];
 
-    let owned_active_cells = function_space
+    let owned_active_cells: Vec<usize> = function_space
         .support_cells()
         .iter()
         .filter(|index| {
@@ -73,9 +68,8 @@ where
                 Ownership::Owned
             )
         })
-        .sorted_by_key(|index| grid.entity(tdim, *index).unwrap().global_index())
-        .map(|index| grid.entity(tdim, *index).unwrap().global_index())
         .copied()
+        .sorted_by_key(|&index| grid.entity(tdim, index).unwrap().global_index())
         .collect_vec();
 
     // Number of cells. We are only interested in owned cells.
@@ -171,9 +165,21 @@ where
     }
 
     if return_transpose {
-        DistributedCsrMatrix::from_aij(range_layout, domain_layout, &cols, &rows, &data).into()
+        Operator::from(DistributedCsrMatrix::from_aij(
+            range_layout,
+            domain_layout,
+            &cols,
+            &rows,
+            &data,
+        ))
     } else {
-        DistributedCsrMatrix::from_aij(domain_layout, range_layout, &rows, &cols, &data).into()
+        Operator::from(DistributedCsrMatrix::from_aij(
+            domain_layout,
+            range_layout,
+            &rows,
+            &cols,
+            &data,
+        ))
     }
 }
 
@@ -235,17 +241,19 @@ impl<
 
         let active_cells = active_cells
             .iter()
-            .sorted_by_key(|&index| grid.entity(tdim, index).unwrap().global_index())
             .copied()
+            .sorted_by_key(|&index| grid.entity(tdim, index).unwrap().global_index())
             .collect_vec();
 
-        let owned_active_cells = active_cells
+        let owned_active_cells: Vec<usize> = active_cells
+            .iter()
             .filter(|index| {
                 matches!(
                     grid.entity(tdim, **index).unwrap().ownership(),
                     Ownership::Owned
                 )
             })
+            .copied()
             .collect_vec();
 
         let n_cells = owned_active_cells.len();
@@ -273,8 +281,8 @@ impl<
         // We also want to create a map from global cell indices to owned active cell indices.
         let mut global_cell_index_to_owned_active_cell_index = HashMap::<usize, usize>::default();
 
-        for (owned_active_index, cell_index) in owned_active_cells.iter().enumerate() {
-            let cell = grid.entity(tdim, *cell_index).unwrap();
+        for (owned_active_index, &cell_index) in owned_active_cells.iter().enumerate() {
+            let cell = grid.entity(tdim, cell_index).unwrap();
             global_cell_index_to_owned_active_cell_index
                 .insert(cell.global_index(), owned_active_index);
         }

src/greens_function_evaluators/dense_evaluator.rs

@@ -1,19 +1,21 @@
 //! A dense evaluator for Green's functions.
 
+use std::rc::Rc;
+
 use green_kernels::{traits::DistributedKernelEvaluator, types::GreenKernelEvalType};
 use mpi::traits::{Communicator, Equivalence};
+use num::One;
 use rlst::{
-    operator::interface::DistributedArrayVectorSpace, rlst_dynamic_array1, AsApply, Element,
-    IndexLayout, OperatorBase, RawAccess, RawAccessMut, RlstScalar,
+    operator::{interface::DistributedArrayVectorSpace, zero_element},
+    rlst_dynamic_array1, AsApply, Element, IndexLayout, OperatorBase, RawAccess, RawAccessMut,
+    RlstScalar,
 };
 
 /// Wrapper for a dense Green's function evaluator.
 pub struct DenseEvaluator<
     'a,
     C: Communicator,
     T: RlstScalar + Equivalence,
-    SourceLayout: IndexLayout<Comm = C>,
-    TargetLayout: IndexLayout<Comm = C>,
     K: DistributedKernelEvaluator<T = T>,
 > where
     T::Real: Equivalence,
@@ -23,17 +25,12 @@ pub struct DenseEvaluator<
     eval_mode: GreenKernelEvalType,
     use_multithreaded: bool,
     kernel: K,
-    domain_space: &'a DistributedArrayVectorSpace<'a, SourceLayout, T>,
-    range_space: &'a DistributedArrayVectorSpace<'a, TargetLayout, T>,
+    domain_space: Rc<DistributedArrayVectorSpace<'a, C, T>>,
+    range_space: Rc<DistributedArrayVectorSpace<'a, C, T>>,
 }
 
-impl<
-        C: Communicator,
-        T: RlstScalar + Equivalence,
-        SourceLayout: IndexLayout<Comm = C>,
-        TargetLayout: IndexLayout<Comm = C>,
-        K: DistributedKernelEvaluator<T = T>,
-    > std::fmt::Debug for DenseEvaluator<'_, C, T, SourceLayout, TargetLayout, K>
+impl<C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator<T = T>>
+    std::fmt::Debug for DenseEvaluator<'_, C, T, K>
 where
     T::Real: Equivalence,
 {
@@ -47,14 +44,8 @@ where
     }
 }
 
-impl<
-        'a,
-        C: Communicator,
-        T: RlstScalar + Equivalence,
-        SourceLayout: IndexLayout<Comm = C>,
-        TargetLayout: IndexLayout<Comm = C>,
-        K: DistributedKernelEvaluator<T = T>,
-    > DenseEvaluator<'a, C, T, SourceLayout, TargetLayout, K>
+impl<'a, C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator<T = T>>
+    DenseEvaluator<'a, C, T, K>
 where
     T::Real: Equivalence,
 {
@@ -65,8 +56,8 @@ where
         eval_mode: GreenKernelEvalType,
         use_multithreaded: bool,
         kernel: K,
-        domain_space: &'a DistributedArrayVectorSpace<'a, SourceLayout, T>,
-        range_space: &'a DistributedArrayVectorSpace<'a, TargetLayout, T>,
+        domain_space: Rc<DistributedArrayVectorSpace<'a, C, T>>,
+        range_space: Rc<DistributedArrayVectorSpace<'a, C, T>>,
     ) -> Self {
         // We want that both layouts have the same communicator.
         assert!(std::ptr::addr_eq(domain_space.comm(), range_space.comm()));
@@ -112,37 +103,26 @@ where
     }
 }
 
-impl<
-        'a,
-        C: Communicator,
-        T: RlstScalar + Equivalence,
-        SourceLayout: IndexLayout<Comm = C>,
-        TargetLayout: IndexLayout<Comm = C>,
-        K: DistributedKernelEvaluator<T = T>,
-    > OperatorBase for DenseEvaluator<'a, C, T, SourceLayout, TargetLayout, K>
+impl<'a, C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator<T = T>>
+    OperatorBase for DenseEvaluator<'a, C, T, K>
 where
     T::Real: Equivalence,
 {
-    type Domain = DistributedArrayVectorSpace<'a, SourceLayout, T>;
+    type Domain = DistributedArrayVectorSpace<'a, C, T>;
 
-    type Range = DistributedArrayVectorSpace<'a, TargetLayout, T>;
+    type Range = DistributedArrayVectorSpace<'a, C, T>;
 
-    fn domain(&self) -> &Self::Domain {
-        self.domain_space
+    fn domain(&self) -> Rc<Self::Domain> {
+        self.domain_space.clone()
     }
 
-    fn range(&self) -> &Self::Range {
-        self.range_space
+    fn range(&self) -> Rc<Self::Range> {
+        self.range_space.clone()
     }
 }
 
-impl<
-        C: Communicator,
-        T: RlstScalar + Equivalence,
-        SourceLayout: IndexLayout<Comm = C>,
-        TargetLayout: IndexLayout<Comm = C>,
-        K: DistributedKernelEvaluator<T = T>,
-    > AsApply for DenseEvaluator<'_, C, T, SourceLayout, TargetLayout, K>
+impl<C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator<T = T>> AsApply
+    for DenseEvaluator<'_, C, T, K>
 where
     T::Real: Equivalence,
 {
@@ -152,7 +132,7 @@ where
         x: &<Self::Domain as rlst::LinearSpace>::E,
         beta: <Self::Range as rlst::LinearSpace>::F,
         y: &mut <Self::Range as rlst::LinearSpace>::E,
-    ) -> rlst::RlstResult<()> {
+    ) {
         y.scale_inplace(beta);
         let mut charges = rlst_dynamic_array1!(
             T,
@@ -170,7 +150,14 @@ where
             self.use_multithreaded,
             self.domain_space.comm(), // domain space and range space have the same communicator
         );
+    }
 
-        Ok(())
+    fn apply(
+        &self,
+        x: &<Self::Domain as rlst::LinearSpace>::E,
+    ) -> <Self::Range as rlst::LinearSpace>::E {
+        let mut y = zero_element(self.range());
+        self.apply_extended(<T as One>::one(), x, T::zero(), &mut y);
+        y
     }
 }