diff --git a/examples/assembly.rs b/examples/assembly.rs index 8dfd129a..9b792ca7 100644 --- a/examples/assembly.rs +++ b/examples/assembly.rs @@ -25,7 +25,7 @@ fn main() { ); // Assemble the single layer Laplace boundary operator. - let matrix = laplace_single_layer(&space, &space, &options); + let matrix = laplace_single_layer(&options).assemble(&space, &space); // Print the entries of the matrix println!("Lagrange single layer matrix"); diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 274661e0..bea418e6 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -27,8 +27,8 @@ use ndgrid::traits::{Entity, Grid, Topology}; use ndgrid::types::Ownership; use rayon::prelude::*; use rlst::{ - rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, MatrixInverse, - RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, + rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, DynamicArray, + MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, }; use std::collections::HashMap; @@ -122,26 +122,20 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: } /// Assemble into a dense matrix. - pub fn assemble< - Space: FunctionSpace, - Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut, - >( + pub fn assemble>( &self, - output: &mut Array2, trial_space: &Space, test_space: &Space, - ) { + ) -> DynamicArray { let test_colouring = test_space.cell_colouring(); let trial_colouring = trial_space.cell_colouring(); if !trial_space.is_serial() || !test_space.is_serial() { panic!("Dense assembly can only be used for function spaces stored in serial"); } - if output.shape()[0] != test_space.global_size() - || output.shape()[1] != trial_space.global_size() - { - panic!("Matrix has wrong shape"); - } + + let mut output = + rlst_dynamic_array2!(T, [test_space.global_size(), trial_space.global_size()]); let output_raw = RawData2D { data: output.data_mut().as_mut_ptr(), @@ -164,6 +158,8 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { *output.get_mut([*i, *j]).unwrap() += *value; } + + output } /// Create new Boundary assembler diff --git a/src/helmholtz.rs b/src/helmholtz.rs index 10488d83..a790bf18 100644 --- a/src/helmholtz.rs +++ b/src/helmholtz.rs @@ -22,109 +22,63 @@ pub mod assembler { }; /// Assembler for the Helmholtz single layer operator. - pub fn helmholtz_single_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( + pub fn helmholtz_single_layer + MatrixInverse>( wavenumber: T::Real, - trial_space: &Space, - test_space: &Space, options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Helmholtz3dKernel> { let kernel = KernelEvaluator::new( Helmholtz3dKernel::new(wavenumber), GreenKernelEvalType::Value, ); - let assembler = - BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); - - assembler.assemble(&mut output, trial_space, test_space); - - output + BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0) } /// Assembler for the Helmholtz double layer operator. - pub fn helmholtz_double_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( + pub fn helmholtz_double_layer + MatrixInverse>( wavenumber: T::Real, - trial_space: &Space, - test_space: &Space, options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Helmholtz3dKernel> { let kernel = KernelEvaluator::new( Helmholtz3dKernel::new(wavenumber), GreenKernelEvalType::ValueDeriv, ); - let assembler = - BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); - - assembler.assemble(&mut output, trial_space, test_space); - - output + BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0) } /// Assembler for the Helmholtz adjoint double layer operator. - pub fn helmholtz_adjoint_double_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( + pub fn helmholtz_adjoint_double_layer + MatrixInverse>( wavenumber: T::Real, - trial_space: &Space, - test_space: &Space, options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Helmholtz3dKernel> { let kernel = KernelEvaluator::new( Helmholtz3dKernel::new(wavenumber), GreenKernelEvalType::ValueDeriv, ); - let assembler = BoundaryAssembler::new( + BoundaryAssembler::new( AdjointDoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0, - ); - - assembler.assemble(&mut output, trial_space, test_space); - - output + ) } /// Assembler for the Helmholtz hypersingular operator. - pub fn helmholtz_hypersingular< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( + pub fn helmholtz_hypersingular + MatrixInverse>( wavenumber: T::Real, - trial_space: &Space, - test_space: &Space, options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler< + T, + BoundaryIntegrandSum< + T, + HypersingularCurlCurlBoundaryIntegrand, + BoundaryIntegrandTimesScalar>, + >, + Helmholtz3dKernel, + > { let kernel = KernelEvaluator::new( Helmholtz3dKernel::new(wavenumber), GreenKernelEvalType::ValueDeriv, @@ -138,10 +92,6 @@ pub mod assembler { ), ); - let assembler = BoundaryAssembler::new(integrand, kernel, options, 4, 1); - - assembler.assemble(&mut output, trial_space, test_space); - - output + BoundaryAssembler::new(integrand, kernel, options, 4, 1) } } diff --git a/src/laplace.rs b/src/laplace.rs index e97fea6c..488010d6 100644 --- a/src/laplace.rs +++ b/src/laplace.rs @@ -20,108 +20,49 @@ pub mod assembler { }; /// Assembler for the Laplace single layer operator. - pub fn laplace_single_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( - trial_space: &Space, - test_space: &Space, + pub fn laplace_single_layer + MatrixInverse>( options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Laplace3dKernel> { let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value); - - let assembler = - BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); - - assembler.assemble(&mut output, trial_space, test_space); - - output + BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0) } /// Assembler for the Laplace double layer operator. - pub fn laplace_double_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( - trial_space: &Space, - test_space: &Space, + pub fn laplace_double_layer + MatrixInverse>( options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Laplace3dKernel> { let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); - let assembler = - BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); - - assembler.assemble(&mut output, trial_space, test_space); - - output + BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0) } /// Assembler for the Laplace adjoint double layer operator. - pub fn laplace_adjoint_double_layer< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( - trial_space: &Space, - test_space: &Space, + pub fn laplace_adjoint_double_layer + MatrixInverse>( options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Laplace3dKernel> { let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); - let assembler = BoundaryAssembler::new( + BoundaryAssembler::new( AdjointDoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0, - ); - - assembler.assemble(&mut output, trial_space, test_space); - - output + ) } /// Assembler for the Laplace hypersingular operator. - pub fn laplace_hypersingular< - T: RlstScalar + MatrixInverse, - Space: FunctionSpace + Sync, - >( - trial_space: &Space, - test_space: &Space, + pub fn laplace_hypersingular + MatrixInverse>( options: &BoundaryAssemblerOptions, - ) -> DynamicArray { - let nrows = trial_space.global_size(); - let ncols = test_space.global_size(); - - let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); - + ) -> BoundaryAssembler, Laplace3dKernel> { let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); - let assembler = BoundaryAssembler::new( + BoundaryAssembler::new( HypersingularCurlCurlBoundaryIntegrand::new(), kernel, options, 4, 1, - ); - - assembler.assemble(&mut output, trial_space, test_space); - - output + ) } } diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs index b421ee15..d9628b1c 100644 --- a/tests/compare_to_bempp_cl.rs +++ b/tests/compare_to_bempp_cl.rs @@ -21,7 +21,7 @@ fn test_laplace_single_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = laplace_single_layer(&space, &space, &options); + let matrix = laplace_single_layer(&options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -41,7 +41,7 @@ fn test_laplace_double_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = laplace_double_layer(&space, &space, &options); + let matrix = laplace_double_layer(&options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -60,7 +60,7 @@ fn test_laplace_adjoint_double_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = laplace_adjoint_double_layer(&space, &space, &options); + let matrix = laplace_adjoint_double_layer(&options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -80,7 +80,7 @@ fn test_laplace_hypersingular_p1_p1() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = laplace_hypersingular(&space, &space, &options); + let matrix = laplace_hypersingular(&options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -106,7 +106,7 @@ fn test_helmholtz_single_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = helmholtz_single_layer(3.0, &space, &space, &options); + let matrix = helmholtz_single_layer(3.0, &options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -126,7 +126,7 @@ fn test_helmholtz_double_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = helmholtz_double_layer(3.0, &space, &space, &options); + let matrix = helmholtz_double_layer(3.0, &options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -145,7 +145,7 @@ fn test_helmholtz_adjoint_double_layer_dp0_dp0() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = helmholtz_adjoint_double_layer(3.0, &space, &space, &options); + let matrix = helmholtz_adjoint_double_layer(3.0, &options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip] @@ -165,7 +165,7 @@ fn test_helmholtz_hypersingular_p1_p1() { let space = SerialFunctionSpace::new(&grid, &element); let options = BoundaryAssemblerOptions::default(); - let matrix = helmholtz_hypersingular(3.0, &space, &space, &options); + let matrix = helmholtz_hypersingular(3.0, &options).assemble(&space, &space); // Compare to result from bempp-cl #[rustfmt::skip]