diff --git a/Cargo.toml b/Cargo.toml index 333f9235..83f9b2f7 100644 --- 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 index 847874cf..71a96540 100644 --- 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 index 9e3295d1..4cd4f3b3 100644 --- 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 index d41c3ffb..b445fc7f 100644 --- 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 index 230d8fe8..8229e383 100644 --- 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::::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 index d5e7c732..45100c11 100644 --- a/src/boundary_assemblers.rs +++ b/src/boundary_assemblers.rs @@ -129,7 +129,10 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, 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, 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, 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, diff --git a/src/evaluator_tools.rs b/src/evaluator_tools.rs index b667b0a9..50d43cdc 100644 --- a/src/evaluator_tools.rs +++ b/src/evaluator_tools.rs @@ -1,10 +1,6 @@ //! Various helper functions to support evaluators. -use std::{ - collections::{HashMap, HashSet}, - marker::PhantomData, - rc::Rc, -}; +use std::{collections::HashSet, marker::PhantomData, rc::Rc}; use bempp_distributed_tools::{index_embedding::IndexEmbedding, Global2LocalDataMapper}; use green_kernels::{traits::Kernel, types::GreenKernelEvalType}; @@ -25,14 +21,13 @@ use rlst::{ Operator, }, rlst_array_from_slice2, rlst_dynamic_array1, rlst_dynamic_array2, rlst_dynamic_array3, - rlst_dynamic_array4, Array, AsApply, DefaultIterator, DistributedCsrMatrix, DistributedVector, - DynamicArray, Element, IndexLayout, OperatorBase, RawAccess, RawAccessMut, RlstScalar, Shape, - UnsafeRandomAccessByValue, UnsafeRandomAccessMut, + rlst_dynamic_array4, AsApply, DefaultIterator, DistributedCsrMatrix, DistributedVector, + Element, IndexLayout, OperatorBase, RawAccess, RawAccessMut, RlstScalar, }; use rayon::prelude::*; -use crate::function::FunctionSpaceTrait; +use crate::function::{FunctionSpaceTrait, LocalFunctionSpaceTrait}; /// Create a linear operator from the map of a basis to points. The points are sorted by global /// index of the corresponding cells in the support. The output space is the space obtained through @@ -52,33 +47,31 @@ where T::Real: Equivalence, { // Get the grid. - let grid = function_space.parallel_grid(); + let grid = function_space.grid(); // Get the rank - let rank = grid.comm().rank(); - // Topological dimension of the grid. - let tdim = grid.topology_dim(); + let tdim = grid.local_grid().topology_dim(); // The method is currently restricted to single element type grids. // So let's make sure that this is the case. - assert_eq!(grid.entity_types(tdim).len(), 1); + assert_eq!(grid.local_grid().entity_types(tdim).len(), 1); - let reference_cell = grid.entity_types(tdim)[0]; + let reference_cell = grid.local_grid().entity_types(tdim)[0]; let owned_support_cells: Vec = function_space + .local_space() .support_cells() .iter() .filter(|index| { matches!( - grid.entity(tdim, **index).unwrap().ownership(), + grid.local_grid().entity(tdim, **index).unwrap().ownership(), Ownership::Owned ) }) .copied() - .sorted_by_key(|&index| grid.entity(tdim, index).unwrap().global_index()) .collect_vec(); // Number of cells. We are only interested in owned cells. @@ -97,13 +90,8 @@ where // Create an index embedding from support cells to all local cells. - let n_owned_cells = grid - .entity_iter(tdim) - .filter(|e| matches!(e.ownership(), Ownership::Owned)) - .count(); - let index_embedding = - IndexEmbedding::new(grid.cell_index_layout(), &owned_support_cells, grid.comm()); + IndexEmbedding::new(grid.cell_layout(), &owned_support_cells, grid.comm()); // Create the index layouts. @@ -125,24 +113,33 @@ where let mut table = rlst_dynamic_array4!( T, function_space + .local_space() .element(reference_cell) .tabulate_array_shape(0, n_quadrature_points) ); function_space + .local_space() .element(reference_cell) .tabulate(&quadrature_points, 0, &mut table); // We have tabulated the basis functions on the reference element. Now need // the map to physical elements. - let geometry_evaluator = grid.geometry_map(reference_cell, quadrature_points.data()); + let geometry_evaluator = grid + .local_grid() + .geometry_map(reference_cell, quadrature_points.data()); // The following arrays hold the jacobians, their determinants and the normals. - let mut jacobians = - rlst_dynamic_array3![T::Real, [grid.geometry_dim(), tdim, n_quadrature_points]]; + let mut jacobians = rlst_dynamic_array3![ + T::Real, + [grid.local_grid().geometry_dim(), tdim, n_quadrature_points] + ]; let mut jdets = rlst_dynamic_array1![T::Real, [n_quadrature_points]]; - let mut normals = rlst_dynamic_array2![T::Real, [grid.geometry_dim(), n_quadrature_points]]; + let mut normals = rlst_dynamic_array2![ + T::Real, + [grid.local_grid().geometry_dim(), n_quadrature_points] + ]; // Now iterate through the cells of the grid, get the attached dofs and evaluate the geometry map. @@ -152,7 +149,7 @@ where let mut data = Vec::::default(); for (embedded_index, &cell_index) in owned_support_cells.iter().enumerate() { - let cell = grid.entity(tdim, cell_index).unwrap(); + let cell = grid.local_grid().entity(tdim, cell_index).unwrap(); // We need to get the global embeddeded index of the owned support cell. // First we get the global index of the owned support cell. @@ -170,10 +167,15 @@ where ); // Get the global dofs of the cell. let global_dofs = function_space + .local_space() .cell_dofs(cell_index) .unwrap() .iter() - .map(|local_dof_index| function_space.global_dof_index(*local_dof_index)) + .map(|local_dof_index| { + function_space + .local_space() + .global_dof_index(*local_dof_index) + }) .collect::>(); for (qindex, (jdet, qweight)) in izip!(jdets.iter(), quadrature_weights.iter()).enumerate() @@ -210,7 +212,7 @@ pub struct NeighbourEvaluator< 'a, T: RlstScalar + Equivalence, K: Kernel, - GridImpl: ParallelGrid, + GridImpl: ParallelGrid, C: Communicator, > { eval_points: Vec, @@ -230,7 +232,7 @@ impl< 'a, T: RlstScalar + Equivalence, K: Kernel, - GridImpl: ParallelGrid, + GridImpl: ParallelGrid, C: Communicator, > NeighbourEvaluator<'a, T, K, GridImpl, C> { @@ -242,16 +244,14 @@ impl< support_cells: &[usize], grid: &'a GridImpl, ) -> Self { - let tdim = grid.topology_dim(); + let tdim = grid.local_grid().topology_dim(); let comm = grid.comm(); // The method is currently restricted to single element type grids. // So let's make sure that this is the case. - assert_eq!(grid.entity_types(tdim).len(), 1); - - let reference_cell = grid.entity_types(tdim)[0]; + assert_eq!(grid.local_grid().entity_types(tdim).len(), 1); // Get the number of points assert_eq!(eval_points.len() % tdim, 0); @@ -269,17 +269,11 @@ impl< let support_cells = support_cells .iter() .copied() - .sorted_by_key(|&index| grid.entity(tdim, index).unwrap().global_index()) - .collect_vec(); - - let support_cell_owners = support_cells - .iter() - .map(|index| { - if let Ownership::Ghost(owner, _) = grid.entity(tdim, *index).unwrap().ownership() { - owner - } else { - comm.rank() as usize - } + .sorted_by_key(|&index| { + grid.local_grid() + .entity(tdim, index) + .unwrap() + .global_index() }) .collect_vec(); @@ -287,7 +281,10 @@ impl< .iter() .filter(|cell_index| { matches!( - grid.entity(tdim, **cell_index).unwrap().ownership(), + grid.local_grid() + .entity(tdim, **cell_index) + .unwrap() + .ownership(), Ownership::Owned ) }) @@ -298,7 +295,7 @@ impl< // easier we embed the owned support cells into the set of all owned cells. let index_embedding = - IndexEmbedding::new(grid.cell_index_layout(), &owned_support_cells, grid.comm()); + IndexEmbedding::new(grid.cell_layout(), &owned_support_cells, grid.comm()); let n_owned_support_cells = owned_support_cells.len(); @@ -312,15 +309,18 @@ impl< // We now setup the data mapper with respect to all owned cells on each process. let global_to_local_mapper = Global2LocalDataMapper::new( - grid.cell_index_layout(), + grid.cell_layout(), &support_cells .iter() - .map(|index| grid.entity(tdim, *index).unwrap().global_index()) + .map(|index| { + grid.local_grid() + .entity(tdim, *index) + .unwrap() + .global_index() + }) .collect_vec(), ); - // Finally, we need a map that tells us for each support cell index what the corresponding position in the `support_cells` vector is. - Self { eval_points: eval_points.to_vec(), n_points, @@ -342,11 +342,8 @@ impl< // and then map this vector globally across all processes to the active cells // on each process. - let mut all_cells_vec = vec![ - T::zero(); - self.grid.cell_index_layout().number_of_local_indices() - * self.n_points - ]; + let mut all_cells_vec = + vec![T::zero(); self.grid.cell_layout().number_of_local_indices() * self.n_points]; self.index_embedding .embed_data(x.local().data(), &mut all_cells_vec, self.n_points); @@ -361,7 +358,7 @@ impl< impl< T: RlstScalar + Equivalence, K: Kernel, - GridImpl: ParallelGrid, + GridImpl: ParallelGrid, C: Communicator, > std::fmt::Debug for NeighbourEvaluator<'_, T, K, GridImpl, C> { @@ -379,7 +376,7 @@ impl< 'a, T: RlstScalar + Equivalence, K: Kernel, - GridImpl: ParallelGrid, + GridImpl: ParallelGrid, C: Communicator, > OperatorBase for NeighbourEvaluator<'a, T, K, GridImpl, C> { @@ -399,7 +396,7 @@ impl< impl< T: RlstScalar + Equivalence, K: Kernel, - GridImpl: ParallelGrid, + GridImpl: ParallelGrid, C: Communicator, > AsApply for NeighbourEvaluator<'_, T, K, GridImpl, C> where @@ -418,8 +415,8 @@ where ) { // We need to iterate through the elements. - let tdim = self.grid.topology_dim(); - let gdim = self.grid.geometry_dim(); + let tdim = self.grid.local_grid().topology_dim(); + let gdim = self.grid.local_grid().geometry_dim(); // We need the chunk size of the targets. This is the chunk size of the domain multiplied // with the range component count of the kernel. @@ -429,7 +426,7 @@ where let n_points = self.n_points; // We need the reference cell - let reference_cell = self.grid.entity_types(tdim)[0]; + let reference_cell = self.grid.local_grid().entity_types(tdim)[0]; // We now need to communicate the data. The `charges` vector contains all charges for the active cells on this process. diff --git a/src/function.rs b/src/function.rs index b108ec74..c412c5a4 100644 --- a/src/function.rs +++ b/src/function.rs @@ -2,7 +2,8 @@ //mod function_space; -use bempp_distributed_tools::Global2LocalDataMapper; +use bempp_distributed_tools::{DataPermutation, Global2LocalDataMapper}; +use itertools::{izip, Itertools}; use mpi::request::WaitGuard; use mpi::traits::{Communicator, Destination, Equivalence, Source}; use ndelement::ciarlet::CiarletElement; @@ -12,12 +13,13 @@ use ndgrid::traits::Grid; use ndgrid::traits::ParallelGrid; use ndgrid::traits::{Entity, Topology}; use ndgrid::types::Ownership; +use num::Zero; +use rlst::operator::interface::DistributedArrayVectorSpace; use rlst::{ - rlst_dynamic_array2, rlst_dynamic_array4, DistributedVector, DynamicArray, IndexLayout, - MatrixInverse, RawAccess, RawAccessMut, RlstScalar, + rlst_array_from_slice2, rlst_array_from_slice_mut2, rlst_dynamic_array4, AsApply, IndexLayout, + MatrixInverse, OperatorBase, RawAccess, RawAccessMut, RlstScalar, }; use std::collections::HashMap; -use std::marker::PhantomData; use std::rc::Rc; type DofList = Vec>; @@ -50,6 +52,37 @@ pub trait LocalFunctionSpaceTrait { /// Get the local DOF numbers associated with a cell fn cell_dofs(&self, cell: usize) -> Option<&[usize]>; + /// Get an array of all cell dofs for each cell. + /// + /// Returns a flat array of all cell dofs for all lcoal cells (owned and ghost), ordered by + /// global index of the cell. + fn all_cell_dofs(&self) -> Vec { + self.support_cells() + .iter() + .flat_map(|&entity_index| self.cell_dofs(entity_index).unwrap()) + .copied() + .collect_vec() + } + + /// Get an array of all cell dofs for each owned cell. + /// + /// Returns a flat array of all cell dofs for all owned lcoal cells, ordered by + /// global index of the cell. + fn owned_cell_dofs(&self) -> Vec { + let tdim = self.grid().topology_dim(); + self.support_cells() + .iter() + .filter(|&&entity_index| { + matches!( + self.grid().entity(tdim, entity_index).unwrap().ownership(), + Ownership::Owned + ) + }) + .flat_map(|&entity_index| self.cell_dofs(entity_index).unwrap()) + .copied() + .collect_vec() + } + /// Get the local DOF numbers associated with a cell /// /// # Safety @@ -67,30 +100,32 @@ pub trait LocalFunctionSpaceTrait { /// 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 + /// The vector of support cells is sorted in ascending order by global cell index and may contain /// ghost cells who are not owned by the current process. fn support_cells(&self) -> &[usize]; } /// A function space -pub trait FunctionSpaceTrait: LocalFunctionSpaceTrait { +pub trait FunctionSpaceTrait { /// Communicator type C: Communicator; + /// Data type + type T: RlstScalar + Equivalence; + + /// Definition of a local grid + type LocalGrid: Grid::Real, EntityDescriptor = ReferenceCellType>; + /// The grid type - type ParallelGrid: ParallelGrid< - Self::C, - LocalGrid = Self::LocalGrid, + type Grid: ParallelGrid< T = ::Real, - EntityDescriptor = ReferenceCellType, - >; - /// Local Function Space - type LocalFunctionSpace: LocalFunctionSpaceTrait< - T = Self::T, + C = Self::C, LocalGrid = Self::LocalGrid, - FiniteElement = Self::FiniteElement, >; + /// Local Function Space + type LocalFunctionSpace: LocalFunctionSpaceTrait; + /// Get the communicator fn comm(&self) -> &Self::C; @@ -98,81 +133,14 @@ pub trait FunctionSpaceTrait: LocalFunctionSpaceTrait { fn local_space(&self) -> &Self::LocalFunctionSpace; /// Return the associated parallel grid - fn parallel_grid(&self) -> &Self::ParallelGrid; + fn grid(&self) -> &Self::Grid; /// Return the index layout associated with the function space. fn index_layout(&self) -> Rc>; - /// Evaluate a function space with a given coefficient vector on all support cells. - /// - /// The return value is a Hash map that maps from the cell id to a two-dimensional array - /// with shape `(value_size, npoints)`, where `value_size` is the number of range components - /// of the basis functions and `npoints` is the number of evaluation points. - /// The coefficient vector must have the same length as the `local_size` of the function space. - fn evaluate( - &self, - coefficients: &DistributedVector<'_, Self::C, Self::T>, - eval_points: &[::Real], - ) -> HashMap> - where - Self::T: Equivalence, - { - // Currently only supports single element grids. - - let tdim = self.grid().topology_dim(); - assert_eq!(self.grid().entity_types(tdim).len(), 1); - - // Get the number of eval points - assert_eq!(eval_points.len() % tdim, 0); - let npoints = eval_points.len() / tdim; - let eval_points = { - let mut tmp = rlst_dynamic_array2!(::Real, [tdim, npoints]); - tmp.data_mut().copy_from_slice(eval_points); - tmp - }; - - let reference_cell = self.grid().entity_types(tdim)[0]; - - // Get the finite element of the space - - let element = self.element(reference_cell); - - // And now tabulate the basis functions at the eval points - - let dims = element.tabulate_array_shape(0, npoints); - let mut basis_values = rlst_dynamic_array4!(Self::T, dims); - element.tabulate(&eval_points, 0, &mut basis_values); - - // We have tabulated the elements. We are going to store the results in a Hashmap. - // The keys are the insertion ids of the owned cells on the process. The values - // are two-dimensional arrays of dimension (value_size, npoints) where value_size - // is the dimension of the range of the basis functions. - - let mut result = HashMap::>::new(); - - for cell_index in self.support_cells() { - let cell = self.grid().entity(tdim, *cell_index).unwrap(); - let cell_dofs = self.cell_dofs(*cell_index).unwrap(); - // Let's instantiate the result array for the cell. - let mut cell_result = rlst_dynamic_array2!(Self::T, [element.value_size(), npoints]); - // Now we loop over the basis functions and points and add the contributions to the result. - - for dof in cell_dofs.iter() { - for point_index in 0..dims[1] { - for value_index in 0..dims[3] { - for basis_index in 0..dims[2] { - cell_result[[value_index, point_index]] += coefficients.local().data() - [*dof] - * basis_values[[0, point_index, basis_index, value_index]]; - } - } - } - } - // Store the result in the HashMap - result.insert(cell.id().unwrap(), cell_result); - } - - result + /// Number of global indices + fn global_dof_count(&self) -> usize { + *self.index_layout().counts().last().unwrap() } } @@ -216,6 +184,11 @@ impl< .iter() .map(|&i| grid.entity_count(i)) .sum(); + + let tdim = grid.topology_dim(); + let support_cells = (0..cell_count) + .sorted_by_key(|entity_index| grid.entity(tdim, *entity_index).unwrap().global_index()) + .collect_vec(); Self { grid, elements, @@ -226,7 +199,7 @@ impl< global_dof_numbers, ownership, // At the moment all spaces are global. - support_cells: (0..cell_count).collect(), + support_cells, } } } @@ -355,26 +328,20 @@ impl< } /// Implementation of a general function space. -pub struct FunctionSpace< - 'a, - C: Communicator, - T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid, -> { +pub struct FunctionSpace<'a, T: RlstScalar + MatrixInverse + Equivalence, GridImpl: ParallelGrid> +where + GridImpl::LocalGrid: Grid, +{ grid: &'a GridImpl, local_space: LocalFunctionSpace<'a, T, GridImpl::LocalGrid>, - index_layout: Rc>, - _marker: PhantomData<(C, T)>, + index_layout: Rc>, } -impl< - 'a, - C: Communicator, - T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid, - > FunctionSpace<'a, C, T, GridImpl> +impl<'a, T: RlstScalar + MatrixInverse + Equivalence, GridImpl: ParallelGrid> + FunctionSpace<'a, T, GridImpl> where - GridImpl::LocalGrid: Sync, + GridImpl::LocalGrid: Grid + Sync, + T::Real: MatrixInverse, { /// Create new function space pub fn new( @@ -394,7 +361,10 @@ where assign_dofs(rank as usize, grid.local_grid(), e_family); let mut elements = HashMap::new(); - for cell in grid.entity_types(grid.topology_dim()) { + for cell in grid + .local_grid() + .entity_types(grid.local_grid().topology_dim()) + { elements.insert(*cell, e_family.element(*cell)); } @@ -512,7 +482,6 @@ where ownership, ), index_layout: Rc::new(IndexLayout::from_local_counts(number_of_owned_dofs, comm)), - _marker: PhantomData, } // Self { @@ -529,79 +498,20 @@ where } } -impl< - C: Communicator, - T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid, - > LocalFunctionSpaceTrait for FunctionSpace<'_, C, T, GridImpl> -{ - type T = T; - - type LocalGrid = GridImpl::LocalGrid; - - type FiniteElement = CiarletElement; - - fn grid(&self) -> &Self::LocalGrid { - self.local_space.grid() - } - - fn element(&self, cell_type: ReferenceCellType) -> &Self::FiniteElement { - self.local_space.element(cell_type) - } - - fn get_local_dof_numbers(&self, entity_dim: usize, entity_number: usize) -> &[usize] { - self.local_space - .get_local_dof_numbers(entity_dim, entity_number) - } - - fn local_size(&self) -> usize { - self.local_space.local_size() - } - - fn global_size(&self) -> usize { - self.local_space.global_size() - } - - fn cell_dofs(&self, cell: usize) -> Option<&[usize]> { - self.local_space.cell_dofs(cell) - } - - unsafe fn cell_dofs_unchecked(&self, cell: usize) -> &[usize] { - self.local_space.cell_dofs_unchecked(cell) - } - - fn cell_colouring(&self) -> HashMap>> { - self.local_space.cell_colouring() - } - - fn global_dof_index(&self, local_dof_index: usize) -> usize { - self.local_space.global_dof_index(local_dof_index) - } - - fn ownership(&self, local_dof_index: usize) -> Ownership { - // 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< - 'a, - C: Communicator, - T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid, - > FunctionSpaceTrait for FunctionSpace<'a, C, T, GridImpl> +impl<'a, T: RlstScalar + MatrixInverse + Equivalence, GridImpl: ParallelGrid> + FunctionSpaceTrait for FunctionSpace<'a, T, GridImpl> +where + GridImpl::LocalGrid: Grid, { - fn comm(&self) -> &C { + fn comm(&self) -> &GridImpl::C { self.grid.comm() } - type C = C; + type T = T; + type C = GridImpl::C; - type ParallelGrid = GridImpl; + type Grid = GridImpl; + type LocalGrid = GridImpl::LocalGrid; type LocalFunctionSpace = LocalFunctionSpace<'a, T, GridImpl::LocalGrid>; @@ -613,7 +523,7 @@ impl< self.index_layout.clone() } - fn parallel_grid(&self) -> &Self::ParallelGrid { + fn grid(&self) -> &Self::Grid { self.grid } } @@ -701,24 +611,177 @@ pub fn assign_dofs< (cell_dofs, entity_dofs, size, owner_data) } +/// Evaluator to evaluate functions on the grid pub struct SpaceEvaluator<'a, Space: FunctionSpaceTrait> { space: &'a Space, global_to_local_mapper: Global2LocalDataMapper<'a, Space::C>, eval_points: Vec<::Real>, - support_cell_to_position: HashMap, + n_eval_points: usize, + id_mapper: Option>, + domain_space: Rc>, + range_space: Rc>, +} + +impl<'a, Space: FunctionSpaceTrait> SpaceEvaluator<'a, Space> { + /// Create a new space evaluator + pub fn new( + space: &'a Space, + eval_points: &[::Real], + map_to_ids: bool, + ) -> Self { + let global_to_local_mapper = Global2LocalDataMapper::new( + space.index_layout(), + &space.local_space().owned_cell_dofs(), + ); + + let tdim = space.grid().local_grid().topology_dim(); + assert_eq!(space.grid().local_grid().entity_types(tdim).len(), 1); + + let reference_cell = space.grid().local_grid().entity_types(tdim)[0]; + + let range_component_count = space.local_space().element(reference_cell).value_size(); + + let n_eval_points = { + let tdim = space.grid().local_grid().topology_dim(); + assert_eq!(eval_points.len() % tdim, 0); + eval_points.len() / tdim + }; + + let id_mapper = if map_to_ids { + let tdim = space.grid().local_grid().topology_dim(); + let owned_ids = space + .grid() + .local_grid() + .entity_iter(tdim) + .filter(|entity| matches!(entity.ownership(), Ownership::Owned)) + .map(|entity| entity.id().unwrap()) + .collect_vec(); + + Some(DataPermutation::new(space.grid().cell_layout(), &owned_ids)) + } else { + None + }; + + let range_layout = Rc::new(IndexLayout::from_local_counts( + space.grid().local_grid().cell_count() * n_eval_points * range_component_count, + space.comm(), + )); + + let domain_space = DistributedArrayVectorSpace::from_index_layout(space.index_layout()); + let range_space = DistributedArrayVectorSpace::from_index_layout(range_layout); + + Self { + space, + global_to_local_mapper, + eval_points: eval_points.to_vec(), + n_eval_points, + id_mapper, + domain_space, + range_space, + } + } +} + +impl std::fmt::Debug for SpaceEvaluator<'_, Space> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!( + f, + "SpaceEvaluatorMap on rank from {} dofs to {} cells with {} evaluation points", + self.space.global_dof_count(), + self.space.grid().global_cell_count(), + self.n_eval_points, + ) + } +} + +impl<'a, Space: FunctionSpaceTrait> OperatorBase for SpaceEvaluator<'a, Space> { + type Domain = DistributedArrayVectorSpace<'a, Space::C, Space::T>; + + type Range = DistributedArrayVectorSpace<'a, Space::C, Space::T>; + + fn domain(&self) -> Rc { + self.domain_space.clone() + } + + fn range(&self) -> Rc { + self.range_space.clone() + } } -// impl<'a, Space: FunctionSpaceTrait> SpaceEvaluator<'a, Space> { -// pub fn new(space: &'a Space, eval_points: &[::Real]) -> Self { -// // We need to iterate through and get the ownership and the associated global dof of each dof index. - -// let global_to_local_mapper = -// Global2LocalDataMapper::new(space.index_layout(), support_cells); -// let eval_points = vec![]; -// Self { -// space, -// global_to_local_mapper, -// eval_points, -// } -// } -// } +impl AsApply for SpaceEvaluator<'_, Space> { + fn apply_extended< + ContainerIn: rlst::ElementContainer::E>, + ContainerOut: rlst::ElementContainerMut::E>, + >( + &self, + alpha: ::F, + x: rlst::Element, + beta: ::F, + mut y: rlst::Element, + ) { + y *= beta; + + let tdim = self.space.local_space().grid().topology_dim(); + let gdim = self.space.local_space().grid().geometry_dim(); + + let reference_cell = self.space.local_space().grid().entity_types(tdim)[0]; + + // Get the finite element of the space + + let element = self.space.local_space().element(reference_cell); + + // We now need to map the coefficient vector x to the cellwise coefficients. + // For that we use the data mapper + + let cell_coeffs = self + .global_to_local_mapper + .map_data(x.view().local().data(), 1); + + // And now tabulate the basis functions at the eval points + + let dims = element.tabulate_array_shape(0, self.n_eval_points); + let mut basis_values = rlst_dynamic_array4!(::T, dims); + let eval_points_array = + rlst_array_from_slice2!(&self.eval_points, [gdim, self.n_eval_points]); + element.tabulate(&eval_points_array, 0, &mut basis_values); + + // We now go through each cell and evaluate the local coefficients in the cell + + for (res, basis_coeffs) in izip!( + y.view_mut() + .local_mut() + .data_mut() + .chunks_mut(element.value_size() * self.n_eval_points), + cell_coeffs.chunks(element.dim()) + ) { + // Now we loop over the basis functions and points and add the contributions to the result. + + let mut res = + rlst_array_from_slice_mut2!(res, [element.value_size(), self.n_eval_points]); + + for (basis_index, &basis_coeff) in basis_coeffs.iter().enumerate() { + for point_index in 0..dims[1] { + for value_index in 0..dims[3] { + res[[value_index, point_index]] += alpha + * basis_coeff + * basis_values[[0, point_index, basis_index, value_index]]; + } + } + } + } + + // If we have an id mapper than use that to map back to the ordering given by the ids and not by + // the global dofs. + + if let Some(id_mapper) = self.id_mapper.as_ref() { + let mut id_ordered_data = + vec![::zero(); y.view().index_layout().number_of_local_indices()]; + id_mapper.backward_permute(y.view().local().data(), &mut id_ordered_data); + + y.view_mut() + .local_mut() + .data_mut() + .copy_from_slice(&id_ordered_data); + } + } +} diff --git a/src/greens_function_evaluators/kifmm_evaluator.rs b/src/greens_function_evaluators/kifmm_evaluator.rs index 738bbed3..0f392bd9 100644 --- a/src/greens_function_evaluators/kifmm_evaluator.rs +++ b/src/greens_function_evaluators/kifmm_evaluator.rs @@ -103,19 +103,31 @@ where )); let simple_comm = domain_space.index_layout().comm().duplicate(); + + let builder = MultiNodeBuilder::new(false) + .tree( + &simple_comm, + sources, + targets, + local_tree_depth as u64, + global_tree_depth as u64, + true, + SortKind::Samplesort { n_samples: 100 }, + ) + .unwrap(); + + let n_sources = builder + .tree + .as_ref() + .unwrap() + .source_tree + .global_indices + .len(); + let cell = RefCell::new( - MultiNodeBuilder::new(false) - .tree( - &simple_comm, - sources, - targets, - local_tree_depth as u64, - global_tree_depth as u64, - true, - SortKind::Samplesort { n_samples: 100 }, - ) - .unwrap() + builder .parameters( + &vec![T::zero(); n_sources], expansion_order, Laplace3dKernel::::default(), FftFieldTranslation::::new(None), diff --git a/src/shapes.rs b/src/shapes.rs index cb837652..39dad5c6 100644 --- a/src/shapes.rs +++ b/src/shapes.rs @@ -7,7 +7,7 @@ use ndelement::{ciarlet::CiarletElement, types::ReferenceCellType}; use ndgrid::{ traits::{Builder, ParallelBuilder}, types::RealScalar, - ParallelGrid, SingleElementGrid, SingleElementGridBuilder, + ParallelGridImpl, SingleElementGrid, SingleElementGridBuilder, }; use num::Float; @@ -20,7 +20,7 @@ pub fn regular_sphere( refinement_level: u32, degree: usize, comm: &C, -) -> ParallelGrid>> { +) -> ParallelGridImpl>> { if comm.rank() == 0 { let mut b = SingleElementGridBuilder::new_with_capacity( 3, @@ -120,7 +120,7 @@ pub fn regular_sphere( pub fn screen_triangles( ncells: usize, comm: &C, -) -> ParallelGrid>> { +) -> ParallelGridImpl>> { if ncells == 0 { panic!("Cannot create a grid with 0 cells"); } @@ -178,7 +178,7 @@ pub fn screen_triangles( pub fn screen_quadrilaterals( ncells: usize, comm: &C, -) -> ParallelGrid>> { +) -> ParallelGridImpl>> { if ncells == 0 { panic!("Cannot create a grid with 0 cells"); } diff --git a/tests/test_space.rs b/tests/test_space.rs index d736de63..9e090165 100644 --- a/tests/test_space.rs +++ b/tests/test_space.rs @@ -1,4 +1,4 @@ -use bempp::function::{assign_dofs, FunctionSpace, LocalFunctionSpaceTrait}; +use bempp::function::{assign_dofs, FunctionSpace, FunctionSpaceTrait, LocalFunctionSpaceTrait}; use bempp::shapes::{regular_sphere, screen_triangles}; use mpi::traits::Communicator; use ndelement::ciarlet::{LagrangeElementFamily, RaviartThomasElementFamily}; @@ -14,15 +14,12 @@ static MPI_UNIVERSE: LazyLock = std::sync::LazyLock::new(|| { .0 }); -fn run_test< - C: Communicator, - GridImpl: ParallelGrid, ->( +fn run_test>( grid: &GridImpl, degree: usize, continuity: Continuity, ) where - GridImpl::LocalGrid: Sync, + GridImpl::LocalGrid: Grid + Sync, { let family = LagrangeElementFamily::::new(degree, continuity); let (cell_dofs, entity_dofs, size, owner_data) = assign_dofs(0, grid.local_grid(), &family); @@ -47,15 +44,12 @@ fn run_test< } } -fn run_test_rt< - C: Communicator, - GridImpl: ParallelGrid, ->( +fn run_test_rt>( grid: &GridImpl, degree: usize, continuity: Continuity, ) where - GridImpl::LocalGrid: Sync, + GridImpl::LocalGrid: Grid + Sync, { let family = RaviartThomasElementFamily::::new(degree, continuity); let (cell_dofs, entity_dofs, size, owner_data) = assign_dofs(0, grid.local_grid(), &family); @@ -153,10 +147,13 @@ fn test_dofmap_lagrange0() { //let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = FunctionSpace::new(&grid, &element); - assert_eq!(space.local_size(), space.global_size()); assert_eq!( - space.local_size(), - grid.entity_count(ReferenceCellType::Triangle) + space.local_space().local_size(), + space.local_space().global_size() + ); + assert_eq!( + space.local_space().local_size(), + grid.local_grid().entity_count(ReferenceCellType::Triangle) ); } @@ -167,10 +164,13 @@ fn test_dofmap_lagrange1() { let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = FunctionSpace::new(&grid, &element); - assert_eq!(space.local_size(), space.global_size()); assert_eq!( - space.local_size(), - grid.entity_count(ReferenceCellType::Point) + space.local_space().local_size(), + space.local_space().global_size() + ); + assert_eq!( + space.local_space().local_size(), + grid.local_grid().entity_count(ReferenceCellType::Point) ); } @@ -181,11 +181,14 @@ fn test_dofmap_lagrange2() { let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(2, Continuity::Standard); let space = FunctionSpace::new(&grid, &element); - assert_eq!(space.local_size(), space.global_size()); assert_eq!( - space.local_size(), - grid.entity_count(ReferenceCellType::Point) - + grid.entity_count(ReferenceCellType::Interval) + space.local_space().local_size(), + space.local_space().global_size() + ); + assert_eq!( + space.local_space().local_size(), + grid.local_grid().entity_count(ReferenceCellType::Point) + + grid.local_grid().entity_count(ReferenceCellType::Interval) ); } @@ -196,13 +199,16 @@ fn test_colouring_p1() { let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = FunctionSpace::new(&grid, &element); - let colouring = &space.cell_colouring()[&ReferenceCellType::Triangle]; - let cells = grid.entity_iter(2).collect::>(); + let colouring = &space.local_space().cell_colouring()[&ReferenceCellType::Triangle]; + let cells = grid.local_grid().entity_iter(2).collect::>(); let mut n = 0; for i in colouring { n += i.len() } - assert_eq!(n, grid.entity_count(ReferenceCellType::Triangle)); + assert_eq!( + n, + grid.local_grid().entity_count(ReferenceCellType::Triangle) + ); for (i, ci) in colouring.iter().enumerate() { for (j, cj) in colouring.iter().enumerate() { if i != j { @@ -236,12 +242,15 @@ fn test_colouring_dp0() { let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = FunctionSpace::new(&grid, &element); - let colouring = &space.cell_colouring()[&ReferenceCellType::Triangle]; + let colouring = &space.local_space().cell_colouring()[&ReferenceCellType::Triangle]; let mut n = 0; for i in colouring { n += i.len() } - assert_eq!(n, grid.entity_count(ReferenceCellType::Triangle)); + assert_eq!( + n, + grid.local_grid().entity_count(ReferenceCellType::Triangle) + ); for (i, ci) in colouring.iter().enumerate() { for (j, cj) in colouring.iter().enumerate() { if i != j { @@ -263,12 +272,15 @@ fn test_colouring_rt1() { let grid = regular_sphere::(2, 1, &comm); let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = FunctionSpace::new(&grid, &element); - let colouring = &space.cell_colouring()[&ReferenceCellType::Triangle]; + let colouring = &space.local_space().cell_colouring()[&ReferenceCellType::Triangle]; let mut n = 0; for i in colouring { n += i.len() } - assert_eq!(n, grid.entity_count(ReferenceCellType::Triangle)); + assert_eq!( + n, + grid.local_grid().entity_count(ReferenceCellType::Triangle) + ); for (i, ci) in colouring.iter().enumerate() { for (j, cj) in colouring.iter().enumerate() { if i != j { @@ -285,12 +297,14 @@ fn test_colouring_rt1() { for cell1 in ci { if cell0 != cell1 { for e0 in grid + .local_grid() .entity(2, *cell0) .unwrap() .topology() .sub_entity_iter(1) { for e1 in grid + .local_grid() .entity(2, *cell1) .unwrap() .topology()