From e19af349c4a89744920abe8eb74cba13c1cb98c1 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 22 Nov 2023 11:42:33 +0000 Subject: [PATCH] update to use tensor rlst branch --- Cargo.toml | 8 +- bem/Cargo.toml | 9 +- bem/benches/assembly_benchmark.rs | 11 +- bem/examples/assemble.rs | 2 +- bem/src/assembly.rs | 22 +- bem/src/assembly/batched.rs | 78 +-- bem/src/assembly/cl_kernel.rs | 2 +- element/Cargo.toml | 9 +- element/src/element.rs | 43 +- element/src/element/lagrange.rs | 58 +- element/src/element/raviart_thomas.rs | 26 +- element/src/polynomials.rs | 118 ++-- field/Cargo.toml | 6 +- field/src/fft.rs | 2 +- field/src/field.rs | 32 +- field/src/types.rs | 23 +- fmm/Cargo.toml | 2 +- grid/Cargo.toml | 6 +- grid/examples/curved_cells.rs | 2 +- grid/src/grid.rs | 592 +++++++++---------- grid/src/io.rs | 4 +- grid/src/shapes.rs | 29 +- kernel/Cargo.toml | 6 +- kernel/src/laplace_3d.rs | 105 ++-- tools/Cargo.toml | 3 +- tools/src/arrays.rs | 27 +- traits/Cargo.toml | 4 +- traits/src/element.rs | 2 +- traits/src/grid.rs | 26 +- tree/Cargo.toml | 4 +- tree/examples/test_adaptive.rs | 2 +- tree/examples/test_uniform.rs | 2 +- tree/src/implementations/helpers.rs | 8 +- tree/src/implementations/impl_domain.rs | 4 +- tree/src/implementations/impl_morton.rs | 6 +- tree/src/implementations/impl_single_node.rs | 2 +- 36 files changed, 661 insertions(+), 624 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 466558d0..571fb724 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,8 +3,8 @@ members = [ "bem", "element", - "fmm", - "field", +# "fmm", +# "field", "grid", "hyksort", "quadrature", @@ -19,8 +19,8 @@ resolver = "2" default-members = [ "bem", "element", - "fmm", - "field", +# "fmm", +# "field", "grid", "hyksort", "quadrature", diff --git a/bem/Cargo.toml b/bem/Cargo.toml index f9fc79f6..0ad72832 100644 --- a/bem/Cargo.toml +++ b/bem/Cargo.toml @@ -31,10 +31,11 @@ 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" } +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-blis-src = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-algorithms = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } [dev-dependencies] criterion = { version = "0.3", features = ["html_reports"]} diff --git a/bem/benches/assembly_benchmark.rs b/bem/benches/assembly_benchmark.rs index 9f2c9de3..4ac154ae 100644 --- a/bem/benches/assembly_benchmark.rs +++ b/bem/benches/assembly_benchmark.rs @@ -24,7 +24,7 @@ pub fn full_assembly_benchmark(c: &mut Criterion) { ); let space = SerialFunctionSpace::new(&grid, &element); - let mut matrix = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size())); + let mut matrix = zero_matrix([space.dofmap().global_size(), space.dofmap().global_size()]); group.bench_function( &format!( @@ -62,7 +62,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) { ); let space = SerialFunctionSpace::new(&grid, &element); - let mut matrix = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size())); + let mut matrix = zero_matrix([space.dofmap().global_size(), space.dofmap().global_size()]); let colouring = space.compute_cell_colouring(); @@ -130,7 +130,7 @@ pub fn cl_benchmark(c: &mut Criterion) { ); let space = SerialFunctionSpace::new(&grid, &element); - let mut matrix = zero_matrix((space.dofmap().global_size(), space.dofmap().global_size())); + let mut matrix = zero_matrix([space.dofmap().global_size(), space.dofmap().global_size()]); /* group.bench_function( @@ -181,5 +181,6 @@ pub fn cl_benchmark(c: &mut Criterion) { // criterion_group!(benches, full_assembly_benchmark, assembly_parts_benchmark); criterion_group!(benches, assembly_parts_benchmark); -criterion_group!(cl_benches, cl_benchmark); -criterion_main!(benches, cl_benches); +//criterion_group!(cl_benches, cl_benchmark); +criterion_main!(benches); +//criterion_main!(benches, cl_benches); diff --git a/bem/examples/assemble.rs b/bem/examples/assemble.rs index f6c27ba2..3af85f69 100644 --- a/bem/examples/assemble.rs +++ b/bem/examples/assemble.rs @@ -34,7 +34,7 @@ fn main() { space0.dofmap().global_size() ); let mut matrix = - zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); + zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); println!("Assembling dense matrix (complex)"); assemble_batched( diff --git a/bem/src/assembly.rs b/bem/src/assembly.rs index 37c3e626..65b322d7 100644 --- a/bem/src/assembly.rs +++ b/bem/src/assembly.rs @@ -97,7 +97,7 @@ mod test { use bempp_traits::element::{Continuity, ElementFamily}; // use num::complex::Complex; use bempp_traits::bem::FunctionSpace; - use rlst_dense::RandomAccessByRef; + use rlst_common::traits::RandomAccessByRef; #[test] fn test_laplace_single_layer() { @@ -118,7 +118,7 @@ mod test { let space1 = SerialFunctionSpace::new(&grid, &element1); let mut matrix = - zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); + zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); batched::assemble( &mut matrix, &Laplace3dKernel::new(), @@ -129,7 +129,7 @@ mod test { ); let mut matrix2 = - zero_matrix::((space1.dofmap().global_size(), space0.dofmap().global_size())); + zero_matrix::([space1.dofmap().global_size(), space0.dofmap().global_size()]); assemble_batched( &mut matrix2, @@ -142,8 +142,8 @@ mod test { for i in 0..space1.dofmap().global_size() { for j in 0..space0.dofmap().global_size() { assert_relative_eq!( - *matrix.get(i, j).unwrap(), - *matrix2.get(i, j).unwrap(), + *matrix.get([i, j]).unwrap(), + *matrix2.get([i, j]).unwrap(), epsilon = 0.0001 ); } @@ -163,7 +163,7 @@ mod test { let colouring = space.compute_cell_colouring(); let mut matrix = - zero_matrix::((space.dofmap().global_size(), space.dofmap().global_size())); + zero_matrix::([space.dofmap().global_size(), space.dofmap().global_size()]); batched::assemble_nonsingular::<16, 16>( &mut matrix, &laplace_3d::Laplace3dKernel::new(), @@ -176,7 +176,7 @@ mod test { 128, ); let mut matrix2 = - zero_matrix::((space.dofmap().global_size(), space.dofmap().global_size())); + zero_matrix::([space.dofmap().global_size(), space.dofmap().global_size()]); cl_kernel::assemble( &mut matrix, &Laplace3dKernel::new(), @@ -188,16 +188,16 @@ mod test { for i in 0..5 { for j in 0..5 { - println!("{} {}", *matrix.get(i, j).unwrap(), - *matrix2.get(i, j).unwrap()); + println!("{} {}", *matrix.get([i, j]).unwrap(), + *matrix2.get([i, j]).unwrap()); } println!(); } for i in 0..space.dofmap().global_size() { for j in 0..space.dofmap().global_size() { assert_relative_eq!( - *matrix.get(i, j).unwrap(), - *matrix2.get(i, j).unwrap(), + *matrix.get([i, j]).unwrap(), + *matrix2.get([i, j]).unwrap(), epsilon = 0.0001 ); } diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index 4ad480ce..031c44e4 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -13,7 +13,7 @@ use bempp_traits::kernel::Kernel; use bempp_traits::types::EvalType; use bempp_traits::types::Scalar; use rayon::prelude::*; -use rlst_dense::{RandomAccessByRef, RawAccess, RawAccessMut, Shape}; +use rlst_common::traits::{RandomAccessByRef, RawAccess, RawAccessMut, Shape}; fn get_quadrature_rule( test_celltype: ReferenceCellType, @@ -70,7 +70,7 @@ fn get_quadrature_rule( pub struct RawData2D { pub data: *mut T, - pub shape: (usize, usize), + pub shape: [usize; 2], } unsafe impl Sync for RawData2D {} @@ -95,12 +95,12 @@ fn assemble_batch_singular<'a>( 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 = 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)); + 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 = 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]; @@ -149,13 +149,13 @@ fn assemble_batch_singular<'a>( let mut sum = 0.0; for (index, wt) in weights.iter().enumerate() { - let mut test_row = vec![0.0; test_mapped_pts.shape().1]; + 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(); + *ti = *test_mapped_pts.get([index, i]).unwrap(); } - let mut trial_row = vec![0.0; trial_mapped_pts.shape().1]; + 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(); + *ti = *trial_mapped_pts.get([index, i]).unwrap(); } kernel.assemble_st(EvalType::Value, &test_row, &trial_row, &mut k); @@ -168,7 +168,7 @@ fn assemble_batch_singular<'a>( } unsafe { *output.data.offset( - (*test_dof + output.shape.0 * *trial_dof) + (*test_dof + output.shape[0] * *trial_dof) .try_into() .unwrap(), ) += sum; @@ -204,11 +204,11 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz 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 test_normals = zero_matrix([NPTS_TEST, 3]); + let mut trial_normals = zero_matrix([NPTS_TRIAL, 3]); - let mut test_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; - let mut trial_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TRIAL, 3)]; + let mut test_mapped_pts = rlst_dense::rlst_dynamic_array2![f64, [NPTS_TEST, 3]]; + let mut trial_mapped_pts = rlst_dense::rlst_dynamic_array2![f64, [NPTS_TRIAL, 3]]; let test_element = test_grid.geometry().element(test_cells[0]); let trial_element = trial_grid.geometry().element(trial_cells[0]); @@ -300,7 +300,7 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz if !neighbour { unsafe { *output.data.offset( - (*test_dof + output.shape.0 * *trial_dof) + (*test_dof + output.shape[0] * *trial_dof) .try_into() .unwrap(), ) += sum; @@ -365,8 +365,8 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> 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() + if output.shape()[0] != test_space.dofmap().global_size() + || output.shape()[1] != trial_space.dofmap().global_size() { panic!("Matrix has wrong shape"); } @@ -377,10 +377,10 @@ pub fn assemble_nonsingular<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usize> // TODO: pass cell types into this function let qrule_test = simplex_rule(ReferenceCellType::Triangle, NPTS_TEST).unwrap(); - let qpoints_test = transpose_to_matrix(&qrule_test.points, (NPTS_TEST, 2)); + let qpoints_test = transpose_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 = transpose_to_matrix(&qrule_trial.points, (NPTS_TRIAL, 2)); + let qpoints_trial = transpose_to_matrix(&qrule_trial.points, [NPTS_TRIAL, 2]); let qweights_trial = qrule_trial.weights; let mut test_table = @@ -474,8 +474,8 @@ pub fn assemble_singular<'a>( 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() + if output.shape()[0] != test_space.dofmap().global_size() + || output.shape()[1] != trial_space.dofmap().global_size() { panic!("Matrix has wrong shape"); } @@ -530,21 +530,21 @@ pub fn assemble_singular<'a>( npoints, ); - let points = transpose_to_matrix(&qrule.trial_points, (qrule.npoints, 2)); + let points = transpose_to_matrix(&qrule.trial_points, [qrule.npoints, 2]); let mut table = Array4D::::new( trial_space .element() - .tabulate_array_shape(0, points.shape().0), + .tabulate_array_shape(0, points.shape()[0]), ); trial_space.element().tabulate(&points, 0, &mut table); trial_points.push(points); trial_tables.push(table); - let points = transpose_to_matrix(&qrule.test_points, (qrule.npoints, 2)); + let points = transpose_to_matrix(&qrule.test_points, [qrule.npoints, 2]); let mut table = Array4D::::new( test_space .element() - .tabulate_array_shape(0, points.shape().0), + .tabulate_array_shape(0, points.shape()[0]), ); test_space.element().tabulate(&points, 0, &mut table); test_points.push(points); @@ -638,7 +638,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = zero_matrix::((ndofs, ndofs)); + let mut matrix = zero_matrix::([ndofs, ndofs]); assemble( &mut matrix, &Laplace3dKernel::new(), @@ -654,7 +654,7 @@ mod test { 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); + assert_relative_eq!(*matrix.get([i, j]).unwrap(), entry, epsilon = 1e-3); } } } @@ -673,7 +673,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = Array2D::::new([ndofs, ndofs]); assemble( &mut matrix, &green::LaplaceGreenDyKernel {}, @@ -707,7 +707,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = Array2D::::new([ndofs, ndofs]); assemble( &mut matrix, &green::LaplaceGreenDxKernel {}, @@ -741,7 +741,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = Array2D::::new([ndofs, ndofs]); laplace_hypersingular_assemble(&mut matrix, &space, &space); for i in 0..ndofs { @@ -764,7 +764,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = Array2D::::new([ndofs, ndofs]); laplace_hypersingular_assemble(&mut matrix, &space, &space); @@ -798,7 +798,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::::new((ndofs, ndofs)); + let mut matrix = Array2D::::new([ndofs, ndofs]); assemble( &mut matrix, &green::HelmholtzGreenKernel { k: 3.0 }, @@ -831,7 +831,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::>::new((ndofs, ndofs)); + let mut matrix = Array2D::>::new([ndofs, ndofs]); assemble( &mut matrix, &green::HelmholtzGreenKernel { k: 3.0 }, @@ -865,7 +865,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::>::new((ndofs, ndofs)); + let mut matrix = Array2D::>::new([ndofs, ndofs]); assemble( &mut matrix, &green::HelmholtzGreenDyKernel { k: 3.0 }, @@ -900,7 +900,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::>::new((ndofs, ndofs)); + let mut matrix = Array2D::>::new([ndofs, ndofs]); assemble( &mut matrix, &green::HelmholtzGreenDxKernel { k: 3.0 }, @@ -935,7 +935,7 @@ mod test { let ndofs = space.dofmap().global_size(); - let mut matrix = Array2D::>::new((ndofs, ndofs)); + let mut matrix = Array2D::>::new([ndofs, ndofs]); helmholtz_hypersingular_assemble(&mut matrix, &space, &space, 3.0); diff --git a/bem/src/assembly/cl_kernel.rs b/bem/src/assembly/cl_kernel.rs index 483cf998..eedae2c7 100644 --- a/bem/src/assembly/cl_kernel.rs +++ b/bem/src/assembly/cl_kernel.rs @@ -14,7 +14,7 @@ use bempp_traits::kernel::Kernel; use bempp_traits::types::EvalType; use bempp_traits::types::Scalar; use rayon::prelude::*; -use rlst_dense::{RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, Shape}; +use rlst_common::traits::{RandomAccessByRef, RandomAccessMut, RawAccess, RawAccessMut, Shape}; fn get_corners<'a>(grid: &impl Grid<'a>, index: &usize, corners: &mut Vec>) { let v = grid.geometry().cell_vertices(*index).unwrap(); diff --git a/element/Cargo.toml b/element/Cargo.toml index e7afe14f..849c18b2 100644 --- a/element/Cargo.toml +++ b/element/Cargo.toml @@ -26,8 +26,9 @@ bempp-quadrature = { path = "../quadrature" } paste = "1.*" libc = "0.2" approx = "0.5" -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" } +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-blis-src = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-algorithms = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } diff --git a/element/src/element.rs b/element/src/element.rs index 32269d33..ff702211 100644 --- a/element/src/element.rs +++ b/element/src/element.rs @@ -6,9 +6,9 @@ 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::{rlst_dynamic_mat, RandomAccessByRef, RandomAccessMut, Shape}; +use rlst_dense::linalg::Trans; +use rlst_dense::rlst_dynamic_array2; +use rlst_common::traits::{ RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessMut}; pub mod lagrange; pub mod raviart_thomas; @@ -71,20 +71,20 @@ impl CiarletElement { let new_pts = if continuity == Continuity::Discontinuous { let mut new_pts: [Vec>; 4] = [vec![], vec![], vec![], vec![]]; let mut pn = 0; - let mut all_pts = rlst_dynamic_mat![f64, (npts, tdim)]; + let mut all_pts = rlst_dynamic_array2![f64, [npts, tdim]]; for (i, pts_i) in interpolation_points.iter().take(tdim).enumerate() { for _pts in pts_i { - new_pts[i].push(rlst_dynamic_mat![f64, (0, tdim)]); + new_pts[i].push(rlst_dynamic_array2![f64, [0, tdim]]); } } for pts_i in interpolation_points.iter() { for pts in pts_i { - for j in 0..pts.shape().0 { + for j in 0..pts.shape()[0] { for k in 0..tdim { - *all_pts.get_mut(pn + j, k).unwrap() = *pts.get(j, k).unwrap(); + *all_pts.get_mut([pn + j, k]).unwrap() = *pts.get([j, k]).unwrap(); } } - pn += pts.shape().0; + pn += pts.shape()[0]; } } new_pts[tdim].push(all_pts); @@ -129,8 +129,8 @@ impl CiarletElement { let mut dof = 0; for d in 0..4 { for (e, pts) in new_pts[d].iter().enumerate() { - if pts.shape().0 > 0 { - let mut table = Array3D::::new((1, pdim, pts.shape().0)); + if pts.shape()[0] > 0 { + let mut table = Array3D::::new((1, pdim, pts.shape()[0])); tabulate_legendre_polynomials(cell_type, pts, highest_degree, 0, &mut table); let mat = &new_wts[d][e]; for i in 0..mat.shape().0 { @@ -138,7 +138,7 @@ impl CiarletElement { for l in 0..pdim { let value = d_matrix.get_mut(j, l, dof + i).unwrap(); *value = 0.0; - for k in 0..pts.shape().0 { + for k in 0..pts.shape()[0] { *value += *mat.get(i, j, k).unwrap() * *table.get(0, l, k).unwrap(); } @@ -150,11 +150,11 @@ impl CiarletElement { } } - let mut dual_matrix = rlst_dense::rlst_dynamic_mat!(f64, (dim, dim)); + let mut dual_matrix = rlst_dense::rlst_dynamic_array2!(f64, [dim, dim]); for i in 0..dim { for j in 0..dim { - let entry = dual_matrix.get_mut(i, j).unwrap(); + let entry = dual_matrix.get_mut([i, j]).unwrap(); *entry = 0.0; for k in 0..value_size { for l in 0..pdim { @@ -165,7 +165,12 @@ impl CiarletElement { } } - let inverse = dual_matrix.linalg().inverse().unwrap(); + let mut ident = rlst_dense::rlst_dynamic_array2!(f64, [dim, dim]); + for i in 0..dim { unsafe { + *ident.get_unchecked_mut([i, i]) = 1.0; + }} + let lu = dual_matrix.into_lu().unwrap(); + let inverse = lu.solve(Trans::NoTrans, ident).unwrap(); let mut coefficients = Array3D::::new((dim, value_size, pdim)); for i in 0..dim { @@ -173,7 +178,7 @@ impl CiarletElement { for j in 0..value_size { for k in 0..pdim { *coefficients.get_mut(i, j, k).unwrap() += - *inverse.get(i, l).unwrap() * *polynomial_coeffs.get(l, j, k).unwrap() + *inverse.get([i, l]).unwrap() * *polynomial_coeffs.get(l, j, k).unwrap() } } } @@ -188,9 +193,9 @@ impl CiarletElement { let mut dof = 0; for i in 0..4 { for pts in &new_pts[i] { - let dofs: Vec = (dof..dof + pts.shape().0).collect(); + let dofs: Vec = (dof..dof + pts.shape()[0]).collect(); entity_dofs[i].add_row(&dofs); - dof += pts.shape().0; + dof += pts.shape()[0]; } } CiarletElement { @@ -240,7 +245,7 @@ impl FiniteElement for CiarletElement { fn dim(&self) -> usize { self.dim } - fn tabulate + Shape>( + fn tabulate + Shape<2>>( &self, points: &T, nderivs: usize, @@ -261,7 +266,7 @@ impl FiniteElement for CiarletElement { ); for d in 0..table.shape().0 { - for p in 0..points.shape().0 { + for p in 0..points.shape()[0] { for j in 0..self.value_size { for b in 0..self.dim { let value = data.get_mut(d, p, b, j).unwrap(); diff --git a/element/src/element/lagrange.rs b/element/src/element/lagrange.rs index 1ecf923d..8351d9fb 100644 --- a/element/src/element/lagrange.rs +++ b/element/src/element/lagrange.rs @@ -6,7 +6,7 @@ use bempp_tools::arrays::{to_matrix, zero_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; -use rlst_dense::RandomAccessMut; +use rlst_common::traits::RandomAccessMut; /// Create a Lagrange element pub fn create( @@ -30,24 +30,24 @@ pub fn create( } for d in 0..tdim { for _e in 0..cell.entity_count(d) { - x[d].push(zero_matrix((0, tdim))); + x[d].push(zero_matrix([0, tdim])); m[d].push(Array3D::::new((0, 1, 0))); } } - x[tdim].push(to_matrix(&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 = zero_matrix((1, tdim)); + let mut pts = zero_matrix([1, tdim]); for i in 0..tdim { - *pts.get_mut(0, i).unwrap() = cell.vertices()[e * tdim + i]; + *pts.get_mut([0, i]).unwrap() = cell.vertices()[e * tdim + i]; } x[0].push(pts); m[0].push(Array3D::::from_data(vec![1.0], (1, 1, 1))); } for e in 0..cell.entity_count(1) { - let mut pts = zero_matrix((degree - 1, tdim)); + let mut pts = zero_matrix([degree - 1, tdim]); let mut ident = vec![0.0; (degree - 1).pow(2)]; let vn0 = cell.edges()[2 * e]; let vn1 = cell.edges()[2 * e + 1]; @@ -56,7 +56,7 @@ pub fn create( for i in 1..degree { ident[(i - 1) * degree] = 1.0; for j in 0..tdim { - *pts.get_mut(i - 1, j).unwrap() = + *pts.get_mut([i - 1, j]).unwrap() = v0[j] + i as f64 / degree as f64 * (v1[j] - v0[j]); } } @@ -80,7 +80,7 @@ pub fn create( } else { panic!("Unsupported face type"); }; - let mut pts = zero_matrix((npts, tdim)); + let mut pts = zero_matrix([npts, tdim]); let mut ident = vec![0.0; npts.pow(2)]; let vn0 = cell.faces()[start]; @@ -96,7 +96,7 @@ pub fn create( for i0 in 1..degree { for i1 in 1..degree - i0 { for j in 0..tdim { - *pts.get_mut(n, j).unwrap() = v0[j] + *pts.get_mut([n, j]).unwrap() = v0[j] + i0 as f64 / degree as f64 * (v1[j] - v0[j]) + i1 as f64 / degree as f64 * (v2[j] - v0[j]); } @@ -109,7 +109,7 @@ pub fn create( for i0 in 1..degree { for i1 in 1..degree { for j in 0..tdim { - *pts.get_mut(n, j).unwrap() = v0[j] + *pts.get_mut([n, j]).unwrap() = v0[j] + i0 as f64 / degree as f64 * (v1[j] - v0[j]) + i1 as f64 / degree as f64 * (v2[j] - v0[j]); } @@ -150,7 +150,7 @@ mod test { use bempp_tools::arrays::Array4D; use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; - use rlst_dense::RandomAccessByRef; + use rlst_common::traits::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -187,7 +187,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 = to_matrix(&[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 { @@ -201,15 +201,15 @@ 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 = to_matrix(&[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 { assert_relative_eq!( *data.get(0, pt, 0, 0).unwrap(), - 1.0 - *points.get(pt, 0).unwrap() + 1.0 - *points.get([pt, 0]).unwrap() ); - assert_relative_eq!(*data.get(0, pt, 1, 0).unwrap(), *points.get(pt, 0).unwrap()); + assert_relative_eq!(*data.get(0, pt, 1, 0).unwrap(), *points.get([pt, 0]).unwrap()); } check_dofs(e); } @@ -221,7 +221,7 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); @@ -238,17 +238,17 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); for pt in 0..6 { assert_relative_eq!( *data.get(0, pt, 0, 0).unwrap(), - 1.0 - *points.get(pt, 0).unwrap() - *points.get(pt, 1).unwrap() + 1.0 - *points.get([pt, 0]).unwrap() - *points.get([pt, 1]).unwrap() ); - assert_relative_eq!(*data.get(0, pt, 1, 0).unwrap(), *points.get(pt, 0).unwrap()); - assert_relative_eq!(*data.get(0, pt, 2, 0).unwrap(), *points.get(pt, 1).unwrap()); + assert_relative_eq!(*data.get(0, pt, 1, 0).unwrap(), *points.get([pt, 0]).unwrap()); + assert_relative_eq!(*data.get(0, pt, 2, 0).unwrap(), *points.get([pt, 1]).unwrap()); } check_dofs(e); } @@ -319,7 +319,7 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); @@ -336,26 +336,26 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); for pt in 0..6 { assert_relative_eq!( *data.get(0, pt, 0, 0).unwrap(), - (1.0 - *points.get(pt, 0).unwrap()) * (1.0 - *points.get(pt, 1).unwrap()) + (1.0 - *points.get([pt, 0]).unwrap()) * (1.0 - *points.get([pt, 1]).unwrap()) ); assert_relative_eq!( *data.get(0, pt, 1, 0).unwrap(), - *points.get(pt, 0).unwrap() * (1.0 - *points.get(pt, 1).unwrap()) + *points.get([pt, 0]).unwrap() * (1.0 - *points.get([pt, 1]).unwrap()) ); assert_relative_eq!( *data.get(0, pt, 2, 0).unwrap(), - (1.0 - *points.get(pt, 0).unwrap()) * *points.get(pt, 1).unwrap() + (1.0 - *points.get([pt, 0]).unwrap()) * *points.get([pt, 1]).unwrap() ); assert_relative_eq!( *data.get(0, pt, 3, 0).unwrap(), - *points.get(pt, 0).unwrap() * *points.get(pt, 1).unwrap() + *points.get([pt, 0]).unwrap() * *points.get([pt, 1]).unwrap() ); } check_dofs(e); @@ -368,13 +368,13 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 1.0, 0.25, 0.3, 0.0, 0.0, 1.0, 1.0, 0.5, 0.2], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); for pt in 0..6 { - let x = *points.get(pt, 0).unwrap(); - let y = *points.get(pt, 1).unwrap(); + let x = *points.get([pt, 0]).unwrap(); + let y = *points.get([pt, 1]).unwrap(); assert_relative_eq!( *data.get(0, pt, 0, 0).unwrap(), (1.0 - x) * (1.0 - 2.0 * x) * (1.0 - y) * (1.0 - 2.0 * y) diff --git a/element/src/element/raviart_thomas.rs b/element/src/element/raviart_thomas.rs index 7bb60551..22785aaf 100644 --- a/element/src/element/raviart_thomas.rs +++ b/element/src/element/raviart_thomas.rs @@ -6,7 +6,7 @@ use bempp_tools::arrays::{zero_matrix, Array3D}; use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; use bempp_traits::element::{Continuity, ElementFamily, MapType}; -use rlst_dense::RandomAccessMut; +use rlst_common::traits::RandomAccessMut; /// Create a Raviart-Thomas element pub fn create( @@ -46,19 +46,19 @@ 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(zero_matrix((0, tdim))); + x[0].push(zero_matrix([0, tdim])); m[0].push(Array3D::::new((0, 2, 0))); } for e in 0..cell.entity_count(1) { - let mut pts = zero_matrix((1, tdim)); + let mut pts = zero_matrix([1, tdim]); let mut mat = vec![0.0; 2]; let vn0 = cell.edges()[2 * e]; let vn1 = cell.edges()[2 * e + 1]; let v0 = &cell.vertices()[vn0 * tdim..(vn0 + 1) * tdim]; let v1 = &cell.vertices()[vn1 * tdim..(vn1 + 1) * tdim]; for i in 0..tdim { - *pts.get_mut(0, i).unwrap() = (v0[i] + v1[i]) / 2.0; + *pts.get_mut([0, i]).unwrap() = (v0[i] + v1[i]) / 2.0; } mat[0] = v0[1] - v1[1]; mat[1] = v1[0] - v0[0]; @@ -67,7 +67,7 @@ pub fn create( } for _e in 0..cell.entity_count(2) { - x[2].push(zero_matrix((0, tdim))); + x[2].push(zero_matrix([0, tdim])); m[2].push(Array3D::::new((0, 2, 0))); } @@ -93,7 +93,7 @@ mod test { use bempp_tools::arrays::{to_matrix, Array4D}; use bempp_traits::arrays::Array4DAccess; use bempp_traits::element::FiniteElement; - use rlst_dense::RandomAccessByRef; + use rlst_common::traits::RandomAccessByRef; fn check_dofs(e: impl FiniteElement) { let cell_dim = match e.cell_type() { @@ -132,31 +132,31 @@ mod test { let mut data = Array4D::::new(e.tabulate_array_shape(0, 6)); let points = to_matrix( &[0.0, 1.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.5], - (6, 2), + [6, 2], ); e.tabulate(&points, 0, &mut data); for pt in 0..6 { assert_relative_eq!( *data.get(0, pt, 0, 0).unwrap(), - -*points.get(pt, 0).unwrap() + -*points.get([pt, 0]).unwrap() ); assert_relative_eq!( *data.get(0, pt, 0, 1).unwrap(), - -*points.get(pt, 1).unwrap() + -*points.get([pt, 1]).unwrap() ); assert_relative_eq!( *data.get(0, pt, 1, 0).unwrap(), - *points.get(pt, 0).unwrap() - 1.0 + *points.get([pt, 0]).unwrap() - 1.0 ); - assert_relative_eq!(*data.get(0, pt, 1, 1).unwrap(), *points.get(pt, 1).unwrap()); + assert_relative_eq!(*data.get(0, pt, 1, 1).unwrap(), *points.get([pt, 1]).unwrap()); assert_relative_eq!( *data.get(0, pt, 2, 0).unwrap(), - -*points.get(pt, 0).unwrap() + -*points.get([pt, 0]).unwrap() ); assert_relative_eq!( *data.get(0, pt, 2, 1).unwrap(), - 1.0 - *points.get(pt, 1).unwrap() + 1.0 - *points.get([pt, 1]).unwrap() ); } check_dofs(e); diff --git a/element/src/polynomials.rs b/element/src/polynomials.rs index 677c46d1..0b542348 100644 --- a/element/src/polynomials.rs +++ b/element/src/polynomials.rs @@ -2,10 +2,10 @@ use bempp_traits::arrays::Array3DAccess; use bempp_traits::cell::ReferenceCellType; -use rlst_dense::{RandomAccessByRef, Shape}; +use rlst_common::traits::{RandomAccessByRef, Shape}; /// Tabulate orthonormal polynomials on a interval -fn tabulate_legendre_polynomials_interval + Shape>( +fn tabulate_legendre_polynomials_interval + Shape<2>>( points: &T, degree: usize, derivatives: usize, @@ -13,8 +13,8 @@ fn tabulate_legendre_polynomials_interval + Sha ) { assert_eq!(data.shape().0, derivatives + 1); assert_eq!(data.shape().1, degree + 1); - assert_eq!(data.shape().2, points.shape().0); - assert_eq!(points.shape().1, 1); + assert_eq!(data.shape().2, points.shape()[0]); + assert_eq!(points.shape()[1], 1); for i in 0..data.shape().2 { *data.get_mut(0, 0, i).unwrap() = 1.0; @@ -31,7 +31,7 @@ fn tabulate_legendre_polynomials_interval + Sha let b = (a + 1.0) * ((2.0 * p as f64 + 1.0) / (2.0 * p as f64 - 1.0)).sqrt(); for i in 0..data.shape().2 { *data.get_mut(k, p, i).unwrap() = - (points.get(i, 0).unwrap() * 2.0 - 1.0) * data.get(k, p - 1, i).unwrap() * b; + (points.get([i, 0]).unwrap() * 2.0 - 1.0) * data.get(k, p - 1, i).unwrap() * b; } if p > 1 { let c = a * ((2.0 * p as f64 + 1.0) / (2.0 * p as f64 - 3.0)).sqrt(); @@ -58,7 +58,7 @@ fn quad_index(i: usize, j: usize, n: usize) -> usize { } /// Tabulate orthonormal polynomials on a quadrilateral -fn tabulate_legendre_polynomials_quadrilateral + Shape>( +fn tabulate_legendre_polynomials_quadrilateral + Shape<2>>( points: &T, degree: usize, derivatives: usize, @@ -66,8 +66,8 @@ fn tabulate_legendre_polynomials_quadrilateral ) { assert_eq!(data.shape().0, (derivatives + 1) * (derivatives + 2) / 2); assert_eq!(data.shape().1, (degree + 1) * (degree + 1)); - assert_eq!(data.shape().2, points.shape().0); - assert_eq!(points.shape().1, 2); + assert_eq!(data.shape().2, points.shape()[0]); + assert_eq!(points.shape()[1], 2); for i in 0..data.shape().2 { *data @@ -91,7 +91,7 @@ fn tabulate_legendre_polynomials_quadrilateral for i in 0..data.shape().2 { *data .get_mut(tri_index(k, 0), quad_index(p, 0, degree), i) - .unwrap() = (points.get(i, 0).unwrap() * 2.0 - 1.0) + .unwrap() = (points.get([i, 0]).unwrap() * 2.0 - 1.0) * data .get(tri_index(k, 0), quad_index(p - 1, 0, degree), i) .unwrap() @@ -139,7 +139,7 @@ fn tabulate_legendre_polynomials_quadrilateral for i in 0..data.shape().2 { *data .get_mut(tri_index(0, k), quad_index(0, p, degree), i) - .unwrap() = (points.get(i, 1).unwrap() * 2.0 - 1.0) + .unwrap() = (points.get([i, 1]).unwrap() * 2.0 - 1.0) * data .get(tri_index(0, k), quad_index(0, p - 1, degree), i) .unwrap() @@ -192,7 +192,7 @@ fn tabulate_legendre_polynomials_quadrilateral } } /// Tabulate orthonormal polynomials on a triangle -fn tabulate_legendre_polynomials_triangle + Shape>( +fn tabulate_legendre_polynomials_triangle + Shape<2>>( points: &T, degree: usize, derivatives: usize, @@ -200,8 +200,8 @@ fn tabulate_legendre_polynomials_triangle + Sha ) { assert_eq!(data.shape().0, (derivatives + 1) * (derivatives + 2) / 2); assert_eq!(data.shape().1, (degree + 1) * (degree + 2) / 2); - assert_eq!(data.shape().2, points.shape().0); - assert_eq!(points.shape().1, 2); + assert_eq!(data.shape().2, points.shape()[0]); + assert_eq!(points.shape()[1], 2); for i in 0..data.shape().2 { *data.get_mut(tri_index(0, 0), tri_index(0, 0), i).unwrap() = f64::sqrt(2.0); @@ -221,7 +221,7 @@ fn tabulate_legendre_polynomials_triangle + Sha f64::sqrt((p as f64 + 0.5) * (p as f64 + 1.0) / ((p as f64 - 0.5) * p as f64)); for i in 0..data.shape().2 { *data.get_mut(tri_index(kx, ky), tri_index(0, p), i).unwrap() = - (*points.get(i, 0).unwrap() * 2.0 + *points.get(i, 1).unwrap() - 1.0) + (*points.get([i, 0]).unwrap() * 2.0 + *points.get([i, 1]).unwrap() - 1.0) * *data.get(tri_index(kx, ky), tri_index(0, p - 1), i).unwrap() * a * scale1; @@ -252,7 +252,7 @@ fn tabulate_legendre_polynomials_triangle + Sha / f64::sqrt((p as f64 - 1.5) * (p as f64 - 1.0)); for i in 0..data.shape().2 { - let b = 1.0 - *points.get(i, 1).unwrap(); + let b = 1.0 - *points.get([i, 1]).unwrap(); *data.get_mut(tri_index(kx, ky), tri_index(0, p), i).unwrap() -= b * b * *data.get(tri_index(kx, ky), tri_index(0, p - 2), i).unwrap() @@ -263,7 +263,7 @@ fn tabulate_legendre_polynomials_triangle + Sha for i in 0..data.shape().2 { *data.get_mut(tri_index(kx, ky), tri_index(0, p), i).unwrap() -= 2.0 * ky as f64 - * (*points.get(i, 1).unwrap() - 1.0) + * (*points.get([i, 1]).unwrap() - 1.0) * *data .get(tri_index(kx, ky - 1), tri_index(0, p - 2), i) .unwrap() @@ -291,7 +291,7 @@ fn tabulate_legendre_polynomials_triangle + Sha *data.get_mut(tri_index(kx, ky), tri_index(1, p), i).unwrap() = *data.get(tri_index(kx, ky), tri_index(0, p), i).unwrap() * scale3 - * ((*points.get(i, 1).unwrap() * 2.0 - 1.0) * (1.5 + p as f64) + * ((*points.get([i, 1]).unwrap() * 2.0 - 1.0) * (1.5 + p as f64) + 0.5 + p as f64); } @@ -321,7 +321,7 @@ fn tabulate_legendre_polynomials_triangle + Sha .unwrap() = *data.get_mut(tri_index(kx, ky), tri_index(q, p), i).unwrap() * scale4 - * ((*points.get(i, 1).unwrap() * 2.0 - 1.0) * a1 + a2) + * ((*points.get([i, 1]).unwrap() * 2.0 - 1.0) * a1 + a2) - *data .get_mut(tri_index(kx, ky), tri_index(q - 1, p), i) .unwrap() @@ -369,7 +369,7 @@ pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usi } } -pub fn legendre_shape + Shape>( +pub fn legendre_shape + Shape<2>>( cell_type: ReferenceCellType, points: &T, degree: usize, @@ -378,12 +378,12 @@ pub fn legendre_shape + Shape>( ( derivative_count(cell_type, derivatives), polynomial_count(cell_type, degree), - points.shape().0, + points.shape()[0], ) } /// Tabulate orthonormal polynomials -pub fn tabulate_legendre_polynomials + Shape>( +pub fn tabulate_legendre_polynomials + Shape<2>>( cell_type: ReferenceCellType, points: &T, degree: usize, @@ -412,14 +412,14 @@ mod test { use approx::*; use bempp_quadrature::simplex_rules::simplex_rule; use bempp_tools::arrays::{transpose_to_matrix, zero_matrix, Array3D}; - use rlst_dense::RandomAccessMut; + use rlst_common::traits::RandomAccessMut; #[test] fn test_legendre_interval() { let degree = 6; let rule = simplex_rule(ReferenceCellType::Interval, degree + 1).unwrap(); - let points = transpose_to_matrix(&rule.points, (rule.npoints, 1)); + let points = transpose_to_matrix(&rule.points, [rule.npoints, 1]); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Interval, @@ -450,7 +450,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Triangle, 79).unwrap(); - let points = transpose_to_matrix(&rule.points, (rule.npoints, 2)); + let points = transpose_to_matrix(&rule.points, [rule.npoints, 2]); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Triangle, @@ -481,7 +481,7 @@ mod test { let degree = 5; let rule = simplex_rule(ReferenceCellType::Quadrilateral, 85).unwrap(); - let points = transpose_to_matrix(&rule.points, (rule.npoints, 2)); + let points = transpose_to_matrix(&rule.points, [rule.npoints, 2]); let mut data = Array3D::::new(legendre_shape( ReferenceCellType::Quadrilateral, @@ -518,10 +518,10 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix((20, 1)); + let mut points = zero_matrix([20, 1]); for i in 0..10 { - *points.get_mut(2 * i, 0).unwrap() = i as f64 / 10.0; - *points.get_mut(2 * i + 1, 0).unwrap() = points.get(2 * i, 0).unwrap() + epsilon; + *points.get_mut([2 * i, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([2 * i + 1, 0]).unwrap() = points.get([2 * i, 0]).unwrap() + epsilon; } let mut data = Array3D::::new(legendre_shape( @@ -533,7 +533,7 @@ mod test { tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 1, &mut data); for i in 0..degree + 1 { - for k in 0..points.shape().0 / 2 { + for k in 0..points.shape()[0] / 2 { assert_relative_eq!( *data.get(1, i, 2 * k).unwrap(), (data.get(0, i, 2 * k + 1).unwrap() - data.get(0, i, 2 * k).unwrap()) / epsilon, @@ -548,18 +548,18 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix((165, 2)); + let mut points = zero_matrix([165, 2]); let mut index = 0; for i in 0..10 { for j in 0..10 - i { - *points.get_mut(3 * index, 0).unwrap() = i as f64 / 10.0; - *points.get_mut(3 * index, 1).unwrap() = j as f64 / 10.0; - *points.get_mut(3 * index + 1, 0).unwrap() = - *points.get(3 * index, 0).unwrap() + epsilon; - *points.get_mut(3 * index + 1, 1).unwrap() = *points.get(3 * index, 1).unwrap(); - *points.get_mut(3 * index + 2, 0).unwrap() = *points.get(3 * index, 0).unwrap(); - *points.get_mut(3 * index + 2, 1).unwrap() = - *points.get(3 * index, 1).unwrap() + epsilon; + *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0; + *points.get_mut([3 * index + 1, 0]).unwrap() = + *points.get([3 * index, 0]).unwrap() + epsilon; + *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap(); + *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap(); + *points.get_mut([3 * index + 2, 1]).unwrap() = + *points.get([3 * index, 1]).unwrap() + epsilon; index += 1; } } @@ -573,7 +573,7 @@ mod test { tabulate_legendre_polynomials(ReferenceCellType::Triangle, &points, degree, 1, &mut data); for i in 0..degree + 1 { - for k in 0..points.shape().0 / 3 { + for k in 0..points.shape()[0] / 3 { assert_relative_eq!( *data.get(1, i, 3 * k).unwrap(), (data.get(0, i, 3 * k + 1).unwrap() - data.get(0, i, 3 * k).unwrap()) / epsilon, @@ -593,18 +593,18 @@ mod test { let degree = 6; let epsilon = 1e-10; - let mut points = zero_matrix((300, 2)); + let mut points = zero_matrix([300, 2]); for i in 0..10 { for j in 0..10 { let index = 10 * i + j; - *points.get_mut(3 * index, 0).unwrap() = i as f64 / 10.0; - *points.get_mut(3 * index, 1).unwrap() = j as f64 / 10.0; - *points.get_mut(3 * index + 1, 0).unwrap() = - *points.get(3 * index, 0).unwrap() + epsilon; - *points.get_mut(3 * index + 1, 1).unwrap() = *points.get(3 * index, 1).unwrap(); - *points.get_mut(3 * index + 2, 0).unwrap() = *points.get(3 * index, 0).unwrap(); - *points.get_mut(3 * index + 2, 1).unwrap() = - *points.get(3 * index, 1).unwrap() + epsilon; + *points.get_mut([3 * index, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([3 * index, 1]).unwrap() = j as f64 / 10.0; + *points.get_mut([3 * index + 1, 0]).unwrap() = + *points.get([3 * index, 0]).unwrap() + epsilon; + *points.get_mut([3 * index + 1, 1]).unwrap() = *points.get([3 * index, 1]).unwrap(); + *points.get_mut([3 * index + 2, 0]).unwrap() = *points.get([3 * index, 0]).unwrap(); + *points.get_mut([3 * index + 2, 1]).unwrap() = + *points.get([3 * index, 1]).unwrap() + epsilon; } } @@ -623,7 +623,7 @@ mod test { ); for i in 0..degree + 1 { - for k in 0..points.shape().0 / 3 { + for k in 0..points.shape()[0] / 3 { assert_relative_eq!( *data.get(1, i, 3 * k).unwrap(), (data.get(0, i, 3 * k + 1).unwrap() - data.get(0, i, 3 * k).unwrap()) / epsilon, @@ -642,9 +642,9 @@ mod test { fn test_legendre_interval_against_known_polynomials() { let degree = 3; - let mut points = zero_matrix((11, 1)); + let mut points = zero_matrix([11, 1]); for i in 0..11 { - *points.get_mut(i, 0).unwrap() = i as f64 / 10.0; + *points.get_mut([i, 0]).unwrap() = i as f64 / 10.0; } let mut data = Array3D::::new(legendre_shape( @@ -655,8 +655,8 @@ mod test { )); tabulate_legendre_polynomials(ReferenceCellType::Interval, &points, degree, 3, &mut data); - for k in 0..points.shape().0 { - let x = *points.get(k, 0).unwrap(); + for k in 0..points.shape()[0] { + let x = *points.get([k, 0]).unwrap(); // 0 => 1 assert_relative_eq!(*data.get(0, 0, k).unwrap(), 1.0, epsilon = 1e-12); @@ -724,11 +724,11 @@ mod test { fn test_legendre_quadrilateral_against_known_polynomials() { let degree = 2; - let mut points = zero_matrix((121, 2)); + let mut points = zero_matrix([121, 2]); for i in 0..11 { for j in 0..11 { - *points.get_mut(11 * i + j, 0).unwrap() = i as f64 / 10.0; - *points.get_mut(11 * i + j, 1).unwrap() = j as f64 / 10.0; + *points.get_mut([11 * i + j, 0]).unwrap() = i as f64 / 10.0; + *points.get_mut([11 * i + j, 1]).unwrap() = j as f64 / 10.0; } } @@ -746,9 +746,9 @@ mod test { &mut data, ); - for k in 0..points.shape().0 { - let x = *points.get(k, 0).unwrap(); - let y = *points.get(k, 1).unwrap(); + for k in 0..points.shape()[0] { + let x = *points.get([k, 0]).unwrap(); + let y = *points.get([k, 1]).unwrap(); // 0 => 1 assert_relative_eq!(*data.get(0, 0, k).unwrap(), 1.0, epsilon = 1e-12); diff --git a/field/Cargo.toml b/field/Cargo.toml index 29c198b7..0cf07533 100644 --- a/field/Cargo.toml +++ b/field/Cargo.toml @@ -28,11 +28,13 @@ bempp-kernel = { path = "../kernel" } bempp-tools = { path = "../tools" } itertools = "0.10" num = "0.4" -rlst = { git = "https://github.com/linalg-rs/rlst.git" } +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } fftw = {git = "https://github.com/skailasa/fftw.git" } cauchy = "0.4.*" dashmap = {version = "5.5.0", features=["rayon"]} rayon = "1.7" [dev-dependencies] -approx_eq = "0.1.8" \ No newline at end of file +approx_eq = "0.1.8" diff --git a/field/src/fft.rs b/field/src/fft.rs index 021df0c6..8efc0d6f 100644 --- a/field/src/fft.rs +++ b/field/src/fft.rs @@ -4,7 +4,7 @@ use num::Complex; use rayon::prelude::*; use crate::types::{FftMatrixc32, FftMatrixc64, FftMatrixf32, FftMatrixf64}; -use rlst::dense::RawAccessMut; +use rlst_common::traits::RawAccessMut; pub trait Fft where diff --git a/field/src/field.rs b/field/src/field.rs index 3f7a9acc..9874eb2a 100644 --- a/field/src/field.rs +++ b/field/src/field.rs @@ -4,17 +4,11 @@ use itertools::Itertools; use num::Zero; use num::{Complex, Float}; use rlst::dense::LayoutType; -use rlst::{ - algorithms::{ - linalg::{DenseMatrixLinAlgBuilder, LinAlg}, - traits::svd::{Mode, Svd}, - }, - common::traits::{Eval, Transpose}, - dense::{ - rlst_dynamic_mat, rlst_pointer_mat, Dot, Dynamic, MultiplyAdd, RawAccess, RawAccessMut, +use rlst_common::traits::{Eval, Transpose}; +use rlst_dense::{ + rlst_dynamic_array2, rlst_pointer_mat, Dot, Dynamic, MultiplyAdd, RawAccess, RawAccessMut, Shape, VectorContainer, - }, -}; + }; use std::collections::{HashMap, HashSet}; use bempp_tools::Array3D; @@ -68,8 +62,8 @@ where let ncols = self.ncoeffs(order); let ntransfer_vectors = self.transfer_vectors.len(); - let mut se2tc_fat = rlst_dynamic_mat![T, (nrows, ncols * ntransfer_vectors)]; - let mut se2tc_thin = rlst_dynamic_mat![T, (nrows * ntransfer_vectors, ncols)]; + let mut se2tc_fat = rlst_dynamic_array2![T, [nrows, ncols * ntransfer_vectors]]; + let mut se2tc_thin = rlst_dynamic_array2![T, [nrows * ntransfer_vectors, ncols]]; for (i, t) in self.transfer_vectors.iter().enumerate() { let source_equivalent_surface = t.source.compute_surface(&domain, order, self.alpha); @@ -78,7 +72,7 @@ where let target_check_surface = t.target.compute_surface(&domain, order, self.alpha); let ntargets = target_check_surface.len() / self.kernel.space_dimension(); - let mut tmp_gram = rlst_dynamic_mat![T, (ntargets, nsources)]; + let mut tmp_gram = rlst_dynamic_array2![T, [ntargets, nsources]]; self.kernel.assemble_st( EvalType::Value, @@ -111,7 +105,7 @@ where let vt = vt.unwrap(); // Keep 'k' singular values - let mut sigma_mat = rlst_dynamic_mat![T, (self.k, self.k)]; + let mut sigma_mat = rlst_dynamic_array2![T, [self.k, self.k]]; for i in 0..self.k { sigma_mat[[i, i]] = T::from(sigma[i]).unwrap() } @@ -140,7 +134,7 @@ where let s_block = unsafe { rlst_pointer_mat!['static, T, s_block.as_ptr(), (nst, self.k), (1, nst)] }; - let mut c = rlst_dynamic_mat![T, (self.k, self.k * ntransfer_vectors)]; + let mut c = rlst_dynamic_array2![T, [self.k, self.k * ntransfer_vectors]]; for i in 0..self.transfer_vectors.len() { let top_left = (0, i * ncols); @@ -549,7 +543,7 @@ mod test { use bempp_kernel::laplace_3d::Laplace3dKernel; use cauchy::{c32, c64}; use num::complex::Complex; - use rlst::dense::RandomAccessMut; + use rlst_common::traits::RandomAccessMut; #[test] pub fn test_svd_operator_data() { @@ -636,7 +630,7 @@ mod test { // Some expansion data let ncoeffs = 6 * (order - 1).pow(2) + 2; - let mut multipole = rlst_dynamic_mat![f64, (ncoeffs, 1)]; + let mut multipole = rlst_dynamic_array2![f64, [ncoeffs, 1]]; for i in 0..ncoeffs { *multipole.get_mut(i, 0).unwrap() = i as f64; @@ -721,7 +715,7 @@ mod test { // Some expansion data998 let ncoeffs = 6 * (order - 1).pow(2) + 2; - let mut multipole = rlst_dynamic_mat![f64, (ncoeffs, 1)]; + let mut multipole = rlst_dynamic_array2![f64, [ncoeffs, 1]]; for i in 0..ncoeffs { *multipole.get_mut(i, 0).unwrap() = i as f64; @@ -885,7 +879,7 @@ mod test { // Some expansion data let ncoeffs = 6 * (order - 1).pow(2) + 2; - let mut multipole = rlst_dynamic_mat![f64, (ncoeffs, 1)]; + let mut multipole = rlst_dynamic_array2![f64, [ncoeffs, 1]]; for i in 0..ncoeffs { *multipole.get_mut(i, 0).unwrap() = i as f64; diff --git a/field/src/types.rs b/field/src/types.rs index 2723b5a1..9aec7bbb 100644 --- a/field/src/types.rs +++ b/field/src/types.rs @@ -3,12 +3,9 @@ use std::collections::HashMap; use cauchy::{c32, c64}; use num::{Complex, Float}; -use rlst::{ - common::traits::{Eval, NewLikeSelf}, - dense::{ - base_matrix::BaseMatrix, data_container::VectorContainer, matrix::Matrix, rlst_dynamic_mat, - Dynamic, - }, +use rlst_common::traits::{Eval, NewLikeSelf}; +use rlst_dense::{ + base_array::BaseArray, data_container::VectorContainer, array::Array, rlst_dynamic_array2, }; use bempp_traits::kernel::Kernel; @@ -16,7 +13,7 @@ use bempp_traits::types::Scalar; use bempp_tree::types::morton::MortonKey; /// Simple type alias for a 2D `Matrix` -pub type SvdM2lEntry = Matrix, Dynamic>, Dynamic>; +pub type SvdM2lEntry = Array, 2>, 2>; /// Simple type alias for pre-computed FFT of green's function evaluations computed for each transfer vector in a box's halo /// Each index corresponds to a halo position, and contains 64 convolutions, one for each of a box's siblings with each child @@ -117,7 +114,7 @@ where T: Scalar, { fn default() -> Self { - let tmp = rlst_dynamic_mat![T, (1, 1)]; + let tmp = rlst_dynamic_array2![T, [1, 1]]; SvdM2lOperatorData { u: tmp.new_like_self().eval(), @@ -128,16 +125,16 @@ where } /// Type alias for real coefficients for into FFTW wrappers -pub type FftMatrixf64 = Matrix, Dynamic>, Dynamic>; +pub type FftMatrixf64 = Array, 2>, 2>; /// Type alias for real coefficients for into FFTW wrappers -pub type FftMatrixf32 = Matrix, Dynamic>, Dynamic>; +pub type FftMatrixf32 = Array, 2>, 2>; /// Type alias for complex coefficients for FFTW wrappers -pub type FftMatrixc64 = Matrix, Dynamic>, Dynamic>; +pub type FftMatrixc64 = Array, 2>, 2>; /// Type alias for complex coefficients for FFTW wrappers -pub type FftMatrixc32 = Matrix, Dynamic>, Dynamic>; +pub type FftMatrixc32 = Array, 2>, 2>; /// Type alias for real coefficients for into FFTW wrappers -pub type FftMatrix = Matrix, Dynamic>, Dynamic>; +pub type FftMatrix = Array, 2>, 2>; diff --git a/fmm/Cargo.toml b/fmm/Cargo.toml index 60a5e9f7..1e8951ba 100644 --- a/fmm/Cargo.toml +++ b/fmm/Cargo.toml @@ -35,6 +35,6 @@ rand = "0.8.*" float-cmp = "0.9.0" num_cpus = "1" num = "0.4" -rlst = { git = "https://github.com/linalg-rs/rlst.git" } +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } fftw = {git = "https://github.com/skailasa/fftw.git" } rayon = "1.7" diff --git a/grid/Cargo.toml b/grid/Cargo.toml index f540f853..15fdf529 100644 --- a/grid/Cargo.toml +++ b/grid/Cargo.toml @@ -26,6 +26,6 @@ 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-proc-macro = { git = "https://github.com/linalg-rs/rlst.git" } -rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-proc-macro = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } diff --git a/grid/examples/curved_cells.rs b/grid/examples/curved_cells.rs index 85e2df25..7d7b4e9a 100644 --- a/grid/examples/curved_cells.rs +++ b/grid/examples/curved_cells.rs @@ -14,7 +14,7 @@ fn main() { 0.25, 0.5, 0.5, 0.5, 0.5, 0.75, 1.0, 1.0, 1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ], - (13, 3), + [13, 3], ), AdjacencyList::from_data( vec![2, 7, 12, 0, 2, 9, 11, 1, 4, 6, 10, 5, 2, 7, 11, 8, 6, 3], diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 69bdcde6..cadc733a 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -7,11 +7,10 @@ use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; use bempp_traits::grid::{Geometry, GeometryEvaluator, Grid, Ownership, Topology}; use itertools::izip; -use rlst_dense::{ - rlst_static_mat, RandomAccessByRef, RandomAccessMut, Shape, SizeIdentifier, +use rlst_common::traits::{RandomAccessByRef, RandomAccessMut, Shape, UnsafeRandomAccessByRef, UnsafeRandomAccessMut, }; -use rlst_proc_macro::rlst_static_size; +use rlst_proc_macro::rlst_static_array; use std::cell::RefCell; use std::ptr; @@ -30,13 +29,13 @@ impl<'a> EvaluatorTdim2Gdim3<'a> { element: &impl FiniteElement, points: &'a Mat, ) -> Self { - let npts = points.shape().0; - assert_eq!(points.shape().1, 2); + let npts = points.shape()[0]; + assert_eq!(points.shape()[1], 2); assert_eq!(geometry.dim(), 3); let mut table = Array4D::::new(element.tabulate_array_shape(1, npts)); element.tabulate(points, 1, &mut table); - let axes = RefCell::new(zero_matrix((2, 3))); - let js = RefCell::new(zero_matrix((npts, 6))); + let axes = RefCell::new(zero_matrix([2, 3])); + let js = RefCell::new(zero_matrix([npts, 6])); Self { geometry, @@ -58,7 +57,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { for i in 0..3 { for p in 0..self.npts { unsafe { - *points.get_unchecked_mut(p, i) = 0.0; + *points.get_unchecked_mut([p, i]) = 0.0; } } } @@ -67,7 +66,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { for j in 0..3 { for p in 0..self.npts { unsafe { - *points.get_unchecked_mut(p, j) += + *points.get_unchecked_mut([p, j]) += *self.geometry.coordinate_unchecked(v, j) * *self.table.get(0, p, i, 0).unwrap(); } @@ -82,7 +81,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { for i in 0..2 { for j in 0..3 { unsafe { - *axes.get_unchecked_mut(i, j) = 0.0; + *axes.get_unchecked_mut([i, j]) = 0.0; } } } @@ -90,30 +89,30 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { let v = unsafe { *self.geometry.cells.get_unchecked(cell_index, i) }; for j in 0..3 { unsafe { - *axes.get_unchecked_mut(0, j) += *self.geometry.coordinate_unchecked(v, j) + *axes.get_unchecked_mut([0, j]) += *self.geometry.coordinate_unchecked(v, j) * self.table.get(1, p, i, 0).unwrap(); - *axes.get_unchecked_mut(1, j) += *self.geometry.coordinate_unchecked(v, j) + *axes.get_unchecked_mut([1, j]) += *self.geometry.coordinate_unchecked(v, j) * self.table.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 = (*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)) + *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 = (*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(); - *normals.get_unchecked_mut(p, 0) /= size; - *normals.get_unchecked_mut(p, 1) /= size; - *normals.get_unchecked_mut(p, 2) /= size; + *normals.get_unchecked_mut([p, 0]) /= size; + *normals.get_unchecked_mut([p, 1]) /= size; + *normals.get_unchecked_mut([p, 2]) /= size; } } } @@ -122,7 +121,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { for i in 0..6 { for p in 0..self.npts { unsafe { - *jacobians.get_unchecked_mut(p, i) = 0.0; + *jacobians.get_unchecked_mut([p, i]) = 0.0; } } } @@ -132,7 +131,7 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { for k in 0..2 { for p in 0..self.npts { unsafe { - *jacobians.get_unchecked_mut(p, k + 2 * j) += + *jacobians.get_unchecked_mut([p, k + 2 * j]) += *self.geometry.coordinate_unchecked(v, j) * self.table.get(k + 1, p, i, 0).unwrap(); } @@ -147,15 +146,15 @@ impl<'a> GeometryEvaluator, Mat> for EvaluatorTdim2Gdim3<'a> { self.compute_jacobians(cell_index, &mut js); for (p, jdet) in jdets.iter_mut().enumerate() { unsafe { - *jdet = ((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)) + *jdet = ((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(); } @@ -181,12 +180,12 @@ impl<'a> LinearSimplexEvaluatorTdim2Gdim3<'a> { element: &impl FiniteElement, points: &'a Mat, ) -> Self { - let npts = points.shape().0; - assert_eq!(points.shape().1, 2); + let npts = points.shape()[0]; + assert_eq!(points.shape()[1], 2); assert_eq!(geometry.dim(), 3); let mut table = Array4D::::new(element.tabulate_array_shape(1, npts)); element.tabulate(points, 1, &mut table); - let axes = RefCell::new(zero_matrix((2, 3))); + let axes = RefCell::new(zero_matrix([2, 3])); let js = RefCell::new(vec![0.0; 6]); Self { @@ -232,7 +231,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd for i in 0..3 { for p in 0..self.npts { unsafe { - *points.get_unchecked_mut(p, i) = 0.0; + *points.get_unchecked_mut([p, i]) = 0.0; } } } @@ -241,7 +240,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd for j in 0..3 { for p in 0..self.npts { unsafe { - *points.get_unchecked_mut(p, j) += + *points.get_unchecked_mut([p, j]) += *self.geometry.coordinate_unchecked(v, j) * *self.table.get(0, p, i, 0).unwrap(); } @@ -255,7 +254,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd for j in 0..3 { for i in 0..2 { unsafe { - *axes.get_unchecked_mut(i, j) = 0.0; + *axes.get_unchecked_mut([i, j]) = 0.0; } } } @@ -264,34 +263,34 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd for j in 0..3 { for k in 0..2 { unsafe { - *axes.get_unchecked_mut(k, j) += *self.geometry.coordinate_unchecked(v, j) + *axes.get_unchecked_mut([k, j]) += *self.geometry.coordinate_unchecked(v, j) * self.table.get(k + 1, 0, i, 0).unwrap(); } } } } unsafe { - *normals.get_unchecked_mut(0, 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(0, 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(0, 2) = *axes.get_unchecked(0, 0) - * *axes.get_unchecked(1, 1) - - *axes.get_unchecked(0, 1) * *axes.get_unchecked(1, 0); - let size = ((*normals.get_unchecked(0, 0)).powi(2) - + (*normals.get_unchecked(0, 1)).powi(2) - + (*normals.get_unchecked(0, 2)).powi(2)) + *normals.get_unchecked_mut([0, 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([0, 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([0, 2]) = *axes.get_unchecked([0, 0]) + * *axes.get_unchecked([1, 1]) + - *axes.get_unchecked([0, 1]) * *axes.get_unchecked([1, 0]); + let size = ((*normals.get_unchecked([0, 0])).powi(2) + + (*normals.get_unchecked([0, 1])).powi(2) + + (*normals.get_unchecked([0, 2])).powi(2)) .sqrt(); - *normals.get_unchecked_mut(0, 0) /= size; - *normals.get_unchecked_mut(0, 1) /= size; - *normals.get_unchecked_mut(0, 2) /= size; + *normals.get_unchecked_mut([0, 0]) /= size; + *normals.get_unchecked_mut([0, 1]) /= size; + *normals.get_unchecked_mut([0, 2]) /= size; } for i in 0..3 { for p in 1..self.npts { unsafe { - *normals.get_unchecked_mut(p, i) = *normals.get_unchecked(0, i); + *normals.get_unchecked_mut([p, i]) = *normals.get_unchecked([0, i]); } } } @@ -303,7 +302,7 @@ impl<'a> GeometryEvaluator, Mat> for LinearSimplexEvaluatorTdim2Gd for i in 0..6 { for p in 0..self.npts { unsafe { - *jacobians.get_unchecked_mut(p, i) = *js.get_unchecked(i); + *jacobians.get_unchecked_mut([p, i]) = *js.get_unchecked(i); } } } @@ -344,9 +343,6 @@ pub struct SerialGeometry { index_map: Vec, } -#[rlst_static_size(2, 3)] -struct TwoByThree; - fn element_from_npts(cell_type: ReferenceCellType, npts: usize) -> CiarletElement { create_element( ElementFamily::Lagrange, @@ -420,14 +416,14 @@ impl SerialGeometry { impl SerialGeometry { unsafe fn coordinate_unchecked(&self, point_index: usize, coord_index: usize) -> &f64 { - self.coordinates.get_unchecked(point_index, coord_index) + self.coordinates.get_unchecked([point_index, coord_index]) } } impl Geometry for SerialGeometry { type T = Mat; type TMut = Mat; fn dim(&self) -> usize { - self.coordinates.shape().1 + self.coordinates.shape()[1] } fn coordinate(&self, point_index: usize, coord_index: usize) -> Option<&f64> { @@ -439,7 +435,7 @@ impl Geometry for SerialGeometry { } fn point_count(&self) -> usize { - self.coordinates.shape().0 + self.coordinates.shape()[0] } fn cell_vertices(&self, index: usize) -> Option<&[usize]> { @@ -465,20 +461,20 @@ impl Geometry for SerialGeometry { } fn compute_points< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, cell: usize, physical_points: &mut TMut, ) { - let npts = points.shape().0; + let npts = points.shape()[0]; let gdim = self.dim(); - if physical_points.shape().0 != npts { + if physical_points.shape()[0] != npts { panic!("physical_points has wrong number of rows."); } - if physical_points.shape().1 != gdim { + if physical_points.shape()[1] != gdim { panic!("physical_points has wrong number of columns."); } let element = self.element(cell); @@ -486,90 +482,90 @@ impl Geometry for SerialGeometry { element.tabulate(points, 0, &mut data); for p in 0..npts { for i in 0..gdim { - *physical_points.get_mut(p, i).unwrap() = 0.0; + *physical_points.get_mut([p, i]).unwrap() = 0.0; } } for i in 0..data.shape().2 { let v = *self.cells.get(cell, i).unwrap(); for j in 0..gdim { for p in 0..npts { - *physical_points.get_mut(p, j).unwrap() += + *physical_points.get_mut([p, j]).unwrap() += *self.coordinate(v, j).unwrap() * data.get(0, p, i, 0).unwrap(); } } } } fn compute_normals< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, cell: usize, normals: &mut TMut, ) { - let npts = points.shape().0; - let tdim = points.shape().1; + let npts = points.shape()[0]; + let tdim = points.shape()[1]; let gdim = self.dim(); if gdim != 3 { unimplemented!("normals currently only implemented for 2D cells embedded in 3D."); } - if normals.shape().0 != npts { + if normals.shape()[0] != npts { panic!("normals has wrong number of columns."); } - if normals.shape().1 != gdim { + if normals.shape()[1] != gdim { panic!("normals has wrong number of rows."); } let element = self.element(cell); let mut data = Array4D::::new(element.tabulate_array_shape(1, npts)); - let mut axes = rlst_static_mat![f64, TwoByThree]; + let mut axes = rlst_static_array![f64, 2, 3]; element.tabulate(points, 1, &mut data); for p in 0..npts { for i in 0..tdim { for j in 0..gdim { - *axes.get_mut(i, j).unwrap() = 0.0; + *axes.get_mut([i, j]).unwrap() = 0.0; } } for i in 0..data.shape().2 { let v = *self.cells.get(cell, i).unwrap(); for j in 0..gdim { - *axes.get_mut(0, j).unwrap() += + *axes.get_mut([0, j]).unwrap() += *self.coordinate(v, j).unwrap() * data.get(1, p, i, 0).unwrap(); - *axes.get_mut(1, j).unwrap() += + *axes.get_mut([1, j]).unwrap() += *self.coordinate(v, j).unwrap() * data.get(2, p, i, 0).unwrap(); } } - *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()) + *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; + *normals.get_mut([p, 0]).unwrap() /= size; + *normals.get_mut([p, 1]).unwrap() /= size; + *normals.get_mut([p, 2]).unwrap() /= size; } } fn compute_jacobians< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, cell: usize, jacobians: &mut TMut, ) { - let npts = points.shape().0; - let tdim = points.shape().1; + let npts = points.shape()[0]; + let tdim = points.shape()[1]; let gdim = self.dim(); - if jacobians.shape().0 != npts { + if jacobians.shape()[0] != npts { panic!("jacobians has wrong number of rows."); } - if jacobians.shape().1 != gdim * tdim { + if jacobians.shape()[1] != gdim * tdim { panic!("jacobians has wrong number of columns."); } let element = self.element(cell); @@ -577,8 +573,8 @@ impl Geometry for SerialGeometry { let tdim = data.shape().0 - 1; element.tabulate(points, 1, &mut data); for p in 0..npts { - for i in 0..jacobians.shape().1 { - *jacobians.get_mut(p, i).unwrap() = 0.0; + for i in 0..jacobians.shape()[1] { + *jacobians.get_mut([p, i]).unwrap() = 0.0; } } for i in 0..data.shape().2 { @@ -586,38 +582,38 @@ impl Geometry for SerialGeometry { for p in 0..npts { for j in 0..gdim { for k in 0..tdim { - *jacobians.get_mut(p, k + tdim * j).unwrap() += + *jacobians.get_mut([p, k + tdim * j]).unwrap() += *self.coordinate(v, j).unwrap() * data.get(k + 1, p, i, 0).unwrap(); } } } } } - fn compute_jacobian_determinants + Shape>( + fn compute_jacobian_determinants + Shape<2>>( &self, points: &T, cell: usize, jacobian_determinants: &mut [f64], ) { - let npts = points.shape().0; - let tdim = points.shape().1; + let npts = points.shape()[0]; + let tdim = points.shape()[1]; let gdim = self.dim(); - if points.shape().0 != jacobian_determinants.len() { + if points.shape()[0] != jacobian_determinants.len() { panic!("jacobian_determinants has wrong length."); } - let mut js = zero_matrix((npts, gdim * tdim)); + let mut js = zero_matrix([npts, gdim * tdim]); self.compute_jacobians(points, cell, &mut js); for (p, jdet) in jacobian_determinants.iter_mut().enumerate() { *jdet = match tdim { 1 => match gdim { - 1 => *js.get(p, 0).unwrap(), + 1 => *js.get([p, 0]).unwrap(), 2 => { - ((*js.get(p, 0).unwrap()).powi(2) + (*js.get(p, 1).unwrap()).powi(2)).sqrt() + ((*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)) + 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."); @@ -625,18 +621,18 @@ impl Geometry for SerialGeometry { }, 2 => match gdim { 2 => { - *js.get(p, 0).unwrap() * *js.get(p, 3).unwrap() - - *js.get(p, 1).unwrap() * *js.get(p, 2).unwrap() + *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()) + 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(), _ => { @@ -645,15 +641,15 @@ impl Geometry for SerialGeometry { }, 3 => match gdim { 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()) + *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."); @@ -666,21 +662,21 @@ impl Geometry for SerialGeometry { } } fn compute_jacobian_inverses< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, cell: usize, jacobian_inverses: &mut TMut, ) { - let npts = points.shape().0; + let npts = points.shape()[0]; let gdim = self.dim(); - let tdim = points.shape().1; - if jacobian_inverses.shape().0 != npts { + let tdim = points.shape()[1]; + if jacobian_inverses.shape()[0] != npts { panic!("jacobian_inverses has wrong number of rows."); } - if jacobian_inverses.shape().1 != gdim * tdim { + if jacobian_inverses.shape()[1] != gdim * tdim { panic!("jacobian_inverses has wrong number of columns."); } let element = self.element(cell); @@ -689,13 +685,13 @@ impl Geometry for SerialGeometry { && element.degree() == 1 { // Map is affine - let mut js = zero_matrix((npts, gdim * tdim)); + let mut js = zero_matrix([npts, gdim * tdim]); self.compute_jacobians(points, cell, &mut js); for p in 0..npts { if tdim == 1 { if gdim == 1 { - *jacobian_inverses.get_mut(p, 0).unwrap() = 1.0 / *js.get(p, 0).unwrap(); + *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 { @@ -705,63 +701,63 @@ impl Geometry for SerialGeometry { } } else if tdim == 2 { if gdim == 2 { - 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; + 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 = (*js.get(p, 3).unwrap() * *js.get(p, 4).unwrap() - - *js.get(p, 2).unwrap() * *js.get(p, 5).unwrap()) + 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()) + + (*js.get([p, 5]).unwrap() * *js.get([p, 0]).unwrap() + - *js.get([p, 4]).unwrap() * *js.get([p, 1]).unwrap()) .powi(2) - + (*js.get(p, 1).unwrap() * *js.get(p, 2).unwrap() - - *js.get(p, 0).unwrap() * *js.get(p, 3).unwrap()) + + (*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())) + *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())) + *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())) + *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())) + *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())) + *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())) + *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."); @@ -1136,7 +1132,7 @@ mod test { #[test] fn test_connectivity() { let g = SerialGrid::new( - to_matrix(&[0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0], (4, 2)), + to_matrix(&[0.0, 1.0, 0.0, 1.0, 0.0, 0.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], ); @@ -1242,7 +1238,7 @@ mod test { 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, ], - (6, 3), + [6, 3], ), AdjacencyList::from_data( vec![ @@ -1268,7 +1264,7 @@ mod test { 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, ], - (9, 2), + [9, 2], ), AdjacencyList::from_data( vec![ @@ -1294,7 +1290,7 @@ mod test { 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, ], - (9, 2), + [9, 2], ), AdjacencyList::from_data( vec![0, 1, 4, 0, 4, 3, 1, 2, 4, 5, 3, 4, 7, 3, 7, 6, 4, 5, 7, 8], @@ -1325,7 +1321,7 @@ mod test { &[ 0.0, 1.0, s, 0.0, -s, -1.0, -s, 0.0, s, 0.0, 0.0, s, 1.0, s, 0.0, -s, -1.0, -s, ], - (9, 2), + [9, 2], ), AdjacencyList::from_data(vec![4, 8, 2, 1, 3, 0, 4, 6, 8, 7, 0, 5], vec![0, 6, 12]), vec![ReferenceCellType::Triangle, ReferenceCellType::Triangle], @@ -1347,7 +1343,7 @@ mod test { 0.25, 0.5, 0.5, 0.5, 0.5, 0.75, 1.0, 1.0, 1.0, -0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ], - (13, 3), + [13, 3], ), AdjacencyList::from_data( vec![2, 7, 12, 0, 2, 9, 11, 1, 4, 6, 10, 5, 2, 7, 11, 8, 6, 3], @@ -1375,7 +1371,7 @@ mod test { 2.0, 4.0, 5.0, 0.0, -1.0, 0.0, 2.0, 2.0, 3.0, 0.0, 0.0, -1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, ], - (6, 3), + [6, 3], ), AdjacencyList::from_data(vec![0, 1, 2, 3, 4, 5], vec![0, 3, 6]), vec![ReferenceCellType::Triangle, ReferenceCellType::Triangle], @@ -1385,159 +1381,159 @@ mod test { let points = to_matrix( &[0.2, 0.5, 1.0 / 3.0, 0.15, 0.0, 0.5, 1.0 / 3.0, 0.3], - (4, 2), + [4, 2], ); // Test compute_points - let mut physical_points = zero_matrix((points.shape().0, 3)); + let mut physical_points = zero_matrix([points.shape()[0], 3]); g.geometry() .compute_points(&points, 0, &mut physical_points); assert_relative_eq!( - *physical_points.get(0, 0).unwrap(), + *physical_points.get([0, 0]).unwrap(), 2.4, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(0, 1).unwrap(), + *physical_points.get([0, 1]).unwrap(), 2.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(0, 2).unwrap(), + *physical_points.get([0, 2]).unwrap(), 0.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 0).unwrap(), + *physical_points.get([1, 0]).unwrap(), 4.5, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 1).unwrap(), + *physical_points.get([1, 1]).unwrap(), 2.5, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 2).unwrap(), + *physical_points.get([1, 2]).unwrap(), 0.5, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 0).unwrap(), + *physical_points.get([2, 0]).unwrap(), 11.0 / 3.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 1).unwrap(), + *physical_points.get([2, 1]).unwrap(), 7.0 / 3.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 2).unwrap(), + *physical_points.get([2, 2]).unwrap(), 1.0 / 3.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 0).unwrap(), + *physical_points.get([3, 0]).unwrap(), 3.2, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 1).unwrap(), + *physical_points.get([3, 1]).unwrap(), 2.3, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 2).unwrap(), + *physical_points.get([3, 2]).unwrap(), 0.3, max_relative = 1e-14 ); g.geometry() .compute_points(&points, 1, &mut physical_points); assert_relative_eq!( - *physical_points.get(0, 0).unwrap(), + *physical_points.get([0, 0]).unwrap(), -0.2, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(0, 1).unwrap(), + *physical_points.get([0, 1]).unwrap(), 0.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(0, 2).unwrap(), + *physical_points.get([0, 2]).unwrap(), 1.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 0).unwrap(), + *physical_points.get([1, 0]).unwrap(), -0.5, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 1).unwrap(), + *physical_points.get([1, 1]).unwrap(), -0.5, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(1, 2).unwrap(), + *physical_points.get([1, 2]).unwrap(), 1.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 0).unwrap(), + *physical_points.get([2, 0]).unwrap(), -1.0 / 3.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 1).unwrap(), + *physical_points.get([2, 1]).unwrap(), -1.0 / 3.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(2, 2).unwrap(), + *physical_points.get([2, 2]).unwrap(), 1.0, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 0).unwrap(), + *physical_points.get([3, 0]).unwrap(), -0.15, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 1).unwrap(), + *physical_points.get([3, 1]).unwrap(), -0.3, max_relative = 1e-14 ); assert_relative_eq!( - *physical_points.get(3, 2).unwrap(), + *physical_points.get([3, 2]).unwrap(), 1.0, max_relative = 1e-14 ); // Test compute_jacobians - let mut jacobians = zero_matrix((points.shape().0, 6)); + let mut jacobians = zero_matrix([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); - assert_relative_eq!(*jacobians.get(i, 1).unwrap(), 3.0, max_relative = 1e-14); - // assert_relative_eq!(*jacobians.get(i, 2).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 3).unwrap(), 1.0, max_relative = 1e-14); - // assert_relative_eq!(*jacobians.get(i, 4).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 5).unwrap(), 1.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 0]).unwrap(), 2.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 1]).unwrap(), 3.0, max_relative = 1e-14); + // assert_relative_eq!(*jacobians.get([i, 2]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 3]).unwrap(), 1.0, max_relative = 1e-14); + // assert_relative_eq!(*jacobians.get([i, 4]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 5]).unwrap(), 1.0, max_relative = 1e-14); } g.geometry().compute_jacobians(&points, 1, &mut jacobians); for i in 0..3 { - assert_relative_eq!(*jacobians.get(i, 0).unwrap(), -1.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 1).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 2).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 3).unwrap(), -1.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 4).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jacobians.get(i, 5).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 0]).unwrap(), -1.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 1]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 2]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 3]).unwrap(), -1.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 4]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jacobians.get([i, 5]).unwrap(), 0.0, max_relative = 1e-14); } // test compute_jacobian_determinants - let mut dets = vec![0.0; points.shape().0]; + let mut dets = vec![0.0; points.shape()[0]]; g.geometry() .compute_jacobian_determinants(&points, 0, &mut dets); for d in &dets { @@ -1550,26 +1546,26 @@ mod test { } // Test compute_jacobian_inverses - let mut jinvs = zero_matrix((points.shape().0, 6)); + let mut jinvs = zero_matrix([points.shape()[0], 6]); g.geometry() .compute_jacobian_inverses(&points, 0, &mut jinvs); for i in 0..3 { - assert_relative_eq!(*jinvs.get(i, 0).unwrap(), 0.5, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 1).unwrap(), -0.75, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 2).unwrap(), -0.75, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 3).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 4).unwrap(), 0.5, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 5).unwrap(), 0.5, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 0]).unwrap(), 0.5, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 1]).unwrap(), -0.75, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 2]).unwrap(), -0.75, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 3]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 4]).unwrap(), 0.5, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 5]).unwrap(), 0.5, max_relative = 1e-14); } g.geometry() .compute_jacobian_inverses(&points, 1, &mut jinvs); for i in 0..3 { - assert_relative_eq!(*jinvs.get(i, 0).unwrap(), -1.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 1).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 2).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 3).unwrap(), 0.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 4).unwrap(), -1.0, max_relative = 1e-14); - assert_relative_eq!(*jinvs.get(i, 5).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 0]).unwrap(), -1.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 1]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 2]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 3]).unwrap(), 0.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 4]).unwrap(), -1.0, max_relative = 1e-14); + assert_relative_eq!(*jinvs.get([i, 5]).unwrap(), 0.0, max_relative = 1e-14); } } @@ -1580,7 +1576,7 @@ mod test { &[ 0.0, 1.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, ], - (5, 3), + [5, 3], ), AdjacencyList::from_data(vec![0, 1, 2, 1, 3, 2, 2, 3, 4], vec![0, 3, 6, 9]), vec![ @@ -1592,25 +1588,25 @@ mod test { ], ); - let pt = to_matrix(&[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 = zero_matrix((1, 3)); + let mut normal = zero_matrix([1, 3]); g.geometry().compute_normals(&pt, 0, &mut normal); - assert_relative_eq!(*normal.get(0, 0).unwrap(), 0.0); - assert_relative_eq!(*normal.get(0, 1).unwrap(), -1.0); - assert_relative_eq!(*normal.get(0, 2).unwrap(), 0.0); + assert_relative_eq!(*normal.get([0, 0]).unwrap(), 0.0); + assert_relative_eq!(*normal.get([0, 1]).unwrap(), -1.0); + assert_relative_eq!(*normal.get([0, 2]).unwrap(), 0.0); g.geometry().compute_normals(&pt, 1, &mut normal); let a = f64::sqrt(1.0 / 3.0); - assert_relative_eq!(*normal.get(0, 0).unwrap(), a); - assert_relative_eq!(*normal.get(0, 1).unwrap(), a); - assert_relative_eq!(*normal.get(0, 2).unwrap(), a); + assert_relative_eq!(*normal.get([0, 0]).unwrap(), a); + assert_relative_eq!(*normal.get([0, 1]).unwrap(), a); + assert_relative_eq!(*normal.get([0, 2]).unwrap(), a); g.geometry().compute_normals(&pt, 2, &mut normal); - assert_relative_eq!(*normal.get(0, 0).unwrap(), 0.0); - assert_relative_eq!(*normal.get(0, 1).unwrap(), 0.0); - assert_relative_eq!(*normal.get(0, 2).unwrap(), 1.0); + assert_relative_eq!(*normal.get([0, 0]).unwrap(), 0.0); + assert_relative_eq!(*normal.get([0, 1]).unwrap(), 0.0); + assert_relative_eq!(*normal.get([0, 2]).unwrap(), 1.0); // Test a curved quadrilateral cell let curved_g = SerialGrid::new( @@ -1619,67 +1615,67 @@ mod test { -1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, -1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, ], - (9, 3), + [9, 3], ), AdjacencyList::from_data(vec![0, 1, 2, 3, 4, 5, 6, 7, 8], vec![0, 9]), vec![ReferenceCellType::Quadrilateral], ); - let points = to_matrix(&[0.0, 0.2, 0.5, 0.7, 1.0, 0.0, 0.3, 0.9, 1.0, 0.3], (5, 2)); - let mut normals = zero_matrix((5, 3)); + let points = to_matrix(&[0.0, 0.2, 0.5, 0.7, 1.0, 0.0, 0.3, 0.9, 1.0, 0.3], [5, 2]); + let mut normals = zero_matrix([5, 3]); curved_g .geometry() .compute_normals(&points, 0, &mut normals); assert_relative_eq!( - *normals.get(0, 0).unwrap(), + *normals.get([0, 0]).unwrap(), 2.0 * f64::sqrt(1.0 / 5.0), epsilon = 1e-12 ); - assert_relative_eq!(*normals.get(0, 1).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([0, 1]).unwrap(), 0.0, epsilon = 1e-12); assert_relative_eq!( - *normals.get(0, 2).unwrap(), + *normals.get([0, 2]).unwrap(), f64::sqrt(1.0 / 5.0), epsilon = 1e-12 ); assert_relative_eq!( - *normals.get(1, 0).unwrap(), + *normals.get([1, 0]).unwrap(), 1.2 * f64::sqrt(1.0 / 2.44), epsilon = 1e-12 ); - assert_relative_eq!(*normals.get(1, 1).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([1, 1]).unwrap(), 0.0, epsilon = 1e-12); assert_relative_eq!( - *normals.get(1, 2).unwrap(), + *normals.get([1, 2]).unwrap(), f64::sqrt(1.0 / 2.44), epsilon = 1e-12 ); - assert_relative_eq!(*normals.get(2, 0).unwrap(), 0.0, epsilon = 1e-12); - assert_relative_eq!(*normals.get(2, 1).unwrap(), 0.0, epsilon = 1e-12); - assert_relative_eq!(*normals.get(2, 2).unwrap(), 1.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([2, 0]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([2, 1]).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([2, 2]).unwrap(), 1.0, epsilon = 1e-12); assert_relative_eq!( - *normals.get(3, 0).unwrap(), + *normals.get([3, 0]).unwrap(), -0.8 * f64::sqrt(1.0 / 1.64), epsilon = 1e-12 ); - assert_relative_eq!(*normals.get(3, 1).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([3, 1]).unwrap(), 0.0, epsilon = 1e-12); assert_relative_eq!( - *normals.get(3, 2).unwrap(), + *normals.get([3, 2]).unwrap(), f64::sqrt(1.0 / 1.64), epsilon = 1e-12 ); assert_relative_eq!( - *normals.get(4, 0).unwrap(), + *normals.get([4, 0]).unwrap(), -2.0 * f64::sqrt(1.0 / 5.0), epsilon = 1e-12 ); - assert_relative_eq!(*normals.get(4, 1).unwrap(), 0.0, epsilon = 1e-12); + assert_relative_eq!(*normals.get([4, 1]).unwrap(), 0.0, epsilon = 1e-12); assert_relative_eq!( - *normals.get(4, 2).unwrap(), + *normals.get([4, 2]).unwrap(), f64::sqrt(1.0 / 5.0), epsilon = 1e-12 ); @@ -1694,19 +1690,19 @@ mod test { 1, Continuity::Continuous, ); - let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], (3, 2)); + let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); - let mut points0 = zero_matrix((3, 3)); - let mut points1 = zero_matrix((3, 3)); + let mut points0 = zero_matrix([3, 3]); + let mut points1 = zero_matrix([3, 3]); for c in 0..grid.geometry().cell_count() { grid.geometry().compute_points(&pts, c, &mut points0); e.compute_points(c, &mut points1); for i in 0..3 { for j in 0..3 { assert_relative_eq!( - *points0.get(i, j).unwrap(), - *points1.get(i, j).unwrap(), + *points0.get([i, j]).unwrap(), + *points1.get([i, j]).unwrap(), epsilon = 1e-12 ); } @@ -1723,19 +1719,19 @@ mod test { 1, Continuity::Continuous, ); - let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], (3, 2)); + let pts = to_matrix(&[0.1, 0.2, 0.6, 0.1, 0.4, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); - let mut normals0 = zero_matrix((3, 3)); - let mut normals1 = zero_matrix((3, 3)); + let mut normals0 = zero_matrix([3, 3]); + let mut normals1 = zero_matrix([3, 3]); for c in 0..grid.geometry().cell_count() { grid.geometry().compute_normals(&pts, c, &mut normals0); e.compute_normals(c, &mut normals1); for i in 0..3 { for j in 0..3 { assert_relative_eq!( - *normals0.get(i, j).unwrap(), - *normals1.get(i, j).unwrap(), + *normals0.get([i, j]).unwrap(), + *normals1.get([i, j]).unwrap(), epsilon = 1e-12 ); } @@ -1752,19 +1748,19 @@ mod test { 1, Continuity::Continuous, ); - let pts = to_matrix(&[0.1, 0.1, 0.2, 0.4, 0.6, 0.2], (3, 2)); + let pts = to_matrix(&[0.1, 0.1, 0.2, 0.4, 0.6, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); - let mut jacobians0 = zero_matrix((3, 6)); - let mut jacobians1 = zero_matrix((3, 6)); + let mut jacobians0 = zero_matrix([3, 6]); + let mut jacobians1 = zero_matrix([3, 6]); for c in 0..grid.geometry().cell_count() { grid.geometry().compute_jacobians(&pts, c, &mut jacobians0); e.compute_jacobians(c, &mut jacobians1); for i in 0..3 { for j in 0..6 { assert_relative_eq!( - *jacobians0.get(i, j).unwrap(), - *jacobians1.get(i, j).unwrap(), + *jacobians0.get([i, j]).unwrap(), + *jacobians1.get([i, j]).unwrap(), epsilon = 1e-12 ); } @@ -1781,7 +1777,7 @@ mod test { 1, Continuity::Continuous, ); - let pts = to_matrix(&[0.1, 0.1, 0.2, 0.4, 0.6, 0.2], (3, 2)); + let pts = to_matrix(&[0.1, 0.1, 0.2, 0.4, 0.6, 0.2], [3, 2]); let e = grid.geometry().get_evaluator(&element, &pts); let mut jacobian_determinants0 = vec![0.0; 3]; diff --git a/grid/src/io.rs b/grid/src/io.rs index 7014e12e..ab4e3c02 100644 --- a/grid/src/io.rs +++ b/grid/src/io.rs @@ -139,7 +139,7 @@ mod test { 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, ], - (9, 2), + [9, 2], ), AdjacencyList::from_data( vec![0, 1, 3, 4, 3, 4, 6, 7, 1, 2, 4, 5, 4, 5, 7, 8], @@ -158,7 +158,7 @@ mod test { 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, ], - (9, 2), + [9, 2], ), AdjacencyList::from_data( vec![0, 1, 4, 0, 4, 3, 1, 2, 4, 5, 3, 4, 7, 3, 7, 6, 4, 5, 7, 8], diff --git a/grid/src/shapes.rs b/grid/src/shapes.rs index 3905b2f5..dcd19ab6 100644 --- a/grid/src/shapes.rs +++ b/grid/src/shapes.rs @@ -2,11 +2,11 @@ use crate::grid::SerialGrid; use bempp_element::cell::Triangle; -use bempp_tools::arrays::{to_matrix, AdjacencyList}; +use bempp_tools::arrays::{zero_matrix, 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}; +use rlst_common::traits::{RandomAccessByRef, RandomAccessMut}; /// Create a regular sphere /// @@ -20,7 +20,7 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, ], - (6, 3), + [6, 3], ), AdjacencyList::from_data( vec![ @@ -35,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 = to_matrix(&vec![0.0; nvertices * 3], (nvertices, 3)); + let mut coordinates = zero_matrix([nvertices, 3]); let mut cells = AdjacencyList::::new(); for i in 0..ncells_old { @@ -48,7 +48,7 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { let gv = g.geometry().cell_vertices(gi).unwrap(); for (j, gv_j) in gv.iter().enumerate() { for k in 0..g.geometry().dim() { - *coordinates.get_mut(tv[j], k).unwrap() = + *coordinates.get_mut([tv[j], k]).unwrap() = *g.geometry().coordinate(*gv_j, k).unwrap(); } } @@ -57,15 +57,15 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { let vs = ref_e.connectivity(1, j, 0).unwrap(); let pt = (0..3) .map(|k| { - (*coordinates.get(tv[vs[0]], k).unwrap() - + *coordinates.get(tv[vs[1]], k).unwrap()) + (*coordinates.get([tv[vs[0]], k]).unwrap() + + *coordinates.get([tv[vs[1]], k]).unwrap()) / 2.0 }) .collect::>(); let norm = (pt[0].powi(2) + pt[1].powi(2) + pt[2].powi(2)).sqrt(); for (k, pt_k) in pt.iter().enumerate() { - *coordinates.get_mut(nvertices_old + tedges_j, k).unwrap() = *pt_k / norm; + *coordinates.get_mut([nvertices_old + tedges_j, k]).unwrap() = *pt_k / norm; } } @@ -92,7 +92,6 @@ pub fn regular_sphere(refinement_level: usize) -> SerialGrid { #[cfg(test)] mod test { use crate::shapes::*; - use bempp_tools::arrays::zero_matrix; #[test] fn test_regular_sphere_0() { @@ -110,17 +109,17 @@ mod test { fn test_normal_is_outward() { for i in 0..3 { let g = regular_sphere(i); - let points = to_matrix(&[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 = zero_matrix((1, 3)); - let mut normal = zero_matrix((1, 3)); + let mut mapped_pt = zero_matrix([1, 3]); + let mut normal = zero_matrix([1, 3]); for i in 0..g.geometry().cell_count() { g.geometry().compute_points(&points, i, &mut mapped_pt); g.geometry().compute_normals(&points, i, &mut normal); - let dot = *mapped_pt.get(0, 0).unwrap() * *normal.get(0, 0).unwrap() - + *mapped_pt.get(0, 1).unwrap() * *normal.get(0, 1).unwrap() - + *mapped_pt.get(0, 2).unwrap() * *normal.get(0, 2).unwrap(); + let dot = *mapped_pt.get([0, 0]).unwrap() * *normal.get([0, 0]).unwrap() + + *mapped_pt.get([0, 1]).unwrap() * *normal.get([0, 1]).unwrap() + + *mapped_pt.get([0, 2]).unwrap() * *normal.get([0, 2]).unwrap(); assert!(dot > 0.0); } } diff --git a/kernel/Cargo.toml b/kernel/Cargo.toml index c7f1f8a8..d06967b1 100644 --- a/kernel/Cargo.toml +++ b/kernel/Cargo.toml @@ -28,7 +28,9 @@ approx = "0.5" rayon = "1.7" num = "0.4" num_cpus = "1" -rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rand = "0.8.5" [dev-dependencies] -rlst = { git = "https://github.com/linalg-rs/rlst.git" } \ No newline at end of file +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } diff --git a/kernel/src/laplace_3d.rs b/kernel/src/laplace_3d.rs index bfc95c9f..ea7dd394 100644 --- a/kernel/src/laplace_3d.rs +++ b/kernel/src/laplace_3d.rs @@ -364,12 +364,42 @@ fn laplace_component_count(eval_type: EvalType) -> usize { mod test { use super::*; + use rand::prelude::*; use approx::assert_relative_eq; use bempp_traits::types::Scalar; - use rlst; - use rlst::common::traits::{Copy, Eval, Transpose}; - use rlst::dense::traits::*; - use rlst_dense; + use rlst_common::traits::{RawAccess, RawAccessMut, RandomAccessByRef, RandomAccessMut, Shape, FillFrom}; + use rlst_dense::rlst_dynamic_array2; + use bempp_tools::arrays::Mat; + + fn copy(m_in: &Mat) -> Mat { + let mut m = rlst_dynamic_array2!(f64, m_in.shape()); + for i in 0..m_in.shape()[0] { + for j in 0..m_in.shape()[1] { + *m.get_mut([i, j]).unwrap() = *m_in.get([i, j]).unwrap(); + } + } + m + } + + fn rand_mat(shape: [usize; 2]) -> Mat { + let mut m = rlst_dynamic_array2!(f64, shape); + let mut rng = rand::thread_rng(); + for i in 0..shape[0] { + for j in 0..shape[1] { + *m.get_mut([i, j]).unwrap() = rng.gen() + } + } + m + } + + fn rand_vec(size: usize) -> Mat { + let mut v = rlst_dynamic_array2!(f64, [size, 1]); + let mut rng = rand::thread_rng(); + for i in 0..size { + *v.get_mut([i, 0]).unwrap() = rng.gen(); + } + v + } #[test] fn test_laplace_3d() { @@ -378,10 +408,10 @@ mod test { let nsources = 5; let ntargets = 3; - let sources = rlst::dense::rlst_rand_mat![f64, (nsources, 3)]; - let targets = rlst::dense::rlst_rand_mat![f64, (ntargets, 3)]; - let charges = rlst::dense::rlst_rand_col_vec![f64, nsources]; - let mut green_value = rlst::dense::rlst_rand_col_vec![f64, ntargets]; + let sources = rand_mat([nsources, 3]); + let targets = rand_mat([ntargets, 3]); + let charges = rand_vec(nsources); + let mut green_value = rand_vec(ntargets); Laplace3dKernel::::default().evaluate_st( EvalType::Value, @@ -405,9 +435,9 @@ mod test { assert_relative_eq!(green_value[[target_index, 0]], expected, epsilon = 1E-12); } - let mut targets_x_eps = targets.copy(); - let mut targets_y_eps = targets.copy(); - let mut targets_z_eps = targets.copy(); + let mut targets_x_eps = copy(&targets); + let mut targets_y_eps = copy(&targets); + let mut targets_z_eps = copy(&targets); for index in 0..ntargets { targets_x_eps[[index, 0]] += eps; @@ -415,7 +445,7 @@ mod test { targets_z_eps[[index, 2]] += eps; } - let mut expected = rlst_dense::rlst_dynamic_mat!(f64, (4, ntargets)); + let mut expected = rlst_dynamic_array2!(f64, [4, ntargets]); Laplace3dKernel::::default().evaluate_st( EvalType::ValueDeriv, @@ -425,7 +455,7 @@ mod test { expected.data_mut(), ); - let mut green_value_x_eps = rlst::dense::rlst_col_vec![f64, ntargets]; + let mut green_value_x_eps = rlst_dynamic_array2![f64, [ntargets, 1]]; Laplace3dKernel::::default().evaluate_st( EvalType::Value, @@ -435,7 +465,7 @@ mod test { green_value_x_eps.data_mut(), ); - let mut green_value_y_eps = rlst::dense::rlst_col_vec![f64, ntargets]; + let mut green_value_y_eps = rlst_dynamic_array2![f64, [ntargets, 1]]; Laplace3dKernel::::default().evaluate_st( EvalType::Value, @@ -444,7 +474,7 @@ mod test { charges.data(), green_value_y_eps.data_mut(), ); - let mut green_value_z_eps = rlst::dense::rlst_col_vec![f64, ntargets]; + let mut green_value_z_eps = rlst_dynamic_array2![f64, [ntargets, 1]]; Laplace3dKernel::::default().evaluate_st( EvalType::Value, @@ -454,9 +484,16 @@ mod test { green_value_z_eps.data_mut(), ); - let x_deriv = ((green_value_x_eps - &green_value) * (1.0 / eps)).eval(); - let y_deriv = ((green_value_y_eps - &green_value) * (1.0 / eps)).eval(); - let z_deriv = ((green_value_z_eps - &green_value) * (1.0 / eps)).eval(); + let gv0 = copy(&green_value); + let gv1 = copy(&green_value); + let gv2 = copy(&green_value); + + let mut x_deriv = rlst_dynamic_array2![f64, [ntargets, 1]]; + let mut y_deriv = rlst_dynamic_array2![f64, [ntargets, 1]]; + let mut z_deriv = rlst_dynamic_array2![f64, [ntargets, 1]]; + x_deriv.fill_from((green_value_x_eps - gv0) * (1.0 / eps)); + y_deriv.fill_from((green_value_y_eps - gv1) * (1.0 / eps)); + z_deriv.fill_from((green_value_z_eps - gv2) * (1.0 / eps)); for target_index in 0..ntargets { assert_relative_eq!( @@ -489,25 +526,26 @@ mod test { let nsources = 3; let ntargets = 5; - let sources = rlst::dense::rlst_rand_mat![f64, (nsources, 3)]; - let targets = rlst::dense::rlst_rand_mat![f64, (ntargets, 3)]; - let mut green_value = rlst_dense::rlst_dynamic_mat!(f64, (nsources, ntargets)); + let sources = rand_mat([nsources, 3]); + let targets = rand_mat([ntargets, 3]); + let mut green_value_t = rlst_dynamic_array2!(f64, [nsources, ntargets]); Laplace3dKernel::::default().assemble_st( EvalType::Value, sources.data(), targets.data(), - green_value.data_mut(), + green_value_t.data_mut(), ); // The matrix needs to be transposed so that the first row corresponds to the first target, // second row to the second target and so on. - let green_value = green_value.transpose().eval(); + let mut green_value = rlst_dynamic_array2!(f64, [ntargets, nsources]); + green_value.fill_from(green_value_t.transpose()); for charge_index in 0..nsources { - let mut charges = rlst::dense::rlst_col_vec![f64, nsources]; - let mut expected = rlst::dense::rlst_col_vec![f64, ntargets]; + let mut charges = rlst_dynamic_array2![f64, [nsources, 1]]; + let mut expected = rlst_dynamic_array2![f64, [ntargets, 1]]; charges[[charge_index, 0]] = 1.0; Laplace3dKernel::::default().evaluate_st( @@ -527,22 +565,23 @@ mod test { } } - let mut green_value_deriv = rlst_dense::rlst_dynamic_mat!(f64, (nsources, 4 * ntargets)); + let mut green_value_deriv_t = rlst_dynamic_array2!(f64, [nsources, 4 * ntargets]); Laplace3dKernel::::default().assemble_st( EvalType::ValueDeriv, sources.data(), targets.data(), - green_value_deriv.data_mut(), + green_value_deriv_t.data_mut(), ); // The matrix needs to be transposed so that the first row corresponds to the first target, etc. - let green_value_deriv = green_value_deriv.transpose().eval(); + let mut green_value_deriv = rlst_dynamic_array2!(f64, [4 * ntargets, nsources]); + green_value_deriv.fill_from(green_value_deriv_t.transpose()); for charge_index in 0..nsources { - let mut charges = rlst::dense::rlst_col_vec![f64, nsources]; - let mut expected = rlst_dense::rlst_dynamic_mat!(f64, (4, ntargets)); + let mut charges = rlst_dynamic_array2![f64, [nsources, 1]]; + let mut expected = rlst_dynamic_array2!(f64, [4, ntargets]); charges[[charge_index, 0]] = 1.0; @@ -571,9 +610,9 @@ mod test { let nsources = 3; let ntargets = 5; - let sources = rlst::dense::rlst_rand_mat![f64, (nsources, 3)]; - let targets = rlst::dense::rlst_rand_mat![f64, (ntargets, 3)]; - let mut green_value_deriv = rlst_dense::rlst_dynamic_mat!(f64, (nsources, 4 * ntargets)); + let sources = rand_mat([nsources, 3]); + let targets = rand_mat([ntargets, 3]); + let mut green_value_deriv = rlst_dynamic_array2!(f64, [nsources, 4 * ntargets]); Laplace3dKernel::::default().assemble_st( EvalType::ValueDeriv, diff --git a/tools/Cargo.toml b/tools/Cargo.toml index ebc5b177..5e0d128d 100644 --- a/tools/Cargo.toml +++ b/tools/Cargo.toml @@ -24,4 +24,5 @@ libc = "0.2" num = "0.4" bempp-traits = { path = "../traits"} rayon = "1.7" -rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index be110f6e..9d3ba1e0 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -1,37 +1,34 @@ //! Containers to store multi-dimensional data use bempp_traits::arrays::{AdjacencyListAccess, Array3DAccess, Array4DAccess}; use num::Num; -use rlst_dense::{operations::transpose::Scalar, rlst_dynamic_mat, UnsafeRandomAccessMut}; +use rlst_common::{types::Scalar, traits::{UnsafeRandomAccessMut}}; +use rlst_dense::{rlst_dynamic_array2, array::Array, base_array::BaseArray, data_container::VectorContainer}; use std::clone::Clone; -pub type Mat = rlst_dense::Matrix< - T, - rlst_dense::base_matrix::BaseMatrix, rlst_dense::Dynamic>, - rlst_dense::Dynamic, ->; +pub type Mat = Array, 2>, 2>; -pub fn to_matrix(data: &[T], shape: (usize, usize)) -> Mat { - let mut mat = rlst_dynamic_mat![T, shape]; +pub fn to_matrix(data: &[T], shape: [usize; 2]) -> Mat { + let mut mat = rlst_dynamic_array2![T, shape]; for (i, d) in data.iter().enumerate() { unsafe { - *mat.get_unchecked_mut(i % shape.0, i / shape.0) = *d; + *mat.get_unchecked_mut([i % shape[0], i / shape[0]]) = *d; } } mat } -pub fn transpose_to_matrix(data: &[T], shape: (usize, usize)) -> Mat { - let mut mat = rlst_dynamic_mat![T, shape]; +pub fn transpose_to_matrix(data: &[T], shape: [usize; 2]) -> Mat { + let mut mat = rlst_dynamic_array2![T, shape]; for (i, d) in data.iter().enumerate() { unsafe { - *mat.get_unchecked_mut(i / shape.1, i % shape.1) = *d; + *mat.get_unchecked_mut([i / shape[1], i % shape[1]]) = *d; } } mat } -pub fn zero_matrix(shape: (usize, usize)) -> Mat { - rlst_dynamic_mat![T, shape] +pub fn zero_matrix(shape: [usize; 2]) -> Mat { + rlst_dynamic_array2![T, shape] } /// A three-dimensional rectangular array @@ -47,7 +44,7 @@ impl Array3D { /// Create an array from a data vector pub fn new(shape: (usize, usize, usize)) -> Self { Self { - data: vec![T::zero(); shape.0 * shape.1 * shape.2], + data: vec![T::zero(); shape.0 * shape. 1 * shape.2], shape, } } diff --git a/traits/Cargo.toml b/traits/Cargo.toml index d9af6c03..d5f89577 100644 --- a/traits/Cargo.toml +++ b/traits/Cargo.toml @@ -22,5 +22,5 @@ crate-type = ["lib", "cdylib"] cauchy="0.4.*" thiserror="1.*" num = "0.4" -rlst-common = { git = "https://github.com/linalg-rs/rlst.git" } -rlst-dense = { git = "https://github.com/linalg-rs/rlst.git" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } diff --git a/traits/src/element.rs b/traits/src/element.rs index 91bd4600..515e5508 100644 --- a/traits/src/element.rs +++ b/traits/src/element.rs @@ -71,7 +71,7 @@ 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 + Shape>( + fn tabulate + Shape<2>>( &self, points: &T, nderivs: usize, diff --git a/traits/src/grid.rs b/traits/src/grid.rs index 22893077..53c8273b 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -13,8 +13,8 @@ pub enum Ownership { } pub trait GeometryEvaluator< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, > { /// The points on the reference cell that this evaluator computes information at @@ -40,8 +40,8 @@ pub trait Geometry { //! Grid geometry //! //! Grid geometry provides information about the physical locations of mesh points in space - type T: RandomAccessByRef + Shape; - type TMut: RandomAccessByRef + RandomAccessMut + Shape; + type T: RandomAccessByRef<2, Item = f64> + Shape<2>; + type TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>; /// The geometric dimension fn dim(&self) -> usize; @@ -70,8 +70,8 @@ pub trait Geometry { /// Compute the physical coordinates of a set of points in a given cell fn compute_points< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, @@ -81,8 +81,8 @@ pub trait Geometry { /// Compute the normals to a set of points in a given cell fn compute_normals< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, @@ -94,8 +94,8 @@ pub trait Geometry { /// /// The input points should be given using coordinates on the reference element fn compute_jacobians< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, @@ -106,7 +106,7 @@ pub trait Geometry { /// 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 + Shape>( + fn compute_jacobian_determinants + Shape<2>>( &self, points: &T, cell: usize, @@ -117,8 +117,8 @@ pub trait Geometry { /// /// The input points should be given using coordinates on the reference element fn compute_jacobian_inverses< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, + T: RandomAccessByRef<2, Item = f64> + Shape<2>, + TMut: RandomAccessByRef<2, Item = f64> + RandomAccessMut<2, Item = f64> + Shape<2>, >( &self, points: &T, diff --git a/tree/Cargo.toml b/tree/Cargo.toml index 081da0bb..76690f7d 100644 --- a/tree/Cargo.toml +++ b/tree/Cargo.toml @@ -21,7 +21,9 @@ memoffset = "0.6" rand = "0.8.*" hyksort = { path = "../hyksort", optional = true } bempp-traits = { path = "../traits" } -rlst = { git = "https://github.com/linalg-rs/rlst.git" } +rlst = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-dense = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } +rlst-common = { git = "https://github.com/linalg-rs/rlst.git", branch="tensor" } num = "0.4" [features] diff --git a/tree/examples/test_adaptive.rs b/tree/examples/test_adaptive.rs index 65998399..bb380f39 100644 --- a/tree/examples/test_adaptive.rs +++ b/tree/examples/test_adaptive.rs @@ -10,7 +10,7 @@ use bempp_tree::implementations::helpers::points_fixture; use bempp_tree::types::{domain::Domain, morton::MortonKey, multi_node::MultiNodeTree}; use bempp_traits::types::Scalar; -use rlst::dense::RawAccess; +use rlst_common::traits::RawAccess; use rand::distributions::uniform::SampleUniform; diff --git a/tree/examples/test_uniform.rs b/tree/examples/test_uniform.rs index fc5c86c4..e45372e2 100644 --- a/tree/examples/test_uniform.rs +++ b/tree/examples/test_uniform.rs @@ -13,7 +13,7 @@ use rand::distributions::uniform::SampleUniform; use bempp_tree::types::{domain::Domain, morton::MortonKey, multi_node::MultiNodeTree}; use bempp_tree::implementations::helpers::points_fixture; -use rlst::dense::RawAccess; +use rlst_common::traits::RawAccess; use num::traits::Float; diff --git a/tree/src/implementations/helpers.rs b/tree/src/implementations/helpers.rs index 9b3d8183..4d40a30b 100644 --- a/tree/src/implementations/helpers.rs +++ b/tree/src/implementations/helpers.rs @@ -8,12 +8,12 @@ use num::Float; use rand::prelude::*; use rand::SeedableRng; -use rlst::dense::{base_matrix::BaseMatrix, rlst_dynamic_mat, Dynamic, Matrix, VectorContainer}; +use rlst_dense::{base_array::BaseArray, array::Array, rlst_dynamic_array2, data_container::VectorContainer}; use crate::types::morton::MortonKey; /// Alias for an rlst container for point data. -pub type PointsMat = Matrix, Dynamic>, Dynamic>; +pub type PointsMat = Array, 2>, 2>; /// Points fixture for testing, uniformly samples in each axis from min to max. /// @@ -36,7 +36,7 @@ pub fn points_fixture Domain { #[cfg(test)] mod test { use bempp_traits::types::Scalar; - use rlst::dense::{RawAccess, Shape}; + use rlst_common::traits::{RawAccess, Shape}; use crate::implementations::helpers::{points_fixture, points_fixture_col, PointsMat}; @@ -80,7 +80,7 @@ mod test { assert!(domain.diameter.iter().all(|&x| x == domain.diameter[0])); // Test that all local points are contained within the local domain - let npoints = points.shape().0; + let npoints = points.shape()[0]; for i in 0..npoints { let point = [points[[i, 0]], points[[i, 1]], points[[i, 2]]]; diff --git a/tree/src/implementations/impl_morton.rs b/tree/src/implementations/impl_morton.rs index b2765282..f5eaaf00 100644 --- a/tree/src/implementations/impl_morton.rs +++ b/tree/src/implementations/impl_morton.rs @@ -1011,7 +1011,7 @@ impl MortonKeyInterface for MortonKey { #[cfg(test)] mod test { use itertools::Itertools; - use rlst::dense::{RawAccess, Shape}; + use rlst_common::traits::{RawAccess, Shape}; use std::vec; use crate::implementations::helpers::points_fixture; @@ -1193,7 +1193,7 @@ mod test { let mut keys: Vec = Vec::new(); - for i in 0..points.shape().0 { + for i in 0..points.shape()[0] { let point = [points[[i, 0]], points[[i, 1]], points[[i, 2]]]; keys.push(MortonKey::from_point(&point, &domain, DEEPEST_LEVEL)); @@ -1466,7 +1466,7 @@ mod test { let mut keys = Vec::new(); - for i in 0..points.shape().0 { + for i in 0..points.shape()[0] { let point = [points[[i, 0]], points[[i, 1]], points[[i, 2]]]; keys.push(MortonKey::from_point(&point, &domain, DEEPEST_LEVEL)) } diff --git a/tree/src/implementations/impl_single_node.rs b/tree/src/implementations/impl_single_node.rs index b9c6a1a3..73de9696 100644 --- a/tree/src/implementations/impl_single_node.rs +++ b/tree/src/implementations/impl_single_node.rs @@ -575,7 +575,7 @@ where #[cfg(test)] mod test { - use rlst::dense::RawAccess; + use rlst_common::traits::RawAccess; use crate::implementations::helpers::{points_fixture, points_fixture_col};