diff --git a/bem/Cargo.toml b/bem/Cargo.toml index b202493c..f9fc79f6 100644 --- a/bem/Cargo.toml +++ b/bem/Cargo.toml @@ -24,9 +24,22 @@ bempp-tools = { path = "../tools"} bempp-traits = { path = "../traits"} bempp-element = { path = "../element"} bempp-grid = { path = "../grid"} +bempp-kernel = { path = "../kernel"} bempp-quadrature = { path = "../quadrature"} approx = "0.5" itertools = "0.10" mpi = { version = "0.6.*", optional = true } num = "0.4" rayon = "1.7" +rlst = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-blis-src = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-algorithms = { git = "https://github.com/linalg-rs/rlst.git" } + +[dev-dependencies] +criterion = { version = "0.3", features = ["html_reports"]} + +[[bench]] +name = "assembly_benchmark" +harness = false + diff --git a/bem/benches/assembly_benchmark.rs b/bem/benches/assembly_benchmark.rs new file mode 100644 index 00000000..38bf1986 --- /dev/null +++ b/bem/benches/assembly_benchmark.rs @@ -0,0 +1,121 @@ +use bempp_bem::assembly::{assemble_batched, batched, BoundaryOperator, PDEType}; +use bempp_bem::function_space::SerialFunctionSpace; +use bempp_element::element::create_element; +use bempp_grid::shapes::regular_sphere; +use bempp_kernel::laplace_3d; +use bempp_tools::arrays::zero_matrix; +use bempp_traits::bem::DofMap; +use bempp_traits::bem::FunctionSpace; +use bempp_traits::cell::ReferenceCellType; +use bempp_traits::element::{Continuity, ElementFamily}; +use criterion::{criterion_group, criterion_main, Criterion}; + +pub fn full_assembly_benchmark(c: &mut Criterion) { + let mut group = c.benchmark_group("assembly"); + group.sample_size(20); + + for i in 3..5 { + 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 = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size())); + + group.bench_function( + &format!( + "Full assembly of {}x{} matrix", + space.dofmap().global_size(), + space.dofmap().global_size() + ), + |b| { + b.iter(|| { + assemble_batched( + &mut matrix, + BoundaryOperator::SingleLayer, + PDEType::Laplace, + &space, + &space, + ) + }) + }, + ); + } + group.finish(); +} + +pub fn assembly_parts_benchmark(c: &mut Criterion) { + let mut group = c.benchmark_group("assembly"); + group.sample_size(20); + + for i in 3..5 { + 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 = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size())); + + let colouring = space.compute_cell_colouring(); + + /* + group.bench_function( + &format!( + "Assembly of singular terms of {}x{} matrix", + space.dofmap().global_size(), + space.dofmap().global_size() + ), + |b| { + b.iter(|| { + batched::assemble_singular( + &mut matrix, + &laplace_3d::Laplace3dKernel::new(), + false, + false, + &space, + &space, + &colouring, + &colouring, + 128, + ) + }) + }, + ); + */ + group.bench_function( + &format!( + "Assembly of non-singular terms of {}x{} matrix", + space.dofmap().global_size(), + space.dofmap().global_size() + ), + |b| { + b.iter(|| { + batched::assemble_nonsingular::<16, 16>( + &mut matrix, + &laplace_3d::Laplace3dKernel::new(), + false, + false, + &space, + &space, + &colouring, + &colouring, + 128, + ) + }) + }, + ); + } + group.finish(); +} + +// criterion_group!(benches, full_assembly_benchmark, assembly_parts_benchmark); +criterion_group!(benches, assembly_parts_benchmark); +criterion_main!(benches); diff --git a/bem/examples/assemble.rs b/bem/examples/assemble.rs index 255f2bae..f6c27ba2 100644 --- a/bem/examples/assemble.rs +++ b/bem/examples/assemble.rs @@ -1,13 +1,12 @@ -use bempp_bem::assembly::{assemble_dense, BoundaryOperator, PDEType}; +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_tools::arrays::zero_matrix; 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; fn main() { println!("Creating grid"); @@ -34,16 +33,14 @@ fn main() { space1.dofmap().global_size(), space0.dofmap().global_size() ); - let mut matrix = Array2D::>::new(( - space1.dofmap().global_size(), - space0.dofmap().global_size(), - )); + let mut matrix = + zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); println!("Assembling dense matrix (complex)"); - assemble_dense( + assemble_batched( &mut matrix, BoundaryOperator::SingleLayer, - PDEType::Helmholtz(5.0), + PDEType::Laplace, &space0, &space1, ); diff --git a/bem/examples/assembly_timing.rs b/bem/examples/assembly_timing.rs deleted file mode 100644 index e84258d5..00000000 --- a/bem/examples/assembly_timing.rs +++ /dev/null @@ -1,44 +0,0 @@ -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 0091fe2a..9ca2f193 100644 --- a/bem/src/assembly.rs +++ b/bem/src/assembly.rs @@ -1,10 +1,7 @@ pub mod batched; -pub mod dense; use crate::function_space::SerialFunctionSpace; -use crate::green; -use crate::green::Scalar; -use bempp_tools::arrays::Array2D; -use bempp_traits::bem::FunctionSpace; +use bempp_kernel::laplace_3d; +use bempp_tools::arrays::Mat; #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] #[repr(u8)] @@ -24,101 +21,10 @@ pub enum PDEType { Helmholtz(f64), } -// TODO: template over float type - -/// Assemble an operator into a dense matrix -pub fn assemble_dense<'a, T: Scalar>( - // TODO: ouput should be `&mut impl ArrayAccess2D` once such a trait exists - output: &mut Array2D, - operator: BoundaryOperator, - pde: PDEType, - trial_space: &impl FunctionSpace<'a>, - test_space: &impl FunctionSpace<'a>, -) { - match pde { - PDEType::Laplace => match operator { - BoundaryOperator::SingleLayer => { - dense::assemble( - output, - &green::LaplaceGreenKernel {}, - false, - false, - trial_space, - test_space, - ); - } - BoundaryOperator::DoubleLayer => { - dense::assemble( - output, - &green::LaplaceGreenDyKernel {}, - false, - true, - trial_space, - test_space, - ); - } - BoundaryOperator::AdjointDoubleLayer => { - dense::assemble( - output, - &green::LaplaceGreenDxKernel {}, - true, - false, - trial_space, - test_space, - ); - } - BoundaryOperator::Hypersingular => { - dense::laplace_hypersingular_assemble(output, trial_space, test_space); - } - _ => { - panic!("Invalid operator"); - } - }, - PDEType::Helmholtz(k) => match operator { - BoundaryOperator::SingleLayer => { - dense::assemble( - output, - &green::HelmholtzGreenKernel { k }, - false, - false, - trial_space, - test_space, - ); - } - BoundaryOperator::DoubleLayer => { - dense::assemble( - output, - &green::HelmholtzGreenDyKernel { k }, - false, - true, - trial_space, - test_space, - ); - } - BoundaryOperator::AdjointDoubleLayer => { - dense::assemble( - output, - &green::HelmholtzGreenDxKernel { k }, - true, - false, - trial_space, - test_space, - ); - } - BoundaryOperator::Hypersingular => { - dense::helmholtz_hypersingular_assemble(output, trial_space, test_space, k); - } - _ => { - panic!("Invalid operator"); - } - }, - }; -} - /// Assemble an operator into a dense matrix using batched parallelisation -pub fn assemble_batched<'a, T: Scalar + Copy + Sync>( +pub fn assemble_batched<'a>( // TODO: ouput should be `&mut impl ArrayAccess2D` once such a trait exists - output: &mut Array2D, + output: &mut Mat, operator: BoundaryOperator, pde: PDEType, trial_space: &SerialFunctionSpace<'a>, @@ -129,96 +35,39 @@ pub fn assemble_batched<'a, T: Scalar + Copy + Sync>( BoundaryOperator::SingleLayer => { batched::assemble( output, - &green::LaplaceGreenKernel {}, + &laplace_3d::Laplace3dKernel::new(), 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"); } }, + _ => { + panic!("Invalid PDE"); + } }; } #[cfg(test)] mod test { - use crate::assembly::dense; + use crate::assembly::batched; use crate::assembly::*; use crate::function_space::SerialFunctionSpace; - use crate::green::{HelmholtzGreenKernel, LaplaceGreenKernel}; use approx::*; use bempp_element::element::create_element; use bempp_grid::shapes::regular_sphere; - use bempp_tools::arrays::Array2D; - use bempp_traits::arrays::Array2DAccess; + use bempp_kernel::laplace_3d::Laplace3dKernel; + use bempp_tools::arrays::zero_matrix; use bempp_traits::bem::DofMap; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily}; - use num::complex::Complex; + // use num::complex::Complex; + use bempp_traits::bem::FunctionSpace; + use rlst_dense::RandomAccessByRef; #[test] fn test_laplace_single_layer() { @@ -239,10 +88,10 @@ mod test { let space1 = SerialFunctionSpace::new(&grid, &element1); let mut matrix = - Array2D::::new((space1.dofmap().global_size(), space0.dofmap().global_size())); - dense::assemble( + zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); + batched::assemble( &mut matrix, - &LaplaceGreenKernel {}, + &Laplace3dKernel::new(), false, false, &space0, @@ -250,9 +99,9 @@ mod test { ); let mut matrix2 = - Array2D::::new((space1.dofmap().global_size(), space0.dofmap().global_size())); + zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); - assemble_dense( + assemble_batched( &mut matrix2, BoundaryOperator::SingleLayer, PDEType::Laplace, @@ -270,64 +119,4 @@ mod test { } } } - - #[test] - fn test_helmholtz_single_layer() { - let grid = regular_sphere(1); - let element0 = create_element( - ElementFamily::Lagrange, - ReferenceCellType::Triangle, - 0, - Continuity::Discontinuous, - ); - let element1 = create_element( - ElementFamily::Lagrange, - ReferenceCellType::Triangle, - 1, - Continuity::Continuous, - ); - let space0 = SerialFunctionSpace::new(&grid, &element0); - let space1 = SerialFunctionSpace::new(&grid, &element1); - - let mut matrix = Array2D::>::new(( - space1.dofmap().global_size(), - space0.dofmap().global_size(), - )); - dense::assemble( - &mut matrix, - &HelmholtzGreenKernel { k: 2.5 }, - false, - false, - &space0, - &space1, - ); - - let mut matrix2 = Array2D::>::new(( - space1.dofmap().global_size(), - space0.dofmap().global_size(), - )); - - assemble_dense( - &mut matrix2, - BoundaryOperator::SingleLayer, - PDEType::Helmholtz(2.5), - &space0, - &space1, - ); - - for i in 0..space1.dofmap().global_size() { - for j in 0..space0.dofmap().global_size() { - assert_relative_eq!( - matrix.get(i, j).unwrap().re, - matrix2.get(i, j).unwrap().re, - epsilon = 0.0001 - ); - assert_relative_eq!( - matrix.get(i, j).unwrap().im, - matrix2.get(i, j).unwrap().im, - epsilon = 0.0001 - ); - } - } - } } diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index 6d7c913e..ff9ecbf8 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -1,20 +1,19 @@ use crate::function_space::SerialFunctionSpace; -use crate::green::{ - //HelmholtzGreenHypersingularTermKernel, HelmholtzGreenKernel, LaplaceGreenKernel, - Scalar, - SingularKernel, -}; +use bempp_kernel::traits::Kernel; +use bempp_kernel::types::EvalType; 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_tools::arrays::{to_matrix, zero_matrix, Array4D, Mat}; +use bempp_traits::arrays::{AdjacencyListAccess, 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 bempp_traits::types::Scalar; use rayon::prelude::*; +use rlst_dense::{RandomAccessByRef, RawAccess, RawAccessMut, Shape, UnsafeRandomAccessMut}; fn get_quadrature_rule( test_celltype: ReferenceCellType, @@ -69,37 +68,39 @@ fn get_quadrature_rule( } } -struct RawData2D { +pub struct RawData2D { pub data: *mut T, pub shape: (usize, usize), } unsafe impl Sync for RawData2D {} +// TODO: use T not f64 #[allow(clippy::too_many_arguments)] -fn assemble_batch_singular<'a, T: Scalar + Clone + Copy>( - output: &RawData2D, - kernel: &impl SingularKernel, +fn assemble_batch_singular<'a>( + output: &RawData2D, + kernel: &impl Kernel, 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, + trial_points: &Mat, + test_points: &Mat, weights: &[f64], trial_table: &Array4D, test_table: &Array4D, ) -> usize { let grid = test_space.grid(); + let mut k = vec![0.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)); + let mut test_mapped_pts = zero_matrix((test_points.shape().0, 3)); + let mut trial_mapped_pts = zero_matrix((trial_points.shape().0, 3)); + let mut test_normals = zero_matrix((test_points.shape().0, 3)); + let mut trial_normals = zero_matrix((trial_points.shape().0, 3)); for (test_cell, trial_cell) in cell_pairs { let test_cell_tindex = grid.topology().index_map()[*test_cell]; @@ -145,20 +146,25 @@ fn assemble_batch_singular<'a, T: Scalar + Clone + Copy>( .iter() .enumerate() { - let mut sum = T::zero(); + let mut sum = 0.0; 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) } + let mut test_row = vec![0.0; test_mapped_pts.shape().1]; + for (i, ti) in test_row.iter_mut().enumerate() { + *ti = *test_mapped_pts.get(index, i).unwrap(); + } + let mut trial_row = vec![0.0; trial_mapped_pts.shape().1]; + for (i, ti) in trial_row.iter_mut().enumerate() { + *ti = *trial_mapped_pts.get(index, i).unwrap(); + } + + kernel.assemble_st(EvalType::Value, &test_row, &trial_row, &mut k); + sum += k[0] + * (wt + * test_table.get(0, index, test_i, 0).unwrap() * test_jdet[index] - * unsafe { trial_table.get_unchecked(0, index, trial_i, 0) } - * trial_jdet[index], - ); + * trial_table.get(0, index, trial_i, 0).unwrap() + * trial_jdet[index]); } unsafe { *output.data.offset( @@ -174,18 +180,18 @@ fn assemble_batch_singular<'a, T: Scalar + Clone + Copy>( } #[allow(clippy::too_many_arguments)] -fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( - output: &RawData2D, - kernel: &impl SingularKernel, +fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize>( + output: &RawData2D, + kernel: &impl Kernel, 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_points: &Mat, trial_weights: &[f64], - test_points: &Array2D, + test_points: &Mat, test_weights: &[f64], trial_table: &Array4D, test_table: &Array4D, @@ -195,27 +201,58 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( 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)); + let mut k = vec![0.0; NPTS_TEST * NPTS_TRIAL]; + let mut test_jdet = vec![0.0; NPTS_TEST]; + let mut trial_jdet = vec![0.0; NPTS_TRIAL]; + let mut test_normals = zero_matrix((NPTS_TEST, 3)); + let mut trial_normals = zero_matrix((NPTS_TRIAL, 3)); + + // let mut rlst_test_normals = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; + // let mut rlst_trial_normals = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; + let mut rlst_test_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; + let mut rlst_trial_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TRIAL, 3)]; + + let test_element = test_grid.geometry().element(test_cells[0]); + + let mut test_data = Array4D::::new(test_element.tabulate_array_shape(0, NPTS_TEST)); + test_element.tabulate(test_points, 0, &mut test_data); + let gdim = test_grid.geometry().dim(); + + // TODO: move this to grid.get_compute_points_function(test_points) + let test_compute_points = |cell: usize, pts: &mut Mat| { + for p in 0..NPTS_TEST { + for i in 0..gdim { + unsafe { + *pts.get_unchecked_mut(p, i) = 0.0; + } + } + } + let vertices = test_grid.geometry().cell_vertices(cell).unwrap(); + for (i, n) in vertices.iter().enumerate() { + let pt = test_grid.geometry().point(*n).unwrap(); + for p in 0..NPTS_TEST { + for (j, pt_j) in pt.iter().enumerate() { + unsafe { + *pts.get_unchecked_mut(p, j) += + *pt_j * *test_data.get_unchecked(0, p, i, 0); + } + } + } + } + }; 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(); + let test_vertices = unsafe { test_c20.row_unchecked(test_cell_tindex) }; 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); + test_compute_points(test_cell_gindex, &mut rlst_test_mapped_pts); + if needs_test_normal { test_grid .geometry() @@ -225,7 +262,7 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( 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(); + let trial_vertices = unsafe { trial_c20.row_unchecked(trial_cell_tindex) }; trial_grid.geometry().compute_jacobian_determinants( trial_points, @@ -235,7 +272,7 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( trial_grid.geometry().compute_points( trial_points, trial_cell_gindex, - &mut trial_mapped_pts, + &mut rlst_trial_mapped_pts, ); if needs_trial_normal { trial_grid.geometry().compute_normals( @@ -245,6 +282,13 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( ); } + kernel.assemble_st( + EvalType::Value, + rlst_test_mapped_pts.data(), + rlst_trial_mapped_pts.data(), + &mut k, + ); + for (test_i, test_dof) in test_space .dofmap() .cell_dofs(test_cell_tindex) @@ -259,25 +303,17 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( .iter() .enumerate() { - let mut sum = T::zero(); + let mut sum = 0.0; 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 + sum += k[test_index * trial_weights.len() + trial_index] + * (test_wt * trial_wt - * unsafe { test_table.get_unchecked(0, test_index, test_i, 0) } + * test_table.get(0, test_index, test_i, 0).unwrap() * test_jdet[test_index] - * unsafe { - trial_table.get_unchecked(0, trial_index, trial_i, 0) - } - * trial_jdet[test_index], - ); + * trial_table.get(0, trial_index, trial_i, 0).unwrap() + * trial_jdet[test_index]); } } // TODO: should we write into a result array, then copy into output after this loop? @@ -304,19 +340,55 @@ fn assemble_batch_nonadjacent<'a, T: Scalar + Clone + Copy>( 1 } -pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( - output: &mut Array2D, - kernel: &impl SingularKernel, +pub fn assemble<'a>( + output: &mut Mat, + kernel: &impl Kernel, 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"); - // } + let test_colouring = test_space.compute_cell_colouring(); + let trial_colouring = trial_space.compute_cell_colouring(); + // TODO: make these configurable + let blocksize = 128; + + assemble_nonsingular::<16, 16>( + output, + kernel, + needs_trial_normal, + needs_test_normal, + trial_space, + test_space, + &trial_colouring, + &test_colouring, + blocksize, + ); + assemble_singular( + output, + kernel, + needs_trial_normal, + needs_test_normal, + trial_space, + test_space, + &trial_colouring, + &test_colouring, + blocksize, + ); +} + +#[allow(clippy::too_many_arguments)] +pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize>( + output: &mut Mat, + kernel: &impl Kernel, + needs_trial_normal: bool, + needs_test_normal: bool, + trial_space: &SerialFunctionSpace<'a>, + test_space: &SerialFunctionSpace<'a>, + trial_colouring: &Vec>, + test_colouring: &Vec>, + blocksize: usize, +) { if !trial_space.is_serial() || !test_space.is_serial() { panic!("Dense assemble can only be used for function spaces stored in serial"); } @@ -326,46 +398,37 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( 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 qrule_test = simplex_rule(ReferenceCellType::Triangle, NPTS_TEST).unwrap(); + let qpoints_test = to_matrix(&qrule_test.points, (NPTS_TEST, 2)); + let qweights_test = qrule_test.weights; + let qrule_trial = simplex_rule(ReferenceCellType::Triangle, NPTS_TRIAL).unwrap(); + let qpoints_trial = to_matrix(&qrule_trial.points, (NPTS_TRIAL, 2)); + let qweights_trial = qrule_trial.weights; + + let mut test_table = + Array4D::::new(test_space.element().tabulate_array_shape(0, NPTS_TEST)); + test_space + .element() + .tabulate(&qpoints_test, 0, &mut test_table); - let mut trial_table = Array4D::::new( - trial_space - .element() - .tabulate_array_shape(0, qpoints.shape().0), - ); + let mut trial_table = + Array4D::::new(trial_space.element().tabulate_array_shape(0, NPTS_TRIAL)); trial_space .element() - .tabulate(&qpoints, 0, &mut trial_table); + .tabulate(&qpoints_test, 0, &mut trial_table); let output_raw = RawData2D { - data: output.data.as_mut_ptr(), - shape: *output.shape(), + data: output.data_mut().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 { + 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![]; @@ -391,11 +454,11 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( test_start = test_end } - let numthreads = test_cells.len(); - let r: usize = (0..numthreads) + let numtasks = test_cells.len(); + let r: usize = (0..numtasks) .into_par_iter() .map(&|t| { - assemble_batch_nonadjacent( + assemble_batch_nonadjacent::( &output_raw, kernel, needs_trial_normal, @@ -404,27 +467,58 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( trial_cells[t], test_space, test_cells[t], - &qpoints, - &qweights, - &qpoints, - &qweights, + &qpoints_trial, + &qweights_trial, + &qpoints_test, + &qweights_test, &trial_table, &test_table, ) }) .sum(); - assert_eq!(r, numthreads); + assert_eq!(r, numtasks); } } +} +#[allow(clippy::too_many_arguments)] +pub fn assemble_singular<'a>( + output: &mut Mat, + kernel: &impl Kernel, + needs_trial_normal: bool, + needs_test_normal: bool, + trial_space: &SerialFunctionSpace<'a>, + test_space: &SerialFunctionSpace<'a>, + trial_colouring: &Vec>, + test_colouring: &Vec>, + blocksize: usize, +) { if test_space.grid() != trial_space.grid() { // If the test and trial grids are different, there are no neighbouring triangles return; } + 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"); + } + + // 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 = 4; + let output_raw = RawData2D { + data: output.data_mut().as_mut_ptr(), + shape: output.shape(), + }; + let grid = test_space.grid(); let c20 = grid.topology().connectivity(2, 0); @@ -463,7 +557,7 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( npoints, ); - let points = Array2D::from_data(qrule.trial_points, (qrule.npoints, 2)); + let points = to_matrix(&qrule.trial_points, (qrule.npoints, 2)); let mut table = Array4D::::new( trial_space .element() @@ -473,7 +567,7 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( trial_points.push(points); trial_tables.push(table); - let points = Array2D::from_data(qrule.test_points, (qrule.npoints, 2)); + let points = to_matrix(&qrule.test_points, (qrule.npoints, 2)); let mut table = Array4D::::new( test_space .element() @@ -485,8 +579,8 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( qweights.push(qrule.weights); } - for test_c in &test_colouring { - for trial_c in &trial_colouring { + 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]; @@ -522,8 +616,8 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( start = end; } - let numthreads = cell_blocks.len(); - let r: usize = (0..numthreads) + let numtasks = cell_blocks.len(); + let r: usize = (0..numtasks) .into_par_iter() .map(&|t| { assemble_batch_singular( @@ -542,7 +636,7 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( ) }) .sum(); - assert_eq!(r, numthreads); + assert_eq!(r, numtasks); } } } @@ -550,20 +644,17 @@ pub fn assemble<'a, T: Scalar + Clone + Copy + Sync>( #[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_kernel::laplace_3d::Laplace3dKernel; 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); + fn test_laplace_single_layer_dp0_dp0() { + let grid = regular_sphere(0); let element = create_element( ElementFamily::Lagrange, ReferenceCellType::Triangle, @@ -574,39 +665,147 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = zero_matrix::((ndofs, ndofs)); assemble( &mut matrix, - &green::LaplaceGreenKernel {}, + &Laplace3dKernel::new(), false, false, &space, &space, ); - let mut dmat = Array2D::::new((ndofs, ndofs)); - dense::assemble( - &mut dmat, - &green::LaplaceGreenKernel {}, + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-3); + } + } + } + /* + + #[test] + fn test_laplace_double_layer_dp0_dp0() { + let grid = regular_sphere(0); + 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::LaplaceGreenDyKernel {}, + true, false, + &space, + &space, + ); + + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); + } + } + } + + #[test] + fn test_laplace_adjoint_double_layer_dp0_dp0() { + let grid = regular_sphere(0); + 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::LaplaceGreenDxKernel {}, false, + true, &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(), - ); + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); } - }*/ - for i in 0..matrix.shape().0 { - for j in 0..matrix.shape().1 { + } + } + + #[test] + fn test_laplace_hypersingular_dp0_dp0() { + let grid = regular_sphere(0); + 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)); + laplace_hypersingular_assemble(&mut matrix, &space, &space); + + for i in 0..ndofs { + for j in 0..ndofs { + assert_relative_eq!(*matrix.get(i, j).unwrap(), 0.0, epsilon = 1e-4); + } + } + } + + #[test] + fn test_laplace_hypersingular_p1_p1() { + let grid = regular_sphere(0); + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 1, + Continuity::Continuous, + ); + let space = SerialFunctionSpace::new(&grid, &element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::::new((ndofs, ndofs)); + + laplace_hypersingular_assemble(&mut matrix, &space, &space); + + // Compare to result from bempp-cl + #[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]; + + for (i, pi) in perm.iter().enumerate() { + for (j, pj) in perm.iter().enumerate() { assert_relative_eq!( *matrix.get(i, j).unwrap(), - *dmat.get(i, j).unwrap(), + from_cl[*pi][*pj], epsilon = 1e-4 ); } @@ -614,8 +813,8 @@ mod test { } #[test] - fn test_laplace_single_layer_dp0_dp0_smaller() { - let grid = regular_sphere(1); + fn test_helmholtz_single_layer_real_dp0_dp0() { + let grid = regular_sphere(0); let element = create_element( ElementFamily::Lagrange, ReferenceCellType::Triangle, @@ -629,31 +828,164 @@ mod test { let mut matrix = Array2D::::new((ndofs, ndofs)); assemble( &mut matrix, - &green::LaplaceGreenKernel {}, + &green::HelmholtzGreenKernel { k: 3.0 }, false, false, &space, &space, ); - let mut dmat = Array2D::::new((ndofs, ndofs)); - dense::assemble( - &mut dmat, - &green::LaplaceGreenKernel {}, + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); + } + } + } + #[test] + fn test_helmholtz_single_layer_complex_dp0_dp0() { + let grid = regular_sphere(0); + 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::HelmholtzGreenKernel { k: 3.0 }, false, false, &space, &space, ); - for i in 0..matrix.shape().0 { - for j in 0..matrix.shape().1 { + // Compare to result from bempp-cl + #[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); + assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); + } + } + } + + #[test] + fn test_helmholtz_double_layer_dp0_dp0() { + let grid = regular_sphere(0); + 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::HelmholtzGreenDyKernel { k: 3.0 }, + true, + false, + &space, + &space, + ); + + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(matrix.get(i, j).unwrap().re, entry.re, epsilon = 1e-4); + assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); + } + } + } + + #[test] + fn test_helmholtz_adjoint_double_layer_dp0_dp0() { + let grid = regular_sphere(0); + 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::HelmholtzGreenDxKernel { k: 3.0 }, + false, + true, + &space, + &space, + ); + + // Compare to result from bempp-cl + #[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() { + assert_relative_eq!(matrix.get(i, j).unwrap().re, entry.re, epsilon = 1e-4); + assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); + } + } + } + + #[test] + fn test_helmholtz_hypersingular_p1_p1() { + let grid = regular_sphere(0); + let element = create_element( + ElementFamily::Lagrange, + ReferenceCellType::Triangle, + 1, + Continuity::Continuous, + ); + let space = SerialFunctionSpace::new(&grid, &element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::>::new((ndofs, ndofs)); + + helmholtz_hypersingular_assemble(&mut matrix, &space, &space, 3.0); + + // Compare to result from bempp-cl + #[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]; + + for (i, pi) in perm.iter().enumerate() { + for (j, pj) in perm.iter().enumerate() { assert_relative_eq!( - *matrix.get(i, j).unwrap(), - *dmat.get(i, j).unwrap(), + matrix.get(i, j).unwrap().re, + from_cl[*pi][*pj].re, + epsilon = 1e-3 + ); + assert_relative_eq!( + matrix.get(i, j).unwrap().im, + from_cl[*pi][*pj].im, epsilon = 1e-3 ); } } } + */ } diff --git a/bem/src/assembly/dense.rs b/bem/src/assembly/dense.rs deleted file mode 100644 index 888239ea..00000000 --- a/bem/src/assembly/dense.rs +++ /dev/null @@ -1,813 +0,0 @@ -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::{available_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}; - -fn get_quadrature_rule( - test_celltype: ReferenceCellType, - trial_celltype: ReferenceCellType, - pairs: Vec<(usize, usize)>, - npoints: usize, -) -> TestTrialNumericalQuadratureDefinition { - if pairs.is_empty() { - // Standard rules - let mut npoints_test = 10 * npoints * npoints; - for p in available_rules(test_celltype) { - if p >= npoints * npoints && p < npoints_test { - npoints_test = p; - } - } - let mut npoints_trial = 10 * npoints * npoints; - for p in available_rules(trial_celltype) { - if p >= npoints * npoints && p < npoints_trial { - npoints_trial = p; - } - } - let test_rule = simplex_rule(test_celltype, npoints_test).unwrap(); - let trial_rule = simplex_rule(trial_celltype, npoints_trial).unwrap(); - if test_rule.dim != trial_rule.dim { - unimplemented!("Quadrature with different dimension cells not supported"); - } - if test_rule.order != trial_rule.order { - unimplemented!("Quadrature with different trial and test orders not supported"); - } - let dim = test_rule.dim; - let npts = test_rule.npoints * trial_rule.npoints; - let mut test_points = Vec::::with_capacity(dim * npts); - let mut trial_points = Vec::::with_capacity(dim * npts); - let mut weights = Vec::::with_capacity(npts); - - for test_i in 0..test_rule.npoints { - for trial_i in 0..trial_rule.npoints { - for d in 0..dim { - test_points.push(test_rule.points[dim * test_i + d]); - trial_points.push(trial_rule.points[dim * trial_i + d]); - } - weights.push(test_rule.weights[test_i] * trial_rule.weights[trial_i]); - } - } - - TestTrialNumericalQuadratureDefinition { - dim, - order: test_rule.order, - npoints: npts, - weights, - test_points, - trial_points, - } - } 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, - }, - 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, - }, - npoints, - ) - .unwrap() - } - } -} - -pub fn assemble<'a, T: Scalar>( - output: &mut Array2D, - kernel: &impl SingularKernel, - needs_trial_normal: bool, - needs_test_normal: bool, - trial_space: &impl FunctionSpace<'a>, - test_space: &impl FunctionSpace<'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: allow user to configure this - let npoints = 4; - - let grid = trial_space.grid(); - let c20 = grid.topology().connectivity(2, 0); - - for test_cell in 0..grid.geometry().cell_count() { - let test_cell_tindex = grid.topology().index_map()[test_cell]; - let test_cell_gindex = grid.geometry().index_map()[test_cell]; - let test_vertices = c20.row(test_cell_tindex).unwrap(); - - 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 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)); - } - } - } - let rule = get_quadrature_rule( - grid.topology().cell_type(test_cell_tindex).unwrap(), - grid.topology().cell_type(trial_cell_tindex).unwrap(), - pairs, - npoints, - ); - - let test_points = Array2D::from_data(rule.test_points, (rule.npoints, 2)); - let trial_points = Array2D::from_data(rule.trial_points, (rule.npoints, 2)); - let mut test_table = - Array4D::::new(test_space.element().tabulate_array_shape(0, rule.npoints)); - let mut trial_table = - Array4D::::new(trial_space.element().tabulate_array_shape(0, rule.npoints)); - - test_space - .element() - .tabulate(&test_points, 0, &mut test_table); - trial_space - .element() - .tabulate(&trial_points, 0, &mut trial_table); - - let mut test_jdet = vec![0.0; rule.npoints]; - let mut trial_jdet = vec![0.0; rule.npoints]; - - grid.geometry().compute_jacobian_determinants( - &test_points, - test_cell_gindex, - &mut test_jdet, - ); - grid.geometry().compute_jacobian_determinants( - &trial_points, - trial_cell_gindex, - &mut trial_jdet, - ); - - let mut test_mapped_pts = Array2D::::new((rule.npoints, 3)); - let mut trial_mapped_pts = Array2D::::new((rule.npoints, 3)); - let mut test_normals = Array2D::::new((rule.npoints, 3)); - let mut trial_normals = Array2D::::new((rule.npoints, 3)); - - grid.geometry() - .compute_points(&test_points, test_cell_gindex, &mut test_mapped_pts); - grid.geometry() - .compute_points(&trial_points, trial_cell_gindex, &mut trial_mapped_pts); - if needs_test_normal { - grid.geometry() - .compute_normals(&test_points, test_cell_gindex, &mut test_normals); - } - 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 in 0..rule.npoints { - 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( - rule.weights[index] - * 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], - ); - } - *output.get_mut(*test_dof, *trial_dof).unwrap() += sum; - } - } - } - } -} - -pub fn curl_curl_assemble<'a, T: Scalar>( - output: &mut Array2D, - kernel: &impl SingularKernel, - trial_space: &impl FunctionSpace<'a>, - test_space: &impl FunctionSpace<'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"); - } - - let npoints = 4; - - let grid = trial_space.grid(); - let c20 = grid.topology().connectivity(2, 0); - - for test_cell in 0..grid.geometry().cell_count() { - let test_cell_tindex = grid.topology().index_map()[test_cell]; - let test_cell_gindex = grid.geometry().index_map()[test_cell]; - let test_vertices = c20.row(test_cell_tindex).unwrap(); - - 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 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)); - } - } - } - let rule = get_quadrature_rule( - grid.topology().cell_type(test_cell_tindex).unwrap(), - grid.topology().cell_type(trial_cell_tindex).unwrap(), - pairs, - npoints, - ); - let test_points = Array2D::from_data(rule.test_points, (rule.npoints, 2)); - let trial_points = Array2D::from_data(rule.trial_points, (rule.npoints, 2)); - let mut test_table = - Array4D::::new(test_space.element().tabulate_array_shape(1, rule.npoints)); - let mut trial_table = - Array4D::::new(trial_space.element().tabulate_array_shape(1, rule.npoints)); - - test_space - .element() - .tabulate(&test_points, 1, &mut test_table); - trial_space - .element() - .tabulate(&trial_points, 1, &mut trial_table); - - let mut test_jdet = vec![0.0; rule.npoints]; - let mut trial_jdet = vec![0.0; rule.npoints]; - let mut test_jinv = Array2D::::new((rule.npoints, 6)); - let mut trial_jinv = Array2D::::new((rule.npoints, 6)); - let mut test_mapped_pts = Array2D::::new((rule.npoints, 3)); - let mut trial_mapped_pts = Array2D::::new((rule.npoints, 3)); - let mut test_normals = Array2D::::new((rule.npoints, 3)); - let mut trial_normals = Array2D::::new((rule.npoints, 3)); - - grid.geometry().compute_jacobian_determinants( - &test_points, - test_cell_gindex, - &mut test_jdet, - ); - grid.geometry().compute_jacobian_determinants( - &trial_points, - trial_cell_gindex, - &mut trial_jdet, - ); - grid.geometry().compute_jacobian_inverses( - &test_points, - test_cell_gindex, - &mut test_jinv, - ); - grid.geometry().compute_jacobian_inverses( - &trial_points, - trial_cell_gindex, - &mut trial_jinv, - ); - grid.geometry() - .compute_points(&test_points, test_cell_gindex, &mut test_mapped_pts); - grid.geometry() - .compute_points(&trial_points, trial_cell_gindex, &mut trial_mapped_pts); - grid.geometry() - .compute_normals(&test_points, test_cell_gindex, &mut test_normals); - 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 in 0..rule.npoints { - let g0 = ( - unsafe { - *trial_jinv.get_unchecked(index, 0) - * *trial_table.get_unchecked(1, index, trial_i, 0) - + *trial_jinv.get_unchecked(index, 3) - * *trial_table.get_unchecked(2, index, trial_i, 0) - }, - unsafe { - *trial_jinv.get_unchecked(index, 1) - * *trial_table.get_unchecked(1, index, trial_i, 0) - + *trial_jinv.get_unchecked(index, 4) - * *trial_table.get_unchecked(2, index, trial_i, 0) - }, - unsafe { - *trial_jinv.get_unchecked(index, 2) - * *trial_table.get_unchecked(1, index, trial_i, 0) - + *trial_jinv.get_unchecked(index, 5) - * *trial_table.get_unchecked(2, index, trial_i, 0) - }, - ); - let g1 = ( - unsafe { - *test_jinv.get_unchecked(index, 0) - * *test_table.get_unchecked(1, index, test_i, 0) - + *test_jinv.get_unchecked(index, 3) - * *test_table.get_unchecked(2, index, test_i, 0) - }, - unsafe { - *test_jinv.get_unchecked(index, 1) - * *test_table.get_unchecked(1, index, test_i, 0) - + *test_jinv.get_unchecked(index, 4) - * *test_table.get_unchecked(2, index, test_i, 0) - }, - unsafe { - *test_jinv.get_unchecked(index, 2) - * *test_table.get_unchecked(1, index, test_i, 0) - + *test_jinv.get_unchecked(index, 5) - * *test_table.get_unchecked(2, index, test_i, 0) - }, - ); - let n0 = ( - unsafe { *trial_normals.get_unchecked(index, 0) }, - unsafe { *trial_normals.get_unchecked(index, 1) }, - unsafe { *trial_normals.get_unchecked(index, 2) }, - ); - let n1 = ( - unsafe { *test_normals.get_unchecked(index, 0) }, - unsafe { *test_normals.get_unchecked(index, 1) }, - unsafe { *test_normals.get_unchecked(index, 2) }, - ); - - let dot_curls = (g0.0 * g1.0 + g0.1 * g1.1 + g0.2 * g1.2) - * (n0.0 * n1.0 + n0.1 * n1.1 + n0.2 * n1.2) - - (g0.0 * n1.0 + g0.1 * n1.1 + g0.2 * n1.2) - * (n0.0 * g1.0 + n0.1 * g1.1 + n0.2 * g1.2); - - 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( - rule.weights[index] * dot_curls * test_jdet[index] * trial_jdet[index], - ); - } - *output.get_mut(*test_dof, *trial_dof).unwrap() += sum; - } - } - } - } -} - -pub fn laplace_hypersingular_assemble<'a, T: Scalar>( - output: &mut Array2D, - trial_space: &impl FunctionSpace<'a>, - test_space: &impl FunctionSpace<'a>, -) { - curl_curl_assemble(output, &LaplaceGreenKernel {}, trial_space, test_space); -} - -pub fn helmholtz_hypersingular_assemble<'a, T: Scalar>( - output: &mut Array2D, - trial_space: &impl FunctionSpace<'a>, - test_space: &impl FunctionSpace<'a>, - k: f64, -) { - curl_curl_assemble(output, &HelmholtzGreenKernel { k }, trial_space, test_space); - assemble( - output, - &HelmholtzGreenHypersingularTermKernel { k }, - true, - true, - trial_space, - test_space, - ); -} - -#[cfg(test)] -mod test { - 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; - - #[test] - fn test_laplace_single_layer_dp0_dp0() { - let grid = regular_sphere(0); - 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, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); - } - } - } - - #[test] - fn test_laplace_double_layer_dp0_dp0() { - let grid = regular_sphere(0); - 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::LaplaceGreenDyKernel {}, - true, - false, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); - } - } - } - - #[test] - fn test_laplace_adjoint_double_layer_dp0_dp0() { - let grid = regular_sphere(0); - 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::LaplaceGreenDxKernel {}, - false, - true, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); - } - } - } - - #[test] - fn test_laplace_hypersingular_dp0_dp0() { - let grid = regular_sphere(0); - 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)); - laplace_hypersingular_assemble(&mut matrix, &space, &space); - - for i in 0..ndofs { - for j in 0..ndofs { - assert_relative_eq!(*matrix.get(i, j).unwrap(), 0.0, epsilon = 1e-4); - } - } - } - - #[test] - fn test_laplace_hypersingular_p1_p1() { - let grid = regular_sphere(0); - let element = create_element( - ElementFamily::Lagrange, - ReferenceCellType::Triangle, - 1, - Continuity::Continuous, - ); - let space = SerialFunctionSpace::new(&grid, &element); - - let ndofs = space.dofmap().global_size(); - - let mut matrix = Array2D::::new((ndofs, ndofs)); - - laplace_hypersingular_assemble(&mut matrix, &space, &space); - - // Compare to result from bempp-cl - #[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]; - - for (i, pi) in perm.iter().enumerate() { - for (j, pj) in perm.iter().enumerate() { - assert_relative_eq!( - *matrix.get(i, j).unwrap(), - from_cl[*pi][*pj], - epsilon = 1e-4 - ); - } - } - } - - #[test] - fn test_helmholtz_single_layer_real_dp0_dp0() { - let grid = regular_sphere(0); - 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::HelmholtzGreenKernel { k: 3.0 }, - false, - false, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); - } - } - } - #[test] - fn test_helmholtz_single_layer_complex_dp0_dp0() { - let grid = regular_sphere(0); - 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::HelmholtzGreenKernel { k: 3.0 }, - false, - false, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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); - assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); - } - } - } - - #[test] - fn test_helmholtz_double_layer_dp0_dp0() { - let grid = regular_sphere(0); - 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::HelmholtzGreenDyKernel { k: 3.0 }, - true, - false, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(matrix.get(i, j).unwrap().re, entry.re, epsilon = 1e-4); - assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); - } - } - } - - #[test] - fn test_helmholtz_adjoint_double_layer_dp0_dp0() { - let grid = regular_sphere(0); - 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::HelmholtzGreenDxKernel { k: 3.0 }, - false, - true, - &space, - &space, - ); - - // Compare to result from bempp-cl - #[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() { - assert_relative_eq!(matrix.get(i, j).unwrap().re, entry.re, epsilon = 1e-4); - assert_relative_eq!(matrix.get(i, j).unwrap().im, entry.im, epsilon = 1e-4); - } - } - } - - #[test] - fn test_helmholtz_hypersingular_p1_p1() { - let grid = regular_sphere(0); - let element = create_element( - ElementFamily::Lagrange, - ReferenceCellType::Triangle, - 1, - Continuity::Continuous, - ); - let space = SerialFunctionSpace::new(&grid, &element); - - let ndofs = space.dofmap().global_size(); - - let mut matrix = Array2D::>::new((ndofs, ndofs)); - - helmholtz_hypersingular_assemble(&mut matrix, &space, &space, 3.0); - - // Compare to result from bempp-cl - #[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]; - - for (i, pi) in perm.iter().enumerate() { - for (j, pj) in perm.iter().enumerate() { - assert_relative_eq!( - matrix.get(i, j).unwrap().re, - from_cl[*pi][*pj].re, - epsilon = 1e-3 - ); - assert_relative_eq!( - matrix.get(i, j).unwrap().im, - from_cl[*pi][*pj].im, - epsilon = 1e-3 - ); - } - } - } -} diff --git a/bem/src/dofmap.rs b/bem/src/dofmap.rs index 7c8f50ca..726ecfd3 100644 --- a/bem/src/dofmap.rs +++ b/bem/src/dofmap.rs @@ -1,12 +1,12 @@ -use bempp_tools::arrays::{AdjacencyList, Array2D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess}; +use bempp_tools::arrays::AdjacencyList; +use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::bem::DofMap; use bempp_traits::element::FiniteElement; use bempp_traits::grid::{Grid, Topology}; pub struct SerialDofMap { entity_dofs: [AdjacencyList; 4], - cell_dofs: Array2D, + cell_dofs: AdjacencyList, size: usize, } @@ -18,8 +18,14 @@ impl SerialDofMap { for d in 0..tdim + 1 { entity_dofs_data.push(vec![vec![]; grid.topology().entity_count(d)]); } - let mut cell_dofs = - Array2D::::new((grid.topology().entity_count(tdim), element.dim())); + let mut offsets = vec![]; + for i in 0..grid.topology().entity_count(tdim) + 1 { + offsets.push(i * element.dim()); + } + let mut cell_dofs = AdjacencyList::from_data( + vec![0; grid.topology().entity_count(tdim) * element.dim()], + offsets, + ); for cell in 0..grid.topology().entity_count(tdim) { for (d, ed_data) in entity_dofs_data.iter_mut().enumerate() { for (i, e) in unsafe { diff --git a/bem/src/green.rs b/bem/src/green.rs deleted file mode 100644 index 2f3d215b..00000000 --- a/bem/src/green.rs +++ /dev/null @@ -1,175 +0,0 @@ -use num::complex::Complex; -use num::Num; - -pub trait Scalar: Num + std::ops::AddAssign { - /// Get 1 over 4*pi as this scalar type - fn inv_4pi() -> Self; - /// Get -1 over 4*pi as this scalar type - fn neg_inv_4pi() -> Self; - /// Get the distance between x and y as this scalar type - fn dist(x: &[f64], y: &[f64]) -> Self; - /// Get the square of the distance between x and y as this scalar type - fn dist_squared(x: &[f64], y: &[f64]) -> Self; - /// Get the cube of the distance between x and y as this scalar type - fn dist_cubed(x: &[f64], y: &[f64]) -> Self; - /// Get x.y as this scalar type - fn dot(x: &[f64], y: &[f64]) -> Self; - /// Get (x-y).n as this scalar type - fn subdot(x: &[f64], y: &[f64], n: &[f64]) -> Self; - /// Convert a f64 to this type - fn from_f64(v: f64) -> Self; - /// Get e^(i*x) as this scalar type - fn eix(x: f64) -> Self; - /// Get i*e^(i*x) as this scalar type - fn ieix(x: f64) -> Self; -} - -impl Scalar for f64 { - fn inv_4pi() -> Self { - 0.25 * std::f64::consts::FRAC_1_PI - } - fn neg_inv_4pi() -> Self { - -0.25 * std::f64::consts::FRAC_1_PI - } - fn dist(x: &[f64], y: &[f64]) -> Self { - f64::sqrt( - (x[0] - y[0]) * (x[0] - y[0]) - + (x[1] - y[1]) * (x[1] - y[1]) - + (x[2] - y[2]) * (x[2] - y[2]), - ) - } - fn dist_squared(x: &[f64], y: &[f64]) -> Self { - (x[0] - y[0]) * (x[0] - y[0]) - + (x[1] - y[1]) * (x[1] - y[1]) - + (x[2] - y[2]) * (x[2] - y[2]) - } - fn dist_cubed(x: &[f64], y: &[f64]) -> Self { - let sq = Self::dist_squared(x, y); - sq * f64::sqrt(sq) - } - fn dot(x: &[f64], y: &[f64]) -> Self { - x[0] * y[0] + x[1] * y[1] + x[2] * y[2] - } - fn subdot(x: &[f64], y: &[f64], n: &[f64]) -> Self { - (x[0] - y[0]) * n[0] + (x[1] - y[1]) * n[1] + (x[2] - y[2]) * n[2] - } - fn from_f64(v: f64) -> Self { - v - } - fn eix(x: f64) -> Self { - Complex::::eix(x).re - } - fn ieix(x: f64) -> Self { - -Complex::::eix(x).im - } -} - -impl Scalar for Complex { - fn inv_4pi() -> Self { - Complex::::new(f64::inv_4pi(), 0.0) - } - fn neg_inv_4pi() -> Self { - Complex::::new(f64::neg_inv_4pi(), 0.0) - } - fn dist(x: &[f64], y: &[f64]) -> Self { - Complex::::new(f64::dist(x, y), 0.0) - } - fn dist_squared(x: &[f64], y: &[f64]) -> Self { - Complex::::new(f64::dist_squared(x, y), 0.0) - } - fn dist_cubed(x: &[f64], y: &[f64]) -> Self { - Complex::::new(f64::dist_cubed(x, y), 0.0) - } - fn dot(x: &[f64], y: &[f64]) -> Self { - Complex::::new(f64::dot(x, y), 0.0) - } - fn subdot(x: &[f64], y: &[f64], n: &[f64]) -> Self { - Complex::::new(f64::subdot(x, y, n), 0.0) - } - fn eix(x: f64) -> Self { - Complex::::new(0.0, x).exp() - } - fn ieix(x: f64) -> Self { - Complex::::new(0.0, 1.0) * Complex::::new(0.0, x).exp() - } - fn from_f64(v: f64) -> Self { - Complex::::new(v, 0.0) - } -} - -pub trait SingularKernel: Sync { - fn eval(&self, x: &[f64], y: &[f64], nx: &[f64], ny: &[f64]) -> T; -} - -pub struct LaplaceGreenKernel {} -impl SingularKernel for LaplaceGreenKernel { - fn eval(&self, x: &[f64], y: &[f64], _nx: &[f64], _ny: &[f64]) -> T { - T::inv_4pi() / T::dist(x, y) - } -} - -pub struct LaplaceGreenDxKernel {} -impl SingularKernel for LaplaceGreenDxKernel { - fn eval(&self, x: &[f64], y: &[f64], nx: &[f64], _ny: &[f64]) -> T { - T::inv_4pi() * T::subdot(y, x, nx) / T::dist_cubed(x, y) - } -} - -pub struct LaplaceGreenDyKernel {} -impl SingularKernel for LaplaceGreenDyKernel { - fn eval(&self, x: &[f64], y: &[f64], _nx: &[f64], ny: &[f64]) -> T { - T::inv_4pi() * T::subdot(x, y, ny) / T::dist_cubed(x, y) - } -} - -pub struct HelmholtzGreenKernel { - pub k: f64, -} -impl SingularKernel for HelmholtzGreenKernel { - fn eval(&self, x: &[f64], y: &[f64], _nx: &[f64], _ny: &[f64]) -> T { - let dist = f64::dist(x, y); - T::inv_4pi() * T::eix(self.k * dist) / T::from_f64(dist) - } -} - -pub struct HelmholtzGreenDxKernel { - pub k: f64, -} -impl SingularKernel for HelmholtzGreenDxKernel { - fn eval(&self, x: &[f64], y: &[f64], nx: &[f64], _ny: &[f64]) -> T { - let sq = f64::dist_squared(x, y); - let dist = sq.sqrt(); - T::inv_4pi() - * T::subdot(x, y, nx) - * (T::from_f64(self.k) * T::ieix(self.k * dist) - - T::eix(self.k * dist) / T::from_f64(dist)) - / T::from_f64(sq) - } -} - -pub struct HelmholtzGreenDyKernel { - pub k: f64, -} -impl SingularKernel for HelmholtzGreenDyKernel { - fn eval(&self, x: &[f64], y: &[f64], _nx: &[f64], ny: &[f64]) -> T { - let sq = f64::dist_squared(x, y); - let dist = sq.sqrt(); - T::inv_4pi() - * T::subdot(y, x, ny) - * (T::from_f64(self.k) * T::ieix(self.k * dist) - - T::eix(self.k * dist) / T::from_f64(dist)) - / T::from_f64(sq) - } -} - -pub struct HelmholtzGreenHypersingularTermKernel { - pub k: f64, -} -impl SingularKernel for HelmholtzGreenHypersingularTermKernel { - fn eval(&self, x: &[f64], y: &[f64], nx: &[f64], ny: &[f64]) -> T { - let dist = f64::dist(x, y); - T::from_f64(self.k) * T::from_f64(self.k) * T::neg_inv_4pi() * T::eix(self.k * dist) - / T::from_f64(dist) - * T::dot(nx, ny) - } -} diff --git a/bem/src/lib.rs b/bem/src/lib.rs index a022006d..d068780c 100644 --- a/bem/src/lib.rs +++ b/bem/src/lib.rs @@ -4,4 +4,3 @@ pub mod assembly; pub mod dofmap; pub mod function_space; -pub mod green; diff --git a/element/src/element.rs b/element/src/element.rs index 503c6205..32269d33 100644 --- a/element/src/element.rs +++ b/element/src/element.rs @@ -2,13 +2,13 @@ use crate::cell::create_cell; use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials}; -use bempp_tools::arrays::{AdjacencyList, Array2D, Array3D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess}; +use bempp_tools::arrays::{AdjacencyList, Array3D, Mat}; +use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess}; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement, MapType}; use rlst_algorithms::linalg::LinAlg; use rlst_algorithms::traits::inverse::Inverse; -use rlst_dense::{RandomAccessByRef, RandomAccessMut}; +use rlst_dense::{rlst_dynamic_mat, RandomAccessByRef, RandomAccessMut, Shape}; pub mod lagrange; pub mod raviart_thomas; @@ -37,7 +37,7 @@ impl CiarletElement { degree: usize, value_shape: Vec, polynomial_coeffs: Array3D, - interpolation_points: [Vec>; 4], + interpolation_points: [Vec>; 4], interpolation_weights: [Vec>; 4], map_type: MapType, continuity: Continuity, @@ -69,12 +69,12 @@ impl CiarletElement { } let new_pts = if continuity == Continuity::Discontinuous { - let mut new_pts = [vec![], vec![], vec![], vec![]]; + let mut new_pts: [Vec>; 4] = [vec![], vec![], vec![], vec![]]; let mut pn = 0; - let mut all_pts = Array2D::::new((npts, tdim)); + let mut all_pts = rlst_dynamic_mat![f64, (npts, tdim)]; for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { for _pts in pts_i { - new_pts[i].push(Array2D::::new((0, tdim))); + new_pts[i].push(rlst_dynamic_mat![f64, (0, tdim)]); } } for pts_i in interpolation_points.iter() { @@ -240,9 +240,9 @@ impl FiniteElement for CiarletElement { fn dim(&self) -> usize { self.dim } - fn tabulate<'a>( + fn tabulate + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, nderivs: usize, data: &mut impl Array4DAccess, ) { diff --git a/element/src/element/lagrange.rs b/element/src/element/lagrange.rs index d41d0e86..e683a181 100644 --- a/element/src/element/lagrange.rs +++ b/element/src/element/lagrange.rs @@ -2,7 +2,7 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::{Array2D, Array3D}; +use bempp_tools::arrays::{to_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; @@ -29,18 +29,18 @@ pub fn create( } for d in 0..tdim { for _e in 0..cell.entity_count(d) { - x[d].push(Array2D::::new((0, tdim))); + x[d].push(to_matrix(&[], (0, tdim))); m[d].push(Array3D::::new((0, 1, 0))); } } - x[tdim].push(Array2D::::from_data(cell.midpoint(), (1, tdim))); + x[tdim].push(to_matrix(&cell.midpoint(), (1, tdim))); m[tdim].push(Array3D::::from_data(vec![1.0], (1, 1, 1))); } else { // TODO: GLL points for e in 0..cell.entity_count(0) { let mut pts = vec![0.0; tdim]; pts.copy_from_slice(&cell.vertices()[e * tdim..(e + 1) * tdim]); - x[0].push(Array2D::::from_data(pts, (1, tdim))); + x[0].push(to_matrix(&pts, (1, tdim))); m[0].push(Array3D::::from_data(vec![1.0], (1, 1, 1))); } for e in 0..cell.entity_count(1) { @@ -56,7 +56,7 @@ pub fn create( pts[(i - 1) * tdim + j] = v0[j] + i as f64 / degree as f64 * (v1[j] - v0[j]); } } - x[1].push(Array2D::::from_data(pts, (degree - 1, tdim))); + x[1].push(to_matrix(&pts, (degree - 1, tdim))); m[1].push(Array3D::::from_data( ident, (degree - 1, 1, degree - 1), @@ -119,7 +119,7 @@ pub fn create( for i in 0..npts { ident[i * npts + i] = 1.0; } - x[2].push(Array2D::::from_data(pts, (npts, tdim))); + x[2].push(to_matrix(&pts, (npts, tdim))); m[2].push(Array3D::::from_data(ident, (npts, 1, npts))); start += nvertices; } @@ -143,9 +143,10 @@ mod test { use crate::cell::*; use crate::element::lagrange::*; use approx::*; - use bempp_tools::arrays::{Array2D, Array4D}; - use bempp_traits::arrays::{Array2DAccess, Array4DAccess}; + use bempp_tools::arrays::Array4D; + use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; + use rlst_dense::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -182,7 +183,7 @@ mod test { let e = create(ReferenceCellType::Interval, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 4)); - let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1)); + let points = to_matrix(&[0.0, 0.2, 0.4, 1.0], (4, 1)); e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -196,7 +197,7 @@ mod test { let e = create(ReferenceCellType::Interval, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 4)); - let points = Array2D::from_data(vec![0.0, 0.2, 0.4, 1.0], (4, 1)); + let points = to_matrix(&[0.0, 0.2, 0.4, 1.0], (4, 1)); e.tabulate(&points, 0, &mut data); for pt in 0..4 { @@ -214,8 +215,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 0, Continuity::Discontinuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -231,8 +232,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -312,8 +313,8 @@ mod test { ); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -329,8 +330,8 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 1, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); @@ -361,8 +362,8 @@ mod test { let e = create(ReferenceCellType::Quadrilateral, 2, Continuity::Continuous); assert_eq!(e.value_size(), 1); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.25, 0.5, 0.3, 0.2], (6, 2), ); e.tabulate(&points, 0, &mut data); diff --git a/element/src/element/raviart_thomas.rs b/element/src/element/raviart_thomas.rs index 2964bc7a..2875af65 100644 --- a/element/src/element/raviart_thomas.rs +++ b/element/src/element/raviart_thomas.rs @@ -2,7 +2,7 @@ use crate::element::{create_cell, CiarletElement}; use crate::polynomials::polynomial_count; -use bempp_tools::arrays::{Array2D, Array3D}; +use bempp_tools::arrays::{to_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; @@ -45,7 +45,7 @@ pub fn create( let mut x = [vec![], vec![], vec![], vec![]]; let mut m = [vec![], vec![], vec![], vec![]]; for _e in 0..cell.entity_count(0) { - x[0].push(Array2D::::new((0, tdim))); + x[0].push(to_matrix(&[], (0, tdim))); m[0].push(Array3D::::new((0, 2, 0))); } @@ -61,12 +61,12 @@ pub fn create( } mat[0] = v0[1] - v1[1]; mat[1] = v1[0] - v0[0]; - x[1].push(Array2D::::from_data(pts, (1, tdim))); + x[1].push(to_matrix(&pts, (1, tdim))); m[1].push(Array3D::::from_data(mat, (1, 2, 1))); } for _e in 0..cell.entity_count(2) { - x[2].push(Array2D::::new((0, tdim))); + x[2].push(to_matrix(&[], (0, tdim))); m[2].push(Array3D::::new((0, 2, 0))); } @@ -89,9 +89,10 @@ mod test { use crate::cell::*; use crate::element::raviart_thomas::*; use approx::*; - use bempp_tools::arrays::{Array2D, Array4D}; - use bempp_traits::arrays::{Array2DAccess, Array4DAccess}; + use bempp_tools::arrays::Array4D; + use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; + use rlst_dense::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -128,8 +129,8 @@ mod test { let e = create(ReferenceCellType::Triangle, 1, Continuity::Continuous); assert_eq!(e.value_size(), 2); let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); - let points = Array2D::from_data( - vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], + let points = to_matrix( + &[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5], (6, 2), ); e.tabulate(&points, 0, &mut data); diff --git a/element/src/polynomials.rs b/element/src/polynomials.rs index 2c9c60c8..b4d4144e 100644 --- a/element/src/polynomials.rs +++ b/element/src/polynomials.rs @@ -1,11 +1,12 @@ //! Orthonormal polynomials -use bempp_traits::arrays::{Array2DAccess, Array3DAccess}; +use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; +use rlst_dense::{RandomAccessByRef, Shape}; /// Tabulate orthonormal polynomials on a interval -fn tabulate_legendre_polynomials_interval<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_interval + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -57,8 +58,8 @@ fn quad_index(i: usize, j: usize, n: usize) -> usize { } /// Tabulate orthonormal polynomials on a quadrilateral -fn tabulate_legendre_polynomials_quadrilateral<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_quadrilateral + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -191,8 +192,8 @@ fn tabulate_legendre_polynomials_quadrilateral<'a>( } } /// Tabulate orthonormal polynomials on a triangle -fn tabulate_legendre_polynomials_triangle<'a>( - points: &impl Array2DAccess<'a, f64>, +fn tabulate_legendre_polynomials_triangle + Shape>( + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -368,9 +369,9 @@ pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usi } } -pub fn legendre_shape<'a>( +pub fn legendre_shape + Shape>( cell_type: ReferenceCellType, - points: &impl Array2DAccess<'a, f64>, + points: &T, degree: usize, derivatives: usize, ) -> (usize, usize, usize) { @@ -382,9 +383,9 @@ pub fn legendre_shape<'a>( } /// Tabulate orthonormal polynomials -pub fn tabulate_legendre_polynomials<'a>( +pub fn tabulate_legendre_polynomials + Shape>( cell_type: ReferenceCellType, - points: &impl Array2DAccess<'a, f64>, + points: &T, degree: usize, derivatives: usize, data: &mut impl Array3DAccess, @@ -410,14 +411,13 @@ mod test { use crate::polynomials::*; use approx::*; use bempp_quadrature::simplex_rules::simplex_rule; - use bempp_tools::arrays::{Array2D, Array3D}; - + use bempp_tools::arrays::{to_matrix, Array3D}; #[test] fn test_legendre_interval() { let degree = 6; let rule = simplex_rule(ReferenceCellType::Interval, degree + 1).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 1)); + let points = to_matrix(&rule.points, (rule.npoints, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -448,7 +448,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Triangle, 79).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 2)); + let points = to_matrix(&rule.points, (rule.npoints, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Triangle, @@ -479,7 +479,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Quadrilateral, 85).unwrap(); - let points = Array2D::from_data(rule.points, (rule.npoints, 2)); + let points = to_matrix(&rule.points, (rule.npoints, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, @@ -521,7 +521,7 @@ mod test { p[2 * i] = i as f64 / 10.0; p[2 * i + 1] = p[2 * i] + epsilon; } - let points = Array2D::from_data(p, (20, 1)); + let points = to_matrix(&p, (20, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -560,7 +560,7 @@ mod test { index += 1; } } - let points = Array2D::from_data(p, (165, 2)); + let points = to_matrix(&p, (165, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Triangle, @@ -603,7 +603,7 @@ mod test { p[6 * index + 5] = p[6 * index + 1] + epsilon; } } - let points = Array2D::from_data(p, (300, 2)); + let points = to_matrix(&p, (300, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, @@ -643,7 +643,7 @@ mod test { for (i, pi) in p.iter_mut().enumerate() { *pi = i as f64 / 10.0; } - let points = Array2D::from_data(p, (11, 1)); + let points = to_matrix(&p, (11, 1)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -729,7 +729,7 @@ mod test { p[2 * (11 * i + j) + 1] = j as f64 / 10.0; } } - let points = Array2D::from_data(p, (121, 2)); + let points = to_matrix(&p, (121, 2)); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, diff --git a/grid/Cargo.toml b/grid/Cargo.toml index 38295d16..96d3287b 100644 --- a/grid/Cargo.toml +++ b/grid/Cargo.toml @@ -26,3 +26,5 @@ bempp-element = { path = "../element"} approx = "0.5" itertools = "0.10" mpi = { version = "0.6.*", optional = true } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } diff --git a/grid/examples/curved_cells.rs b/grid/examples/curved_cells.rs index 4ef609e3..2e142761 100644 --- a/grid/examples/curved_cells.rs +++ b/grid/examples/curved_cells.rs @@ -1,6 +1,6 @@ use bempp_grid::grid::SerialGrid; use bempp_grid::io::export_as_gmsh; -use bempp_tools::arrays::{AdjacencyList, Array2D}; +use bempp_tools::arrays::{to_matrix, AdjacencyList}; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::grid::{Geometry, Grid, Topology}; @@ -8,8 +8,8 @@ use bempp_traits::grid::{Geometry, Grid, Topology}; fn main() { // Create a grid of three cells. One cell is a curved triangle, one cell is a flat triangle, the other is a curved quadrilateral let grid = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &vec![ 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 1.5, 0.25, 0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 2.0, 0.5, 0.0, 1.5, 0.75, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0, 2.0, -0.5, 0.0, diff --git a/grid/examples/parallel_grid.rs b/grid/examples/parallel_grid.rs index fbcd81af..89868074 100644 --- a/grid/examples/parallel_grid.rs +++ b/grid/examples/parallel_grid.rs @@ -1,20 +1,21 @@ //? mpirun -n {{NPROCESSES}} --features "mpi" +#[cfg(feature = "mpi")] +use approx::*; #[cfg(feature = "mpi")] use bempp_grid::parallel_grid::ParallelGrid; #[cfg(feature = "mpi")] -use bempp_tools::arrays::{AdjacencyList, Array2D}; +use bempp_tools::arrays::{zero_matrix, AdjacencyList}; #[cfg(feature = "mpi")] -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess}; +use bempp_traits::arrays::AdjacencyListAccess; #[cfg(feature = "mpi")] use bempp_traits::cell::ReferenceCellType; #[cfg(feature = "mpi")] use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; #[cfg(feature = "mpi")] use mpi::{environment::Universe, request::WaitGuard, topology::Communicator, traits::*}; - #[cfg(feature = "mpi")] -use approx::*; +use rlst_dense::RandomAccessMut; #[cfg(feature = "mpi")] fn test_parallel_grid() { @@ -28,7 +29,7 @@ fn test_parallel_grid() { let n = 10; let grid = if rank == 0 { - let mut pts = Array2D::new((n * n, 3)); + let mut pts = zero_matrix((n * n, 3)); let mut i = 0; for y in 0..n { for x in 0..n { diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 19eca8c0..3973bc10 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -1,18 +1,18 @@ //! A serial implementation of a grid use bempp_element::cell; use bempp_element::element::{create_element, CiarletElement}; -use bempp_tools::arrays::{AdjacencyList, Array2D, Array4D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array4DAccess}; +use bempp_tools::arrays::{to_matrix, AdjacencyList, Array4D, Mat}; +use bempp_traits::arrays::{AdjacencyListAccess, Array4DAccess}; 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 rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape}; use std::ptr; - /// Geometry of a serial grid pub struct SerialGeometry { coordinate_elements: Vec, - coordinates: Array2D, + coordinates: Mat, cells: AdjacencyList, element_changes: Vec, index_map: Vec, @@ -35,7 +35,7 @@ fn element_from_npts(cell_type: ReferenceCellType, npts: usize) -> CiarletElemen impl SerialGeometry { pub fn new( - coordinates: Array2D, + coordinates: Mat, cells: &AdjacencyList, cell_types: &[ReferenceCellType], ) -> Self { @@ -94,8 +94,16 @@ impl Geometry for SerialGeometry { self.coordinates.shape().1 } - fn point(&self, i: usize) -> Option<&[f64]> { - self.coordinates.row(i) + fn point(&self, index: usize) -> Option> { + if index > self.point_count() { + None + } else { + let mut pt = vec![0.0; self.dim()]; + for (i, p) in pt.iter_mut().enumerate() { + *p = *self.coordinates.get(index, i).unwrap(); + } + Some(pt) + } } fn point_count(&self) -> usize { @@ -111,11 +119,44 @@ impl Geometry for SerialGeometry { fn index_map(&self) -> &[usize] { &self.index_map } - fn compute_points<'a>( + /* + fn get_compute_points_function + Shape>( + &self, + element: &impl FiniteElement, + points: &T, + ) -> Box { + let npts = points.shape().0; + let mut table = Array4D::::new(element.tabulate_array_shape(0, npts)); + element.tabulate(points, 0, &mut table); + let gdim = self.dim(); + + Box::new(|cell: usize, pts: &mut T| { + for p in 0..npts { + for i in 0..gdim { + *pts.get_mut(p, i).unwrap() = 0.0; + } + } + let vertices = self.cell_vertices(cell).unwrap(); + for (i, n) in vertices.iter().enumerate() { + let pt = self.point(*n).unwrap(); + for p in 0..points.shape().0 { + for (j, pt_j) in pt.iter().enumerate() { + *pts.get_mut(p, j).unwrap() += + *pt_j * *table.get(0, p, i, 0); + } + } + } + }) + } + */ + fn compute_points< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - physical_points: &mut impl Array2DAccess<'a, f64>, + physical_points: &mut TMut, ) { let gdim = self.dim(); if points.shape().0 != physical_points.shape().0 { @@ -129,31 +170,27 @@ impl Geometry for SerialGeometry { element.tabulate(points, 0, &mut data); for p in 0..points.shape().0 { for i in 0..physical_points.shape().1 { - unsafe { - *physical_points.get_unchecked_mut(p, i) = 0.0; - } + *physical_points.get_mut(p, i).unwrap() = 0.0; } } for i in 0..data.shape().2 { - let pt = unsafe { - self.coordinates - .row_unchecked(*self.cells.get_unchecked(cell, i)) - }; + let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); for p in 0..points.shape().0 { for (j, pt_j) in pt.iter().enumerate() { - unsafe { - *physical_points.get_unchecked_mut(p, j) += - *pt_j * data.get_unchecked(0, p, i, 0); - } + *physical_points.get_mut(p, j).unwrap() += + *pt_j * data.get(0, p, i, 0).unwrap(); } } } } - fn compute_normals<'a>( + fn compute_normals< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - normals: &mut impl Array2DAccess<'a, f64>, + normals: &mut TMut, ) { let gdim = self.dim(); if gdim != 3 { @@ -167,57 +204,44 @@ impl Geometry for SerialGeometry { } let element = self.element(cell); let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? - let mut axes = Array2D::::new((2, 3)); + let mut axes = to_matrix(&[0.0; 6], (2, 3)); element.tabulate(points, 1, &mut data); for p in 0..points.shape().0 { for i in 0..axes.shape().0 { for j in 0..axes.shape().1 { - unsafe { - *axes.get_unchecked_mut(i, j) = 0.0; - } + *axes.get_mut(i, j).unwrap() = 0.0; } } for i in 0..data.shape().2 { - let pt = unsafe { - self.coordinates - .row_unchecked(*self.cells.get_unchecked(cell, i)) - }; + let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); for (j, pt_j) in pt.iter().enumerate() { - unsafe { - *axes.get_unchecked_mut(0, j) += *pt_j * data.get_unchecked(1, p, i, 0); - *axes.get_unchecked_mut(1, j) += *pt_j * data.get_unchecked(2, p, i, 0); - } + *axes.get_mut(0, j).unwrap() += *pt_j * data.get(1, p, i, 0).unwrap(); + *axes.get_mut(1, j).unwrap() += *pt_j * data.get(2, p, i, 0).unwrap(); } } - unsafe { - *normals.get_unchecked_mut(p, 0) = *axes.get_unchecked(0, 1) - * *axes.get_unchecked(1, 2) - - *axes.get_unchecked(0, 2) * *axes.get_unchecked(1, 1); - *normals.get_unchecked_mut(p, 1) = *axes.get_unchecked(0, 2) - * *axes.get_unchecked(1, 0) - - *axes.get_unchecked(0, 0) * *axes.get_unchecked(1, 2); - *normals.get_unchecked_mut(p, 2) = *axes.get_unchecked(0, 0) - * *axes.get_unchecked(1, 1) - - *axes.get_unchecked(0, 1) * *axes.get_unchecked(1, 0); - } - let size = unsafe { - (*normals.get_unchecked(p, 0) * *normals.get_unchecked(p, 0) - + *normals.get_unchecked(p, 1) * *normals.get_unchecked(p, 1) - + *normals.get_unchecked(p, 2) * *normals.get_unchecked(p, 2)) - .sqrt() - }; - unsafe { - *normals.get_unchecked_mut(p, 0) /= size; - *normals.get_unchecked_mut(p, 1) /= size; - *normals.get_unchecked_mut(p, 2) /= size; - } + *normals.get_mut(p, 0).unwrap() = *axes.get(0, 1).unwrap() * *axes.get(1, 2).unwrap() + - *axes.get(0, 2).unwrap() * *axes.get(1, 1).unwrap(); + *normals.get_mut(p, 1).unwrap() = *axes.get(0, 2).unwrap() * *axes.get(1, 0).unwrap() + - *axes.get(0, 0).unwrap() * *axes.get(1, 2).unwrap(); + *normals.get_mut(p, 2).unwrap() = *axes.get(0, 0).unwrap() * *axes.get(1, 1).unwrap() + - *axes.get(0, 1).unwrap() * *axes.get(1, 0).unwrap(); + let size = (*normals.get(p, 0).unwrap() * *normals.get(p, 0).unwrap() + + *normals.get(p, 1).unwrap() * *normals.get(p, 1).unwrap() + + *normals.get(p, 2).unwrap() * *normals.get(p, 2).unwrap()) + .sqrt(); + *normals.get_mut(p, 0).unwrap() /= size; + *normals.get_mut(p, 1).unwrap() /= size; + *normals.get_mut(p, 2).unwrap() /= size; } } - fn compute_jacobians<'a>( + fn compute_jacobians< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobians: &mut impl Array2DAccess<'a, f64>, + jacobians: &mut TMut, ) { let gdim = self.dim(); let tdim = points.shape().1; @@ -233,31 +257,24 @@ impl Geometry for SerialGeometry { element.tabulate(points, 1, &mut data); for p in 0..points.shape().0 { for i in 0..jacobians.shape().1 { - unsafe { - *jacobians.get_unchecked_mut(p, i) = 0.0; - } + *jacobians.get_mut(p, i).unwrap() = 0.0; } } for i in 0..data.shape().2 { - let pt = unsafe { - self.coordinates - .row_unchecked(*self.cells.get_unchecked(cell, i)) - }; + let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); for p in 0..points.shape().0 { for (j, pt_j) in pt.iter().enumerate() { for k in 0..tdim { - unsafe { - *jacobians.get_unchecked_mut(p, k + tdim * j) += - *pt_j * data.get_unchecked(k + 1, p, i, 0); - } + *jacobians.get_mut(p, k + tdim * j).unwrap() += + *pt_j * data.get(k + 1, p, i, 0).unwrap(); } } } } } - fn compute_jacobian_determinants<'a>( + fn compute_jacobian_determinants + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, jacobian_determinants: &mut [f64], ) { @@ -266,62 +283,60 @@ impl Geometry for SerialGeometry { if points.shape().0 != jacobian_determinants.len() { panic!("jacobian_determinants has wrong length."); } - let mut js = Array2D::::new((points.shape().0, gdim * tdim)); // TODO: Memory is assigned here. Can we avoid this? + let mut js = to_matrix( + &vec![0.0; points.shape().0 * gdim * tdim], + (points.shape().0, gdim * tdim), + ); // TODO: Memory is assigned here. Can we avoid this? self.compute_jacobians(points, cell, &mut js); // TODO: is it faster if we move this for inside the match statement? for (p, jdet) in jacobian_determinants.iter_mut().enumerate() { *jdet = match tdim { 1 => match gdim { - 1 => unsafe { *js.get_unchecked(p, 0) }, - 2 => unsafe { - ((*js.get_unchecked(p, 0)).powi(2) + (*js.get_unchecked(p, 1)).powi(2)) - .sqrt() - }, - 3 => unsafe { - ((*js.get_unchecked(p, 0)).powi(2) - + (*js.get_unchecked(p, 1)).powi(2) - + (*js.get_unchecked(p, 2)).powi(2)) - .sqrt() - }, + 1 => *js.get(p, 0).unwrap(), + 2 => { + ((*js.get(p, 0).unwrap()).powi(2) + (*js.get(p, 1).unwrap()).powi(2)).sqrt() + } + 3 => ((*js.get(p, 0).unwrap()).powi(2) + + (*js.get(p, 1).unwrap()).powi(2) + + (*js.get(p, 2).unwrap()).powi(2)) + .sqrt(), _ => { panic!("Unsupported dimensions."); } }, 2 => match gdim { - 2 => unsafe { - *js.get_unchecked(p, 0) * *js.get_unchecked(p, 3) - - *js.get_unchecked(p, 1) * *js.get_unchecked(p, 2) - }, - 3 => unsafe { - (((*js.get_unchecked(p, 0)).powi(2) - + (*js.get_unchecked(p, 2)).powi(2) - + (*js.get_unchecked(p, 4)).powi(2)) - * ((*js.get_unchecked(p, 1)).powi(2) - + (*js.get_unchecked(p, 3)).powi(2) - + (*js.get_unchecked(p, 5)).powi(2)) - - (*js.get_unchecked(p, 0) * *js.get_unchecked(p, 1) - + *js.get_unchecked(p, 2) * *js.get_unchecked(p, 3) - + *js.get_unchecked(p, 4) * *js.get_unchecked(p, 5)) - .powi(2)) - .sqrt() - }, + 2 => { + *js.get(p, 0).unwrap() * *js.get(p, 3).unwrap() + - *js.get(p, 1).unwrap() * *js.get(p, 2).unwrap() + } + 3 => (((*js.get(p, 0).unwrap()).powi(2) + + (*js.get(p, 2).unwrap()).powi(2) + + (*js.get(p, 4).unwrap()).powi(2)) + * ((*js.get(p, 1).unwrap()).powi(2) + + (*js.get(p, 3).unwrap()).powi(2) + + (*js.get(p, 5).unwrap()).powi(2)) + - (*js.get(p, 0).unwrap() * *js.get(p, 1).unwrap() + + *js.get(p, 2).unwrap() * *js.get(p, 3).unwrap() + + *js.get(p, 4).unwrap() * *js.get(p, 5).unwrap()) + .powi(2)) + .sqrt(), _ => { panic!("Unsupported dimensions."); } }, 3 => match gdim { - 3 => unsafe { - *js.get_unchecked(p, 0) - * (*js.get_unchecked(p, 4) * *js.get_unchecked(p, 8) - - *js.get_unchecked(p, 5) * *js.get_unchecked(p, 7)) - - *js.get_unchecked(p, 1) - * (*js.get_unchecked(p, 3) * *js.get_unchecked(p, 8) - - *js.get_unchecked(p, 5) * *js.get_unchecked(p, 6)) - + *js.get_unchecked(p, 2) - * (*js.get_unchecked(p, 3) * *js.get_unchecked(p, 7) - - *js.get_unchecked(p, 4) * *js.get_unchecked(p, 6)) - }, + 3 => { + *js.get(p, 0).unwrap() + * (*js.get(p, 4).unwrap() * *js.get(p, 8).unwrap() + - *js.get(p, 5).unwrap() * *js.get(p, 7).unwrap()) + - *js.get(p, 1).unwrap() + * (*js.get(p, 3).unwrap() * *js.get(p, 8).unwrap() + - *js.get(p, 5).unwrap() * *js.get(p, 6).unwrap()) + + *js.get(p, 2).unwrap() + * (*js.get(p, 3).unwrap() * *js.get(p, 7).unwrap() + - *js.get(p, 4).unwrap() * *js.get(p, 6).unwrap()) + } _ => { panic!("Unsupported dimensions."); } @@ -332,11 +347,14 @@ impl Geometry for SerialGeometry { } } } - fn compute_jacobian_inverses<'a>( + fn compute_jacobian_inverses< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobian_inverses: &mut impl Array2DAccess<'a, f64>, + jacobian_inverses: &mut TMut, ) { let gdim = self.dim(); let tdim = points.shape().1; @@ -352,17 +370,17 @@ impl Geometry for SerialGeometry { && element.degree() == 1 { // Map is affine - let mut js = Array2D::::new((points.shape().0, gdim * tdim)); // TODO: Memory is assigned here. Can we avoid this? + let mut js = to_matrix( + &vec![0.0; points.shape().0 * gdim * tdim], + (points.shape().0, gdim * tdim), + ); // TODO: Memory is assigned here. Can we avoid this? self.compute_jacobians(points, cell, &mut js); // TODO: is it faster if we move this for inside the if statement? for p in 0..points.shape().0 { if tdim == 1 { if gdim == 1 { - unsafe { - *jacobian_inverses.get_unchecked_mut(p, 0) = - 1.0 / *js.get_unchecked(p, 0); - } + *jacobian_inverses.get_mut(p, 0).unwrap() = 1.0 / *js.get(p, 0).unwrap(); } else if gdim == 2 { unimplemented!("Inverse jacobian for this dimension not implemented yet."); } else if gdim == 3 { @@ -372,76 +390,64 @@ impl Geometry for SerialGeometry { } } else if tdim == 2 { if gdim == 2 { - let det = unsafe { - *js.get_unchecked(p, 0) * *js.get_unchecked(p, 3) - - *js.get_unchecked(p, 1) * *js.get_unchecked(p, 2) - }; - unsafe { - *jacobian_inverses.get_unchecked_mut(p, 0) = - js.get_unchecked(p, 3) / det; - *jacobian_inverses.get_unchecked_mut(p, 1) = - -js.get_unchecked(p, 1) / det; - *jacobian_inverses.get_unchecked_mut(p, 2) = - -js.get_unchecked(p, 2) / det; - *jacobian_inverses.get_unchecked_mut(p, 3) = - js.get_unchecked(p, 0) / det; - } + let det = *js.get(p, 0).unwrap() * *js.get(p, 3).unwrap() + - *js.get(p, 1).unwrap() * *js.get(p, 2).unwrap(); + *jacobian_inverses.get_mut(p, 0).unwrap() = js.get(p, 3).unwrap() / det; + *jacobian_inverses.get_mut(p, 1).unwrap() = -js.get(p, 1).unwrap() / det; + *jacobian_inverses.get_mut(p, 2).unwrap() = -js.get(p, 2).unwrap() / det; + *jacobian_inverses.get_mut(p, 3).unwrap() = js.get(p, 0).unwrap() / det; } else if gdim == 3 { - let c = unsafe { - (*js.get_unchecked(p, 3) * *js.get_unchecked(p, 4) - - *js.get_unchecked(p, 2) * *js.get_unchecked(p, 5)) + let c = (*js.get(p, 3).unwrap() * *js.get(p, 4).unwrap() + - *js.get(p, 2).unwrap() * *js.get(p, 5).unwrap()) + .powi(2) + + (*js.get(p, 5).unwrap() * *js.get(p, 0).unwrap() + - *js.get(p, 4).unwrap() * *js.get(p, 1).unwrap()) .powi(2) - + (*js.get_unchecked(p, 5) * *js.get_unchecked(p, 0) - - *js.get_unchecked(p, 4) * *js.get_unchecked(p, 1)) - .powi(2) - + (*js.get_unchecked(p, 1) * *js.get_unchecked(p, 2) - - *js.get_unchecked(p, 0) * *js.get_unchecked(p, 3)) - .powi(2) - }; - unsafe { - *jacobian_inverses.get_unchecked_mut(p, 0) = (*js.get_unchecked(p, 0) - * ((*js.get_unchecked(p, 5)).powi(2) - + (*js.get_unchecked(p, 3)).powi(2)) - - *js.get_unchecked(p, 1) - * (*js.get_unchecked(p, 2) * *js.get_unchecked(p, 3) - + *js.get_unchecked(p, 4) * *js.get_unchecked(p, 5))) - / c; - *jacobian_inverses.get_unchecked_mut(p, 1) = (*js.get_unchecked(p, 2) - * ((*js.get_unchecked(p, 1)).powi(2) - + (*js.get_unchecked(p, 5)).powi(2)) - - *js.get_unchecked(p, 3) - * (*js.get_unchecked(p, 4) * *js.get_unchecked(p, 5) - + *js.get_unchecked(p, 0) * *js.get_unchecked(p, 1))) - / c; - *jacobian_inverses.get_unchecked_mut(p, 2) = (*js.get_unchecked(p, 4) - * ((*js.get_unchecked(p, 3)).powi(2) - + (*js.get_unchecked(p, 1)).powi(2)) - - *js.get_unchecked(p, 5) - * (*js.get_unchecked(p, 0) * *js.get_unchecked(p, 1) - + *js.get_unchecked(p, 2) * *js.get_unchecked(p, 3))) - / c; - *jacobian_inverses.get_unchecked_mut(p, 3) = (*js.get_unchecked(p, 1) - * ((*js.get_unchecked(p, 4)).powi(2) - + (*js.get_unchecked(p, 2)).powi(2)) - - *js.get_unchecked(p, 0) - * (*js.get_unchecked(p, 2) * *js.get_unchecked(p, 3) - + *js.get_unchecked(p, 4) * *js.get_unchecked(p, 5))) - / c; - *jacobian_inverses.get_unchecked_mut(p, 4) = (*js.get_unchecked(p, 3) - * ((*js.get_unchecked(p, 0)).powi(2) - + (*js.get_unchecked(p, 4)).powi(2)) - - *js.get_unchecked(p, 2) - * (*js.get_unchecked(p, 4) * *js.get_unchecked(p, 5) - + *js.get_unchecked(p, 0) * *js.get_unchecked(p, 1))) - / c; - *jacobian_inverses.get_unchecked_mut(p, 5) = (*js.get_unchecked(p, 5) - * ((*js.get_unchecked(p, 2)).powi(2) - + (*js.get_unchecked(p, 0)).powi(2)) - - *js.get_unchecked(p, 4) - * (*js.get_unchecked(p, 0) * *js.get_unchecked(p, 1) - + *js.get_unchecked(p, 2) * *js.get_unchecked(p, 3))) - / c; - } + + (*js.get(p, 1).unwrap() * *js.get(p, 2).unwrap() + - *js.get(p, 0).unwrap() * *js.get(p, 3).unwrap()) + .powi(2); + *jacobian_inverses.get_mut(p, 0).unwrap() = (*js.get(p, 0).unwrap() + * ((*js.get(p, 5).unwrap()).powi(2) + + (*js.get(p, 3).unwrap()).powi(2)) + - *js.get(p, 1).unwrap() + * (*js.get(p, 2).unwrap() * *js.get(p, 3).unwrap() + + *js.get(p, 4).unwrap() * *js.get(p, 5).unwrap())) + / c; + *jacobian_inverses.get_mut(p, 1).unwrap() = (*js.get(p, 2).unwrap() + * ((*js.get(p, 1).unwrap()).powi(2) + + (*js.get(p, 5).unwrap()).powi(2)) + - *js.get(p, 3).unwrap() + * (*js.get(p, 4).unwrap() * *js.get(p, 5).unwrap() + + *js.get(p, 0).unwrap() * *js.get(p, 1).unwrap())) + / c; + *jacobian_inverses.get_mut(p, 2).unwrap() = (*js.get(p, 4).unwrap() + * ((*js.get(p, 3).unwrap()).powi(2) + + (*js.get(p, 1).unwrap()).powi(2)) + - *js.get(p, 5).unwrap() + * (*js.get(p, 0).unwrap() * *js.get(p, 1).unwrap() + + *js.get(p, 2).unwrap() * *js.get(p, 3).unwrap())) + / c; + *jacobian_inverses.get_mut(p, 3).unwrap() = (*js.get(p, 1).unwrap() + * ((*js.get(p, 4).unwrap()).powi(2) + + (*js.get(p, 2).unwrap()).powi(2)) + - *js.get(p, 0).unwrap() + * (*js.get(p, 2).unwrap() * *js.get(p, 3).unwrap() + + *js.get(p, 4).unwrap() * *js.get(p, 5).unwrap())) + / c; + *jacobian_inverses.get_mut(p, 4).unwrap() = (*js.get(p, 3).unwrap() + * ((*js.get(p, 0).unwrap()).powi(2) + + (*js.get(p, 4).unwrap()).powi(2)) + - *js.get(p, 2).unwrap() + * (*js.get(p, 4).unwrap() * *js.get(p, 5).unwrap() + + *js.get(p, 0).unwrap() * *js.get(p, 1).unwrap())) + / c; + *jacobian_inverses.get_mut(p, 5).unwrap() = (*js.get(p, 5).unwrap() + * ((*js.get(p, 2).unwrap()).powi(2) + + (*js.get(p, 0).unwrap()).powi(2)) + - *js.get(p, 4).unwrap() + * (*js.get(p, 0).unwrap() * *js.get(p, 1).unwrap() + + *js.get(p, 2).unwrap() * *js.get(p, 3).unwrap())) + / c; } else { panic!("Unsupported dimensions."); } @@ -541,7 +547,7 @@ impl SerialTopology { starts[i + 1] }; for c in cstart..cend { - let cell = unsafe { cells.row_unchecked(c) }; + let cell = cells.row(c).unwrap(); for e in &ref_entities { let vertices = e.iter().map(|x| cell[*x]).collect::>(); let mut found = false; @@ -646,8 +652,7 @@ impl SerialTopology { starts[i + 1] }; for c in cstart..cend { - for (e, t) in izip!(unsafe { cell_to_entities0.row_unchecked(c) }, &etypes) - { + for (e, t) in izip!(cell_to_entities0.row(c).unwrap(), &etypes) { sub_cell_types[*e] = *t; } } @@ -731,7 +736,7 @@ impl Topology<'_> for SerialTopology { } fn cell(&self, index: usize) -> Option<&[usize]> { if index < self.entity_count(self.dim) { - Some(unsafe { self.connectivity(self.dim, 0).row_unchecked(index) }) + Some(self.connectivity(self.dim, 0).row(index).unwrap()) } else { None } @@ -770,7 +775,7 @@ pub struct SerialGrid { impl SerialGrid { pub fn new( - coordinates: Array2D, + coordinates: Mat, cells: AdjacencyList, cell_types: Vec, ) -> Self { @@ -814,7 +819,7 @@ mod test { #[test] fn test_connectivity() { let g = SerialGrid::new( - Array2D::from_data(vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0], (4, 2)), + to_matrix(&[0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0], (4, 2)), AdjacencyList::from_data(vec![0, 1, 2, 2, 1, 3], vec![0, 3, 6]), vec![ReferenceCellType::Triangle; 2], ); @@ -915,8 +920,8 @@ mod test { #[test] fn test_serial_triangle_grid_octahedron() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, ], @@ -941,8 +946,8 @@ mod test { #[test] fn test_serial_triangle_grid_screen() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0, 1.0, ], @@ -967,8 +972,8 @@ mod test { #[test] fn test_serial_mixed_grid_screen() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0, 1.0, ], @@ -999,8 +1004,8 @@ mod test { fn test_higher_order_grid() { let s = 1.0 / (2.0_f64).sqrt(); let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 1.0, 0.0, s, s, 0.0, 1.0, -s, s, -1.0, 0.0, -s, -s, 0.0, -1.0, s, -s, ], (9, 2), @@ -1019,8 +1024,8 @@ mod test { #[test] fn test_higher_order_mixed_grid() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 1.5, 0.25, 0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 2.0, 0.5, 0.0, 1.5, 0.75, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0, 2.0, -0.5, 0.0, @@ -1048,8 +1053,8 @@ mod test { #[test] fn test_points_and_jacobians() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 2.0, 2.0, 0.0, 4.0, 2.0, 0.0, 5.0, 3.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 1.0, 0.0, -1.0, 1.0, ], @@ -1061,13 +1066,14 @@ mod test { assert_eq!(g.topology().dim(), 2); assert_eq!(g.geometry().dim(), 3); - let points = Array2D::from_data( - vec![0.2, 0.0, 0.5, 0.5, 1.0 / 3.0, 1.0 / 3.0, 0.15, 0.3], + let points = to_matrix( + &[0.2, 0.0, 0.5, 0.5, 1.0 / 3.0, 1.0 / 3.0, 0.15, 0.3], (4, 2), ); // Test compute_points - let mut physical_points = Array2D::new((points.shape().0, 3)); + let mut physical_points = + to_matrix(&vec![0.0; points.shape().0 * 3], (points.shape().0, 3)); g.geometry() .compute_points(&points, 0, &mut physical_points); assert_relative_eq!( @@ -1194,7 +1200,7 @@ mod test { ); // Test compute_jacobians - let mut jacobians = Array2D::new((points.shape().0, 6)); + let mut jacobians = to_matrix(&vec![0.0; points.shape().0 * 6], (points.shape().0, 6)); g.geometry().compute_jacobians(&points, 0, &mut jacobians); for i in 0..3 { assert_relative_eq!(*jacobians.get(i, 0).unwrap(), 2.0, max_relative = 1e-14); @@ -1228,7 +1234,7 @@ mod test { } // Test compute_jacobian_inverses - let mut jinvs = Array2D::new((points.shape().0, 6)); + let mut jinvs = to_matrix(&vec![0.0; points.shape().0 * 6], (points.shape().0, 6)); g.geometry() .compute_jacobian_inverses(&points, 0, &mut jinvs); for i in 0..3 { @@ -1254,8 +1260,8 @@ mod test { #[test] fn test_normals() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 1.0, ], (5, 3), @@ -1270,9 +1276,9 @@ mod test { ], ); - let pt = Array2D::from_data(vec![1.0 / 3.0, 1.0 / 3.0], (1, 2)); + let pt = to_matrix(&[1.0 / 3.0, 1.0 / 3.0], (1, 2)); - let mut normal = Array2D::::new((1, 3)); + let mut normal = to_matrix(&[0.0; 3], (1, 3)); g.geometry().compute_normals(&pt, 0, &mut normal); assert_relative_eq!(*normal.get(0, 0).unwrap(), 0.0); @@ -1292,8 +1298,8 @@ mod test { // Test a curved quadrilateral cell let curved_g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, -1.0, 0.0, -1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, ], @@ -1303,11 +1309,8 @@ mod test { vec![ReferenceCellType::Quadrilateral], ); - let points = Array2D::from_data( - vec![0.0, 0.0, 0.2, 0.3, 0.5, 0.9, 0.7, 1.0, 1.0, 0.3], - (5, 2), - ); - let mut normals = Array2D::::new((5, 3)); + let points = to_matrix(&[0.0, 0.0, 0.2, 0.3, 0.5, 0.9, 0.7, 1.0, 1.0, 0.3], (5, 2)); + let mut normals = to_matrix(&[0.0; 15], (5, 3)); curved_g .geometry() diff --git a/grid/src/io.rs b/grid/src/io.rs index 4d52ca44..b9eb0fb4 100644 --- a/grid/src/io.rs +++ b/grid/src/io.rs @@ -123,7 +123,7 @@ mod test { use crate::grid::SerialGrid; use crate::io::*; use crate::shapes::regular_sphere; - use bempp_tools::arrays::{AdjacencyList, Array2D}; + use bempp_tools::arrays::{to_matrix, AdjacencyList}; use bempp_traits::cell::ReferenceCellType; #[test] @@ -135,8 +135,8 @@ mod test { #[test] fn test_gmsh_output_quads() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0, 1.0, ], @@ -154,8 +154,8 @@ mod test { #[test] fn test_gmsh_output_mixed_cell_type() { let g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 0.5, 0.0, 1.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.5, 1.0, 1.0, 1.0, ], diff --git a/grid/src/parallel_grid.rs b/grid/src/parallel_grid.rs index fd7f318f..4edc62da 100644 --- a/grid/src/parallel_grid.rs +++ b/grid/src/parallel_grid.rs @@ -1,11 +1,12 @@ //! A parallel implementation of a grid use crate::grid::{SerialGeometry, SerialTopology}; use bempp_element::element::CiarletElement; -use bempp_tools::arrays::{AdjacencyList, Array2D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess}; +use bempp_tools::arrays::{zero_matrix, AdjacencyList, Mat}; +use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; use mpi::{request::WaitGuard, topology::Communicator, traits::*}; +use rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape}; /// Geometry of a parallel grid pub struct ParallelGeometry<'a, C: Communicator> { @@ -16,7 +17,7 @@ pub struct ParallelGeometry<'a, C: Communicator> { impl<'a, C: Communicator> ParallelGeometry<'a, C> { pub fn new( comm: &'a C, - coordinates: Array2D, + coordinates: Mat, cells: &AdjacencyList, cell_types: &Vec, ) -> Self { @@ -44,7 +45,7 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { self.serial_geometry.dim() } - fn point(&self, i: usize) -> Option<&[f64]> { + fn point(&self, i: usize) -> Option> { self.serial_geometry.point(i) } @@ -61,46 +62,58 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { fn index_map(&self) -> &[usize] { self.serial_geometry.index_map() } - fn compute_points<'b>( + fn compute_points< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'b, f64>, + points: &T, cell: usize, - physical_points: &mut impl Array2DAccess<'b, f64>, + physical_points: &mut TMut, ) { self.serial_geometry .compute_points(points, cell, physical_points) } - fn compute_normals<'b>( + fn compute_normals< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'b, f64>, + points: &T, cell: usize, - normals: &mut impl Array2DAccess<'b, f64>, + normals: &mut TMut, ) { self.serial_geometry.compute_points(points, cell, normals) } - fn compute_jacobians<'b>( + fn compute_jacobians< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'b, f64>, + points: &T, cell: usize, - jacobians: &mut impl Array2DAccess<'b, f64>, + jacobians: &mut TMut, ) { self.serial_geometry .compute_jacobians(points, cell, jacobians) } - fn compute_jacobian_determinants<'b>( + fn compute_jacobian_determinants + Shape>( &self, - points: &impl Array2DAccess<'b, f64>, + points: &T, cell: usize, jacobian_determinants: &mut [f64], ) { self.serial_geometry .compute_jacobian_determinants(points, cell, jacobian_determinants) } - fn compute_jacobian_inverses<'b>( + fn compute_jacobian_inverses< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'b, f64>, + points: &T, cell: usize, - jacobian_inverses: &mut impl Array2DAccess<'b, f64>, + jacobian_inverses: &mut TMut, ) { self.serial_geometry .compute_jacobian_inverses(points, cell, jacobian_inverses) @@ -179,7 +192,7 @@ pub struct ParallelGrid<'a, C: Communicator> { impl<'a, C: Communicator> ParallelGrid<'a, C> { pub fn new( comm: &'a C, - coordinates: Array2D, + coordinates: Mat, cells: AdjacencyList, cell_types: Vec, cell_owners: Vec, @@ -214,8 +227,8 @@ impl<'a, C: Communicator> ParallelGrid<'a, C> { vertex_indices_per_proc[p].push(*v); vertex_owners_per_proc[p].push(vertex_owners[*v].0 as usize); vertex_local_indices_per_proc[p].push(vertex_owners[*v].1); - for x in unsafe { coordinates.row_unchecked(*v) } { - coordinates_per_proc[p].push(*x) + for i in 0..coordinates.shape().1 { + coordinates_per_proc[p].push(*coordinates.get(*v, i).unwrap()) } } } @@ -237,8 +250,8 @@ impl<'a, C: Communicator> ParallelGrid<'a, C> { vertex_indices_per_proc[p].push(*v); vertex_owners_per_proc[p].push(vertex_owners[*v].0 as usize); vertex_local_indices_per_proc[p].push(vertex_owners[*v].1); - for x in unsafe { coordinates.row_unchecked(*v) } { - coordinates_per_proc[p].push(*x) + for i in 0..coordinates.shape().1 { + coordinates_per_proc[p].push(*coordinates.get(*v, i).unwrap()) } } cells_per_proc[p].push( @@ -379,7 +392,7 @@ impl<'a, C: Communicator> ParallelGrid<'a, C> { panic!("Unsupported cell type"); }); } - let mut coordinates = Array2D::::new((vertex_indices.len(), gdim)); + let mut coordinates = zero_matrix((vertex_indices.len(), gdim)); for i in 0..vertex_indices.len() { for j in 0..gdim { *coordinates.get_mut(i, j).unwrap() = flat_coordinates[i * gdim + j]; diff --git a/grid/src/shapes.rs b/grid/src/shapes.rs index 1311b3bd..d9822177 100644 --- a/grid/src/shapes.rs +++ b/grid/src/shapes.rs @@ -2,10 +2,11 @@ use crate::grid::SerialGrid; use bempp_element::cell::Triangle; -use bempp_tools::arrays::{AdjacencyList, Array2D}; -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess}; +use bempp_tools::arrays::{to_matrix, AdjacencyList}; +use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::grid::{Geometry, Grid, Topology}; +use rlst_dense::{RandomAccessByRef, RandomAccessMut}; /// Create a regular sphere /// @@ -14,8 +15,8 @@ use bempp_traits::grid::{Geometry, Grid, Topology}; /// each edge). The new points are then scaled so that they are a distance of 1 from the origin. pub fn regular_sphere(refinement_level: usize) -> SerialGrid { let mut g = SerialGrid::new( - Array2D::from_data( - vec![ + to_matrix( + &[ 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, ], @@ -34,7 +35,7 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { let nvertices_old = g.topology().entity_count(0); let ncells_old = g.topology().entity_count(2); let nvertices = g.topology().entity_count(0) + g.topology().entity_count(1); - let mut coordinates = Array2D::::new((nvertices, 3)); + let mut coordinates = to_matrix(&vec![0.0; nvertices * 3], (nvertices, 3)); let mut cells = AdjacencyList::::new(); for i in 0..ncells_old { @@ -108,10 +109,10 @@ mod test { fn test_normal_is_outward() { for i in 0..3 { let g = regular_sphere(i); - let points = Array2D::from_data(vec![1.0 / 3.0, 1.0 / 3.0], (1, 2)); + let points = to_matrix(&[1.0 / 3.0, 1.0 / 3.0], (1, 2)); - let mut mapped_pt = Array2D::::new((1, 3)); - let mut normal = Array2D::::new((1, 3)); + let mut mapped_pt = to_matrix(&[0.0; 3], (1, 3)); + let mut normal = to_matrix(&[0.0; 3], (1, 3)); for i in 0..g.geometry().cell_count() { g.geometry().compute_points(&points, i, &mut mapped_pt); diff --git a/kernel/src/traits.rs b/kernel/src/traits.rs index eb722366..78fd74cb 100644 --- a/kernel/src/traits.rs +++ b/kernel/src/traits.rs @@ -4,7 +4,7 @@ use crate::types::KernelType; use bempp_traits::types::Scalar; /// Interface to evaluating Green's functions for given sources and targets. -pub trait Kernel { +pub trait Kernel: Sync { type T: Scalar; /// Single threaded evaluation of Green's functions. diff --git a/python/test/test_dependencies.py b/python/test/test_dependencies.py index 170b7c12..1d27b43a 100644 --- a/python/test/test_dependencies.py +++ b/python/test/test_dependencies.py @@ -25,6 +25,7 @@ def test_dependencies(): cargos = walk_dirs(root_dir, "Cargo.toml") deps = {} + errors = [] for c in cargos: with open(c, "rb") as f: data = tomllib.load(f) @@ -38,6 +39,8 @@ def test_dependencies(): if d not in deps: deps[d] = (info, c) elif deps[d][0] != info: - raise ValueError( + errors.append( f"Version of {d} in {c} ({info}) does not agree " f"with version in {deps[d][1]} ({deps[d][0]})") + if len(errors) > 0: + raise ValueError("\n".join(errors)) diff --git a/tools/Cargo.toml b/tools/Cargo.toml index 0f877963..ebc5b177 100644 --- a/tools/Cargo.toml +++ b/tools/Cargo.toml @@ -24,3 +24,4 @@ libc = "0.2" num = "0.4" bempp-traits = { path = "../traits"} rayon = "1.7" +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } diff --git a/tools/examples/array2d.rs b/tools/examples/array2d.rs deleted file mode 100644 index d594c9e0..00000000 --- a/tools/examples/array2d.rs +++ /dev/null @@ -1,52 +0,0 @@ -use bempp_tools::arrays::Array2D; -use bempp_traits::arrays::Array2DAccess; - -fn main() { - // Create an 2D array of zero integers of a given size - // A 2D array is a two-dimensional storage format that allows items to be accessed using two indices (a row and a column) - let mut a = Array2D::::new((5, 3)); - - // Set values in the first row of the array - *a.get_mut(0, 0).unwrap() = 3; - *a.get_mut(0, 1).unwrap() = 1; - *a.get_mut(0, 2).unwrap() = 4; - - // Print values from the first column of the array - println!( - "The first column of a is {} {} {} {} {}", - *a.get(0, 0).unwrap(), - *a.get(1, 0).unwrap(), - *a.get(2, 0).unwrap(), - *a.get(3, 0).unwrap(), - *a.get(4, 0).unwrap() - ); - - // Items can also be accessed without bound checking using the unchecked methods - println!( - "The entry in the first row and first column of a is {}", - unsafe { *a.get_unchecked(0, 0) } - ); - - // Print values from the first row of the array - // This can be done in two ways - println!( - "The first row of a is {} {} {}", - *a.get(0, 0).unwrap(), - *a.get(0, 1).unwrap(), - *a.get(0, 2).unwrap() - ); - let row0 = a.row(0).unwrap(); - println!("The first row of a is {} {} {}", row0[0], row0[1], row0[2]); - - // Print the shape of the array - println!("a has {} rows and {} columns", a.shape().0, a.shape().1); - - // Iterate through the rows of a - for (i, row) in a.iter_rows().enumerate() { - println!( - "The sum of the values in row {} of a is {}", - i, - row.iter().sum::() - ); - } -} diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index 008e12ef..dfd36c83 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -1,95 +1,27 @@ //! Containers to store multi-dimensional data -use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess}; +use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess}; use num::Num; +use rlst_dense::{operations::transpose::Scalar, rlst_dynamic_mat, UnsafeRandomAccessMut}; use std::clone::Clone; -/// A two-dimensional rectangular array -#[derive(Clone)] -pub struct Array2D { - /// The data in the array, in row-major order - pub data: Vec, - /// The shape of the array - shape: (usize, usize), -} - -impl Array2D { - /// Create an array from a data vector - pub fn new(shape: (usize, usize)) -> Self { - Self { - data: vec![T::zero(); shape.0 * shape.1], - shape, +pub type Mat = rlst_dense::Matrix< + T, + rlst_dense::base_matrix::BaseMatrix, rlst_dense::Dynamic>, + rlst_dense::Dynamic, +>; + +pub fn to_matrix(data: &[T], shape: (usize, usize)) -> Mat { + let mut mat = rlst_dynamic_mat![T, shape]; + for (i, d) in data.iter().enumerate() { + unsafe { + *mat.get_unchecked_mut(i / shape.1, i % shape.1) = *d; } } - - /// Create an array from a data vector - pub fn from_data(data: Vec, shape: (usize, usize)) -> Self { - assert_eq!(data.len(), shape.0 * shape.1); - Self { data, shape } - } + mat } -impl<'a, T: Num + 'a> Array2DAccess<'a, T> for Array2D { - type I = Array2DRowIterator<'a, T>; - - fn get(&self, index0: usize, index1: usize) -> Option<&T> { - if index0 >= self.shape.0 || index1 >= self.shape.1 { - None - } else { - unsafe { Some(self.get_unchecked(index0, index1)) } - } - } - /// Get a mutable item from the array - fn get_mut(&mut self, index0: usize, index1: usize) -> Option<&mut T> { - if index0 >= self.shape.0 || index1 >= self.shape.1 { - None - } else { - unsafe { Some(self.get_unchecked_mut(index0, index1)) } - } - } - /// Get a row of the array - fn row(&self, index: usize) -> Option<&[T]> { - if index >= self.shape.0 { - None - } else { - unsafe { Some(self.row_unchecked(index)) } - } - } - /// Get an item from the array without checking bounds - unsafe fn get_unchecked(&self, index0: usize, index1: usize) -> &T { - self.data.get_unchecked(index0 * self.shape.1 + index1) - } - /// Get a mutable item from the array without checking bounds - unsafe fn get_unchecked_mut(&mut self, index0: usize, index1: usize) -> &mut T { - self.data.get_unchecked_mut(index0 * self.shape.1 + index1) - } - /// Get a row of the array without checking bounds - unsafe fn row_unchecked(&self, index: usize) -> &[T] { - &self.data[index * self.shape.1..(index + 1) * self.shape.1] - } - /// Get the shape of the array - fn shape(&self) -> &(usize, usize) { - &self.shape - } - /// Iterate through the rows - fn iter_rows(&'a self) -> Self::I { - Array2DRowIterator:: { - array: self, - index: 0, - } - } -} - -pub struct Array2DRowIterator<'a, T: Num> { - array: &'a Array2D, - index: usize, -} - -impl<'a, T: Num> Iterator for Array2DRowIterator<'a, T> { - type Item = &'a [T]; - fn next(&mut self) -> Option { - self.index += 1; - self.array.row(self.index - 1) - } +pub fn zero_matrix(shape: (usize, usize)) -> Mat { + rlst_dynamic_mat![T, shape] } /// A three-dimensional rectangular array @@ -338,57 +270,6 @@ impl<'a, T: Num> Iterator for AdjacencyListRowIterator<'a, T> { mod test { use crate::arrays::*; - #[test] - fn test_array_2d() { - let mut arr = Array2D::from_data(vec![1, 2, 3, 4, 5, 6], (2, 3)); - assert_eq!(*arr.get(0, 0).unwrap(), 1); - assert_eq!(*arr.get(0, 1).unwrap(), 2); - assert_eq!(*arr.get(0, 2).unwrap(), 3); - assert_eq!(*arr.get(1, 0).unwrap(), 4); - assert_eq!(*arr.get(1, 1).unwrap(), 5); - assert_eq!(*arr.get(1, 2).unwrap(), 6); - - let row1 = arr.row(1).unwrap(); - assert_eq!(row1.len(), 3); - assert_eq!(row1[0], 4); - assert_eq!(row1[1], 5); - assert_eq!(row1[2], 6); - - *arr.get_mut(1, 2).unwrap() = 7; - assert_eq!(*arr.get(1, 2).unwrap(), 7); - - for (index, row) in arr.iter_rows().enumerate() { - assert_eq!(*arr.get(index, 0).unwrap(), row[0]); - } - - let mut arr2 = Array2D::from_data(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0], (2, 3)); - assert_eq!(*arr2.get(0, 0).unwrap(), 1.0); - assert_eq!(*arr2.get(0, 1).unwrap(), 2.0); - assert_eq!(*arr2.get(0, 2).unwrap(), 3.0); - assert_eq!(*arr2.get(1, 0).unwrap(), 4.0); - assert_eq!(*arr2.get(1, 1).unwrap(), 5.0); - assert_eq!(*arr2.get(1, 2).unwrap(), 6.0); - - let row1 = arr2.row(1).unwrap(); - assert_eq!(row1.len(), 3); - assert_eq!(row1[0], 4.0); - assert_eq!(row1[1], 5.0); - assert_eq!(row1[2], 6.0); - - *arr2.get_mut(1, 2).unwrap() = 7.; - assert_eq!(*arr2.get(1, 2).unwrap(), 7.0); - - let mut arr3 = Array2D::::new((4, 5)); - assert_eq!(*arr3.get(1, 2).unwrap(), 0); - *arr3.get_mut(1, 2).unwrap() = 5; - assert_eq!(*arr3.get(1, 2).unwrap(), 5); - - let mut arr4 = Array2D::::new((4, 5)); - assert_eq!(*arr4.get(1, 2).unwrap(), 0.0); - *arr4.get_mut(1, 2).unwrap() = 5.0; - assert_eq!(*arr4.get(1, 2).unwrap(), 5.0); - } - #[test] fn test_adjacency_list() { let mut arr = AdjacencyList::from_data(vec![1, 2, 3, 4, 5, 6], vec![0, 2, 3, 6]); diff --git a/traits/Cargo.toml b/traits/Cargo.toml index 9b642c10..b6df2c86 100644 --- a/traits/Cargo.toml +++ b/traits/Cargo.toml @@ -22,3 +22,4 @@ crate-type = ["lib", "cdylib"] cauchy="0.4.*" thiserror="1.*" num = "0.4" +rlst-common = { git = "https://github.com/linalg-rs/rlst.git" } diff --git a/traits/src/arrays.rs b/traits/src/arrays.rs index a9625cbb..4af2a01c 100644 --- a/traits/src/arrays.rs +++ b/traits/src/arrays.rs @@ -1,71 +1,6 @@ //! Containers to store multi-dimensional data use num::Num; -pub trait Array1DAccess<'a, T: Num> { - type I: Iterator; - - /// Get an item from the array - fn get(&self, index: usize) -> Option<&T>; - - /// Get a mutable item from the array - fn get_mut(&mut self, index: usize) -> Option<&mut T>; - - /// Get an item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked(&self, index: usize) -> &T; - - /// Get a mutable item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked_mut(&mut self, index: usize) -> &mut T; - - /// Get the shape of the array - fn shape(&self) -> &usize; - - /// Iterate through the values - fn iter(&'a self) -> Self::I; -} - -pub trait Array2DAccess<'a, T: Num> { - type I: Iterator; - - /// Get an item from the array - fn get(&self, index0: usize, index1: usize) -> Option<&T>; - - /// Get a mutable item from the array - fn get_mut(&mut self, index0: usize, index1: usize) -> Option<&mut T>; - - /// Get a row of the array - fn row(&self, index: usize) -> Option<&[T]>; - - /// Get an item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked(&self, index0: usize, index1: usize) -> &T; - - /// Get a mutable item from the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn get_unchecked_mut(&mut self, index0: usize, index1: usize) -> &mut T; - - /// Get a row of the array without checking bounds - /// - /// # Safety - /// This function does not perform bound checks - unsafe fn row_unchecked(&self, index: usize) -> &[T]; - - /// Get the shape of the array - fn shape(&self) -> &(usize, usize); - - /// Iterate through the rows - fn iter_rows(&'a self) -> Self::I; -} - pub trait Array3DAccess { /// Get an item from the array fn get(&self, index0: usize, index1: usize, index2: usize) -> Option<&T>; diff --git a/traits/src/element.rs b/traits/src/element.rs index 24c27fd9..91bd4600 100644 --- a/traits/src/element.rs +++ b/traits/src/element.rs @@ -1,7 +1,8 @@ //! Finite element definitions -use crate::arrays::{Array2DAccess, Array4DAccess}; +use crate::arrays::Array4DAccess; use crate::cell::ReferenceCellType; +use rlst_common::traits::{RandomAccessByRef, Shape}; /// The family of an element #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] @@ -70,9 +71,9 @@ pub trait FiniteElement { fn value_size(&self) -> usize; /// Tabulate the values of the basis functions and their derivatives at a set of points - fn tabulate<'a>( + fn tabulate + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, nderivs: usize, data: &mut impl Array4DAccess, ); diff --git a/traits/src/grid.rs b/traits/src/grid.rs index d3c02315..47a85df7 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -1,7 +1,9 @@ //! Geometry and topology definitions -use crate::arrays::{AdjacencyListAccess, Array2DAccess}; +use crate::arrays::AdjacencyListAccess; use crate::cell::ReferenceCellType; +// use crate::element::FiniteElement; +use rlst_common::traits::{RandomAccessByRef, RandomAccessMut, Shape}; /// The ownership of a mesh entity #[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] @@ -18,8 +20,8 @@ pub trait Geometry { /// The geometric dimension fn dim(&self) -> usize; - /// Get the point with the index `i` - fn point(&self, i: usize) -> Option<&[f64]>; + /// Get the coordinates of a point + fn point(&self, index: usize) -> Option>; // TODO: Make this Option<&[f64]> /// The number of points stored in the geometry fn point_count(&self) -> usize; @@ -33,38 +35,58 @@ pub trait Geometry { /// Return the index map from the input cell numbers to the storage numbers fn index_map(&self) -> &[usize]; + /* + /// Compute the physical coordinates of a set of points in a given cell + fn get_compute_points_function< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &self, + element: &impl FiniteElement, + points: &T, + ) -> Box; + */ /// Compute the physical coordinates of a set of points in a given cell - fn compute_points<'a>( + fn compute_points< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - physical_points: &mut impl Array2DAccess<'a, f64>, + physical_points: &mut TMut, ); /// Compute the normals to a set of points in a given cell - fn compute_normals<'a>( + fn compute_normals< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - normals: &mut impl Array2DAccess<'a, f64>, + normals: &mut TMut, ); /// Evaluate the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobians<'a>( + fn compute_jacobians< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobians: &mut impl Array2DAccess<'a, f64>, + jacobians: &mut TMut, ); /// Evaluate the determinand of the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobian_determinants<'a>( + fn compute_jacobian_determinants + Shape>( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, jacobian_determinants: &mut [f64], ); @@ -72,11 +94,14 @@ pub trait Geometry { /// Evaluate the jacobian inverse at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn compute_jacobian_inverses<'a>( + fn compute_jacobian_inverses< + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( &self, - points: &impl Array2DAccess<'a, f64>, + points: &T, cell: usize, - jacobian_inverses: &mut impl Array2DAccess<'a, f64>, + jacobian_inverses: &mut TMut, ); }