WIP: Refactor

src/evaluator_tools.rs

@@ -2,7 +2,7 @@
 
 use std::{collections::HashMap, marker::PhantomData};
 
-use bempp_distributed_tools::GhostCommunicator;
+use bempp_distributed_tools::{GhostCommunicator, IndexLayoutFromLocalCounts};
 use green_kernels::{traits::Kernel, types::GreenKernelEvalType};
 use itertools::{izip, Itertools};
 use mpi::traits::{Communicator, Equivalence};
@@ -14,8 +14,13 @@ use ndgrid::{
 
 use rayon::prelude::*;
 use rlst::{
-    operator::interface::{
-        distributed_sparse_operator::DistributedCsrMatrixOperator, DistributedArrayVectorSpace,
+    operator::{
+        interface::{
+            distributed_array_vector_space,
+            distributed_sparse_operator::DistributedCsrMatrixOperatorImpl,
+            DistributedArrayVectorSpace,
+        },
+        Operator,
     },
     rlst_array_from_slice2, rlst_dynamic_array1, rlst_dynamic_array2, rlst_dynamic_array3,
     rlst_dynamic_array4, Array, AsApply, DefaultIterator, DistributedCsrMatrix, DistributedVector,
@@ -35,11 +40,9 @@ pub fn basis_to_point_map<
     Space: FunctionSpaceTrait<T = T>,
 >(
     function_space: &Space,
-    domain_space: &'a DistributedArrayVectorSpace<'a, DomainLayout, T>,
-    range_space: &'a DistributedArrayVectorSpace<'a, RangeLayout, T>,
     quadrature_points: &[T::Real],
     quadrature_weights: &[T::Real],
-) -> DistributedCsrMatrixOperator<'a, DomainLayout, RangeLayout, T, C>
+) -> Operator<DistributedCsrMatrixOperatorImpl<'a, DomainLayout, RangeLayout, T, C>>
 where
     T::Real: Equivalence,
 {
@@ -69,14 +72,10 @@ where
     assert_eq!(quadrature_points.len() % tdim, 0);
     assert_eq!(quadrature_points.len() / tdim, n_quadrature_points);
 
-    // Check that domain space and function space are compatible.
-
-    let n_domain_dofs = domain_space.index_layout().number_of_global_indices();
-    let n_range_dofs = range_space.index_layout().number_of_global_indices();
-
-    assert_eq!(function_space.local_size(), n_domain_dofs,);
+    // Assign number of domain and range dofs.
 
-    assert_eq!(n_cells * n_quadrature_points, n_range_dofs);
+    let n_domain_dofs = function_space.local_size();
+    let n_range_dofs = n_cells * n_quadrature_points;
 
     // All the dimensions are OK. Let's get to work. We need to iterate through the elements,
     // get the attached global dofs and the corresponding jacobian map.

src/function.rs

@@ -2,6 +2,7 @@
 
 //mod function_space;
 
+use bempp_distributed_tools::IndexLayoutFromLocalCounts;
 use mpi::request::WaitGuard;
 use mpi::traits::{Communicator, Destination, Source};
 use ndelement::ciarlet::CiarletElement;
@@ -13,6 +14,7 @@ use ndgrid::{traits::Grid, types::Ownership};
 use rlst::{MatrixInverse, RlstScalar};
 use std::collections::HashMap;
 use std::marker::PhantomData;
+use std::ops::Index;
 
 type DofList = Vec<Vec<usize>>;
 type OwnerData = Vec<(usize, usize, usize, usize)>;
@@ -58,6 +60,12 @@ pub trait LocalFunctionSpaceTrait {
 
     /// Get ownership of a local DOF
     fn ownership(&self, local_dof_index: usize) -> Ownership;
+
+    /// Get the local indices of the support cells associated with this space.
+    ///
+    /// The vector of support cells is sorted in ascending order and may contain
+    /// ghost cells who are not owned by the current process.
+    fn support_cells(&self) -> &[usize];
 }
 
 /// A function space
@@ -79,6 +87,9 @@ pub trait FunctionSpaceTrait: LocalFunctionSpaceTrait {
 
     /// Get the local function space
     fn local_space(&self) -> &Self::LocalFunctionSpace;
+
+    /// Return the index layout associated with the function space.
+    fn index_layout(&self) -> &IndexLayoutFromLocalCounts<'_, Self::C>;
 }
 
 /// Definition of a local function space.
@@ -95,6 +106,7 @@ pub struct LocalFunctionSpace<
     global_size: usize,
     global_dof_numbers: Vec<usize>,
     ownership: Vec<Ownership>,
+    support_cells: Vec<usize>,
 }
 
 impl<
@@ -124,6 +136,7 @@ impl<
             global_size,
             global_dof_numbers,
             ownership,
+            support_cells: (0..local_size).collect(),
         }
     }
 }
@@ -245,6 +258,10 @@ impl<
     fn ownership(&self, local_dof_index: usize) -> Ownership {
         self.ownership[local_dof_index]
     }
+
+    fn support_cells(&self) -> &[usize] {
+        self.support_cells.as_slice()
+    }
 }
 
 /// Implementation of a general function space.
@@ -256,6 +273,7 @@ pub struct FunctionSpace<
 > {
     grid: &'a GridImpl,
     local_space: LocalFunctionSpace<'a, T, GridImpl::LocalGrid>,
+    index_layout: IndexLayoutFromLocalCounts<'a, C>,
     _marker: PhantomData<(C, T)>,
 }
 
@@ -401,6 +419,7 @@ where
                 global_dof_numbers,
                 ownership,
             ),
+            index_layout: IndexLayoutFromLocalCounts::new(dofmap_size, comm),
             _marker: PhantomData,
         }
 
@@ -471,6 +490,10 @@ impl<
         // Syntactical workaround as rust-analyzer mixed up this ownership with entity ownership.
         LocalFunctionSpaceTrait::ownership(&self.local_space, local_dof_index)
     }
+
+    fn support_cells(&self) -> &[usize] {
+        self.local_space.support_cells()
+    }
 }
 
 impl<
@@ -493,6 +516,10 @@ impl<
     fn local_space(&self) -> &Self::LocalFunctionSpace {
         &self.local_space
     }
+
+    fn index_layout(&self) -> &IndexLayoutFromLocalCounts<'_, Self::C> {
+        &self.index_layout
+    }
 }
 
 /// Assign DOFs to entities.