diff --git a/src/boundary_assemblers.rs b/src/boundary_assemblers.rs index d90ce997..536ce9b7 100644 --- a/src/boundary_assemblers.rs +++ b/src/boundary_assemblers.rs @@ -121,27 +121,15 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: BoundaryAssembler<'o, T, Integrand, K> { /// Assemble the singular part into a CSR matrix. - pub fn assemble_singular<'a, C: Communicator, Space: FunctionSpaceTrait>( + pub fn assemble_singular<'a, C: Communicator, Space: FunctionSpaceTrait>( &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()); @@ -149,8 +137,8 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: // 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, diff --git a/src/evaluator_tools.rs b/src/evaluator_tools.rs index 4c583158..41080e8d 100644 --- a/src/evaluator_tools.rs +++ b/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, ->( - function_space: &Space, +pub fn basis_to_point_map<'a, T: RlstScalar + Equivalence, Space: FunctionSpaceTrait>( + 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 = 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 = 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::::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); } diff --git a/src/greens_function_evaluators/dense_evaluator.rs b/src/greens_function_evaluators/dense_evaluator.rs index 6260d479..392d6bc6 100644 --- a/src/greens_function_evaluators/dense_evaluator.rs +++ b/src/greens_function_evaluators/dense_evaluator.rs @@ -1,10 +1,14 @@ //! 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. @@ -12,8 +16,6 @@ pub struct DenseEvaluator< 'a, C: Communicator, T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, K: DistributedKernelEvaluator, > 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>, + range_space: Rc>, } -impl< - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - K: DistributedKernelEvaluator, - > std::fmt::Debug for DenseEvaluator<'_, C, T, SourceLayout, TargetLayout, K> +impl> + 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, - TargetLayout: IndexLayout, - K: DistributedKernelEvaluator, - > DenseEvaluator<'a, C, T, SourceLayout, TargetLayout, K> +impl<'a, C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator> + 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>, + range_space: Rc>, ) -> 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, - TargetLayout: IndexLayout, - K: DistributedKernelEvaluator, - > OperatorBase for DenseEvaluator<'a, C, T, SourceLayout, TargetLayout, K> +impl<'a, C: Communicator, T: RlstScalar + Equivalence, K: DistributedKernelEvaluator> + 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_space.clone() } - fn range(&self) -> &Self::Range { - self.range_space + fn range(&self) -> Rc { + self.range_space.clone() } } -impl< - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - K: DistributedKernelEvaluator, - > AsApply for DenseEvaluator<'_, C, T, SourceLayout, TargetLayout, K> +impl> AsApply + for DenseEvaluator<'_, C, T, K> where T::Real: Equivalence, { @@ -152,7 +132,7 @@ where x: &::E, beta: ::F, y: &mut ::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: &::E, + ) -> ::E { + let mut y = zero_element(self.range()); + self.apply_extended(::one(), x, T::zero(), &mut y); + y } } diff --git a/src/greens_function_evaluators/kifmm_evaluator.rs b/src/greens_function_evaluators/kifmm_evaluator.rs index 86720cee..6388cf59 100644 --- a/src/greens_function_evaluators/kifmm_evaluator.rs +++ b/src/greens_function_evaluators/kifmm_evaluator.rs @@ -1,6 +1,6 @@ //! Interface to kifmm library -use std::cell::RefCell; +use std::{cell::RefCell, rc::Rc}; use bempp_distributed_tools::{self, permutation::DataPermutation}; @@ -24,38 +24,29 @@ use kifmm::{ use mpi::traits::{Communicator, Equivalence}; use num::Float; use rlst::{ - operator::interface::DistributedArrayVectorSpace, rlst_dynamic_array1, AsApply, Element, - IndexLayout, MatrixSvd, OperatorBase, RawAccess, RawAccessMut, RlstScalar, + operator::{interface::DistributedArrayVectorSpace, zero_element}, + rlst_dynamic_array1, AsApply, Element, IndexLayout, MatrixSvd, OperatorBase, RawAccess, + RawAccessMut, RlstScalar, }; /// This structure instantiates an FMM evaluator. -pub struct KiFmmEvaluator< - 'a, - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, -> where +pub struct KiFmmEvaluator<'a, C: Communicator, T: RlstScalar + Equivalence> +where T::Real: Equivalence, T: DftType::ComplexType>, T: Dft + AsComplex + Epsilon + MatrixSvd + Float, KiFmmMulti, FftFieldTranslation>: SourceToTargetTranslationMetadata, { - domain_space: &'a DistributedArrayVectorSpace<'a, SourceLayout, T>, - range_space: &'a DistributedArrayVectorSpace<'a, TargetLayout, T>, - source_permutation: DataPermutation<'a, SourceLayout>, - target_permutation: DataPermutation<'a, TargetLayout>, + domain_space: Rc>, + range_space: Rc>, + source_permutation: DataPermutation<'a, C>, + target_permutation: DataPermutation<'a, C>, n_permuted_sources: usize, n_permuted_targets: usize, fmm: RefCell, FftFieldTranslation>>, } -impl< - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - > std::fmt::Debug for KiFmmEvaluator<'_, C, T, SourceLayout, TargetLayout> +impl std::fmt::Debug for KiFmmEvaluator<'_, C, T> where T::Real: Equivalence, T: DftType::ComplexType>, @@ -72,13 +63,7 @@ where } } -impl< - 'a, - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - > KiFmmEvaluator<'a, C, T, SourceLayout, TargetLayout> +impl<'a, C: Communicator, T: RlstScalar + Equivalence> KiFmmEvaluator<'a, C, T> where T::Real: Equivalence, T: DftType::ComplexType>, @@ -92,8 +77,8 @@ where local_tree_depth: usize, global_tree_depth: usize, expansion_order: usize, - domain_space: &'a DistributedArrayVectorSpace<'a, SourceLayout, T>, - range_space: &'a DistributedArrayVectorSpace<'a, TargetLayout, T>, + domain_space: Rc>, + range_space: Rc>, ) -> Self { // We want that both layouts have the same communicator. assert!(std::ptr::addr_eq(domain_space.comm(), range_space.comm())); @@ -169,38 +154,28 @@ where } } -impl< - 'a, - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - > OperatorBase for KiFmmEvaluator<'a, C, T, SourceLayout, TargetLayout> +impl<'a, C: Communicator, T: RlstScalar + Equivalence> OperatorBase + for KiFmmEvaluator<'a, C, T> where T::Real: Equivalence, T: DftType::ComplexType>, T: Dft + AsComplex + Epsilon + MatrixSvd + Float, KiFmmMulti, FftFieldTranslation>: SourceToTargetTranslationMetadata, { - 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_space.clone() } - fn range(&self) -> &Self::Range { - self.range_space + fn range(&self) -> Rc { + self.range_space.clone() } } -impl< - C: Communicator, - T: RlstScalar + Equivalence, - SourceLayout: IndexLayout, - TargetLayout: IndexLayout, - > AsApply for KiFmmEvaluator<'_, C, T, SourceLayout, TargetLayout> +impl + Equivalence> AsApply for KiFmmEvaluator<'_, C, T> where T::Real: Equivalence, T: DftType::ComplexType>, @@ -219,7 +194,7 @@ where x: &::E, beta: ::F, y: &mut ::E, - ) -> rlst::RlstResult<()> { + ) { let mut x_permuted = rlst_dynamic_array1![T, [self.n_permuted_sources]]; let mut y_permuted = rlst_dynamic_array1![T, [self.n_permuted_targets]]; @@ -262,7 +237,14 @@ where y.scale_inplace(beta); // Now add the result. y.view_mut().local_mut().sum_into(y_result.r()); + } - Ok(()) + fn apply( + &self, + x: &::E, + ) -> ::E { + let mut y = zero_element(self.range()); + self.apply_extended(T::one(), x, T::zero(), &mut y); + y } }