From 4590040ef20eb57933dcca67c1d0aaf773afc902 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 10 Oct 2023 10:00:52 +0100 Subject: [PATCH] First batched assembler (#122) * format data from bempp-cl to take up less space * start working on batched assembly * workin on batched assembly * Working on batched assembly * clippy * Use rayon in assembly * move memory assignment, evaluation of basis functions out of function * assemble the batches before moving on to next colour (!) * start planning out singular assembly * formatting * add example for timing assembly * cty * update parallel grid * batch the singular terms too * only make raw output once * run tests with --release on CI too --- .github/workflows/run-tests.yml | 4 + .github/workflows/run-weekly-tests.yml | 5 + bem/Cargo.toml | 1 + bem/examples/assembly_timing.rs | 44 ++ bem/src/assembly.rs | 92 ++++ bem/src/assembly/batched.rs | 659 +++++++++++++++++++++++ bem/src/assembly/dense.rs | 718 +------------------------ bem/src/function_space.rs | 15 +- bem/src/green.rs | 2 +- grid/src/grid.rs | 465 ++++++---------- grid/src/parallel_grid.rs | 24 +- tools/src/arrays.rs | 3 +- traits/src/grid.rs | 34 +- 13 files changed, 1015 insertions(+), 1051 deletions(-) create mode 100644 bem/examples/assembly_timing.rs create mode 100644 bem/src/assembly/batched.rs diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml index 36ec4561..85eb5ef3 100644 --- a/.github/workflows/run-tests.yml +++ b/.github/workflows/run-tests.yml @@ -44,8 +44,12 @@ jobs: - name: Run unit tests run: cargo test --lib --features "strict" + - name: Run unit tests (release) + run: cargo test --lib --release --features "strict" - name: Run unit tests (with mpi enabled) run: cargo test --lib --features "mpi,strict" + - name: Run unit tests (release with mpi enabled) + run: cargo test --lib --release --features "mpi,strict" - name: Run tests run: cargo test --examples --release --features "mpi,strict" - name: Run examples diff --git a/.github/workflows/run-weekly-tests.yml b/.github/workflows/run-weekly-tests.yml index 302e2ca7..c5cdede4 100644 --- a/.github/workflows/run-weekly-tests.yml +++ b/.github/workflows/run-weekly-tests.yml @@ -36,10 +36,15 @@ jobs: - name: Build rust library (release with mpi) run: cargo build --release --features "strict,mpi" + - name: Run unit tests run: cargo test --lib --features "strict" + - name: Run unit tests (release) + run: cargo test --lib --release --features "strict" - name: Run unit tests (with mpi enabled) run: cargo test --lib --features "mpi,strict" + - name: Run unit tests (release with mpi enabled) + run: cargo test --lib --release --features "mpi,strict" - name: Run tests run: cargo test --examples --release --features "mpi,strict" - name: Run examples diff --git a/bem/Cargo.toml b/bem/Cargo.toml index d7f78fab..b202493c 100644 --- a/bem/Cargo.toml +++ b/bem/Cargo.toml @@ -29,3 +29,4 @@ approx = "0.5" itertools = "0.10" mpi = { version = "0.6.*", optional = true } num = "0.4" +rayon = "1.7" diff --git a/bem/examples/assembly_timing.rs b/bem/examples/assembly_timing.rs new file mode 100644 index 00000000..e84258d5 --- /dev/null +++ b/bem/examples/assembly_timing.rs @@ -0,0 +1,44 @@ +use bempp_bem::assembly::{assemble_batched, BoundaryOperator, PDEType}; +use bempp_bem::function_space::SerialFunctionSpace; +use bempp_element::element::create_element; +use bempp_grid::shapes::regular_sphere; +use bempp_tools::arrays::Array2D; +use bempp_traits::bem::DofMap; +use bempp_traits::bem::FunctionSpace; +use bempp_traits::cell::ReferenceCellType; +use bempp_traits::element::{Continuity, ElementFamily}; +use num::complex::Complex; +use std::time::Instant; + +fn main() { + for i in 0..5 { + let now = Instant::now(); + let grid = regular_sphere(i); + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 0, + Continuity::Discontinuous, + ); + + let space = SerialFunctionSpace::new(&grid, &element); + let mut matrix = Array2D::>::new(( + space.dofmap().global_size(), + space.dofmap().global_size(), + )); + + assemble_batched( + &mut matrix, + BoundaryOperator::SingleLayer, + PDEType::Helmholtz(5.0), + &space, + &space, + ); + + println!( + "{} {}", + space.dofmap().global_size(), + now.elapsed().as_millis() + ) + } +} diff --git a/bem/src/assembly.rs b/bem/src/assembly.rs index 1161c41e..0091fe2a 100644 --- a/bem/src/assembly.rs +++ b/bem/src/assembly.rs @@ -1,4 +1,6 @@ +pub mod batched; pub mod dense; +use crate::function_space::SerialFunctionSpace; use crate::green; use crate::green::Scalar; use bempp_tools::arrays::Array2D; @@ -112,6 +114,96 @@ pub fn assemble_dense<'a, T: Scalar>( }, }; } + +/// Assemble an operator into a dense matrix using batched parallelisation +pub fn assemble_batched<'a, T: Scalar + Copy + Sync>( + // TODO: ouput should be `&mut impl ArrayAccess2D` once such a trait exists + output: &mut Array2D, + operator: BoundaryOperator, + pde: PDEType, + trial_space: &SerialFunctionSpace<'a>, + test_space: &SerialFunctionSpace<'a>, +) { + match pde { + PDEType::Laplace => match operator { + BoundaryOperator::SingleLayer => { + batched::assemble( + output, + &green::LaplaceGreenKernel {}, + false, + false, + trial_space, + test_space, + ); + } + BoundaryOperator::DoubleLayer => { + batched::assemble( + output, + &green::LaplaceGreenDyKernel {}, + false, + true, + trial_space, + test_space, + ); + } + BoundaryOperator::AdjointDoubleLayer => { + batched::assemble( + output, + &green::LaplaceGreenDxKernel {}, + true, + false, + trial_space, + test_space, + ); + } + //BoundaryOperator::Hypersingular => { + // batched::laplace_hypersingular_assemble(output, trial_space, test_space); + //} + _ => { + panic!("Invalid operator"); + } + }, + PDEType::Helmholtz(k) => match operator { + BoundaryOperator::SingleLayer => { + batched::assemble( + output, + &green::HelmholtzGreenKernel { k }, + false, + false, + trial_space, + test_space, + ); + } + BoundaryOperator::DoubleLayer => { + batched::assemble( + output, + &green::HelmholtzGreenDyKernel { k }, + false, + true, + trial_space, + test_space, + ); + } + BoundaryOperator::AdjointDoubleLayer => { + batched::assemble( + output, + &green::HelmholtzGreenDxKernel { k }, + true, + false, + trial_space, + test_space, + ); + } + //BoundaryOperator::Hypersingular => { + // batched::helmholtz_hypersingular_assemble(output, trial_space, test_space, k); + //} + _ => { + panic!("Invalid operator"); + } + }, + }; +} + #[cfg(test)] mod test { use crate::assembly::dense; diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs new file mode 100644 index 00000000..6d7c913e --- /dev/null +++ b/bem/src/assembly/batched.rs @@ -0,0 +1,659 @@ +use crate::function_space::SerialFunctionSpace; +use crate::green::{ + //HelmholtzGreenHypersingularTermKernel, HelmholtzGreenKernel, LaplaceGreenKernel, + Scalar, + SingularKernel, +}; +use bempp_quadrature::duffy::quadrilateral::quadrilateral_duffy; +use bempp_quadrature::duffy::triangle::triangle_duffy; +use bempp_quadrature::simplex_rules::simplex_rule; +use bempp_quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; +use bempp_tools::arrays::{Array2D, Array4D}; +use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array4DAccess}; +use bempp_traits::bem::{DofMap, FunctionSpace}; +use bempp_traits::cell::ReferenceCellType; +use bempp_traits::element::FiniteElement; +use bempp_traits::grid::{Geometry, Grid, Topology}; +use rayon::prelude::*; + +fn get_quadrature_rule( + test_celltype: ReferenceCellType, + trial_celltype: ReferenceCellType, + pairs: &[(usize, usize)], + npoints: usize, +) -> TestTrialNumericalQuadratureDefinition { + if pairs.is_empty() { + panic!("Non-singular rule"); + } else { + // Singular rules + if test_celltype == ReferenceCellType::Triangle { + if trial_celltype != ReferenceCellType::Triangle { + unimplemented!("Mixed meshes not yet supported"); + } + triangle_duffy( + &CellToCellConnectivity { + connectivity_dimension: if pairs.len() == 1 { + 0 + } else if pairs.len() == 2 { + 1 + } else { + 2 + }, + local_indices: pairs.to_vec(), + }, + npoints, + ) + .unwrap() + } else { + if test_celltype != ReferenceCellType::Quadrilateral { + unimplemented!("Only triangles and quadrilaterals are currently supported"); + } + if trial_celltype != ReferenceCellType::Quadrilateral { + unimplemented!("Mixed meshes not yet supported"); + } + quadrilateral_duffy( + &CellToCellConnectivity { + connectivity_dimension: if pairs.len() == 1 { + 0 + } else if pairs.len() == 2 { + 1 + } else { + 2 + }, + local_indices: pairs.to_vec(), + }, + npoints, + ) + .unwrap() + } + } +} + +struct RawData2D { + pub data: *mut T, + pub shape: (usize, usize), +} + +unsafe impl Sync for RawData2D {} + +#[allow(clippy::too_many_arguments)] +fn assemble_batch_singular<'a, T: Scalar + Clone + Copy>( + output: &RawData2D, + kernel: &impl SingularKernel, + needs_trial_normal: bool, + needs_test_normal: bool, + trial_space: &SerialFunctionSpace<'a>, + test_space: &SerialFunctionSpace<'a>, + cell_pairs: &[(usize, usize)], + trial_points: &Array2D, + test_points: &Array2D, + weights: &[f64], + trial_table: &Array4D, + test_table: &Array4D, +) -> usize { + let grid = test_space.grid(); + + // Memory assignment to be moved elsewhere as passed into here mutable? + let mut test_jdet = vec![0.0; test_points.shape().0]; + let mut trial_jdet = vec![0.0; trial_points.shape().0]; + let mut test_mapped_pts = Array2D::::new((test_points.shape().0, 3)); + let mut trial_mapped_pts = Array2D::::new((trial_points.shape().0, 3)); + let mut test_normals = Array2D::::new((test_points.shape().0, 3)); + let mut trial_normals = Array2D::::new((trial_points.shape().0, 3)); + + for (test_cell, trial_cell) in cell_pairs { + let test_cell_tindex = grid.topology().index_map()[*test_cell]; + let test_cell_gindex = grid.geometry().index_map()[*test_cell]; + let trial_cell_tindex = grid.topology().index_map()[*trial_cell]; + let trial_cell_gindex = grid.geometry().index_map()[*trial_cell]; + + grid.geometry().compute_jacobian_determinants( + test_points, + test_cell_gindex, + &mut test_jdet, + ); + grid.geometry() + .compute_points(test_points, test_cell_gindex, &mut test_mapped_pts); + if needs_test_normal { + grid.geometry() + .compute_normals(test_points, test_cell_gindex, &mut test_normals); + } + + grid.geometry().compute_jacobian_determinants( + trial_points, + trial_cell_gindex, + &mut trial_jdet, + ); + grid.geometry() + .compute_points(trial_points, trial_cell_gindex, &mut trial_mapped_pts); + if needs_trial_normal { + grid.geometry() + .compute_normals(trial_points, trial_cell_gindex, &mut trial_normals); + } + + for (test_i, test_dof) in test_space + .dofmap() + .cell_dofs(test_cell_tindex) + .unwrap() + .iter() + .enumerate() + { + for (trial_i, trial_dof) in trial_space + .dofmap() + .cell_dofs(trial_cell_tindex) + .unwrap() + .iter() + .enumerate() + { + let mut sum = T::zero(); + + for (index, wt) in weights.iter().enumerate() { + sum += kernel.eval::( + unsafe { test_mapped_pts.row_unchecked(index) }, + unsafe { trial_mapped_pts.row_unchecked(index) }, + unsafe { test_normals.row_unchecked(index) }, + unsafe { trial_normals.row_unchecked(index) }, + ) * T::from_f64( + wt * unsafe { test_table.get_unchecked(0, index, test_i, 0) } + * test_jdet[index] + * unsafe { trial_table.get_unchecked(0, index, trial_i, 0) } + * trial_jdet[index], + ); + } + unsafe { + *output.data.offset( + (*test_dof + output.shape.0 * *trial_dof) + .try_into() + .unwrap(), + ) += sum; + } + } + } + } + 1 +} + +#[allow(clippy::too_many_arguments)] +fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( + output: &RawData2D, + kernel: &impl SingularKernel, + needs_trial_normal: bool, + needs_test_normal: bool, + trial_space: &SerialFunctionSpace<'a>, + trial_cells: &[usize], + test_space: &SerialFunctionSpace<'a>, + test_cells: &[usize], + trial_points: &Array2D, + trial_weights: &[f64], + test_points: &Array2D, + test_weights: &[f64], + trial_table: &Array4D, + test_table: &Array4D, +) -> usize { + let test_grid = test_space.grid(); + let test_c20 = test_grid.topology().connectivity(2, 0); + let trial_grid = trial_space.grid(); + let trial_c20 = trial_grid.topology().connectivity(2, 0); + + // Memory assignment to be moved elsewhere as passed into here mutable? + let mut test_jdet = vec![0.0; test_points.shape().0]; + let mut trial_jdet = vec![0.0; trial_points.shape().0]; + let mut test_mapped_pts = Array2D::::new((test_points.shape().0, 3)); + let mut trial_mapped_pts = Array2D::::new((trial_points.shape().0, 3)); + let mut test_normals = Array2D::::new((test_points.shape().0, 3)); + let mut trial_normals = Array2D::::new((trial_points.shape().0, 3)); + + for test_cell in test_cells { + let test_cell_tindex = test_grid.topology().index_map()[*test_cell]; + let test_cell_gindex = test_grid.geometry().index_map()[*test_cell]; + let test_vertices = test_c20.row(test_cell_tindex).unwrap(); + + test_grid.geometry().compute_jacobian_determinants( + test_points, + test_cell_gindex, + &mut test_jdet, + ); + test_grid + .geometry() + .compute_points(test_points, test_cell_gindex, &mut test_mapped_pts); + if needs_test_normal { + test_grid + .geometry() + .compute_normals(test_points, test_cell_gindex, &mut test_normals); + } + + for trial_cell in trial_cells { + let trial_cell_tindex = trial_grid.topology().index_map()[*trial_cell]; + let trial_cell_gindex = trial_grid.geometry().index_map()[*trial_cell]; + let trial_vertices = trial_c20.row(trial_cell_tindex).unwrap(); + + trial_grid.geometry().compute_jacobian_determinants( + trial_points, + trial_cell_gindex, + &mut trial_jdet, + ); + trial_grid.geometry().compute_points( + trial_points, + trial_cell_gindex, + &mut trial_mapped_pts, + ); + if needs_trial_normal { + trial_grid.geometry().compute_normals( + trial_points, + trial_cell_gindex, + &mut trial_normals, + ); + } + + for (test_i, test_dof) in test_space + .dofmap() + .cell_dofs(test_cell_tindex) + .unwrap() + .iter() + .enumerate() + { + for (trial_i, trial_dof) in trial_space + .dofmap() + .cell_dofs(trial_cell_tindex) + .unwrap() + .iter() + .enumerate() + { + let mut sum = T::zero(); + + for (test_index, test_wt) in test_weights.iter().enumerate() { + for (trial_index, trial_wt) in trial_weights.iter().enumerate() { + sum += kernel.eval::( + unsafe { test_mapped_pts.row_unchecked(test_index) }, + unsafe { trial_mapped_pts.row_unchecked(trial_index) }, + unsafe { test_normals.row_unchecked(test_index) }, + unsafe { trial_normals.row_unchecked(trial_index) }, + ) * T::from_f64( + test_wt + * trial_wt + * unsafe { test_table.get_unchecked(0, test_index, test_i, 0) } + * test_jdet[test_index] + * unsafe { + trial_table.get_unchecked(0, trial_index, trial_i, 0) + } + * trial_jdet[test_index], + ); + } + } + // TODO: should we write into a result array, then copy into output after this loop? + let mut neighbour = false; + for v in test_vertices { + if trial_vertices.contains(v) { + neighbour = true; + break; + } + } + if !neighbour { + unsafe { + *output.data.offset( + (*test_dof + output.shape.0 * *trial_dof) + .try_into() + .unwrap(), + ) += sum; + } + } + } + } + } + } + 1 +} + +pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( + output: &mut Array2D, + kernel: &impl SingularKernel, + needs_trial_normal: bool, + needs_test_normal: bool, + trial_space: &SerialFunctionSpace<'a>, + test_space: &SerialFunctionSpace<'a>, +) { + // Note: currently assumes that the two grids are the same + // TODO: implement == and != for grids, then add: + // if *trial_space.grid() != *test_space.grid() { + // unimplemented!("Assembling operators with spaces on different grids not yet supported"); + // } + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assemble can only be used for function spaces stored in serial"); + } + if output.shape().0 != test_space.dofmap().global_size() + || output.shape().1 != trial_space.dofmap().global_size() + { + panic!("Matrix has wrong shape"); + } + + // TODO: make these configurable + let blocksize = 128; + + // Size of this might not be known at compile time + // let test_dofs_per_cell = 1; + // let trial_dofs_per_cell = 1; + + // TODO: allow user to configure this + let npoints = 16; + + // TODO: pass cell types into this function + let qrule = simplex_rule(ReferenceCellType::Triangle, npoints).unwrap(); + let qpoints = Array2D::from_data(qrule.points, (qrule.npoints, 2)); + let qweights = qrule.weights; + + let mut test_table = Array4D::::new( + test_space + .element() + .tabulate_array_shape(0, qpoints.shape().0), + ); + test_space.element().tabulate(&qpoints, 0, &mut test_table); + + let mut trial_table = Array4D::::new( + trial_space + .element() + .tabulate_array_shape(0, qpoints.shape().0), + ); + trial_space + .element() + .tabulate(&qpoints, 0, &mut trial_table); + + let output_raw = RawData2D { + data: output.data.as_mut_ptr(), + shape: *output.shape(), + }; + + let test_colouring = test_space.compute_cell_colouring(); + let trial_colouring = trial_space.compute_cell_colouring(); + for test_c in &test_colouring { + for trial_c in &trial_colouring { + let mut test_cells: Vec<&[usize]> = vec![]; + let mut trial_cells: Vec<&[usize]> = vec![]; + + let mut test_start = 0; + while test_start < test_c.len() { + let test_end = if test_start + blocksize < test_c.len() { + test_start + blocksize + } else { + test_c.len() + }; + + let mut trial_start = 0; + while trial_start < trial_c.len() { + let trial_end = if trial_start + blocksize < trial_c.len() { + trial_start + blocksize + } else { + trial_c.len() + }; + test_cells.push(&test_c[test_start..test_end]); + trial_cells.push(&trial_c[trial_start..trial_end]); + trial_start = trial_end; + } + test_start = test_end + } + + let numthreads = test_cells.len(); + let r: usize = (0..numthreads) + .into_par_iter() + .map(&|t| { + assemble_batch_nonadjacent( + &output_raw, + kernel, + needs_trial_normal, + needs_test_normal, + trial_space, + trial_cells[t], + test_space, + test_cells[t], + &qpoints, + &qweights, + &qpoints, + &qweights, + &trial_table, + &test_table, + ) + }) + .sum(); + assert_eq!(r, numthreads); + } + } + + if test_space.grid() != trial_space.grid() { + // If the test and trial grids are different, there are no neighbouring triangles + return; + } + + // TODO: allow user to configure this + let npoints = 4; + + let grid = test_space.grid(); + let c20 = grid.topology().connectivity(2, 0); + + let mut possible_pairs = vec![]; + // Vertex-adjacent + for i in 0..3 { + for j in 0..3 { + possible_pairs.push(vec![(i, j)]); + } + } + // edge-adjacent + for i in 0..3 { + for j in i + 1..3 { + for k in 0..3 { + for l in 0..3 { + if k != l { + possible_pairs.push(vec![(i, k), (j, l)]); + } + } + } + } + } + // Same cell + possible_pairs.push(vec![(0, 0), (1, 1), (2, 2)]); + + let mut qweights = vec![]; + let mut trial_points = vec![]; + let mut test_points = vec![]; + let mut trial_tables = vec![]; + let mut test_tables = vec![]; + for pairs in &possible_pairs { + let qrule = get_quadrature_rule( + ReferenceCellType::Triangle, + ReferenceCellType::Triangle, + pairs, + npoints, + ); + + let points = Array2D::from_data(qrule.trial_points, (qrule.npoints, 2)); + let mut table = Array4D::::new( + trial_space + .element() + .tabulate_array_shape(0, points.shape().0), + ); + trial_space.element().tabulate(&points, 0, &mut table); + trial_points.push(points); + trial_tables.push(table); + + let points = Array2D::from_data(qrule.test_points, (qrule.npoints, 2)); + let mut table = Array4D::::new( + test_space + .element() + .tabulate_array_shape(0, points.shape().0), + ); + test_space.element().tabulate(&points, 0, &mut table); + test_points.push(points); + test_tables.push(table); + qweights.push(qrule.weights); + } + + for test_c in &test_colouring { + for trial_c in &trial_colouring { + let mut cell_pairs: Vec> = vec![vec![]; possible_pairs.len()]; + for test_cell in test_c { + let test_cell_tindex = grid.topology().index_map()[*test_cell]; + let test_vertices = c20.row(test_cell_tindex).unwrap(); + for trial_cell in trial_c { + let trial_cell_tindex = grid.topology().index_map()[*trial_cell]; + let trial_vertices = c20.row(trial_cell_tindex).unwrap(); + + let mut pairs = vec![]; + for (test_i, test_v) in test_vertices.iter().enumerate() { + for (trial_i, trial_v) in trial_vertices.iter().enumerate() { + if test_v == trial_v { + pairs.push((test_i, trial_i)); + } + } + } + if !pairs.is_empty() { + cell_pairs[possible_pairs.iter().position(|r| *r == pairs).unwrap()] + .push((*trial_cell, *test_cell)) + } + } + } + for (i, cells) in cell_pairs.iter().enumerate() { + let mut start = 0; + let mut cell_blocks = vec![]; + while start < cells.len() { + let end = if start + blocksize < cells.len() { + start + blocksize + } else { + cells.len() + }; + cell_blocks.push(&cells[start..end]); + start = end; + } + + let numthreads = cell_blocks.len(); + let r: usize = (0..numthreads) + .into_par_iter() + .map(&|t| { + assemble_batch_singular( + &output_raw, + kernel, + needs_trial_normal, + needs_test_normal, + trial_space, + test_space, + cell_blocks[t], + &trial_points[i], + &test_points[i], + &qweights[i], + &trial_tables[i], + &test_tables[i], + ) + }) + .sum(); + assert_eq!(r, numthreads); + } + } + } +} +#[cfg(test)] +mod test { + use crate::assembly::batched::*; + use crate::assembly::dense; + use crate::function_space::SerialFunctionSpace; + use crate::green; + use approx::*; + use bempp_element::element::create_element; + use bempp_grid::shapes::regular_sphere; + use bempp_traits::cell::ReferenceCellType; + use bempp_traits::element::{Continuity, ElementFamily}; + // use num::complex::Complex; + + #[cfg_attr(debug_assertions, ignore)] + #[test] + fn test_laplace_single_layer_dp0_dp0_larger() { + let grid = regular_sphere(2); + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 0, + Continuity::Discontinuous, + ); + let space = SerialFunctionSpace::new(&grid, &element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::::new((ndofs, ndofs)); + assemble( + &mut matrix, + &green::LaplaceGreenKernel {}, + false, + false, + &space, + &space, + ); + + let mut dmat = Array2D::::new((ndofs, ndofs)); + dense::assemble( + &mut dmat, + &green::LaplaceGreenKernel {}, + false, + false, + &space, + &space, + ); + + /*for i in 0..matrix.shape().0 { + for j in 0..matrix.shape().1 { + println!("{} {}", + *matrix.get(i, j).unwrap(), + *dmat.get(i, j).unwrap(), + ); + } + }*/ + for i in 0..matrix.shape().0 { + for j in 0..matrix.shape().1 { + assert_relative_eq!( + *matrix.get(i, j).unwrap(), + *dmat.get(i, j).unwrap(), + epsilon = 1e-4 + ); + } + } + } + + #[test] + fn test_laplace_single_layer_dp0_dp0_smaller() { + let grid = regular_sphere(1); + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 0, + Continuity::Discontinuous, + ); + let space = SerialFunctionSpace::new(&grid, &element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::::new((ndofs, ndofs)); + assemble( + &mut matrix, + &green::LaplaceGreenKernel {}, + false, + false, + &space, + &space, + ); + + let mut dmat = Array2D::::new((ndofs, ndofs)); + dense::assemble( + &mut dmat, + &green::LaplaceGreenKernel {}, + false, + false, + &space, + &space, + ); + + for i in 0..matrix.shape().0 { + for j in 0..matrix.shape().1 { + assert_relative_eq!( + *matrix.get(i, j).unwrap(), + *dmat.get(i, j).unwrap(), + epsilon = 1e-3 + ); + } + } + } +} diff --git a/bem/src/assembly/dense.rs b/bem/src/assembly/dense.rs index 3f9e6c36..888239ea 100644 --- a/bem/src/assembly/dense.rs +++ b/bem/src/assembly/dense.rs @@ -143,24 +143,11 @@ pub fn assemble<'a, T: Scalar>( let test_cell_gindex = grid.geometry().index_map()[test_cell]; let test_vertices = c20.row(test_cell_tindex).unwrap(); - let mut npoints_test_cell = 10 * npoints * npoints; - for p in available_rules(grid.topology().cell_type(test_cell_tindex).unwrap()) { - if p >= npoints * npoints && p < npoints_test_cell { - npoints_test_cell = p; - } - } for trial_cell in 0..grid.geometry().cell_count() { let trial_cell_tindex = grid.topology().index_map()[trial_cell]; let trial_cell_gindex = grid.geometry().index_map()[trial_cell]; let trial_vertices = c20.row(trial_cell_tindex).unwrap(); - let mut npoints_trial_cell = 10 * npoints * npoints; - for p in available_rules(grid.topology().cell_type(trial_cell_tindex).unwrap()) { - if p >= npoints * npoints && p < npoints_trial_cell { - npoints_trial_cell = p; - } - } - let mut pairs = vec![]; for (test_i, test_v) in test_vertices.iter().enumerate() { for (trial_i, trial_v) in trial_vertices.iter().enumerate() { @@ -292,24 +279,11 @@ pub fn curl_curl_assemble<'a, T: Scalar>( let test_cell_gindex = grid.geometry().index_map()[test_cell]; let test_vertices = c20.row(test_cell_tindex).unwrap(); - let mut npoints_test_cell = 10 * npoints * npoints; - for p in available_rules(grid.topology().cell_type(test_cell_tindex).unwrap()) { - if p >= npoints * npoints && p < npoints_test_cell { - npoints_test_cell = p; - } - } for trial_cell in 0..grid.geometry().cell_count() { let trial_cell_tindex = grid.topology().index_map()[trial_cell]; let trial_cell_gindex = grid.geometry().index_map()[trial_cell]; let trial_vertices = c20.row(trial_cell_tindex).unwrap(); - let mut npoints_trial_cell = 10 * npoints * npoints; - for p in available_rules(grid.topology().cell_type(trial_cell_tindex).unwrap()) { - if p >= npoints * npoints && p < npoints_trial_cell { - npoints_trial_cell = p; - } - } - let mut pairs = vec![]; for (test_i, test_v) in test_vertices.iter().enumerate() { for (trial_i, trial_v) in trial_vertices.iter().enumerate() { @@ -526,88 +500,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - 0.1854538822982487, - 0.08755414595678074, - 0.05963897421514472, - 0.08755414595678074, - 0.08755414595678074, - 0.05963897421514473, - 0.04670742127454548, - 0.05963897421514472, - ], - vec![ - 0.08755414595678074, - 0.1854538822982487, - 0.08755414595678074, - 0.05963897421514472, - 0.05963897421514472, - 0.08755414595678074, - 0.05963897421514473, - 0.04670742127454548, - ], - vec![ - 0.05963897421514472, - 0.08755414595678074, - 0.1854538822982487, - 0.08755414595678074, - 0.04670742127454548, - 0.05963897421514472, - 0.08755414595678074, - 0.05963897421514473, - ], - vec![ - 0.08755414595678074, - 0.05963897421514472, - 0.08755414595678074, - 0.1854538822982487, - 0.05963897421514473, - 0.04670742127454548, - 0.05963897421514472, - 0.08755414595678074, - ], - vec![ - 0.08755414595678074, - 0.05963897421514472, - 0.046707421274545476, - 0.05963897421514473, - 0.1854538822982487, - 0.08755414595678074, - 0.05963897421514472, - 0.08755414595678074, - ], - vec![ - 0.05963897421514473, - 0.08755414595678074, - 0.05963897421514472, - 0.046707421274545476, - 0.08755414595678074, - 0.1854538822982487, - 0.08755414595678074, - 0.05963897421514472, - ], - vec![ - 0.046707421274545476, - 0.05963897421514473, - 0.08755414595678074, - 0.05963897421514472, - 0.05963897421514472, - 0.08755414595678074, - 0.1854538822982487, - 0.08755414595678074, - ], - vec![ - 0.05963897421514472, - 0.046707421274545476, - 0.05963897421514473, - 0.08755414595678074, - 0.08755414595678074, - 0.05963897421514472, - 0.08755414595678074, - 0.1854538822982487, - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.08755414595678074, 0.05963897421514473, 0.04670742127454548, 0.05963897421514472], vec![0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.05963897421514472, 0.08755414595678074, 0.05963897421514473, 0.04670742127454548], vec![0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.04670742127454548, 0.05963897421514472, 0.08755414595678074, 0.05963897421514473], vec![0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.05963897421514473, 0.04670742127454548, 0.05963897421514472, 0.08755414595678074], vec![0.08755414595678074, 0.05963897421514472, 0.046707421274545476, 0.05963897421514473, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074], vec![0.05963897421514473, 0.08755414595678074, 0.05963897421514472, 0.046707421274545476, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472], vec![0.046707421274545476, 0.05963897421514473, 0.08755414595678074, 0.05963897421514472, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074], vec![0.05963897421514472, 0.046707421274545476, 0.05963897421514473, 0.08755414595678074, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -640,88 +534,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - -1.9658941517361406e-33, - -0.08477786720045567, - -0.048343860959178774, - -0.08477786720045567, - -0.08477786720045566, - -0.048343860959178774, - -0.033625570841778946, - -0.04834386095917877, - ], - vec![ - -0.08477786720045567, - -1.9658941517361406e-33, - -0.08477786720045567, - -0.048343860959178774, - -0.04834386095917877, - -0.08477786720045566, - -0.048343860959178774, - -0.033625570841778946, - ], - vec![ - -0.048343860959178774, - -0.08477786720045567, - -1.9658941517361406e-33, - -0.08477786720045567, - -0.033625570841778946, - -0.04834386095917877, - -0.08477786720045566, - -0.048343860959178774, - ], - vec![ - -0.08477786720045567, - -0.048343860959178774, - -0.08477786720045567, - -1.9658941517361406e-33, - -0.048343860959178774, - -0.033625570841778946, - -0.04834386095917877, - -0.08477786720045566, - ], - vec![ - -0.08477786720045566, - -0.04834386095917877, - -0.033625570841778946, - -0.04834386095917877, - 4.910045345075783e-33, - -0.08477786720045566, - -0.048343860959178774, - -0.08477786720045566, - ], - vec![ - -0.04834386095917877, - -0.08477786720045566, - -0.04834386095917877, - -0.033625570841778946, - -0.08477786720045566, - 4.910045345075783e-33, - -0.08477786720045566, - -0.048343860959178774, - ], - vec![ - -0.033625570841778946, - -0.04834386095917877, - -0.08477786720045566, - -0.04834386095917877, - -0.048343860959178774, - -0.08477786720045566, - 4.910045345075783e-33, - -0.08477786720045566, - ], - vec![ - -0.04834386095917877, - -0.033625570841778946, - -0.04834386095917877, - -0.08477786720045566, - -0.08477786720045566, - -0.048343860959178774, - -0.08477786720045566, - 4.910045345075783e-33, - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![-1.9658941517361406e-33, -0.08477786720045567, -0.048343860959178774, -0.08477786720045567, -0.08477786720045566, -0.048343860959178774, -0.033625570841778946, -0.04834386095917877], vec![-0.08477786720045567, -1.9658941517361406e-33, -0.08477786720045567, -0.048343860959178774, -0.04834386095917877, -0.08477786720045566, -0.048343860959178774, -0.033625570841778946], vec![-0.048343860959178774, -0.08477786720045567, -1.9658941517361406e-33, -0.08477786720045567, -0.033625570841778946, -0.04834386095917877, -0.08477786720045566, -0.048343860959178774], vec![-0.08477786720045567, -0.048343860959178774, -0.08477786720045567, -1.9658941517361406e-33, -0.048343860959178774, -0.033625570841778946, -0.04834386095917877, -0.08477786720045566], vec![-0.08477786720045566, -0.04834386095917877, -0.033625570841778946, -0.04834386095917877, 4.910045345075783e-33, -0.08477786720045566, -0.048343860959178774, -0.08477786720045566], vec![-0.04834386095917877, -0.08477786720045566, -0.04834386095917877, -0.033625570841778946, -0.08477786720045566, 4.910045345075783e-33, -0.08477786720045566, -0.048343860959178774], vec![-0.033625570841778946, -0.04834386095917877, -0.08477786720045566, -0.04834386095917877, -0.048343860959178774, -0.08477786720045566, 4.910045345075783e-33, -0.08477786720045566], vec![-0.04834386095917877, -0.033625570841778946, -0.04834386095917877, -0.08477786720045566, -0.08477786720045566, -0.048343860959178774, -0.08477786720045566, 4.910045345075783e-33]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -754,88 +568,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - 1.9658941517361406e-33, - -0.08478435261011981, - -0.048343860959178774, - -0.0847843526101198, - -0.08478435261011981, - -0.04834386095917877, - -0.033625570841778946, - -0.048343860959178774, - ], - vec![ - -0.0847843526101198, - 1.9658941517361406e-33, - -0.08478435261011981, - -0.048343860959178774, - -0.048343860959178774, - -0.08478435261011981, - -0.04834386095917877, - -0.033625570841778946, - ], - vec![ - -0.048343860959178774, - -0.0847843526101198, - 1.9658941517361406e-33, - -0.08478435261011981, - -0.033625570841778946, - -0.048343860959178774, - -0.08478435261011981, - -0.04834386095917877, - ], - vec![ - -0.08478435261011981, - -0.048343860959178774, - -0.0847843526101198, - 1.9658941517361406e-33, - -0.04834386095917877, - -0.033625570841778946, - -0.048343860959178774, - -0.08478435261011981, - ], - vec![ - -0.0847843526101198, - -0.04834386095917877, - -0.033625570841778946, - -0.04834386095917877, - -4.910045345075783e-33, - -0.0847843526101198, - -0.048343860959178774, - -0.08478435261011981, - ], - vec![ - -0.04834386095917877, - -0.0847843526101198, - -0.04834386095917877, - -0.033625570841778946, - -0.08478435261011981, - -4.910045345075783e-33, - -0.0847843526101198, - -0.048343860959178774, - ], - vec![ - -0.033625570841778946, - -0.04834386095917877, - -0.0847843526101198, - -0.04834386095917877, - -0.048343860959178774, - -0.08478435261011981, - -4.910045345075783e-33, - -0.0847843526101198, - ], - vec![ - -0.04834386095917877, - -0.033625570841778946, - -0.04834386095917877, - -0.0847843526101198, - -0.0847843526101198, - -0.048343860959178774, - -0.08478435261011981, - -4.910045345075783e-33, - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![1.9658941517361406e-33, -0.08478435261011981, -0.048343860959178774, -0.0847843526101198, -0.08478435261011981, -0.04834386095917877, -0.033625570841778946, -0.048343860959178774], vec![-0.0847843526101198, 1.9658941517361406e-33, -0.08478435261011981, -0.048343860959178774, -0.048343860959178774, -0.08478435261011981, -0.04834386095917877, -0.033625570841778946], vec![-0.048343860959178774, -0.0847843526101198, 1.9658941517361406e-33, -0.08478435261011981, -0.033625570841778946, -0.048343860959178774, -0.08478435261011981, -0.04834386095917877], vec![-0.08478435261011981, -0.048343860959178774, -0.0847843526101198, 1.9658941517361406e-33, -0.04834386095917877, -0.033625570841778946, -0.048343860959178774, -0.08478435261011981], vec![-0.0847843526101198, -0.04834386095917877, -0.033625570841778946, -0.04834386095917877, -4.910045345075783e-33, -0.0847843526101198, -0.048343860959178774, -0.08478435261011981], vec![-0.04834386095917877, -0.0847843526101198, -0.04834386095917877, -0.033625570841778946, -0.08478435261011981, -4.910045345075783e-33, -0.0847843526101198, -0.048343860959178774], vec![-0.033625570841778946, -0.04834386095917877, -0.0847843526101198, -0.04834386095917877, -0.048343860959178774, -0.08478435261011981, -4.910045345075783e-33, -0.0847843526101198], vec![-0.04834386095917877, -0.033625570841778946, -0.04834386095917877, -0.0847843526101198, -0.0847843526101198, -0.048343860959178774, -0.08478435261011981, -4.910045345075783e-33]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -885,56 +619,8 @@ mod test { laplace_hypersingular_assemble(&mut matrix, &space, &space); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - 0.33550642155494004, - -0.10892459915262698, - -0.05664545560057827, - -0.05664545560057828, - -0.0566454556005783, - -0.05664545560057828, - ], - vec![ - -0.10892459915262698, - 0.33550642155494004, - -0.05664545560057828, - -0.05664545560057827, - -0.05664545560057828, - -0.05664545560057829, - ], - vec![ - -0.05664545560057828, - -0.05664545560057827, - 0.33550642155494004, - -0.10892459915262698, - -0.056645455600578286, - -0.05664545560057829, - ], - vec![ - -0.05664545560057827, - -0.05664545560057828, - -0.10892459915262698, - 0.33550642155494004, - -0.05664545560057828, - -0.056645455600578286, - ], - vec![ - -0.05664545560057829, - -0.0566454556005783, - -0.05664545560057829, - -0.05664545560057829, - 0.33550642155494004, - -0.10892459915262698, - ], - vec![ - -0.05664545560057829, - -0.05664545560057831, - -0.05664545560057829, - -0.05664545560057829, - -0.10892459915262698, - 0.33550642155494004, - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![0.33550642155494004, -0.10892459915262698, -0.05664545560057827, -0.05664545560057828, -0.0566454556005783, -0.05664545560057828], vec![-0.10892459915262698, 0.33550642155494004, -0.05664545560057828, -0.05664545560057827, -0.05664545560057828, -0.05664545560057829], vec![-0.05664545560057828, -0.05664545560057827, 0.33550642155494004, -0.10892459915262698, -0.056645455600578286, -0.05664545560057829], vec![-0.05664545560057827, -0.05664545560057828, -0.10892459915262698, 0.33550642155494004, -0.05664545560057828, -0.056645455600578286], vec![-0.05664545560057829, -0.0566454556005783, -0.05664545560057829, -0.05664545560057829, 0.33550642155494004, -0.10892459915262698], vec![-0.05664545560057829, -0.05664545560057831, -0.05664545560057829, -0.05664545560057829, -0.10892459915262698, 0.33550642155494004]]; let perm = [0, 5, 2, 4, 3, 1]; @@ -973,88 +659,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - 0.08742460357596939, - -0.02332791148192136, - -0.04211947809894265, - -0.02332791148192136, - -0.023327911481921364, - -0.042119478098942634, - -0.03447046598405515, - -0.04211947809894265, - ], - vec![ - -0.023327911481921364, - 0.08742460357596939, - -0.02332791148192136, - -0.04211947809894265, - -0.04211947809894265, - -0.02332791148192136, - -0.042119478098942634, - -0.03447046598405515, - ], - vec![ - -0.04211947809894265, - -0.02332791148192136, - 0.08742460357596939, - -0.02332791148192136, - -0.03447046598405515, - -0.04211947809894265, - -0.023327911481921364, - -0.042119478098942634, - ], - vec![ - -0.02332791148192136, - -0.04211947809894265, - -0.023327911481921364, - 0.08742460357596939, - -0.042119478098942634, - -0.03447046598405515, - -0.04211947809894265, - -0.02332791148192136, - ], - vec![ - -0.023327911481921364, - -0.04211947809894265, - -0.03447046598405515, - -0.042119478098942634, - 0.08742460357596939, - -0.02332791148192136, - -0.04211947809894265, - -0.023327911481921364, - ], - vec![ - -0.042119478098942634, - -0.02332791148192136, - -0.04211947809894265, - -0.034470465984055156, - -0.02332791148192136, - 0.08742460357596939, - -0.023327911481921364, - -0.04211947809894265, - ], - vec![ - -0.03447046598405515, - -0.042119478098942634, - -0.023327911481921364, - -0.04211947809894265, - -0.04211947809894265, - -0.023327911481921364, - 0.08742460357596939, - -0.02332791148192136, - ], - vec![ - -0.04211947809894265, - -0.034470465984055156, - -0.042119478098942634, - -0.02332791148192136, - -0.023327911481921364, - -0.04211947809894265, - -0.02332791148192136, - 0.08742460357596939, - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![0.08742460357596939, -0.02332791148192136, -0.04211947809894265, -0.02332791148192136, -0.023327911481921364, -0.042119478098942634, -0.03447046598405515, -0.04211947809894265], vec![-0.023327911481921364, 0.08742460357596939, -0.02332791148192136, -0.04211947809894265, -0.04211947809894265, -0.02332791148192136, -0.042119478098942634, -0.03447046598405515], vec![-0.04211947809894265, -0.02332791148192136, 0.08742460357596939, -0.02332791148192136, -0.03447046598405515, -0.04211947809894265, -0.023327911481921364, -0.042119478098942634], vec![-0.02332791148192136, -0.04211947809894265, -0.023327911481921364, 0.08742460357596939, -0.042119478098942634, -0.03447046598405515, -0.04211947809894265, -0.02332791148192136], vec![-0.023327911481921364, -0.04211947809894265, -0.03447046598405515, -0.042119478098942634, 0.08742460357596939, -0.02332791148192136, -0.04211947809894265, -0.023327911481921364], vec![-0.042119478098942634, -0.02332791148192136, -0.04211947809894265, -0.034470465984055156, -0.02332791148192136, 0.08742460357596939, -0.023327911481921364, -0.04211947809894265], vec![-0.03447046598405515, -0.042119478098942634, -0.023327911481921364, -0.04211947809894265, -0.04211947809894265, -0.023327911481921364, 0.08742460357596939, -0.02332791148192136], vec![-0.04211947809894265, -0.034470465984055156, -0.042119478098942634, -0.02332791148192136, -0.023327911481921364, -0.04211947809894265, -0.02332791148192136, 0.08742460357596939]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -1086,88 +692,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - Complex::new(0.08742460357596939, 0.11004203436820102), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(-0.04211947809894265, 0.003720159902487029), - Complex::new(-0.02332791148192136, 0.04919102584271125), - Complex::new(-0.023327911481921364, 0.04919102584271124), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.03447046598405515, -0.02816544680626108), - Complex::new(-0.04211947809894265, 0.0037201599024870254), - ], - vec![ - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(0.08742460357596939, 0.11004203436820104), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.04211947809894265, 0.0037201599024870254), - Complex::new(-0.02332791148192136, 0.04919102584271125), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.03447046598405515, -0.028165446806261072), - ], - vec![ - Complex::new(-0.04211947809894265, 0.003720159902487029), - Complex::new(-0.02332791148192136, 0.04919102584271125), - Complex::new(0.08742460357596939, 0.11004203436820102), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(-0.03447046598405515, -0.02816544680626108), - Complex::new(-0.04211947809894265, 0.0037201599024870254), - Complex::new(-0.023327911481921364, 0.04919102584271124), - Complex::new(-0.042119478098942634, 0.003720159902487025), - ], - vec![ - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(0.08742460357596939, 0.11004203436820104), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.03447046598405515, -0.028165446806261072), - Complex::new(-0.04211947809894265, 0.0037201599024870254), - Complex::new(-0.02332791148192136, 0.04919102584271125), - ], - vec![ - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.03447046598405515, -0.02816544680626108), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(0.08742460357596939, 0.11004203436820104), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(-0.04211947809894265, 0.0037201599024870267), - Complex::new(-0.023327911481921364, 0.04919102584271125), - ], - vec![ - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.02332791148192136, 0.04919102584271125), - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.034470465984055156, -0.028165446806261075), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(0.08742460357596939, 0.11004203436820104), - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(-0.04211947809894265, 0.0037201599024870237), - ], - vec![ - Complex::new(-0.03447046598405515, -0.02816544680626108), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.04211947809894265, 0.0037201599024870267), - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(0.08742460357596939, 0.11004203436820104), - Complex::new(-0.02332791148192136, 0.04919102584271124), - ], - vec![ - Complex::new(-0.04211947809894265, 0.0037201599024870263), - Complex::new(-0.034470465984055156, -0.028165446806261075), - Complex::new(-0.042119478098942634, 0.003720159902487025), - Complex::new(-0.02332791148192136, 0.04919102584271125), - Complex::new(-0.023327911481921364, 0.04919102584271125), - Complex::new(-0.04211947809894265, 0.0037201599024870237), - Complex::new(-0.02332791148192136, 0.04919102584271124), - Complex::new(0.08742460357596939, 0.11004203436820104), - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![Complex::new(0.08742460357596939, 0.11004203436820102), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(-0.04211947809894265, 0.003720159902487029), Complex::new(-0.02332791148192136, 0.04919102584271125), Complex::new(-0.023327911481921364, 0.04919102584271124), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.03447046598405515, -0.02816544680626108), Complex::new(-0.04211947809894265, 0.0037201599024870254)], vec![Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(0.08742460357596939, 0.11004203436820104), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.04211947809894265, 0.0037201599024870254), Complex::new(-0.02332791148192136, 0.04919102584271125), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.03447046598405515, -0.028165446806261072)], vec![Complex::new(-0.04211947809894265, 0.003720159902487029), Complex::new(-0.02332791148192136, 0.04919102584271125), Complex::new(0.08742460357596939, 0.11004203436820102), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(-0.03447046598405515, -0.02816544680626108), Complex::new(-0.04211947809894265, 0.0037201599024870254), Complex::new(-0.023327911481921364, 0.04919102584271124), Complex::new(-0.042119478098942634, 0.003720159902487025)], vec![Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(0.08742460357596939, 0.11004203436820104), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.03447046598405515, -0.028165446806261072), Complex::new(-0.04211947809894265, 0.0037201599024870254), Complex::new(-0.02332791148192136, 0.04919102584271125)], vec![Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.03447046598405515, -0.02816544680626108), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(0.08742460357596939, 0.11004203436820104), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(-0.04211947809894265, 0.0037201599024870267), Complex::new(-0.023327911481921364, 0.04919102584271125)], vec![Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.02332791148192136, 0.04919102584271125), Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.034470465984055156, -0.028165446806261075), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(0.08742460357596939, 0.11004203436820104), Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(-0.04211947809894265, 0.0037201599024870237)], vec![Complex::new(-0.03447046598405515, -0.02816544680626108), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.04211947809894265, 0.0037201599024870267), Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(0.08742460357596939, 0.11004203436820104), Complex::new(-0.02332791148192136, 0.04919102584271124)], vec![Complex::new(-0.04211947809894265, 0.0037201599024870263), Complex::new(-0.034470465984055156, -0.028165446806261075), Complex::new(-0.042119478098942634, 0.003720159902487025), Complex::new(-0.02332791148192136, 0.04919102584271125), Complex::new(-0.023327911481921364, 0.04919102584271125), Complex::new(-0.04211947809894265, 0.0037201599024870237), Complex::new(-0.02332791148192136, 0.04919102584271124), Complex::new(0.08742460357596939, 0.11004203436820104)]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { assert_relative_eq!(matrix.get(i, j).unwrap().re, entry.re, epsilon = 1e-4); @@ -1200,88 +726,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - Complex::new(-1.025266688854119e-33, -7.550086433767158e-36), - Complex::new(-0.07902626473768169, -0.08184681047051735), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(-0.07902626473768169, -0.08184681047051737), - Complex::new(0.01906923918000323, -0.10276858786959302), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.019069239180003215, -0.10276858786959299), - ], - vec![ - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(-1.025266688854119e-33, 1.0291684702482414e-35), - Complex::new(-0.0790262647376817, -0.08184681047051737), - Complex::new(0.019069239180003212, -0.10276858786959299), - Complex::new(0.019069239180003212, -0.10276858786959298), - Complex::new(-0.07902626473768168, -0.08184681047051737), - Complex::new(0.01906923918000323, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - ], - vec![ - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(-1.025266688854119e-33, -7.550086433767158e-36), - Complex::new(-0.07902626473768169, -0.08184681047051735), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.019069239180003215, -0.10276858786959299), - Complex::new(-0.07902626473768169, -0.08184681047051737), - Complex::new(0.01906923918000323, -0.10276858786959302), - ], - vec![ - Complex::new(-0.0790262647376817, -0.08184681047051737), - Complex::new(0.019069239180003212, -0.10276858786959299), - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(-1.025266688854119e-33, 1.0291684702482414e-35), - Complex::new(0.01906923918000323, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(0.019069239180003212, -0.10276858786959298), - Complex::new(-0.07902626473768168, -0.08184681047051737), - ], - vec![ - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(0.019069239180003215, -0.10276858786959298), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.01906923918000323, -0.10276858786959299), - Complex::new(5.00373588753262e-33, -1.8116810507789718e-36), - Complex::new(-0.07902626473768169, -0.08184681047051735), - Complex::new(0.019069239180003212, -0.10276858786959299), - Complex::new(-0.07902626473768169, -0.08184681047051737), - ], - vec![ - Complex::new(0.019069239180003222, -0.10276858786959299), - Complex::new(-0.07902626473768173, -0.08184681047051737), - Complex::new(0.01906923918000322, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(-0.07902626473768169, -0.08184681047051735), - Complex::new(7.314851820797302e-33, -1.088140415641433e-35), - Complex::new(-0.07902626473768169, -0.08184681047051737), - Complex::new(0.01906923918000322, -0.10276858786959299), - ], - vec![ - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.01906923918000323, -0.10276858786959299), - Complex::new(-0.07902626473768172, -0.08184681047051737), - Complex::new(0.019069239180003215, -0.10276858786959298), - Complex::new(0.019069239180003212, -0.10276858786959299), - Complex::new(-0.07902626473768169, -0.08184681047051737), - Complex::new(5.00373588753262e-33, -1.8116810507789718e-36), - Complex::new(-0.07902626473768169, -0.08184681047051735), - ], - vec![ - Complex::new(0.01906923918000322, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(0.019069239180003222, -0.10276858786959299), - Complex::new(-0.07902626473768173, -0.08184681047051737), - Complex::new(-0.07902626473768169, -0.08184681047051737), - Complex::new(0.01906923918000322, -0.10276858786959299), - Complex::new(-0.07902626473768169, -0.08184681047051735), - Complex::new(7.314851820797302e-33, -1.088140415641433e-35), - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![Complex::new(-1.025266688854119e-33, -7.550086433767158e-36), Complex::new(-0.07902626473768169, -0.08184681047051735), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(-0.07902626473768169, -0.08184681047051737), Complex::new(0.01906923918000323, -0.10276858786959302), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.019069239180003215, -0.10276858786959299)], vec![Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(-1.025266688854119e-33, 1.0291684702482414e-35), Complex::new(-0.0790262647376817, -0.08184681047051737), Complex::new(0.019069239180003212, -0.10276858786959299), Complex::new(0.019069239180003212, -0.10276858786959298), Complex::new(-0.07902626473768168, -0.08184681047051737), Complex::new(0.01906923918000323, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506)], vec![Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(-1.025266688854119e-33, -7.550086433767158e-36), Complex::new(-0.07902626473768169, -0.08184681047051735), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.019069239180003215, -0.10276858786959299), Complex::new(-0.07902626473768169, -0.08184681047051737), Complex::new(0.01906923918000323, -0.10276858786959302)], vec![Complex::new(-0.0790262647376817, -0.08184681047051737), Complex::new(0.019069239180003212, -0.10276858786959299), Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(-1.025266688854119e-33, 1.0291684702482414e-35), Complex::new(0.01906923918000323, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(0.019069239180003212, -0.10276858786959298), Complex::new(-0.07902626473768168, -0.08184681047051737)], vec![Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(0.019069239180003215, -0.10276858786959298), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.01906923918000323, -0.10276858786959299), Complex::new(5.00373588753262e-33, -1.8116810507789718e-36), Complex::new(-0.07902626473768169, -0.08184681047051735), Complex::new(0.019069239180003212, -0.10276858786959299), Complex::new(-0.07902626473768169, -0.08184681047051737)], vec![Complex::new(0.019069239180003222, -0.10276858786959299), Complex::new(-0.07902626473768173, -0.08184681047051737), Complex::new(0.01906923918000322, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(-0.07902626473768169, -0.08184681047051735), Complex::new(7.314851820797302e-33, -1.088140415641433e-35), Complex::new(-0.07902626473768169, -0.08184681047051737), Complex::new(0.01906923918000322, -0.10276858786959299)], vec![Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.01906923918000323, -0.10276858786959299), Complex::new(-0.07902626473768172, -0.08184681047051737), Complex::new(0.019069239180003215, -0.10276858786959298), Complex::new(0.019069239180003212, -0.10276858786959299), Complex::new(-0.07902626473768169, -0.08184681047051737), Complex::new(5.00373588753262e-33, -1.8116810507789718e-36), Complex::new(-0.07902626473768169, -0.08184681047051735)], vec![Complex::new(0.01906923918000322, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(0.019069239180003222, -0.10276858786959299), Complex::new(-0.07902626473768173, -0.08184681047051737), Complex::new(-0.07902626473768169, -0.08184681047051737), Complex::new(0.01906923918000322, -0.10276858786959299), Complex::new(-0.07902626473768169, -0.08184681047051735), Complex::new(7.314851820797302e-33, -1.088140415641433e-35)]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -1315,88 +761,8 @@ mod test { ); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - Complex::new(1.025266688854119e-33, 7.550086433767158e-36), - Complex::new(-0.079034545070751, -0.08184700030244885), - Complex::new(0.019069239180003205, -0.10276858786959298), - Complex::new(-0.07903454507075097, -0.08184700030244886), - Complex::new(-0.07903454507075099, -0.08184700030244887), - Complex::new(0.01906923918000323, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.019069239180003212, -0.10276858786959298), - ], - vec![ - Complex::new(-0.07903454507075097, -0.08184700030244885), - Complex::new(1.025266688854119e-33, -1.0291684702482414e-35), - Complex::new(-0.079034545070751, -0.08184700030244887), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244887), - Complex::new(0.019069239180003233, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - ], - vec![ - Complex::new(0.019069239180003205, -0.10276858786959298), - Complex::new(-0.07903454507075097, -0.08184700030244886), - Complex::new(1.025266688854119e-33, 7.550086433767158e-36), - Complex::new(-0.079034545070751, -0.08184700030244885), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.019069239180003212, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244887), - Complex::new(0.01906923918000323, -0.10276858786959299), - ], - vec![ - Complex::new(-0.079034545070751, -0.08184700030244887), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07903454507075097, -0.08184700030244885), - Complex::new(1.025266688854119e-33, -1.0291684702482414e-35), - Complex::new(0.019069239180003233, -0.10276858786959299), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244887), - ], - vec![ - Complex::new(-0.07903454507075099, -0.08184700030244887), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.01906923918000323, -0.10276858786959302), - Complex::new(-5.00373588753262e-33, 1.8116810507789718e-36), - Complex::new(-0.07903454507075099, -0.08184700030244885), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244886), - ], - vec![ - Complex::new(0.019069239180003233, -0.10276858786959302), - Complex::new(-0.07903454507075099, -0.08184700030244886), - Complex::new(0.019069239180003212, -0.10276858786959298), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(-0.07903454507075099, -0.08184700030244885), - Complex::new(-7.314851820797302e-33, 1.088140415641433e-35), - Complex::new(-0.07903454507075099, -0.08184700030244886), - Complex::new(0.019069239180003215, -0.10276858786959298), - ], - vec![ - Complex::new(0.10089706509966115, -0.07681163409722505), - Complex::new(0.01906923918000323, -0.10276858786959302), - Complex::new(-0.07903454507075099, -0.08184700030244887), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(0.01906923918000321, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244886), - Complex::new(-5.00373588753262e-33, 1.8116810507789718e-36), - Complex::new(-0.07903454507075099, -0.08184700030244885), - ], - vec![ - Complex::new(0.019069239180003212, -0.10276858786959298), - Complex::new(0.10089706509966115, -0.07681163409722506), - Complex::new(0.019069239180003233, -0.10276858786959302), - Complex::new(-0.07903454507075099, -0.08184700030244886), - Complex::new(-0.07903454507075099, -0.08184700030244886), - Complex::new(0.019069239180003215, -0.10276858786959298), - Complex::new(-0.07903454507075099, -0.08184700030244885), - Complex::new(-7.314851820797302e-33, 1.088140415641433e-35), - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![Complex::new(1.025266688854119e-33, 7.550086433767158e-36), Complex::new(-0.079034545070751, -0.08184700030244885), Complex::new(0.019069239180003205, -0.10276858786959298), Complex::new(-0.07903454507075097, -0.08184700030244886), Complex::new(-0.07903454507075099, -0.08184700030244887), Complex::new(0.01906923918000323, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.019069239180003212, -0.10276858786959298)], vec![Complex::new(-0.07903454507075097, -0.08184700030244885), Complex::new(1.025266688854119e-33, -1.0291684702482414e-35), Complex::new(-0.079034545070751, -0.08184700030244887), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244887), Complex::new(0.019069239180003233, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506)], vec![Complex::new(0.019069239180003205, -0.10276858786959298), Complex::new(-0.07903454507075097, -0.08184700030244886), Complex::new(1.025266688854119e-33, 7.550086433767158e-36), Complex::new(-0.079034545070751, -0.08184700030244885), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.019069239180003212, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244887), Complex::new(0.01906923918000323, -0.10276858786959299)], vec![Complex::new(-0.079034545070751, -0.08184700030244887), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07903454507075097, -0.08184700030244885), Complex::new(1.025266688854119e-33, -1.0291684702482414e-35), Complex::new(0.019069239180003233, -0.10276858786959299), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244887)], vec![Complex::new(-0.07903454507075099, -0.08184700030244887), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.01906923918000323, -0.10276858786959302), Complex::new(-5.00373588753262e-33, 1.8116810507789718e-36), Complex::new(-0.07903454507075099, -0.08184700030244885), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244886)], vec![Complex::new(0.019069239180003233, -0.10276858786959302), Complex::new(-0.07903454507075099, -0.08184700030244886), Complex::new(0.019069239180003212, -0.10276858786959298), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(-0.07903454507075099, -0.08184700030244885), Complex::new(-7.314851820797302e-33, 1.088140415641433e-35), Complex::new(-0.07903454507075099, -0.08184700030244886), Complex::new(0.019069239180003215, -0.10276858786959298)], vec![Complex::new(0.10089706509966115, -0.07681163409722505), Complex::new(0.01906923918000323, -0.10276858786959302), Complex::new(-0.07903454507075099, -0.08184700030244887), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(0.01906923918000321, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244886), Complex::new(-5.00373588753262e-33, 1.8116810507789718e-36), Complex::new(-0.07903454507075099, -0.08184700030244885)], vec![Complex::new(0.019069239180003212, -0.10276858786959298), Complex::new(0.10089706509966115, -0.07681163409722506), Complex::new(0.019069239180003233, -0.10276858786959302), Complex::new(-0.07903454507075099, -0.08184700030244886), Complex::new(-0.07903454507075099, -0.08184700030244886), Complex::new(0.019069239180003215, -0.10276858786959298), Complex::new(-0.07903454507075099, -0.08184700030244885), Complex::new(-7.314851820797302e-33, 1.088140415641433e-35)]]; for (i, row) in from_cl.iter().enumerate() { for (j, entry) in row.iter().enumerate() { @@ -1424,56 +790,8 @@ mod test { helmholtz_hypersingular_assemble(&mut matrix, &space, &space, 3.0); // Compare to result from bempp-cl - let from_cl = vec![ - vec![ - Complex::new(-0.24054975187128322, -0.37234907871793793), - Complex::new(-0.2018803657726846, -0.3708486980714607), - Complex::new(-0.31151549914430937, -0.36517694339435425), - Complex::new(-0.31146604913280734, -0.3652407688678574), - Complex::new(-0.3114620814217625, -0.36524076431695807), - Complex::new(-0.311434147468966, -0.36530056813389983), - ], - vec![ - Complex::new(-0.2018803657726846, -0.3708486980714607), - Complex::new(-0.24054975187128322, -0.3723490787179379), - Complex::new(-0.31146604913280734, -0.3652407688678574), - Complex::new(-0.31151549914430937, -0.36517694339435425), - Complex::new(-0.3114620814217625, -0.36524076431695807), - Complex::new(-0.311434147468966, -0.36530056813389983), - ], - vec![ - Complex::new(-0.31146604913280734, -0.3652407688678574), - Complex::new(-0.31151549914430937, -0.36517694339435425), - Complex::new(-0.24054975187128322, -0.3723490787179379), - Complex::new(-0.2018803657726846, -0.3708486980714607), - Complex::new(-0.31146208142176246, -0.36524076431695807), - Complex::new(-0.31143414746896597, -0.36530056813389983), - ], - vec![ - Complex::new(-0.31151549914430937, -0.36517694339435425), - Complex::new(-0.31146604913280734, -0.3652407688678574), - Complex::new(-0.2018803657726846, -0.3708486980714607), - Complex::new(-0.24054975187128322, -0.3723490787179379), - Complex::new(-0.3114620814217625, -0.36524076431695807), - Complex::new(-0.311434147468966, -0.36530056813389983), - ], - vec![ - Complex::new(-0.31146208142176257, -0.36524076431695807), - Complex::new(-0.3114620814217625, -0.3652407643169581), - Complex::new(-0.3114620814217625, -0.3652407643169581), - Complex::new(-0.3114620814217625, -0.3652407643169581), - Complex::new(-0.24056452443903534, -0.37231826606213236), - Complex::new(-0.20188036577268464, -0.37084869807146076), - ], - vec![ - Complex::new(-0.3114335658086867, -0.36530052927274986), - Complex::new(-0.31143356580868675, -0.36530052927274986), - Complex::new(-0.3114335658086867, -0.36530052927274986), - Complex::new(-0.3114335658086867, -0.36530052927274986), - Complex::new(-0.2018803657726846, -0.37084869807146076), - Complex::new(-0.2402983805938184, -0.37203286968364935), - ], - ]; + #[rustfmt::skip] + let from_cl = vec![vec![Complex::new(-0.24054975187128322, -0.37234907871793793), Complex::new(-0.2018803657726846, -0.3708486980714607), Complex::new(-0.31151549914430937, -0.36517694339435425), Complex::new(-0.31146604913280734, -0.3652407688678574), Complex::new(-0.3114620814217625, -0.36524076431695807), Complex::new(-0.311434147468966, -0.36530056813389983)], vec![Complex::new(-0.2018803657726846, -0.3708486980714607), Complex::new(-0.24054975187128322, -0.3723490787179379), Complex::new(-0.31146604913280734, -0.3652407688678574), Complex::new(-0.31151549914430937, -0.36517694339435425), Complex::new(-0.3114620814217625, -0.36524076431695807), Complex::new(-0.311434147468966, -0.36530056813389983)], vec![Complex::new(-0.31146604913280734, -0.3652407688678574), Complex::new(-0.31151549914430937, -0.36517694339435425), Complex::new(-0.24054975187128322, -0.3723490787179379), Complex::new(-0.2018803657726846, -0.3708486980714607), Complex::new(-0.31146208142176246, -0.36524076431695807), Complex::new(-0.31143414746896597, -0.36530056813389983)], vec![Complex::new(-0.31151549914430937, -0.36517694339435425), Complex::new(-0.31146604913280734, -0.3652407688678574), Complex::new(-0.2018803657726846, -0.3708486980714607), Complex::new(-0.24054975187128322, -0.3723490787179379), Complex::new(-0.3114620814217625, -0.36524076431695807), Complex::new(-0.311434147468966, -0.36530056813389983)], vec![Complex::new(-0.31146208142176257, -0.36524076431695807), Complex::new(-0.3114620814217625, -0.3652407643169581), Complex::new(-0.3114620814217625, -0.3652407643169581), Complex::new(-0.3114620814217625, -0.3652407643169581), Complex::new(-0.24056452443903534, -0.37231826606213236), Complex::new(-0.20188036577268464, -0.37084869807146076)], vec![Complex::new(-0.3114335658086867, -0.36530052927274986), Complex::new(-0.31143356580868675, -0.36530052927274986), Complex::new(-0.3114335658086867, -0.36530052927274986), Complex::new(-0.3114335658086867, -0.36530052927274986), Complex::new(-0.2018803657726846, -0.37084869807146076), Complex::new(-0.2402983805938184, -0.37203286968364935)]]; let perm = [0, 5, 2, 4, 3, 1]; diff --git a/bem/src/function_space.rs b/bem/src/function_space.rs index 2136949a..bbe1b993 100644 --- a/bem/src/function_space.rs +++ b/bem/src/function_space.rs @@ -1,18 +1,19 @@ use crate::dofmap::SerialDofMap; +use bempp_element::element::CiarletElement; use bempp_grid::grid::SerialGrid; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::bem::FunctionSpace; use bempp_traits::element::FiniteElement; use bempp_traits::grid::{Grid, Topology}; -pub struct SerialFunctionSpace<'a, E: FiniteElement> { +pub struct SerialFunctionSpace<'a> { grid: &'a SerialGrid, - element: &'a E, + element: &'a CiarletElement, dofmap: SerialDofMap, } -impl<'a, E: FiniteElement> SerialFunctionSpace<'a, E> { - pub fn new(grid: &'a SerialGrid, element: &'a E) -> Self { +impl<'a> SerialFunctionSpace<'a> { + pub fn new(grid: &'a SerialGrid, element: &'a CiarletElement) -> Self { let dofmap = SerialDofMap::new(grid, element); Self { grid, @@ -67,10 +68,10 @@ impl<'a, E: FiniteElement> SerialFunctionSpace<'a, E> { } } -impl<'a, E: FiniteElement> FunctionSpace<'a> for SerialFunctionSpace<'a, E> { +impl<'a> FunctionSpace<'a> for SerialFunctionSpace<'a> { type DofMap = SerialDofMap; type Grid = SerialGrid; - type FiniteElement = E; + type FiniteElement = CiarletElement; fn dofmap(&self) -> &Self::DofMap { &self.dofmap @@ -78,7 +79,7 @@ impl<'a, E: FiniteElement> FunctionSpace<'a> for SerialFunctionSpace<'a, E> { fn grid(&self) -> &Self::Grid { self.grid } - fn element(&self) -> &E { + fn element(&self) -> &CiarletElement { self.element } } diff --git a/bem/src/green.rs b/bem/src/green.rs index 7bb7e189..2f3d215b 100644 --- a/bem/src/green.rs +++ b/bem/src/green.rs @@ -97,7 +97,7 @@ impl Scalar for Complex { } } -pub trait SingularKernel { +pub trait SingularKernel: Sync { fn eval(&self, x: &[f64], y: &[f64], nx: &[f64], ny: &[f64]) -> T; } diff --git a/grid/src/grid.rs b/grid/src/grid.rs index cb180573..19eca8c0 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -7,7 +7,7 @@ use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; use itertools::izip; -use std::cell::{Ref, RefCell}; +use std::ptr; /// Geometry of a serial grid pub struct SerialGeometry { @@ -463,14 +463,10 @@ impl Geometry for SerialGeometry { /// Topology of a serial grid pub struct SerialTopology { dim: usize, - connectivity: Vec>>>, + connectivity: Vec>>, index_map: Vec, starts: Vec, cell_types: Vec, - facet_adjacent_cells: RefCell>, - ridge_adjacent_cells: RefCell>, - peak_adjacent_cells: RefCell>, - nonadjacent_cells: RefCell>, } fn get_reference_cell(cell_type: ReferenceCellType) -> Box { @@ -484,6 +480,8 @@ fn get_reference_cell(cell_type: ReferenceCellType) -> Box { } } +unsafe impl Sync for SerialTopology {} + impl SerialTopology { pub fn new(cells: &AdjacencyList, cell_types: &[ReferenceCellType]) -> Self { let mut index_map = vec![]; @@ -491,19 +489,22 @@ impl SerialTopology { let mut starts = vec![]; let mut cell_types_new = vec![]; let dim = get_reference_cell(cell_types[0]).dim(); + let mut connectivity = vec![]; for i in 0..dim + 1 { connectivity.push(vec![]); for _j in 0..dim + 1 { - connectivity[i].push(RefCell::new(AdjacencyList::::new())); + connectivity[i].push(AdjacencyList::::new()); } } + + // dim0 = dim, dim1 = 0 for c in cell_types { if dim != get_reference_cell(*c).dim() { panic!("Grids with cells of mixed topological dimension not supported."); } if !cell_types_new.contains(c) { - starts.push(connectivity[dim][0].borrow().num_rows()); + starts.push(connectivity[dim][0].num_rows()); cell_types_new.push(*c); let n = get_reference_cell(*c).vertex_count(); for (i, cell) in cells.iter_rows().enumerate() { @@ -517,133 +518,185 @@ impl SerialTopology { } row.push(vertices.iter().position(|&r| r == *v).unwrap()); } - connectivity[dim][0].borrow_mut().add_row(&row); + connectivity[dim][0].add_row(&row); } } } } - let facet_adjacent_cells = RefCell::new(vec![]); - let ridge_adjacent_cells = RefCell::new(vec![]); - let peak_adjacent_cells = RefCell::new(vec![]); - let nonadjacent_cells = RefCell::new(vec![]); + // dim1 == 0 + for dim0 in 1..dim { + let mut cty = AdjacencyList::::new(); + let cells = &connectivity[dim][0]; + for (i, cell_type) in cell_types_new.iter().enumerate() { + let ref_cell = get_reference_cell(*cell_type); + let ref_entities = (0..ref_cell.entity_count(dim0)) + .map(|x| ref_cell.connectivity(dim0, x, 0).unwrap()) + .collect::>>(); - Self { - dim, - connectivity, - index_map, - starts, - cell_types: cell_types_new, - facet_adjacent_cells, - ridge_adjacent_cells, - peak_adjacent_cells, - nonadjacent_cells, + let cstart = starts[i]; + let cend = if i == starts.len() - 1 { + connectivity[2][0].num_rows() + } else { + starts[i + 1] + }; + for c in cstart..cend { + let cell = unsafe { cells.row_unchecked(c) }; + for e in &ref_entities { + let vertices = e.iter().map(|x| cell[*x]).collect::>(); + let mut found = false; + for entity in cty.iter_rows() { + if all_equal(entity, &vertices) { + found = true; + break; + } + } + if !found { + cty.add_row(&vertices); + } + } + } + } + connectivity[dim0][0] = cty; } - } - fn compute_adjacent_cells(&self) { - let tdim = self.dim(); - if tdim != 2 { - panic!("Adjacent cell computation only implemented for 2D cells."); + // dim0 == dim1 == 0 + let mut nvertices = 0; + let mut cty = AdjacencyList::::new(); + let cells = &connectivity[dim][0]; + for cell in cells.iter_rows() { + for j in cell { + if *j >= nvertices { + nvertices = *j + 1; + } + } + } + for i in 0..nvertices { + cty.add_row(&[i]); + } + connectivity[0][0] = cty; + + // dim0 == dim1 + for (dim0, c) in connectivity.iter_mut().enumerate().skip(1) { + for i in 0..c[0].num_rows() { + c[dim0].add_row(&[i]); + } } - let cells_to_vertices = self.connectivity(tdim, 0); - // Compute cells adjacent via a facet (ie edge in 2D) - for cells in self.connectivity(tdim - 1, tdim).iter_rows() { - if cells.len() == 2 { - let tindex0 = self.index_map()[cells[0]]; - let vertices0 = cells_to_vertices.row(tindex0).unwrap(); - let tindex1 = self.index_map()[cells[1]]; - let vertices1 = cells_to_vertices.row(tindex1).unwrap(); - - let mut pairs: Vec = vec![0, 0, 0, 0]; - let mut n = 0; - for (i0, v0) in vertices0.iter().enumerate() { - for (i1, v1) in vertices1.iter().enumerate() { - if v0 == v1 { - pairs[2 * n] = i0.try_into().unwrap(); - pairs[2 * n + 1] = i1.try_into().unwrap(); - n += 1; + + // dim0 == dim + for dim1 in 1..dim + 1 { + let mut cty = AdjacencyList::::new(); + let entities0 = &connectivity[dim][0]; + let entities1 = &connectivity[dim1][0]; + + let mut sub_cell_types = vec![ReferenceCellType::Point; entities0.num_rows()]; + for (i, cell_type) in cell_types_new.iter().enumerate() { + let ref_cell = get_reference_cell(*cell_type); + let etypes = ref_cell.entity_types(dim); + + let cstart = starts[i]; + let cend = if i == starts.len() - 1 { + connectivity[2][0].num_rows() + } else { + starts[i + 1] + }; + for t in sub_cell_types.iter_mut().skip(cstart).take(cend) { + *t = etypes[0]; + } + } + for (ei, entity0) in entities0.iter_rows().enumerate() { + let entity = get_reference_cell(sub_cell_types[ei]); + let mut row = vec![]; + for i in 0..entity.entity_count(dim1) { + let vertices = entity + .connectivity(dim1, i, 0) + .unwrap() + .iter() + .map(|x| entity0[*x]) + .collect::>(); + for (j, entity1) in entities1.iter_rows().enumerate() { + if all_equal(&vertices, entity1) { + row.push(j); + break; } } } - self.facet_adjacent_cells.borrow_mut().push(( - cells[0], - cells[1], - pairs[0] << 6 | pairs[1] << 4 | pairs[2] << 2 | pairs[3], - )); - self.facet_adjacent_cells.borrow_mut().push(( - cells[1], - cells[0], - if pairs[1] < pairs[3] { - pairs[1] << 6 | pairs[0] << 4 | pairs[3] << 2 | pairs[2] - } else { - pairs[3] << 6 | pairs[2] << 4 | pairs[1] << 2 | pairs[0] - }, - )); + cty.add_row(&row); } + connectivity[dim][dim1] = cty } - // Compute cells adjacent via a ridge (ie vertex in 2D) - for cells in self.connectivity(tdim - 2, tdim).iter_rows() { - for c0 in cells { - let tindex0 = self.index_map()[*c0]; - let vertices0 = cells_to_vertices.row(tindex0).unwrap(); - for c1 in cells { - if *c0 < *c1 { - let tindex1 = self.index_map()[*c1]; - let vertices1 = cells_to_vertices.row(tindex1).unwrap(); - - let mut pairs: Vec = vec![0, 0]; - let mut n = 0; - for (i0, v0) in vertices0.iter().enumerate() { - for (i1, v1) in vertices1.iter().enumerate() { - if v0 == v1 { - pairs[0] = i0.try_into().unwrap(); - pairs[1] = i1.try_into().unwrap(); - n += 1; - } - } + // dim1 < dim0 + for dim1 in 1..dim + 1 { + for dim0 in dim1 + 1..dim { + let mut cty = AdjacencyList::::new(); + let entities0 = &connectivity[dim0][0]; + let entities1 = &connectivity[dim1][0]; + let cell_to_entities0 = &connectivity[dim][dim0]; + + let mut sub_cell_types = vec![ReferenceCellType::Point; entities0.num_rows()]; + for (i, cell_type) in cell_types_new.iter().enumerate() { + let ref_cell = get_reference_cell(*cell_type); + let etypes = ref_cell.entity_types(dim0); + + let cstart = starts[i]; + let cend = if i == starts.len() - 1 { + connectivity[2][0].num_rows() + } else { + starts[i + 1] + }; + for c in cstart..cend { + for (e, t) in izip!(unsafe { cell_to_entities0.row_unchecked(c) }, &etypes) + { + sub_cell_types[*e] = *t; } - if n == 1 { - self.ridge_adjacent_cells.borrow_mut().push(( - *c0, - *c1, - pairs[0] << 2 | pairs[1], - )); - self.ridge_adjacent_cells.borrow_mut().push(( - *c1, - *c0, - pairs[1] << 2 | pairs[0], - )); + } + } + for (ei, entity0) in entities0.iter_rows().enumerate() { + let entity = get_reference_cell(sub_cell_types[ei]); + let mut row = vec![]; + for i in 0..entity.entity_count(dim1) { + let vertices = entity + .connectivity(dim1, i, 0) + .unwrap() + .iter() + .map(|x| entity0[*x]) + .collect::>(); + for (j, entity1) in entities1.iter_rows().enumerate() { + if all_equal(&vertices, entity1) { + row.push(j); + break; + } } } + cty.add_row(&row); } + connectivity[dim0][dim1] = cty; } } - // Compute non-adjacent cells - for c0 in 0..cells_to_vertices.num_rows() { - let tindex0 = self.index_map()[c0]; - let vertices0 = cells_to_vertices.row(tindex0).unwrap(); - for c1 in 0..cells_to_vertices.num_rows() { - if c0 < c1 { - let tindex1 = self.index_map()[c1]; - let vertices1 = cells_to_vertices.row(tindex1).unwrap(); - let mut n = 0; - for v0 in vertices0 { - for v1 in vertices1 { - if v0 == v1 { - n += 1; - } - } - } - if n == 0 { - self.nonadjacent_cells.borrow_mut().push((c0, c1)); - self.nonadjacent_cells.borrow_mut().push((c1, c0)); + // dim1 > dim0 + for dim1 in 1..dim + 1 { + for dim0 in 0..dim1 { + let mut data = vec![vec![]; connectivity[dim0][0].num_rows()]; + for (i, row) in connectivity[dim1][dim0].iter_rows().enumerate() { + for v in row { + data[*v].push(i); } } + for row in data { + connectivity[dim0][dim1].add_row(&row); + } } } + + Self { + dim, + connectivity, + index_map, + starts, + cell_types: cell_types_new, + } } } @@ -676,11 +729,9 @@ impl Topology<'_> for SerialTopology { fn entity_count(&self, dim: usize) -> usize { self.connectivity(dim, 0).num_rows() } - fn cell(&self, index: usize) -> Option> { + fn cell(&self, index: usize) -> Option<&[usize]> { if index < self.entity_count(self.dim) { - Some(Ref::map(self.connectivity(self.dim, 0), |x| unsafe { - x.row_unchecked(index) - })) + Some(unsafe { self.connectivity(self.dim, 0).row_unchecked(index) }) } else { None } @@ -688,7 +739,7 @@ impl Topology<'_> for SerialTopology { fn cell_type(&self, index: usize) -> Option { for (i, start) in self.starts.iter().enumerate() { let end = if i == self.starts.len() - 1 { - self.connectivity[2][0].borrow().num_rows() + self.connectivity[2][0].num_rows() } else { self.starts[i + 1] }; @@ -698,165 +749,17 @@ impl Topology<'_> for SerialTopology { } None } - fn create_connectivity(&self, dim0: usize, dim1: usize) { + + fn connectivity(&self, dim0: usize, dim1: usize) -> &Self::Connectivity { if dim0 > self.dim() || dim1 > self.dim() { panic!("Dimension of connectivity should be higher than the topological dimension"); } - - if self.connectivity[dim0][dim1].borrow().num_rows() > 0 { - return; - } - - if dim0 < dim1 { - self.create_connectivity(dim0, 0); - self.create_connectivity(dim1, dim0); - let mut data = vec![vec![]; self.connectivity[dim0][0].borrow().num_rows()]; - for (i, row) in self.connectivity[dim1][dim0] - .borrow() - .iter_rows() - .enumerate() - { - for v in row { - data[*v].push(i); - } - } - for row in data { - self.connectivity[dim0][dim1].borrow_mut().add_row(&row); - } - } else if dim0 == dim1 { - if dim0 == 0 { - let mut nvertices = 0; - let cells = &self.connectivity[self.dim()][0].borrow(); - for cell in cells.iter_rows() { - for j in cell { - if *j >= nvertices { - nvertices = *j + 1; - } - } - } - for i in 0..nvertices { - self.connectivity[0][0].borrow_mut().add_row(&[i]); - } - } else { - self.create_connectivity(dim0, 0); - for i in 0..self.connectivity[dim0][0].borrow().num_rows() { - self.connectivity[dim0][dim0].borrow_mut().add_row(&[i]); - } - } - } else if dim1 == 0 { - let cells = &self.connectivity[self.dim()][0].borrow(); - for (i, cell_type) in self.cell_types.iter().enumerate() { - let ref_cell = get_reference_cell(*cell_type); - let ref_entities = (0..ref_cell.entity_count(dim0)) - .map(|x| ref_cell.connectivity(dim0, x, 0).unwrap()) - .collect::>>(); - - let cstart = self.starts[i]; - let cend = if i == self.starts.len() - 1 { - self.connectivity[2][0].borrow().num_rows() - } else { - self.starts[i + 1] - }; - for c in cstart..cend { - let cell = unsafe { cells.row_unchecked(c) }; - for e in &ref_entities { - let vertices = e.iter().map(|x| cell[*x]).collect::>(); - let mut found = false; - for entity in self.connectivity[dim0][0].borrow().iter_rows() { - if all_equal(entity, &vertices) { - found = true; - break; - } - } - if !found { - self.connectivity[dim0][0].borrow_mut().add_row(&vertices); - } - } - } - } - } else { - self.create_connectivity(dim0, 0); - self.create_connectivity(dim1, 0); - self.create_connectivity(self.dim(), dim0); - let entities0 = &self.connectivity[dim0][0].borrow(); - let entities1 = &self.connectivity[dim1][0].borrow(); - let cell_to_entities0 = &self.connectivity[self.dim()][dim0].borrow(); - - let mut cell_types = vec![ReferenceCellType::Point; entities0.num_rows()]; - for (i, cell_type) in self.cell_types.iter().enumerate() { - let ref_cell = get_reference_cell(*cell_type); - let etypes = ref_cell.entity_types(dim0); - - let cstart = self.starts[i]; - let cend = if i == self.starts.len() - 1 { - self.connectivity[2][0].borrow().num_rows() - } else { - self.starts[i + 1] - }; - for c in cstart..cend { - for (e, t) in izip!(unsafe { cell_to_entities0.row_unchecked(c) }, &etypes) { - cell_types[*e] = *t; - } - } - } - for (ei, entity0) in entities0.iter_rows().enumerate() { - let entity = get_reference_cell(cell_types[ei]); - let mut row = vec![]; - for i in 0..entity.entity_count(dim1) { - let vertices = entity - .connectivity(dim1, i, 0) - .unwrap() - .iter() - .map(|x| entity0[*x]) - .collect::>(); - for (j, entity1) in entities1.iter_rows().enumerate() { - if all_equal(&vertices, entity1) { - row.push(j); - break; - } - } - } - self.connectivity[dim0][dim1].borrow_mut().add_row(&row); - } - } - } - - fn connectivity(&self, dim0: usize, dim1: usize) -> Ref { - self.create_connectivity(dim0, dim1); - self.connectivity[dim0][dim1].borrow() + &self.connectivity[dim0][dim1] } fn entity_ownership(&self, _dim: usize, _index: usize) -> Ownership { Ownership::Owned } - - fn facet_adjacent_cells(&self) -> Ref> { - if self.facet_adjacent_cells.borrow().len() == 0 { - self.compute_adjacent_cells() - } - self.facet_adjacent_cells.borrow() - } - - fn ridge_adjacent_cells(&self) -> Ref> { - if self.ridge_adjacent_cells.borrow().len() == 0 { - self.compute_adjacent_cells() - } - self.ridge_adjacent_cells.borrow() - } - - fn peak_adjacent_cells(&self) -> Ref> { - if self.peak_adjacent_cells.borrow().len() == 0 { - self.compute_adjacent_cells() - } - self.peak_adjacent_cells.borrow() - } - - fn nonadjacent_cells(&self) -> Ref> { - if self.nonadjacent_cells.borrow().len() == 0 { - self.compute_adjacent_cells() - } - self.nonadjacent_cells.borrow() - } } /// Serial grid @@ -871,8 +774,9 @@ impl SerialGrid { cells: AdjacencyList, cell_types: Vec, ) -> Self { + let topology = SerialTopology::new(&cells, &cell_types); Self { - topology: SerialTopology::new(&cells, &cell_types), + topology, geometry: SerialGeometry::new(coordinates, &cells, &cell_types), } } @@ -895,33 +799,18 @@ impl Grid<'_> for SerialGrid { } } +impl PartialEq for SerialGrid { + fn eq(&self, other: &Self) -> bool { + ptr::eq(self, other) + } +} +impl Eq for SerialGrid {} + #[cfg(test)] mod test { use crate::grid::*; use approx::*; - #[test] - fn test_adjacent_cells() { - let g = SerialGrid::new( - Array2D::from_data( - vec![ - 0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 1.0, 0.0, 2.0, 1.0, 2.0, - 2.0, 2.0, - ], - (9, 2), - ), - AdjacencyList::from_data( - vec![ - 0, 1, 4, 0, 4, 3, 1, 2, 5, 1, 5, 4, 3, 4, 7, 3, 7, 6, 4, 5, 8, 4, 8, 7, - ], - vec![0, 3, 6, 9, 12, 15, 18, 21, 24], - ), - vec![ReferenceCellType::Triangle; 8], - ); - assert_eq!(g.topology().facet_adjacent_cells().len(), 8 * 2); - assert_eq!(g.topology().ridge_adjacent_cells().len(), (29 - 8 * 2) * 2); - assert_eq!(g.topology().peak_adjacent_cells().len(), 0); - } #[test] fn test_connectivity() { let g = SerialGrid::new( diff --git a/grid/src/parallel_grid.rs b/grid/src/parallel_grid.rs index 9706c457..fd7f318f 100644 --- a/grid/src/parallel_grid.rs +++ b/grid/src/parallel_grid.rs @@ -6,7 +6,6 @@ use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess}; use bempp_traits::cell::ReferenceCellType; use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; use mpi::{request::WaitGuard, topology::Communicator, traits::*}; -use std::cell::Ref; /// Geometry of a parallel grid pub struct ParallelGeometry<'a, C: Communicator> { @@ -151,17 +150,14 @@ impl<'a, C: Communicator> Topology<'a> for ParallelTopology<'a, C> { fn entity_count(&self, dim: usize) -> usize { self.serial_topology.entity_count(dim) } - fn cell(&self, index: usize) -> Option> { + fn cell(&self, index: usize) -> Option<&[usize]> { self.serial_topology.cell(index) } fn cell_type(&self, index: usize) -> Option { self.serial_topology.cell_type(index) } - fn create_connectivity(&self, dim0: usize, dim1: usize) { - self.serial_topology.create_connectivity(dim0, dim1) - } - fn connectivity(&self, dim0: usize, dim1: usize) -> Ref { + fn connectivity(&self, dim0: usize, dim1: usize) -> &Self::Connectivity { self.serial_topology.connectivity(dim0, dim1) } @@ -171,22 +167,6 @@ impl<'a, C: Communicator> Topology<'a> for ParallelTopology<'a, C> { } self.ownership[dim][index] } - - fn facet_adjacent_cells(&self) -> Ref> { - self.serial_topology.facet_adjacent_cells() - } - - fn ridge_adjacent_cells(&self) -> Ref> { - self.serial_topology.ridge_adjacent_cells() - } - - fn peak_adjacent_cells(&self) -> Ref> { - self.serial_topology.peak_adjacent_cells() - } - - fn nonadjacent_cells(&self) -> Ref> { - self.serial_topology.nonadjacent_cells() - } } /// Parallel grid diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index 3e13efe0..008e12ef 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -4,9 +4,10 @@ use num::Num; use std::clone::Clone; /// A two-dimensional rectangular array +#[derive(Clone)] pub struct Array2D { /// The data in the array, in row-major order - data: Vec, + pub data: Vec, /// The shape of the array shape: (usize, usize), } diff --git a/traits/src/grid.rs b/traits/src/grid.rs index 059e5126..d3c02315 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -2,7 +2,6 @@ use crate::arrays::{AdjacencyListAccess, Array2DAccess}; use crate::cell::ReferenceCellType; -use std::cell::Ref; /// The ownership of a mesh entity #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] @@ -98,45 +97,16 @@ pub trait Topology<'a> { fn entity_count(&self, dim: usize) -> usize; /// The indices of the vertices that from cell with index `index` - fn cell(&self, index: usize) -> Option>; + fn cell(&self, index: usize) -> Option<&[usize]>; /// The indices of the vertices that from cell with index `index` fn cell_type(&self, index: usize) -> Option; - /// Create the connectivity of entities of dimension `dim0` to entities of dimension `dim1` - /// - /// If this function is called multiple times, it will do nothing after the first call - fn create_connectivity(&self, dim0: usize, dim1: usize); - - /// Create the connectivity information for all dimensions - fn create_connectivity_all(&self) { - for dim0 in 0..self.dim() { - for dim1 in 0..self.dim() { - self.create_connectivity(dim0, dim1); - } - } - } - /// Get the connectivity of entities of dimension `dim0` to entities of dimension `dim1` - fn connectivity(&self, dim0: usize, dim1: usize) -> Ref; + fn connectivity(&self, dim0: usize, dim1: usize) -> &Self::Connectivity; /// Get the ownership of a mesh entity fn entity_ownership(&self, dim: usize, index: usize) -> Ownership; - - /// Get pairs of cells that are adjacent via a facet (entity of dimension tdim - 1). - /// Entries in returned vector are (cell0, cell1, connectivity_type_id) - fn facet_adjacent_cells(&self) -> Ref>; - - /// Get pairs of cells that are adjacent via a ridge (entity of dimension tdim - 2). - /// Entries in returned vector are (cell0, cell1, connectivity_type_id) - fn ridge_adjacent_cells(&self) -> Ref>; - - /// Get pairs of cells that are adjacent via a peak (entity of dimension tdim - 3). - /// Entries in returned vector are (cell0, cell1, connectivity_type_id) - fn peak_adjacent_cells(&self) -> Ref>; - - /// Get pairs of cells that are not adjacent - fn nonadjacent_cells(&self) -> Ref>; } pub trait Grid<'a> {