From d3d1db63dc3b15879b13120277bb13aef6203ed4 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 25 Oct 2023 16:32:04 +0100 Subject: [PATCH 1/5] tabulate once before the loop instead of during every iterations --- bem/src/assembly/batched.rs | 84 ++++--------- grid/Cargo.toml | 1 + grid/src/grid.rs | 229 ++++++++++++++++++++++++++++++------ traits/src/grid.rs | 56 +++++++-- 4 files changed, 261 insertions(+), 109 deletions(-) diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index ff9ecbf8..6633efbf 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -13,7 +13,7 @@ use bempp_traits::element::FiniteElement; use bempp_traits::grid::{Geometry, Grid, Topology}; use bempp_traits::types::Scalar; use rayon::prelude::*; -use rlst_dense::{RandomAccessByRef, RawAccess, RawAccessMut, Shape, UnsafeRandomAccessMut}; +use rlst_dense::{RandomAccessByRef, RawAccess, RawAccessMut, Shape}; fn get_quadrature_rule( test_celltype: ReferenceCellType, @@ -209,36 +209,27 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz // let mut rlst_test_normals = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; // let mut rlst_trial_normals = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; - let mut rlst_test_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TEST, 3)]; - let mut rlst_trial_mapped_pts = rlst_dense::rlst_dynamic_mat![f64, (NPTS_TRIAL, 3)]; + let 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 test_element = test_grid.geometry().element(test_cells[0]); + let trial_element = trial_grid.geometry().element(trial_cells[0]); - let mut test_data = Array4D::::new(test_element.tabulate_array_shape(0, NPTS_TEST)); - test_element.tabulate(test_points, 0, &mut test_data); - let gdim = test_grid.geometry().dim(); + let test_compute_points = test_grid.geometry().get_compute_points_function(test_element, test_points); + let trial_compute_points = trial_grid.geometry().get_compute_points_function(trial_element, trial_points); - // TODO: move this to grid.get_compute_points_function(test_points) - let test_compute_points = |cell: usize, pts: &mut Mat| { - for p in 0..NPTS_TEST { - for i in 0..gdim { - unsafe { - *pts.get_unchecked_mut(p, i) = 0.0; - } - } - } - let vertices = test_grid.geometry().cell_vertices(cell).unwrap(); - for (i, n) in vertices.iter().enumerate() { - let pt = test_grid.geometry().point(*n).unwrap(); - for p in 0..NPTS_TEST { - for (j, pt_j) in pt.iter().enumerate() { - unsafe { - *pts.get_unchecked_mut(p, j) += - *pt_j * *test_data.get_unchecked(0, p, i, 0); - } - } - } - } + let mut test_compute_jdets = test_grid.geometry().get_compute_jacobian_determinants_function(test_element, test_points); + let mut trial_compute_jdets = trial_grid.geometry().get_compute_jacobian_determinants_function(trial_element, trial_points); + + let mut test_compute_normals = if needs_test_normal { + test_grid.geometry().get_compute_normals_function(test_element, test_points) + } else { + Box::new(|_i: usize, _m: &mut Mat| {}) + }; + let mut trial_compute_normals = if needs_trial_normal { + trial_grid.geometry().get_compute_normals_function(trial_element, trial_points) + } else { + Box::new(|_i: usize, _m: &mut Mat| {}) }; for test_cell in test_cells { @@ -246,46 +237,23 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz let test_cell_gindex = test_grid.geometry().index_map()[*test_cell]; let test_vertices = unsafe { test_c20.row_unchecked(test_cell_tindex) }; - test_grid.geometry().compute_jacobian_determinants( - test_points, - test_cell_gindex, - &mut test_jdet, - ); - test_compute_points(test_cell_gindex, &mut rlst_test_mapped_pts); - - if needs_test_normal { - test_grid - .geometry() - .compute_normals(test_points, test_cell_gindex, &mut test_normals); - } + test_compute_jdets(test_cell_gindex, &mut test_jdet); + test_compute_points(test_cell_gindex, &mut test_mapped_pts); + test_compute_normals(test_cell_gindex, &mut test_normals); for trial_cell in trial_cells { let trial_cell_tindex = trial_grid.topology().index_map()[*trial_cell]; let trial_cell_gindex = trial_grid.geometry().index_map()[*trial_cell]; let trial_vertices = unsafe { trial_c20.row_unchecked(trial_cell_tindex) }; - trial_grid.geometry().compute_jacobian_determinants( - trial_points, - trial_cell_gindex, - &mut trial_jdet, - ); - trial_grid.geometry().compute_points( - trial_points, - trial_cell_gindex, - &mut rlst_trial_mapped_pts, - ); - if needs_trial_normal { - trial_grid.geometry().compute_normals( - trial_points, - trial_cell_gindex, - &mut trial_normals, - ); - } + trial_compute_jdets(test_cell_gindex, &mut trial_jdet); + trial_compute_points(trial_cell_gindex, &mut trial_mapped_pts); + trial_compute_normals(trial_cell_gindex, &mut trial_normals); kernel.assemble_st( EvalType::Value, - rlst_test_mapped_pts.data(), - rlst_trial_mapped_pts.data(), + test_mapped_pts.data(), + trial_mapped_pts.data(), &mut k, ); diff --git a/grid/Cargo.toml b/grid/Cargo.toml index 96d3287b..f540f853 100644 --- a/grid/Cargo.toml +++ b/grid/Cargo.toml @@ -27,4 +27,5 @@ 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" } diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 3973bc10..4d800e05 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -1,14 +1,16 @@ //! A serial implementation of a grid use bempp_element::cell; use bempp_element::element::{create_element, CiarletElement}; -use bempp_tools::arrays::{to_matrix, AdjacencyList, Array4D, Mat}; +use bempp_tools::arrays::{zero_matrix, AdjacencyList, Array4D, Mat}; use bempp_traits::arrays::{AdjacencyListAccess, Array4DAccess}; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; use itertools::izip; -use rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape}; +use rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape, rlst_static_mat, SizeIdentifier, RawAccess}; +use rlst_proc_macro::rlst_static_size; use std::ptr; + /// Geometry of a serial grid pub struct SerialGeometry { coordinate_elements: Vec, @@ -18,6 +20,10 @@ 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, @@ -119,36 +125,37 @@ impl Geometry for SerialGeometry { fn index_map(&self) -> &[usize] { &self.index_map } - /* - fn get_compute_points_function + Shape>( - &self, - element: &impl FiniteElement, - points: &T, - ) -> Box { - let npts = points.shape().0; - let mut table = Array4D::::new(element.tabulate_array_shape(0, npts)); - element.tabulate(points, 0, &mut table); - let gdim = self.dim(); - - Box::new(|cell: usize, pts: &mut T| { - for p in 0..npts { - for i in 0..gdim { - *pts.get_mut(p, i).unwrap() = 0.0; - } + fn get_compute_points_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box { + let npts = points.shape().0; + let mut table = Array4D::::new(element.tabulate_array_shape(0, npts)); + element.tabulate(points, 0, &mut table); + let gdim = self.dim(); + + Box::new(move |cell: usize, pts: &mut TMut| { + for p in 0..npts { + for i in 0..gdim { + *pts.get_mut(p, i).unwrap() = 0.0; } - let vertices = self.cell_vertices(cell).unwrap(); - for (i, n) in vertices.iter().enumerate() { - let pt = self.point(*n).unwrap(); - for p in 0..points.shape().0 { - for (j, pt_j) in pt.iter().enumerate() { - *pts.get_mut(p, j).unwrap() += - *pt_j * *table.get(0, p, i, 0); - } + } + let vertices = self.cell_vertices(cell).unwrap(); + for (i, n) in vertices.iter().enumerate() { + let pt = self.point(*n).unwrap(); + for p in 0..points.shape().0 { + for (j, pt_j) in pt.iter().enumerate() { + *pts.get_mut(p, j).unwrap() += + *pt_j * *table.get(0, p, i, 0).unwrap(); } } - }) - } - */ + } + }) + } fn compute_points< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -175,14 +182,55 @@ impl Geometry for SerialGeometry { } for i in 0..data.shape().2 { let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); - for p in 0..points.shape().0 { - for (j, pt_j) in pt.iter().enumerate() { + for (j, pt_j) in pt.iter().enumerate() { + for p in 0..points.shape().0 { *physical_points.get_mut(p, j).unwrap() += *pt_j * data.get(0, p, i, 0).unwrap(); } } } } + fn get_compute_normals_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box { + let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? + let mut axes = rlst_static_mat![f64, TwoByThree]; + element.tabulate(points, 1, &mut data); + Box::new(move |cell: usize, normals: &mut TMut| { + for p in 0..points.shape().0 { + for i in 0..axes.shape().0 { + for j in 0..axes.shape().1 { + *axes.get_mut(i, j).unwrap() = 0.0; + } + } + for i in 0..data.shape().2 { + let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); + for (j, pt_j) in pt.iter().enumerate() { + *axes.get_mut(0, j).unwrap() += *pt_j * data.get(1, p, i, 0).unwrap(); + *axes.get_mut(1, j).unwrap() += *pt_j * data.get(2, p, i, 0).unwrap(); + } + } + *normals.get_mut(p, 0).unwrap() = *axes.get(0, 1).unwrap() * *axes.get(1, 2).unwrap() + - *axes.get(0, 2).unwrap() * *axes.get(1, 1).unwrap(); + *normals.get_mut(p, 1).unwrap() = *axes.get(0, 2).unwrap() * *axes.get(1, 0).unwrap() + - *axes.get(0, 0).unwrap() * *axes.get(1, 2).unwrap(); + *normals.get_mut(p, 2).unwrap() = *axes.get(0, 0).unwrap() * *axes.get(1, 1).unwrap() + - *axes.get(0, 1).unwrap() * *axes.get(1, 0).unwrap(); + let size = (*normals.get(p, 0).unwrap() * *normals.get(p, 0).unwrap() + + *normals.get(p, 1).unwrap() * *normals.get(p, 1).unwrap() + + *normals.get(p, 2).unwrap() * *normals.get(p, 2).unwrap()) + .sqrt(); + *normals.get_mut(p, 0).unwrap() /= size; + *normals.get_mut(p, 1).unwrap() /= size; + *normals.get_mut(p, 2).unwrap() /= size; + } + }) + } fn compute_normals< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -204,7 +252,7 @@ impl Geometry for SerialGeometry { } let element = self.element(cell); let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? - let mut axes = to_matrix(&[0.0; 6], (2, 3)); + let mut axes = rlst_static_mat![f64, TwoByThree]; element.tabulate(points, 1, &mut data); for p in 0..points.shape().0 { for i in 0..axes.shape().0 { @@ -234,6 +282,37 @@ impl Geometry for SerialGeometry { *normals.get_mut(p, 2).unwrap() /= size; } } + fn get_compute_jacobians_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box { + let tdim = points.shape().1; + let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? + element.tabulate(points, 1, &mut data); + + Box::new(move |cell: usize, jacobians: &mut TMut| { + for p in 0..points.shape().0 { + for i in 0..jacobians.shape().0 { + *jacobians.get_mut(i, p).unwrap() = 0.0; + } + } + for i in 0..data.shape().2 { + let pt = self.point(*self.cells.get(cell, i).unwrap()).unwrap(); + for p in 0..points.shape().0 { + for (j, pt_j) in pt.iter().enumerate() { + for k in 0..tdim { + *jacobians.get_mut(k + tdim * j, p).unwrap() += + *pt_j * data.get(k + 1, p, i, 0).unwrap(); + } + } + } + } + }) + } fn compute_jacobians< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -272,6 +351,82 @@ impl Geometry for SerialGeometry { } } } + fn get_compute_jacobian_determinants_function<'a, + T: RandomAccessByRef + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box { + + let gdim = self.dim(); + let tdim = points.shape().1; + let mut js = zero_matrix((gdim * tdim, points.shape().0)); + + let det = match tdim { + 1 => match gdim { + 1 => |x: &[f64]| x[0], + 2 => { |x: &[f64]| + ((x[0]).powi(2) + (x[1]).powi(2)).sqrt() + } + 3 => |x: &[f64]| ((x[0]).powi(2) + + (x[1]).powi(2) + + (x[2]).powi(2)) + .sqrt(), + _ => { + panic!("Unsupported dimensions."); + } + }, + 2 => match gdim { + 2 => { |x: &[f64]| + x[0] * x[3] + - x[1] * x[2] + } + 3 => |x: &[f64]| (((x[0]).powi(2) + + (x[2]).powi(2) + + (x.get(4).unwrap()).powi(2)) + * ((x[1]).powi(2) + + (x[3]).powi(2) + + (x[5]).powi(2)) + - (x[0] * x[1] + + x[2] * x[3] + + x.get(4).unwrap() * x[5]) + .powi(2)) + .sqrt(), + _ => { + panic!("Unsupported dimensions."); + } + }, + 3 => match gdim { + 3 => { |x: &[f64]| + x[0] + * (x.get(4).unwrap() * x[8] + - x[5] * x[7]) + - x[1] + * (x[3] * x[8] + - x[5] * x[6]) + + x[2] + * (x[3] * x[7] + - x.get(4).unwrap() * x[6]) + } + _ => { + panic!("Unsupported dimensions."); + } + }, + _ => { + panic!("Unsupported dimensions."); + } + }; + + let compute_jacobians = self.get_compute_jacobians_function(element, points); + + Box::new(move |cell: usize, jacobian_determinants: &mut [f64]| { + compute_jacobians(cell, &mut js); + for (p, jdet) in jacobian_determinants.iter_mut().enumerate() { + *jdet = det(&js.data()[tdim * gdim * p.. tdim * gdim * (p+1)]); + } + }) + } fn compute_jacobian_determinants + Shape>( &self, points: &T, @@ -283,9 +438,7 @@ impl Geometry for SerialGeometry { if points.shape().0 != jacobian_determinants.len() { panic!("jacobian_determinants has wrong length."); } - let mut js = to_matrix( - &vec![0.0; points.shape().0 * gdim * tdim], - (points.shape().0, gdim * tdim), + let mut js = zero_matrix((points.shape().0, gdim * tdim), ); // TODO: Memory is assigned here. Can we avoid this? self.compute_jacobians(points, cell, &mut js); @@ -370,10 +523,7 @@ impl Geometry for SerialGeometry { && element.degree() == 1 { // Map is affine - let mut js = to_matrix( - &vec![0.0; points.shape().0 * gdim * tdim], - (points.shape().0, gdim * tdim), - ); // TODO: Memory is assigned here. Can we avoid this? + let mut js = zero_matrix((points.shape().0, gdim * tdim)); // TODO: Memory is assigned here. Can we avoid this? self.compute_jacobians(points, cell, &mut js); // TODO: is it faster if we move this for inside the if statement? @@ -815,6 +965,7 @@ impl Eq for SerialGrid {} mod test { use crate::grid::*; use approx::*; + use bempp_tools::arrays::to_matrix; #[test] fn test_connectivity() { diff --git a/traits/src/grid.rs b/traits/src/grid.rs index 47a85df7..b8001088 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -2,7 +2,7 @@ use crate::arrays::AdjacencyListAccess; use crate::cell::ReferenceCellType; -// use crate::element::FiniteElement; +use crate::element::FiniteElement; use rlst_common::traits::{RandomAccessByRef, RandomAccessMut, Shape}; /// The ownership of a mesh entity @@ -35,17 +35,16 @@ pub trait Geometry { /// Return the index map from the input cell numbers to the storage numbers fn index_map(&self) -> &[usize]; - /* - /// Compute the physical coordinates of a set of points in a given cell - fn get_compute_points_function< - T: RandomAccessByRef + Shape, - TMut: RandomAccessByRef + RandomAccessMut + Shape, - >( - &self, - element: &impl FiniteElement, - points: &T, - ) -> Box; - */ + /// Get function that computes the physical coordinates of a set of points in a given cell + fn get_compute_points_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box; + /// Compute the physical coordinates of a set of points in a given cell fn compute_points< T: RandomAccessByRef + Shape, @@ -57,6 +56,16 @@ pub trait Geometry { physical_points: &mut TMut, ); + /// Get function that computes the normals to a set of points in a given cell + fn get_compute_normals_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box; + /// Compute the normals to a set of points in a given cell fn compute_normals< T: RandomAccessByRef + Shape, @@ -68,6 +77,18 @@ pub trait Geometry { normals: &mut TMut, ); + /// Get function that evaluates the jacobian at a set of points in a given cell + /// + /// The input points should be given using coordinates on the reference element + fn get_compute_jacobians_function<'a, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box; + /// Evaluate the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element @@ -81,6 +102,17 @@ pub trait Geometry { jacobians: &mut TMut, ); + /// Get function that evaluates 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 get_compute_jacobian_determinants_function<'a, + T: RandomAccessByRef + Shape, + >( + &'a self, + element: &impl FiniteElement, + points: &'a T, + ) -> Box; + /// 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 From da8775e1b7cc6351485690f41043e67425f190d6 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 25 Oct 2023 16:32:27 +0100 Subject: [PATCH 2/5] fmt --- bem/src/assembly/batched.rs | 26 +++++++++--- grid/src/grid.rs | 85 +++++++++++++++---------------------- traits/src/grid.rs | 13 +++--- 3 files changed, 60 insertions(+), 64 deletions(-) diff --git a/bem/src/assembly/batched.rs b/bem/src/assembly/batched.rs index 6633efbf..a0c21890 100644 --- a/bem/src/assembly/batched.rs +++ b/bem/src/assembly/batched.rs @@ -215,19 +215,31 @@ fn assemble_batch_nonadjacent<'a, const NPTS_TEST: usize, const NPTS_TRIAL: usiz let test_element = test_grid.geometry().element(test_cells[0]); let trial_element = trial_grid.geometry().element(trial_cells[0]); - let test_compute_points = test_grid.geometry().get_compute_points_function(test_element, test_points); - let trial_compute_points = trial_grid.geometry().get_compute_points_function(trial_element, trial_points); - - let mut test_compute_jdets = test_grid.geometry().get_compute_jacobian_determinants_function(test_element, test_points); - let mut trial_compute_jdets = trial_grid.geometry().get_compute_jacobian_determinants_function(trial_element, trial_points); + let test_compute_points = test_grid + .geometry() + .get_compute_points_function(test_element, test_points); + let trial_compute_points = trial_grid + .geometry() + .get_compute_points_function(trial_element, trial_points); + + let mut test_compute_jdets = test_grid + .geometry() + .get_compute_jacobian_determinants_function(test_element, test_points); + let mut trial_compute_jdets = trial_grid + .geometry() + .get_compute_jacobian_determinants_function(trial_element, trial_points); let mut test_compute_normals = if needs_test_normal { - test_grid.geometry().get_compute_normals_function(test_element, test_points) + test_grid + .geometry() + .get_compute_normals_function(test_element, test_points) } else { Box::new(|_i: usize, _m: &mut Mat| {}) }; let mut trial_compute_normals = if needs_trial_normal { - trial_grid.geometry().get_compute_normals_function(trial_element, trial_points) + trial_grid + .geometry() + .get_compute_normals_function(trial_element, trial_points) } else { Box::new(|_i: usize, _m: &mut Mat| {}) }; diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 4d800e05..6c99d5b4 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -7,7 +7,9 @@ use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; use itertools::izip; -use rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape, rlst_static_mat, SizeIdentifier, RawAccess}; +use rlst_dense::{ + rlst_static_mat, RandomAccessByRef, RandomAccessMut, RawAccess, Shape, SizeIdentifier, +}; use rlst_proc_macro::rlst_static_size; use std::ptr; @@ -23,7 +25,6 @@ pub struct SerialGeometry { #[rlst_static_size(2, 3)] struct TwoByThree; - fn element_from_npts(cell_type: ReferenceCellType, npts: usize) -> CiarletElement { create_element( ElementFamily::Lagrange, @@ -125,7 +126,8 @@ impl Geometry for SerialGeometry { fn index_map(&self) -> &[usize] { &self.index_map } - fn get_compute_points_function<'a, + fn get_compute_points_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -141,7 +143,7 @@ impl Geometry for SerialGeometry { Box::new(move |cell: usize, pts: &mut TMut| { for p in 0..npts { for i in 0..gdim { - *pts.get_mut(p, i).unwrap() = 0.0; + *pts.get_mut(p, i).unwrap() = 0.0; } } let vertices = self.cell_vertices(cell).unwrap(); @@ -149,8 +151,7 @@ impl Geometry for SerialGeometry { let pt = self.point(*n).unwrap(); for p in 0..points.shape().0 { for (j, pt_j) in pt.iter().enumerate() { - *pts.get_mut(p, j).unwrap() += - *pt_j * *table.get(0, p, i, 0).unwrap(); + *pts.get_mut(p, j).unwrap() += *pt_j * *table.get(0, p, i, 0).unwrap(); } } } @@ -190,7 +191,8 @@ impl Geometry for SerialGeometry { } } } - fn get_compute_normals_function<'a, + fn get_compute_normals_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -215,11 +217,14 @@ impl Geometry for SerialGeometry { *axes.get_mut(1, j).unwrap() += *pt_j * data.get(2, p, i, 0).unwrap(); } } - *normals.get_mut(p, 0).unwrap() = *axes.get(0, 1).unwrap() * *axes.get(1, 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() + *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() + *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() @@ -282,7 +287,8 @@ impl Geometry for SerialGeometry { *normals.get_mut(p, 2).unwrap() /= size; } } - fn get_compute_jacobians_function<'a, + fn get_compute_jacobians_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -351,14 +357,11 @@ impl Geometry for SerialGeometry { } } } - fn get_compute_jacobian_determinants_function<'a, - T: RandomAccessByRef + Shape, - >( + fn get_compute_jacobian_determinants_function<'a, T: RandomAccessByRef + Shape>( &'a self, element: &impl FiniteElement, points: &'a T, ) -> Box { - let gdim = self.dim(); let tdim = points.shape().1; let mut js = zero_matrix((gdim * tdim, points.shape().0)); @@ -366,49 +369,30 @@ impl Geometry for SerialGeometry { let det = match tdim { 1 => match gdim { 1 => |x: &[f64]| x[0], - 2 => { |x: &[f64]| - ((x[0]).powi(2) + (x[1]).powi(2)).sqrt() - } - 3 => |x: &[f64]| ((x[0]).powi(2) - + (x[1]).powi(2) - + (x[2]).powi(2)) - .sqrt(), + 2 => |x: &[f64]| ((x[0]).powi(2) + (x[1]).powi(2)).sqrt(), + 3 => |x: &[f64]| ((x[0]).powi(2) + (x[1]).powi(2) + (x[2]).powi(2)).sqrt(), _ => { panic!("Unsupported dimensions."); } }, 2 => match gdim { - 2 => { |x: &[f64]| - x[0] * x[3] - - x[1] * x[2] - } - 3 => |x: &[f64]| (((x[0]).powi(2) - + (x[2]).powi(2) - + (x.get(4).unwrap()).powi(2)) - * ((x[1]).powi(2) - + (x[3]).powi(2) - + (x[5]).powi(2)) - - (x[0] * x[1] - + x[2] * x[3] - + x.get(4).unwrap() * x[5]) - .powi(2)) - .sqrt(), + 2 => |x: &[f64]| x[0] * x[3] - x[1] * x[2], + 3 => |x: &[f64]| { + (((x[0]).powi(2) + (x[2]).powi(2) + (x.get(4).unwrap()).powi(2)) + * ((x[1]).powi(2) + (x[3]).powi(2) + (x[5]).powi(2)) + - (x[0] * x[1] + x[2] * x[3] + x.get(4).unwrap() * x[5]).powi(2)) + .sqrt() + }, _ => { panic!("Unsupported dimensions."); } }, 3 => match gdim { - 3 => { |x: &[f64]| - x[0] - * (x.get(4).unwrap() * x[8] - - x[5] * x[7]) - - x[1] - * (x[3] * x[8] - - x[5] * x[6]) - + x[2] - * (x[3] * x[7] - - x.get(4).unwrap() * x[6]) - } + 3 => |x: &[f64]| { + x[0] * (x.get(4).unwrap() * x[8] - x[5] * x[7]) + - x[1] * (x[3] * x[8] - x[5] * x[6]) + + x[2] * (x[3] * x[7] - x.get(4).unwrap() * x[6]) + }, _ => { panic!("Unsupported dimensions."); } @@ -423,7 +407,7 @@ impl Geometry for SerialGeometry { Box::new(move |cell: usize, jacobian_determinants: &mut [f64]| { compute_jacobians(cell, &mut js); for (p, jdet) in jacobian_determinants.iter_mut().enumerate() { - *jdet = det(&js.data()[tdim * gdim * p.. tdim * gdim * (p+1)]); + *jdet = det(&js.data()[tdim * gdim * p..tdim * gdim * (p + 1)]); } }) } @@ -438,8 +422,7 @@ impl Geometry for SerialGeometry { if points.shape().0 != jacobian_determinants.len() { panic!("jacobian_determinants has wrong length."); } - let mut js = zero_matrix((points.shape().0, gdim * tdim), - ); // TODO: Memory is assigned here. Can we avoid this? + let mut js = zero_matrix((points.shape().0, gdim * tdim)); // TODO: Memory is assigned here. Can we avoid this? self.compute_jacobians(points, cell, &mut js); // TODO: is it faster if we move this for inside the match statement? diff --git a/traits/src/grid.rs b/traits/src/grid.rs index b8001088..c23009cb 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -36,7 +36,8 @@ pub trait Geometry { fn index_map(&self) -> &[usize]; /// Get function that computes the physical coordinates of a set of points in a given cell - fn get_compute_points_function<'a, + fn get_compute_points_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -57,7 +58,8 @@ pub trait Geometry { ); /// Get function that computes the normals to a set of points in a given cell - fn get_compute_normals_function<'a, + fn get_compute_normals_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -80,7 +82,8 @@ pub trait Geometry { /// Get function that evaluates the jacobian at a set of points in a given cell /// /// The input points should be given using coordinates on the reference element - fn get_compute_jacobians_function<'a, + fn get_compute_jacobians_function< + 'a, T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, >( @@ -105,9 +108,7 @@ pub trait Geometry { /// Get function that evaluates 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 get_compute_jacobian_determinants_function<'a, - T: RandomAccessByRef + Shape, - >( + fn get_compute_jacobian_determinants_function<'a, T: RandomAccessByRef + Shape>( &'a self, element: &impl FiniteElement, points: &'a T, From 60ee344b692969cf5182355b4ffe3316f47ca7a7 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 25 Oct 2023 16:39:27 +0100 Subject: [PATCH 3/5] clippy --- grid/src/grid.rs | 10 +++++----- traits/src/grid.rs | 11 +++++++---- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/grid/src/grid.rs b/grid/src/grid.rs index 6c99d5b4..a780c122 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -5,7 +5,7 @@ use bempp_tools::arrays::{zero_matrix, AdjacencyList, Array4D, Mat}; use bempp_traits::arrays::{AdjacencyListAccess, Array4DAccess}; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; -use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; +use bempp_traits::grid::{Geometry, Grid, Ownership, Topology, GeomF, GeomFMut}; use itertools::izip; use rlst_dense::{ rlst_static_mat, RandomAccessByRef, RandomAccessMut, RawAccess, Shape, SizeIdentifier, @@ -134,7 +134,7 @@ impl Geometry for SerialGeometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box { + ) -> GeomF<'a, TMut> { let npts = points.shape().0; let mut table = Array4D::::new(element.tabulate_array_shape(0, npts)); element.tabulate(points, 0, &mut table); @@ -199,7 +199,7 @@ impl Geometry for SerialGeometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box { + ) -> GeomFMut<'a, TMut> { let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? let mut axes = rlst_static_mat![f64, TwoByThree]; element.tabulate(points, 1, &mut data); @@ -295,7 +295,7 @@ impl Geometry for SerialGeometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box { + ) -> GeomF<'a, TMut> { let tdim = points.shape().1; let mut data = Array4D::::new(element.tabulate_array_shape(1, points.shape().0)); // TODO: Memory is assigned here. Can we avoid this? element.tabulate(points, 1, &mut data); @@ -361,7 +361,7 @@ impl Geometry for SerialGeometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box { + ) -> GeomFMut<'a, [f64]> { let gdim = self.dim(); let tdim = points.shape().1; let mut js = zero_matrix((gdim * tdim, points.shape().0)); diff --git a/traits/src/grid.rs b/traits/src/grid.rs index c23009cb..f9c39fe5 100644 --- a/traits/src/grid.rs +++ b/traits/src/grid.rs @@ -12,6 +12,9 @@ pub enum Ownership { Ghost(usize, usize), } +pub type GeomF<'a, T> = Box; +pub type GeomFMut<'a, T> = Box; + pub trait Geometry { //! Grid geometry //! @@ -44,7 +47,7 @@ pub trait Geometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box; + ) -> GeomF<'a, TMut>; /// Compute the physical coordinates of a set of points in a given cell fn compute_points< @@ -66,7 +69,7 @@ pub trait Geometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box; + ) -> GeomFMut<'a, TMut>; /// Compute the normals to a set of points in a given cell fn compute_normals< @@ -90,7 +93,7 @@ pub trait Geometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box; + ) -> GeomF<'a, TMut>; /// Evaluate the jacobian at a set of points in a given cell /// @@ -112,7 +115,7 @@ pub trait Geometry { &'a self, element: &impl FiniteElement, points: &'a T, - ) -> Box; + ) -> GeomFMut<[f64]>; /// Evaluate the determinand of the jacobian at a set of points in a given cell /// From 635d7076edd7110efe7265dd3bec8a7068406c09 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 25 Oct 2023 17:08:04 +0100 Subject: [PATCH 4/5] fmt --- grid/src/grid.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/grid/src/grid.rs b/grid/src/grid.rs index a780c122..9e733bc1 100644 --- a/grid/src/grid.rs +++ b/grid/src/grid.rs @@ -5,7 +5,7 @@ use bempp_tools::arrays::{zero_matrix, AdjacencyList, Array4D, Mat}; use bempp_traits::arrays::{AdjacencyListAccess, Array4DAccess}; use bempp_traits::cell::{ReferenceCell, ReferenceCellType}; use bempp_traits::element::{Continuity, ElementFamily, FiniteElement}; -use bempp_traits::grid::{Geometry, Grid, Ownership, Topology, GeomF, GeomFMut}; +use bempp_traits::grid::{GeomF, GeomFMut, Geometry, Grid, Ownership, Topology}; use itertools::izip; use rlst_dense::{ rlst_static_mat, RandomAccessByRef, RandomAccessMut, RawAccess, Shape, SizeIdentifier, From 33d2a0fac0b86989886531b8d1c2fb935a0252f9 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 25 Oct 2023 17:12:06 +0100 Subject: [PATCH 5/5] fix parallel grid --- grid/src/parallel_grid.rs | 47 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/grid/src/parallel_grid.rs b/grid/src/parallel_grid.rs index 4edc62da..f561edd0 100644 --- a/grid/src/parallel_grid.rs +++ b/grid/src/parallel_grid.rs @@ -4,7 +4,8 @@ use bempp_element::element::CiarletElement; use bempp_tools::arrays::{zero_matrix, AdjacencyList, Mat}; use bempp_traits::arrays::AdjacencyListAccess; use bempp_traits::cell::ReferenceCellType; -use bempp_traits::grid::{Geometry, Grid, Ownership, Topology}; +use bempp_traits::element::FiniteElement; +use bempp_traits::grid::{GeomF, GeomFMut, Geometry, Grid, Ownership, Topology}; use mpi::{request::WaitGuard, topology::Communicator, traits::*}; use rlst_dense::{RandomAccessByRef, RandomAccessMut, Shape}; @@ -62,6 +63,18 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { fn index_map(&self) -> &[usize] { self.serial_geometry.index_map() } + fn get_compute_points_function< + 'b, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'b self, + element: &impl FiniteElement, + points: &'b T, + ) -> GeomF<'b, TMut> { + self.serial_geometry + .get_compute_points_function(element, points) + } fn compute_points< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -74,6 +87,18 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { self.serial_geometry .compute_points(points, cell, physical_points) } + fn get_compute_normals_function< + 'b, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'b self, + element: &impl FiniteElement, + points: &'b T, + ) -> GeomFMut<'b, TMut> { + self.serial_geometry + .get_compute_normals_function(element, points) + } fn compute_normals< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -85,6 +110,18 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { ) { self.serial_geometry.compute_points(points, cell, normals) } + fn get_compute_jacobians_function< + 'b, + T: RandomAccessByRef + Shape, + TMut: RandomAccessByRef + RandomAccessMut + Shape, + >( + &'b self, + element: &impl FiniteElement, + points: &'b T, + ) -> GeomF<'b, TMut> { + self.serial_geometry + .get_compute_jacobians_function(element, points) + } fn compute_jacobians< T: RandomAccessByRef + Shape, TMut: RandomAccessByRef + RandomAccessMut + Shape, @@ -97,6 +134,14 @@ impl<'a, C: Communicator> Geometry for ParallelGeometry<'a, C> { self.serial_geometry .compute_jacobians(points, cell, jacobians) } + fn get_compute_jacobian_determinants_function<'b, T: RandomAccessByRef + Shape>( + &'b self, + element: &impl FiniteElement, + points: &'b T, + ) -> GeomFMut<'b, [f64]> { + self.serial_geometry + .get_compute_jacobian_determinants_function(element, points) + } fn compute_jacobian_determinants + Shape>( &self, points: &T,