Skip to content

Commit

Permalink
WIP: Refactor assembler
Browse files Browse the repository at this point in the history
  • Loading branch information
@tbetcke
tbetcke committed Nov 5, 2024
1 parent 821cad8 commit 269e6fa
Show file tree
Hide file tree
Showing 5 changed files with 53 additions and 166 deletions.
2 changes: 1 addition & 1 deletion examples/assembly.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down
22 changes: 9 additions & 13 deletions src/assembly/boundary/assemblers.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -122,26 +122,20 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, K:
}

/// Assemble into a dense matrix.
pub fn assemble<
Space: FunctionSpace<T = T>,
Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut<Item = T>,
>(
pub fn assemble<Space: FunctionSpace<T = T>>(
&self,
output: &mut Array2,
trial_space: &Space,
test_space: &Space,
) {
) -> DynamicArray<T, 2> {
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(),
Expand All @@ -164,6 +158,8 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand<T = T>, 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
Expand Down
92 changes: 21 additions & 71 deletions src/helmholtz.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,109 +22,63 @@ pub mod assembler {
};

/// Assembler for the Helmholtz single layer operator.
pub fn helmholtz_single_layer<
T: RlstScalar<Complex = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
pub fn helmholtz_single_layer<T: RlstScalar<Complex = T> + MatrixInverse>(
wavenumber: T::Real,
trial_space: &Space,
test_space: &Space,
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, SingleLayerBoundaryIntegrand<T>, Helmholtz3dKernel<T>> {
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<Complex = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
pub fn helmholtz_double_layer<T: RlstScalar<Complex = T> + MatrixInverse>(
wavenumber: T::Real,
trial_space: &Space,
test_space: &Space,
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, DoubleLayerBoundaryIntegrand<T>, Helmholtz3dKernel<T>> {
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<Complex = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
pub fn helmholtz_adjoint_double_layer<T: RlstScalar<Complex = T> + MatrixInverse>(
wavenumber: T::Real,
trial_space: &Space,
test_space: &Space,
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, AdjointDoubleLayerBoundaryIntegrand<T>, Helmholtz3dKernel<T>> {
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<Complex = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
pub fn helmholtz_hypersingular<T: RlstScalar<Complex = T> + MatrixInverse>(
wavenumber: T::Real,
trial_space: &Space,
test_space: &Space,
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<

Check failure on line 73 in src/helmholtz.rs

View workflow job for this annotation

GitHub Actions / Rust style checks 

very complex type used. Consider factoring parts into `type` definitions
T,
BoundaryIntegrandSum<
T,
HypersingularCurlCurlBoundaryIntegrand<T>,
BoundaryIntegrandTimesScalar<T, HypersingularNormalNormalBoundaryIntegrand<T>>,
>,
Helmholtz3dKernel<T>,
> {
let kernel = KernelEvaluator::new(
Helmholtz3dKernel::new(wavenumber),
GreenKernelEvalType::ValueDeriv,
Expand All @@ -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)
}
}
87 changes: 14 additions & 73 deletions src/laplace.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,108 +20,49 @@ pub mod assembler {
};

/// Assembler for the Laplace single layer operator.
pub fn laplace_single_layer<
T: RlstScalar<Real = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
trial_space: &Space,
test_space: &Space,
pub fn laplace_single_layer<T: RlstScalar<Real = T> + MatrixInverse>(
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, SingleLayerBoundaryIntegrand<T>, Laplace3dKernel<T>> {
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<Real = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
trial_space: &Space,
test_space: &Space,
pub fn laplace_double_layer<T: RlstScalar<Real = T> + MatrixInverse>(
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, DoubleLayerBoundaryIntegrand<T>, Laplace3dKernel<T>> {
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<Real = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
trial_space: &Space,
test_space: &Space,
pub fn laplace_adjoint_double_layer<T: RlstScalar<Real = T> + MatrixInverse>(
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, AdjointDoubleLayerBoundaryIntegrand<T>, Laplace3dKernel<T>> {
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<Real = T> + MatrixInverse,
Space: FunctionSpace<T = T> + Sync,
>(
trial_space: &Space,
test_space: &Space,
pub fn laplace_hypersingular<T: RlstScalar<Real = T> + MatrixInverse>(
options: &BoundaryAssemblerOptions,
) -> DynamicArray<T, 2> {
let nrows = trial_space.global_size();
let ncols = test_space.global_size();

let mut output = rlst_dynamic_array2!(T, [nrows, ncols]);

) -> BoundaryAssembler<T, HypersingularCurlCurlBoundaryIntegrand<T>, Laplace3dKernel<T>> {
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
)
}
}
16 changes: 8 additions & 8 deletions tests/compare_to_bempp_cl.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand Down

0 comments on commit 269e6fa

Please sign in to comment.